Werner Bauer

Pranav Chandramouli

Long Li
email: long.li@inria.fr

Etienne Mémin

Stochastic representation of mesoscale eddy effects in

Keywords: Stochastic modelling, Mesoscale eddies, Geostrophic turbulence, Wind-driven circulation

come   L'archive ouverte pluridisciplinaire

Introduction

Mesoscale eddies contain a significant proportion of ocean energy and have an important impact on large-scale circulations. They are found everywhere in the ocean, and are particularly intensive in the western boundary currents like the Gulf Stream and the Antarctic Circumpolar Current. Unfortunately, to fully resolve these eddies in numerical simulations, a horizontal resolution of ∼ 10km is required, which is far too expensive for a large ensemble of realizations or simulations over a long time duration. Neglecting mesoscale eddy effects may lead to strong errors in the evolution of the large-scale dynamics. Therefore, they need to be properly modeled or parametrized.

A classical parametrization approach is to introduce eddy viscosity in coarse models to mimic the action of the computationally unresolved scales while simultaneously ensuring numerical stability by avoiding pile up of energy at the cutoff scale. The explicit dissipation mechanism is often represented either by a harmonic or biharmonic friction term with uniform coefficient, or through functional operators [START_REF] Smagorinsky | General circulation experiments with the primitive equa-840 tions[END_REF][START_REF] Leith | Atmospheric predictability and two-dimensional turbulence[END_REF][START_REF] Griffies | Biharmonic friction with a Smagorinskylike viscosity for use in large-scale eddy-permitting ocean models[END_REF] that depend on the resolved flow. A more widely adopted approach in ocean modeling is to include the Gent-McWilliams parametrization 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] in addition to eddy viscosity, to model the potential energy flux by smoothing the neutral surface height. However, encoding only large-scale dissipation in coarse models often leads to an excessive decreasing of the resolved kinetic energy [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]).

An alternative approach is based on stochastic parametrization [START_REF] Berloff | Random-forcing model of the mesoscale oceanic eddies[END_REF][START_REF] Grooms | Stochastic superparameterization in quasigeostrophic turbulence[END_REF][START_REF] Porta Mana | Toward a stochastic parameterization of ocean mesoscale eddies[END_REF]Cooper and Zanna, 2015;[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], which aims to introduce energy backscattering across scales. These models provide a marked benefit in improving the internal ocean variability, which can be paramount in ensemble forecasting and data assimilation. As a matter of fact, it is well known that models with poor variability usually lead to very low spread of the ensemble [START_REF] Karspeck | An ensemble ajustment Kalman filter for the CCSM4 ocean component[END_REF]. Hence, assimilation systems tend to be over-confident in the model as compared to the observations. However, to overcome numerical instability brought by introducing random forcing, specific tuning parameters are often included in these parameterized models. The success of such tuning methods often do not extend into new flow regimes.

Stochastic parameterization techniques have been proposed for reduced order climate models based on rigorous homogenization techniques [START_REF] Frank | Stochastic homogenization for an energy conserving multi-scale toy model of the atmosphere[END_REF][START_REF] Franzke | Low-order stochastic mode reduction for a realistic barotropic model climate[END_REF][START_REF] Franzke | Low-order stochastic mode reduction for a prototype atmospheric GCM[END_REF][START_REF] Franzke | Stochastic climate theory and modeling[END_REF][START_REF] Gottwald | Stochastic climate theory[END_REF][START_REF] Majda | Models for stochastic climate prediction[END_REF]. These models rely on a scale-separation principle and introduce a linear stochastic Ansatz model with damping terms for the nonlinear small-scale evolution equation. The resulting homogenized dynamics are cubic with correlated additive and multiplicative (CAM) noises. In the absence of scale-separation, the system usually becomes non-Markovian and incorporates memory terms, as shown in the Mori-Zwanzig equation [START_REF] Givon | Extracting macroscopic dynamics: model problems and algorithms[END_REF][START_REF] Gottwald | Stochastic climate theory[END_REF].

Alternatively, [START_REF] Mémin | Fluid flow dynamics under location uncertainty[END_REF] proposed a consistent stochastic framework defined from physical conservation laws. This derivation keeps the full nonlinearity of the system yet relies on a strong temporal scale-separation principle. Within this framework, the Lagrangian velocity is decomposed into a smooth component and a highly oscillating random field. A stochastic transport principle is subsequently derived using stochastic calculus. Notably, the resulting evolution of a random tracer includes a multiplicative random forcing, a heterogeneous diffusion and an advection correction due to inhomogeneity of the random flow component. With these additional terms, a remarkable energy conservation property along time for any realization of the advected tracer still holds (Resseguier et al., 2017a). This stochastic transport principle has been used as a fundamental tool to derive stochastic representations of large-scale geophysical dynamics (Resseguier et al., 2017a,b,c;[START_REF] Chapron | Large-scale flows Cotter[END_REF] in which the missing contributions of unresolved processes are explicitly taken into account.

Similar approaches based on the same decomposition have been also recently proposed by [START_REF] Holm | Variational principles for stochastic fluid dynamics[END_REF]; [START_REF] Gugole | Numerical development and evaluation of an energy conserving conceptual stochastic climate model[END_REF].

The performance of such a random model has been evaluated and analyzed in terms of uncertainty quantification and ensemble forecasting [START_REF] Resseguier | New trends in ensemble forecast strategy: uncertainty quantification for coarsegrid computational fluid dynamics[END_REF] for a surface quasi-geostrophic (SQG) flow. A more efficient spread is produced by the proposed model compared to a deterministic model with perturbed initial condition. As discussed above, this ability is essential for data assimilation applications. Recently, a stochastic barotropic quasi-geostrophic (QG) model has been proposed [START_REF] Bauer | Deciphering the role of small-scale inhomogeneity on geophysical flow structuration: a stochastic approach[END_REF] within this setting to study the structuration effect of the random field on the large-scale flow. Numerical results illustrate that, encoding an inhomogeneous random component into a propagating monochromatic Rossby wave, induces the formation of extra large vortices.

In the present work, the performance of this stochastic barotropic model is assessed for the numerical simulation of an idealized double-gyre wind forcing within an enclosed shallow basin at midlatitude. The wind-driven circulation is a classical simplified problem in oceanography [START_REF] Vallis | Atmospheric and oceanic fluid dynamics: fundamentals and large-scale circulation[END_REF], which produces qualitatively realistic patterns of mesoscale eddies in approximate geostrophic equilibrium. A particular circulation [START_REF] Greatbatch | Four-gyre circulation in a barotropic model with double-gyre wind forcing[END_REF] living in a highly turbulent regime under weak dissipation of potential enstrophy leads to a stationary four-gyre structure in a long-time average sense. We focus here on the ability of the proposed stochastic models to accurately represent at a coarse resolution the four first statistical moments (mean, variance, skewness and kurtosis) of the flow. Comparing this statistical distribution through its four moments to that predicted by the eddy-resolving data enables us to qualify and quantify the accuracy of our stochastic representation of mesoscale eddy effects on large-scale circulation.

The paper is structured as follows. Section 2 briefly reviews the barotropic wind-driven model in adimensional form. Section 3 focuses on the stochastic transport principles and the derived stochastic barotopic vorticity equation.

Section 4 details the data-driven approaches adopted for the modeling of the random small-scale velocity field. Section 5 discusses the numerical results and their long-term statistics. Finally, Section 6 concludes this work and gives some outlook for future research.

Barotropic vorticity equation

In this work, we use a single-layer QG formulation to study the wind-driven circulation in an oceanic basin following [START_REF] Vallis | Atmospheric and oceanic fluid dynamics: fundamentals and large-scale circulation[END_REF]. Under this regime, the dimensional barotropic vorticity equation (BVE) can be written as:

∂ω ∂t + J(ψ, ω) + β ∂ψ ∂x = F + D, (2.1a) ∇ 2 ψ = ω, (2.1b)
where

ω = k • ∇ × u = ∂ x v -∂ y u
is the relative (or kinematic) vorticity (henceforth, referred to as vorticity) with k = [0, 0, 1] T . The oceanic geostrophic velocity u can be defined by a stream function

ψ such that u = ∇ ⊥ ψ = [-∂ y ψ, ∂ x ψ] T .
100

The nonlinear advection is transformed into a Jacobian operator which is defined as J(ψ, ω) = ∂ x ψ∂ y ω -∂ y ψ∂ x ω. The linear term β∂ x ψ describes the advection of β-planetary vorticity. An active tracer in this case is given by the potential vorticity (PV) defined as q = ω + βy.

On the right-hand side (RHS) of (2.1a), F = k • ∇ × τ /(ρH) is a forcing which adds vorticity into the gyres, due to the wind stress τ over the ocean surface, where ρ and H are respectively (resp.) the basic fluid density and depth of the basin. An idealized double-gyre wind stress [START_REF] Greatbatch | Four-gyre circulation in a barotropic model with double-gyre wind forcing[END_REF], defined only in zonal direction, is used in this work within the basin

Ω = [0, L] × [-L, L], that is 110 τ = [τ 0 cos( πy L ), 0] T , (2.2)
where τ 0 is the magnitude of the wind. This form of wind stress [START_REF] San | Approximate deconvolution large eddy simulation of a barotropic ocean circulation model[END_REF][START_REF] San | Approximate deconvolution large eddy simulation of a stratified two-layer quasigeostrophic ocean model[END_REF] represents the meridional profile of easterly trade winds, mid-latitude westerlies, and polar easterlies from south to north over the ocean basin.

The boundary layer friction D can be interpreted either as a linear drag for the Ekman layer as presented in the Stommel problem [START_REF] Stommel | The westward intensification of wind-driven ocean currents[END_REF], an eddy viscosity term as presented in the Munk problem [START_REF] Munk | On the wind-driven ocean circulation[END_REF], or a combination of the two [START_REF] Fox-Kemper | Reevaluating the roles of eddies in multiple barotropic wind-driven gyres[END_REF]. In this work, we are more interested in the Munk model, by assuming that the ocean has a flat-bottom. The eddy viscosity that we will discuss in the following will be either harmonic D = ν 2 ∇ 2 ω or biharmonic D = -ν 4 ∇ 4 ω, with a uniform coefficient ν 2 (of unit m 2 s -1 ) or ν 4 120 (of unit m 4 s -1 ).

To simplify the problem, one may scale the equation (2.1a), by comparing each term to the dominant β-effect [START_REF] Vallis | Atmospheric and oceanic fluid dynamics: fundamentals and large-scale circulation[END_REF]. The leading order is given by the Sverdrup balance between the rotation and wind forcing, i.e. β∂ x ψ ≈ |F |, which provides a characteristic size of velocity:

V = τ 0 ρH π βL .
(2.3a)

This leads to the following scaling of time, vorticity and stream function:

t = L V t , ω = V L ω , ψ = V Lψ , (2.3b)
where the variables with prime symbol ( ) are adimensionalized.

The thickness of the Munk boundary layer can be then quantified by the balance between the β-effect and friction [START_REF] Munk | On the wind-driven ocean circulation[END_REF]. For instance, β∂ x ψ ≈ ν 2 ∇ 2 ω gives us a harmonic-boundary-layer scale, that is

130 δ 2 = ν 2 β 1/3 . (2.3c)
Similarly, β∂ x ψ ≈ ν 4 ∇ 4 ω gives us a biharmonic-boundary-layer scale:

δ 4 = ν 4 β 1/5
.

(2.3d)

The nonlinear advection term J(ψ, ω) is smaller than the linear terms. Nevertheless, the nonlinear effect may still be important in the boundary layer, especially in the western one. To measure its strength, one may define a β-Rossby number (denoted as R β ) as the ratio of the size of the nonlinear term to the β-effect:

R β = V βL 2 .
(2.3e)

Using these scaling numbers (2.3a)-(2.3e) for (2.1a), the dimensional BVE reduces to its adimensional form as:

∂ω ∂t + J(ψ , ω ) + 1 R β ∂ψ ∂x = 1 R β sin(πy ) + 1 R β D, (2.4) with D = (δ 2 /L) 3 ∇ 2 ω or D = -(δ 4 /L) 5 ∇ 4 ω resulting from (2.3c) or (2.3d),
respectively. The adimensional PV is written as q = R β ω + y , and the Poisson 140 equation (2.1b) is invariant under this adimensionalization, i.e. ∇ 2 ψ = ω . For the sake of readability, in the following we drop the prime for all the adimensional variables.

To close the problem, we need one initial condition -that will be discussed in section 5.1 -and two boundary conditions. The first boundary condition is imposed by the no-normal-flow condition due to the forcing form:

ψ| ∂Ω = 0, i.e. u| x=0,L = v| y=-L,L = 0, (2.5a)
where ∂Ω denotes the basin's boundary. The second one depends on the chosen eddy viscosity form. For a harmonic friction, i.e. D = (δ 2 /L) δ 4 /L).

Stochastic barotropic vorticity equation

This section provides a stochastic representation of the barotropic QG flow.

We start by introducing the stochastic Lagrangian flow (X ∈ Ω ⊂ R d , d = 2, 3)

given by [START_REF] Mémin | Fluid flow dynamics under location uncertainty[END_REF]: given by

dX t = u(X t , t)dt + σ(X t ,
Q(x, y, t, s) = E σ(x, t)dB t σ(y, s)dB s T = δ(t -s)dt Ω σ(x, z, t) σT (y, z, s)dz, (3.2b)
where E stands for the expectation and the last equality ensues from Itô isometry [START_REF] Da Prato | Stochastic equations in infinite dimensions[END_REF]. The variance (or auto-covariance) tensor, denoted as a, is defined by the diagonal components of the covariance per unit of time, a(x, t) = Q(x, x, t, t)/dt = σσ T (x, t), which has the unit of a diffusion tensor (m 2 s -1 ). In addition, the density of the turbulent kinetic energy (TKE) under this framework can be defined by 1 2 tr(a)/dt that has a unit of m 2 s -2 . The previous representation (3.2) is a general way to define the small-scale flow. In particular, the fact that σ is Hilbert-Schmidt, ensures that the covariance operator per unit of time, Q/dt, admits an orthogonal eigenfunction basis {Φ n (•, t)} n∈N weighted by the eigenvalues Λ n ≥ 0 such that n∈N Λ n < ∞.

Therefore, one may equivalently define the small-scale flow based on the following spectral decomposition (Da Prato and Zabczyk, 2014):

σ(x, t)dB t = n∈N Φ n (x,
Φ T m (x, t) = n∈N Φ n (x, t)Φ T n (x, t), (3.3b)
where δ n,m denotes the Kronecker symbol.

Hereafter, the rate of change of a random scalar process θ, within a volume V transported by the stochastic flow (3.1), can be deduced from the Itô-Wenzell 190 theorem [START_REF] Kunita | Stochastic flows and stochastic differential equations[END_REF]. Under the incompressible assumption for the smallscale flow, ∇•σ = 0, it can be written in Eulerian coordinates as d

V(t) θ(x, t)dx = V(t) (D t θ + θ ∇• (u -u s ))dx, (3.4a) D t θ = d t θ + (u -u s ) •∇ θdt + σdB t •∇ θ - 1 2 ∇• (a∇θ)dt, (3.4b)
where D t is introduced as a stochastic transport operator (Resseguier et al., 2017a). Note that d t θ = θ t+dt -θ t stands for the forward time-increment of the scalar θ at a fixed point x. The turbophoresis term, u s = 1 2 ∇• a, accounting for the effect of statistical inhomogeneity of the small-scale field on the largescale current, is referred to as the Itô-Stokes drift in [START_REF] Bauer | Deciphering the role of small-scale inhomogeneity on geophysical flow structuration: a stochastic approach[END_REF]. This term was shown to play a crucial role in the transition from the viscous layer regime to the logarithmic layer regime in wall bounded turbulent flows [START_REF] Pinier | Stochastic flow approach to model the mean velocity profile of wall-bounded flows[END_REF]. It can be considered as a generalization of the Stokes drift, which 200 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]. As shown in [START_REF] Mémin | Fluid flow dynamics under location uncertainty[END_REF], under a spatially heterogeneous and temporally non-stationary random field in general, the last term in (3.4b) plays a role similar to the functional eddy viscosity as introduced in many largescale circulation models [START_REF] Smagorinsky | General circulation experiments with the primitive equa-840 tions[END_REF][START_REF] Gent | Isopycnal mixing in ocean circulation models[END_REF]. In particular, for a homogeneous, isotropic and stationary random field, in which the variance tensor a becomes a 0 I d , the diffusive term boils down immediately to a harmonic friction term, 1 2 a 0 ∇ 2 θ, with a uniform coefficient a 0 to be specified. In order to ensure an isochoric flow, an incompressibility constraint on the corrected large-scale drift, ∇•(u-u s ) = 0, is additionally required. A stochastic 210 transport equation of the extensive tracer θ is directly deduced from (3.4a), D t θ = 0.

(3.5a)

In Resseguier et al. (2017a), it is shown that those incompressibility constraints enable us to establish an energy conservation property:

d dt Ω 1 2 θ 2 dx = 1 2 Ω θ ∇• (a∇θ)dx
Energy loss by diffusion

+ 1 2 Ω (∇θ) T a∇θdx
Energy intake by noise = 0, (3.5b) in which, for any realization of the random tracer, the global energy brought by the small-scale flow is exactly compensated by that dissipated by its diffusive contribution (within ideal boundary conditions). Note that the energyincreasing term arises from Itô integration by part formula.

The derivation of the stochastic barotropic vorticity equation, fully detailed in [START_REF] Bauer | Deciphering the role of small-scale inhomogeneity on geophysical flow structuration: a stochastic approach[END_REF], follows a similar strategy as in the classical framework.

10

The main steps of the derivation procedure are: first, the three-dimensional stochastic mass and momentum equations are obtained by applying the stochastic transport principle (3.4a); then, a two-dimensional stochastic rotating shallow water system is deduced from the classical hydrostatic assumption; subsequently, substituting the unknown variables, written as a power series of (small)

Rossby number, into the dimensionless equations, we get the asymptotic solutions of each order. Introducing the wind forcing, the eddy viscosity, and assuming an infinite Rossby radius of deformation (poor height stratification), the dimensional stochastic barotropic vorticity equation (SBVE) reads

d t ω + J(ψ, ω)dt + k • ∇ × dM t + β ∂ ∂x (ψdt + ϕdB t ) = (F + D)dt, (3.6a) dM t = (σdB t -u s dt) •∇ u - 1 2 ∇• (a∇u)dt. (3.6b)
The process dM t gathers the additional momentum terms introduced in the stochastic transport equation (3.4b). Due to the geostrophic balance and the Doob-Meyer decomposition theorem [START_REF] Kunita | Stochastic flows and stochastic differential equations[END_REF], the small-scale flow is defined from a random stream function ϕdB t as σdB t = ∇ ⊥ ϕdB t [START_REF] Bauer | Deciphering the role of small-scale inhomogeneity on geophysical flow structuration: a stochastic approach[END_REF]. The curl of such a process can be expanded as

k • ∇ × dM t = J(ϕdB t , ω) -u s •∇ ωdt - 1 2 ∇• (a∇ω)dt + dS t , (3.6c) dS t = i=1,2 J(σdB i t -u i s dt, u i ) - 1 2 ∇• (∂x ⊥ i a∇u i )dt. (3.6d)
where dS t stands for the source/sink process of the vorticity, due to the rotating interactions between the strains of the large and small scale flows (Resseguier et al., 2017b). The first term in dS t has a similar form as the additional term introduced in the barotropic Leray α-model studied in [START_REF] Holm | Modeling mesoscale turbulence in the barotropic double-gyre circulation[END_REF].

In addition, we highlight from [START_REF] Bauer | Deciphering the role of small-scale inhomogeneity on geophysical flow structuration: a stochastic approach[END_REF] that, without any forcing and damping, the proposed model preserves the total energy (which reduces in this work to the kinetic energy) of the large-scale flow within ideal boundary conditions. A more compact form of SBVE (3.8) can be obtained under Stratonovich stochastic integrals. In the following we give its expression in an adimensionalized form.

Besides the scaling numbers given in Section 2, we need to scale the variance tensor, a = Aa , to precise the strength of uncertainty included in the SBVE.

As mentioned above, since a has the unit of a diffusion tensor m 2 s -1 , one may consider that A is proportional to V L up to a factor , i.e. A = V L. Hereafter, this scaling number, , can be related to the ratio between the TKE, A/T σ , and the mean kinetic energy (MKE), V 2 , and proportional to the ratio between the small-scale correlation time, T σ , and the large-scale one, T (Resseguier et al., 2017b). This reads:

= T σ T TKE MKE . (3.7a)
This leads to the following scaling of variance tensor and small-scale flow:

a = V La , σdB t = √ LσdB t . (3.7b)
The greater this scaling number the larger the variance tensor and the stronger the uncertainty. Furthermore, as interpreted in Resseguier et al. (2017b) and [START_REF] Bauer | Deciphering the role of small-scale inhomogeneity on geophysical flow structuration: a stochastic approach[END_REF], the geostrophic balance is valid only for weak ( 1) to moderate ( ∼ 1) uncertainty in the stochastic case. Beyond this scaling the geostrophic balance is eventually modified and includes correction terms to isobaric velocities. In the present work, only moderate uncertainty is adopted.

Under such an assumption, the final dimensionless SBVE in Stratonovich notation is written as To close the problem, we assume that the small-scale component σdB t and the Itô-Stokes drift u s have the same boundary conditions as the large-scale current u , given in (2.5a) and (2.5d).

d t ω+J ψdt+ √ ϕdB t , ω -u s •∇ωdt+ 1 R β ∂ ∂x (ψdt+ √ ϕdB t ) = 1 R β (F +D)dt-dS t , (3.8a) dS t = i=1,2 J( √ σdB i t -u i s dt, u i ). ( 3
It can be remarked that cancelling the source term (3.8b) and the Itô-Stokes drift in (3.8), we obtain a stochastic potential vorticity equation that corresponds exactly to the model proposed in [START_REF] Cotter | Numerically modeling stochastic Lie transport in fluid dynamics[END_REF], built upon imposing a strong circulation conservation constraint [START_REF] Holm | Variational principles for stochastic fluid dynamics[END_REF]. By defini-

Data-driven modeling of uncertainty

In order to perform a numerical simulation of the SBVE (3.8), the uncertainty field σdB t has to be a priori modeled. This results from (3.3) to construct the eigenfunction basis of the spatial covariance. In practice, we work with a finite set of eigenfunctions of the small-scale Eulerian velocity rather than with the Lagrangian displacement. Data-driven approaches are presented in this section to estimate these empirical basis functions. The first method is based on the proper orthogonal decomposition (POD) method where the covariance is assumed to be quasi-stationary. Moreover, we propose in Section 4.2 a second approach which introduces time-dependent weight coefficients into the spectral decomposition.

Pre-processing of data

In order to estimate the basis functions for coarse SBVE model using (highresolution) eddy-resolving data, a coarse-grainning procedure is required. To this end, a collection of stream function snapshots {ψ HR (x, t i )} i=1,...,Nt , provided by a high-resolution simulation of the BVE (2.4) with grid spacing ∆ HR , are first filtered to a coarser resolution of grid spacing ∆ LR through a Gaussian filter:

ψ HR (x, t i ) = 6 π∆ 2 Ω exp - 6(x -y) 2 ∆ 2 ψ HR (y, t i )dy, (4.1)
with width ∆ = 2∆ LR /∆ HR . The filtered snapshots ψ HR , are subsequently subsampled to give the reference data ψ o (also referred to as observation data in the following) at the coarse resolution. The reference velocity snapshots

{u o (x, t i )} i=1,...,Nt are then deduced from ∇ ⊥ ψ o .

POD method

Applying the snapshot POD procedure [START_REF] Sirovich | Turbulence and the dynamics of coherent structures, part I: Coherent structures[END_REF] (given in appendix A) for the fluctuations u o = u o -u t o (where the overbar ( t ) denotes a temporal average), enables us to build a set of (mutually) orthonormal spatial modes (of unit ms -1 ) {φ i } i=1,...,Nt , and a set of orthogonal temporal modes {b i (t j )} i,j=1,...,Nt associated with a set of decaying eigenvalues {λ i } i=1,...,Nt . In addition, we suppose that such a set of empirical eigenfunctions has a complete (or direct) decomposition [START_REF] Mémin | Fluid flow dynamics under location uncertainty[END_REF]Resseguier et al., 2017d) such that the fluctuations u of the large-scale current lives in a subspace spanned by {φ i } i=1,...,M0-1 , and the small-scale random drift σdB t /∆t with a sufficiently small time step ∆t lives in the residual subspace spanned by

{φ i } i=M0,...,M1 with M 0 < M 1 ≤ N t such that 1 ∆t σ(x)dB t ≈ M1 m=M0 λ m φ m (x)ξ m , (4.2a)
where ξ m are i.i.d. standard Gaussian variables. The corresponding variance tensor is then given by

1 ∆t a(x) ≈ M1 m=M0 λ m φ m (x)φ T m (x). (4.2b)
Therefore, such a POD approach depends only on two parameters: M 0 and M 1 .

The choice of these parameters depends on the energy ratio γ 0 , respectively γ 1 , with 0 < γ 0 < γ 1 < 1, that needs to be captured by the largest, respectively the smallest, spatial scales of the random flow component. More precisely, let us first introduce the so-called relative information content (RIC) of the eigen decomposition:

RIC(m) = m i=1 λ i Nt i=1 λ i , m = 1, . . . , N t . (4.3a)
Suppose that the largest structure of the random flow is required to contain the ratio γ 0 of the total energy of the fluctuations, the first truncated mode is then determined by

M 0 = min{ m | RIC(m) ≥ γ 0 }, (4.3b)
as shown in Figure 1. Analogously, the last truncated mode M 1 can be found with a given ratio γ 1 for the smallest structure. In practice, this latter can be Now, the problem boils down to choose adequately the ratio γ 0 . We propose to estimate it by comparing the kinetic energy spectrums, between the ensemble of observation data {u o (x, t i )} i=1,...,Nt and an extra collection of snapshots {u LR (x, t i )} i=1,...,Nt , obtained from a simulation of the BVE at the coarse resolution ∆ LR . The parameter γ 0 is approximated by the proportion of the partial energy, accumulated up to the first wavenumber κ 0 for which the two temporally averaged spectrums start to deviate (c.f. Figure 2):

γ 0 ≈ κ≤κ0 Êo t (κ) κ≤κc Êo t (κ) , (4.3c)
where Êo denotes the instantaneous kinetic energy spectral density of the observations, and κ c = π/∆ LR stands for the theoretical effective cutoff.

340

Note that both the free-slip boundary conditions and the divergence-free constraint imposed in the previous section, are well-satisfied for the parameterized random velocity (4.2a). Indeed, the proposed spatial modes are represented as a linear combination of the instantaneous observed velocity fields (see (6.1d)).

Mode matching method

The previous POD procedure is an efficient off-line learning method, yet it relies on a strong stationary assumption, and thus leads to a sequence of random velocity fields with no temporal connection with the resolved dynamics. In the following, we propose a novel approach that introduces a time-dependent weight coefficient α m (t) in the POD representation. In this approach, the instantaneous random velocity at each time t is now defined as

1 ∆t σ(x, t)dB t ≈ M1 m=M0 λ m φ m (x) √ α m (t)ξ m , (4.4a)
with the corresponding variance tensor given by

1 ∆t a(x, t) ≈ M1 m=M0 λ m φ m (x)φ T m (x)α m (t). (4.4b)
Indeed, such a weighting provides an energy re-distribution of the spatial modes at each time step. The weighting principle proposed here, consists of selecting from the reference data the set of time instances that match to the large-scale structure of the current simulation. To be more specific, let us consider a current velocity field u l (x, t) at a given time t of the SBVE simulation. The projection coefficient b l 1 of the current fluctuation u l on the first spatial mode φ 1 is defined by

b l 1 (t) = u l (•, t), φ 1 Ω , (4.5a)
where the fluctuation u l at one position are obtained by subtracting a local 360 average of the current field around that position, and where u, v Ω = Ω u•vdx denotes the L 2 (Ω)-inner product. As illustrated in Figure 3, a collection of matching instants is constructed by identifying the current projection b l 1 to the time series of the first temporal mode {b 1 (s)} s=1,...,Nt subject to a consistent condition of its time increments:

S(t) = s |b 1 (s) -b l 1 (t)| ≤ c; ∆ s [b 1 ]∆ t [b l 1 ] ≥ 0 , (4.5b)
where c is a sufficiently small threshold and ∆ t

[b] = b(t) -b(t -∆t) stands
for the temporal variation of b at time t. This aims at selecting the events corresponding to the same projection coefficient and the same sign of the time increment. The weight coefficient α m for each mode m = M 0 , . . . , M 1 is then fixed from the sample variance:

370 α m (t) = 1 |S(t)| -1 s∈S(t) (b m (s) -µ m (t)) 2 , µ m (t) = 1 |S(t)| s∈S(t) b m (s), (4.5c)
where |S| stands for the sample size. These time dependent coefficients allows us to slave a set of modes on some dominant modes. Note that in the present study, we work only with the first mode, however this technique could be extended to a vector of dominant modes in order to select more complex turbulent events. Let us also outline that the boundary conditions and the divergence-free constraint of the random flow (4.4a) remain valid with this weighting method.

Numerical results

In this section, we discuss and compare the respective numerical simulations of the BVE (2.4) and the SBVE (3.8). The main motivation here is to numerically assess if the proposed random model reproduces well the long-term 380 statistics of the high resolution (eddy-resolving) simulation.

All the models have been discretized with the same numerical schemes. As detailed in Appendix B, a staggered Arakawa C-grid [START_REF] Arakawa | Computational design of the basic dynamical processes of the UCLA general circulation model[END_REF] has been considered. In that respect, the nonlinear Jacobian terms in the governing equations are discretized using Arakawa's 9-points conservative scheme (Arakawa and Lamb, 1981). To invert the Poisson equation (2.1b) associated to the stream function, an efficient discrete sine transform solver [START_REF] Press | Numerical recipes 3rd edition: The art of scientific computing[END_REF]) is adopted. For the time-stepping, a strong stability preserving 3rd order Runge Kutta scheme [START_REF] Gottlieb | On higher order strong stability preserving Runge-Kutta and multi-step time discretizations[END_REF] with a Courant-Friedrichs-Lewy (CFL) number of 1/3 is considered for BVE. As further detailed in Apendix B, for the 390 SBVE we used a similar time integration scheme.

Model configurations and simulations

In all the configurations we fix the basin length to L = 1 and the Rossby number to R β = 0.06 2 . For the SBVE simulations, the uncertainty strength parameter has been fixed to ε = 1. For the high resolution eddy-resolving model, a regular mesh with 256 × 512 cells with uniform grid spacing ∆ HR = 0.004 and a five times wider harmonic boundary layer δ 2 = 0.02 have been used. We consider a quiescent state as the initial condition, that is ψ(x, t = 0) = 0, ∀x ∈ Ω. For such an initial condition, the dominant Sverdrup balance between the forcing and rotation leads to a symmetric PV field during a short period. As the nonlinear inertial term becomes more and more important, a symmetry breaking phenomena occurs (at t ≈ 2), which can be observed from the time series of the global kinetic energy in figure 4. This so-called spin-up period is then followed by a dissipation stage (up to t ≈ 5) of the very high enstrophy that has been produced during the spin-up. Immediately after, the flow dynamics becomes rapidly turbulent. The three subsequent snapshots in Figure 4 illustrate this vigorous eddying nature.

At coarse resolution, the subgrid dissipation model is defined through a biharmonic friction term with a grid-dependent uniform coefficient. We choose such a simple, yet commonly used, eddy-viscosity scheme to single out the effects of the proposed random model and the sub-grid dissipation. Besides, as shown in Appendix C, our model provides a very useful technique to estimate the uniform coefficient δ 4 from the high-resolution data ω o = ∇ 2 ψ o . The idea consists in fixing the amplitude of a specific noise (with a corresponding noise diffusion of biharmonic form) such that its energy matches to the observed turbulent kinetic energy. This simple estimation is a very interesting by-product of our stochastic setting. For instance in our case, the estimated values of δ 4 at coarse resolutions 64 × 128, 32 × 64 and 16 × 32 are, respectively, 0.026, 0.040 and 0.049. In practice, all these values have shown to be very good estimates for the subgrid dissipation.

The numerical simulations of the SBVE are performed using both the POD (denoted as SBVE POD ) technique and Mode Matching (denoted as SBVE MM ) approach. In both simulations, the spatial modes for different coarse resolutions are trained during the same period consisting of 6000 snapshots. The energy proportion parameter γ 0 of the first truncated mode for the random flow

E(t) = 1
2 Ω (u 2 + v 2 )dx and the global enstrophy is defined by Z(t) = 1 2 Ω ω 2 dx. The plots show their graph normalized by their temporal maxima. is, resp., estimated at 0.95 and 0.92 at resolution 64 × 128 and 32 × 64 (same for 16 × 32). As shown in Figure 5, by introducing randomness into the initial symmetric double-gyre circulation, the symmetry breaking state is reached much earlier for the SBVE simulations, than for the BVE. Hereafter, in order to compare the different models and to reduce the spin-up errors, we use the coarse-grained version of one specific eddy-resolving snapshot (after t = 5) as the initial condition for all coarse model runs. In other words, the BVE and the SBVE at each coarse resolution are simulated from the very same initial field, in which the spin-up period is accounted for at the eddy-resolving resolution. An instantaneous illustration of the small-scale random stream function, denoted as ψ r = 1 ∆t ϕdB t , and the Itô-Stokes stream function ψ s , is shown in Figure 6. It appears that both ψ r and ψ s based on MM are stronger and more regular than those based on POD.

Long-term prediction of statistics

Although we are working in a turbulent regime, the statistics of the largescale tracers ψ and q tend to reach a statistical steady state equilibrium. As shown in [START_REF] Greatbatch | Four-gyre circulation in a barotropic model with double-gyre wind forcing[END_REF], a robust four-gyre structure is characterized in time-averaged circulation, as long as the dissipation is sufficiently weak. Here, a weak dissipation means that the boundary layer size δ 2 or δ 4 has a smaller order than the so-called Rhines scale R β [START_REF] Vallis | Atmospheric and oceanic fluid dynamics: fundamentals and large-scale circulation[END_REF]. However, this does not indicate that the flow dynamics are under resolved. Note that in under resolved simulations, the contour lines of the averaged tracers would be oscillating. On the other hand, increasing the explicit dissipation up to the order of Rhines scale, would result in a conventional double-gyre.

In this work, apart from the mean structure, we are also interested in the eddy energy distributions and higher order moments of the tracers, such as skewness and kurtosis. These two standard moments of a probability distribution characterize the asymmetry and extreme events, respectively. They are particularly informative when the distribution is non-Gaussian. The associated large-scale velocity field is indicated here by the black arrows. Note that these velocity values are located on the PV-grid (see figure 15) through linear interpolations. More precisely, the first four central moments of ψ are defined by 3, 4, (5.1a) where the superscript ( k ) denotes the power, while the subscript ( k ) denotes the order of the moment order. Similarly, the central moments of q as function of the prognostic variable ω are defined by

m 1 [ψ] = ψ t , m k [ψ] = (ψ -m 1 [ψ]) k t , k = 2,
m 1 [q] = R β m 1 [ω] + y, m k [q] = R k β m k [ω], k = 2, 3, 4.
(5.1b)

The skewness s (resp. kurtosis k) of ψ reduces to

s[ψ] = m 3 [ψ] (m 2 [ψ]) 3/2 , k[ψ] = m 4 [ψ] (m 2 [ψ]) 2 -3, (5.1c)
where algebraic manipulations ensure that the kurtosis of the Gaussian distri-460 bution is zero. The skewness (resp. kurtosis) of q is given by

s[q] = s[ω], k[q] = k[ω].
(5.1d)

We remark from (5.1) that the skewness and kurtosis of both tracers ψ and q are not defined at boundaries, since the second moments are zero there. In addition, the eddy kinetic energy (EKE) and the eddy potential enstrophy (EPE) are provided through second order moments by:

EKE = 1 2 (m 2 [u] + m 2 [v]), EPE = 1 2 m 2 [q].
(5.1e)

In the following, theses statistics are computed for both BVE and SBVE at resolution 64 × 128, 32 × 64 and 16 × 32. Before discussing the results, the convergence of each statistic at each resolution is quantified. This can be done by progressively increasing the time interval, and computing a global error of the statistics between two adjacent intervals. More precisely, let us consider a point-wise statistic f obtained for a sufficiently long interval [t 0 , t 1 ] (where t 1 depends on the resolution considered) with a uniform partition of increment δt.

We propose to measure the convergence by a relative global error ˜ between the subintervals [t 0 , t] and [t 0 , t -δt]:

˜ (f t ) = f t -f t-δt 2 f t1 2 , (5.2)
where and t 0 is a fixed time after the spin-up (set to t 0 = 20 in this work, c.f. Figure 4) and the time increment has been fixed to δt = 0.1. A statistic is considered to be converged, as soon as the time series of relative global errors reaches a stable low error level. As shown in Figure 7, we observe that the convergence to an error less than 1% for resolutions 256 × 512, 64 × 128, 32 × 64, and 16 × 32 is reached approximatively after the time 140, 250, 350 and 500, respectively. We note that the coarser the resolution, the longer it takes to get converged statistics. This is even more pronounced for higher moments. This is likely due to higher values of the turbulent viscosity which prevent the flow to visit freely its attractor and enforce it to stay for a much longer time in the attraction basin of the equilibrium points [START_REF] Chapron | Large-scale flows Cotter[END_REF]. Note also that as observed therein, the convergence time for SBVE is shorter for all resolutions studied here (not shown). Therefore, we choose to use for all simulations the slowest convergence time (i.e. the one computed for BVE).

• 2 = •, • Ω stands for the L 2 (Ω)-
Hereafter, we focus on the comparisons of the statistics obtained for the different coarse models. To build a reference (REF) for each resolution, we directly subsample the statistics computed on the eddy-resolving data -i.e. we do not smooth them in order not to lower their energy. Figure 8 shows that at the coarsest resolution 16 × 32, the four-gyre structure is captured for both models, nonlinearity is produced such that a stronger distortion of the PV field between inner and outer gyres is observed. From Figure 9, we observe that compared to BVE, SBVE MM produces higher eddy energy in the front between the outer and inner gyres, and higher eddy enstrophy in the region between the two inner gyres. However, both coarse models do not produce enough energy flux in the western and eastern boundary layers. In particular, the too low tracers' variance in the eastern boundary layers leads to markedly higher skewness and kurtosis than those observed in the reference. Nevertheless, the introduction of randomness enables us to increase the internal variability of the tracers. For In order to provide a more quantitative comparison, we propose here a global performance index, measured by the root mean squared error (RMSE) with an aposteriori normalization to ensure a similar error level of the different statistics.

Given a statistic f with reference f REF , the normalized RMSE is defined as

RMSE(f ) = 1 |Ω| f -f REF 2 max x∈Ω |f REF (x)| . (5.3)
t ) -1 .8 -1 .4 -1 -1 -0.6 -0 .6 -0.2 -0.2 -0.2 -0 .2 0.2 0. 2 0.2 0.2 0.6 0.6 0.6 0.6 1 1 1 .4 1.8 0 0.5 1 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 -2 .2 -1 . 8 -1 . 4 -1 -0.6 -0. 6 -0.6 -0 .2 -0.2 -0.2 -0 .2 0.2 0 .2 0.2 0.2 0. 6 0.6 0.6 1 1 1 .4 1 .8 2.2 0 0.5 1 -1 .8 -1 .4 -1 -0.6 -0.6 -0.2 -0 .2 -0.2 0.2 0 .2 0.2 0.6 0.6 1 1 1 . 4 1 . 8 2 .2 0 0.5 1 -2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5
Average of potential vorticity (q t ) -0 . 6 In each panel, the first one is REF, the second one is SBVEMM and the third one is BVE.

Eddy kinetic energy (EKE) brought by properly incorporating, into large-scale simulations the effects of the small-scale flow component through its statistical distribution. From (NS for No Itô-Stokes drift), the terms related to u s are dropped in (3.8). In the second one, denoted as SBVE CP (CP for Circulation Preserving), source term dS t is removed in addition to the Itô-Stokes drift terms. This second version corresponds to the model described in [START_REF] Cotter | Numerically modeling stochastic Lie transport in fluid dynamics[END_REF]Cotter et al. ( , 2018)), for which there is no energy conservation due to the absence of the stochastic source term, c.f. [START_REF] Bauer | Deciphering the role of small-scale inhomogeneity on geophysical flow structuration: a stochastic approach[END_REF]. These two models are evaluated at the resolution 32 × 64 550 with the same parameters R β = 0.06 2 and δ 4 = 0.040 as before and the same POD noise.

Note that the advection by the Itô-Stokes drift has no effect on the resolved energy. However, as it can be observed in In addition to the discussions above, it is also important to show if the SBVE on coarse mesh can reproduce the temporal correlation behaviors of the reference [START_REF] Gugole | Numerical development and evaluation of an energy conserving conceptual stochastic climate model[END_REF]. To this end, the autocorrelation functions (ACF) for the time series of the global stream function is adopted. More precisely, this ACF is defined as

ACF(τ ) = Ψ(t) -Ψ t Ψ(t + τ ) -Ψ t t σ 2 Ψ , (5.4)
where τ stands for a time-lag, Ψ(t) = 1 |Ω| Ω ψ(x, t)dx is the global stream function at time t, and σ Ψ is the (temporal) standard deviation of Ψ. Figure 14 shows that compared to the BVE at each coarse resolution, both SBVE POD and SBVE MM capture better the ACF of the reference. For instance, they have smaller decorrelation time scales compared to the BVE. Besides, the best results are provided by the mode matching method, which is consistent with our previous conclusions.

Conclusions

The approach explored in this work consists in a stochastic representation of (2017). The encouraging results presented here inspire us to implement the proposed stochastic approach on more complex flow, and to test more physical parameterizations for the small-scale random flow. Two subsequent projects on the study of Q-GCM [START_REF] Hogg | A quasigeostrophic coupled model (Q-GCM)[END_REF] and NEMO (NEMO team, 2016) are already in progress. In particular, we aim to parametrize the noise on the isopycnal surfaces [START_REF] Gent | Isopycnal mixing in ocean circulation models[END_REF], such that the transfer of the available potential energy to the resolved kinetic energy can be efficiently achieved. Successes of other stochastic parameterizations [START_REF] Grooms | Stochastic superparameterization in a quasigeostrophic model of the Antarctic Circumpolar Current[END_REF][START_REF] Gugole | Numerical development and evaluation of an energy conserving conceptual stochastic climate model[END_REF][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] provide some confidence that the backscatter and jet enforcement can be reproduced with some success under the proposed stochastic framework. Besides, the 620 feasibility for the application of the mode matching strategy on these models will be analyzed. All these efforts aim at progressively going toward the study of data-driven stochastic IPCC-class climate models and to confirm that relevant stochastic flow models contribute to improve them.

v •,0 = v •,N = 0, v 0,• = v 1,• , v M +1,• = v M,• . (6.2c)
The same discrete representations apply to u s and σdB t .

Discretized differential operators can now be consistently built within such a specific staggered grid. They are based on the following first derivative ap-660 proximations:

∂ h x θ
i+ 1 2

,j = θ i+1,j -θ i,j ∆x , ∂ h y θ i,j+ 1 2 = θ i,j+1 -θ i,j ∆y , (6.3) which remains valid whether θ be defined in p-, u-and v-grid. As such, the velocity u (resp. for σdB t ) can be derived from a given stream function ψ (resp.

for ϕdB t ), by applying the discretized perpendicular gradient ∇ ⊥ h = [-∂ h y , ∂ h x ] T . Subsequently, the vorticity is given by ω = (∇ ⊥ h ) T u = (∇ ⊥ h ) T ∇ ⊥ h ψ = ∇ 2 h ψ, where the discretized Laplacian operator reads: ω i,j = ψ i+1,j -2ψ i,j + ψ i-1,j ∆x 2 + ψ i,j+1 -2ψ i,j + ψ i,j-1 ∆y 2 , (6.4)

for any interior (i, j) points of the p-grid. Conversely, the stream function ψ

(and similarly for ϕdB t ) can be re-constructed from a current vorticity ω and the inverse of the Laplacian operator (∇ 2 h ) -1 , expressed in practice in the Fourier domain through an efficient discrete Fourier transform solver [START_REF] Press | Numerical recipes 3rd edition: The art of scientific computing[END_REF]. To discretize the nonlinear jacobian terms in the BVE and in the SBVE, we employ Arakawa's 9-points conservative scheme (Arakawa and Lamb, 1981).

Actually, such a discretization can be interpreted through interpolated derivatives on the staggered grid. For instance, the advection of ω by ψ can be written as 680

J h (ψ, ω) = 1 3 ∂ h x ψ x ∂ h y ω y -∂ h y ψ y ∂ h x ω x + 1 3 ∂ h x (ψ∂ h y ω y ) x -∂ h y (ψ∂ h x ω x ) y + 1 3 ∂ h y (ω∂ h x ψ x ) y -∂ h x (ω∂ h y ψ y )
x , (6.6a)

where (θ x ) i+1/2,j = (θ i+1,j + θ i,j )/2 and (θ y ) i,j+1/2 = (θ i,j+1 + θ i,j )/2 stand for central interpolations between two neighboring points in x-and y-directions respectively. Such a discretized operator is applied in the very same way on the other advection terms of the SBVE (3.8) associated to ϕdB t . The source terms (3.8b) are otherwise discretized as J h (U y , u y ) + J h (V x , v x ), (6.6b)

where U and V denote the two components of U = σdB t -u s dt.

The stochastic RK3 scheme of [START_REF] Cotter | Numerically modeling stochastic Lie transport in fluid dynamics[END_REF] is given by ω (1) = ω n + f ∆tψ n , ω n + g ϕdB n t , ∆tu n s , ω n + h σdB n t , ∆tu n s , u n , (6.7a) 1) +f ∆tψ (1) , ω (1) +g ϕdB n t , ∆tu n s , ω (1) +h σdB n t , ∆tu n s , u (1) , (6.7b) 2) +f ∆tψ (2) , ω (2) +g ϕdB n t , ∆tu n s , ω (2) +h σdB n t , ∆tu n s , u (2) , (6.7c) where f (∆tψ, ω) = -∆tJ(ψ, ω)+ ∆t R β (F +D-∂ψ ∂x ), g(ϕdB t , ∆tu s , ω) = -J(ϕdB t , ω)+ ∆tu s •∇ ω -1 R β ∂ ∂x ϕdB t and h(σdB t , ∆tu s , u) = -J(σdB t -∆tu s , u).

ω (2) = 3 4 ω n + 1 4 ω (
ω n+1 = 1 3 ω n + 2 3 ω (
We assume that there exists an additional isotropic random field living at the unresolved (sub-grid) scales -i.e. not represented at the considered resolution scale. Thus, the global contribution of its variance tensor a 0 I 2 to the enstrophy dissipation can be expressed by Arakawa, A., Lamb, V.R., 1981. A potential enstrophy and energy conserving

1 2 Ω ω o a 0 ∇ 2 ω o dx = - a 0 2 Ω ∇ω o 2 dx, ( 6

  .8b) where d t ω = ω t+dt/2 -ω t-dt/2 stands for the central time-increment and where the prime symbols have been dropped for all the adimensional variables. Note that compared to the explicit Itô form (3.6), the diffusive terms are now implicit in the Stratonovich time-integral. Switching from Itô to Stratonovich integral allows us to benefit from the advantages of both integral representations: the Itô flow (3.1) allows us to keep a zero mean noise term (whereas it is not true for Stratonovich convention) and provides a way to explain more easily the different physical contributions of the noise terms. The Stratonovich representation permits the use of the classical chain rule differentiation and leads to more efficient numerical implementation, in which the diffusive contribution is implicitly taken into account (Cotter et al., 2019). The advantages and limitations of Itô and Stratonovich formulations in the context of fluid flow dynamics together with their relationship are detailed in Bauer et al. (2020).

Figure 1 :

 1 Figure 1: Illustration of the spatial modes truncation for the random velocity, within the spectrum of the corresponding eigenvalues.

Figure 2 :

 2 Figure 2: Illustration of the time-averaged kinetic energy spectrums. The wavenumber κ 0 is searched as the first point where the observation and the coarse-simulation BVE derivate, in order to estimate γ 0 from (4.3c).

Figure 3 :

 3 Figure 3: Illustration of mode matching principle: Selection of a sample set of time based on (4.5b) corresponding to potential events matching the large-scale configuration of the current simulation.

Figure 4 :

 4 Figure 4: Instantaneous snapshots of PV and time series of the global energy and enstrophy, provided by the eddy-resolving BVE at resolution 256 × 512. The global energy is defined by

Figure 5 :

 5 Figure 5: Instantaneous snapshots of PV provided by different models at resolution 64 × 128.

Figure 6 :

 6 Figure 6: Instantaneous snapshots of the small-scale random stream functions ψ POD r , ψ MMr

Figure 7 :

 7 Figure 7: Time series of the relative errors of the statistics by progressively increasing the time interval. In each row, the left plot shows the statistical errors of the stream function (or velocity), and the right one shows that of the PV. In each column, the results correspond, from top to bottom, to resolutions 256 × 512, 64 × 128, 32 × 64, and 16 × 32. Note that in both cases, the first (adimensioned time) interval on which we compute the statistics is set to be [20, 80]; this interval is progressively augmented with a time step of 0.1. The Y -axis values describe the converging percentage of one statistic w.r.t. its global (over the spatial domain) value performed at previous instant.

  instance, as shown in Figures 10 and 11, the region with extreme values of 510 skewness and kurtosis is significantly reduced for SBVE MM when compared to BVE. As the resolution increases, it can be noticed from Figures 12 and 13 that the local structures of the PV statistics provided by SBVE MM , qualitatively converges to the reference.

Figure 8 :

 8 Figure 8: Contour plots of the time-average fields at resolution 16 × 32. The top three plots depict the SF with contour interval (CI) of 0.2, and the bottom three show PV with CI of 0.1.

Figure 9 :Figure 10 :Figure 11 :

 91011 Figure 9: Contour plots of the time-variance fields at resolution 16 × 32. The top three plots depict EKE with CI of 30, and the bottom three show EPE with CI of 0.025. In each panel, the first one is REF, the second one is SBVEMM and the third one is BVE.

Figure 14 :

 14 Figure 14: Comparison of the autocorrelation functions (ACF) of the global stream function between different models, at resolution 16 × 32, 32 × 64 and 64 × 128. All the ACFs are calculated from t = 20 to t = 100 using 8001 snapshots.

  mesoscale eddy effects on large-scale ocean circulation. The main result demonstrates that the large-scale flow can be simulated by a coarse-resolution model composed of a multiplicative random forcing, a heterogeneous diffusion and an advection correction. All these ingredients allow us to correctly backscatter, dissipate and distribute the large-scale energy. Such a random model, built from 590 classical conservation laws, provides here an explicit eddy representation for a single-layered QG model. Under this regime, additional vorticity sources arise from the interaction of the strains between the small-scale random component and the large-scale current. These terms are important in conserving the global energy of the resolved scales. Numerically, the spatial correlation of the random fields in the coarse model is defined from the coherent structures of an eddy-resolving simulation. In order to quantify the accuracy of the proposed random model, a statistical analysis of the flow tracers has been performed. As expected, compared to a classical coarse model, the proposed stochastic model better represents the nonlinearity at the 600 resolved scales while properly dissipating the unresolved scales, leading hence to a balanced correction of excessive dissipation and the continuous increase of internal variability. As a result, it reproduces better on a coarse mesh, the local structures of the distribution of eddy-resolving tracers. Further analysis showed that the vorticity sources are important in locally strengthening the eddy flux of PV. Although the idealized barotropic model used in this work cannot describe quantitatively the real ocean, they do in fact produce qualitatively realistic patterns of large-scale flow in the major basins of the world, as illustrated in Vallis

  vertical cell interfaces respectively (they are referred to as u-points and v-points respectively). Considering M × N cells, then the p-grid has (M + 1) × (N + 1) points, with homogeneous Dirichlet boundary values defined asψ 0,• = ψ M,• = ψ •,0 = ψ •,N = 0. (6.2a)The same boundary condition is imposed to ω and ϕdB t . The u-grid has a dimension of (M + 1) × (N + 2) points together with free-slip boundary values imposed asu 0,• = u M,• = 0, u •,0 = u •,1 , u •,N +1 = u •,N ,(6.2b)while the v-grid has (M +2)×(N +1) points with the free-slip boundary values:

670Figure 15 :

 15 Figure 15: Illustration of the staggered grid.

Figure 16 :

 16 Figure 16: Estimated values of the biharmonic boundary layer size δ4 at different coarse resolutions, using the eddy-resolving data.

  norm, and f t (x), ∀x ∈ Ω, denotes the local-in-time point-wise statistics associated to the interval [t 0 , t]. In practice, we initiate this procedure from a reasonable intermediate instant t c ∈ [t 0 , t 1 ],

Table 1

 1 compares the results of the different models at the coarsest resolution 16 × 32. The proposed stochastic model shows a clear improvement of all the 520 statistics w.r.t. the references. This improvement is particularly noticeable for the higher moments. For instance, compared to BVE, SBVE MM has 35.87% and 39.26% less errors in skewness and kurtosis of the stream function (SF), respectively. The mode matching strategy, SBVE MM , performs better than the POD strategy, SBVE POD , for all moments, although the latter already reduces the BVE error of the first and second moments (with an improvement of 9, 7% for the SF mean and 12, 6% for EKE). Both SBVE MM and SBVE POD reach very similar errors in terms of EKE and EPE (with an improvement above

	Average of stream function (ψ

Table 1 :

 1 Table2and 3 we see that these RMSEs improvements still hold as the resolution is increased. The improvements at resolution 64 × 128 in terms of EKE and EPE are still noticeable (25%). The third order moment of SF continues to improve (20%) while for the fourth order moments the improvement is less significant.Both SBVE MM and SBVE POD improve also the first order moments at resolution 32 × 64 (at almost the same rate as for the coarsest resolution) and 64 × 128 (with a smaller decreasing of the errors). This latter has by definition a lower noise level. This illustrates that even for weak noise levels the stochastic systems 540 lead to better results than the deterministic version. Comparison of the normalized RMSEs between different models at resolution 16 × 32 with R β = 0.06 2 and δ 4 = 0.049 fixed. The lowest errors are highlighted in bold.

	Model	RMSE	ψ	t	q t	EKE	EPE	s[ψ]	s[q]	k[ψ]	k[q]
	BVE	0.245	0.091	0.111	0.148	0.499	0.406	0.782	0.806
	SBVE POD	0.221	0.082 0.097 0.132 0.489	0.390	0.624	0.758
	SBVE MM	0.197 0.075 0.098 0.131 0.320 0.325 0.475 0.631
	Model	RMSE	ψ	t	q t	EKE	EPE	s[ψ]	s[q]	k[ψ]	k[q]
	BVE	0.108	0.061	0.073	0.122	0.190	0.166	0.218	0.155
	SBVE POD	0.094 0.056 0.064	0.116 0.161 0.146 0.182 0.122
	SBVE MM	0.089 0.055 0.058 0.107 0.161 0.136 0.181 0.106

Table 2 :

 2 Comparison of the normalized RMSEs between different models at resolution 32 × 64 with R β = 0.06 2 and δ 4 = 0.040 fixed. The lowest errors are highlighted in bold.We analyze now the individual effects of the Itô-Stokes drift and the additional vorticity sources on the accuracy of the statistics. To that end, two particular versions of SBVE have been run. In the first one, denoted as SBVE NS

	Model	RMSE	ψ	t	q t	EKE	EPE	s[ψ]	s[q]	k[ψ]	k[q]
	BVE	0.075	0.028	0.036	0.055	0.087	0.039	0.068	0.035
	SBVE POD	0.073 0.024 0.034	0.047	0.080	0.036 0.061 0.031
	SBVE MM	0.069 0.023 0.027 0.041 0.068 0.034 0.061 0.029

Table 3 :

 3 Comparison of the normalized RMSEs between different models at resolution 64×128 with R β = 0.06 2 and δ 4 = 0.026 fixed. The lowest errors are highlighted in bold.

  Table4from the comparison between SBVE NS and SBVE POD , its inclusion improves all the SF moments (with a marked decrease of errors in the first and third moments). These improvements outlines the importance of taking into account properly the inhomogeneity of the small-scale component as captured by the Itô-Stokes drift. In contrast, the Itô-Stokes drift plays no role in terms of the PV mean and EPE (and only leads to marginally better third and fourth order PV moments). Finally SBVE CP

		order PV moments. For the barotropic regime studied here, global energy con-
		servation together with the Itô-Stokes correction as considered in this stochastic
	570	framework, provides more accurate long-term statistics, than models in which
		these two features are not taken into account.
		Model	RMSE	ψ	t	q t	EKE	EPE	s[ψ]	s[q]	k[ψ]	k[q]
		BVE	0.108	0.061	0.073	0.122	0.190	0.166	0.218	0.155
		SBVE POD	0.094 0.056 0.064 0.116 0.161 0.146 0.182 0.122
		SBVE NS POD	0.100 0.056 0.067 0.116 0.185	0.148	0.191	0.130
		SBVE CP POD	0.104 0.056 0.068 0.115 0.185	0.156	0.208	0.138
	560
		provides intermediate results between the traditional large-scale deterministic
		model and the proposed stochastic model.
		The influence of the stochastic source term can be appreciated comparing
		SBVE NS and SBVE CP . This term, which guarantees the conservation of the
		global energy, enables us to improve the SF mean and kurtosis as well as the
		PV skewness and kurtosis. The association of both the Itô-Stokes drift and the
		stochastic source terms improves the four SF moments and the third and fourth

Table 4 :

 4 Comparison of the normalized RMSEs between the proposed stochastic model SBVEPOD, a version without the Itô-Stokes drift (SBVE NS POD ) and a circulation preserving version without both the Itô-Stokes drift and the stochastic source term (SBVE CP POD ) at resolution 32 × 64 with R β = 0.06 2 and δ 4 = 0.040 fixed. The lowest errors are highlighted in bold.

tion, in the absence of forcing, the resulting model preserves potential vorticity while model (3.8) conserves the global energy. We will see, however, that model (3.8) enables us to reproduce more accurately potential vorticity and enstrophy statistics, highlighting in this setting the importance of energy conservation.

defined empirically. For instance, we fix it as γ 1 = 0.999 in this work.

Acknowledgements

The authors acknowledge the support of the "Laboratoires d'Excellence" CominLabs, Lebesgue and Mer through the SEACS project and of the ERC EU project 856408 -STUOD.

Average of potential vorticity (q t ) Eddy potential enstrophy (EPE) 10% for both quantities) and SBVE MM is more efficient in reducing errors in the third and fourth moments. These results highlight the benefits that are Appendix A. Snapshot POD This section describes briefly the snapshot POD method [START_REF] Sirovich | Turbulence and the dynamics of coherent structures, part I: Coherent structures[END_REF].

The eigenvalues and their associated eigenfunctions can be estimated from the following eigenvalues problem:

where Λ = (λ i ) i=1,...,Nt is the set of decaying eigenvalues, i.e.

eigenvectors. The temporal modes {b i } i=1,...,Nt are then normalized such that (6.1c) where δ ij denotes for the Kronecker symbol here. The spatial modes {φ i } i=1,...,Nt

given by

are orthonormal:

And, from this spectral decomposition, each snapshot can be reconstructed by