V Resseguier
email: valentin.resseguier@inria.fr

E Mémin

B Chapron

Geophysical flows under

Keywords: stochastic sub-grid parameterization, uncertainty quantification, ensemble forecasts

Introduction

Ensemble forecasting and filtering are widely used in geophysical sciences for forecasting and climate projection. In practice, dynamical models are randomized through their initial conditions and a Gaussian error model, and are generally found to be underdispersive [START_REF] Mitchell | Data assimilation in slow-fast systems using homogenized climate models[END_REF][START_REF] Gottwald | The role of additive and multiplicative noise in filtering complex dynamical systems[END_REF][START_REF] Berner | Model uncertainty in a mesoscale ensemble prediction system: Stochastic versus multiphysics representations[END_REF][START_REF] Snyder | Performance bounds for particle filters using the optimal proposal[END_REF] with a low variance. As a consequence, errors are underestimated and observations are hardly taken into account. Corrections are considered by incorporating inflation procedures or hyperprior to increase the variance of ensemble Kalman filters [START_REF] Anderson | A Monte Carlo implementation of the nonlinear filtering problem to produce ensemble assimilations and forecasts[END_REF][START_REF] Bocquet | Expanding the validity of the ensemble kalman filter without the intrinsic need for inflation[END_REF]. However, such corrections do not provide an accurate spatial localization of the errors.

Another difficulty of ensemble methods lies in the huge dimensions of the involved state spaces. For obvious computational reasons, ensembles for geophysical applications appear constrained and limited to small sizes. It thus becomes primordial to build strategies to best track the most likely dynamical events. From this point of view, ensemble simulations and stochastic dynamics have clear advantages over the deterministic models.

The simplest random models are defined from Langevin equations with linear damping and additive isotropic Gaussian noise, as, for instance, the linear inverse models (e.g. [START_REF] Penland | A balance condition for stochastic numerical models with application to the el nino-southern oscillation[END_REF][START_REF] Penland | The optimal growth of tropical sea surface temperature anomalies[END_REF], or the Eddy-Damped Quasi Normal Markovian (EDQNM) models (e.g. [START_REF] Orszag | Analytical theories of turbulence[END_REF][START_REF] Leith | Atmospheric predictability and two-dimensional turbulence[END_REF][START_REF] Jeffrey R Chasnov | Simulation of the kolmogorov inertial subrange using an improved subgrid model[END_REF]. Among other empirical stochastic models, the Stochastic Kinetic Energy Backscatter (SKEBS) [START_REF] Shutts | A kinetic energy backscatter algorithm for use in ensemble prediction systems[END_REF][START_REF] Berner | A spectral stochastic kinetic energy backscatter scheme and its impact on flow-dependent predictability in the ecmwf ensemble prediction system[END_REF][START_REF] Berner | Model uncertainty in a mesoscale ensemble prediction system: Stochastic versus multiphysics representations[END_REF] and the Stochastic Perturbed Physics Tendency scheme (SPPT) [START_REF] Buizza | Stochastic representation of model uncertainties in the ecmwf ensemble prediction system[END_REF] introduce correlated multiplicative noises. SPPT and SKEBS methods have been successfully applied in operational weather forecast centers [START_REF] Franzke | Stochastic climate theory and modeling[END_REF]. To target highly non-Gaussian distribution of fluid dynamics properties, an attractive path is to infer randomness from physics [START_REF] Berner | Stochastic parameterization: Towards a new view of weather and climate models[END_REF]. For this purpose, the time-scale separation assumption is convenient. [START_REF] Hasselmann | Stochastic climate models. part i: theory[END_REF] already relied on it for geophysical fluid dynamics. This assumption is the foundation of averaging and homogenization theories [START_REF] Kurtz | A limit theorem for perturbed operator semigroups with applications to random evolutions[END_REF][START_REF] Gc Papanicolaou | Asymptotic theory of mixing stochastic ordinary differential equations[END_REF][START_REF] Givon | Extracting macroscopic dynamics: model problems and algorithms[END_REF][START_REF] Gottwald | Homogenization for deterministic maps and multiplicative noise[END_REF][START_REF] Mitchell | Data assimilation in slow-fast systems using homogenized climate models[END_REF][START_REF] Gottwald | The role of additive and multiplicative noise in filtering complex dynamical systems[END_REF][START_REF] Franzke | Stochastic climate theory and modeling[END_REF][START_REF] Gottwald | Stochastic climate theory[END_REF]. A successful application of homogenization theory in geophysics is the MTV algorithms (Majda et al., 1999[START_REF] Majda | A mathematical framework for stochastic climate models[END_REF][START_REF] Franzke | Low-order stochastic mode reduction for a realistic barotropic model climate[END_REF][START_REF] Majda | An applied mathematics perspective on stochastic modelling for climate[END_REF]. The homogenized dynamics is cubic with correlated additive and multiplicative (CAM) noises. This noise structure is able to produce intermittency and extreme events. In practice, the non-linearity of the small-scale equation (fast dynamics) is conveniently replaced by a noise and a damping terms before the homogenization procedure. Noise statistics are estimated from data, with Gaussian assumptions.

In Resseguier et al. (2017a), following [START_REF] Mémin | Fluid flow dynamics under location uncertainty[END_REF], another approach has been considered to help derive models under location uncertainty based on stochastic calculus and the Ito-Wentzell formula [START_REF] Kunita | Stochastic flows and stochastic differential equations[END_REF]. [START_REF] Mikulevicius | Stochastic Navier-Stokes equations for turbulent flows[END_REF] and [START_REF] Flandoli | The interaction between noise and transport mechanisms in pdes[END_REF] already introduce this methodology. Yet, their works mostly focused on pure mathematical aims: existence and uniqueness of SPDE solutions. For our more practical purpose, the large-scale is understood as sub-sampled in time, and the remaining small-scale velocity component is then considered as uncorrelated in time.

Starting with the definition of the revised transport under location uncertainty (section 2), developments are then carried out to derive and analyze the stochastic versions of Quasi-Geostrophy (QG) and Surface Quasi-Geostrophy (SQG) models with a moderate influence of sub-grid terms (section 3). Numerical results highlight the potential of these models under location uncertainty, especially for ensemble forecast (Section 4).

Models under location uncertainty

This section briefly outlines main theoretical results discussed in Resseguier et al. (2017a). The velocity is decomposed between a possibly random large-scale component, w, and a time-uncorrelated component, σ Ḃ. The latter is Gaussian, correlated in space with possible inhomogeneities and anisotropy. Hereafter, this unresolved velocity component will further be assumed to be solenoidal. To parameterize those spatial correlations, we apply an infinite-dimensional linear operator, σ, to a d-dimensional space-time white noise1 , Ḃ.

In time, the velocity is irregular. The material derivative, D t , is then changed. In most cases, it coincides with the stochastic transport operator, D t , defined for every field, Θ, as follows:

D t Θ = d t Θ = Θ(x,t+dt)-Θ(x,t) Time increment + (w dt + σdB t ) • ∇Θ Advection -∇ • 1 2 a∇Θ Diffusion dt, (1)
where the time increment term d t Θ is used in place of the partial time derivative, as Θ is in general non-differentiable. The diffusion coefficient matrix, a, is solely defined by the one-point one-time covariance of the unresolved displacement per unit of time:

a = σσ T = E σdB t (σdB t ) T dt , (2)
and the modified drift is given by

w = w -1 2 (∇ • a) T . (3)
For a divergent small-scale velocity, this drift would involve an additional component (see equation ( 4) of Resseguier et al. (2017a)). With this modified material derivative (1), the transport equations under location uncertainty involve three new terms: a modification of the large-scale advection (w instead of w), an inhomogeneous and anisotropic diffusion and a multiplicative noise. This random forcing is directly related to the advection by the unresolved velocity.

For incompressible flows (∇•w = 0), the energy of any tracer, Θ, is conserved for each realization:

d Ω Θ 2 = 0, ( 4
)
where Ω is the spatial domain. This still holds for active tracers. The diffusion dissipates as much energy as the multiplicative noise is injecting it in the system. In particular, the (ensemble) mean of the energy, E Ω Θ 2 , is conserved. This results ensures a constant balance between the energy of the mean and the (ensemble) variance. The energy fluxes in these stochastic models are more thoroughly described in Resseguier et al. (2017a).

A random version of the Reynolds transport theorem can further be derived [START_REF] Mémin | Fluid flow dynamics under location uncertainty[END_REF]Resseguier et al., 2017a). From this theorem, usual conservation of mechanics (mass, linear momentum, energy and amount of substance) can be expressed in a stochastic sense. Random Navier-Stokes and Boussinesq models can then be derived. This last model describes the stochastic transports of velocity and density anomaly, as well as incompressibility conditions.

Mesoscales under moderate uncertainty

To simplify the stochastic Boussinesq model of Resseguier et al. (2017a), Quasi-Geostrophic (QG) models are developed for large horizontal length scales, L, such as:

1 Bu = Fr Ro 2 = L L d 2 ∼ 1 and 1 Ro = Lf 0 U 1, ( 5
)
where U is the horizontal velocity scale, L d = N h f is the Rossby deformation radius, N is the stratification (Brunt-Väisälä frequency) and h is the characteristic vertical length scale. The Rossby deformation radius explicitly defines the mesoscale range, over which both kinetic and buoyancy effects are important, and strongly interact. In the following, both differential operators Del, ∇, and Laplacian, ∆, represent 2D operators.

Specific scaling assumptions

Hereafter, we explicit scaling assumptions to derive the non-dimensional version of the stochastic Boussinesq model.

Quadratic variation scaling

Besides traditional ones, another dimensionless number, Υ , is introduced to relate the large-scale kinetic energy to the energy dissipation due to the horizontal small-scale random component. In the following, σ H• stands for the horizontal component of σ, a H for σ H• σ T H• and A u for its scaling. The new dimensionless number is defined by:

Υ = U L A u = U 2 A u /T . ( 6
)
This number compares horizontal advective and diffusive terms in the momentum and buoyancy equations. This number can also be related to the ratio between the Mean Kinetic Energy (MKE), U 2 , and the Turbulent Kinetic Energy (TKE), A u /T σ , where T σ is the small-scale correlation time. This reads:

Υ = 1 M KE T KE , (7)
where = T σ /T is the ratio of the small-scale to the large-scale correlation times. This parameter, , is central in homogenization and averaging methods (Majda et al., 1999;[START_REF] Givon | Extracting macroscopic dynamics: model problems and algorithms[END_REF][START_REF] Gottwald | Homogenization for deterministic maps and multiplicative noise[END_REF]. The number Ro/Υ can then be stated to measure the ratio between sub-grid terms and the Coriolis force. In the usual deterministic case and the limit of small Rossby number, the predominant terms of the horizontal momentum equation then correspond to the geostrophic balance. In the stochastic case, this balance also applies from weak (Υ 1) to moderate (Υ ∼ 1) uncertainty. However, if Υ/Ro is close enough to O(1), this geostrophic balance is modified due to the diffusion effects introduced by the small-scale random velocity. Hereafter, developments focus on the moderate uncertainty case. Resseguier et al. (2017b) deals with the strong uncertainty case.

To evaluate Υ for a given flow at a given scale, eddy viscosity or diffusivity values help the determination of A u . [START_REF] Boccaletti | Mixed layer instabilities and restratification[END_REF] give some examples of canonical values. Then, the typical resolved velocity and length scale lead to Υ . If no canonical values are known, absolute diffusivity or similar mixing diagnoses could be measured [START_REF] Keating | Diagnosing lateral mixing in the upper ocean with virtual tracers: Spatial and temporal resolution dependence[END_REF] as a proxy of the variance tensor.

Vertical unresolved velocity

The scaling to compare vertical to horizontal unresolved velocities is also considered:

(σdB t ) z (σdB t ) H ∼ Ro Bu D, (8)
where D = h L is the aspect ratio and the subscript H indicates horizontal coordinates. This scaling can be derived from the ω-equation [START_REF] Giordani | Advanced insights into sources of vertical velocity in the ocean[END_REF]. For any velocity u = (u H , w) T , which scales as (U, U, W) T , this equation reads

f 2 0 ∂ 2 z w + N 2 ∆w = ∇ • Q ≈ -∇ • (∇u T H ∇b) ≈ -f 0 ∇ • ∇u T H ∂ z u ⊥ H , (9)
where b stands for the buoyancy variable and Q for the so-called Q-vector. In its non-dimentional version, the ω-equation reads:

W U ∂ 2 z w + Bu∆w ≈ DRo∇ • Q. ( 10
)
At planetary scales, Burger number is small and the rotation dominates the stratification, W/U ∼ DRo. At smaller scales, with a larger Burger number, the stratification dominates the rotation, W/U ∼ DRo/Bu. For the small-scale velocity σ Ḃ, the latter is thus more relevant. Note that the angle between the small-scale component and the horizontal one can be assumed to be constrained by the angle between the isopycnical and the horizontal plane. Invoked to describe baroclinic instabilities theory, this statement helps to specify the anisotropy of the eddy diffusivity [START_REF] Vallis | Atmospheric and oceanic fluid dynamics: fundamentals and large-scale circulation[END_REF]. The argument of the orientation of the eddies activity with isentropic surfaces and the related mixing is also supported by several other authors [START_REF] Gent | Isopycnal mixing in ocean circulation models[END_REF][START_REF] Pierrehumbert | Global chaotic mixing on isentropic surfaces[END_REF].

In the case of QG models, the large and small Burger scaling cases lead to the same result: the unresolved velocity is mainly horizontal.

(σdB t ) z (σdB t ) H ∼ Ro Bu D D. (11)
This is consistent with the assumption of a large stratification, i.e. flat isopycnicals, if we admit that the eddies activity appears preferentially along the isentropic surfaces. As a consequence, the terms (σdB t ) z ∂ z scale as Ro/Bu (σdB t ) H • ∇. In the QG approximation, the scaling of the diffusion and effective advection terms including σ z• are one to two orders smaller (in power of Ro/Bu) than terms involving σ H• . For any function ξ, the vertical diffusion

∂ z (σ z• σ T z• /2 ∂ z ξ) is one order smaller than the horizontal-vertical diffusion term ∇ • (σ H• σ T z• /2 ∂ z ξ)
and two orders smaller than the horizontal diffusion term ∇ • (σ H• σ T H• /2 ∇ξ).

Beta effect

At mid-latitudes, the related term, given by β = ∂ y f , is much smaller than the constant part of the Coriolis frequency. Nevertheless, it can govern a large part of the relative vorticity at large scales.

The following scaling is thus chosen [START_REF] Vallis | Atmospheric and oceanic fluid dynamics: fundamentals and large-scale circulation[END_REF]:

βy ∼ ∇ ⊥ • u ∼ U L = Rof 0 . ( 12
)

Stratified Quasi-Geostrophic model under moderate uncertainty

The moderate uncertainty case corresponds to Υ ∼ 1. Horizontal advective terms and horizontal sub-grid terms are comparable.

Following similar principles as those used to derive the deterministic stratified QG model [START_REF] Vallis | Atmospheric and oceanic fluid dynamics: fundamentals and large-scale circulation[END_REF], a stochastic QG model can be derived (see Appendix B). This QG solution corresponds to the limit of the Boussinesq solution when the Rossby number goes to zero. The resulting potential vorticity (PV), Q, is then found to be conserved, along the horizontal random flow, up to three source terms:

D H t Q = 1 2 i,j∈H ∂ 2 ij ∇ ⊥ a ij • u dt -1 2 ∇ • (∇ • (a H f )) T dt -tr [SJ S σdBt ] , (13)
where the QG PV is:

Q = ∆ψ + f + f 0 N 2 ∂ 2 z ψ, (14)
ψ is the streamfunction, J = 0 -1 1 0 is the π/2 rotation matrix,

S = 1 2 ∇u T + (∇u T ) T and S σdBt = 1 2 ∇(σdB t ) T H + (∇(σdB t ) T H ) T ( 15
)
denotes the strain rate tensor of the horizontal resolved and unresolved velocities, u and (σdB t ) H , respectively. To interpret the source terms, we rather focus on the material derivative of the PV:

D H t Q = ∇ • (α∇ψ)dt -1 2 ∇ • (∇ • (a H f )) T dt -tr [SJ S σdBt ] , (16)
with

α T = i ∇σ T Hi 2 , (17)
which can be decomposed into a symmetric part, positive or negative diffusion of the stream function, and an anti-symmetric part, skew diffusion advection of the stream function. Compared to the traditional QG model, this system includes two smooth (continuous) source/sink terms that depend on the variance tensor, and a random forcing term. The first source term in ( 16) is correlated in time and may decrease or increase the PV energy. This term is due to the spatial variations of both the diffusion coefficient and the drift correction. The second term takes into account interactions between the Coriolis frequency, including beta effects, and inhomogeneous sub-grid eddies. The last source term in ( 16) is a noise term, encoding the interactions between the resolved and the unresolved strain rate tensors. Uncorrelated in time, this noise increases the potential enstrophy along time.

To further understand this source term, let us denote Ξ and Λ the eigenvalues associated with the stable directions (i.e. negative eigenvalue) of the strain rate tensors of the large-scale flow, S, and of the small-scale flow, S σdBt respectively. We note θ, the angle between these two stable directions

-tr [SJ S σdBt ] = 2 ΞΛ >0 sin (2θ) . ( 18
)
The detailed derivation is provided in Appendix B. This random source vanishes when the stable directions of u and (σdB t ) H are aligned or orthogonal. It is maximum and positive (respectively minimum and negative) when there is an angle of π/4 (respectively -π/4) between those directions. Around the local position x, stable and unstable directions of the large-scale velocity define 2 axes and 4 quadrants. As understood, the strain rate tensor does not depend on the local vorticity.

Yet, an hyperbolic deformation will almost resemble a positive vorticity in the upper-left and bottom-right quadrants, and a negative vorticity in the upper-right and bottom-left quadrants. For θ = π/4, the stable direction of the small-scale velocity aligns along the upper-left to bottom-right direction. The small-scale velocity then compresses the flow in this direction and dilates the flow in the orthogonal direction (upper-right to bottom-left). The quadrants associated with a seemingly positive (resp. negative) vorticity are brought closer (resp. farther) to x. Accordingly, the vorticity increases at x. For -π/4, the vorticity would decrease.

Note the dB t factor has been omitted in the right-hand side of equation ( 18). This term remains a linear function of the uncorrelated noise z → dB t (z). Whatever the angle between the stable directions, the source term always has a zero (ensemble) mean and increases the enstrophy since it is a term in dB t . Equation ( 18) could then be used to define the horizontal inhomogeneous smallscale component of the velocity. If the conservation of PV is a strong constraint, this component can indeed be defined to ensure that its stable direction is always along or orthogonal to the stable direction of u.

A two-layer model could also be deduced from equation ( 13) or ( 16). This would help identifying the stochastic parameterization effects on the barotropic and baroclinic modes. In particular, the particular forms of the operator σ able to trigger barotropization effects can bemore efficiently studied.

In the stochastic QG model, the stream function ψ is related to the buoyancy, b, the pressure, p , and the velocity, u, by the usual relations:

b = f 0 ∂ z ψ, p = ρ b f 0 ψ and u = ∇ ⊥ ψ, (19)
where ρ b is the mean (background) density. The horizontal noise term, (σdB t ) H , appearing in both the horizontal stochastic material derivative and in the 2 × 2 horizontal variance tensor, a H , is in geostrophic balance with a pressure component uncorrelated in time. Due to their scaling, the vertical noise and its variance are neglected in the final equations.

For homogeneous turbulence conditions, the transport of PV ( 16) simplifies. The variance tensor becomes constant, the first two source terms disappear, to give

D H t Q = D H t Q = -tr [SJ S σdBt ] . (20)
The transport of the PV (equation ( 13) or ( 20)) determines the dynamics of the fluid interior. Boundary conditions are then necessary to specify completely the dynamics.

Surface Quasi-Geostrophic model under moderate uncertainty

A classical choice considers a vanishing solution in the deep ocean and a buoyancy transport at the surface [START_REF] Vallis | Atmospheric and oceanic fluid dynamics: fundamentals and large-scale circulation[END_REF][START_REF] Lapeyre | Dynamics of the upper oceanic layers in terms of surface quasigeostrophy theory[END_REF]:

ψ -→ z→-∞ 0 and D H t b |z=0 = D H t b |z=0 = 0. ( 21
)
Assuming zero PV in the interior but keeping these boundary conditions leads to the Surface Quasi-Geostrophic model (SQG) [START_REF] Blumen | Uniform potential vorticity flow: part i. theory of wave interactions and twodimensional turbulence[END_REF][START_REF] Held | Surface quasi-geostrophic dynamics[END_REF][START_REF] Lapeyre | Dynamics of the upper oceanic layers in terms of surface quasigeostrophy theory[END_REF][START_REF] Constantin | Formation of strong fronts in the 2-d quasigeostrophic thermal active scalar[END_REF][START_REF] Constantin | Front formation in an active scalar equation[END_REF][START_REF] Constantin | New numerical results for the surface quasi-geostrophic equation[END_REF]. Under the stochastic framework, the derivation is similar. The PV is indeed identical to the classical one (see equation ( 14)), assuming zero PV in the interior and vanishing solution as z → -∞ unsurprisingly yields the same SQG relationship:

b = N k ψ. ( 22
)
The top boundary condition, equation ( 21), provides an evolution equation, namely the horizontal transport of surface buoyancy, in the stochastic sense:

D H t b = 0. ( 23
)
The time-uncorrelated component of the velocity, σ Ḃ, is divergence-free. Its inhomogeneous and anisotropic spatial covariance has then to be specified. The time-correlated component of the velocity is also divergence-free, with a stream function specified by the SQG relation ( 22). The buoyancy is randomly advected, and the resulting smooth velocity component is random as well.

Summary

For simplified models, stochastic versions are derived for scaling assumptions related to the sub-grid terms. For moderate uncertainty, the PV is transported along the random flow up to three source terms. The first one, smooth in time, is due to spatial variations of the inhomogeneous diffusion and the drift correction. The second one, also smooth, encodes the interaction between inhomogeneous turbulence and Coriolis frequency. These terms disappear for an homogeneous turbulence. The last term, a time-uncorrelated multiplicative noise, involves the large-scale and the small-scale strain rate tensors. It is a source of potential enstrophy and its instantaneous value depends on the angle between the large-scale and small-scale stable directions. Assuming zero PV in the interior, a SQG model follows from this QG model.

Numerical results

We focus on this SQG M U model (3.3). A high-resolution deterministic SQG simulation provides a reference. The MATLAB R codes are available online (http://vressegu.github.io/sqgmu). Numerical results are analyzed in terms of the resolution gains (when a single realization is simulated) and the potential for ensemble forecasting in estimating spatial and spectral reconstruction errors (for an ensemble of realizations).

Test flow

The initial conditions defining the test flow, Figure 1, consist of a spatially smooth buoyancy field with two warm elliptical anticyclones and two cold elliptical cyclones given by:

b(x, t = 0) = F x - 250 km 250 km + F x - 750 km 250 km -F x - 250 km 750 km -F x - 750 km 750 km , (24a)
with

F (x) = B 0 exp -1 2 x 2 σ 2 x + y 2 σ 2 y and σ x = 67 km, σ y = 133 km. (24b)
The size of the vortices is of order of the Rossby radius L d . The buoyancy and the stratification have been set with: B 0 = 10 -3 m.s -2 and N = 3f 0 . The Coriolis frequency is set to 1.028 × 10 -4 s -1 (45 • N). Periodic boundaries conditions are considered. The deterministic high-resolution SQG reference model is associated with a spatial mesh grid of 512 2 points, whereas the low-resolution (deterministic or stochastic) SQG models are run on 128 2 points. The simulations have been performed through a pseudo-spectral code in space. As for the temporal discrete scheme the deterministic simulation relies on a fourth-order Runge-Kutta scheme, whereas the stochastic ones are based on an Euler-Maruyama scheme [START_REF] Kloeden | Numerical Solution of Stochastic Differential Equations[END_REF]. For our application, the weak precision of this scheme is balanced by the use of a small time step. In all the simulations (deterministic and random, high-resolution and low-resolution), a standard hyperviscosity model is used:

D t b = α hv ∆ 4 b dt, ( 25
)
with a coefficient α hv = (5 × 10 29 m 8 .s -1 )M -8

x where M x denotes the meshgrid size (i.e. 128 or 512).

Figure 1 displays the high-resolution buoyancy field at t = 0, 5, 10, 13, 15, 16, 20 and 30 days. During the first ten days, the vortices turn with slight deformation. Vortices of the same sign have their tails that draw closer. This creates high shears around four saddle points located at (x, y) = (0, 250), (500, 250), (0, 750) and (500, 750) (in km). A strong non-linearity in the neighborhood of a saddle point has been identified to become a major source of instability [START_REF] Constantin | Formation of strong fronts in the 2-d quasigeostrophic thermal active scalar[END_REF][START_REF] Constantin | Front formation in an active scalar equation[END_REF][START_REF] Constantin | New numerical results for the surface quasi-geostrophic equation[END_REF]. In our case, this effect is weak but yields an effective creation of turbulence 10 days later. Shears create long and fine filaments, wrapping around the vortices until the 15 th day. At this time, the filaments become unstable, break and a so-called "pearl-necklace" appears, characteristic of the SQG model, days 17-18 in the simulation. These small vortices are then ejected from their orbits. Between days 17 th and 25 th , they interact with the large vortices, the filaments and other small vortices, to create a fully-developed SQG turbulence orbiting around the four large vortices.

Simulation of the random velocity

To simulate the SQG MU model ( 22)-( 23), the covariance of the unresolved velocity σ Ḃ must be specified. As this unresolved velocity field is assumed divergence-free, we introduce the following stream function linear operator, ψ σ , and its kernel, ψσ :

σ H (x)dB t = ∇ ⊥ ψ σ (x)dB t , = Ω dz ∇ ⊥ x ψσ (x, z)dB t (z). ( 26
)
As such, a single cylindrical Wiener process, B t , is sufficient to sample our Gaussian process. This is specific to two-dimensional domains. In 3D, a vector of 3 independent I d -cylindrical Wiener processes, and a projection operator on the divergence-free vector space or a curl must be considered to simulate an isotropic small-scale velocity [START_REF] Mémin | Fluid flow dynamics under location uncertainty[END_REF]. For a divergent unresolved velocity, equation ( 26) can additionally involve the gradient of a random potential, ∇ ψσ dB t . Then, similar to the Kraichnan's model, a solenoidal homogeneous field can be considered: [START_REF] Kraichnan | Small-scale structure of a scalar field convected by turbulence[END_REF][START_REF] Kraichnan | Anomalous scaling of a randomly advected passive scalar[END_REF][START_REF] Gawedzky | Anomalous scaling of the passive scalar[END_REF][START_REF] Majda | Simplified models for turbulent diffusion: Theory, numerical modelling, and physical phenomena[END_REF]:

High-resolution buoyancy

σ H (x)dB t = Ω dz ∇ ⊥ x ψσ (x -z)dB t (z) = ∇ ⊥ ψσ dB t (x). ( 27
)
where denotes a convolution. Although spatially inhomogeneous field would be more physically relevant, homogeneity greatly simplifies the random field simulation. Indeed, homogeneity in physical space implies independence between the Fourier modes

σ H Ḃ(k) = ik ⊥ ψσ (k) Ḃ(k), ( 28
) in the half-space k ∈ (R × R + * ) ∪ (R + × {0}).
Thus, the small-scale velocity can be conveniently specified from its omnidirectional spectrum:

k → 1 µ(Ω) E [0,2π] dθ k k σ H Ḃ(k) 2 = 2π ∆t k 3 ψσ ( k ) 2 , ( 29
)
where µ(Ω) is the surface of the spatial domain Ω, θ k is the angle of the wave-vector k and ∆t the simulation time-step. Consistent with SQG turbulence, the omni-directional spectrum slope, denoted s, is fixed to -5/3. For 2D Euler equations, the slope would be set to -3. If the small scales spectrum slope is unknown, the spectrum slope of the resolve scales -estimated on linemay enable to specify s through a scale similarity assumption. The unresolved velocity should be energetic only where the dynamics cannot be properly resolved. Consequently, we apply to the spectrum a smooth band-pass filter, f BP , which has non-zero values between two wavenumbers κ min and κ max . The parameter κ min is inversely related to the spatial correlation length of the unresolved component. In practice, we set κ max to the theoretical resolution, π/∆x, and κ min to the effective resolution (hereafter κ min = κ max /2). Figure 2 illustrates this spectrum specification.

The small scales' energy is specified by the diffusion coefficient a H and the simulation time step:

E σ H Ḃ σ H Ḃ T = 1 ∆t a H = 1 ∆t a H 0 0 a H . ( 30
)
The diagonal structure of the variance tensor is due both to incompressiblity and isotropy. The scalar variance tensor, a H , is similar to an eddy viscosity coefficient. So, a typical value of eddy viscosity used in practice is a good proxy to setup this parameter. Otherwise, this parameter can be tuned. For this paper, it is set to 9 m 2 .s -1 . The time step depends itself, through the CFL conditions, on both the spatial resolution and the maximum magnitude of the resolved velocity. Finally, equation ( 28) writes:

σ H Ḃ(k) = A √ ∆t ik ⊥ f BP ( k ) k -α dB t √ ∆t (k) with s = 3 -2α = - 5 3 , ( 31
)
where A is a constant to ensure E σ H Ḃ 2 = 2 a H ∆t (see equation ( 30) above), dB t is the spatial Fourier transform of dB t , with dBt √ ∆t , a discrete scalar white noise process of unit variance in space and time. To sample the small-scale velocity, we first sample dB t / √ ∆t, to get dB t / √ ∆t, and finally σ H Ḃ(k) with the above equation.

Resolution gain on a single simulation

In Figure 3, the buoyancy field and its spectrum for low resolution SQG MU and deterministic SQG simulations are displayed for the day 17 th . For the spectrum plots (right column), the slope -5/3 is superimposed. While the spectrum tail of the SQG model falls slightly before the stochastic one, the most significant gain is observed in the spatial domain, i.e. in the phase of the tracer. Indeed, the SQG MU buoyancy field exhibits pearl-necklaces, only obtained at higher resolution. The low-resolved SQG simulation only generates smooth and stable filaments. Though small-scale energy distribution remains similar for both low-resolved models, the phase of the stochastic tracer is more accurate. This may seem surprising since the unresolved velocity, σ H Ḃ, is defined in a loose way, through its spectrum, without prescribing the nature of its phase. However, the noise is multiplicative, and the random forcing, -(σ H Ḃ) • ∇b, does implicitly take into account the tracer phase.

Note, within the stochastic framework, the diffusion coefficient is explicitly related to the noise variance. If the small-scale velocity is set to a magnitude three times smaller than the one prescribed by the diffusion coefficient a H /2, the tracer field becomes quickly too smooth (see Figure 4). Conversely, if the small-scale velocity is set to a magnitude three times larger than dictated by the stochastic transport model, the tracer field becomes rapidly too noisy. This is visible both in the spatial and Fourier spaces (Figure 4). The stochastic transport model thus imposes a correct balance between noise and diffusion. ) and its spectrum (m 2 .s -4 /(rad.m -1 )) at the 17 th day of advection for SQG MU at resolution 128 2 (top), SQG at resolution 512 2 (middle) and at resolution 128 2 (bottom). Unlike SQG MU , the low-resolved SQG simulation diffuses the "pearl necklaces", noticeable only at higher resolution (Colour online). ) and its spectrum (m2 .s -4 /(rad.m -1 )) at the t = 17 th day of advection for the SQG MU model with a small-scale velocity component three times weaker than the one prescribed by the diffusion coefficient a H (top), with the correct amount of small-scale energy (middle) and a small-scale velocity three times higher than the model diffusion. If the prescribed balance between noise and diffusion is not met the tracer field becomes quickly too smooth or too noisy (Colour online). and filtering. Straightforwardly, an ensemble of independently randomly forced realizations {b (i) |i = 1, • • • , N e } of tracer b can be simulated according to the SPDE (23). The probability density function and all the statistical moments of the simulated tracer can then be approximated. For instance, the (ensemble) mean of the buoyancy is a spatio-temporal field defined by:

Ensemble forecasts

E(b)(x, t) ≈ Ê(b)(x, t) = 1 N e Ne i=1 b (i) (x, t), (32)
where N e denotes the ensemble size. This is in essence a Monte-Carlo Markov Chain (MCMC) simulation. The ensemble size is deliberately kept small 2 in order to assess the proposed stochastic framework skills.

We compare the ensemble bias with the estimated error provided by the ensemble itself. The bias corresponds to the discrepancy between the tracer ensemble mean and the SQG simulation at high resolution3 (512 2 ).

Our reference is deterministic since the initial condition is perfectly known and the target dynamics is deterministic, as the real ocean dynamics. The partial knowledge of initial conditions is a complementary issue not addressed in this paper. The reference being deterministic, the bias represents both the error of the mean and the mean of the error:

E{b} -b ref = E{ }, (33)
where = b -b ref stands for the (random) error. We denote by e the absolute value of this bias.

Another error metric could be the Root Mean Square Error (RMSE), E{ 2 }. Yet slightly larger, it is found to have similar spatial and spectral distributions (not shown).

The estimated error, denoted est , is set to 1.96 times the ensemble standard deviation. This specific value corresponds to the (Gaussian) 95% confidence interval. Although the tracer distribution is not Gaussian, this value provides an accurate conventional error estimate:

2 est (x, t) = (1.96) 2 V ar(b) = (1.96) 2 1 N e -1 Ne i=1 b (i) -Ê(b) 2 . ( 34
)
As this error depends on time and space, several comparisons are performed at several distinct times in both the spatial and Fourier domains. In Figures 5 and6, the absolute value of spatial fields ( 33) and (34) (i.e. e and est ) are compared at days 10, 13, 15, 17, 20 and 25. As obtained, the SQG M U model enables the ensemble to predict the positions and the amplitudes of its own errors with a very good accuracy.

To compare the spread-error consistency of the proposed model, a more classical type of random simulation is considered. An ensemble of the same size is initialized with random perturbations of the initial conditions (24a). The perturbations are assumed to be homogeneous, isotropic, Gaussian and are sampled from a (-5/3) spectrum restricted to the small spatial scales, as shown in Figure 7. Then, the ensemble is forecast with the deterministic SQG model.

Figures 8 and9 represent the spectrum of the errors. The blue and red lines with crosses stand for the spectrum of the bias absolute value, e, of the SQG M U with deterministic initial conditions and of the SQG model with random initial conditions, respectively. | f ( κ)

| 2 κ( r a d . m -1 )

D e t e r mi n i st i c an d r an d om p ar t s of t h e i n i t i al c on d i t i on

Figure 7: Spectrum (m 2 .s -4 /(rad.m -1 )), at the initial time, of the mean buoyancy, in blue, spectrum of its random perturbation, in red, and slope -5/3 in black. The initial perturbation is restricted to a narrow spectral band. This random initial condition has been used to simulate an ensemble with the deterministic SQG model (Colour online).

Deterministic and stochastic models have close distribution of errors over the scales, although the SQG M U ensemble mean generally leads to lower errors than the SQG ensemble mean.

The blue line with circles denotes the spectrum of the SQG M U ensemble estimated error, est . As a benchmark, we superimposed the spectrum of the same estimator, est , but simulated with the usual model (red curve with circles). This estimation is dramatically underestimated. It is generally one order of magnitude smaller that the real error. To reduce this drawback, a solution would be to multiply by 10 the perturbations of the initial condition. However, this solution introduces strong errors on the realizations (not shown). Their small-scale errors are generally one order of magnitude larger than the ones of our model. These realizations of the deterministic model remain far from the reference for about ten days. On the contrary, the SQG M U predicts the correct spectral distribution of errors at each time, except at very small-scales, and each of its realizations are accurate as shown in the previous subsection. Let us note however that most of the errors are concentrated at large scales.

SQG M U thus appears to provide a relevant ensemble of realizations, as it enables us to estimate the amplitude of its own error with a good accuracy both in the spatial and spectral domains.

With such an ensemble of realizations, it is now possible to analyze the spatio-temporal evolution of the statistical moments. In Figure 10, we plotted the ensemble tracer mean and variance for t = 17, 20 and 30 days of advection. As expected, the mean field is more smooth than the realizations (see Figure 4 for comparison at t = 17 days). One realization provides a more realistic field than the mean from a topological point of view. Indeed, the realization exhibits physically relevant smallscale structures. Nevertheless, those structures have uncertain shapes and positions. Therefore, on average, the mean field is closer (in the sense of the norm • 2 L 2 (Ω) ) to the reference. Besides, those uncertain small-scale structures, forgotten by the mean field, are visible in the variance. The variance becomes significant after 10 days of advection, near the stretched saddle points. The Spread-error consistency in the Fourier domain S p e c t r u m of t h e e r r or s an d i t s e st i mat i on )) of the low-resolution SQG model with random initial conditions (red) and of the SQG MU model at the same resolution (blue), at (from top to bottom) t = 10, 13 and 15 days of advection. The reference is the usual SQG model at resolution 512 2 -adequately filtered and subsampled. The low-resolution deterministic model with random initial conditions underestimates the error by at least one order of magnitude whereas our estimation is very precise except at small scales (Colour online). S p e c t r u m of t h e e r r or s an d i t s e st i mat i on )) of the low-resolution SQG model with random initial conditions (red) and of the SQG MU model at the same resolution (blue), at (from top to bottom) t = 17, 20 and 25 days of advection. The reference is the usual SQG model at resolution 512 2 -adequately filtered and subsampled. The low-resolution deterministic model with random initial conditions underestimates the error by at least one order of magnitude whereas our estimation is very precise except at small scales (Colour online).

strong tracer gradients create strong multiplicative noises. Indeed, strong large-scale gradients involve smaller scales, and thus interact with the small-scale velocity σ Ḃ. Then, at t = 17 days, the filament instabilities are triggered by the unresolved velocity stretching effects. The appearance of "pearl necklaces" and the underlying motions of those small-scale eddies are mainly determined by the action of the unresolved velocity component. In consequence, these structures are associated with a high uncertainty in their shapes and locations. Hence, they appear naturally on the variance field. At t = 20, those sources of variance remain and mushroom-like structures also develop near (x, y) = (0, 100), (500, 100), (0, 900) and (500, 900) (in km). The evolution of these fronts are uncertain, and also show up in the variance field. On the day 30 th , these random structures are transported by the zonal jets which are located at y = 0 and y = 500 km.

The empirical moments of order 3 and 4 can also be evaluated with the ensemble. A high 4 th order moment directly relates to the occurrence of extreme events, which is very relevant for dynamical analysis. The point-wise 4-th order moment is centered and normalized to obtain the so-called kurtosis:

m 4 = E (b -E(b)) 4 E (b -E(b)) 2 2 . ( 35
)
The excess kurtosis, m 4 -3 highlights deviations from Gaussianity. In particular, positive values figure the existence of fat-tail distribution. On the right column of Figure 11, the logarithm of the excess kurtosis is displayed for several distinct times. Negative values of the excess kurtosis (which indicates a flatter peak around the mean) have been set to zero. The "pearl necklaces", identified in the variance plots, engender fat-tailed distribution at days t = 17 and 20. The small eddies of a "pearl necklace" have similar vorticity and are close to each other, creating high shears between them. A given eddy can be ejected from the necklace by its closest neighbors, and led up to the north or south down. In such a case, the eddy reaches a zone of the space, neither warm nor cold, with weak variability (e.g. with both local mean and variance being low compared to eddy's temperature). This brings extreme tracer values in statistical homogeneous areas. Finally, the random structures, associated with extreme events are trapped in the zonal jets.

The point-wise moment of order 3 marks the asymmetry of the point-wise tracer distribution. The skewness is the third-order moment of the centered and normalized tracer:

m 3 = E (b -E(b)) 3 E (b -E(b)) 2 3/2 . ( 36
)
Considering the interpretation of excess-kurtosis, the skewness identifies the predominant occurrence of cold (resp. warm) extreme events, associated with the cold (resp. warm) "pearl-necklaces".

Conclusion

Models under location uncertainty involve a velocity partially time-uncorrelated. Accordingly, the material derivative, the interpretation of conservation laws, and the usual fluid dynamics models are modified. In this paper, the random Boussinesq model is approximated by the so-called QG equations. In our random framework, the approximation depends on sub-grid terms scaling. With moderate turbulent dissipation, the PV is randomly transported in the fluid interior up to three source/sink terms. Two of them are smooth in time and cancel out for homogeneous turbulence.

The last forcing term -a random enstrophy source -is related to the angle between stable directions of resolved and unresolved velocities. Similarly to the deterministic case, a uniform PV yields a randomized SQG model, called SQG M U , where the buoyancy is transported in the stochastic sense. Simulation results are considered for the SQG M U model which is a good representation of the transport under location uncertainty. As such, results are believed to hold for any fluid dynamics models under location uncertainty. As found, SQG M U better resolves small-scale tracer structures than a usual SQG model simulated at the same resolution. The prescribed balance between noise and diffusion has also been confirmed. As further highlighted, an ensemble of simulations was able to estimate the amplitude and the position of its own errors in both spatial and spectral domains. This result suggests that the proposed randomized dynamics should be well suited for filtering and other data assimilation methods. On the contrary, a deterministic model with randomized initial conditions, either creates strong errors in its realizations (one order of magnitude larger than the unperturbed deterministic dynamics), or underestimates its own errors (one order of magnitude too low). A MATLAB R code simulating the SQG M U model is available online (http://vressegu.github.io/sqgmu).

As a discussion, we can address the problem of uncertainty quantification (UQ) of an unresolved dynamics from an opposite point of view as the usual setting. Instead of specifying a form for the sub-grid velocity, we can wonder what is the optimal form of SPDE for UQ in fluid dynamics. As demonstrated, randomization of initial conditions is far from being sufficient to quantify uncertainty. Therefore, a random forcing is needed to inject randomness at each time step. The simplest choice is a forcing uncorrelated in time. Otherwise, additional stochastic equations need to be simulated to sample a time-correlated process. This is not desirable in high dimension and the correlation time of the process is often small anyway [START_REF] Berner | Model uncertainty in a mesoscale ensemble prediction system: Stochastic versus multiphysics representations[END_REF]. A forcing uncorrelated in time is a source of energy. So, to be physically acceptable, the SPDE should involve a dissipative term to exactly compensate this source, even in non-stationary regime. The simplest choices of dissipation are diffusion and linear drag. For small-scale processes, the first is more suitable. Now, what is the form of a noise which brings as much energy as a diffusion removes? The proposed approach constitutes a suitable solution toward this goal.

To further improve the accuracy of the UQ, spatial inhomogeneity of the variance tensor a can be introduced from data or from additional models, as discussed in Resseguier et al. (2017a). This inhomogeneity may reduce possible spurious oscillations of tracer stable isolines. Such oscillations are visible on Figure 3 on the sides of the largest vortices. The assumption of time decorrelation may also be a limitation. Nevertheless, as shown by the numerical simulations, the method already achieve very good outcomes with an homogeneous noise component and no real time-scale separation between the resolved and unresolved velocities. Note in particular that since the noise is multiplicative, the random forcing is inhomogeneous even for homogeneous small-scale velocity. Resseguier et al. (2017b) focuses on a system with a clear time-scale separation between the meso and sub-meso scale dynamics to explore the consequences of the QG assumptions under a strong uncertainty assumption (Υ ∼ Ro). A zero PV directly appears in the fluid interior and the horizontal velocity becomes divergent. This divergence provides a simple diagnosis of the frontolysis on warm sides of fronts and frontogenesis on cold sides of fronts.

Future works shall also focus on the potential benefits of the stochastic transport for data assimilation issues. As foreseen, the proposed stochastic formalism opens new horizons for ensemble forecasting techniques and other UQ based dynamical approaches (e.g. [START_REF] Ubelmann | Dynamic interpolation of sea surface height and potential applications for future high-resolution altimetry mapping[END_REF]. This stochastic setup has also been used to characterize chaotic transitions associated with breaking symmetries, also demonstrating interesting perspectives in that context. large-scale for the 0-th order terms. The buoyancy variable does not involve any white noise term, and the small-scale random velocity is thus almost constant along z, as

∂ z u 0 = ∇ ⊥ b 0 and ∂ z (σdB t ) H = O RoD 2 . ( 54
)
Accordingly the variance tensor scales as:

∀i, j ∈ H, ∂ z a ij = O Ro 2 D 4 , (55)
which is negligible in all equations, and the uncertain random field solely depends on the horizontal coordinates. Since Ro/Bu ∼ Ro, the 1-st order term of the buoyancy equation must be kept to describe the evolution of b 0 :

1 Bu D H 0t b 0 + w 1 dt - 1 Υ (∇ • a Hz ) T dt + 1 Υ 1/2 (σdB t ) z = 0, (56)
where, for all functions h,

D H 0t h = d t h + ∇h • u 0 dt + 1 Υ 1/2 (σdB t ) H - 1 2Υ i,j∈H ∂ 2 ij (a ij h) dt. ( 57
)
Taking the derivative along z leads to:

1 Bu   D H 0t ∂ z b 0 + ∇b 0 • ∂ z u 0 dt + 1 Υ 1/2 (σdB t ) H - 1 2Υ i,j∈H ∂ 2 ij (∂ z a ij b 0 ) dt   + ∂ z w 1 dt - 1 Υ (∇ • ∂ z a Hz ) T dt + 1 Υ 1/2 ∂ z (σdB t ) z = 0. (58)
The introduction of the thermal wind equations (54) and incompressibility conditions (47-49) helps simplifying this equation as:

1 Bu D H 0t ∂ z b 0 -∇ • u 1 dt + Ro Bu -1 1 Υ ∇ • (∇ • a H ) T dt - Ro Bu -1 1 Υ 1/2 ∇ • (σdB t ) H = 0. ( 59
)
Note the factor (Ro/Bu) -1 appears. It comes from the incompressible conditions ( 48) and ( 49), leading ∇ • (σdB t ) H and ∇ • (∇ • a H )

T dt to both scale as Ro/Bu. The hydrostatic balance at 0-order links the buoyancy to the pressure, and then to the stream function

∂ z b 0 = ∂ 2 z φ 0 = ∂ 2 z ψ 0 . (60)
The 1-st order term of the vertical velocity is not known. Yet, the system can be closed using the vorticity equation at order 1:

∇ ⊥ • D H 0t u 0 + ∇•u 1 + Ro Bu -1 ∇•(σdB t ) H + ∇(βy) • (u 0 dt + (σdB t ) H ) = 0, ( 61
)
where the divergence terms come from the constant Coriolis term.

Again, factors (Ro/Bu) -1 compensate the order of magnitude of ∇•(σdB t ) H and ∇•(∇ • a H ) T dt. Then,

D H 0t (∆ψ 0 ) + ∇•u 1 dt + 1 Υ 1/2 Ro Bu -1 ∇•(σdB t ) H + β v 0 dt + 1 Υ 1/2 (σdB t ) y + 1 Υ 1/2 tr ∇ ⊥ (σdB t ) T H ∇u T 0 - 1 2Υ i,j∈H ∂ 2 ij ∇ ⊥ a ij • u 0 dt = 0. ( 62
)
To make appear the transport of PV, we note that:

D H 0t (1 + βy) = -β v 0 dt + 1 Υ 1/2 (σdB t ) y + ∇•a Hy βdt - Ro Bu -1 1 2Υ ∇ • (∇ • a H ) T dt. ( 63
)
Then, using ( 59), ( 60) and ( 62), we get:

D H 0t ∆ψ 0 + 1 + βy + 1 Bu ∂ 2 z ψ 0 = -∇•a Hy βdt - Ro Bu -1 1 2Υ ∇ • (∇ • a H ) T dt -tr ∇ ⊥ (σdB t ) T H ∇u T 0 + 1 2Υ i,j∈H ∂ 2 ij ∇ ⊥ a ij • u 0 dt. ( 64
)
We recall that coefficients (Ro/Bu) -1 are still present since

∇•(σdB t ) H ∼ ∇ • (∇ • a H ) T dt ∼ Ro Bu . (65)
If we rewrite the equation with dimensional quantities, the evolution equation for u 0 = lim Ro→0 u is obtained (dropping the index 0 for clarity):

D H t Q = -tr ∇ ⊥ (σdB t ) T H ∇u T + 1 2 i,j∈H ∂ 2 ij ∇ ⊥ a ij • u dt -1 2 ∇ • (∇ • (a H f )) T dt, ( 66
)
where Q is the QG potential vorticity:

Q = ∆ψ + f + 1 N 2 ∂ 2 z ψ. (67)
Note, (51) provides the geostrophic balance for the small-scale velocity component. To express the material derivative of Q, the noise term is expanded:

-tr ∇ ⊥ (σdB t ) T H ∇u T = - k,j∈H ∂ 2 kj ψ∂ k σ j• dB t . ( 68
)
According to Resseguier et al. (2017a), the difference between the material derivative, D t Q, and the stochastic transport operator D t Q, is a function of the time-uncorrelated forcing:

D t Q = f 1 dt + h T 1 dB t , D t Q = f 2 dt + h T 2 dB t , ⇐⇒ f 2 = f 1 + tr (σ T ∇) h T 1 , h 1 = h 2 . ( 69
)
The expression of h 1 is given by equation (66) and the above formulas give:

D t Q -D t Q = i∈H σ i• ∂ i   - j,k∈H ∂ k σ j• ∂ 2 kj ψ   T , (70)
=i,j,k∈H

σ i• ∂ 2 ik σ T j• ∂ 2 kj ψ + σ i• ∂ k σ T j• ∂ 3 ijk ψ . ( 71
)
With the use of the small-scale incompressibility, we obtain:

1 2 i,j∈H ∂ 2 ij ∇ ⊥ a ij • u = i,j,k∈H ∂ j σ i• ∂ 2 ik σ T j• ∂ k ψ + ∂ j σ i• ∂ k σ T j• ∂ 2 ik ψ + σ i• ∂ 2 ik σ T j• ∂ 2 jk ψ + σ i• ∂ k σ T j• ∂ 3 ijk ψ . ( 72
)
From ( 71) and ( 72), it yields:

D t Q -   D t Q -1 2 i,j∈H ∂ 2 ij ∇ ⊥ a ij • u   = i,j,k∈H ∂ j σ i• ∂ 2 ik σ T j• ∂ k ψ + ∂ j σ i• ∂ k σ T j• ∂ 2 ik ψ . (73)
Denoting, α, the following matrix

α ij = k∈H ∂ k σ i• ∂ j σ T k• = k∈H ∂ k (σ i• ∂ j σ T k• ), (74)
we have ∇ • (α∇ψ) = i,j,k∈H

∂ j σ i• ∂ 2 ik σ T j• ∂ k ψ + ∂ j σ i• ∂ k σ T j• ∂ 2 ik ψ , (75)
= D t Q -   D t Q -1 2 i,j∈H ∂ 2 ij ∇ ⊥ a ij • u   , (76)
and the material derivative of the PV finally reads:

D H t Q = ∇ • (α∇ψ)dt -1 2 ∇ • (∇ • (a H f )) T dt -tr ∇ ⊥ (σdB t ) T H ∇u T . ( 77
)
To note, the transpose of the matrix α has a compact expression:

α T = p ∇σ T Hp 2 . ( 78
)
To better assess the role of the random source term (the last term of (77)), it is decomposed in terms of symmetric and anti-symmetric parts of the small-scale/large-scale deformation tensors. Let us denote S and S σdBt the symmetric parts of ∇u T and ∇(σdB t ) T H , respectively. Associated with divergence-free velocities, these symmetric parts, so-called strain rate tensors, have zero trace.

34

Terms -1 2 ωJ and -1 2 ω σdBt J will stand for the anti-symmetric parts, where J = 0 -1 1 0 is the π/2 rotation. The factors ω and ω σdBt are the large-scale and the small-scale components of the vorticity, respectively. Using J J = -I d and tr[M N ] = tr[N M ] yields:

-tr ∇ ⊥ (σdB t ) T H ∇u T = -tr J S σdBt -1 2 ω σdBt J S -1 2 ωJ ,

This term thus only depends on the strain rate tensors of u and (σdB t ) H . The PV transport can thus be rewritten as:

D H t Q = ∇ • (α∇ψ)dt -1 2 ∇ • (∇ • (a H f )) T dt -tr [SJ S σdBt ] . (82)
The noise term can be further expressed using the stable directions of the flows defined by u and (σdB t ) H , respectively. In the following, we will omit writing the dB t factor. The two strain rate tensors are decomposed in orthogonal basis:

S = V ΞV T = 2 p=1 Ξ pp V •p V T •p and S σdBt = W ΛW T , (83)
where V T •p V •q = W T •p W •q = δ pq , Ξ 11 = -Ξ 22 < 0 and Λ 11 = -Λ 22 < 0.

-tr [SJ S σdBt ] = -2 p,q=1

Ξ pp Λ qq tr V •p V T •p J W •q W T •q , (84) = - 2 p,q=1 Ξ pp Λ qq V T •p W •q V T •p J W •q , (85)
where θ pq = (V •p , W •q ) is the angle between V •p and W •q . Using the relations between the eigenvalues and the orthogonality of the eigenvectors, it finally comes:

-tr [SJ S σdBt ] =

Figure 1 :

 1 Figure 1: Buoyancy (m.s -2 ) at t = 0, 5, 10, 13, 15, 16, 20, 30 days of advection for the usual SQG model at resolution 512 2 (Colour online).

  Figure2: Spectrum (m 2 .s -2 /(rad.m -1 )) of w, at the initial time, in blue, spectrum of σ H Ḃ (up to a multiplicative constant), in red, and slope -5/3 in black. In the simulation performed, σ H Ḃ is restricted to a narrow spectral band. Thus, this velocity component almost only acts near the resolution cutoff, where the large-scale component, w, has a low energy (Colour online).

Figure 3 :

 3 Figure3: Buoyancy (m.s -2 ) and its spectrum (m 2 .s -4 /(rad.m -1 )) at the 17 th day of advection for SQG MU at resolution 128 2 (top), SQG at resolution 512 2 (middle) and at resolution 128 2 (bottom). Unlike SQG MU , the low-resolved SQG simulation diffuses the "pearl necklaces", noticeable only at higher resolution (Colour online).

Figure 4 :

 4 Figure4: Buoyancy (m.s -2 ) and its spectrum (m 2 .s -4 /(rad.m -1 )) at the t = 17 th day of advection for the SQG MU model with a small-scale velocity component three times weaker than the one prescribed by the diffusion coefficient a H (top), with the correct amount of small-scale energy (middle) and a small-scale velocity three times higher than the model diffusion. If the prescribed balance between noise and diffusion is not met the tracer field becomes quickly too smooth or too noisy (Colour online).

Figure 5 :Figure 6 :

 56 Figure 5: Buoyancy bias absolute value, e = | E{b} -b ref |, (m.s -2 ) of the SQG MU model (left) and its estimation, est , (1.96× the standard deviation of the ensemble) (right) at resolution 128 2 at (from top to bottom) t = 10, 13 and 15 days of advection. The reference is the usual SQG model at resolution 512 2 -adequately filtered and subsampled (Colour online). 16

  d . m -1S p e c t r u m of t h e e r r or s an d i t s e st i mat i on d . m -1 S p e c t r u m of t h e e r r or s an d i t s e st i mat i on

Figure 8 :

 8 Figure 8: Spectrum of the buoyancy bias absolute value, e = | E{b} -b ref |, (lines with crosses) and spectrum of the estimated error, est , (1.96× the standard deviation of the ensemble) (lines with circles) (m 2 .s -4 /(rad.m -1)) of the low-resolution SQG model with random initial conditions (red) and of the SQG MU model at the same resolution (blue), at (from top to bottom) t = 10, 13 and 15 days of advection. The reference is the usual SQG model at resolution 512 2 -adequately filtered and subsampled. The low-resolution deterministic model with random initial conditions underestimates the error by at least one order of magnitude whereas our estimation is very precise except at small scales (Colour online).

  d . m -1S p e c t r u m of t h e e r r or s an d i t s e st i mat i on d . m -1 S p e c t r u m of t h e e r r or s an d i t s e st i mat i on

Figure 9 :

 9 Figure 9: Spectrum of the buoyancy bias absolute value, e = | E{b} -b ref |, (lines with crosses) and spectrum of the estimated error, est , (1.96× the standard deviation of the ensemble) (lines with circles) (m 2 .s -4 /(rad.m -1)) of the low-resolution SQG model with random initial conditions (red) and of the SQG MU model at the same resolution (blue), at (from top to bottom) t = 17, 20 and 25 days of advection. The reference is the usual SQG model at resolution 512 2 -adequately filtered and subsampled. The low-resolution deterministic model with random initial conditions underestimates the error by at least one order of magnitude whereas our estimation is very precise except at small scales (Colour online).

1Figure 10 :

 10 Figure10: Point-wise mean (left) and standard deviation (right) of the buoyancy (m.s -2 ) at t = 17, and 30 days of advection for SQG M U model at resolution 128 2 . The moments are computed through MCMC simulations. The mean field is smoother than the individual realizations. Areas of higher variance appear first near the stretched saddle points. Then, at t = 17 days, the filament instabilities are triggered by the unresolved velocity component. The appearance of "pearl necklaces" can be observed. At t = 20, mushrooms-like structures also develop in the variance field near (x, y) = (0, 100), (500, 100), (0, 900) and (500, 900) (in km). At t = 30 days, these random structures are transported by the zonal jets (Colour online).

3Figure 11 :

 11 Figure 11: Point-wise skewness, and logarithm of the excess kurtosis of the buoyancy at t = 17, 20 and 30 days of advection for SQG MU model at resolution 128 2 . The moments are computed through MCMC simulations. Negative excess kurtosis is set to 0. The point-wise law of the tracer is clearly non-Gaussian. The "pearl necklace" events identified in the variance plots leads to fattailed distributions with skewness at t = 17 and 20 days. The random structures, associated with fat tails are then trapped in the zonal jets (Colour online).

=

  -tr [SJ S σdBt ] -1 2 ω σdBt tr [S] [SJ S σdBt ] .

ΞΞ

  pp Λ qq cos(θ pq ) cos θ pq + pp Λ qq sin(2θ pq ),

  1 2 Ξ 11 Λ 11 sin (2θ 11 ) -sin 2 θ 11 -

		π 2	-sin 2 θ 11 +	π 2	+ sin (2θ 11 ) ,
	= 2 Ξ 11 Λ 11	sin (2θ 11 ) .			(88)
	>0

Formally each coefficient of (t → Bt) is a cylindrical I d -Wiener process (see Da Prato and Zabczyk (1992) and[START_REF] Prévôt | A concise course on stochastic partial differential equations[END_REF] for more information on infinite dimensional Wiener processes and cylindrical I d -Wiener processes).

All the random simulations are performed with 200 -128 2 mesh-size -realizations.

Note this simulation is afterward spatially filtered and subsampled to the same resolution as the ensemble

Acknowledgments

The authors thank Guillaume Lapeyre, Aurélien Ponte, Jeroen Molemaker, Guillaume Roulet and Jonathan Gula for helpful discussions. We also acknowledge the support of the ESA DUE GlobCurrent project (contract no. 4000109513/13/I-LG), the "Laboratoires d'Excellence" Comin-Labs, Lebesgue and Mer (grant ANR-10-LABX-19-01) through the SEACS project.

A Non-dimensional Boussinesq equations

To derive a non-dimensional version of the Boussinesq equations under location uncertainty (Resseguier et al., 2017a), we scale the horizontal coordinates xh = Lx h , the vertical coordinate z = hz, the aspect ratio D = h/L between the vertical and horizontal length scales. A characteristic time t = T t corresponds to the horizontal advection time U/L with horizontal velocity ũ = U u. A vertical velocity w = (h/L)U w is deduced from the divergence-free condition. We further take a scaled buoyancy b = Bb, pressure φ = Φφ (with the density scaled pressures φ = p /ρ b and d t φ σ = d t p σ /ρ b ), and the earth rotation f * = f k. For the uncertainty variables, we consider a horizontal uncertainty ãH = A u a H corresponding to the horizontal 2 × 2 variance tensor; a vertical uncertainty vector ãzz = A w a zz and a horizontal-vertical uncertainty vector ãHz = √ A u A w a Hz related to the variance between the vertical and horizontal velocity components. The resulting non-dimensional Boussinesq system under location uncertainty becomes:

Nondimensional Boussinesq equations under location uncertainty

Momentum equations

Buoyancy equation

Effective drift

Here, we do not separate the time-correlated components and the time-uncorrelated components in the momentum equations. The terms in O (Ro/Bu) and O (Ro/Bu) 2 are related to the time-uncorrelated vertical velocity. These terms are too small to appear in the final QG model (Bu = O (1) in QG approximation) and not explicitly shown. We only make appear the big O approximations. Traditional non-dimensional numbers are introduced : the Rossby number R o = U/(f 0 L) with f 0 the average Coriolis frequency; the Froude number (Fr = U/(N h)), ratio between the advective time to the buoyancy time; Eu, the Euler number, ratio between the pressure force and the inertial forces, Γ = Bh/U 2 = D 2 BT /W the ratio between the mean potential energy to the mean kinetic energy. To scale the buoyancy equation, the ratio between the buoyancy advection and the stratification term has also been introduced:

Besides those traditional dimensionless numbers, this system introduces Υ , relating the largescale kinetic energy to the energy dissipated by the unresolved component:

B QG model under moderate uncertainty

Hereafter, we consider the QG approximation (Ro 1 and Bu ∼ 1), for Υ ∼ 1. We focus on solutions of the Boussinesq model with Rossby number going to zero. To derive the evolution equations corresponding to this limit, the solution of the non-dimentional Boussinesq model (Appendix A) is developed as a power series of the Rossby number:

According to the horizontal momentum equation (37a), the scaling of the pressure still corresponds to the usual geostrophic balance. This sets the Euler number as:

For the ocean, the aspect ratio, D = H/L, is small and D 2 1. As a consequence,

Therefore, the inertial and diffusion terms are negligible in the vertical momentum equation. The hydrostatic assumption is still valid. This leads to the classical QG scaling of the buoyancy equation:

In the following, the subscript H is omitted for the differential operators Del, ∇, and Laplacian, ∆. They all represent 2D operators. Only keeping terms of order 0 and 1, we get the following system:

The thermodynamic equation ( 46) at 0 order leads to :

and then, by the large-scale incompressibility equation ( 47), the 0-order horizontal velocity is divergence-free. Following the scaling assumption, the horizontal small-scale velocity is also divergencefree (48). The horizontal momentum equation ( 44) at the 0-th order leads to:

where time-correlated and time-uncorrelated components have been separated by the mean of uniqueness of the semi-martingale decomposition [START_REF] Kunita | Stochastic flows and stochastic differential equations[END_REF]. Being divergent-free, both components can be expressed with two stream functions ψ 0 and d t ψ σ :

exactly corresponding to the dimensionless pressure terms:

Deriving these equations along z and introducing the hydrostatic equilibrium (45) -decomposed between correlated and uncorrelated components -yields the classical thermal wind balance at