Long Li

Bruno Deremble

Noé Lahaye

Etienne Mémin

Stochastic data-driven parameterization of unresolved mesoscale eddies

Introduction

Ocean mesoscale eddies contain a large proportion of energy and have an important impact on largescale circulations. They are found everywhere and are particularly intensive in the western boundary currents such as the Gulf Stream and the Kuroshio. However, these eddies are generally not resolved in ocean general circulation models, since the deformation radius in the ocean is at most of the order of 100 km. In particular, the effects of the mesoscale eddies need to be properly parameterized in coarse-resolution ocean models.

The most successful parameterization of mesoscale effect is based on the eddy-induced advection scheme [START_REF] Gent | Isopycnal mixing in ocean circulation models[END_REF][START_REF] Gent | Parameterising eddy-induced tracer transports in ocean circulation models[END_REF][START_REF] Griffies | The Gent-McWilliams skew flux[END_REF], which consists in mimicking the impact of the baroclinic instability by flattening the isopycnal surfaces to transfer the available potential energy of the resolved flow to the eddy kinetic energy at the sub-grid scales. However, such scheme does not account for the backscattering of kinetic energy from small to large scales [START_REF] Bachman | The GM+E closure: A framework for coupling backscatter with the Gent and McWilliams parameterization[END_REF]. Classical eddy viscosity models are introduced in coarse models to mimic the mixing action of the unresolved scales. The associated energy dissipation mechanism is often represented by some functional operators [START_REF] Leith | Atmospheric predictability and two-dimensional turbulence[END_REF][START_REF] Griffies | Biharmonic friction with a Smagorinsky-like viscosity for use in large-scale eddy-permitting ocean models[END_REF][START_REF] Bachman | A scale-aware subgrid model for quasi-geostrophic turbulence[END_REF] that depend on the resolved flow. These approaches seem to improve the effective resolution of eddy-permitting models. Nevertheless, they still lead to excessive decreases of the resolved kinetic energy in large-scale ocean models. The lack of variability due to energy missing is particularly detrimental in ensemble forecasting and data assimilation applications [START_REF] Karspeck | An ensemble ajustment Kalman filter for the CCSM4 ocean component[END_REF][START_REF] Franzke | Stochastic climate theory and modeling[END_REF].

To address these shortcomings, stochastic parameterizations are becoming more and more popular and several schemes have been devised and studied for ocean circulation models. These schemes introduce energy backscattering across scales and enable to increase the internal ocean variability. For instance, [START_REF] Berloff | Random-forcing model of the mesoscale oceanic eddies[END_REF] developed a random-forcing model based on the dynamical decomposition of the flow into a large-scale component and an eddy components; [START_REF] Jansen | Parameterizing subgrid-scale eddy effects using energetically consistent backscatter[END_REF] proposed a stochastic forcing in the form of a negative Laplacian viscosity to inject dissipated energy back to the resolved flow; [START_REF] Grooms | Stochastic superparameterization in quasigeostrophic turbulence[END_REF]; [START_REF] Grooms | Stochastic superparameterization in a quasigeostrophic model of the Antarctic Circumpolar Current[END_REF] proposed a stochastic superparameterization scheme that includes stochastic Reynolds stress terms to backscatter kinetic energy; Porta Mana and [START_REF] Porta Mana | Toward a stochastic parameterization of ocean mesoscale eddies[END_REF]; [START_REF] Zanna | Scale-aware deterministic and stochastic parametrizations of eddy-mean flow interaction[END_REF] proposed an eddy parameterization based on a non-Newtonian stress which depends on the partially resolved scales and their variability.

We propose here a specific stochastic representation of the unresolved flow dynamics based on the location uncertainty (LU) framework [START_REF] Mémin | Fluid flow dynamics under location uncertainty[END_REF]. This representation introduces a random unresolved flow component and relies on a stochastic transport principle. Similar approaches based on the same decomposition have also been recently proposed by [START_REF] Holm | Variational principles for stochastic fluid dynamics[END_REF]; [START_REF] Cotter | Numerically modeling stochastic Lie transport in fluid dynamics[END_REF][START_REF] Cotter | Modelling uncertainty using stochastic transport noise in a 2-layer quasi-geostrophic model[END_REF]; [START_REF] Gugole | Numerical development and evaluation of an energy conserving conceptual stochastic climate model[END_REF]. This kind of random models enable to consider systems that are less dissipative than their deterministic counterparts. Nevertheless, the ability of such models to represent faithfully the uncertainties associated to the actual unresolved small scales highly depends on the definition of the random component and on its evolution in time. Unsurprisingly, stationarity/timevarying and homogeneity/inhomogeneity characteristics of the unresolved flow component have strong influences on the numerical results. For instance, Bauer et al. (2020a) illustrates that the noise inhomogeneity induces a structuration of the large-scale flow reminiscent to the action of the vortex force associated to surface wave-induced Stokes drift; Bauer et al. (2020b) shows that the introduction of inhomogeneous noise into the barotropic quasi-geostrophic (QG) model enables us to reproduce accurately, on a coarse mesh, the high order statistics of the eddy-resolving data; a stochastic shallow water model preserving the resolved total energy has been proposed by [START_REF] Brecht | Rotating shallow water flow under location uncertainty with a structure-preserving discretization[END_REF] and the results show that the stochastic parameterization provides a good trade-off between model error representation and ensemble spread.

In this work, we investigate the LU formulation for the baroclinic QG model [START_REF] Hogg | Mechanisms of decadal variability of the winddriven ocean circulation[END_REF]) within an idealized double-gyre circulation configuration. In particular, we focus on the reproduction of the eastwards jet as well as the prediction of the low-frequency variability for the proposed random model on a coarse resolution grid, for cases in which the baroclinic instability is not resolved. Another, important aspect, concerns the ability to include in the noise representation a stationary drift component associated to the temporal mean of the high-resolution fluctuations. As shown in this paper such stationary drift can be elegantly introduced in the noise through Girsanov transformation. In addition, a projection method has been proposed to update in time the noise along the iso-surfaces of the stratification.

The remainder of this paper is structured as follows. Section 2 describes the framework of LU and presents the derived stochastic QG model. Section 3 explains the discretization of the proposed QG model and the parameterization method for the unresolved noise. Section 4 discusses the numerical results with some statistical diagnosis and energetic analysis. In Section 5 we draw some conclusions and provide an outlook for future works. In the Appendices we demonstrate the energy conservation of the random system and explore the energy conversion between the ensemble-mean and the ensembleeddy components.

Continuous models

In this section, we first review the general setting of Location Uncertainty (LU) [START_REF] Mémin | Fluid flow dynamics under location uncertainty[END_REF], we then present the resulting stochastic Boussinesq quasi-geostrophic (QG) model for a continuously stratified ocean, finally we complete the random formulation by including a time-correlated drift into the unresolved flow component through the Girsanov transformation.

Stochastic flow

The evolution of Lagrangian particle trajectory (X t ) under LU is described by the following stochastic differential equation (SDE):

dX t (x) = v X t (x), t dt + σ X t (x), t dB t , X 0 (x) = x ∈ D, (2.1)
where v denotes the time-smooth resolved velocity that is both spatially and temporally correlated, σdB t stands for the fast oscillating unresolved flow component (also called noise in the following) that is correlated in space yet uncorrelated in time, and D ⊂ R d (d = 2 or 3) is a bounded spatial domain.

Mathematically, {B t } 0≤t≤T is an I d -cylindrical Brownian motion (Da Prato and Zabczyk, 2014) on a filtered probability space (Ω, F, {F t } 0≤t≤T , P)1 and takes values in the Hilbert space H := (L 2 (D)) d . Informally, one can consider its time-derivative as a space-time "white noise" (in a distribution sense) with independent components (B i t ) i=1,...,d . The spatial structure of the unresolved flow component is modeled by the correlation operator, σ, which is not necessarily deterministic. In the most general case it is random and time-dependent with some regularity conditions. More precisely, for each (ω, t) ∈ Ω × [0, T ], σ(•, t)[•] is assumed to be a Hilbert-Schmidt integral operator on H with a bounded matrix kernel σ = (σ ij ) i,j=1,...,d such that

σ(x, t) f = D σ(x, y, t)f (y) dy, f ∈ H, x ∈ D.
(2.2)

The composition of σ(•, t)[•] and its adjoint σ * (•, t)[•] is a trace class operator on H and admits eigenfunctions ξ n (•, t) with eigenvalues λ n (t) satisfying n∈N λ n (t) < +∞. Then, the noise can be equally defined by the following spectral decomposition:

σ(x, t) dB t = n∈N λ 1/2 n (t)ξ n (x, t) dβ n (t), (2.3)
where β n are independent (one-dimensional) standard Brownian motions.

In addition, we assume that the operator-space-valued process {σ(

•, t)[•]} 0≤t≤T is stochastically integrable, i.e. P T 0 n∈N λ n (t) dt < +∞ = 1. As such, the stochastic integral { t 0 σ(•, s) dB s } 0≤t≤T
is a H-valued Gaussian process of zero mean and of bounded (global) variance under the probability measure P:

E P t 0 σ(•, s) dB s = 0, E P t 0 σ(•, s) dB s 2 H < +∞.
(2.4)

Moreover, the point-wise (x ∈ D) and path-wise (ω ∈ Ω) strength of the unresolved flow component at each instant (t ∈ [0, T ]) is measured by the matrix kernel of the composite operator σσ * , and denoted by a(x, t), namely

a(x, t) := D σ(x, y, t) σT (y, x, t) dy = n∈N λ n (t) ξ n ξ T n (x, t).
(2.5)

As shown in Bauer et al. (2020a), the process t 0 a(x, s) ds actually corresponds to the quadratic variation (Da Prato and Zabczyk, 2014) of t 0 σ(x, s) dB s , and it is a continuous random process of finite variation (hence time differentiable). In the particular case of a nonrandom σ, a reduces then to the local variance of the noise according to the Itô isometry [START_REF] Da Prato | Stochastic equations in infinite dimensions[END_REF]. Physically, the symmetric non-negative tensor a represents the (possibly random) friction coefficients of the unresolved fluid motions and the eigenvalues λ n have the unit of m 2 /s.

Stochastic transport

The evolution law of a random tracer Θ with extensive property (e.g. temperature, salinity, buoyancy) transported by the stochastic flow, Θ(X t+δt , t + δt) = Θ(X t , t) with δt an infinitely small time variation, is derived by [START_REF] Mémin | Fluid flow dynamics under location uncertainty[END_REF]; Bauer et al. (2020a) using the generalized Itô formula [START_REF] Kunita | Stochastic flows and stochastic differential equations[END_REF]. It is described by the following stochastic partial differential equation (SPDE):

D t Θ = d t Θ + (v ⋆ dt + σdB t ) • ∇Θ - 1 2 ∇ • (a∇Θ) dt = 0, (2.6)
in which D t is introduced as a stochastic transport operator and d t Θ(x) := Θ(x, t + δt) -Θ(x, t) stands for the (forward) time-increment of Θ at a fixed point x ∈ D. This SPDE encompasses physically meaningful terms. For instance, the third term describes the tracer's advection by the unresolved flow component. As shown in Resseguier et al. (2017a); Bauer et al. (2020a), the resulting multiplicative noise σdB t • ∇Θ continuously backscatters random energy to the system through the quadratic variation 1 /2 (∇Θ) • (a∇Θ) of the random tracer. The last term in (2.6) depicts tracer's diffusion due to the mixing of the unresolved scales. In particular, under specific noise definitions [START_REF] Mémin | Fluid flow dynamics under location uncertainty[END_REF], the resulting diffusion can be connected to the functional eddy viscosity as introduced in many large-scale circulation models [START_REF] Smagorinsky | General circulation experiments with the primitive equations[END_REF][START_REF] Redi | Oceanic isopycnal mixing by coordinate rotation[END_REF].

As an additional feature of interest, there exists an effective advection velocity v ⋆ in (2.6) which is defined as

v ⋆ = v - 1 2 ∇ • a + σ * (∇ • σ). (2.7)
This statistical eddy-induced velocity drift captures the action of inhomogeneity of the random field on the transported tracer and the possible divergence of the unresolved flow component. Bauer et al. (2020a) shows that the turbophoresis term 1 /2 ∇ • a can be interpreted as a generalization of the Stokes drift, which occurs, for example, in the Langmuir circulation [START_REF] Craik | A rational model for Langmuir circulations[END_REF][START_REF] Leibovich | On wave-current interaction theories of Langmuir circulations[END_REF]. Many useful properties of the stochastic transport operator D t have been explored by Resseguier et al. (2017a[START_REF] Resseguier | New trends in ensemble forecast strategy: uncertainty quantification for coarse-grid computational fluid dynamics[END_REF]; [START_REF] Li | Stochastic modeling and numerical simulation of ocean dynamics[END_REF]. In particular, if a random tracer is transported by the incompressible stochastic flow under suitable boundary conditions, then the path-wise p-th moment (p ≥ 1) of the tracer is materially and integrally invariant, namely

D t 1 p Θ p = 0, d t D 1 p Θ p dx = 0.
(2.8)

Stochastic QG model

The derivation of the stochastic geophysical models under the LU framework follows almost exactly the same path as the deterministic derivation [START_REF] Vallis | Atmospheric and oceanic fluid dynamics: fundamentals and large-scale circulation[END_REF]. In particular, a stochastic Euler-Boussinesq model has been derived by Resseguier et al. (2017b); Bauer et al. (2020a) for large-scale atmospheric and oceanic circulations, and reads:

Momentum equations

D t v + f k × (u dt + σ h dB t ) = (b dt + db σ t )k -∇(p dt + dp σ t ), (2.9a) Continuity equations ∇ • v ⋆ = 0, ∇ • σdB t = 0, (2.9b)
Thermodynamic equation

D t b + N 2 (w ⋆ dt + σ z dB t ) = 1 2 ∇ • (a •z N 2 ) dt, (2.9c)
where v = (u, w) T , σdB t = (σ h dB t , σ z dB t ) T , k = (0, 0, 1) T , f = f 0 + βy is the Coriolis parameter using beta-plane approximation, p and dp σ t /dt (in a distribution sense) are the time-smooth component and the fast oscillating noise of the pressure fluctuations rescaled by the background density ρ 0 , N 2 (z) = -g∂ z ρ(z)/ρ 0 is the Brunt-Väisälä frequency with g the gravity value and ρ the stationary density in equilibrium that only depend on height, b(x, t) = -gρ ′ (x, t)/ρ 0 is the resolved buoyancy associated with the density anomaly ρ ′ , db σ t /dt is a zero-mean noise that models the unresolved buoyancy fluctuations and a •z denotes the z-column vector of tensor a. The term db σ t can be seen as a thermodynamics noise coming from the nonlinear mapping law of state applied to the noise term of salinity and temperature transports. It is of the same nature as the noise term introduced in Brankart (2013).

From the system (2.9) a diverse set of approximated models under LU can be obtained through nondimensionalization and asymptotic approach with proper scaling. However, the noise introduces an additional degree of freedom that must be appropriately accounted for [START_REF] Brecht | Rotating shallow water flow under location uncertainty with a structure-preserving discretization[END_REF]Bauer et al., 2020b;Resseguier et al., 2017b,c). The horizontal components of the quadratic variation are first scaled as a h ∼ ϵ U L, where U and L are typical velocity and length scales, and the factor ϵ is proportional to the ratio between the eddy kinetic energy (EKE) and the mean kinetic energy (MKE) and to the ratio between the small-scale correlation time and the large-scale advection time. From the definitions (2.2), (2.3) and (2.5), the horizontal components of the noise can be scaled as σ h dB t ∼ ϵ 1/2 L. The greater this scaling number ϵ, the larger the quadratic variation, hence the stronger the path-wise noise. As shown in Resseguier et al. (2017b,c), using different levels of noise in the stochastic system allows us to model different physical regimes of the large-scale flow.

In this work, we only consider the "moderate" horizontal uncertainty, ϵ = O(1) or ϵ R o = O(R o ) with a small Rossby number R o = U/(f 0 L) ≪ 1, in order to maintain the classical geostrophic balance.

Besides, to keep the assumption of flat isopycnal (small variation of stratification) in the QG theory, Resseguier et al. (2017b) proposed to scale the ratio between the vertical scale and the horizontal scale of the noise as σ z dB t /σ h dB t ∼ (H/L)(R o /B u ) with the height scale H ≪ L and the Burger number B u = (L d /L) 2 = O(1), where L d stands for the scale of deformation radius. Moreover, the unresolved pressure and buoyancy are scaled by dp σ t ∼ Φ = ϵ 1/2 f 0 L 2 and db σ t ∼ Φ/H. All the resolved variables (u, w, p, b, N 2 ) are scaled as in the classical framework [START_REF] Vallis | Atmospheric and oceanic fluid dynamics: fundamentals and large-scale circulation[END_REF] and the variation of Coriolis parameter stays small, i.e. βL < f 0 .

Expanding the resolved variables (v, p, b) in power series of the small Rossby number and truncating the equations at zeroth and first orders, then restoring the dimensions of variables, the following QG system can be derived: Momentum equations

D h t u + f 0 k × u a dt + βyk × (u dt + σ h dB t ) = -∇ h p a dt, (2.10a)
Buoyancy equation

D h t b + N 2 w a dt = 0, (2.10b) Continuity equation ∇ h • u a + ∂ z w a = 0, (2.10c)
where u and p denote the geostrophic components of the resolved velocity and pressure, (u a , w a ) T and p a denote the ageostrophic components, and D h t stands for the horizontal component of the stochastic transport operator, which is defined as

D h t Θ = d t Θ + (u ⋆ dt + σ h dB t ) • ∇ h Θ - 1 2 ∇ h • (a h ∇ h Θ) dt (2.10d) with u ⋆ := u -1 /2 ∇ h • a h and ∇ h = (∂ x , ∂ y ) T .
The geostrophic component and the noise satisfy the following equilibrium:

Geostrophic balances and martingale pressure equilibrium

u = 1 f 0 ∇ ⊥ h p, σ h dB t = 1 f 0 ∇ ⊥ h dp σ t , ∇ h • (u ⋆ -u) = 0, (2.10e) Hydrostatic balances b = ∂ z p, db σ t = ∂ z dp σ t , (2.10f)
where

∇ ⊥ h = (-∂ y , ∂ x ) T .
For sake of conciseness of the notations, we drop the horizontal subscript h in the following. Combining the above equilibrium, the thermal wind balances, ∂ z u ⊥ ∇b and ∂ z σdB t ⊥ ∇db σ t , can be deduced. In this work, we only focus on noise fully satisfying the martingale pressure equilibrium. It means that we model the mesoscale flow as a stochastic component given a large-scale flow. However, such random equilibrium is not necessary in a general setting and could be easily relaxed.

In order to represent the previous random system in terms of potential vorticity (PV) in the same way as in the classical framework [START_REF] Vallis | Atmospheric and oceanic fluid dynamics: fundamentals and large-scale circulation[END_REF], we first take the (horizontal) curl of the firstorder momentum equations (2.10a), then sum the resulting equation with the vertical derivative of the buoyancy equation (2.10b) rescaled by f 0 /N 2 , together with the equilibrium (2.10e), (2.10f) and (2.10c), to obtain finally Stochastic evolution of PV

D t q = i=1,2 -J (u ⋆ dt + σdB t ) i , u i + 1 2 ∇ • ∂ ⊥ x i a∇u i dt -∇ • (a∇f ) dt - f 0 N 2 ∂ z (u ⋆ dt + σdB t ) • ∇b + f 0 2N 2 ∇ • (∂ z a∇b) dt, (2.11a)
From PV to streamfunction

q = ∇ 2 ψ + ∂ z f 0 N 2 b + βy, b = f 0 ∂ z ψ (2.11b) Incompressible constraints u = ∇ ⊥ ψ, ∇ • σdB t = ∇ • (u ⋆ -u) = 0, (2.11c)
where q denotes the PV, ψ := p/f 0 stands for the streamfunction of the resolved flow, J(f, g) = ∂ x f ∂ y g -∂ x g∂ y f stands for the Jacobian operator and

∇ 2 = ∂ 2 xx + ∂ 2
yy is the two-dimensional (2D) Laplacian operator. Note that in the deterministic setting, the sources and sinks of PV in (2.11a) cancels due to the anti-symmetry of the Jacobian operator, J(u i , u i ) = 0, and to the thermal wind balance, ∂ z u ⊥ ∇b. However, it should be noted that the presence of sub-grid terms in a largescale deterministic Boussinesq system leads also in general to additional terms in the RHS of the PV equation. The additional terms in the RHS of (2.11a) are coming exactly from the same sources of difficulty inherent to large-scale flow dynamics representations (i.e. non commutation of curl operator with the sub-grid model, breaking of the Jacobian anti-symmetry as well as possibly non commutation of coarsening processes with derivatives).

Energy conservation and transfers

One important characteristic of the random system (2.10) or (2.11) is that it conserves path-wise (i.e. for each realization) the total energy of the resolved geostrophic flow (under natural boundary conditions), namely

d t D 1 2 |u| 2 + b 2 N 2 dx = D w a b dx - D w a b dx dt = 0.
(2.12)

The demonstration can be found in 5. We remark that this conservation property is consistent with the energy conservation of the stochastic barotropic QG system (Bauer et al., 2020a) and of the stochastic rotating shallow-water system [START_REF] Brecht | Rotating shallow water flow under location uncertainty with a structure-preserving discretization[END_REF].

The energy conversions between the ensemble-mean and ensemble-eddy components are in addition explored in 5. Based on the decompositions u

= E[u] + (u -E[u]) := u + u ′ and b = b + b ′ ,
the mean kinetic energy (MKE), the eddy kinetic energy (EKE), the mean potential energy (MPE) and the eddy potential energy (EPE) are respectively defined by

1 /2⟨|u| 2 ⟩, 1 /2⟨|u ′ | 2 ⟩, 1 /2⟨(b/N ) 2 ⟩ and 1 /2⟨(b ′ /N ) 2 ⟩,
where ⟨•⟩ denotes the integration over domain for simplicity. The energy conversions of the proposed random model (2.10) are summarized in the following diagram.

MKE

EKE MPE EPE -∇u, (u ⋆ ) ′ u ′ F + 1 /2 ∇u, a∇u + a ′ ∇u ′ F + ∇u, (u ⋆ ) ′ u ′ F -1 /2 ∇u, a∇u + a ′ ∇u ′ F -(∇b/N 2 ) • (u ⋆ ) ′ b ′ + 1 /2 (∇b/N 2 ) • (a∇b + a ′ ∇b ′ ) + (∇b/N 2 ) • (u ⋆ ) ′ b ′ -1 /2 (∇b/N 2 ) • (a∇b + a ′ ∇b ′ ) -bw a + bw a -b ′ w ′ a + b ′ w ′ a
Here, [A, B] F = Tr(A T B) stands for the Frobenius inner product of the tensors A and B. The ensemble-mean of the path-wise energy conservation (2.12) can be recovered by summing all the conversion terms.

Compared to the energy budget of the classical deterministic system, the random system provides additional conversion terms between the mean and eddy components. In the general setting with σ random, the classical geostrophic Reynolds stress [START_REF] Bachman | The GM+E closure: A framework for coupling backscatter with the Gent and McWilliams parameterization[END_REF] is modified to (u ⋆ ) ′ u ′ due to the possibly inhomogeneous nature of the noise. This term includes an additional term related to the Ito-Stokes drift that can be positive or negative depending on its alignment with the large-scale velocity component. The fact that a dt = Var(σdB t ) is symmetric and non-negative definite, ensures that the integral ∇u, a∇u F (resp. ∇b, a∇b F ) is always positive, hence provides a positive transfer from MKE to EKE (resp. from MPE to EPE). When the correlation tensor σ is nonrandom, only this term remains in addition to the term related to the classical geostrophic Reynolds stress. As the stochastic system (2.10) does not contain vertical noise component, it does not explicitly modify the conversion mechanisms between KE and PE as observed in the classical system. Nevertheless, the effect of the stochastic terms on this energy transfers is implicitly accounted via the time-integration of the resolved variables b and w a .

Girsanov transformation

The previous formulations (2.10) and (2.11) consist of only a zero-mean and temporally uncorrelated noise. However, this might not be enough in practice and including a mean or time-correlated component of the unresolved velocity field could be of crucial importance to obtain a relevant model. For instance, the eddy parameterization proposed by [START_REF] Zanna | Scale-aware deterministic and stochastic parametrizations of eddy-mean flow interaction[END_REF] is decomposed into a deterministic term based on the negative Laplacian of the PV transport and a zero-mean stochastic term. For the double-gyre circulation configuration, the considered deterministic parameterization allows to reproduce the eastwards jet for the coarse-resolution model, while the additional stochastic terms enhance the gyres circulation and improves the flow variability. Similarly, the random-forcing model proposed by [START_REF] Berloff | Random-forcing model of the mesoscale oceanic eddies[END_REF] consists in a space-time correlated stochastic process to enhance the jet extension.

The slow components of the sub-grid scales can be provided by adequate high-pass filtering of high-resolution data on the coarse grid. We aim in this work at investigating the incorporation of such slow components within the LU framework. However, the derivation of LU models [START_REF] Mémin | Fluid flow dynamics under location uncertainty[END_REF]Resseguier et al., 2017a;Bauer et al., 2020a) relies on the martingale properties of the centered noise and we need hence to properly handle non centred Brownian terms. The Girsanov transformation (Da Prato and Zabczyk, 2014) provides a theoretical tool that fully warrants such a superposition: by a change of the probability measure, the composed noise can be centered with respect to a new probability measure while the additional drift term appears, which pulls back time-correlated subgrid-scale components into the dynamical system. The associated mathematical description is given as follows.

Let Γ t be an H-valued F t -predictable process satisfying

E P exp( 1 2 T 0 ∥Γ t ∥ 2 H dt) < +∞, then the process { B t := B t + t 0 Γ s ds} 0≤t≤T is an H-valued cylindrical Wiener process on (Ω, F, {F t } 0≤t≤T , P)
by a change of the probability measure:

d P dP = exp - T 0 ⟨Γ t , dB t ⟩ H - 1 2 T 0 ∥Γ t ∥ 2 H dt , (2.13)
where ⟨f , g⟩ H = D (f T g)(x) dx denotes the inner product on H = (L 2 (D)) d . Under this new probability measure P, the previous system (2.11) must be adapted to the modified noise σd B t with E P [σd B t ] = 0 by a correction of the effective drift

u ⋆ := u ⋆ -σΓ t . (2.14)
Hereafter, the resulting stochastic system will be adopted with σΓ t referred to as the Girsanov drift.

Numerical models

In this section, we first describe the discretization of the continuous stochastic QG system (2.11) (under a simplification) that will be used later for the numerical simulations. Then, we present numerical methods to estimate the spatial correlation functions of the unresolved flow component based on the spectral decomposition (2.3).

Discrete QG model

We consider here a vertically discretized QG model. This formulation is quite common in geophysical fluid dynamics and the derivation follows the standard methods [START_REF] Vallis | Atmospheric and oceanic fluid dynamics: fundamentals and large-scale circulation[END_REF][START_REF] Mcwilliams | Fundamentals of Geophysical Fluid Dynamics[END_REF]. As illustrated in Figure 1, such model consists of n ocean layers with constant thickness H k and density ρ k in each layer k. The prognostic variables (q k , ψ k , u k , σ k , a k ) are assumed to be layer-averaged quantities.

We first apply finite differences for the buoyancy frequency N 2 and anomaly b on the interface between layers k and k + 1 (for k = 1, . . . , n -1):

b k = f 0 (ψ k -ψ k+1 ) (H k + H k+1 )/2 , N 2 k = g ′ k (H k + H k+1 )/2 , (3.1)
where g ′ k := -g(ρ k -ρ k+1 )/ρ 0 is introduced as the reduced gravity across the interface k. Applying subsequently finite differences for the stratification term ∂ z (f 0 b/N 2 ) with the zero boundary conditions at ocean surface and bottom, b 0 /N 2 0 = b n /N 2 n = 0, the Helmholtz equation (2.11b) can be vertically discretized by

∇ 2 ψ k + f 2 0 H k ψ k-1 -ψ k g ′ k-1 - ψ k -ψ k+1 g ′ k = q k -βy. (3.2a)
In this work, we adopt the dynamical core of the QG coupled model (Q-GCM) developed by [START_REF] Hogg | A quasi-geostrophic coupled model (Q-GCM)[END_REF], which integrates numerically the state variables q and ψ. In order to reduce the inaccuracy of numerical discretization, we do not include the PV sources and sinks (RHS terms of (2.11a)) in this work.

These terms are naturally considered when working in u-b formulations such as in the primitive equations dynamical core [START_REF] Bachman | The GM+E closure: A framework for coupling backscatter with the Gent and McWilliams parameterization[END_REF] or for multi-layers shallow water systems [START_REF] Adcroft | The GFDL global ocean and sea ice model OM4.0: Model description and simulation features[END_REF]. In our case, for the type of noises considered (with moderated amplitude so as to respect the QG scaling performed in this work) these additional noise terms are small but leads to additional numerical difficulties. For stronger noise, with de facto other scalings with potential modification of the hydrostatic balances these additional noise terms have to be considered in the PV equation.

Taking account the numerical dissipation and the vertical entrainment e k (across interfaces), the simplified stochastic evolution of PV reads

D k t q k = (A 2 ∇ 2 -A 4 ∇ 4 )(∇ 2 ψ) + f 0 H k (e k-1 -e k ) dt, (3.2b)
where D k t denotes the k-th layer stochastic transport operator under the probability measure P (hence includes the Girsanov drift described in Section 2.5, ∇ 4 = ∇ 2 (∇ 2 ) is the bi-Laplacian operator associated with the dissipation coefficient A 4 .

As for the spatial discretization, we propose to use a a conservative flux form ∇ • (uq) that corresponds exactly to the 9-points Arakawa Jacobian scheme [START_REF] Arakawa | A potential enstrophy and energy conserving scheme for the shallow water equations[END_REF] of J(ψ, q), in which the zonal and meridional components of the advection flux are defined as 2 /3 u x y q x + 1 /3 u x q d y and 2 /3 v y x q y + 1 /3 v x q d y . Here, • x , • y and • d denote the central averaging operations in the zonal, meridional and diagonal directions respectively. We then apply it for the advection of PV by all the additional drift (such as the noise). On the other hand, we employ simply a central winding scheme to discretize the inhomogeneous diffusion term ∇ • (a∇q). We adopt a stochastic Leapfrog scheme [START_REF] Ewald | Numerical analysis of stochastic schemes in geophysics[END_REF] for the time-stepping of (3.2b). The inversion of the modified Helmholtz equation (2.11b) is carried out with a discrete Fourier transform method [START_REF] Hogg | A quasi-geostrophic coupled model (Q-GCM)[END_REF].

The vertical entrainment e in the classical QG models [START_REF] Berloff | The turbulent ocsillator: A mechanism of low-frequency variability of the wind-driven ocean gyres[END_REF][START_REF] Marshall | A framework for parameterizing eddy potential vorticity fluxes[END_REF][START_REF] Grooms | Stochastic superparameterization in a quasigeostrophic model of the Antarctic Circumpolar Current[END_REF][START_REF] Zanna | Scale-aware deterministic and stochastic parametrizations of eddy-mean flow interaction[END_REF] accounts only for the Ekman pumping from the upper and bottom boundary layers, e 0 = ∇ ⊥ • τ /f 0 and e n = 1 /2 δ ek ∇ 2 ψ n , where τ denotes a given surface wind stress and δ ek stands for a prescribed thickness of the bottom Ekman layer. In the Q-GCM formulation [START_REF] Hogg | Mechanisms of decadal variability of the winddriven ocean circulation[END_REF], an additional entrainment e 1 is considered due to the Ekman pumping of mixed layer temperature. This will be further detailed in Section 4.3.

Parameterizations of unresolved flow

This section describes the numerical methods of the parameterization for the unresolved flow structure. We first review the empirical orthogonal functions (EOF) method which relies on an assumption of quasi-stationarity of the noise covariance. To relax this stationarity constraint, we then propose a projection method enabling to update on-line the noise spatial correlations.

Let {u HR (x, t i )} i=1,...,N be the finite set of velocity snapshots provided by a high-resolution (HR) simulation. We first build the spatial local fluctuations u f (x, t i ) of each snapshot at the coarse-grid points. In particular, for the QG system, one can first perform a high-pass filtering with a 2D Gaussian convolution kernel G on each HR streamfunction ψ HR , to obtain the streamfunction fluctuations, ψ f (x, t i ) = (I -G) ⋆ ψ HR (x, t i ) (only for the coarse-grid points x). Then, the geostrophic velocity fluctuations can be derived by u f = ∇ ⊥ LR ψ f . We next center the data set by u ′ f = u f -u f t (with • t the temporal mean) and perform the EOF procedure [START_REF] Sirovich | Turbulence and the dynamics of coherent structures, part I: Coherent structures[END_REF] to get a set of orthogonal temporal modes {α n } n=1,...,N and orthonormal spatial modes {ϕ n } n=1,...,N satisfying

u ′ f (x, t i ) = N n=1 α n (t i )ϕ n (x), (3.3)
where ⟨ϕ m , ϕ n ⟩ H = δ m,n , α m α n t = λ n δ m,n , H := L 2 (D) and D = A × [-H, 0] denotes the 3D domain with the horizontal area A and the total depth H.

Truncating subsequently the modes (with M ≪ N ) and rescaling them by a small-scale decorrelation time τ , the k-th layer noise and its variance are built by

σ k (x) d B t = τ 1/2 M n=1 λ 1/2 n ϕ n,k (x) dβ n (t), (3.4a) a k (x) = τ M n=1 λ n ϕ n,k (x)ϕ T n,k (x). (3.4b)
Note that this time scale τ is used to match the fact that the noise in (2.1) admits unit of meter. In practice, we often consider the coarse-grid simulation time-step ∆t LR .

In this work, we simply propose to represent the Girsanov drift by the temporal mean of the unresolved eddies, that is

σ k (x)Γ t = N n=1 ⟨u f t , ϕ n ⟩ H ϕ n,k (x). (3.4c)
We remark that in future works this stationary drift term could be generalized to slowly varying fluctuation components or to drift correction with respect to observation [START_REF] Dufée | Stochastic parametrization: An alternative to inflation in ensemble Kalman filters[END_REF]. The previous EOF procedure is an efficient off-line learning method, yet it relies on a strong stationary assumption, and hence leads to a sequence of random flow fields with no connection with the resolved dynamics. In the following, we propose an novel approach that project the EOF-based noise on a resolved state variable. In order to enforce the noise to act only on the resolved momentum (2.10a) and provide efficient backscattering of KE, we consider constraining the noise along the isosurfaces of the stratification, σd B t ⊥ ∇∂ z b with b := f 0 b/N 2 . This projection procedure reads

σ k (x, t) d B t = A 1/2 P k (x, t) σ k (x) d B t , P k := I 2 - ∇(∂ z b) k ∇(∂ z b) k T |∇(∂ z b) k | 2 , (3.5)
where σ k d B t is generated as in (3.4a) and A is a scaling factor to ensure that the noise amplitude is projection invariant. From the definition of the noise quadratic variation (2.5), tensor a defined through (3.4b) is updated as a k = A P k a k P T k . In the next section, we show that the resulting non-stationary noise leads to significant improvement of variability for coarse-resolution simulations.

Numerical results

In this section, we discuss and compare numerical results of several models run within different configurations. The objective consists in improving the variability of low-order models at low Reynolds numbers flow, and in the particular case of unresolved baroclinic instabilities.

Model configurations

We consider here a finite box ocean at mid-latitude driven by an idealized (stationary and symmetric) wind stress τ = (-τ 0 cos(2πy)/L y , 0) T . A mixed horizontal boundary condition is used for the k-th layer streamfunction:

ψ k | ∂A = f k (t) and ∂ 2 n ψ k | ∂A = -(α bc /∆x)∂ n ψ k | ∂A (
same for the 4-th order derivative). Here, f k is a time-dependent function constrained by mass conservation [START_REF] Mcwilliams | A note on a consistent quasigeostrophic model in a multiply connected domain[END_REF], ∆x stands for the horizontal resolution and α bc is a nondimensional coefficient associated to the slip boundary conditions [START_REF] Haidvogel | Boundary current separation in a quasigeostrophic, eddy-resolving ocean circulation model[END_REF].

A deterministic eddy-resolving (∆x = 5 km) model is first simulated and referred to as the references (REF ). It is then compared to several coarse-resolution (∆x = 40 km, 80 km, 120 km) models: the benchmark deterministic model (DET ), two stochastic models with a stationary EOF-based noise (STO-EOF ) and a non-stationary projection noise (STO-EOF-P ). The eddy-resolving model starts from a quiescent initial condition, whereas spin-up conditions downsampled from REF (after 90-years integration) are adopted for all the coarse-resolution models. The common parameters for all the simulations are listed in Snapshots of the relative vorticity on the two top layers are shown in Figure 2. The dynamics of REF model is mainly characterized by the energetic eastward jet with adjacent recirculations and by the fast westward Rossby waves. The former results from the most active mesoscale eddies effect through baroclinic instability. However, this instability cannot be resolved once the horizontal resolution is similar or greater than the baroclinic deformation radius maximum (39 km here). For instance, without any eddy parameterization, the DET (40 km) simulation only generates a smooth field due to an over-dissipation effect.

To calibrate the structure of the unresolved flow (noise and Girsanov drift), the REF data are sampled at each day during the first ten years after the spin-up, this is followed by the procedure presented in Section 3.2 with a fixed 200 km wide Gaussian filter. Considering the computational cost, we employ only 50 EOFs to build the noises (3.4a) for the stochastic simulations at different resolutions. However, the eigenvalues of these truncated EOFs are amplified in order to capture 95 % energy of the full set of modes. An example of the first EOF ( √ λ 1 ϕ 1 ) and the noise variance (a defined in (3.4b)) is illustrated by Figure 3. The noise structure is mainly characterized by the most energetic fluctuations in the jet region.

Including the additional advection of PV by the calibrated unresolved flow component, both STO-EOF and STO-EOF-P models are able to reproduce the eastward jet on the coarse mesh, and the latter seems to enhance fluctuations everywhere. For instance, such result is illustrated at a resolution of 40 km by Figure 2. The ability of these coarse models to reproduce the statistical properties of the REF model will be diagnosed and analyzed more precisely in the next section.

Statistical diagnostics

To investigate the prediction of the long-term statistics, we sample the data (streamfunction) of each run at 15-day intervals over the 120 years (after the spin-up) and compute the statistics over the last 100 years.

We first compare the temporal mean of the modal streamfunctions for both coarse models to that subsampled (in space) of the eddy-resolving model. From Figure 4, we observe that both barotropic and baroclinic modal mean of the REF characterize the eastwards jet, and the two stochastic models (at 40 km) enable to reproduce qualitatively the local structures of both vertical modes predicted by the REF. Conversely, the DET model can only capture the symmetric double-gyre structure.

To quantify the models' statistics at different resolutions, we then employ the spatial root-meansquare error (RMSE) of the temporal mean and the standard deviation (std) between each coarse model and the subsampled REF one. To assess in a single measure the mean and the variance reconstruction, we adopt the Gaussian relative entropy (GRE) proposed by [START_REF] Grooms | Stochastic superparameterization in a quasigeostrophic model of the Antarctic Circumpolar Current[END_REF] that is defined as

GRE = 1 |D| D 1 2 ψ M t -ψ R t 2 σ 2 [ψ M ] + σ 2 [ψ R ] σ 2 [ψ M ] -1 -log σ 2 [ψ R ] σ 2 [ψ M ] dx, (4.1)
where σ denotes the std, ψ M (resp. ψ R ) stands for the streamfunction of the coarse model (resp. for that subsampled from the REF ).

These three global criterion are computed for all the coarse models at different resolutions and the results are summarized in Figure 5. For all the resolutions tested the two random models have lower errors than the DET model both in terms of mean and std. In addition, the STO-EOF-P models enable to reduce the RMSE of the std and the GRE, which means an improvement on the prediction of variability.

We investigate subsequently the models' variability based on the energy analysis. The timeaveraged KE spectra provided by different models are compared in Figure 6. Unsurprisingly, compared to the REF, the DET coarse models produce a severe lack of resolved KE [START_REF] Arbic | On eddy viscosity, energy cascades, and the horizonal resolution of gridded stallite altimeter products[END_REF][START_REF] Kjellsson | The impact of horizontal resolution on energy transfers in global ocean models[END_REF] due to the excessive dissipation without any eddy parameterization. Both stochastic models bring more energy over all the wavenumbers, particularly at large scales. Moreover, when the resolution increases (from 120 km to 40 km), their spectra slopes in the inertial range become more and more closer to that of the REF model.

We compare finally the energy decomposition as described in Section 2 for all the models. However, we focus only on a single realization here and the eddy-mean decomposition of u and b is simply build from a global time-average. From Figure 7, we observe that both random models have higher MKE and MPE than the DET model. In particular, at resolution 40 km, they are almost at the same order as that of the REF model. Compared to the STO-EOF models, STO-EOF-P models produce globally higher EKE and EPE at all the different resolutions. However, they remain much lower than that of the REF ones.

Furthermore, the horizontal density of the EKE and EPE provided by different models are compared in Figure 8. It shows that the STO-EOF-P model improves locally the eddy energy almost everywhere except in the jet region. A more precise parameterization of the unresolved flow is required to explore in future works to improve locally the energy transfers.

In summary, the diagnostic results of this section demonstrate that including the unresolved flow structure inferred from the data into the proposed random models enable to reproduce on the coarse mesh the mean flow of the eddy-resolving simulation. In addition, as performed here through a projection onto the iso-surfaces of vertical stratification, providing an adequate dynamics for the unresolved flow structure is of crucial importance to increase the variability of the coarse models.

Including mixed layer dynamics

In this section, we follow further the Q-GCM formulation [START_REF] Hogg | Mechanisms of decadal variability of the winddriven ocean circulation[END_REF] to embed a mixed layer (of thickness H m = 100 m) in the upper QG layer. This brings an additional diabatic forcing into the PV equation (3.2b) due to the Ekman pumping of sea surface temperature (SST). The vertical velocity across the first ocean interface is e 1 = -1 /2 (∆ m T /∆ 1 T ) e 0 , where ∆ m T = T m -T 1 (resp. ∆ 1 T = T 1 -T 2 ) denotes the temperature difference across the bottom boundary of the mixed layer (resp. of the upper QG layer). Here, the QG potential temperature are fixed to be (287, 282, 276) K.

The additional forcing e 1 evolves in time as the SST does. In the present work, we only considered the geostrophic component of the unresolved noise in the mixed layer. In this case, the stochastic evolution of SST reads

D m t T m + T m ∇ • u m dt = (K 2 ∇ 2 -K 4 ∇ 4 )T m + T m + T 1 2H m e 0 dt, (4.2)
where transport operator under probability measure P (associated to u m ), K 2 = 100 m 2 /s and K 4 = A 4 (see values in Table 2). Note that this shallow-water type equation is derived by vertical integration of a 3D heat equation. Equation (4.2) is numerically discretized by the central winding scheme (in space) and the stochastic Leapfroq scheme (in time).

u m = ∇ ⊥ ψ 1 -τ ⊥ /(f 0 H m )
As shown in Figure 9, the inclusion of the additional forcing coming from the SST evolution generates more small-scale eddies almost everywhere in REF. The improvements of the two stochastic models observed in the previous test case still hold here, namely the zonal jet is reproduced and perturbed. In addition, the variability of the gyres seems also to be enhanced in this case.

Subsequently, an analogue procedure is performed for all the coarse-resolution models in this configuration. For instance, the results of the time-mean streamfunctions at 40 km are provided in Figures 10 anda conclusion consistent with the previous test case is recovered.

As illustrated in Figure 11, the STO-EOF-P models, with a noise along the iso-surfaces of vertical stratification, involves an additional SST random forcing that brings an higher KE backscattering than the stationary noise model at both resolutions. This seems to highlight the importance of the non-stationary characteristic of the noise brought by the projection. As shown in Figure 12, the STO-EOF-P models provide the best results for all the metrics of variability at different resolutions. The local structures of the eddy energy are also better than those recovered in the other coarse models, which is demonstrated by Figure 13.

Conclusions

A stochastic parameterization of the unresolved eddy field is successfully implemented in a well-known QG model. The empirical spatial correlation of the small-scale noise has been first estimated from eddy-resolving simulation data. An additional correction drift that is fully justified in the stochastic setting considered and ensuing from a change of probability measure has been introduced. This non intuitive term seems quite important in the reproduction of the eastward jet within the wind-driven double-gyre circulation. In addition, a projection method has been proposed to constrain the noise living along the iso-surfaces of the vertical stratification. The resulting non-stationary noise model enables us to improve the intrinsic variability of the large-scale resolved flow. This improvement has been demonstrated through some statistical criterion. We have highlighted that both stationary and non-stationary noise models provide improvements at ocean climatic scale in terms of variability metrics, compared to a corresponding deterministic coarse model. At that resolution, the effects of the mesoscale eddies within the ocean basin at midlatitude (of deformation radius around 30 ∼ 40 km) are better represented by the random models, even though the baroclinic instability can not be resolved.

The numerical results presented in this work encourage us to implement the proposed random model on more complex and realistic flow configurations. A following study for a fully coupled model of the ocean and atmosphere [START_REF] Hogg | A quasi-geostrophic coupled model (Q-GCM)[END_REF] is already in progress. Including random forcing into the bottom layer of the atmosphere leads to an unsteady wind stress and to a random diabatic forcing due to Ekman pumping of the sea surface temperature (SST) in both atmospheric and oceanic mixed layers. Thus, it will be interesting to verify the response of the ocean dynamics in terms of variability and energy to these highly variable forcing, and particularly to compare the results to a stationary wind forced ocean dynamics as well as to an unforced dynamics. The stochastic primitive Boussinesq models (2.9) could then be tested. However, in that case the random pressure term in the LU model has to be a priori modeled or computed based on proper scaling of hydrostasy.

Another possible research axis focuses on new parameterization methods for the unresolved flow component. On the one hand, techniques relying on the dynamic mode decomposition [START_REF] Schmid | Dynamic mode decomposition of numerical and experimental data[END_REF][START_REF] Gugole | Spatial covariance modeling for stochastic subgrid-scale parameterizations using dynamic mode decomposition[END_REF], or on the spectral analysis of the Koopman operator [START_REF] Mezić | Analysis of fluid flows via spectral properties of the Koopman operator[END_REF][START_REF] Giannakis | Data-driven spectral decomposition and forecasting of ergodic dynamical system[END_REF] could be explored in practice to estimate the dynamics of the noise and the Girsanov drift from high-resolution data. The objective will consist in evaluating their ability to characterize the long-terms impacts of some small-scale events. Preliminary very encouraging results for the same stochastic QG model as the one studied in the present work can be found in [START_REF] Li | Stochastic parameterization with dynamic mode decomposition[END_REF]. On the other hand, it is also interesting to design adequate noises that provide efficient energy conversion [START_REF] Bachman | The GM+E closure: A framework for coupling backscatter with the Gent and McWilliams parameterization[END_REF][START_REF] Jansen | Toward an energetically consistent, resolution aware parameterization of ocean mesoscale eddies[END_REF] based on the analysis presented in this work. The objective is to release the APE through the noise quadratic variation (2.10b) and to backscatter the unresolved KE to the resolved KE through the advection of momentum by the noise (2.10a).

(2.10b) for the resolved buoyancy b, where the first equality comes from the geostrophic balances (2.10e), the third equality results from both the Stokes theorem with ideal boundary condition of u a and the continuity equation (2.10c), the fourth equality derives from the vertical integration-by-parts formula, and the last equality is based on both the vertical boundary conditions, w a (z = 0) = w a (z = -H) = 0 and the hydrostatic balances (2.10f). Finally, by summing equations (5.5) and ( 5.3), one show that the stochastic system (5.2) conserves the path-wise total energy of the resolved geostrophic flow component, namely

d t D 1 2 |u| 2 + b 2 N 2 dx = D d t |u| 2 2 dx + D d t b 2
2N 2 dx = 0.

(5.6)

Since this conservation property is path-wise, it holds also for their ensemble-mean.

where the first two lines describe the eddy-mean interactions, whereas the last two lines represent the effect of eddy-eddy interactions on the ensemble-eddy component. To evaluate the ensemble-mean of the energy for Θ ′ , the Itô's integration-by-part formula (5.1) is required again. In this case, one can first show that the mean of the joint quadratic variation process for Θ

Figure 1 :

 1 Figure 1: Illustration for the vertical discretization of a three-layered QG system.

Figure 2 :Figure 3 :Figure 4 :

 234 Figure 2: Snapshots of relative vorticity (divided by f0) for the two upper layers (by rows) provided by different simulations (by columns) after 60-years integration.

Figure 5 :

 5 Figure 5: Comparison of variability measures (by sub-figures) for different coarse models (by colors and groups in each sub-figure). Note that the y-axis of the last sub-figure is in log-scale due to poor variability of the DET models.

Figure 6 :

 6 Figure 6: Temporal mean of vertically integrated KE spectra for different models (by colors) at different resolutions (by sub-figures).

Figure 7 :

 7 Figure 7: Comparison of the global MKE, EKE, MPE and EPE (by sub-figures) for different coarse models (by colors and groups in each sub-figure). The corresponding values of the REF model: MKE = 0.23 J/m 3 , EKE = 0.62 J/m 3 , MPE = 4.01 J/m 3 and EPE = 0.34 J/m 3 .

Figure 9 :

 9 Figure 9: Snapshots of relative vorticity (divided by f0) for the two upper layers (by rows) provided by different simulations (by columns) after 60-years integration.

Figure 10 :

 10 Figure 10: Comparison of time-mean contour of barotropic (top row) and first baroclinic (bottom row) streamfunctions for different models (by columns).

Figure 11 :Figure 12 :

 1112 Figure 11: Temporal mean of vertically integrated KE spectra for different models (by colors) at different resolutions (by sub-figures).

d

  t u = -(u ⋆ dt + σdB t ) • ∇u + 1 2 ∇ • (a∇u) + f 0 v a + βyv -∂ x p a dt(5.2a)d t v = -(u ⋆ dt + σdB t ) • ∇v + 1 2 ∇ • (a∇v) -f 0 u a -βyu -∂ y p a dt (5.2b) d t b = -(u ⋆ dt + σdB t ) • ∇b + 1 2 ∇ • (a∇b) dt -N 2 w a dt,(5.2c) together with the geostrophic equilibrium (2.10e), the hydrostatic balances (2.10f) and the continuity equation (2.10c). Note that a dt corresponds to the quadratic variation of the noise σdB t , which is a defined in (2.5). Applying equation (5.1) for the density of APE, b 2 /(2N 2 ), and integrating subsequently over the ocean domain D = A × [0, -H], we have b(u ⋆ dt + σdB t ) • ∇b + 1 2 b∇ • (a∇b) dt -N 2 w a b dt + equality results from the geostrophic balances (2.10e) and the last equality comes from the divergence theorem with ideal boundary conditions, u • n | ∂A = σdB t • n | ∂A = 0 (which leads to u ⋆ • n | ∂A = 0 from equations (2.7) and (2.5)). Similarly to the previous calculation, one can show that 0 v a + βyv -∂ x p a ) dx dt, 0 u a + βyu + ∂ y p a ) dx dt. (5.4b) Summing these two equations, we deduce the evolution of KE: p) • u ⊥ a + ∇ • (up a ) dx dt = -D ∇ ⊥ • (pu ⊥ a ) -p∇ • u a dx dt = -

  ′ reads d⟨Θ ′ , Θ ′ ⟩ t = (∇Θ ′ ) • (a∇Θ ′ ) + (∇Θ ′ ) • (a ′ ∇Θ ′ ) + (∇Θ ′ ) • (a ′ ∇Θ) + (∇Θ) • (a ′ ∇Θ ′ ) + (∇Θ) • (a∇Θ).(5.8b) Subsequently, we can deduce the evolution of the global variance for Θ• (a∇Θ + a ′ ∇Θ ′ ) dx.(5.8c)Applying similarly this result for u and v (resp. for b), we can deduce the evolution of the EKE (resp. EPE), and the results are shown in the diagram in Section 2.4.

Table 1 ,

 1 whereas resolution dependent parameters are presented separately in Table2.

	Parameters	Value	Description
	X × Y	(3840 × 4800) km	Domain size
	H k	(350, 750, 2900) m	Mean layer thickness
	ρ	1000 kg m 3	Density
	g ′ k	(0.025, 0.0125) m s -2	Reduced gravity
	δ ek	2 m	Bottom Ekman layer thickness
	τ 0	2 × 10 -5 m 2 s -2	Wind stress magnitude
	α bc	0.2	Mixed boundary condition coefficient
	f 0	9.375 × 10 -5 s -1	Mean Coriolis parameter
	β	1.754 × 10 -11 (m s) -1	Coriolis parameter gradient
	L d	(39, 22) km	Baroclinic Rossby radii

Table 1 :

 1 Common parameters for all the models.

	∆x (km) ∆t (s) A 4 (m 4 s -1 )
	5	600	2 × 10 9
	40	1200	5 × 10 11
	80	1440	5 × 10 12
	120	1800	1 × 10 13

Table 2 :

 2 Values of grid varying parameters.

Ω is a set of samples, F is a σ-field (a collection of subsets of Ω), {Ft}0≤t≤T is a filtration (a family of sub σ-fields of F that are continuously indexed in time and ordered non-decreasingly), and P is a probability measure.

Acknowledgments

The authors acknowledge the support of the ERC EU project 856408-STUOD. We would like to thank the technical support provided by Charles Deltel and Pranav Chandramouli.

Data Availability Statement

The source codes to reproduce the results are available at https://github.com/matlong/qgcm_lu, all the coarse-resolution data can be found in https://doi.org/10.5281/zenodo.6818751, and the highresolution data can be either reproduced by the restart files or accessed by contacting the corresponding author.

Appendix A. Conservation of path-wise energy

In this section, we show that the stochastic Boussinesq QG system (2.9), without considering the beta effect on the noise, preserves along time the global energy of the resolved geostrophic flow component for each realization. To this end, we first recall the Itô's integration-by-part formula [START_REF] Kunita | Stochastic flows and stochastic differential equations[END_REF]:

where ⟨X, Y ⟩ t

stands for the joint quadratic variation between the processes X and Y with 0 = t n 0 < t n 1 < • • • < t n pn = t a partition of the interval [0, t], and the limit, if it exists, is defined in the sense of convergence in probability. In particular, we have ⟨B, B⟩ t = t.

We now expand equation (2.10a) for the resolved geostrophic velocity u = (u, v) T and equation

Appendix B. Conversion of ensemble energies

In this section, we show briefly the conversions between the ensemble-variance and the energy of the ensemble-mean. To that end, we first decompose each prognostic variable Θ (could be u, v, b or a) of ( 5.2) into a mean component Θ := E[Θ] and an eddy component

The other time-smooth components F dt apart from the geostrophic transport operator are decomposed in the same way. For instance, F = -N 2 (w a + w ′ a ) in the buoyancy equation. However, the noise is only an eddy component since E[σdB t ] = 0. Substituting such decomposition into equations (5.2) and taking subsequently the ensemble-mean of the expanded equation, we obtain

where the second and the third terms on the RHS describe the mean-mean interactions, whereas the last two terms represent the effect of eddy-eddy interactions on the ensemble-mean component. Under natural boundary conditions, on can easily show the time evolution of the global energy for Θ:

(5.7b)

Applying this result for u and v (resp. for b) together with (5.5) (resp. with (5.3)), one can deduce the evolution of the MKE (resp. MPE), and the conversion terms are summarized in the diagram in Section 2.4. Subtracting next (5.7a) from (5.2), we deduce the evolution of Θ ′ : (5.8a)