V Resseguier
email: valentin.resseguier@inria.fr

E Mémin

B Chapron

Geophysical flows under location uncertainty, Part III SQG and frontal dynamics under strong turbulence conditions

Keywords: stochastic subgrid tensor, uncertainty quantification, upper ocean dynamics

published or not. The documents may come

Geophysical flows under location uncertainty, Part III SQG and frontal dynamics under strong turbulence conditions

Introduction

Quasi-Geostrophic (QG) models are standard models to study mesoscale barotropic and baroclinic dynamics. Assuming uniform Potential Vorticity (PV) in the fluid interior, the Surface Quasi-Geostrophic (SQG) model helps describe the surface dynamics [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]Lapeyre and Klein, 2006;[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]. Despite its simplicity, the SQG relation provides a good diagnosis to relate mesoscale surface buoyancy fields to surface and interior velocity fields. Nevertheless, QG and SQG paradigms assume strong rotation and strong stratification (Fr ∼ Ro 1) and thus neglect the submesoscale ageostrophic dynamics. In particular, the QG velocity is horizontal and solenoidal. This structure prevents the emergence and development of realistic submesoscale features such as frontogenesis, restratification, and asymmetry between cyclones and anticyclones (Lapeyre et al., 2006;[START_REF] Klein | Upper ocean turbulence from high-resolution 3D simulations[END_REF]. In contrast, the QG +1 [START_REF] Muraki | The next-order corrections to quasigeostrophic theory[END_REF] and SQG +1 [START_REF] Hakim | A new surface model for cyclone-anticyclone asymmetry[END_REF]) models capture such phenomenon with a (one degree) higher order power series expansions in the Rossby number. This comes with an 2 Models under location uncertainty Hereafter, we briefly outline the main ideas for the derivation of these stochastic models (for a more complete description, see Resseguier et al. (2017a)). This relies on a decomposition of the flow velocity in terms of a large-scale component, w, and a random field uncorrelated in time, σ Ḃ:

dX dt = w + σ Ḃ. ( 1
)
The latter represents the small-scale velocity component. This solenoidal, possibly anisotropic and non-homogeneous random field corresponds to the aliasing effect of the unresolved velocity component. To parametrize its spatial correlations, an infinite-dimensional linear operator, σ, is applied to a space-time white noise, Ḃ. The decomposition (1) leads to a stochastic representation of the Reynolds transport theorem (RTT) and of the material derivative, D t (derivative along the flow (1)). In most cases, this derivative coincides with the stochastic transport operator, D t , defined for every field, Θ, as follow:

D t Θ = d t Θ = Θ(x,t+dt)-Θ(x,t) Time increment + (w dt + σdB t ) • ∇Θ Advection -∇ • 1 2 a∇Θ Diffusion dt, (2)
where the time increment term d t Θ stands instead of the partial time derivative as Θ is 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 , (3)
and the modified drift is given by

w = w -1 2 (∇ • a) T . (4)
For a divergent small-scale velocity, this drift would involve an additional component (see equation ( 4) of Resseguier et al. (2017a)). The stochastic RTT and material derivative involve a diffusive subgrid term, a multiplicative noise and a modified advection drift induced by the small-scale inhomogeneity. This material derivative has a remarkable conservative property. Indeed, for any field, Θ, randomly transported, i.e.

Θ(X(t + ∆t), t + ∆t) = Θ(X(t), t), (5)
Resseguier et al. (2017a) showed that the energy of each realization is conserved:

d dt Ω Θ 2 = 0. ( 6
)
The RTT enables us to express the conservation law of mechanics (linear momentum, energy, mass) with a partially known velocity. Deterministic and random subgrid parametrizations for various geophysical flow dynamics can then directly be obtained. Stochastic Navier-Stokes and Boussinesq models can be derived as discussed by [START_REF] Mémin | Fluid flow dynamics under location uncertainty[END_REF] and Resseguier et al. (2017a). The latter model involves random transports of buoyancy and velocity, together with incompressibility constraints.

Mesoscale flows under strong uncertainty

From the Boussinesq model, the QG assumptions state a strong rotation and a strong stratification. This is of particular interest to study flows at mesoscale, where both kinetic and buoyant dynamics are important. More specifically, we focus on horizontal length scales, L, such as:

1 Bu = Fr Ro 2 = L L d 2 ∼ 1 and 1 Ro = Lf 0 U 1, ( 7
)
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. In the following, both differential operators Del, ∇, and Laplacian, ∆, represent 2D operators. Moreover, σ H• stands for the horizontal component of σ, a H for σ H• σ T H• and A u for its scaling.

Specific scaling assumptions

Similarly to Resseguier et al. (2017b), scalings within the QG framework (7) can authorize the set up of a non-dimensional stochastic Boussinesq model amenable to further simplifications.

Quadratic variation scaling

Models under location uncertainty involve subgrid terms which have also to be scaled. A new dimentionless number, Υ , quantifying the ratio of horizontal advection and horizontal turbulent diffusion is therefore introduced:

Υ = U/L A u /L 2 = U 2 A u /T . ( 8
)
We can also relate it to the ratio of Mean Kinetic Energy (MKE), U 2 , to the Turbulent Kinetic Energy (TKE), A u /T σ , where T σ is the small-scale correlation time. This reads:

Υ = 1 M KE T KE , (9)
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 [START_REF] Majda | Models for stochastic climate prediction[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]. The number Υ/Ro measures the ratio between rotation and horizontal diffusion. For a parameter Υ close or larger than unity, the geostrophic balance still holds (Resseguier et al., 2017b), whereas for Υ ∼ Ro, this balance is modified. Throughout this paper, we focus on this specific scaling. The parameter Υ depends through A u on the flow and on the resolution scale. In order to specify the scaling and the resulting associated model, knowledge of the characteristic horizontal eddy diffusivity or eddy viscosity is needed. Tuning experiences of usual subgrid parametrizations may provide such information, and [START_REF] Boccaletti | Mixed layer instabilities and restratification[END_REF] give some examples of canonical values.

If absence of characteristic values, 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. Small values of Υ are generally relevant for the ocean where the TKE is often one order of magnitude larger than the MKE (Wyrtki et al., 1976;[START_REF] Richardson | Eddy kinetic energy in the north atlantic from surface drifters[END_REF][START_REF] Stammer | Global characteristics of ocean variability estimated from regional topex/poseidon altimeter measurements[END_REF]Vallis, 2006). Note that here the TKE may encompass all the unresolved dynamics down to the Kolmogorov scale.

Vertical unresolved velocity

To scale the vertical unresolved velocity, we consider

(σdB t ) z (σdB t ) H ∼ Ro Bu D, (10)
where D = h/L is the aspect ratio and the subscript H indicates horizontal coordinates. The ω-equation [START_REF] Giordani | Advanced insights into sources of vertical velocity in the ocean[END_REF] justifies such a scaling. 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 , (11)
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. ( 12
)
The Burger number is small at planetary scales where the rotation dominates (W/U ∼ DRo) and is large at submesoscales where the stratification dominates (W/U ∼ DRo/Bu). For the small-scale velocity σ Ḃ, the latter is thus more relevant.

Relations between the isopicnical tilt and mixing give another justification of the scaling (10). Based on baroclinic instabilities theory, anisotropy specifications of eddy diffusivity sometimes rely on this tilt (Vallis, 2006). Moreover, several other authors suggest that the eddy activity and the associated mixing mainly occur along isentropic surfaces [START_REF] Gent | Isopycnal mixing in ocean circulation models[END_REF][START_REF] Pierrehumbert | Global chaotic mixing on isentropic surfaces[END_REF].

For QG dynamics, the Burger number is of order one and the scaling in DRo and in DRo/Bu coincides. In particular, they encode a mainly horizontal unresolved velocity:

(σdB t ) z (σdB t ) H ∼ Ro Bu D D. ( 13
)
This is consistent with the assumption of a large stratification, i.e. flat isopycnicals, if we admit that the activity of eddies preferentially appears 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

The beta effect is weak at mid-latitude mesoscales. Yet, at the first order, it influences the absolute vorticity. So, we choose the same scaling as Vallis ( 2006):

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

Stratified Quasi-Geostrophic model under strong uncertainty

Strong uncertainty condition corresponds to Υ having an order of magnitude close to the Rossby number. More specifically, we assume Ro Υ 1. In this situation, the random eddies have larger energy than the large-scale mean kinetic energy. Accordingly, the diffusion and drift terms are one order of magnitude larger than the advection terms.

In the case of strong ratio Υ , the diffusion is very large and the system is not approximately in geostrophic balance anymore. The large-scale horizontal velocity becomes divergent, and decoupling the system is more tedious. For sake of simplicity, in the following we consider the case of homogeneous and horizontally isotropic turbulence. As a consequence, the variance tensor, a, is constant in space and diagonal:

a =   a H 0 0 0 a H 0 0 0 a z   . ( 15
)

Modified geostrophic balance under strong uncertainty

For horizontal homogeneous turbulence, the large-scale geostrophic balance is modified by the horizontal diffusion, whereas the unresolved velocity is in geostrophic balance:

f × u - a H 2 ∆u = - 1 ρ b ∇p , (16a) f × σ H dB t = - 1 ρ b ∇d t p σ , ( 16b
)
where u is the large-scale horizontal velocity, p the time-correlated component of the pressure, ṗσ = dp σ /dt the time-uncorrelated component, and ρ b is the mean density. For a constant Coriolis frequency, the first equation can be solved in Fourier space. The Helmholtz decomposition of the velocity reads:

u = ∇ ⊥ ψ + ∇ ψ, (17a)
ψ = 1 + k kc 4 2 -1 p ρ b f , (17b) ψ = 1 k 2 c ∆ψ, (17c)
where k c = 2f 0 /a H and the hat accent indicates a horizontal Fourier transform. This solution is derived in Appendix A using geometric power series of matrices. The obtained formula is valid for any right-hand side in equation (16a). For instance, additional forcing such as an Ekman stress could be taken into account. In equation ( 17), the solenoidal component of the velocity, ∇ ⊥ ψ, corresponds to the usual geostrophic velocity multiplied by a low-pass filter (17b). The irrotational (ageostrophic) component of the velocity, ∇ ψ, dilates the anticyclones (maximum of pressure and negative vorticity) and shrinks the cyclones (minimum of pressure and positive vorticity) at small scales. Indeed, according to equation (17c), the divergence of the velocity corresponds to the vorticity Laplacian divided by k 2 c . Naturally, this structure is reminiscent of the Ekman model where divergence and vorticity would be related by a double vertical derivative:

δ = 2 E k ∂ 2 z ζ where δ = ∇•u, ζ = ∇ ⊥ • u, (18)
and E k is the thickness of the Ekman layer. The turbulent diffusion involved in equation ( 17c) is rather horizontal due to the strong stratification assumption (see ( 10)). In the proposed stochastic model, the divergent component and the low-pass filter of the system (17) are parameterized by the spatial cutoff frequency k c , which moves toward larger scales when the diffusion coefficient a H increases. If both the vorticity and the divergence can be measured at large scales, the previous relation should enable to estimate the cutoff frequency k c by fitting terms of equation (17c). Then, the horizontal diffusive coefficient, a H , or the variance of the horizontal small-scale velocity (at the time scale ∆t), a H /∆t, can be deduced.

Modified SQG relation under strong uncertainty

To derive a QG model, we use the other equations of the stochastic Boussinesq model at the 0-order. After some algebra (see Appendix C), we obtain directly a zero PV in the fluid interior:

P V = ∆ + 1 + ∆ 2 k 4 c ∂ z f 0 N 2 ∂ z ψ = 0, ( 19
)
where k c = 2f 0 /a H . The assumptions used here correspond to the same used for a classical QG model (Vallis, 2006), except that the dissipation, due to the noise, is strong. It is a striking result.

Instead of finding a model in the form of a classical QG model, developments, through a strong uncertainty, directly leads to the description of surface dynamics, a SQG model. It means that the subgrid dissipation prevents the development of the interior dynamics. Without this dynamics, no baroclinic instabilities can grow (Lapeyre and Klein, 2006). If the stratification is vertically invariant, this static linear equation can be solved by imposing a vanishing condition in the deep ocean (z → -∞) and a specified boundary value at a given depth (z = η). The horizontal Fourier transform of the solution then reads:

ψ(k, z) = ψ(k, η) exp     N k 2 f 0 1 + k kc 4 2 (z -η)     . ( 20
) At z = η, the modified SQG relation is: b(k, η) = N k 2 1 + k k c 4 2 ψ(k, η), (21)
where b stands for the buoyancy. In the following, we will refer to (21) as the SQG relation under Strong Uncertainty (SQG SU ). For low wave number or moderate uncertainty ( k/k c 2 2 ∼ Ro/Υ 1), we retrieve the standard SQG relation. The expression of the stream function as a function of the buoyancy is expressed as the convolution with a Green function, G SQG = x -1 /(2πN ), and the velocity decays rapidly as the inverse of the square distance to the point vortex center. On the other hand, for very high wave numbers or very large uncertainty ( k/k c 2 2 ∼ Ro/Υ 1), the velocity tends very quickly to zero. For strong uncertainty or small scales (

k/k c 2 2 ∼ Ro/Υ ∼ 1), b = √ 2N k c k 2 2 + O k →kc k k c 2 -1 2 ψ. (22)
Accordingly, we may see the SQG SU relation as an intermediary between two relevant models in geophysics: the SQG dynamics where the tracer (the buoyancy) is proportional to k 2 ψ and a two-dimensional flow dynamics where the tracer (the vorticity) is proportional to k 2 2 ψ. In the latter case, the streamfunction can be expressed as the convolution of the buoyancy with the Green function, G 2D = k c /(2 √ 2πN ) ln x , and the velocity decays slowly as the inverse of the distance to the point vortex center. Nevertheless, contrary to the two-dimensional flow and the SQG models, the 2D velocity u is divergent (see equation ( 17c)). The total horizontal velocity can be computed from the buoyancy, through the Helmholtz decomposition (17a), the modified SQG relation ( 21)

Figure 1: Value of the interior buoyancy created by a warm spot and a cold spot at the surface. The two components of the velocity are also shown. The upwelling and shrinking of the cyclone (cold spot) and the downwelling and spreading of the anticyclone (warm spot) are clearly visible (Colour online).

Figure 2: Value of the interior buoyancy created by a front at the surface. The two components of the velocity are also shown. The divergence effects will strengthen the front on the cold side (frontogenesis) and smooth the front on the warm side (frontolysis). The vertical velocity is here much weaker than in the case of isolated spots (Colour online). and the equation (17c). As derived, the vertical velocity is finite and given by the main balance of the buoyancy equation:

w = f 0 N 2 1 k 2 c ∆b. ( 23
)
Note that this equation is not derived from a non-hydrostatic vertical momentum equation. Equation ( 23) is directly obtained from the thermodynamic equation. It expresses the fact that, under strong stratification and strong horizontal diffusion, the buoyancy anomalies are mainly created by vertical advection. This relation is similar to the result of [START_REF] Garrett | Dynamical aspects of shallow sea fronts[END_REF], except the proportionality coefficient. Indeed, [START_REF] Garrett | Dynamical aspects of shallow sea fronts[END_REF] consider vertical diffusion and neglect the horizontal one. Invoking the thermal wind relation and the stratification structure, vertical variations are then associated with horizontal buoyancy variations. In the present development, the vertical velocity scales as Ro/(Υ Bu) DU ∼ k/k c Figures 1 and 2 show the static link between the 3D velocity and buoyancy for two isolated vortices and a front, respectively. As obtained, the solenoidal component is similar to the classic SQG velocity. In Figure 1, the non-rotational component forces the anticyclone (warm spot) to spread, and the cyclone (cold spot) to shrink. Note that our study focuses on the ocean dynamics. For atmospheric applications, the vertical axis should be inverted and the sign of the temperature anomaly changed [START_REF] Ragone | A study of surface semi-geostrophic turbulence: freely decaying dynamics[END_REF]. In Figure 2, the irrotational component is weak on the warm side of the front, but strongly strengthens the cold side. As modeled in the SQG +1 [START_REF] Hakim | A new surface model for cyclone-anticyclone asymmetry[END_REF] and Surface Semi-Geostrophic (SSG) [START_REF] Badin | Surface semi-geostrophic dynamics in the ocean[END_REF][START_REF] Ragone | A study of surface semi-geostrophic turbulence: freely decaying dynamics[END_REF]) models, a frontolysis (resp. frontogenesis) develops on the warm (resp. cold) side of the front. In Figure 1, a downwelling of warm water and a upwelling of cold water appear. As the vertical velocity comes from the thermodynamic equation and not from the vertical momentum equation, it is the cause of the buoyancy anomaly not its consequence. Whereas the irrotational horizontal component is stronger close to a front than within an eddy, the vertical velocity associated with a front is found much weaker than the one associated with an isolated eddy.

Diagnostic under strong uncertainty

As derived, under strong uncertainty, the eddy diffusion is substantial and modifies the geostrophic balance (17). The velocity becomes divergent and equation (17c) offers a diagnostic of this divergence. This diagnostic states that the divergence should be proportional to the Laplacian of the vorticity:

δ = 1 k 2 c ∆ζ where δ = ∇•u, ζ = ∇ ⊥ • u. (24)
To evaluate the relevance of this diagnostic, outputs of a realistic 3D high-resolution oceanic simulation are used. During winter, the eddy activities are usually stronger, especially close to energetic currents. For this reason, the Gulf-Stream during winter season is a test-bed region for highresolution simulation [START_REF] Gula | Topographic vorticity generation, submesoscale instability and vortex street formation in the gulf stream[END_REF].

Figure 3 shows the temperature of the first and of the 58 th day. Simulations are threedimensional and involve a fine spatial and temporal resolutions. Equation ( 17c) is a surface mesoscale diagnostic valid far from the coasts. Consequently, the surface fields are filtered temporally and spatially. The final time step is one day and the final resulting spatial resolution is 3 km. Figure 3 displays the original surface field and the filtered cropped fields.

Figure 4 compares the reference divergence field to our estimate, the Laplacian of the vorticity. An overall agreement clearly emerges. Nonetheless, the small scales of our estimate are more energetic than the small scales of the real divergence field. For this reason, the spatial fields are further filtered at a resolution of 30 km. Except for some small spots, estimation and reference are similar. In particular, fronts -associated with two length scales: one at sub-mesoscales and one at mesoscales -are highlighted.

Figure 5 specifies the relevance and the limitations of the proposed diagnostic. The spectra of the two fields unveil a very good match at mesoscale range (L > 60km i.e. κ < 10 -4 ), whereas they differ at sub-mesoscales. This difference is certainly not surprising, the estimation being derived for large scale components. Note, the velocity divergent component is far from being zero in the mesoscale range. Compared to the solenoidal component, its spectrum is certainly much flatter and smaller in this range. Nevertheless, the mesoscale divergence is stronger than the sub-mesoscales divergence. The ratio of Fourier transform modulus further confirms the accuracy of our diagnostic at mesoscales and makes clear the difference at sub-mesoscales. The -1 slope may suggest that a fractional diffusion would be preferable to a Laplacian diffusion at those scales.

The complementary analysis is the coherence, which is a measure of the phase relationship between two fields. Specifically, the coherence is the Fourier modes correlation coefficient:

  δ(k) ∆ζ(k) δ(k) ∆ζ(k)   , (25)
where denotes the real part. The coherence is the cosinus of the phase shift, θ, between the two fields. Here, we directly show the phase shift averaged on angular spatial frequencies.

For our estimate the phase-shift is about 0.8 ≈ π/4. It means that a linear transformation of the large-scale vorticity can explain more than half of the divergence. As a comparison, the same analysis was done with the SQG relation, using temperature anomaly instead of buoyancy (not shown). The phase shift was similar.

From Figure 4, one further get a rough estimation for the multiplicative constant of the proposed diagnostic:

k 2 c ≈ 10 -7 m -2 . It suggests a spatial cutoff k -1 c
≈ 3 km and a diffusion coefficient a H /2 ≈ 1000 m 2 .s -1 . This value is canonical, according to [START_REF] Boccaletti | Mixed layer instabilities and restratification[END_REF], which upholds the proposed approach. To confirm the validity of our strong uncertain assumption, it can be evaluated:

Ro Υ ∼ k k c 2 2 ∼ a H 2f 0 κ 2 ∼ 0.1, (26)
with 2π/κ = 60 km.

The unresolved energy can also be estimated. From a mesoscale point of view, motions induced with diurnal cycles can be approximated as delta-correlated processes. Hence, an estimation of the unresolved horizontal velocity amplitude shall follow from a H /∆t ≈ 10 -1 m.s -1 , with ∆t = 1 day. Considering the present simulation, this is consistent with the sub-mesoscale velocity field.

Conclusion

To develop models under location uncertainty, the highly-oscillating unresolved velocity component is assumed to be uncorrelated in time. Consequently, the expression of the material derivative and hence most fluid dynamics models are modified, taking into account an inhomogeneous and anisotropic diffusion, an advection correction and a multiplicative noise. In this work, we simplify a Boussinesq model under location uncertainty assuming strong rotation, stratification, and subgrid turbulence. From this last assumption, the geostrophic balance is modified, and an horizontal divergent velocity explicitly appears. Furthermore, the QG approximation implies a zero PV. In other words, the strong uncertainty prevents interior dynamics at mesoscales. This provides a new derivation of the SQG model from the Boussinesq equations. The ensuing SQG model with divergent velocity is denoted SQG SU . It exhibits physically relevant asymmetry between cold and warm areas, and suggests a diagnostic of the mesoscale divergence from the vorticity, as successfully tested on very high-resolution simulated data.

A more complete model could encompass white noise components for temperature, salinity and density. At mesoscales, a thermal wind relation should relate these time uncorrelated components to the unresolved velocity. Therefore, these additional terms should provide the vertical structure of the unresolved velocity, without increasing the complexity of the parametrization.

Finally, besides solar forcing, the restratification is certainly a complicated process related to frontal dynamics. In the Mixed Layer (ML), the ML instabilities are often triggered by nonhydrostatic motions. They generate very-small-scale baroclinic instabilities and slumpings of the fronts [START_REF] Boccaletti | Mixed layer instabilities and restratification[END_REF]. For such phenomena, subgrid parameterizations are necessary. They must act to horizontally homogenize and restratify the ML. In such a context, the SQG SU model may constitute a simple solution or, at least a first step to develop models under location uncertainty in this direction. To encode the weak stratification of the ML, stochastic Semi-Geostrophic (SG) and Surface Semi-Geostrophic (SSG) models could also be derived. According to our scaling of the vertical unresolved velocity (10), a weaker stratification should then enhance the vertical mixing compared to the SQG SU model. The modified geostrophic balance (17) would involve both horizontal and vertical diffusions. Moreover, since the stratification is weaker in the SG scalings, each term of the buoyancy equation (3D transport, 3D turbulent dissipation and stratification) would have the same scaling. Me an p h ase sh i f t

A Modified geostrophic balance

Under strong horizontal homogeneous turbulence, the large-scale geostrophic balance is modified by the horizontal diffusion:

f × u - a H 2 ∆ H u = ξ, ( 27
)
where u is the resolved horizontal velocity and ∆ H = ∂ 2 x + ∂ 2 y the horizontal Laplacian. On the right-hand side, ξ is the pressure gradient. Let us note that f × u = f J u with J = 0 -1 1 0 and that J T = J -1 = -J . For a constant Coriolis frequency, the previous equation can be solved in the horizontal Fourier space :

û = f J + a H 2 k 2 2 I d -1 ξ = I d - k k c 2 2 J -1 - 1 f ξ⊥ , (28)
with k c = 2f /a H . -1/f ξ ⊥ = -1/f J ξ is the solution without diffusion. Expanding the righthand side operator in Taylor series and using the properties J 2p = (-1) p I d and J 2p+1 = (-1) p J ,

I d - k k c 2 2 J -1 = +∞ p=0 k k c 2 2 J p , (29)
= +∞ p=0 (-1) p k k c 4p 2 I d + +∞ p=0 (-1) p k k c 4p+2 2 J , (30)
= +∞ p=0 - k k c 4 2 p I d + k k c 2 2 J , (31)
= 1 1 + k kc 4 2 I d + k k c 2 2 J . ( 32
)
This leads to the following solution for the modified geostrophic balance:

û = 1 1 + k kc 4 2 - 1 f ξ⊥ + k kc 2 2 1 + k kc 4 2 1 f ξ . (33)

B Non-dimensional Boussinesq equations

To derive a non-dimensional version of the Boussinesq equations under location uncertainty (Resseguier et al., 2017a), each term of the evolution laws is scaled (Resseguier et al., 2017b): 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

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

Υ = U L A u = U 2 A u /T . ( 36
)
C QG model under strong uncertainty

For the case Υ close to the Rossby number, the diffusion term is not negligible anymore and the geostrophic balance is modified. As the terms of the geostrophic balance remain large (Ro Υ 1), the scaling of the pressure can still be done with the Coriolis force. This leads to an Euler number scaling as

Eu ∼ 1 Ro . ( 37
)
Keeping a small aspect ratio D 2 1, we get

Eu D 2 ∼ 1 RoD 2 1 Ro 1 Υ . ( 38
)
As the Rossby number and the ratio Υ are both small in the vertical momentum equation, the inertial terms are dominated by the diffusion term which is itself negligible in front of the pressure term. The hydrostatic balance is hence conserved. The buoyancy scaling still correspond to the thermal winds relation:

Γ ∼ Eu ∼ 1 Ro . (39)
Considering the scaling (σdB t ) z / (σdB t ) H ∼ DRo/Bu for the vertical small-scale velocity, the non-dimensional evolution equations are now given by: equation at the order 0 along z, we obtain

a H 2 ∆∂ z b 0 Bu = ∂ z w 0 = -∇ • u 0 = -∆ ψ0 = - ∆ 2 k 4 c ψ 0 . (49)
The time-correlated part of the 0-th order hydrostatic equation relates the buoyancy to the pressure φ 0 :

a H 2 ∆∂ z b 0 Bu = a H 2 ∆∂ 2 z φ 0 = a H 2 ∆∂ 2 z 1 + ∆ 2 k 4 c ψ 0 . (50)
Gathering these two equations leads to:

∆ + 1 + ∆ 2 k 4 c ∂ z f 0 N 2 ∂ z ψ = 0. ( 51
)
Using the horizontal Fourier transform, it writes:

-k 2 2 + 1 + k k c 4 2 ∂ z f 0 N 2 ∂ z ψ = 0. ( 52
)
Under an uniform stratification, with a fixed value at a specific depth (z = η), and a vanishing condition in the deep ocean (z → -∞), a solution is:

ψ(k, z) = ψ(k, η) exp     N k 2 f 0 1 + k kc 4 2 (z -η)     . (53)
Accordingly, the buoyancy is:

b = ∂ z φ = f 0 1 + k k c 4 2 ∂ z ψ = N k 2 1 + k k c 4 2 ψ. (54)

Figure 3 :Figure 4 :

 34 Figure 3: Temperature (in Celsius degree) for the first (top) and 58 th day (bottom) at high temporal and spatial resolution (∆t = 12h and ∆x = 750m) (left) and after filtering (∆t = 1 day and ∆x = 3km) (right). The black line on the top pictures highlight the region selected for the diagnostic (Colour online).

Figure 5 :

 5 Figure5: From top to bottom: Normalized spectrum of the irrotational velocity component (red) and of our estimate of this component (blue), ratio of the Fourier transform modulus of the divergence to the one of our estimate and phase shift (rad) between the divergence and its estimate. Each of these spectral quantities is averaged on angular spatial frequencies and on the 58 winter days (Colour online).

  G. Vallis. Atmospheric and oceanic fluid dynamics: fundamentals and large-scale circulation. Cambridge University Press, 2006. K. Wyrtki, L. Magaard, and J. Hager. Eddy energy in the oceans. J. Geophys. Res., 81(15): 2641-2646, 1976.

2 DU /Bu. It is prominent at small scales and proportional to the variance tensor, such as the divergent component of the horizontal velocity.

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.

Momentum equations

Buoyancy equation

Effective drift

Here, the time-correlated components and the time-uncorrelated components in the momentum equations have not been separated. 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:

Momentum equations

The operators Del, ∇, and Laplacian, ∆ represent 2D operators. If Ro ∼ Υ , the system is not anymore approximately in geostrophic balance. The large-scale velocity becomes divergent and decoupling the system is more involved. For sake of simplicity, we thus focus on the case of homogeneous and horizontally isotropic turbulence. As a consequence, the variance tensor a is constant in space and diagonal:

The time-correlated components of the horizontal momentum at the 0-th order can be written as:

Then, equation (33) of Appendix A expresses the result in Fourier space. In the physical space, the solution reads:

which is the Helmholtz decomposition of the horizontal velocity u 0 into its rotational and divergent component with a stream function ψ 0 and a velocity potential ψ0 . Differentiating the buoyancy