Stochastic models can be developed to perform ensemble forecasts of geophysical fluid dynamical systems to more efficiently handle subgrid parametrizations. As numerical simulations do not usually resolve all temporal scales, solutions can be found by decomposing the velocity into a smooth resolved component and an unresolved one, uncorrelated in time and inhomogeneous in space. In turn, this new point of view changes the usual interpretation of transport and fundamental conservation laws (mass, momentum and 1st principle). Indeed, following a stochastic version of the Reynolds transport theorem (Mémin, 2014), three terms naturally emerge: a multiplicative noise, an anisotropic and inhomogeneous diffusion, and a drift correction. For instance, neglecting diabatic effects implies the conservation of the temperature T along the flow as : DT/Dt=0, with the material derivative D/Dt, now understood in a stochastic sense. As such, the rigorous derivation of this model relates the random forcing to the subgrid parametrization. This further ensures important properties such as energy conservation and also simplifies the derivation of ensemble forecasts in defining a clear mathematical framework. Using this new framework, stochastic versions of geophysical models have been derived, namely: Navier-Stokes equations in a rotating frame, the Boussinesq approximation, Quasi-Geostrophy (QG) and Surface Quasi-Geostrophy (SQG). Depending on the amount of randomness, the QG approximation can lead to two different models. With moderate uncertainty, the horizontal transport of the Potential Vorticity (PV), in the interior of the fluid, has 3 sources terms. For homogeneous turbulence, two of them, disappear. The remaining term, a noise uncorrelated in time, encodes the interactions between the resolved and the unresolved velocity gradient tensors. The ensuing SQG model is then derived assuming a zero PV in the interior, as in the deterministic case (Held et al., 1995). Yet, the buoyancy transport at the surface of the fluid has to be understood in the stochastic sense. Increasing to strong uncertainty, the QG approximation further leads to a vanishing PV in the interior. As such, a classical SQG relationship remains. However, for that case, an ageostrophic component appears, contributing to intensify asymmetries between cyclones and anticyclones, as in the SQG+ model (Hakim et al., 2002). Frontogenesis occurs on the cold side of fronts, on uncertain locations, and frontolysis smooths the warm side of fronts. As obtained, the flow topology is changed by the noise and the divergent component, and the symmetry breaks more rapidly. Larger the uncertainty, faster the symmetry is broken. Simulating an ensemble of realizations then enables to track the different possible topologies and all the bifurcations of the system.