Takaya Uchida

Bruno Deremble
email: brunoderemble@cnrs.fr

Stephane Popinet

Deterministic model of the eddy dynamics for a midlatitude ocean model

studies in which we aim to improve the large-scale state by parameterizing the net mesoscale feedback onto the large-scale flow.

), Juricke et al. (2019) and Perezhogin (2019) implemented a prognostic equation for the sub-grid (unresolved) eddy energy and achieve the backscattering via a negative viscosity. One notable difference in our method is that while many previous studies have formulated their parameterizations based on a local closure (i.e. relating the eddy momentum/buoyancy flux locally at each grid point to the resolved momentum/buoyancy), we construct our closure by incorporating basin-scale information. This is motivated by the fact that Venaille et al. (2011) and Grooms et al. (2013) have shown that the eddy feedback on the large-scale flow is strongly non-local. We also focus on the sub-grid potential vorticity (PV) equation rather than sub-grid energy within the quasi-geostrophic (QG) framework. The QG framework has been shown to be fruitful in examining the eddy-mean flow interaction and formulating eddy closures (e.g. Marshall et al. 2012; Porta Mana and Zanna 2014; Mak et al. 2016; Berloff 2018). In particular, Berloff et al. (2021) have shown some success in accounting for the non-local eddy feedback by solving for the sub-grid QGPV equation (cf. Eq. 22).

While our approach is similar, here, we propose an alternative strategy to achieve a deterministic closure for the sub-grid PV. This approach of prognostically solving for the sub-grid dynamics is sometimes referred to as super parameterization and has been commonly implemented for atmospheric or oceanic convection (e.g.

ABSTRACT: Mesoscale eddies, although being on the scales of (20-100 km), have a disproportionate role in shaping the mean stratification, which varies on the scale of (1000 km). With the increase in computational power, we are now able to partially resolve the eddies in basin-scale and global ocean simulations, a model resolution often referred to as mesoscale permitting. It is well-known, however, that due to grid-scale numerical viscosity, mesoscale-permitting simulations have less energetic eddies and consequently weaker eddy feedback onto the mean flow. In this study, we run a quasi-geostrophic model at mesoscale-resolving resolution in a double gyre configuration and formulate a deterministic closure for the eddy rectification term of potential vorticity (PV), namely, the eddy PV flux divergence. Our closure successfully reproduces the spatial patterns and magnitude of eddy kinetic and potential energy diagnosed from the mesoscale-resolving model.

One novel point about our approach is that we account for non-local eddy feedbacks onto the mean flow by solving the 'sub-grid' eddy PV equation prognostically in addition to the mean PV.

Introduction

In the field of fluid dynamics and turbulence, formulating a closure for the governing equations has been a long standing problem [START_REF] Smagorinsky | The geography of linear baroclinic instability in earth's oceans[END_REF][START_REF] Launder | Progress in the development of a reynolds-stress turbulence closure[END_REF]. Resolving the flow down to the molecular scale where kinetic energy is dissipated to internal energy due to molecular viscosity is usually not feasible, whether in observations or a numerical model. Particularly in the field of geophysical fluid dynamics (GFD) where the scales of interest span up to (1000 km), resolving the molecular scale is practically unachievable and will remain so for the foreseeable future. Due to the lack of resolution, a numerical model will only solve the governing equations for the "resolved" field, and some work has to be done to account for the "unresolved" field. A significant effort in GFD has been, therefore, to formulate a closure for the unresolved field, i.e. represent the unresolved field with the resolved momentum and/or tracer field (e.g. [START_REF] Mellor | Development of a turbulence closure model for geophysical fluid problems[END_REF][START_REF] Redi | Oceanic isopycnal mixing by coordinate rotation[END_REF][START_REF] Gent | Isopycnal mixing in ocean circulation models[END_REF][START_REF] Bachman | A scale-aware subgrid model for quasigeostrophic turbulence[END_REF].

The ocean component of climate models suffers from this issue of missing the unresolved dynamics because it hardly resolves the mesoscale eddies (horizontal scale of (20-100 km)). This is problematic because the unresolved (small-scale) field not only drains energy from the resolved (large-scale) field but also partially feeds back onto the resolved field via up-scale momentum and buoyancy fluxes, and so modifies the dynamics of the large-scale flow [START_REF] Vallis | Atmospheric and Oceanic Fluid Dynamics[END_REF][START_REF] Lévy | Large-scale impacts of submesoscale dynamics on phytoplankton: Local and remote effects[END_REF][START_REF] Arbic | On eddy viscosity, energy cascades, and the horizontal resolution of gridded satellite altimeter products[END_REF][START_REF] Aluie | Mapping the energy cascade in the north atlantic ocean: The coarse-graining approach[END_REF][START_REF] Ajayi | Diagnosing cross-scale kinetic energy exchanges from two submesoscale permitting ocean models[END_REF]. Modelling studies with varying spatial resolution have shown that only partially resolving the mesoscale results in weaker mesoscale eddies, and consequently weaker feedback onto large-scale flows. It is also well-known that mesoscale eddies exert a strong influence on oceanic jets such as the Gulf Stream (Chassignet and Xu 2017;[START_REF] Kjellsson | The impact of horizontal resolution on energy transfers in global ocean models[END_REF][START_REF] Chassignet | On the importance of high-resolution in large-scale ocean models[END_REF]. Considering the impact of the jets on global tracer transport and air-sea interaction [START_REF] Kelly | Western boundary currents and frontal air-sea interaction: Gulf Stream and Kuroshio Extension[END_REF][START_REF] Tréguier | Meridional transport of salt in the global ocean from an eddy-resolving model[END_REF][START_REF] Jones | Components of upper-ocean salt transport by the gyres and the meridional overturning circulation[END_REF][START_REF] Bellucci | Air-sea interaction over the gulf stream in an ensemble of highresmip present climate simulations[END_REF], improving the representation of the eddy feedback onto the jet has climate implications. Hence, there has been a growing effort to represent the inverse cascade of kinetic energy otherwise lost to grid-scale numerical viscosity at mesoscale-permitting resolution, a process often referred to as energy backscattering parameterizations (e.g. [START_REF] Zanna | Scale-aware deterministic and stochastic parametrizations of eddy-mean flow interaction[END_REF][START_REF] Berloff | Dynamically consistent parameterization of mesoscale eddies. part iii: Deterministic approach[END_REF][START_REF] Jansen | Toward an energetically consistent, resolution aware parameterization of ocean mesoscale eddies[END_REF][START_REF] Bachman | The gm+ e closure: A framework for coupling backscatter with the gent and mcwilliams parameterization[END_REF][START_REF] Josey | Wind stress forcing of the ocean in the SOC climatology: Comparisons with the NCEP NCAR, ECMWF, UWM/COADS, and Hellerman and Rosenstein datasets[END_REF][START_REF] Perezhogin | Deterministic and stochastic parameterizations of kinetic energy backscatter in the nemo ocean model in double-gyre configuration[END_REF]Zanna and Bolton 2020, and references therein). Our study here is in the same realm of parameterization

Model and methods

a. The control run

We adopt the QG framework in order to describe the well-known double gyre circulation in an idealized midlatitude ocean basin. This model is known to capture both the large-scale and smallscale variability of the ocean with a relatively coarse vertical resolution (cf. [START_REF] Berloff | Dynamically consistent parameterization of mesoscale eddies. part i: Simple model[END_REF]. The QG formalism is meant to describe dynamical regimes for a prescribed background stratification 2 and Coriolis parameter . Two ingredients are necessary to reproduce the double gyre pattern: the planetary vorticity must vary with latitude for the western boundary intensification and the wind forcing must be cyclonic in the northern part of the domain and anticyclonic in the southern part of the domain. In order to satisfy the first condition, we work with the -plane approximation such that the Coriolis parameter varies linearly with latitude. This sets the planetary scale = 0 / which is large compared to the deformation scale = / 0 , (with the depth of the ocean and 0 the average value of the Coriolis parameter in the domain). In this formalism, the main dynamical variable is the QG potential vorticity (PV) defined as

= ∇ 2 + Γ def = L , (1)
with the stream function, ∇ 2 the horizontal Laplace operator and

Γ def = 2 0 2 (2)
the vertical stretching operator. The horizontal velocity is defined as

= - and = , (3)
and the buoyancy is defined as

= 0 . ( 4
)
The equation of evolution of PV is

+ ( , ) + = 4 ∇ 4 + ∇ 2 + , (5)
( , ) = - , (6)
the Jacobian operator, which corresponds to the non-linear advective term, 4 the bi-harmonic viscosity, the bottom friction coefficient which parameterizes a bottom Ekman layer (and is thus non-zero in the lower layer only), and the forcing resulting from an Ekman pumping in a thin Ekman layer at the surface and is thus non-zero in the upper layer only. We build the numerical version of this model in the Basilisk framework (Popinet 2015, basilisk.fr).

We solve Eqs. ( 5) and ( 1) in a horizontal square domain with side = 5000 km and of vertical extension = 5000 m. We discretize these equations with 512 × 512 horizontal points (which correspond to a horizontal resolution of slightly less than 10 km) and 4 vertical layers of thickness ℎ 1 = 238 m, ℎ 2 = 476 m, ℎ 3 = 953 m and ℎ 4 = 3333 m (from top to bottom). We adjust the background stratification 2 to mimic the stratification in middle of the subtropical gyre in the North Atlantic such that at each layer interface, we have 2

1.5 = 1.7 × 10 -5 s -2 , 2 2.5 = 1.1 × 10 -5 s -2 , 2 3.5 = 3.2 × 10 -7 s -2
, from top to bottom. The average value of the Coriolis parameter is 0 = 9.3 × 10 -5 s -1 and = 1.7 × 10 -11 m -1 s -1 . For these parameters, the three deformation radii are 1 = 25 km, 2 = 10 km and 3 = 7 km. Note that these deformation radii correspond to the inverse squared eigenvalue of the vertical stretching operator. At this resolution we choose 4 = 6.25 × 10 9 m 4 s -1 , and a spindown timescale = 1/166 days (which corresponds to a bottom Ekman layer of thickness = 7.5 m). We solve the elliptic equation (Eq. 1) with homogeneous Dirichlet boundary conditions = 0 on the sides (corresponding to no flux boundary condition) and homogeneous Neumann boundary conditions = = 0 at the top and bottom boundaries.

The forcing is

= ∇ × τ 0 ℎ 1 , with = 0 sin 3 . ( 7
)
We use a a cubic sine function in the definition of the wind in order to reproduce a narrow midlatitude atmospheric jet. For such a narrow jet, the boundary between the positive and negative area of the wind stress curl pattern is sharper than if we use the traditional cosine shape for the wind pattern. We choose 0 = 0.25 N m -2 which is an acceptable value for the difference between the maximum and minimum value of the wind in the North Atlantic [START_REF] Josey | Wind stress forcing of the ocean in the SOC climatology: Comparisons with the NCEP NCAR, ECMWF, UWM/COADS, and Hellerman and Rosenstein datasets[END_REF]. We have 6 Accepted for publication in Journal of Physical Oceanography. DOI 10.1175/JPO-D-21-0217.1.

Unauthenticated | Downloaded 04/07/22 08:27 AM UTC also kept the wind stress axisymmetric as our interest is on eddy time scales and not low-frequency variability [START_REF] Berloff | The turbulent oscillator: A mechanism of lowfrequency variability of the wind-driven ocean gyres[END_REF].

To integrate the model in time, we first perform a spin up phase of 80 years at low resolution (Δ = 78.13 km) followed by another 80 years at the prescribed resolution (Δ = 9.77 km). After this spin up of 160 years in total, the model is in a statistically steady state (i.e.

≈ 0). We show in Fig. 1, the meridional profile of the wind stress and snapshot of the local Rossby number (i.e.

relative vorticity normalized by 0 ). Except in the region of the separated jet, the local Rossby number is much smaller than unity, consistent with the QG scaling. Henceforth, we refer to this run as the CTRL run. 7), however this is only cosmetic because this additional term does not impact the wind stress curl, which is what ultimately matters in QG dynamics.

b. Mean flow and eddy models

In order to set up the framework for parameterization, we perform a Reynolds decomposition of each dynamical variable as the sum of its mean (denoted with an overbar) and a perturbation about the mean (denoted with a prime) as shown here for the stream function

= + . (8)
We leave the definition of the 'mean' intentionally vague for now to keep the arguments general.

If we use this decomposition in the equation of evolution of PV, we get

( + ) + ( + , + ) + ( + ) = 4 ∇ 4 ( + ) + ∇ 2 ( + ) + , (9)
and taking the mean of this equation gives

+ ( , ) + = 4 ∇ 4 + ∇ 2 -( , ) + . ( 10
)
The term ( , ) = ∇ • (u ) is known as the eddy rectification of the large-scale flow. At this point, it is common in eddy parameterization studies to reinterpret the 'mean' flow as the resolved flow of a coarse resolution model and formulate a closure to model the contributions from the sub-grid flow onto the resolved flow. Namely, the eddy rectification is the sub-grid feedback that many studies seek to parameterize (e.g. [START_REF] Eden | Parameterising meso-scale eddy momentum fluxes based on potential vorticity mixing and a gauge term[END_REF][START_REF] Marshall | A framework for parameterizing eddy potential vorticity fluxes[END_REF][START_REF] Mana | Toward a stochastic parameterization of ocean mesoscale eddies[END_REF][START_REF] Mak | A new gauge-invariant method for diagnosing eddy diffusivities[END_REF][START_REF] Berloff | Dynamically consistent parameterization of mesoscale eddies. part iii: Deterministic approach[END_REF]). The reinterpretation is based on the expectation that the reduction in variability resulting from the averaging operator would mimic the partially resolved variability at mesoscale-permitting resolution.

A well-known approach to parameterize the role of the eddies on the mean flow is to approximate the eddy PV flux by a term that is proportional to the local gradient of the mean PV [START_REF] Gent | Isopycnal mixing in ocean circulation models[END_REF]:

u ≈ GM ∇(Γ ) (11)
with GM an eddy diffusivity coefficient and Γ , the vertical stretching component of PV. In effect, this term corresponds to the diffusion of the thickness of an isopycnal layer or a skew diffusion (Griffies 1998). Hence, the GM parameterization belongs to the class of "down-gradient" parameterizations and its effect is always to flatten isopycnal surfaces.

In the present study, we want to overcome the down-gradient parameterization and we are going to explicitly model the (sub-grid) eddy dynamics with an independent model in order to formulate a parameterization for the eddy rectification. The equation for the (sub-grid) eddy dynamics can be obtained by taking the difference between Eq. ( 9) and Eq. ( 10), given by

+ ( , ) + ( , ) + ( , ) + = 4 ∇ 4 + ∇ 2 + ( , ) + . (12)
Note that the presence of ( , ) in this equation is somewhat cumbersome because in order to simulate the eddy equation, which we propose as the independent model, it requires a priori knowledge of the eddy rectification term (viz. the mean properties of eddy-eddy interaction) as a forcing which renders the eddy model meaningless. As we shall see in section 3a, if we run the eddy model without this term, we get a poor representation of the eddy field. The crux of this paper is a proposition to parameterize this term with a modification to the definition of the mean (see section 3).

c. Mean flow and eddy dynamics in the full model

We first analyze the output of reference run (CTRL) which is a mesoscale-resolving simulation.

We recall that this model solves the full PV equation (Eq. 5): we decompose the output of that simulation into a mean and an eddy flow. We perform this decomposition with a time mean. For the remainder of this study, the averaging operator (•) is defined as a time mean; note that because the forcing is stationary, = and = 0, and so the time mean is similar to an ensemble mean here under the ergodic assumption [START_REF] Galanti | Is turbulence ergodic?[END_REF]). We will consider the mean and eddy flow diagnosed from the CTRL run as the "truth". We will then use these diagnostics to validate the model of the eddy dynamics only (section 3).

The stream function of the full model exhibits a standard double gyre pattern with an strong eddying jet that separate the cyclonic and anticylconic gyres. Such pattern has already been observed and described in numerous studies; we wish however to highlight the mean/eddy decomposition from an energetic perspective. In quasi geostrophy, the total energy is the sum of potential energy

= 1 2 2 2 , (13)
and kinetic energy

= 1 2
( 2 + 2 ) , (14)
and since potential and kinetic energies are quadratic quantities, we write their time average as

= 1 2 2 2 + 1 2 2 2 def = + , (15)
= 1 2
( 2 + 2 ) + 1 2 ( 2 + 2 ) def = + , (16)
with and the potential and kinetic energy of the time mean flow and and the mean potential and kinetic energy of the eddy flow.

We plot in Fig. 2a, a snapshot of the eddy kinetic energy in the upper layer. We find at least two distinct dynamical regimes: (i) the eddying jet with on the order of 0.5 m 2 s -2 (corresponding to a velocity of | | ∼ 1 m s -1 ). The intensity of the jet decreases downstream (eastward). (ii) A region with moderate eddies in the middle of each gyre; the magnitude of these eddies increases from East to West but their overall intensity is order ∼ 0.04 m 2 s -2 (| | ∼ 0.2 m s -1 ). There are other dynamical regions such as quiescent zone with no eddies at all at the same latitude as the jet but near the eastern boundary, and the regions near the northern and southern boundaries.

We plot with the same colorbar the eddy potential energy for the same snapshot (Fig. 2b). We observe that the magnitude of is similar to the magnitude of consistent with the QG scaling. We plot in Figs. 2c and 2d the mean eddy kinetic energy and mean eddy potential energy.

The eddy potential energy and eddy kinetic energy exhibit similar patterns and are maximal in the jet. The maximum value of eddy energy in the jet area reflects the meandering jet. These meanders are strongest near the western boundary and decrease in amplitude moving east.

The energy stored in the mean flow exhibits a radically different pattern than the eddy energy (Figs. 2e and2f). The QG model exhibits the standard result that most of the large-scale energy is stored in the form of potential energy. Note that the colorbar in Fig. 2f

d. Vorticity balance of the mean flow

For sufficiently long integration, the first term in the mean flow (Eq. 10) will eventually vanish.

There is thus a balance between the remaining terms of the mean PV equation. We only focus here on the rectification term that we plot in Fig. 3. We plot in Fig. 3a the raw estimate of ( , ) computed with 500 independent snapshots that are 60 days apart (which corresponds to the eddy decorrelation time scale, not shown). And we plot in Fig. 3b the smoothed version where we average 16 neighboring grid points and linearly interpolate back on the fine grid for visualization purposes. From the latter plot, a large-scale component of this field emerges in the return flow area.

The region of the separated jet exhibits a stronger signal whereas the region near the boundaries also exhibit intense magnitude signal. We emphasize one more time that ( , ) is very important to establish the flow pattern that we described earlier: this term is of the same order of magnitude as the other terms in Eq. ( 10) and if it were absent, the mean flow would be quite different.

It is also important to note that the pattern in Fig. 3a clearly has not converged because when we sum all the terms in Eq. ( 10), viz. ( , ) + ( , ) + -4 ∇ 4 -∇ 2 -, we get a noisy field (similar to Fig. 3a), whereas we should actually get zero everywhere if the model were run long enough ( ∼ 0; not shown). With the purpose of formulating a deterministic model for the eddy rectification term, some spatial smoothing is appropriate in order to filter out stochastic variability. If we admit that the smoothed ( , ) is the deterministic part and that ( , ) should converge towards its smoothed version, we can estimate the number of samples we need for convergence with a maximum of 10% error. Indeed the standard error of the mean is given by /

√

where is the standard deviation of the time series at a given point and the number of samples. If we want the error bar to be 10% of the value of the mean , the 95% confidence interval on the mean for that tolerance is given by 0.1 = 2 / √ such that = 400 2 / 2 . We get an estimate of = 10 5 samples to get this 10% precision for the mean. This corresponds to 10 4 years of simulation which is clearly out of reach in the current setup. We have tested this using the 2740 years of output from [START_REF] Kondrashov | Stochastic modeling of decadal variability in ocean gyres[END_REF] and found the convergence to be very slow (personal communication with Pavel Berloff). The fact that such a long integration is required for accurate statistics is problematic from an eddy closure perspective, namely the eddy statistics of today would depend on the dynamical state of the system thousands of years in the past. This conundrum also highlights the need for a closure for the eddy rectification.

The sub-grid PV model

Our goal is now to see if we can approximate ( , ) with a dynamical equation for the perturbed quantities. Given a mean flow, the eddy or sub-grid PV equation (Eq. 12) can be prognostically solved with the caveat of the 'unknown' eddy rectification forcing ( , ) which appears on the right-hand side of Eq. ( 12). We are going to test two strategies to handle this term: (i) we will simply remove it, and (ii) we are going to propose a spatial filter approach. It is also important to note that the eddy dynamics is driven only by the presence of the barred variables in the eddy equation. In this section, we take the time mean field of CTRL run for these barred variables. With this choice, we will test now if the eddy model is able to reproduce the eddy dynamics described in the previous section.

In the remainder of this paper, we adopt the following convention: we write with a prime (e.g.

), the diagnosed eddy field from the CTRL run, and with a dagger (e.g. † ) the prognostic eddy dynamics that result from the explicit time integration of the sub-grid model (Eq. 12) with the mean flow ( , ) from the CTRL run as the input (Table 1). Namely, where we have replaced the primes with daggers to signify the reinterpretation from eddy to sub-grid. We have also replaced ( , ) by R as we are going to design a parameterization for ( , ). Our aim is to build a sub-grid model for which With the lack of a good predictor for the eddy rectification forcing, we can start by examining the sub-grid model (Eq. 17) without it on the right-hand side (viz. R = 0). We recall that Eq. ( 17) with R = 0 has mostly been used to simulate local turbulence in doubly-periodic patches of the ocean with uniform shear (e.g. [START_REF] Venaille | Baroclinic turbulence in the ocean: Analysis with primitive equation and quasigeostrophic simulations[END_REF][START_REF] Grooms | Mesoscale eddy energy locality in an idealized ocean model[END_REF]), whereas we now apply and solve this equation prognostically in the entire domain with a large-scale flow that varies in space.

† + ( † , † ) + ( , † ) + ( † , ) + † = 4 ∇ 4 † + ∇ 2 † + R , (17
For white noise initial conditions, we can decompose the run in several stages: we first observe a linear growth of the most unstable modes mainly in the jet and near the northern and southern boundary. The duration of this phase is on the same order of magnitude as the inverse linear growth rate (see Appendix A, Fig. A1a). We then enter another transient phase during which a large-scale pattern emerges in the PV field, and after this transient phase, we reach a statistical steady state (i.e. ≈ 0). To illustrate this last regime, we plot in Fig. 4 the mean potential and kinetic energy as well as snapshot of these two fields. There are several important things to notice: first we note that † (Fig. 4d) is very different from (Fig. 2d):

† is maximum along the western boundary and does not really reflect the eddies that were present in the jet in the CTRL run. In fact when we look at a snapshot of potential energy (Fig. 4b), we see that this potential energy field is the sum of a large-scale and small-scale flow.

Everywhere in the domain, the mean kinetic energy in this sub-grid run (Fig. 4c) is weaker than the mean eddy kinetic energy diagnosed from the CTRL run (Fig. 2c), viz. † < . The lower energy levels in eddy kinetic and potential energy is also apparent in the isotropic wavenumber spectra (Fig. 5; compare the black solid and dotted lines). We compute the eddy kinetic and potential energy spectra ( | û| 2 2 and | ˆ | 2 2 2 respectively where ( •) is the Fourier transformed amplitude) of the first layer over the whole domain using the xrft Python package (Uchida et al. 2021b) and taper the fields with the Hann window prior to taking the Fourier transform as is commonly done when computing the spectra [START_REF] Uchida | Seasonality of eddy kinetic energy in an eddy permitting global climate model[END_REF][START_REF] Khatri | Surface ocean enstrophy, kinetic energy fluxes, and spectra from satellite altimetry[END_REF][START_REF] Ajayi | Spatial and temporal variability of the north atlantic eddy field from two kilometric-resolution ocean models[END_REF]). The periodogram is computed every 23 days over the last 580 days of output and then averaged.

In the eddy run, we still see a local kinetic energy ( † ) maximum in the middle of the domain where the mean jet is and we also observe deformation radius size eddies in the rest of the gyre (Fig. 4a). Such difference between † and † where we see larger scale patterns in the former indicates that in this eddy run, energy is stored in the large-scale buoyancy field rather than in small-scale eddies. We interpret these energy maps in the light of the inverse cascade in quasi geostrophy that fluxes energy toward larger scales [START_REF] Charney | Impact of horizontal resolution (1/12 • to 1/50 • ) on gulf stream separation, penetration, and variability[END_REF][START_REF] Vallis | Atmospheric and Oceanic Fluid Dynamics[END_REF]). Because of this inverse cascade, we see the appearance of a large-scale pattern superimposed on top of the prescribed large-scale circulation (i.e. and in Eq. 12). The sum of these two large-scale solutions as we see in Fig. 4d corresponds to a less baroclinically unstable state and hence weaker eddies (see Fig. 4a).

We also plot in Fig. 6a the eddy stream function for the same snapshot as the one plotted in Fig. 2, and in Fig. 6b the sub-grid stream function of the sub-grid model for the same snapshot as in Fig. 4.

This plot confirms the differences already highlighted of a weaker baroclinicity in the eddy run and also shows that large-scale Rossby waves present in the eddy field diagnosed from the CTRL run ( ; Fig. 6a) are not present in the eddy model ( † ; Fig. 6b). This is probably because Rossby F . 5: The isotropic wavenumber spectra taken over the whole domain for kinetic and potential energy in the first layer (a,b). The energies diagnosed from the CTRL run are shown in solid black, from the sub-grid model with no forcing in dotted red (R = 0), and from the sub-grid model with the varying spatial filter approach in solid red lines respectively. The sub-grid models at coarser resolutions (256 and 128 grid points; Appendix C) are shown in green dashed and dotted lines respectively.

the result of the rectification of the large-scale flow by small-scale eddies: the eddies tend to create a flow that opposes the large-scale forcing from the CTRL output ( ). As already noted with the energy diagnostics, the intensity of the eddy activity increases near the central latitude and near the western boundary. Near the central latitude, the eddies tend to form an eastward jet, which is also the opposite of what is observed in the CTRL run (a western boundary current that penetrates into the domain as a westward flowing jet). Although a similar mechanism of the eddies counteracting the mean flow is well-known in the Southern Ocean where the overturning circulation by eddies counter balance the mean Ekman steepening of isopycnals (e.g. Sinha and Abernathey 2016), we conclude that the solution produced by the sub-grid model ( † ) is not a fair reproduction of the eddy dynamics in the CTRL run ( ; Fig. 6). We show in section 3b, however, that we have some success in recovering the eddy dynamics from the dagger fields by parameterizing the eddy rectification forcing.

We now focus on the rectification term ( † , † ) (the mean of second term on the left-hand side of Eq. 17) that emerges in this simulation from the white-noise initial condition and plot this field in Fig. 7. The field is smoothed in a similar manner to as described in section 2d where we purposes. The smoothed ( † , † ) shares many common features with the diagnosed rectification term ( ( , ); Fig. 3): both fields are positive (negative) in the subpolar (subtropical) gyre. The magnitude of this term is intensified in the region of the separated jet with roughly the same alternance of positive and negative pattern. Lastly, the boundary dynamics is also of the same sign. The main difference is that the simulated field ( † , † ) is weaker in magnitude than the diagnosed field (Fig. 3). The agreement in spatial patterns between these two fields is pleasing given the discrepancies of the dynamics in the two simulations (cf. Figs. [START_REF] Sinha | Time scales of southern ocean eddy equilibration[END_REF]4,6).

This experiment suggests that eddy dynamics feedback onto the large-scale dynamics via the inverse cascade. In the eddy model, this feedback on the large-scale potential energy concurs to flatten isopycnal surfaces and effectively shuts off the generation of eddies via baroclinic instability.

We conclude that although the term ( , ) has no impact on the domain-integrated energetics of the eddy flow, it is actually very important to counteract the inverse cascade and prevent the formation of spurious large-scale mode in the eddy flow. Even though the stream function we get from the sub-grid model is different from the diagnosed eddy stream function from the CTRL run, we get at this point a viable candidate for the rectification of the large-scale flow by the eddies ( ( † , † )). This result itself is already a big improvement compared to the regional models of [START_REF] Venaille | Baroclinic turbulence in the ocean: Analysis with primitive equation and quasigeostrophic simulations[END_REF]. We recall that the main difference between the present study and Venaille 2011) is that they used regional model with periodic boundary conditions whereas we run the eddy model for the entire domain. With this strategy we do capture the non-local eddy/mean-flow interaction that is impossible to capture with regional models. In the remainder of this section, we propose a parameterization for R in Eq. ( 17) and show that we can improve the eddy statistics and

( † , † ).

b. Parameterizing the eddy rectification forcing (R)

As noted earlier, the field ( † , † ) converges very slowly, and so cannot be computed in practice as a parameterization of R to run the sub-grid model. In order to parameterize ( † , † ), we propose to use the idea developed by [START_REF] Pedlosky | The equations for geostrophic motion in the ocean[END_REF], [START_REF] Grooms | On the interactions between planetary geostrophy and mesoscale eddies[END_REF] and others to decompose the flow into a large-scale component and a small-scale component. In a similar way to the definition of the time mean, we introduce the spatial scale decomposition for a field as

= + * , (18)
where and * are respectively the large-scale and small-scale components of the field . Based on Pedlosky's scale decomposition, the large-scale flow evolves on a slow time scale and the smallscale flow evolves on a fast time scale [START_REF] Pedlosky | Geophysical Fluid Dynamics[END_REF]. We accomplish such scale decomposition Accepted for publication in Journal of Physical Oceanography. DOI 10.1175/JPO-D-21-0217.1.

Unauthenticated | Downloaded 04/07/22 08:27 AM UTC by enforcing † to remain a small-scale field † = 0 ,

which implies, if we set the initial condition of † | =0 = 0 then † = † * is satisfied for all time.

Because of the equivalence between the slow time scale and large-scale spatial scale, our hope is that enforcing Eq. ( 19) will be equivalent to enforcing † = 0. Note that in the run with R = 0, we clearly did not have † = 0. We use this argument to parameterize R as a damping of the large-scale flow

R = - † , (20)
where † in Eq. ( 17) is relaxed towards zero on a three-day time scale ( ; see Appendix B2 for details on the numerical implementation). With this parameterization of the rectification term, we can already anticipate that the spatial filtering strategy will not work well in the region of the separated jet where there is no clear scale separation between the eddy flow and the mean flow (cf. [START_REF] Jamet | On wind-driven energetics of subtropical gyres[END_REF]. However, as we shall see, this strategy works well in the rest of the domain.

We illustrate the effect of the spatial filter operator (Eq. 18) in Fig. 8 where we plot the same sub-grid stream function as the one used in Fig. 4 along with its large-scale and small-scale component. We do this scale separation by applying a low-pass filter with a discrete wavelet transform (numerical details of the implementation are provided in Appendix B). In Fig. 8, we use a cutoff length scale of = 500 km. In the large-scale pattern, we recognize a cyclonic and anticyclonic gyre, and a weak jet in the middle that we described earlier.

In order to use this spatial filter as a parameterization of R, we performed many tests with either uniform or spatially varying . Henceforth, we present our best results obtained with non-uniform , and now briefly describe the reasons for opting for a non-uniform . We see in Fig. 2a that the patch of high eddy kinetic energy has horizontal dimensions on the order of 1000 km. In the region of the separated jet, there is thus no clear scale separation between the eddy flow and the mean flow. From hereon, when we refer to the sub-grid model (Eq. 17), R is that of described in Eq. ( 20) (i.e. the linear damping of low-pass filtered sub-grid PV). We now run the sub-grid model with the same mean flow ( ) as in section 3a: namely the mean variables from the CTRL run.

We plot the energy diagnostics in Fig. 10. Comparing Fig. 10c,d with Fig. 4c,d, we see that using this R (Eq. 20) in Eq. ( 17) succeeds in increasing the eddy amplitude overall and in particular around the separated jet (compared to the solution with R = 0). The energy levels come closer to the eddy field diagnosed from the CTRL run (Figs. 2c and2d), which is also apparent in the isotropic wavenumber spectra (Fig. 5). We see clear increases in energy from the run with R = 0 and that the varying spatial filter approach captures energy levels close to the diagnosed eddy kinetic and potential energy except for the smallest wavenumbers (largest spatial scales; compare the black solid and red dashed lines in Fig. 5). This is expected as we extract the large-scale component with F . 9: The cut-off length scale ( ) in meters based on the instability length scale.

If our parameterization were perfect, time averaging Eq. ( 17) would return the balance

( † , † ) = R (21)
because the linear terms should vanish. Although the balance in Eq. ( 21) requires there to be a clear scale separation between the eddy and mean flow, we expect the balance to approximately hold, viz. † ∼ † ∼ 0 for a converged simulation. We plot in Fig. 11,( † ,† ) smoothed by 16 neighboring grid points and R. (The difference between Fig. 11a and 7 is in Eq. ( 17) prognostically solved with and without the eddy rectification parameterization on the right-hand side.) We first see that ( † , † ) captures the same patterns as the diagnosed field from the CTRL run ( ( , ); Fig. 3b). We also see improvements compared to the run without the rectification forcing (R = 0).

Along with this visual comparison, we plot in Fig. 12 the joint histogram of ( , ) and ( † , † ).

We see that this joint histogram aligns more around the one-to-one line with the varying spatial filter approach compared to when R = 0. If we now compare ( † , † ) and R, we see that the latter captures the large-scale pattern in the return flow of the gyre but misses the small-scale variability in the separated jet and right at the western boundary. We could have anticipated the lack of small-scale variability in R because of the nature of our filter which only retains the large-scale component of ( † , † ). In the separated jet, the agreement between ( , Lastly, we note that Eq. ( 21) is complimentary to a recent work by Porta [START_REF] Mana | Toward a stochastic parameterization of ocean mesoscale eddies[END_REF] and [START_REF] Grooms | A note on 'toward a stochastic parameterization of ocean mesoscale eddies[END_REF] where they find a local relation ( * , * ) ∇ 2 . We emphasize that by explicitly solving for Eq. ( 17) and parameterizing the eddy rectification forcing with Eq. ( 20), the parameterization incorporates non-local effects as it partially balances the advective term on the left-hand side. Notably, in a recent work, [START_REF] Berloff | On dynamically unresolved oceanic mesoscale motions[END_REF] achieved such non-local closure by diagnosing the eddy rectification forcing term as the mismatch between the left-hand and right-hand side of a coarse-grained PV equation, viz.

)
R + ( , ) + + 4 ∇ 4 + ∇ 2 -, (22)
and then plugging it along with , into the sub-grid equation (Eq. 17). While our approach is similar, the difference is in how the eddy rectification forcing is defined: we define it by applying a low-pass spatial filter to the sub-grid stream function (Appendix B; whereas they use the full PV equation).

Modification of the mean flow due to the eddy rectification term

The procedure described in the previous section demonstrates that the sub-grid model can fairly reproduce the "true" eddy dynamics given a prescribed background flow. There is one caveat, however, which is precisely the specification of this background flow. Indeed, from an eddy parameterization perspective, taking as the mean of high resolution model (as we did so far) is very different than taking from a coarse resolution model which has never seen properly resolved eddies. This is because the unstable modes are very different depending on whether the eddies are 25 Accepted for publication in Journal of Physical Oceanography. DOI 10.1175/JPO-D-21-0217.1.

Unauthenticated | Downloaded 04/07/22 08:27 AM UTC resolved or not and so we expect the eddy dynamics to be function of the background flow. In this section, we first explore the sensitivity of the eddy model with respect to the background flow. We then study how our sub-grid model modifies the mean flow by feeding back onto it via the eddy rectification forcing.

a. Non-eddying full model and mesoscale-resolving sub-grid model

In order to see how the eddy model performs in the more realistic situation where comes from a coarse resolution model, we now run a coarse full QG model (Eq. 5) with the same parameters as CTRL, except we lower the resolution to Δ ≈ 78.13 km, increase the bi-harmonic viscosity to 4 = 6.25 × 10 11 m 4 s -1 , and also use a harmonic viscosity with 2 = 1000 m 2 s -1 . Hereon we call this configuration the REF run. In this coarse resolution model, the flow converges to a stationary state with almost no variability. This mean flow has less potential energy than the CTRL run and the mid-latitude eastward jet is very weak (see Fig. 13). Note that we spun up the coarse full model with white noise initial conditions but without any rectification term.

The sub-grid model itself is still Eq. ( 17 and13a). As a consequence, the eddy rectification of the separated jet in the domain interior is negligible (Fig. 13b).

b. Impact of the rectification on the large-scale flow

In order to see how we can use this eddy parameterization for coarse resolution models, we turn our attention to Eq. ( 10). The only difference between this equation and the full model is the presence of the rectification term ( , ) and the purpose of this study was to propose a closure for this term. We can now use either ( † , † ) or R as an approximation for ( , ) and plug it into Eq. ( 10 10) with R. We first note in Fig. 13b that the magnitude of this term is comparable to the wind stress forcing in the western part of the basin (not shown). Also, compared to the wind forcing, this term has a vertical structure (not shown; whereas the wind forcing is only present in the surface layer). Hence, we expect the rectification to have a significant impact on the mean flow.

When we integrate in time the coarse resolution model with the rectification term, the circulation changes in a few places. We plot in Fig. 14a-b the change in the stream function when we force the coarse model with ( † , † ) and R respectively. Both of these runs undergo very similar changes so it does not really matter which of these term we choose to force the coarse resolution model.

Both runs exhibit a weakening of the western boundary current (patch of color of the opposite sign as the mean circulation). However, the rectification strengthen the separated jet (patch of color of the same sign as the mean circulation).

Since the resolution of the full model is non-eddying, a common eddy parameterization to implement would be the GM parameterization (cf. Eq. 11). We implement it in the QG model and 27 Accepted for publication in Journal of Physical Oceanography. DOI 10.1175/JPO-D-21-0217.1.

Unauthenticated | Downloaded 04/07/22 08:27 AM UTC we use a diffusivity coefficient ( GM = 1000 m 2 s -1 applied only to buoyancy, equivalently the layer thickness in quasi geostrophy; cf. Uchida et al. 2021a). As GM is intended to mimic the baroclinic process of reducing PE, it would tend to further weaken the separated jet, which is what we see over the entire domain (blue in the subtropical and red in the subpolar gyre; Fig. 14c). The two runs with the eddy rectification forcing, on the other hand, tends to sharpen and strengthen the jet upon separation near the western boundary as we see between the meridional extent of 150-350 km (Fig. 14b,c). In other words, our closure captures the energy backscattering from the "sub-grid" eddies onto the coarse full flow as they would if the eddy model were run until it reaches statistical convergence (see the similarity between Fig. 14b,c). The benefit of using R instead of ( † , † ) is that it converges much faster than directly diagnosing ( † , † ), reducing the computational cost by a factor of (10 2 ). We have shown that for a non-eddying resolution, our closure provides a potential path forward to go beyond GM.

Conclusions and discussion

In this study, we have examined the eddy rectification term, which encapsulates the net eddy feedback onto the mean flow, from a quasi-geostrophic (QG) double gyre simulation. In doing so, we decompose the QG potential vorticity (PV) into its mean flow, defined by a time mean, and eddies as the fluctuations about the mean. This paper is an attempt to estimate the rectification term ( , ) based on the knowledge of the mean flow only. For that purpose, we solve an eddy equation that describes the dynamics of the perturbation around that mean flow. Since we solve for the perturbation equation, we now need a closure for non-linear interaction between the perturbation variables as is always the case in closure problems. We have shown that we can use the eddy model (Eq. 17) to diagnose the eddy rectification term without any closure. With R = 0, the eddy model gives a rough estimate for the rectification term diagnosed from the mesoscale-resolving full model, viz. ( † , † ) ∼ ( , ) (Figs. 3b and7). The improvement compared to previous studies for which local closures were developed in a doubly periodic regional model [START_REF] Venaille | Baroclinic turbulence in the ocean: Analysis with primitive equation and quasigeostrophic simulations[END_REF]) is that we solve the eddy model at the basin scale thus allowing non-local eddy feedback. However, the fact that a large-scale component of the sub-grid stream function itself ( † ) emerges opposing the background flow without the eddy rectification forcing, which is not apparent in the eddy stream function diagnosed from the full model ( ), perhaps warrants some attention (Figs. Once the eddy rectification forcing is estimated from the (sub-grid) eddy model (R; Eqs. 10 and 17), we can then use this term in the mean flow model (Eq. ( 10)) as a forcing term on the right-hand side. For a coupled system between the mean flow and sub-grid model, this leads to a process where we march forward in time by: i) re-interpreting the mean flow model as the full model at non-mesoscale-resolving resolutions, ii) feeding the resolved flow to the sub-grid model as the background flow with the parameterization for the eddy rectification forcing (R), and iii) from which we force the full model with the eddy rectification forcing estimated from the eddy model (R). This is similar to other energy backscatter parameterization studies where they solve the (sub-grid) eddy energy equation and take that as a forcing for the resolved momentum equation (e.g. [START_REF] Jansen | Toward an energetically consistent, resolution aware parameterization of ocean mesoscale eddies[END_REF][START_REF] Josey | Wind stress forcing of the ocean in the SOC climatology: Comparisons with the NCEP NCAR, ECMWF, UWM/COADS, and Hellerman and Rosenstein datasets[END_REF][START_REF] Perezhogin | Deterministic and stochastic parameterizations of kinetic energy backscatter in the nemo ocean model in double-gyre configuration[END_REF]. Here, we have formulated a deterministic closure based on PV instead of energy; PV is a more fundamental variable in quasi geostrophy as it is materially conserved while energy is only conserved in the volume integrated sense. Our approach of parameterizing the eddy rectification term via a spatially-filtered eddy stream function is complementary to a recent work by [START_REF] Mana | Toward a stochastic parameterization of ocean mesoscale eddies[END_REF] and [START_REF] Grooms | A note on 'toward a stochastic parameterization of ocean mesoscale eddies[END_REF] where they find a closure for the rectification term in relation to the low-pass filtered PV. One major difference here is that while their closure was local, we have accounted for non-local effects by approximating the eddy rectification forcing prognostically from the eddy model (cf.

Berloff et al. 2021).

As a first step towards a PV-based coupled closure, we have emphasized the importance of solving the sub-grid model explicitly and provided a proof of concept by solving the 'partially' coupled system within the QG framework. We emphasize 'partially' as the eddy rectification forcing we gave the full model at non-eddying resolution was the time mean of the rectification predicted from the sub-grid model (R). This has to do with the fact that we decompose the eddy-mean flow with a temporal averaging. While the temporal averaging was chosen originally to examine the eddy model under a prescribed double-gyre background flow and to allow for commutability between the averaging operator and spatial derivatives, this makes the coupling process and interpretation convoluted in our case. In other words, if the averaging operator were orthogonal to the time dimension, we would have total = coarse + † at each time step where coarse here is the full PV resolved at coarse resolution. In such case, the total eddy kinetic energy would become total = coarse + † , where we would be able to directly compare it with the 30 Accepted for publication in Journal of Physical Oceanography. DOI 10.1175/JPO-D-21-0217.1.

Unauthenticated | Downloaded 04/07/22 08:27 AM UTC eddy kinetic energy from the CTRL run. Nevertheless, we have shown that our time-mean eddy rectification forcing sharpens the jet as the eddies would if they were resolved when the full model is non-eddying (section 4).

We also tested a case where the full model was mesoscale-permitting (Δ = 19.5 km; 4 = 6.25 × 10 11 m 4 s -1 ). The idea was to examine how an eddy model would perform if the full model also partially resolved the eddies. We followed the same procedure as described in section 4: i) run the full model without the rectification (R = 0), ii) diagnose the time-mean rectification (R) from the sub-grid model taking the mean flow from the full model as its background flow, and iii) plug the rectification into the full model as forcing. However, as the full model was already baroclinically unstable and partially resolved eddies, the process led to the full model having weaker eddies in step (iii); the eddies which fed off of the mean flow of the full model in step (ii) resulted in giving a rectification forcing that actually reduced the instability of the full flow in step (iii). In hindsight, this may have been expected as the eddies, if resolved, tend to extract PE from the background flow. For the case where the full model was non-eddying, the resolved flow was never unstable so the reduction in PE upon iteration did not happen.

While we have attempted to design a deterministic super parameterization where one explicitly solves the sub-grid processes, it is possible that we are facing the limit of deterministic closures for the mesoscale-permitting regime. Stochastic and/or machine learning approaches may need to be considered [START_REF] Bauer | Stochastic representation of mesoscale eddy effects in coarse-resolution barotropic models[END_REF][START_REF] Guillaumin | Stochastic-deep learning parameterization of ocean momentum forcing[END_REF][START_REF] Frezat | Physical invariance in neural networks for subgrid-scale scalar flux modeling[END_REF]. Nonetheless, we have shown that our closure improves the eddy model in representing the eddies in comparison to them diagnosed from a mesoscale-resolving full model. Lastly, one may ask how our results can be extended to primitive equation models. In primitive equations, the eddy Ertel PV flux encapsulates the eddy feedback onto the mean flow [START_REF] Young | An exact thickness-weighted average formulation of the Boussinesq equations[END_REF]. In other words, a closure based on Ertel PV may allow one to capture the eddy variability in a primitive eddy model.

As an alternative to our spatial filtering approach, we hypothesize that it is possible to obtain the rectification term through iteratively solving for Eq. ( 12) as the Fixed-Point Theorem would predict.

As we discussed in section 3a, the sub-grid model without any forcing term (R = 0) produces a good first guess of the rectification term, namely the mean of ( † , † ) on the left-hand side of Eq. ( 12) (Fig. 7). The idea is then to re-run the sub-grid model with this first guess as the forcing term (R = ( † , † )) and repeat this iterative procedure until convergence is reached. We already 31 Accepted for publication in Journal of Physical Oceanography. DOI 10.1175/JPO-D-21-0217.1.

Unauthenticated | Downloaded 04/07/22 08:27 AM UTC know that this convergence is extremely slow (order of million of eddy time scale; section 2d) so this process cannot be practically done with the raw estimate of the rectification term but may be possible for its spatially smoothed version. The proof for mathematical convergence of this iterative process is beyond the scope of this study. energy levels similar to the eddy energies diagnosed from the CTRL run in the separated jet region (Fig. 2). We also see this from the wavenumber spectra where in the spatial range of ∼ 300 km, the level of EKE is similar between † and (Fig. 5). Moving to the coarsest resolution, we see that the jet penetration into the gyre deteriorates due to insufficient resolution and high viscosity prohibiting the instabilities to grow (Fig. C1c,d). The lack of energy is apparent in the wavenumber spectra where they fall off too quickly with wavenumber (Fig. 5).

With the numerical viscosity as a tuning parameter, we end this appendix by showing the dependency of the system on it. Figure C2 shows the ratio between domain integrated EKE

5

  Accepted for publication in Journal of Physical Oceanography. DOI 10.1175/JPO-D-21-0217.1. Unauthenticated | Downloaded 04/07/22 08:27 AM UTC with

F . 1 :

 1 A schematic of the four-layer configuration with a rigid lid and flat bottom (a). The background stratification is prescribed at the layer interfaces. (b) The meridional profile of the wind stress and a snapshot of the surface relative vorticity normalized by 0 . Note that the wind pattern takes only positive values: we could have added a term -0 /2 in the definition of the wind in Eq. (

F . 2 :

 2 is extended by a factor 20 compared to the other plots because there is approximately 20 times more potential energy than kinetic energy in the large-scale flow. This result corresponds to the traditional view of the ocean circulation ( ∼ > ), expressed here in the QG framework. In Fig. 2f, we see the bowl shape of the anticylconic gyre in the southern part of the domain (and respectively the dome shape of the cyclonic gyre in the northern part of the domain). Potential energy is maximum in the middle of the gyre where the buoyancy anomaly is maximum. The mean jet is much less energetic as shown in the kinetic energy panel (Fig. 2e). 10 Accepted for publication in Journal of Physical Oceanography. DOI 10.1175/JPO-D-21-0217.1. Snapshots and time-mean of potential energy and kinetic energy diagnosed from the CTRL run. A snapshot of the eddy kinetic and potential are shown in panels (a) and (b). Their time mean in panels (c) and (d). The mean kinetic and potential energy are shown in panels (e) and (f). Units: m 2 s -2 11 Accepted for publication in Journal of Physical Oceanography. DOI 10.1175/JPO-D-21-0217.1. Unauthenticated | Downloaded 04/07/22 08:27 AM UTC

12

  Accepted for publication in Journal of Physical Oceanography. DOI 10.1175/JPO-D-21-0217.1. Unauthenticated | Downloaded 04/07/22 08:27 AM UTC (a) Raw J(ψ , q ) The raw ( , ) and ( , ) smoothed by averaging 16 neighboring grid points and linearly interpolated back on the fine grid (a,b).

  ) 13 Accepted for publication in Journal of Physical Oceanography. DOI 10.1175/JPO-D-21-0217.1. Unauthenticated | Downloaded 04/07/22 08:27 AM UTC

1 :

 1 † , † and ( † , † ) mimic , and ( , ) diagnosed from the CTRL run. Eddy terms about the time mean diagnosed from the CTRL run (•) † Prognostic eddy terms from the sub-grid model , Time mean of the eddy kinetic and potential energy diagnosed from the CTRL run † , † Time mean of the kinetic and potential energy diagnosed from the sub-grid model ( , ), ( † , † ) Eddy rectification diagnosed from the CTRL run and sub-grid model respectively by taking the time mean of the simulation outputs R Eddy rectification forcing (i.e. the target of parameterization) T Definition of the notations. a. No eddy rectification forcing (R = 0)

F . 4 :

 4 waves in the full model are triggered by intense eddies in the meandering jet. Since this model only produces mild eddies, there are no Rossby wave that will emerge in the eddy model. Another possibility is that Rossby waves are excited by the winds ( in Eq. 9), which project themselves 15 Accepted for publication in Journal of Physical Oceanography. DOI 10.1175/JPO-D-21-0217.1. Snapshots and time-mean of kinetic and potential energy diagnosed from the eddy model with no forcing (R = 0). The snapshots are shown in panels (a) and (b), and their time means in (c) and (d) respectively. Units: m 2 s -2 . onto the temporally varying fields of , whereas the sub-grid model ( † ) has no input to excite such waves. The interesting point is that without the eddy rectification forcing, the large-scale pattern in † that emerges corresponds to a the cyclonic gyre (in blue) is in the southern part of the domain and the anticyclonic gyre (red) is in the northern part of the domain (Fig. 6b), which is precisely the opposite from the stream function in the CTRL run. We interpret this large-scale pattern in † as 16 Accepted for publication in Journal of Physical Oceanography. DOI 10.1175/JPO-D-21-0217.1. Unauthenticated | Downloaded 04/07/22 08:27 AM UTC

  average 16 neighboring grid points and linearly interpolate back on the fine grid for visualization 17 Accepted for publication in Journal of Physical Oceanography. DOI 10.1175/JPO-D-21-0217.1. The eddy stream function diagnosed from the CTRL run and sub-grid stream function † simulated from the sub-grid model with no forcing (R = 0) (a,b).

18

  Accepted for publication in Journal of Physical Oceanography. DOI 10.1175/JPO-D-21-0217.1. Unauthenticated | Downloaded 04/07/22 08:27 AM UTC J(ψ † , q † ) (with R = 0) ( † , † ) diagnosed from the sub-grid model without forcing (R = 0), smoothed by averaging 16 neighboring grid points and linearly interpolated back on the fine grid. et al. (

  To a certain extent, this corroborates what we observe in the instability analysis (Appendix A). In Fig. A1, we see that in the region of the separated jet, the most unstable mode has a characteristic length scale = 300 km compared to the most unstable length scale in the return flow which is = 230 km. We use this information to build a filter with non-uniform 20 Accepted for publication in Journal of Physical Oceanography. DOI 10.1175/JPO-D-21-0217.1. Low pass and high pass filtered sub-grid stream function diagnosed from the sub-grid model with R = 0 (a,b). The eddy stream function spatially decomposed is the one in Fig. 6b. length scale in the form of = , and we set = 4.5 to get ∼ (1000 km) in the area of the return flow. We plot in Fig. 9 the final map of which corresponds to a smoothed version of the most unstable length scale (see Appendix A). As desired, has values on the order of 1000 km with a maximum of 1350 km in the region of the separated jet and a minimum of 850 km near the north-east and south-east corners.

  the spatial filter. 21 Accepted for publication in Journal of Physical Oceanography. DOI 10.1175/JPO-D-21-0217.1. Unauthenticated | Downloaded 04/07/22 08:27 AM UTC λ c (m)

F

  and R is poor and we 22 Accepted for publication in Journal of Physical Oceanography. DOI 10.1175/JPO-D-21-0217.1. . 10: Potential energy and kinetic energy diagnosed from the eddy model where R is implemented with the variable length scale filter. The snapshots are shown in panels (a) and (b), and their time means in (c) and (d) respectively. Units: m 2 s -2 . face here the limits of approximating the time average by a low-pass spatial filter (Eq. 20). Reducing the length scale of the filter is problematic because it degrades the quality of the eddy solution (not shown). Nevertheless, even with this bias, the rectification term ( † , † ) compares well with the diagnosed rectification ( ( , ); Figs. 3 and 11). We show the resolution dependence of the sub-grid model given the same background flow from the CTRL run in Appendix C. 23 Accepted for publication in Journal of Physical Oceanography. DOI 10.1175/JPO-D-21-0217.1. Unauthenticated | Downloaded 04/07/22 08:27 AM UTC

F

  (a) J(ψ † , q † ) (varying spatial filter) . 11: ( † , † ) diagnosed from the eddy model with the varying spatial filter approach, smoothed by averaging 16 neighboring grid points and linearly interpolated back on the fine grid, and R = † / (a,b).

  : Joint histogram of the spatially smoothed ( , ) and ( † , † ) for the first and second layer plotted against a logarithmic scaling (the masked out regions have zero values). The left column shows the run with no forcing (R = 0) (a,c), and right the run with the varying spatial filter approach (b,d). The one-to-one line is shown in grey dashed lines. The histograms were computed using the xhistogram Python package[START_REF] Abernathey | xhistogram: Fast, flexible, label-aware histograms for numpy and xarray[END_REF]).

  ) which now takes the time mean of the coarse model as barred variables, and for R we use a spatial filter with uniform cutoff length scale = 1000 km (simply because the unstable modes of the mean flow of REF exhibit an almost uniform pattern for both the instability time scale and the instability length scale). A snapshot of the eddies and diagnosed eddy rectification from the eddy model are shown in Fig. 13. The eddy activity resemble the CTRL run near the western boundary but lacks the signature in the separated jet region (Figs. 2a

  ) to see how it would in turn modify the flow of the coarse resolution model; note that R is barred: R = - † / . Under stationary forcing conditions as we have set up here ( = ), Accepted for publication in Journal of Physical Oceanography. DOI 10.1175/JPO-D-21-0217.1. Unauthenticated | Downloaded 04/07/22 08:27 AM UTC

F

  . 13: Snapshot of the EKE of the eddy model driven by the low resolution background flow, namely the outputs from REF (a). The contours show the time-mean reference stream function from the low resolution REF run. (b) ( † , † ) diagnosed from the high resolution eddy model driven by the low resolution background flow, and then smoothed by 16 neighboring grid points and linearly interpolated back on the fine grid. a converged flow would give ∼ 0. Hence, we gave the eddy rectification forcing as its time mean to a priori remove time dependency. Namely, we replaced ( , ) on the right-hand side of Eq. (

F

  6 and 28 Accepted for publication in Journal of Physical Oceanography. DOI 10.1175/JPO-D-21-0217.1. Unauthenticated | Downloaded 04/07/22 08:27 AM UTC (a) Forced with J(ψ † , q † ) . 14: Color: The difference in stream function between the coarse reference run with R = 0 and coarse runs with eddy closures ([m 2 s -1 ]). Contours: Stream function of the low resolution REF run. The run with ( † , † ) smoothed by 16 neighboring grid points (a), with R (b), and with GM only (c). 8a). We have shown that approximating the eddy rectification forcing with the spatially-filtered eddy PV (R = - † / ) improves the eddy kinetic and potential energy and ( † , † ) (Figs. 5, 10-12). In other words, we have provided a closure to circumvent the necessity to diagnose the mean properties of eddy-eddy interaction from an eddy resolving simulation (section 3). 29 Accepted for publication in Journal of Physical Oceanography. DOI 10.1175/JPO-D-21-0217.1. Unauthenticated | Downloaded 04/07/22 08:27 AM UTC

32

  Accepted for publication in Journal of Physical Oceanography. DOI 10.1175/JPO-D-21-0217.1. Unauthenticated | Downloaded 04/07/22 08:27 AM UTCthe inverse growth rate of the most unstable mode, and , the zonal and meridional wavelength of that most unstable mode, and of that mode. We plot T and in Fig.A1. One first important information from these plots is that the large-scale solution is unstable almost everywhere in the domain (except in the small white area at = 2500 km near the eastern boundary). This was not obvious a priori because we computed the most unstable mode with the same viscosity as the CTRL run and viscosity is known to damp instabilities. We divide the time scale pattern into three distinct dynamical regimes: the western boundary and the intergyre jet which have the fastest growing mode (order 20 days), the return flow near the northern and southern boundary for which the instability time scale is order 60 days, and the rest of the domain for which the instability time scale is greater than 115 days (the colorbar saturates beyond this value). We do not consider the instability with long time scale because such long time scale is much bigger than the eddy time scale and become irrelevant for the eddy dynamics (local instability analysis is probably not relevant in areas with such long time scales). The instability length scale is noisier but overall in the area where T < 115 days, the length scale of the instability is 10 times the deformation radius (consistent with the canonical 2-layer baroclinic instability; Cushman-Roisin and Beckers 2011).

  : Time scale and length scale of the most unstable mode (computed at every fourth grid point) (a,b).

  : The time-mean of kinetic and potential energy diagnosed from the eddy model at coarser resolutions with the varying spatial filter. The energies from the run with 256 grids are shown in panels (a) and (b), and 128 grids in (c) and (d) respectively. Units: m 2 s -2 .

Accepted for publication in Journal of Physical Oceanography. DOI 10.1175/JPO-D-21-0217.1. Unauthenticated | Downloaded 04/07/22 08:27 AM UTC

Acknowledgments. We would like to acknowledge the editor Paola Cessi along with Elizabeth Yankovsky and two other anonymous reviewers for their comments, which led to significant improvements in the manuscript. We wish to thank Antoine Venaille, Laure Saint Raymond and William K. Dewar for their insightful comments and suggestions. This work has been supported by the French national program LEFE/INSU. Uchida acknowledges support from the French 'Make Our Planet Great Again' (MOPGA) initiative managed by the Agence Nationale de la Recherche under the Programme d'Investissement d'Avenir, with the reference ANR-18-MPGA-0002. Accepted for publication in Journal of Physical Oceanography. DOI 10.1175/JPO-D-21-0217.1. Unauthenticated | Downloaded 04/07/22 08:27 AM UTC Accepted for publication in Journal of Physical Oceanography. DOI 10.1175/JPO-D-21-0217.1. Unauthenticated | Downloaded 04/07/22 08:27 AM UTC Accepted for publication in Journal of Physical Oceanography. DOI 10.1175/JPO-D-21-0217.1. Unauthenticated | Downloaded 04/07/22 08:27 AM UTC

Data availability statement. The open-source software for the QG model can be found at github. com/bderembl/msom. It was developed as a module of Basilisk (available at basilisk.fr). Simulation outputs are available upon request. APPENDIX A

Linear stability analysis

In this appendix, we perform a linear stability analysis of the mean state of the CTRL run described in section 2d. Methods to perform such analysis have been reported elsewhere (e.g. [START_REF] Vallis | Atmospheric and Oceanic Fluid Dynamics[END_REF][START_REF] Smagorinsky | The geography of linear baroclinic instability in earth's oceans[END_REF][START_REF] Tulloch | Scales, growth rates, and spectral fluxes of baroclinic instability in the ocean[END_REF][START_REF] Uchida | Seasonality of eddy kinetic energy in an eddy permitting global climate model[END_REF]) and we only recall the main steps here. From the eddy equation (Eq. 12), we drop the non-linear advective term as well as the rectification term and replace by one Fourier component

where stands for complex conjugate. For each Fourier component, we get an equation with four unknown: ( ), , , and , respectively the vertical structure of the Fourier mode, the zonal, meridional, and temporal wave number. We span the ( , ) space in order to find ( ) and , which are the eigenvector and the eigenvalue of the equation. If the imaginary part of is negative, the corresponding mode is exponentially decaying and the solution is stable but if the imaginary part of is positive, the solution is unstable. In the ( , ) space, the most unstable mode corresponds to the solution for which ( ) is maximum. We call

When we compare these plot with Fig. 2c, there does not seem to be an obvious link between the local instability parameter and the observed eddy kinetic energy. The path of the jet has a wider signature in the map than in the instability analysis. The demarcation between the return flow and the rest of the gyre that we observe in Fig. A1a also does not show up in the kinetic energy map. This confirms the conclusion of [START_REF] Grooms | Mesoscale eddy energy locality in an idealized ocean model[END_REF] who showed that the eddies observed at one given location are mostly not locally generated but emanate from areas afar (see also [START_REF] Venaille | Baroclinic turbulence in the ocean: Analysis with primitive equation and quasigeostrophic simulations[END_REF]).

We use these two fields to build the length scale cutoff of the spatial filter. We start by simply setting = . However, we argue against using the raw value of as shown in Fig. A1b as this field is noisy and also because some instabilities are not relevant to the dynamics. The instabilities irrelevant to mesoscale dynamics occur in places where the instability time scale is greater than the advection time scale (which is on the order of 20 days in most of the gyre, not shown). To get rid of the non-relevant unstable modes, we adjust the value of to 225 km everywhere where T > 115 days. We then smooth that field with a Gaussian filter with a standard deviation of 4.5 grid points to get rid of the grid scale variations. Lastly, for each point of the domain, we create a halo of size over which we propagate the value of . We take = 4.5. This is done to let enough space for all instabilities to develop around the formation site. Several halos overlap at one point and so for each point we retain the maximum value of all halos that are present at that point.

We smooth the final map to damp the halo pattern that may have persisted. We plot the final map of in Fig. 9.

APPENDIX B

Numerical implementation

B1. Spatial filter

The discrete wavelet transform bears some resemblance with the multigrid solver. We define a set of grids from the finest model resolution 2 × 2 to the coarsest resolution 2 0 × 2 0 (one grid point). In our high resolution model (512 × 512), there are + 1 = 10 sets of grids. The two key operations in the filtering procedure are:

• The restriction R for which we coarsen a field by averaging 4 neighboring points;

• The prolongation P for which we refine a field by linear interpolation of neighboring points.

Let's suppose a field is defined on a grid of level (2 × 2 ). Then we have

We define the wavelet coefficients at level as ˘ = -P ( -1 ).

(B2)

Hence from the wavelet coefficients, one can reconstruct the field at the finest grid with an iterative procedure. The wavelet coefficients at level hold the information about the structure of the field at length scale of the grid size Δ . To high pass filter a field with a cutoff length scale = Δ , we simply need to set to zero the wavelet coefficients ˘ for < . In the case where varies smoothly in space, we can zero the wavelet coefficients locally only.

B2. Computation of R

We propose to approximate R as a damping term on the large-scale part of † as shown in Eq. ( 20). However, the filtering operation can be numerically expensive. Also, because the largescale component of † grows on a slow time scale, we chose to periodically (every three days)

remove the large-scale component of † in Eq. ( 17). We chose this three-day period because it is comparable to the eddy time scale and was short enough compared to the time needed for large-scale mode to build up observed in Fig. 8a, which is on the order of years. Lastly, we found that removing the large-scale component of † is less efficient than removing the large-scale component of † and then applying the linear operator L to † . With the latter technique, we take the derivative of the filtered field which does not create a spurious large-scale component. When the order of operation is the other way around (first filter † and then invert the elliptic equation (Eq. 1) to compute † ), we observed a spurious large-scale component in † . Hence, every three days, we add the term

to the right-hand side of Eq. ( 17) for only one time step (Δ ) and then set R = 0 the rest of the time.

This is equivalent to keeping R = -L ( † )/ constant for the three-day duration until we update it for the next three days. To see the equivalence, the number of time steps within every three days is = /Δ . Therefore, the cumulative effect of R over the three-day period is

where the left-hand side is what we have in Eq. ( 20).

This time scale separation is similar to ocean models where the barotropic and baroclinic modes are solved with different time stepping (cf. [START_REF] Marshall | A finite-volume, incompressible Navier Stokes model for studies of the ocean on parallel computers[END_REF]. The relaxation by our parameterization damps the large-scale component of † , i.e. † ∼ 0.

APPENDIX C

The sub-grid model at coarser resolution with a prescribed background flow

Given that the prognostic sub-grid model (Eq. 17) solved at mesoscale-resolving resolution is the best our method can achieve (section 3b), we examine the sensitivity of how our closure scales at coarser resolutions. We ran two additional cases of the sub-grid model with the resolution of ∼ 19.5 km and ∼ 39 km (256 and 128 grid points respectively) keeping the parameters identical to the mesoscale-resolving run except for numerical viscosity. As noted earlier, the first deformation radius is around 25 km, so the two resolutions can be considered mesoscale permitting [START_REF] Hallberg | Using a resolution function to regulate parameterizations of oceanic mesoscale eddy effects[END_REF]). The biharmonic viscosities were 4 = (6.25, 31.25) × 10 10 m 4 s -1 respectively. The mean flow and length scale of the spatial filter ( ) were provided by coarse graining them with a 2 × 2 and 4 × 4 box-car filter respectively. While we acknowledge there may be more sophisticated approaches to filter the background flow [START_REF] Aluie | Mapping the energy cascade in the north atlantic ocean: The coarse-graining approach[END_REF][START_REF] Grooms | Diffusion-based smoothers for spatial filtering of gridded geophysical data[END_REF], the box-car filter is the simplest operator that commutes with spatial derivatives, and additional terms owing to non-commutative properties between the filter and derivatives do not arise upon coarse graining the background flow.

We show in Fig. C1 the time mean of the eddy kinetic and potential energies from the two runs at coarser resolutions. Notably, the run with 256 grids and eddy rectification forcing performs better than the highest-resolution eddy model without the forcing (Figs. 4 and C1a,b) with the F . C2: A scatter plot showing the ratio between area integrated † and in the first layer, and its dependence on the biharmonic viscosity. The runs with 256 grids are shown as black crosses and 128 grids as red dots. Values larger than unity indicate that the coarse-resolution sub-grid models are more energetic than the eddies in the mesoscale-resolving CTRL run. diagnosed from the CTRL run and respective mesoscale-permitting sub-grid models plotted against the numerical viscosity. The runs we show in Fig. C1 were taken from the runs with the highest viscosity respectively. As we decrease the viscosity, the level of EKE increases as expected, with the run with 128 grids showing a strong dependency. While the sub-grid model with a prescribed double-gyre background flow could be run stably with small numerical viscosity in respect to its resolution, the poorly resolved instabilities tended to excite Rossby waves in the gyre interior (not shown), which accumulated at the western boundary (the western boundary current is too zonally broad in Fig. C1c). This caused the domain integrated EKE to be larger than that diagnosed from the CTRL run, viz. values larger than unity in