<figure><img src="image_1.png" /></figure>

HAL
open science

# Non-Local Eddy-Mean Kinetic Energy Transfers in Submesoscale-Permitting Ensemble Simulations

Quentin Jamet, Stephanie Leroux, William K Dewar, Thierry Penduff, Julien

Le Sommer, Jean-marc Molines, Jonathan Gula

► To cite this version:

Quentin Jamet, Stephanie Leroux, William K Dewar, Thierry Penduff, Julien Le Sommer, et al. Non-Local Eddy-Mean Kinetic Energy Transfers in Submesoscale-Permitting Ensemble Simulations. Journal of Advances in Modeling Earth Systems, 2022, 14 (10), pp.e2022MS003057. 10.1029/2022MS003057. hal-03835762

HAL Id: hal-03835762

https://cnrs.hal.science/hal-03835762

Submitted on 1 Nov 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
RESEARCH ARTICLE

10.1029/2022MS003057

**Key Points:**

* Ensemble-based eddy-mean decomposition of kinetic energy budget supports the view of an ocean turbulence driven by internal dynamics

* Turbulent fluxes of the cross energy term provide a potentially strong horizontal constraint on eddy-mean flow interactions

* Non-localities are leading order at small scales and should be accounted for in submesoscale parameterizations

**Correspondence to:**

Q. Jamet,
quentin.jamet@univ-grenoble-alpes.fr

**Citation:**

Jamet, Q., Leroux, S., Dewar, W. K., Penduff, T., Le Sommer, J., Molines, J.-M., & Gula, J. (2022). Non-local eddy-mean kinetic energy transfers in submesoscale-permitting ensemble simulations. *Journal of Advances in Modeling Earth Systems*, *14*, e2022MS003057. https://doi.org/10.1029/2022MS003057

Received 24 FEB 2022
Accepted 29 AUG 2022

# Non-Local Eddy-Mean Kinetic Energy Transfers in Submesoscale-Permitting Ensemble Simulations

<figure><img src="image_2.png" /></figure>

Quentin Jamet¹ ID, Stephanie Leroux¹,² ID, William K. Dewar¹,³ ID, Thierry Penduff¹ ID, Julien Le Sommer¹ ID, Jean-Marc Molines¹, and Jonathan Gula⁴,⁵ ID

¹CNRS, IRD, Grenoble INP, IGE, Université Grenoble Alpes, Grenoble, France, ²Ocean Next, Grenoble, France,

³Department of EOAS, Florida State University, Tallahassee, FL, USA, ⁴Laboratoire d'Océanographie Physique et Spatiale (LOPS), CNRS, Ifremer, IRD, IUEM, University Brest, Plouzané, France, ⁵Institut Universitaire de France (IUF), Paris, France

**Abstract** Understanding processes associated with eddy-mean flow interactions helps our interpretation of ocean energetics, and guides the development of parameterizations. Here, we focus on the non-local nature of Kinetic Energy (KE) transfers between mean and turbulent reservoirs. Transfers are interpreted as non-local when the energy extracted from the mean flow does not locally sustain an growth of energy in the turbulent flow, or vice versa. The novelty of our approach is to use ensemble statistics to define the mean and the turbulent flow. Based on KE budget considerations, we first rationalize the eddy-mean separation in the ensemble framework, and discuss the interpretation of a mean flow $\langle u \rangle$ driven by the prescribed (surface and boundary) forcing and a turbulent flow $u'$ driven by non-linear dynamics sensitive to initial conditions. We then analyze 120-day long, 20-member ensemble simulations of the Western Mediterranean basin run at $1/60^\circ$ resolution. Our main contribution is to recognize the prominent contribution of the cross energy term $\langle u_h \rangle \cdot u'_h$ to explain non-local energy transfers, which provides a strong constraint on the horizontal organization of eddy-mean flow KE transfers since the cross energy term vanishes identically for perturbations ($u'_h$) orthogonal to the mean flow ($\langle u_h \rangle$). We also highlight the prominent contribution of vertical turbulent fluxes for energy transfers within the surface mixed layer. Analyzing the scale dependence of non-local energy transfers supports the local approximation usually made in the development of meso-scale, energy-aware parameterizations for non-eddying models, but points out to the necessity of accounting for non-local dynamics in the meso-to-submeso scale range.

**Plain Language Summary** The ocean constantly exchanges energy between its mean and its turbulent reservoirs. However, we are still lacking a clear understanding of eddy-mean flow interactions, which limits our ability to represent them in numerical ocean simulations that require turbulent closures. In particular, it has been recently shown that instabilities of midlatitude jets do not necessarily sustain the growth of turbulent eddies locally. Instead, the energy released by the jet can be transported over significant distances to either sustain turbulence or to reinforce the jet. Here, we analyze model outputs of submesoscale-permitting (horizontal resolution of 1–2 km) ensemble simulations of the Western Mediterranean basin with the view of better understanding this non-local dynamics. Starting from 20 initial conditions perturbed by small, independent perturbations, we analyze the development of the ensemble spread during 120-days long simulations exposed to identical forcing. We investigate the spatiotemporal structure of eddy-mean flow interactions through their kinetic energy expression. Our main contribution is to highlight turbulent fluxes of the cross energy term as a driving mechanism to explain non-local dynamics, a process that need to be accounted for in the development of submesoscale parametrizations.

## 1. Introduction

Meso-scale eddies play a crucial role for the energetic balance of the ocean, providing the main pathway toward dissipative scales (Wunsch & Ferrari, 2004). Understanding how eddies interact with the mean flow thus helps our interpretation of the ocean circulation, and also serves as a basis for the development of robust parameterizations for ocean models. In order to gain insights from the different processes controlling eddy energetics, it is usual and natural to investigate the different terms contributing to the time rate of change of the Eddy Kinetic Energy (EKE) equation (e.g., Dewar & Bane, 1989; Webster, 1961, 1965). From the point of view of parameterization, evaluating the energy levels of meso-scale "eddies" is used to constrain numerical eddy dissipation

© 2022 The Authors. Journal of Advances in Modeling Earth Systems published by Wiley Periodicals LLC on behalf of American Geophysical Union. This is an open access article under the terms of the Creative Commons Attribution-NonCommercial License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited and is not used for commercial purposes.
coefficients, either through mixing length arguments (Cessi, 2008; Eden & Greatbatch, 2008; Jansen et al., 2019) or through Eliassen-Palm eddy stress tensor (Mak et al., 2018; Marshall et al., 2012), thus making dissipative coefficients energy-aware. In this context, the “eddies” are associated with unresolved, sub-grid scale physics that need to be parameterized based on the *mean*, resolved flow. A particularity of eddy-mean kinetic energy transfers lies in the difference in the terms involved in KE budget of the mean and the turbulent flow. That is, changes in the mean flow energetics are subject to the divergence of an eddy stress tensor correlated with the mean flow, while changes in the turbulent flow energetics are subject to a turbulent flux up or down the gradient of the mean flow. Equating the eddy-mean interaction term from these two different perspectives is subject to an assumption of locality, where the energy released by the mean flow at one location is assumed to sustain the growth of eddies at that location (or vice versa for energy backscattering processes). However, recent studies based on Lorenz energy cycles at global (Chen et al., 2014, 2016) and regional (Capó et al., 2019; Kang & Curchitser, 2015) scales have shed light on the strong non-locality of such transfers at small scales. Our interest in this study is to further investigate the spatiotemporal structure of non-local eddy-mean KE transfers by leveraging the recent developments of kilometric-scale resolution ensemble simulations to separate mean and eddies based on ensemble statistics.

An emerging concern for the development of turbulent parameterizations for ocean models is placed on the non-locality of energy transfers. In early work on energy-aware parameterizations for mesoscale turbulence, Cessi (2008) has proposed an improved Gent-McWilliams (Gent & McWilliams, 1990) formulation in which the eddy buoyancy diffusivity was defined as a function of the averaged sub-grid scale turbulent kinetic energy through mixing length arguments. Although globally integrated estimates of sub-grid scale kinetic energy offer interesting properties (Marshall & Adcroft, 2010), it obviously *only* provides an averaged estimate. Other studies have provided more elaborate formulations to account for the spatial organization of mesoscale eddy diffusivity (Ferreira et al., 2005; Groeskamp et al., 2020; Visbeck et al., 1997), but at the expense of severely complicating the prognostic equation of sub-grid scale turbulent kinetic energy that needs to be solved (Eden & Greatbatch, 2008; Jansen et al., 2019; Mak et al., 2018). In practice, the several processes involved in this prognostic equation are usually parameterized through isotropic dissipative operators, mostly due to the lack of better theories. However, Grooms (2017) has recently shown that, while the local approximation is valid for isotropic barotropic turbulence with no mean flow, idealized advection-diffusion models rapidly fail to accurately represent the transport of EKE when a mean flow is present in the problem (arising from the presence of the $\beta$ effect in his case). A potential reason to explain this is associated with the non-locality of the eddy energy transfers, as for instance identified in a wind-driven, two-layer QG model by Grooms et al. (2013); in this simulation, the energy lost by eddies in the separated jet is primarily balanced by imports of energy from remote regions. Non-local kinetic energy reported by Grooms et al. (2013) are associated with various processes, such as wave radiation, advection, or eddy-mean flow interactions. The latter relates the dynamics behind energy transfers between the mean and the turbulent flow, and its leading order contribution has been recently reported by Chen et al. (2014), Kang and Curchitser (2015), and Capó et al. (2019) in realistic simulations. It is thus likely to have important implications for the development of future parameterizations.

There are many ways to define “mean” and “eddies,” the most traditional approach being to use a time averaging. This definition offers several advantages, such as ease in implementation and natural interpretation when dealing with observations. Eddies so defined are however associated with all signals that vary in time, which makes the attribution of processes somehow ambiguous (for instance in disentangling processes associated with hydrodynamic instabilities from those associated with time varying forcing). Coarse-graining or spatial filtering (e.g., Aluie et al., 2018; Grooms et al., 2021, and references therein) offer alternative approaches, which are more intuitive in the context of parameterization. Although the time dimension is retained, such approaches induces some subjectivity in the definition of length scale cutoff, thus the size of the eddies, as well as complexities in dealing with solid boundaries, isotropy and inhomogeneities of the flow structure.

Here, we choose to leverage ensemble simulations to define the “mean” flow as that common to all members (i.e., an ensemble mean), and the “eddies” as the deviation of each member with its ensemble mean. We will argue in the following that this approach offers an unambiguous definition of “eddies” through KE budget considerations; it allows to robustly separate the flow in a part that is controlled by the prescribed forcing (the “mean” flow), and a part that is intrinsically driven by non-linear dynamics (the “eddies”). Ensembles also allows the analysis of the spatiotemporal structure of ocean turbulence and its associated flux of energy. An obvious limitation is associated with the computational resources required to produce such a data set. Here, in order to partially account for the
potential effects of submesoscale dynamics in eddy-mean flow interactions, we have used the newly generated kilometric-scale resolution (1/60°) MEDWEST60 ensemble simulations of Leroux et al. (2021). It is composed of 20 ensemble members subject to small initial conditions uncertainties (usually referred to as micro initial conditions; Stainforth et al., 2007), run for 120-days from the already spun-up oceanic state of eNATL60 simulation (Brodeau et al., 2020), a numerically identical, single simulation run over the whole North Atlantic basin. Analyzing the decorrelation of each ensemble member in this context informs us on the processes controlling the growth of ensemble spread, thus on the spatiotemporal structure of eddy-mean flow interactions.

The paper is organized as follows. In Section 2, we first discuss the eddy-mean decomposition of kinetic energy budget in the context of ensemble simulations, and the present the MEDWEST60 ensemble simulations as well as the diagnostic tools used for their analysis. We then discuss the decorrelation of the turbulent flow from initial conditions, as well as some aspects of the associated kinetic energy budgets in Section 3. In Section 4, we first diagnose the non-local kinetic energy transfers, and then estimate their spatial scale dependence with a view toward parameterization. We finally summarize our results and discuss their implications in Section 5.

## 2. Methods

### 2.1. Kinetic Energy Budget of Ensemble Simulations

Our primary interest is to investigate the kinetic energy budget of the MEDWEST60 submesoscale-permitting ensemble simulations, described in Section 2.2, with a focus on energy transfers between the ensemble mean and the turbulent flow. The momentum equations solved by MEDWEST60 ensemble simulations are the Boussinesq, hydrostatic equations written in flux form:

$$ \partial_t u = -\nabla \cdot \mathbf{u}u + f v - \frac{1}{\rho_0} \partial_x p + \mathbf{D}_u, \quad (1a) $$

$$ \partial_t v = -\nabla \cdot \mathbf{u}v - f u - \frac{1}{\rho_0} \partial_y p + \mathbf{D}_v, \quad (1b) $$

with $\mathbf{u} = (u, v, w)$ the three-dimensional velocity field, $\nabla = (\partial_x, \partial_y, \partial_z)$ the three-dimensional gradient operator, $f = 2\Omega\sin(\phi)$ the Coriolis frequency and $\phi$ the latitude, $p = \int_z^\eta \rho g dz$ the (hydrostatic and surface) pressure field, and $\mathbf{D}_u = \partial_z (\mathbf{A}\partial_z u)$ and $\mathbf{D}_v = \partial_z (\mathbf{A}\partial_z v)$, the viscous effects including both surface wind forcing and bottom drag as surface and bottom boundary conditions, respectively, as well as interior ocean dissipation of momentum, with $\mathbf{A}$ the spatiotemporally varying viscous coefficient computed through the TKE turbulent closure scheme. Horizontal viscous effects are implicitly included in the UBS advective scheme as a biharmonic operator (Shchepetkin & McWilliams, 2005) (see Appendix A for further details).

Following standard practices, an equation for the hydrostatic kinetic energy

$$ K = \frac{\rho_0}{2} (\mathbf{u}_h \cdot \mathbf{u}_h), \quad (2) $$

with $\mathbf{u}_h = (u, v)$ the horizontal component of the velocity field, is obtained by multiplying Equation 1a by $\rho_0 u$ and Equation 1b by $\rho_0 v$, and summing the resulting equations, such that:

$$ \partial_t K = -\nabla \cdot (\mathbf{u}K) - \mathbf{u}_h \cdot \nabla_h p + \partial_z (\mathbf{A}\partial_z K) - \epsilon, \quad (3) $$

with $\nabla_h = (\partial_x, \partial_y)$ the horizontal gradient operator, $\partial_z (\mathbf{A}\partial_z K)$ the work done by vertical viscous forces, and $\epsilon = \rho_0 \mathbf{A} (\partial_z \mathbf{u}_h \cdot \partial_z \mathbf{u}_h)$ the vertical dissipation of kinetic energy. Adding and subtracting $-w \partial_z p = w b$ in Equation 3, and using the continuity equation for Boussinesq fluids $\nabla \cdot \mathbf{u} = 0$, allows the pressure term to be written as the divergence of a flux, and makes explicit the exchange of kinetic energy with potential energy through $wb$:

$$ \partial_t K = -\nabla \cdot (\mathbf{u}K) - \nabla \cdot (\mathbf{u}p) - w b + \partial_z (\mathbf{A}\partial_z K) - \epsilon. \quad (4) $$

In our ensemble simulations, the velocity field simulated by each individual ensemble member obeys this KE equation. It is however possible, from ensemble statistics, to decompose the velocity field as that common to all members, and that specific to each member, and analyze their kinetic energy expression.
For this, we consider the Reynolds decomposition

$$x_n = \langle x \rangle + x'_n, \qquad (5)$$

where the mean operator

$$\langle x \rangle = \frac{1}{N} \sum_{n=1}^{N} x_n. \qquad (6)$$

represents the ensemble mean, with N the size of the ensemble. Following this procedure to decompose the zonal and meridional velocities defining the kinetic energy (Equation 2) leads to:

$$K = \tilde{K} + K^* + \rho_0 \langle \mathbf{u}_h \rangle \cdot \mathbf{u}'_h, \qquad (7)$$

where $\tilde{K} = \frac{\rho_0}{2} (\langle \mathbf{u}_h \rangle \cdot \langle \mathbf{u}_h \rangle)$ and $K^* = \frac{\rho_0}{2} (\mathbf{u}'_h \cdot \mathbf{u}'_h)$. For reasons explained below, we will refer the former quantity ($\tilde{K}$) as the Forced Kinetic Energy (FKE), and the ensemble mean of the latter quantity ($\langle K^* \rangle$) as the Internal Kinetic Energy (IKE), which refer to the kinetic energy of the ensemble mean flow and the kinetic energy of the perturbations, respectively. The notation used here is somehow different from the more classical Mean and Eddy Kinetic Energy (MKE, EKE) terminology used when working with time averages. While these terms are formally the same, the different terminology used here aims at highlighting differences in their interpretation and properties in the context of ensemble simulations. Such differences are further discussed below. Finally, we note that the vector form employed here also emphasizes that, in addition to vanishing identically upon averaging, the cross energy term $\rho_0 \langle \mathbf{u}_h \rangle \cdot \mathbf{u}'_h$ is also zero for turbulent flow orthogonal to the mean flow.

The kinetic energy equation for the mean flow and that for the perturbations are usually derived based on averaged and residual forms of Equation 1a and 1b. Formally, multiplying the ensemble mean equations $\langle (1a) \rangle$ and $\langle (1b) \rangle$ by the ensemble mean zonal and meridional velocities $\rho_0 \langle \mathbf{u} \rangle$ and $\rho_0 \langle \mathbf{v} \rangle$, respectively, and summing the resulting equations, leads to an equation for the Forced Kinetic Energy (FKE) of the form:

$$\partial_t \tilde{K} = -\nabla \cdot (\langle \mathbf{u} \rangle \tilde{K}) - \rho_0 \langle \mathbf{u}_h \rangle \cdot \nabla \cdot (\langle \mathbf{u}' \otimes \mathbf{u}'_h \rangle) - \nabla \cdot (\langle \mathbf{u} \rangle \langle \mathbf{p} \rangle) - \langle \mathbf{w} \rangle \langle \mathbf{b} \rangle + \partial_z (\langle \mathbf{A} \rangle \partial_z \tilde{K}) - \epsilon_{\tilde{K}}, \qquad (8)$$

where $\mathbf{u}' \otimes \mathbf{u}'_h = \mathbf{u}' \mathbf{u}'_h^T$ represents the outer product of the three-dimension velocity field $\mathbf{u}'$ with its horizontal component $\mathbf{u}'_h$, with $\mathbf{u}'_h^T$ the transpose of the latter. The first term on the RHS of Equation 8 is associated with the advection of FKE by the mean flow, and the underlined term is associated with eddy-mean flow interactions. Their respective contribution for the time rate of change of FKE ($\partial_t \tilde{K}$) will be further evaluated in Section 3. The exchange of FKE with forced potential energy is made explicit through the inclusion of $\langle \mathbf{w} \rangle \langle \mathbf{b} \rangle$.

A similar equation is obtained for the Internal Kinetic Energy (IKE) by multiplying the residual equation for the zonal and meridional momentum (Equation 1a)' and (Equation 1b)' by the zonal and meridional velocity perturbations $\rho_0 \mathbf{u}'$ and $\rho_0 \mathbf{v}'$, ensemble averaging and then summing the resulting equations, leading to:

$$\partial_t \langle K^* \rangle = -\nabla \cdot (\langle \mathbf{u} K^* \rangle) - \rho_0 \langle \mathbf{u}' \otimes \mathbf{u}'_h \rangle \cdot \nabla \langle \mathbf{u}_h \rangle - \nabla \cdot (\langle \mathbf{u}' \mathbf{p}' \rangle) - \langle \mathbf{w}' \mathbf{b}' \rangle + \partial_z \langle \mathbf{A}' \partial_z K^* \rangle - \epsilon_{K^*}, \qquad (9)$$

where the first term on the RHS of Equation 9 includes advection of IKE by both the ensemble mean and the turbulent flow, and the underlined term is associated with eddy-mean flow interactions. Again, the exchange of IKE with internal potential energy is made explicit through the inclusion of $\langle \mathbf{w}' \mathbf{b}' \rangle$. The respective contribution of these tree terms for the time rate of change of IKE ($\partial_t \langle K^* \rangle$) will be further evaluated in Section 3. The sum of Equation 8 and 9 leads to an equation for the ensemble mean kinetic energy of the full flow, that is, $\partial_t \langle K \rangle = \partial_t \tilde{K} + \partial_t \langle K^* \rangle$.

Another, yet equivalent, procedure to derive an equation for the ensemble mean kinetic energy of the full flow consists in expanding the different components of Equation 4 following the Reynolds decomposition in the ensemble dimension (Equation 5), then ensemble averaging, leading to:

$$\begin{aligned} \partial_t \langle K \rangle = & -\nabla \cdot (\langle \mathbf{u} \rangle \tilde{K}) - \nabla \cdot (\langle \mathbf{u} K^* \rangle) - \rho_0 \nabla \cdot (\langle \mathbf{u}' \rangle (\langle \mathbf{u}_h \rangle \cdot \langle \mathbf{u}'_h \rangle)) \\ & -\nabla \cdot (\langle \mathbf{u} \rangle \langle \mathbf{p} \rangle) - \nabla \cdot (\langle \mathbf{u}' \mathbf{p}' \rangle) - \langle \mathbf{w} \rangle \langle \mathbf{b} \rangle - \langle \mathbf{w}' \mathbf{b}' \rangle + \partial_z (\langle \mathbf{A} \rangle \partial_z \tilde{K}) + \partial_z (\langle \mathbf{A}' \partial_z K^* \rangle) - \epsilon_{\tilde{K}} - \epsilon_{K^*}, \end{aligned} \qquad (10)$$
where $\epsilon_{\bar{K}}$ and $\epsilon_{K^*}$ represents dissipation of FKE and IKE, respectively. Here, the underlined term emerged from the advection of the cross energy term $\langle \mathbf{u}_h \rangle \cdot \mathbf{u}_h'$ by the perturbations. It reflects that, although the covariance of eddy and mean velocity field vanishes identically upon averaging, its advection by perturbations does not. This is of particular interest because it is associated with kinetic energy transfers between the mean and the turbulent flow, thus plays a critical role in eddy-mean flow interactions. Indeed, following the chain rule, the underlined term in Equation 10 can be decomposed as

$$ \underbrace{-\nabla \cdot \langle \mathbf{u}' (\langle \mathbf{u}_h \rangle \cdot \mathbf{u}_h') \rangle}_{\text{DIVEF}} = \underbrace{-\langle \mathbf{u}_h \rangle \cdot \nabla \cdot \langle \mathbf{u}' \otimes \mathbf{u}_h' \rangle}_{\text{MEC}} - \underbrace{\langle \mathbf{u}' \otimes \mathbf{u}_h' \rangle \cdot \nabla \langle \mathbf{u}_h \rangle}_{\text{EDDYFLX}} \quad (11) $$

where the continuity equation has been used to express the last term of the RHS of Equation 11 in a more conventional way. (Note that the LHS of Equation 11 can be formally expressed with tensor notations as $\nabla \cdot (\langle \mathbf{u}' \otimes \mathbf{u}_h' \rangle \cdot \langle \mathbf{u}_h \rangle)$). The first term of the RHS of Equation 11 is the covariance of the horizontal mean flow with the divergence of the Reynolds stress tensor associated with the FKE equation, and the second term of the RHS of Equation 11 is the eddy momentum fluxes up or down the gradient of the mean flow associated with the IKE equation. Expanding the underlined term in Equation 10 as Equation 11 then leads to an equation for the ensemble mean kinetic energy of the full flow that equates the sum of the FKE and the IKE equation, that is, Equation 8 and 9. In the following, we will refer to the three terms of Equation 11, from left to right, as DIVergence of Eddy Flux (DIVEF), Mean-to-Eddy energy Conversion (MEC), and EDDY momentum FLUX (EDDYFLX). A detailed analysis of their spatiotemporal structure is presented in Section 4.

By volume integration, several components of Equation 10 become statements about fluxes at the boundaries of the volume of integration through the divergence theorem. In ensemble simulations such as those we analyze here, ocean surface and boundary conditions are usually prescribed as ensemble mean conditions, common to all members, such that we can neglect turbulent fluxes at the (surface and open) boundaries. (This assumption, along with bottom turbulent fluxes, are further discussed in Section 2.3.) Several terms of the integrated version of Equation 10 thus vanish, and the domain integrated equation for the ensemble mean kinetic energy of the full flow simplifies to:

$$ \begin{aligned} \partial_t \int_V \langle K \rangle dV = \partial_t \int_V \bar{K} dV + \partial_t \int_V \langle K^* \rangle dV = & \\ & - \int_S (\langle \mathbf{u} \rangle \bar{K}) \cdot \mathbf{n} dS - \int_S (\langle \mathbf{u} \rangle \langle p \rangle) \cdot \mathbf{n} dS - \int_V (\langle \omega \rangle \langle b \rangle + \langle \omega' b' \rangle) dV & (12) \\ & + \int_A (\langle \mathbf{u}_h \rangle \cdot \langle \boldsymbol{\tau} \rangle) dA - \int_B (\langle \mathbf{u}_h \rangle \cdot \langle \mathbf{F} \rangle) dB - \int_V (\epsilon_{\bar{K}} + \epsilon_{K^*}) dV, \end{aligned} $$

where V is the volume of integration, S the surface bounding V, A, and B its ocean surface and bottom part, respectively, and n the normal to the surface S. Here, the work done by surface wind stress and bottom friction ($\int_A (\langle \mathbf{u}_h \rangle \cdot \langle \boldsymbol{\tau} \rangle) dA$ and $\int_B (\langle \mathbf{u}_h \rangle \cdot \langle \mathbf{F} \rangle) dB$ with F the vertical diffusive flux at the bottom boundary, respectively) comes from the volume integration of viscous forces. The time rate of change of kinetic energy within the domain thus reflects the import/export of FKE and the wavefield prescribed at the open boundaries (two first terms), exchanges with potential energy (third term), work associated with prescribed surface forcing (fourth term) and bottom boundary condition (fifth term), and dissipation (last term). We note here that although the transfers of kinetic energy between the mean and the turbulent flow (underlined term in Equation 10) can be locally large, they cancel each other when integrated over the entire basin to satisfy the boundary condition of no turbulent flux of the LHS of Equation 11.

The turbulent version of Equation 12 summarizes as:

$$ \partial_t \int_V \langle K^* \rangle dV = -\rho_0 \int_V (\langle \mathbf{u}' \otimes \mathbf{u}_h' \rangle \cdot \nabla \langle \mathbf{u}_h \rangle) dV - \int_V \langle \omega' b' \rangle dV - \int_V \epsilon_{K^*} dV, \quad (13) $$

where the first term of the RHS of Equation 13 comes from the development of Equation 11. In a basin integrated sense, the time rate of change of IKE as diagnosed through ensemble statistics is thus a balance between exchanges with FKE, exchanges with eddy potential energy, and dissipation (horizontal and vertical component, which are treated as residual when interpreting numerical results, see Section 2.3). It is not directly driven by prescribed forcing, but rather reflects the part of the ocean intrinsic dynamics that develops spontaneously in response to the non-linearity of the system. This provides an energy-budget based rationalization that the ensemble strategy
<figure><img src="image_3.png" /><figcaption>**Figure 1.** (top panels) Vertically integrated MEC ($-\rho_0 \langle \mathbf{u}_h \rangle \cdot \nabla \cdot \langle \mathbf{u}' \otimes \mathbf{u}_h' \rangle$, left panel), EDDYFLX ($-\rho_0 \langle \mathbf{u}' \otimes \mathbf{u}_h' \rangle \cdot \nabla \langle \mathbf{u}_h \rangle$, center panel), and DIVEF ($-\rho_0 \nabla \cdot \langle \mathbf{u}' \cdot (\langle \mathbf{u}_h \rangle \cdot \mathbf{u}_h') \rangle$, right panel) after 60 days of simulation. Their volume integrated values are shown at the bottom right of each panels. (bottom panel) Basin integrated MEC (black), EDDYFLX (red) and DIVEF (blue) for each individual members.</figcaption></figure>

provides an unambiguous definition of the ocean *turbulence*. In the following, we pay a particular attention to the contribution of EDDYFLX for the construction of IKE, and its relation to the mean flow (MEC) through the flux divergence DIVEF.

## 2.2. Model and Simulations

We analyze in this study a subset of the MEDWEST60 ensemble simulations (Leroux et al., 2021). These simulations have been produced to evaluate the predictability of the fine scale dynamics in a typical high-resolution Copernicus Marine Environment Monitoring Service (CMEMS) forecasting model by including the effect of initial and model uncertainties. They are based on a kilometric-scale regional configuration of the Western Mediterranean sea (cf., Figure 1) that uses the same numerical choices as the North Atlantic simulation eNATL60 (Brodeau et al., 2020). Briefly, they are NEMO-v3.6 simulations run at $\frac{1}{60}^\circ$ and with vertical grid spacing of 1 m at the surface and 24 m at depth, for a total of 212 vertical levels in MEDWEST60. The simulations are forced at the surface with 3-hourly ERA-interim (ECMWF) atmospheric reanalysis through the CORE bulk flux formulation (Large & Yeager, 2004), and they partially account for surface ocean current feedbacks (e.g., Renault, Molemaker, McWilliams, et al., 2016), where only 50% of surface currents speed is considered in the computation of the wind stress. The tuning is based on Julien Jouanno's recommendations who performed sensitivity tests on modeled EKE levels with (i.e., 100%) and without (0%) ocean current feedbacks in wind stress formulation, and found 50% as a good compromise to reproduce the level of EKE observed by satellite altimetry. Open boundary conditions are applied at the eastern and western boundaries of the domain with a Flow Relaxation Scheme (FRS) for baroclinic velocities and active tracers (Davies, 1976; Engedahl, 1995) and Flather (1994) radiation scheme for sea-surface height and barotropic velocities. The former is a simple relaxation of model fields toward hourly, externally-specified values over the 12 grid points adjacent to the boundaries. The relaxation time scale ranges from $\tau = 0$ s at the domain edge and increases exponentially to about 30 days at grid point 12. The latter ("Flather") applies radiation conditions on the normal depth-mean transport across the open boundaries, set as prescribed values plus a correction based on sea surface height anomalies at the boundaries that allows
gravity waves generated within the domain to exit through the open boundaries. We note that the prescribed boundary conditions are taken from the eNATL60 North Atlantic experiment run with tidal forcing, such that MEDWEST60 includes tides through boundary conditions in addition to tidal potential forcing.

Among the various ensemble simulations produced in the context of MEDWEST60, we focus here on the 20-member ensemble ENS-CI-GSL19, which has been produced as follows. From the already spun-up (through a 18 months integration) oceanic state of the eNALT60 simulation at 5 February 2010, an ensemble of 20 runs has been produced for 1 day with a stochastic perturbation (Brankart et al., 2015) applied on the horizontal grid of the model to represent uncertainties affecting the smallest scales in the model (for more details, see Leroux et al., 2021). The 20 oceanic states so generated have then been used as initial conditions for the production of a 120-day long, 20-member ensemble where all other components of the simulation (including forcing) are common across all members, and the stochastic perturbations are turned off. Such a procedure is usually referred to as micro initial condition uncertainties (Hawkins et al., 2016; Stainforth et al., 2007), and is meant to allow the growth of dynamically consistent small perturbations.

## 2.3. Diagnostic Considerations

During the production of MEDWEST60 ensemble simulations, prognostic variables of the model (T, S, U, V, SSH), as well as vertical velocity (W), have been saved every hour. Based on hourly averaged model outputs, we have used offline diagnostic tools to recompute the kinetic energy budget of MEDWEST60 simulations by closely following the numerical implementations of NEMO. Relevant details for the present analysis are provided in Appendix A, along with validation. These offline tools, along with the high frequency of model outputs (hourly), provide us with a reliable procedure to accurately (relative errors ~ $\mathcal{O}(10^{-3})$, see Table A1) compute the kinetic energy trends due to advection, thus the terms associated with eddy-mean kinetic energy transfers.

In our kinetic energy budget considerations derived in Section 2.1, we have assumed zero turbulent fluxes conditions at the boundaries of the domain. In practice, however, the computation of surface wind stress partially (50%) accounts for ocean-atmosphere feedback (Renault, Molemaker, McWilliams, et al., 2016), such that the turbulent wind work $\langle \mathbf{u}'_h \cdot \boldsymbol{\tau}' \rangle$ is not strictly zero. Its contribution is however weak ($-0.12$ TJ; $1$ TJ = $10^{12}$ J) as compared to mean wind work ($+5.10$ TJ) over the course of the 120-day long simulation, and is several orders of magnitude smaller than the total IKE production of $+2.27$ PJ ($1$ PJ = $10^{15}$ J) within the domain. Furthermore, turbulent wind work is negative, providing a sink for domain integrated IKE time rate of change, in agreement with the eddy-killing effect (Renault, Molemaker, Gula, et al., 2016). Similar considerations are also relevant for turbulent bottom stress, which damps the production of IKE. Our estimates of surface and bottom velocities ensemble spread suggest the bottom contribution is at least one order of magnitude weaker than the surface contribution. As for the open boundary conditions, the “Flather” scheme allows gravity waves generated within the domain to exit the model through boundaries, thus providing an explicit sink of IKE. In an averaged sense, all members are however expected to exhibit similar levels of energy associated with the development of such waves, such that the spread so induced on model velocities is expected to be weak and can be neglected. We recall that baroclinic velocities are strongly relaxed toward prescribed values at the boundaries. The contribution of surface and boundary turbulent forcing, as well as bottom turbulent stress, for the interpretation of IKE production in our ensemble can then be safely neglected.

Finally, we are primarily interested in diagnosing eddy-mean flow kinetic energy transfers through DIVEF, MEC and EDDYFLX (cf., Equation 11). As detailed above, open boundary conditions ensure that the ensemble spread at the boundaries is controlled, such that the domain integrated eddy fluxes of the cross energy term $\rho_0 V \cdot \langle \mathbf{u}' (\langle \mathbf{u}_h \rangle \cdot \mathbf{u}'_h) \rangle$ is negligible. This implies all the energy released by the ensemble mean flow has been used to sustain the growth of IKE within the domain, which we have tested by computing the volume integrated MEC and EDDYFLX for the full domain, and estimating their divergence DIVEF. We show on top panels of Figure 1 the vertically integrated MEC and EDDYFLX, and their divergence (DIVEF) is obtained by simple summation following Equation 11. Integrated over the full domain, MEC drain $-0.53$ GW of energy out of the ensemble mean flow at that particular time (day 60), and EDDYFLX supply $+0.58$ GW of energy to the turbulent flow. The close balance confirms that our procedure provides reliable estimates of these fluxes, with a ~10% error. The error, of about 0.05 GW, is relatively constant across the 20 ensemble members ($\pm 0.01$ GW, Figure 1, lower panel), suggesting a systematic error in our estimates. We attribute the error to the implicit dissipation of the UBS advective scheme used in MEDWEST60. As detailed in Appendix A, we have performed the
<figure><img src="image_4.png" /><figcaption>**Figure 2.** (Upper panels) Spatial maps of surface currents ensemble mean kinetic energy of the full flow (⟨K⟩; left), kinetic energy of the ensemble mean flow (K̃, FKE; center) and the ensemble mean kinetic energy of the turbulent flow (⟨K*⟩, IKE; right) after 60 days of simulation. (Lower panel) 120-day long time series of these quantities, integrated within the green box. The time series of the kinetic energy of a given member is provided for reference (gray line). Units of the spatial maps are J m⁻³ and those of the time series are terrajoules (1 TJ = 10¹² J). The black box on top left panel is used to validate our recomputation of kinetic energy budgets (cf., Appendix A).</figcaption></figure>

eddy-mean flow decomposition of the advective operator based on a fourth order centered scheme, which is the non-dissipative equivalent of the UBS scheme. The error in our estimates being positive and relatively constant across ensemble members suggests it is associated with dissipation.

In the following sections, we turn our attention to the analysis of the MEDWEST60-ENS-CI-GSL19 ensemble simulations, where we first diagnose the decorrelation of the turbulent flow from its ensemble mean, then evaluate the respective contribution of MEC and EDDYFLX for the kinetic energy budget of the ensemble mean and the turbulent flow, and then analyze their interactions through DIVEF.

# 3. Results

## 3.1. Decorrelation of the Turbulent Flow

Figure 2 provides horizontal maps and time evolution of surface kinetic energy, as well as its ensemble statistical decomposition. From left to right, the upper panels show the ensemble mean surface kinetic energy of the full flow ⟨K⟩, the FKE and the IKE at day 60. Their time evolution over the course of the 120 days, integrated within the green box, are shown on the lower panel. The ensemble mean full kinetic energy ⟨K⟩ exhibits a combination of high and low frequency variations, but remains relatively constant (6–8 TJ; 1 TJ = 10¹² J) over the 120 days, reflecting the already spun-up state of the eNATL60 simulation used to initialize the ensemble. For reference, the level of kinetic energy of a given member is shown in light gray. It exhibits small variations around its ensemble mean equivalent, illustrating that the ensemble mean kinetic energy of the full flow provides a statistical estimate of the energy level of the ensemble. We note that the deviation of the kinetic energy of a single member from the ensemble mean kinetic energy is not to be confused with the separation between the kinetic energy of the ensemble mean flow and that of the perturbations, which is the primary focus of our study.

The spatial pattern of the FKE (K̃) is representative of the relatively well organized flow within the western Mediterranean basin. In the northern half, the FKE exhibits high levels of energy associated with the southwestward flowing Liguro-Provençal current (Millot, 1999; Waldman, 2016). In the southern half, FKE exhibits a very large
import of energy through the strait of Gibraltar (exceeding 2,000 J m⁻³), the development of standing eddies downstream, and an eastward flowing boundary current along the southern boundary of the basin (the Algerian Current, Millot, 1985). Around 5°E, the Algerian Current detaches from the coast, forming a “loop current,” a region of intense meso-scale eddies formation through mixed baroclinic-barotropic instabilities (e.g., Obaton et al., 2000; Poulain et al., 2021). We will focus on the eddy dynamics of this region in the following. Although IKE (⟨K*⟩) is more pronounced in the southern than in the northern part of the domain, it somehow follows the spatial organization of FKE, reflecting the link between the two; turbulent dynamics develop in region of strong currents, which are more prone to instabilities.

The lower panel of Figure 2 illustrates the time evolution of surface FKE and IKE, integrated within the green box, during the 120 days of simulation. At the beginning all ensemble members are in phase, such that IKE is zero and FKE reflects the energy content of the full flow. The latter diverges from the ensemble mean full KE about 1 week after initialization as each ensemble member starts to decorrelate. At the end of the 120 days, FKE has dropped to less than 2 TJ, that is, about one third of its initial energy content. In the same time, the turbulent part of the flow (IKE, ⟨K*⟩) develops and reaches about 5 TJ at the end of the 120 days. The development of IKE exhibits several stages before saturation at about day 80. It is interesting to note that a first increase in IKE is observed from day 6 to day 20, where IKE reaches a first plateau. The 6 days time scale for the turbulent flow to start decorrelating from initial conditions is consistent with time scale reported by Fox-Kemper et al. (2008) and Schubert et al. (2020) in their idealized linear study of mixed layer instability and absorption of submesoscale vortices by mesoscale eddies, respectively. In both studies, time scales shorter than 1 week are associated with the development of submesoscale structures through surface mixed layer instabilities, which then saturate and undergo non-linear interactions to transfer their energy upscale. The 6 days time scale in our ensemble simulations is thus likely associated with similar processes, and suggests non-linear interactions of submesoscale instabilities are responsible for the initial growth of IKE. The other stages of IKE increase are associated with further development of turbulent flow. By comparing the IKE patterns at days 30 and 60 for instance (not shown), it appears that initial IKE development mostly takes place along the mean current, while later on, turbulent structures develop more broadly, contributing to the increase in the integrated IKE level within the green box. Additional spectral estimates of the decorrelation of ensemble members over the first 60 days can be found in Leroux et al. (2021). In what follows, we will focus our analysis on day 60, which is about 20 days before the saturation of IKE. As shown in the following, day 60 exhibits a well organized spatial structure in the eddy-mean flow KE interactions that nicely illustrates non-local processes. Such processes are nonetheless observed all along the 120-day long simulation. The interested reader is referred to the following animation: https://doi.org/10.5281/zenodo.6221153. The 120 days of simulation cover the period 6 February–5 June, and a weakened submesoscale activity associated with spring time is observed toward the end of the simulation. It is thus likely such a seasonal cycle will imprint onto eddy-mean flow kinetic energy transfers, a signature observed for instance by Uchida et al. (2022). The relatively short time duration of MEDWEST60 ensemble does however not allow us to quantify such seasonality.

## 3.2. Kinetic Energy Budget

We now turn our attention to the respective contributions of the advective terms of the FKE and IKE budget, focusing on the “loop current” region. We recall here that many other processes contribute to these budgets, such as wave radiation, dissipation or exchanges with turbulent potential energy (cf., Equation 10). We briefly discuss the contribution of the latter in what follows, but otherwise postpone the analysis of other contributions for further work. Here, we focus our attention on the terms driving kinetic energy transfers between the mean and the turbulent flow. We first discuss the kinetic energy budget of the mean flow and that of the turbulent flow, and estimate the respective contribution of MEC and EDDYFLX.

We show on Figure 3 the vertically integrated time rate of change of FKE (top left panel), as well as advection of FKE by the mean flow (−∇·(uK̃); top right panel) and Mean-to-Eddy Conversion (MEC, −ρ₀⟨u_h⟩·∇·⟨u' ⊗ u'_h⟩; bottom left panel) at day 60. Their vertical distributions within the upper 500 m, horizontally integrated within the green box, appear on the bottom right panel as black, blue and red lines, respectively. Note that all horizontal maps have been integrated down to the ocean floor for consistency, but most of the dynamics is observed within the upper 500 m. The contribution from other processes, such as pressure work, surface forcing and viscous effects, as well as small uncertainties associated with our offline estimates (cf., Appendix A), are shown in green as a residual. We first note that the time rate of change of FKE is dominated by a wave-like horizontal structure,
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Vertically integrated time rate of change of FKE (upper left panel), advection of FKE by the mean flow (upper right panel) and Mean-to-Eddy energy Conversion rate (MEC, lower left panel) in the region of the loop current at day 60, with their volume integrated values within the green box shown at the bottom right of each panels. The vertical distribution of these quantities, within the upper 500 m and horizontally integrated within the green box, are shown on the bottom right panel. The other components of the FKE budget, **including viscous effects**, are shown as a residual (green line).</figcaption></figure>

which exhibits a strong baroclinic signature. The fast (daily) evolution of this signal (not shown) suggests it is associated with the high frequency signal observed in the FKE time series of surface currents (Figure 2, bottom panel). As part of the ensemble mean flow, this signal is likely associated with the forcing, such as high frequency winds and, to a smaller extent, tidal forcing. The time rate of change of FKE integrated within the green box is +0.30 GW. In contrast, both advection of FKE by the mean flow and MEC exhibit very different patterns with smaller scale structures. The former exhibits a multipole-like organization, and has an opposite signature in the upper 50 m (i.e., deeper than the ensemble mean and spatially averaged mixed layer depth of about 30 m) than in the rest of the water column. When integrated over the volume however, mean advection of FKE is two orders of magnitude weaker than the volume integrated time rate of change of FKE. Although MEC exhibit weaker signals locally, its volume integrated contribution is significant (–0.24 GW), with a maximum at about 40 m depth.

Figure 4 shows the equivalent of Figure 3 but for the IKE budget. We first note that the spatial pattern of IKE time rate of change is significantly different from that of FKE, with smaller scale structures. Contribution of advection of IKE by the mean and turbulent flow within the box is weak (+0.03 GW), but exhibits local important contributions for the IKE redistribution. EDDYFLX contribute to +0.25 GW to the budget, which slightly exceeds the time rate of change of IKE of +0.21 GW. The vertical profile of turbulent potential to kinetic energy conversion rate – $\langle w'b' \rangle$ is also shown, with a net contribution within the green box of about +0.20 GW. It is maximum at about 30 m depth and tends toward zero at the surface. Although relatively weak when integrated within the green box (–0.08 GW), the large intensification of the residual near the surface is expected to mostly reflect the action of vertical viscous forces and dissipation.
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Same as Figure 3, but for the IKE budget. The advection of IKE (upper right panel) includes advection by both the mean flow (−∇·⟨u⟩⟨K*⟩) and the turbulent flow (−∇·⟨u'K*⟩). Turbulent potential to kinetic energy conversion rate (−⟨w'b'⟩) is also shown in gray and its net contribution within the green box is of about +0.20 GW. Note the change in amplitude of the colorbar as compared to Figure 3.</figcaption></figure>

Finally, we quantify the contribution of EDDYFLX for construction of the IKE over the course of the 120 days of simulations, and assess its relation with the loss of energy of the mean flow through MEC by computing the volume integrated contribution of both EDDYFLX and MEC within the green box of Figure 2 for the 120 day long simulations. We show on Figure 5 the time series of the two contributions (left panel), as well as their time

<figure><img src="image_7.png" /><figcaption>**Figure 5.** (left) Time series of volume integrated MEC (green), EDDYFLX (red), DIVEF (light blue) and −⟨w'b'⟩ (gray) within the green box of Figure 2, and (right) their time integrated contribution. The 120-day long integrated MEC (EDDYFLX, DIVEF, −⟨w'b'⟩) contribution is −2.12 PJ (+2.41 PJ, +0.30 PJ, +1.38 PJ).</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Vertically integrated MEC ($-\rho_0 \langle \mathbf{u}_h \rangle \cdot \nabla \cdot \langle \mathbf{u}' \otimes \mathbf{u}'_h \rangle$, left panel) EDDYFLX ($-\rho_0 \langle \mathbf{u}' \otimes \mathbf{u}'_h \rangle \cdot \nabla \langle \mathbf{u}_h \rangle$, middle panel) and DIVEF ($-\rho_0 \nabla \cdot \langle \mathbf{u}' \otimes (\langle \mathbf{u}_h \rangle \cdot \mathbf{u}'_h) \rangle$, right panel) after 60 days of simulations within the loop current region. Integrated quantities within the green box are shown on the bottom right insert. Ensemble mean surface currents are shown with arrows, and the black line is the section shown in Figure 7.</figcaption></figure>

integrated estimates (right panel). Starting from zero at the beginning of the simulations where all ensemble members are in phase, EDDYFLX start to inject energy in the turbulent flow after about 5–6 days, in agreement with surface IKE increase discussed in Section 3.1. The rate at which EDDYFLX inject energy in the turbulent flow is of about 0.2 GJ s⁻¹ with time variations as large as ±0.13 GJ s⁻¹. MEC drain energy out of the mean flow with similar rate and temporal variations, leading to a small contribution of DIVEF (light blue line). Over the course of the 120 days of simulation, EDDYFLX and MEC have contributed to +2.41 PJ and −2.12 PJ for the IKE and FKE budget, respectively (Figure 5, right panel). The integrated contribution of DIVEF is small within this region, suggesting that eddy-mean energy transfers associated with the loop current instabilities are mostly local. Also shown on this figure is the contribution of the turbulent potential to kinetic energy conversion rate −⟨w′b'⟩. We first note the very large temporal variations in this term as compared to eddy-mean flow interaction processes, suggesting intense exchanges with turbulent potential energy reservoirs on very short time scales. Its time integrated contribution, however, is of the same order of magnitude than EDDYFLX but slightly weaker, supporting mixed barotropic-baroclinic instability processes for driving the growth of Algerian Eddies as proposed earlier (Obaton et al., 2000; Poulain et al., 2021). It is interesting to compare these estimates to the total IKE and FKE changes. During the 120 days of simulation, the volume integrated IKE within the green box has grown by +0.98 PJ, which is only about a quarter of the total energy injected by EDDYFLX and −⟨w′b'⟩. Similarly, the FKE destruction over the full simulation is −0.91 PJ, which is about half of the energy drained by MEC, highlighting the leading order contribution of other processes for balancing kinetic energy budgets of this region.

# 4. Non-Locality of FKE-IKE Energy Transfers

## 4.1. Diagnosing Non-Local KE Transfers

The patterns and amplitude of MEC and EDDYFLX discussed in the previous section are associated with energy transfers between the mean and the turbulent flow. As discussed in the Introduction and in Section 2.1, eddy-mean flow interactions can either be local, with a negligible contribution of DIVEF (left-hand side of Equation 11), or non-local, with transfers of energy with turbulent processes of remote regions. Dynamically, the amplitude of DIVEF provides an estimate of the level of energy released by the mean flow that *locally* sustains the growth of eddies. Or, vice versa, an estimate of the level of energy released by the eddies that is *locally* backscattered to energize the mean flow. We further analyze this local versus non-local contribution in what follows.

Horizontal maps of vertically integrated MEC, EDDYFLX and DIVEF are shown in Figure 6 at day 60, and their volume integrated values within the green box appear at the bottom right of each panel. Averaged over the box, the energy lost by the mean flow (MEC, −0.24 GW) is used to support eddy growth (EDDYFLX, +0.25 GW), and the divergence of eddy flux is weak (DIVEF, +0.01 GW). That MEC drain −2.12 PJ out from FKE and EDDYFLX inject +2.41 PJ into IKE during the 120 days of simulation, as diagnosed in Section 3.2, also supports the interpretation of a turbulence controlled by local processes in this region. However, the details of these energy transfers are complex, and the radically different spatial structure of MEC and EDDYFLX strongly suggests
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Associated vertical structure of MEC, EDDYFLX, and DIVEF along the cross-stream section of Figure 6. Gray contours represent the ensemble mean current across the section. Dashed green contours on middle and right panels show the main structure of MEC.</figcaption></figure>

eddy-mean flow kinetic energy transfers are non-local at small scales. The spatial scale dependence of non-local KE transfers is further analyzed in Section 4.2.

At day 60, the horizontal structure of MEC (Figure 6, left panel) exhibit alternation of FKE destruction (blue spots) with FKE production (red spot), which tend to organize mostly along the mean flow. In contrast, EDDY-FLX (Figure 6, middle panel) exhibit signals of weaker amplitude, which tend to be more pronounced on the flanks of the flow. This suggests a significant part of the kinetic energy lost by the mean flow at one location is advected further downstream before being re-injected in the mean flow, but little is used to sustain the growth of eddies locally. The connection between MEC and EDDYFLX involves DIVEF, which is associated with eddy flux divergence of the cross energy term $\langle \mathbf{u}_h \rangle \cdot \mathbf{u}_h'$. This term exhibits a rich spatial organization (Figure 6, right panel), with regions of destruction of FKE associated with a divergence of eddy flux, that is, the cross energy term is fluxed out of the control volume by the turbulent flow, and regions of FKE production associated with a convergence of eddy fluxes, that is, the cross energy term is fluxed within the controlled volume by the turbulent flow. The region indicated by the black line is of particular interest because it exhibits a region of production of IKE (red spot of EDDYFLX) to the northeast of the region of FKE destruction. MEC, EDDYFLX and DIVEF vertical cross sections along this line are shown in Figure 7. At the surface, MEC exhibit largest negative values about 10 km away from the core of the mean current, and exhibits a tilted vertical structure. In contrast, the EDDYFLX are largest about 20 km northeastward of the minimum of MEC, a region of strong horizontal mean flow gradient, but exhibits a shallower vertical penetration as compared to MEC. As a result, DIVEF are dominated by a divergence of eddy flux near the core of the mean flow, and a convergence on its flank. Although a direct interpretation of a turbulent flux of the cross energy term $\langle \mathbf{u}_h \rangle \cdot \mathbf{u}_h'$ to connect regions of FKE destruction with regions of IKE production is tempting, we recall here that this term vanishes identically for turbulent flow orthogonal to the mean flow. This suggests that DIVEF is more efficient at transporting energy in the along stream direction than in the across stream direction, providing a strong horizontal constraint for eddy-mean flow interactions. This may well provide a dynamical rationalization to explain the large variations of MEC observed in the along stream direction, where energy extracted from the mean flow would be transported downstream before to be reinjected into the mean flow, but little would actually be transferred to the turbulent flow through EDDYFLX.

Figure 8 shows the horizontal and vertical contribution for the three components involved in eddy-mean flow kinetic energy transfers in the upper ocean layer. We first note that, as expected, vertical fluxes are much weaker than horizontal fluxes. However, while weak at each location, vertical turbulent fluxes are predominately positive in the upper layer, such that their horizontally integrated contribution is of the same order of magnitude as the horizontal turbulent fluxes for the three terms (Figure 9). More interestingly, while the horizontal component of MEC and EDDYFLX tend to oppose each other, the vertical components tend to have the same sign. Indeed, the horizontal contribution of MEC are relatively constant and negative in the upper 100 m and smoothly decreases further below (left panel), while the horizontal contribution of EDDYFLX is negligible at the surface, reaches its maximum at about 30 m and smoothly decreases further below (center panel). In contrast, in both MEC and
<figure><img src="image_10.png" /><figcaption>**Figure 8.** (left) Upper layer MEC, (center) EDDYFLX, and (right) DIVEF at day 60, decomposed into an horizontal (top panels) and a vertical (bottom panels) contribution. Ensemble mean surface currents are shown with arrows.</figcaption></figure>

EDDYFLX, vertical turbulent fluxes are upward in the upper 15 m, reach a maximum downward contribution at the base of the spatially averaged mixed layer (about 30 m), and decrease further below to reach negligible contribution below about 100 m. The balanced DIVEF within the green box (right panel) thus results in a balance between horizontal MEC and EDDYFLX below 100 m, but involves strong contributions from the vertical turbulent fluxes within the upper 100 m, with a prominent downward turbulent flux across the base of the mixed layer. Our results thus highlight the leading order contribution of vertical turbulent fluxes in eddy-mean flow kinetic energy interactions at the base of the mixed layer.

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Vertical profile of horizontally integrated MEC (left), EDDYFLX (center) and DIVEF (right) within the green box of Figure 6. Three-dimensional estimates (black) are decomposed into an horizontal (blue) and vertical (red) contribution. Positive vertical eddy fluxes are oriented upward, and the dashed gray line represent the spatially averaged mixed layer depth at about 30 m.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Spatial correlation of MEC and EDDYFLX as a function of the coarse grained grid size at day 60. Each colored line is associated with a different 729 × 729 (i.e., 3⁶ × 3⁶) grid points box covering a slightly different portion of the full, 883 × 803 grid points domain. The lower left insert indicate the location of each boxes. The black line provides an averaged estimate of the correlation coefficient as a function of the coarse grained grid size.</figcaption></figure>

**4.2. Horizontal Scale Dependence**

Finally, we assess the scale-dependence of non-local kinetic energy transfers. Although at small scales, our results suggest eddy-mean flow interactions are largely non-local, our estimates on larger scales tend toward a local balance (i.e., DIVEF is negligible). It is true for the 3° × 3° green box of Figure 6, as well as for other places in the western Mediterranean basin (not shown), suggesting non-local effects are predominantly small scale features. We have thus computed the spatial correlation *r* between MEC and EDDYFLX as a function of coarse grained grid size (Figure 10). Starting from the initial model grid size at 1/60°, a spatial averaging is performed with the adjacent grid points, that is, a factor 3, up to a grid size of about 4°. This procedure has been performed on four different boxes of 3⁶ × 3⁶ (i.e., 729 × 729) grid points (colored lines) in order to cover the entire 883 × 803 grid points MEDWEST60 domain. The spatial correlation between MEC and EDDYFLX ranges from -0.12 on average at the model grid size to -0.96 at about 4°. This suggests that although non-local at small scales, kinetic energy transfers can be seen as local processes for scales larger than a few hundreds of kilometers. However, correlations lower than -0.5 are found for grid size of about 1/2° and finer, suggesting non-local dynamics would become leading order contribution as soon as mesoscale eddies are (even partially) resolved. It suggests that the processes associated with this non-locality need to be accounted for in the development of submesoscale parametrizations for eddy-permitting/eddy-resolving ocean models.

**5. Conclusion**

In this study, we have investigated the spatiotemporal structure of the kinetic energy transfers between the ensemble mean and the turbulent flow. We have performed our analysis with a kilometric-scale resolution (1/60°), 120-day long, 20-member ensemble simulations of the Western Mediterranean basin (Leroux et al., 2021). We have first introduced the Forced and Internal Kinetic Energy equation (FKE and IKE, respectively) in this framework, and discussed the implications for their interpretation. In particular, the prescribed surface and boundary forcings drive the basin integrated time rate of change of FKE, and the basin integrated time rate of change of IKE reflects the energy of the turbulent flow that develops within the domain through the non-linear dynamics sensitive to initial conditions. We have then quantified the respective contributions of Mean-to-Eddy energy Conversion (MEC, $\langle u_h \rangle \cdot \nabla \cdot \langle u' \otimes u' \rangle$) and the EDDY momentum FLuX (EDDYFLX,
$\langle \mathbf{u}' \otimes \mathbf{u}'_h \rangle \cdot \nabla \langle \mathbf{u}_h \rangle$ in the FKE and IKE budgets during the 120-day long runs. By further analyzing their spatial organization, we have then highlighted the non-locality of the energy transfers between the ensemble mean and the turbulent flow, where non-local processes are associated with energy destruction in one reservoir that does not locally sustain the growth of kinetic energy in the other reservoir, in agreement with previous studies (Capó et al., 2019; Chen et al., 2014; Kang & Curchitser, 2015). We have pointed out that non-local transfers are driven by turbulent fluxes of eddy-mean cross energy term, which are captured by the DIVergence of Eddy Flux (DIVEF, $\nabla \cdot \langle \mathbf{u}' \cdot (\mathbf{u}_h \cdot \mathbf{u}'_h) \rangle$). Our main contribution is to recognize that DIVEF is associated with advection of the cross energy term $\langle \mathbf{u}_h \rangle \cdot \mathbf{u}'_h$ by the turbulent flow, which provides a strong spatial constraint on these transfers since the cross energy term vanishes identically for turbulent flow orthogonal to the mean flow. Finally, we have shown that although weaker than the horizontal component at the model grid size, the vertical eddy fluxes become leading order when horizontally integrated over sufficiently large scales. On average, the three kind of vertical eddy energy fluxes (mean, eddy and cross energy term) are oriented downward at the base of the mixed layer.

Analyzing the scale dependence of these non-local KE transfers, we have shown that, although prevalent at eddy scales, KE transfers tend toward a local balance at non-eddying scale (i.e., >1°). Thus, while our results support approximations usually made in the development of energy-aware parameterizations of meso-scale turbulence (Eden & Greatbatch, 2008; Jansen et al., 2019; Mak et al., 2018), that is, that the growth of sub-grid scale turbulent kinetic energy is locally sustained by a weakening of the kinetic energy of the resolved flow, they point out to the necessity of accounting for non-local dynamics for the development of submesoscale parametrizations. In particular, accounting for such dynamics in eddy-permitting ocean models, such as those that will equibe the next generation climate model, could lead to significant improvements given non-locality has been found to be leading order contribution for scales as large as $\frac{1}{2}$°.

In this direction, the emerging approach of transport under *Location Uncertainty* (LU) for the representation of small scale, stochastic dynamics and its effect on the large scale flow (e.g., Chapron et al., 2018; Mémin, 2014; Resseguier et al., 2017) is an attractive alternative to the mixing length approach. Through a stochastic representation of the transport operator, LU indeed has the potential of providing interesting non-local properties, which will be the focus of future work.

We have performed our analysis based on ensemble simulations, with a view of inferring dynamical processes that need to be accounted for in submesoscale parametrizations. The ensemble approach differs from other time averaging, coarse-graining or spatial filtering methods. Although a comparative analysis between the different approaches is out of the scope of this paper, we want to point out to two potential benefits of ensemble simulations. First, when considering turbulence as the residual from a time averaging, ergodicity of the system is implied, that is, the time averaging is treated as an ensemble averaging. Although such assumption might be valid in the case of steady forcing, its validity is questionable for non-stationary systems. Thus, ensemble simulations may help in examining the response of eddy-mean interactions to changes in the forcing, such as what Uchida et al. (2022) have found for the seasonal variation of Eliassen-Palm fluxes in $\frac{1}{12}$°, 48-ensemble member ensemble simulations of the North Atlantic subtropical gyre. Second, coarse-graining or spatial filtering (e.g., Aluie et al., 2018; Grooms et al., 2021) approaches are subject to the definition of a length scale cut-off, thus to the size of the "eddies." However, it remains unclear how non-local energy transfers would depend on the length scale cut-off. In particular, questions remain on the spectral expression of MEC, EDDYFLX and DIVEF, as well as their respective contributions in fluxing energy up or down scale. We are currently investigating this last point and will report on the results in a dedicated paper.

Finally, we want to discuss the implications of our results for the interpretation of the dynamics of western boundary currents jet extension such as the Gulf Stream. Jamet et al. (2021) have recently shown the leading order contribution of MEC for the energetic balance of the North Atlantic subtropical, wind driven gyre. They concluded that MEC in the Gulf Stream extension region are the primary sink of 26-year mean kinetic energy within the gyre, balancing the energy inputted by the wind in the westerly wind region of the North Atlantic subtropical gyre. However, how this loss of mean kinetic energy interacts with the turbulent flow remains an open question. Some indications of spatial organization of EDDYFLX can be found in previous in-situ and satellite observation analyzes. In their earlier work on Gulf Stream energetics based on in-situ observations, Webster (1961), Webster (1965), Rossby (1987), and Dewar and Bane (1989) have reported on eddy fluxes that are more pronounced on the inshore flank of the Gulf Stream, both along the US coastline and downstream of
Cap Hatteras. Based on satellite observations, Ducet and Le Traon (2001) and Greatbatch et al. (2010) have highlighted a prominent feature of the Gulf Stream, so-called the “double-blade” structure, associated with the turbulent dynamics just downstream of Cape Hatteras. There, the Reynolds stress cross-covariance was found to be maximum on both flanks on the stream, and to exhibit alternation of highs and lows further downstream. This “double-blade” structure suggests that eddy fluxes (EDDYFLX) are more pronounced on the flank of the jet, where large Reynolds stresses $\overline{u'v'}$ are colocalized with a strong horizontal shear of the mean flow $\partial_y \overline{u}$, while mean-to-eddy conversion rates (MEC) would be more pronounced toward the core of the jet, where the cross-stream gradient of Reynolds stresses $\partial_y \overline{u'v'}$ are colocalized with maximum of the mean zonal current $\overline{u}$. We can also find some indications of such a spatial organization of eddy-mean flow interactions in the Lorenz energy cycle based on eddy-resolving numerical simulations of Kang and Curchitser (2015), although further analyses are needed to conclude on this.

## Appendix A: Offline Recomputation of Kinetic Energy Budget

We are interested in analyzing the energetic of the MEDWEST60 ensemble simulations, which have been recently produced (Leroux et al., 2021). We thus developed diagnostic tools to recompute the momentum budget, which kinetic energy builds upon, of these simulations based on the variables saved during the production of these simulations, that is, three-dimensional temperature (T), salinity (S) and velocity (U, V, W), as well as two-dimensional free-surface elevation (SSH). These offline tools are developed as part of the CDFTOOLS diagnostic package for the analysis of NEMO model output (https://github.com/meom-group/CDFTOOLS), which are written in FORTRAN 90 and follow the numerical implementation of the NEMO General Circulation Model (Madec et al., 2017).

As all GCM, NEMO offers different numerical schemes to integrate the Primitive Equations with various levels of approximation. The numerical schemes that have currently been implemented in these tools are those relevant for the analysis of the energetic of the MEDWEST60 ensemble simulations, which are based on the version 3.6 of the NEMO model. This includes: A dynamical vertical coordinate following the free surface elevation, with partial stepping along the ocean floor; the third order upstream biased scheme (UBS, Shchepetkin & McWilliams, 2005) to advect momentum; the TEOS-10 equation of state (Roquet et al., 2015) to compute density; a split-explicit formulation to compute surface pressure gradients (Shchepetkin & McWilliams, 2005), which also accounts for atmospheric surface pressure loading and freshwater air-land-sea fluxes; and an implicit time differencing scheme to compute vertical viscous effects, which include surface wind stress forcing following the CORE bulk flux formulation (Large & Yeager, 2004), bottom friction due to bottom boundary condition, tides, internal waves breaking and other short time scale currents, as well as vertical dissipation of momentum within the water column based on the Turbulent Kinetic Energy (TKE) turbulent closure scheme (Blanke & Delecluse, 1993; Gaspar et al., 1990; Mellor & Yamada, 1982). A full description of these schemes is available online (https://github.com/quentinjamet/CDFTOOLS/tree/cdf_medwest/note_KE_bgt_cdftools.pdf). With shorthands, the full kinetic energy budget can be represented as:

$$NXT = ADV + (HPG + SPG_{1st\ guess}) + SPG_{correction} + ZDF, \quad (A1)$$

where NXT refers to the time rate of change $\partial_t$ (before application of the Asselin filter), ADV to three-dimensional advection, HPG to hydrostatic pressure work, SPG$_{1st\ guess}$ to surface pressure work computed at baroclinic time step due to the rescaled vertical coordinate following free surface elevation, SPG$_{correction}$ to surface pressure work correction associated with the time-splitting scheme of Shchepetkin and McWilliams (2005) which includes atmospheric pressure loading and freshwater fluxes, and ZDF to vertical viscous effects.

## A1. Validation at Model Time Step

In order to insure that our offline recomputation lines up with the online estimates computed by the NEMO model, we have re-run for a short period of time one member of the ensemble and outputted, at the model time step ($\Delta t = 80$ s), momentum and kinetic energy trends, as well as required prognostic variables necessary for their offline recomputation, within the 150 × 150 grid point sub-region (black box on Figure 2). Comparing our offline recomputation with the online estimates provides an robust estimate of the errors. An example is provided on Figure A1 for the three-dimensional advection of kinetic energy within the model upper layer. The errors are relatively small
<figure><img src="image_13.png" /><figcaption>**Figure A1.** Upper layer Kinetic Energy trends associated with three-dimensional advection based on the model outputs (left), its offline recomputation (center), and associated errors (right). The offline recomputation is performed at model time step and accounts for the forward time discretization of the third order upstream biased part of UBS advective scheme. Note the different scale factor used for errors.</figcaption></figure>

(locally four order of magnitude, but five order of magnitude when horizontally averaged within the sub-domain, cf., Table A1), providing strong confidence in the accuracy of these tools. Tests for the other terms of the KE budget have been conducted, providing similar level of accuracy for time rate of change and pressure work (cf., Table A1). Offline estimates of vertical viscous effects are associated with much larger errors, of the order of 10%, and we currently have no estimates for the surface pressure correction associated with the split-explicit scheme.

## A2. Estimation of Errors Due to Time Discretization and Averaging

Based on model time step accuracy estimates, we have quantified the errors associated with time discretization of the different operators, as well as the use of time averaged quantities. We discuss here these implications for the estimates of the advective component of the budget.

The advective operator used in the MEDWEST60 is an upstream biased third order scheme (UBS, Shchepetkin & McWilliams, 2005). This scheme has two component, a second order scheme and a third order biased scheme. While the former is centered in time, the latter is implemented forward in time, that is, it is evaluated with *before* velocities. While this numerical detail provides stability for a GCM, it is not required in the context of *offline* computations and introduces ambiguities about how this should be evaluated when working with time averaged quantities. We thus decided to evaluate the third order biased scheme of the advective operator as centered in time instead. This leads to a growth of the errors made in the recomputation by one order of magnitude (cf., Table A1). When computed based on hourly model outputs, as available from MEDWEST60, the error increases by another order of magnitude to reach $10^{-3}$. Also increased from model time step to hourly model outputs, the accuracy of these *offline* diagnostic tools remains high, providing reliable estimates of the advective operator of the model.

<table><caption>**Table A1**</caption><caption>Order of Magnitude of the Relative Errors of the Offline Estimates for the Different Terms of the Kinetic Energy Budget, Computed as the Spatial Root-Mean-Square Error Normalized by the Spatial Standard Deviation of the Reference, NEMO Outputs</caption><thead><tr><th></th><th>∂<sub>t</sub>K</th><th>= −</th><th>V · (uK)</th><th>−</th><th>u<sub>h</sub> · ∇<sub>h</sub>φ<sub>hyd</sub></th><th>−</th><th>u<sub>h</sub> · ∇<sub>h</sub>φ<sub>surf</sub></th><th>+</th><th>ρ<sub>0</sub>u<sub>h</sub> · D<sup>m</sup></th></tr></thead><tbody><tr><th>Model time step</th><td>10<sup>-3</sup></td><td></td><td>10<sup>-5</sup></td><td></td><td>10<sup>-5</sup></td><td></td><td>−</td><td></td><td>10<sup>-1</sup></td></tr><tr><th>Time discretization</th><td>−</td><td></td><td>10<sup>-4</sup></td><td></td><td>−</td><td></td><td>−</td><td></td><td>10<sup>-1</sup></td></tr><tr><th>Hourly average</th><td>10<sup>-2</sup></td><td></td><td>10<sup>-3</sup></td><td></td><td>10<sup>-3</sup></td><td></td><td>−</td><td></td><td>10<sup>-1</sup></td></tr></tbody></table>

*Note.* The third line stands for the sensitivity of the error associated with the forward time discretization of the third order upstream biased part of UBS advective scheme and in the TKE turbulent closure scheme. We currently have estimates for the surface pressure work correction associated with the split-explicit scheme (third term of the RHS), such that no values are reported on here.
<figure><img src="image_14.png" /><figcaption>**Figure A2.** Time integrated KE trends of the full sub-domain, volume integrated time rate of change (upper left), three-dimensional advection (upper right), pressure work (bottom left) and vertical dissipation (bottom right) based on hourly averaged model outputs (black lines), recomputation based on hourly averaged T, S, U, V, W, η (red dots), and the associated errors (blue lines). Note the scale factor used for errors in the legend panels, which differs for each quantities.</figcaption></figure>

Similar considerations are applied for the vertical viscous effects (i.e., time discretization, hourly model outputs), but the already large error of $10^{-1}$ is found to be unchanged.

Finally, we estimate the evolution in time of these errors by comparing the recomputation made with hourly model outputs with estimates outputted by the model over a time period of 10 days (Figure A2). From these tests, no systematic errors emerged for both time rate of change (upper left panel) and hydrostatic pressure work (bottom left panel). We observe, however, a steady growth in the error made in the recomputation of the advective term (top right panel), reaching about $-20 \times 10^{-3}$ GW h$^{-1}$ at the end of the 10 days of simulation. Finally, the largest errors are observed in the recomputation of the vertical viscous effects (bottom right panel), in agreement with errors reported earlier. We are currently working on improving this recomputation.

### A3. Eddy-Mean Separation

Based on these offline estimates, we explicitly decompose the full equation into mean and eddy contributions. For the zonal momentum advection, it leads to:

$$ \nabla \cdot (\mathbf{u}\mathbf{u}) = \nabla \cdot (\langle\langle \mathbf{u} \rangle\rangle \langle\rangle u) + \nabla \cdot (\langle\langle \mathbf{u} \rangle\rangle u') + \nabla \cdot (\langle\langle \mathbf{u}' \rangle\rangle u) + \nabla \cdot (\langle\langle \mathbf{u}' \rangle\rangle u') \quad (A2) $$
where $\langle \cdot \rangle$ and $\cdot'$ denotes averaging and perturbation, respectively (cf., Section 2.1 for details on the decomposition used in this study). Performing a similar procedure for the advection of meridional momentum, multiplying the former by $\rho_0 (\langle u \rangle + u')$ and the latter by $\rho_0 (\langle v \rangle + v')$ and summing the resulting equations leads to a decomposition of the advection of kinetic energy that accounts for the different contributions that compose the FKE and IKE budgets (Equation 8 and 9, respectively). We note here that in MEDWEST60, the advection of momentum is achieved by the upstream biased third order scheme (UBS, Shchepetkin & McWilliams, 2005). This scheme accounts for the horizontal dissipation of momentum through an implicit formulation which takes the form of a biharmonic operator with an eddy coefficient proportional to the velocity $A_h = -|u|\Delta x^3/12$. The formulation of this implicit dissipation introduces complexities in the eddy-mean decomposition. We thus decided to evaluate the horizontal advection terms using a fourth order finite differencing centered scheme instead, which is the non-dissipative equivalent of the UBS scheme (Jouanno et al., 2016; Madec et al., 2017).

## Data Availability Statement

Further details on the MEDWEST60 ensemble simulations are available at https://zenodo.org/record/4570159 and the NEMO code used for the MEDWEST60 configuration is available at https://github.com/ocean-next/MEDWEST60/tree/main/src_config. Python scripts used to produce the figures of this manuscript are available at https://github.com/quentinjamet/publications-codes/tree/master/Jamet_etal_JAMES2022. Dedicated CDFTOOLS are available at https://github.com/quentinjamet/CDFTOOLS/tree/cdf_medwest.

## Acknowledgments

This work was supported through the French "Make Our Planet Great Again" program managed by the Agence Nationale de la Recherche under the Programme d'Investissement d'Avenir, with the reference ANR-18-MPGA-0002, as well as the NSF grants OCE-1829856 and OCE-2123632. MEDWEST60 ensemble simulations and their analyses have been performed using HPC resources from GENCI-IDRIS through the Grant A008-0111279 allocated in the context of the H2020-IMMERSE projet (European Union Horizon, 2020 research and innovation programme, Grant No. 821926). We want to thank Julien Jouanno and Robin Waldman for their help in the development of offline energy budgets in NEMO, as well as the two anonymous reviewers for their insightful comments and suggestions for the preparation of this manuscript.

## References

Aluie, H., Hecht, M., & Vallis, G. K. (2018). Mapping the energy cascade in the North Atlantic Ocean: The coarse-graining approach. *Journal of Physical Oceanography*, 48(2), 225–244. https://doi.org/10.1175/jpo-d-17-0100.1

Blanke, B., & Delecluse, P. (1993). Variability of the tropical Atlantic Ocean simulated by a general circulation model with two different mixed-layer physics. *Journal of Physical Oceanography*, 23(7), 1363–1388. https://doi.org/10.1175/1520-0485(1993)023<1363:vottao>2.0.co;2

Brankart, J.-M., Candille, G., Garnier, F., Calone, C., Melet, A., Bouttier, P.-A., et al. (2015). A generic approach to explicit simulation of uncertainty in the nemo ocean model. *Geoscientific Model Development*, 8(5), 1285–1297. https://doi.org/10.5194/gmd-8-1285-2015

Brodeau, L., Le Sommer, J., & Albert, A. (2020). Ocean-next/eNATL60: Material describing the set-up and the assessment of NEMO-eNATL60 simulations (Version v1) [Dataset]. Zenodo. https://doi.org/10.5281/zenodo.4032732

Capó, E., Orfila, A., Mason, E., & Ruiz, S. (2019). Energy conversion routes in the Western Mediterranean Sea estimated from eddy–mean flow interactions. *Journal of Physical Oceanography*, 49(1), 247–267. https://doi.org/10.1175/jpo-d-18-0036.1

Cessi, P. (2008). An energy-constrained parameterization of eddy buoyancy flux. *Journal of Physical Oceanography*, 38(8), 1807–1819. https://doi.org/10.1175/2007jpo3812.1

Chapron, B., Dérrian, P., Mémin, E., & Resseguier, V. (2018). Large-scale flows under location uncertainty: A consistent stochastic framework. *Quarterly Journal of the Royal Meteorological Society*, 144(710), 251–260. https://doi.org/10.1002/qj.3198

Chen, R., Flierl, G. R., & Wunsch, C. (2014). A description of local and nonlocal eddy–mean flow interaction in a global eddy-permitting state estimate. *Journal of Physical Oceanography*, 44(9), 2336–2352. https://doi.org/10.1175/jpo-d-14-0009.1

Chen, R., Thompson, A. F., & Flierl, G. R. (2016). Time-dependent eddy-mean energy diagrams and their application to the ocean. *Journal of Physical Oceanography*, 46(9), 2827–2850. https://doi.org/10.1175/jpo-d-16-0012.1

Davies, H. (1976). A lateral boundary formulation for multi-level prediction models. *Quarterly Journal of the Royal Meteorological Society*, 102(432), 405–418. https://doi.org/10.1002/qj.49710243210

Dewar, W. K., & Bane, J. M. (1989). Gulf Stream dynamics. Pad II: Eddy energetics at 73°W. *Journal of Physical Oceanography*, 19(10), 1574–1587. https://doi.org/10.1175/1520-0485(1989)019<1574:gsdpie>2.0.co;2

Ducet, N., & Le Traon, P.-Y. (2001). A comparison of surface eddy kinetic energy and Reynolds stresses in the Gulf Stream and the Kuroshio Current systems from merged TOPEX/Poseidon and ERS-1/2 altimetric data. *Journal of Geophysical Research*, 106(C8), 16603–16622. https://doi.org/10.1029/2000jc000205

Eden, C., & Greatbatch, R. J. (2008). Towards a mesoscale eddy closure. *Ocean Model*, 20(3), 223–239. https://doi.org/10.1016/j.ocemod.2007.09.002

Engedahl, H. (1995). Use of the flow relaxation scheme in a three-dimensional baroclinic ocean model with realistic topography. *Tellus A*, 47(3), 365–382. https://doi.org/10.1034/j.1600-0870.1995.t01-2-00006.x

Ferreira, D., Marshall, J., & Heimbach, P. (2005). Estimating eddy stresses by fitting dynamics to observations using a residual-mean ocean circulation model and its adjoint. *Journal of Physical Oceanography*, 35(10), 1891–1910. https://doi.org/10.1175/jpo2785.1

Flather, R. A. (1994). A storm surge prediction model for the northern Bay of Bengal with application to the cyclone disaster in April 1991. *Journal of Physical Oceanography*, 24(1), 172–190. https://doi.org/10.1175/1520-0485(1994)024<0172:asspmf>2.0.co;2

Fox-Kemper, B., Ferrari, R., & Hallberg, R. (2008). Parameterization of mixed layer eddies. Part I: Theory and diagnosis. *Journal of Physical Oceanography*, 38(6), 1145–1165. https://doi.org/10.1175/2007jpo3792.1

Gaspar, P., Grégoris, Y., & Lefevre, J.-M. (1990). A simple eddy kinetic energy model for simulations of the oceanic vertical mixing: Tests at station papa and long-term upper ocean study site. *Journal of Geophysical Research*, 95(C9), 16179–16193. https://doi.org/10.1029/jc095ic09p16179

Gent, P. R., & McWilliams, J. C. (1990). Isopycnal mixing in ocean circulation models. *Journal of Physical Oceanography*, 20(1), 150–155. https://doi.org/10.1175/1520-0485(1990)020<0150:imiocm>2.0.co;2

Greatbatch, R. J., Zhai, X., Kohlmann, J.-D., & Czeschel, L. (2010). Ocean eddy momentum fluxes at the latitudes of the Gulf Stream and the Kuroshio extensions as revealed by satellite data. *Ocean Model*, 60(3), 617–628. https://doi.org/10.1007/s10236-010-0282-6

Groeskamp, S., LaCasce, J. H., McDougall, T. J., & Rogé, M. (2020). Full-depth global estimates of ocean mesoscale eddy mixing from observations and theory. *Geophysical Research Letters*, 47(18), e2020GL089425. https://doi.org/10.1029/2020GL089425
Grooms, I. (2017). Simulations of eddy kinetic energy transport in barotropic turbulence. *Physical Review Fluids*, 2(11), 113801. https://doi.org/10.1103/physrevfluids.2.113801

Grooms, I., Loose, N., Abernathey, R., Steinberg, J., Bachman, S. D., Marques, G., et al. (2021). Diffusion-based smoothers for spatial filtering of gridded geophysical data. *Journal of Advances in Modeling Earth Systems*, 13(9), e2021MS002552. https://doi.org/10.1029/2021ms002552

Grooms, I., Nadeau, L.-P., & Smith, K. S. (2013). Mesoscale eddy energy locality in an idealized ocean model. *Journal of Physical Oceanography*, 43(9), 1911–1923. https://doi.org/10.1175/jpo-d-13-036.1

Hawkins, E., Smith, R. S., Gregory, J. M., & Stainforth, D. A. (2016). Irreducible uncertainty in near-term climate projections. *Climate Dynamics*, 46(11–12), 3807–3819. https://doi.org/10.1007/s00382-015-2806-8

Jamet, Q., Deremble, B., Wienders, N., Uchida, T., & Dewar, W. (2021). On wind-driven energetics of subtropical gyres. *Journal of Advances in Modeling Earth Systems*, 13(4), e2020MS002329. https://doi.org/10.1029/2020ms002329

Jansen, M. F., Adcroft, A., Khani, S., & Kong, H. (2019). Toward an energetically consistent, resolution aware parameterization of ocean mesoscale eddies. *Journal of Advances in Modeling Earth Systems*, 11(8), 2844–2860. https://doi.org/10.1029/2019ms001750

Jouanno, J., Capet, X., Madec, G., Roullet, G., & Klein, P. (2016). Dissipation of the energy imparted by mid-latitude storms in the Southern Ocean. *Ocean Science*, 12(3), 743–769. https://doi.org/10.5194/os-12-743-2016

Kang, D., & Curchitser, E. N. (2015). Energetics of eddy–mean flow interactions in the Gulf Stream region. *Journal of Physical Oceanography*, 45(4), 1103–1120. https://doi.org/10.1175/jpo-d-14-0200.1

Large, W. G., & Yeager, S. G. (2004). Diurnal to decadal global forcing for ocean and sea-ice models: The data sets and flux climatologies.

Leroux, S., Brankart, J.-M., Albert, A., Brodeau, L., Molines, J.-M., Jamet, Q., et al. (2021). Ensemble quantification of short-term predictability of the ocean fine-scale dynamics: A western mediterranean test-case at kilometric-scale resolution. *Ocean Science Discussions*. https://doi.org/10.5194/os-2022-11

Madec, G., Bourdallé-Badie, R., Bouttier, P.-A., Bricaud, C., Bruciaferri, D., Calvert, D., et al. (2017). NEMO ocean engine.

Mak, J., Maddison, J. R., Marshall, D. P., & Munday, D. R. (2018). Implementation of a geometrically informed and energetically constrained mesoscale eddy parameterization in an ocean circulation model. *Journal of Physical Oceanography*, 48(10), 2363–2382. https://doi.org/10.1175/jpo-d-18-0017.1

Marshall, D. P., & Adcroft, A. J. (2010). Parameterization of ocean eddies: Potential vorticity mixing, energetics and Arnold’s first stability theorem. *Ocean Model*, 32(3–4), 188–204. https://doi.org/10.1016/j.ocemod.2010.02.001

Marshall, D. P., Maddison, J. R., & Berloff, P. S. (2012). A framework for parameterizing eddy potential vorticity fluxes. *Journal of Physical Oceanography*, 42(4), 539–557. https://doi.org/10.1175/jpo-d-11-048.1

Mellor, G. L., & Yamada, T. (1982). Development of a turbulence closure model for geophysical fluid problems. *Reviews of Geophysics*, 20(4), 851–875. https://doi.org/10.1029/rg020i004p00851

Mémin, E. (2014). Fluid flow dynamics under location uncertainty. *Geophysical & Astrophysical Fluid Dynamics*, 108(2), 119–146. https://doi.org/10.1080/03091929.2013.836190

Millot, C. (1985). Some features of the Algerian Current. *Journal of Geophysical Research*, 90(C4), 7169–7176. https://doi.org/10.1029/jc090ic04p07169

Millot, C. (1999). Circulation in the Western Mediterranean Sea. *Journal of Marine Systems*, 20(1–4), 423–442. https://doi.org/10.1016/s0924-7963(98)00078-5

Obaton, D., Millot, C., d'Hières, G. C., & Taupier-Letage, I. (2000). The Algerian Current: Comparisons between in situ and laboratory data sets. *Deep Sea Research Part I: Oceanographic Research Papers*, 47(11), 2159–2190. https://doi.org/10.1016/s0967-0637(00)00014-5

Poulain, P.-M., Centurioni, L., Özgökmen, T., Tarry, D., Pascual, A., Ruiz, S., et al. (2021). On the structure and kinematics of an Algerian Eddy in the southwestern Mediterranean Sea. *Remote Sensing*, 13(15), 3039. https://doi.org/10.3390/rs13153039

Renault, L., Molemaker, M. J., Gula, J., Masson, S., & McWilliams, J. C. (2016). Control and stabilization of the Gulf Stream by oceanic current interaction with the atmosphere. *Journal of Physical Oceanography*, 46(11), 3439–3453. https://doi.org/10.1175/jpo-d-16-0115.1

Renault, L., Molemaker, M. J., McWilliams, J. C., Shchepetkin, A. F., Lemarié, F., Chelton, D., et al. (2016). Modulation of wind work by oceanic current interaction with the atmosphere. *Journal of Physical Oceanography*, 46(6), 1685–1704. https://doi.org/10.1175/jpo-d-15-0232.1

Resseguier, V., Mémin, E., & Chapron, B. (2017). Geophysical flows under location uncertainty, Part I random transport and general models. *Geophysical & Astrophysical Fluid Dynamics*, 111(3), 149–176. https://doi.org/10.1080/03091929.2017.1310210

Roquet, F., Madec, G., McDougall, T. J., & Barker, P. M. (2015). Accurate polynomial expressions for the density and specific volume of seawater using the TEOS-10 standard. *Ocean Modelling*, 90, 29–43. https://doi.org/10.1016/j.ocemod.2015.04.002

Rossby, T. (1987). On the energetics of the Gulf Stream at 73°W. *Journal of Marine Research*, 45(1), 59–82. https://doi.org/10.1357/002224087788400918

Schubert, R., Gula, J., Greatbatch, R. J., Baschek, B., & Biastoch, A. (2020). The submesoscale kinetic energy cascade: Mesoscale absorption of submesoscale mixed layer eddies and frontal downscale fluxes. *Journal of Physical Oceanography*, 50(9), 2573–2589. https://doi.org/10.1175/jpo-d-19-0311.1

Shchepetkin, A. F., & McWilliams, J. C. (2005). The regional oceanic modeling system (ROMS): A split-explicit, free-surface, topography- following-coordinate oceanic model. *Ocean modelling*, 9(4), 347–404. https://doi.org/10.1016/j.ocemod.2004.08.002

Stainforth, D. A., Allen, M. R., Tredger, E. R., & Smith, L. A. (2007). Confidence, uncertainty and decision-support relevance in climate predic- tions. *Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences*, 365(1857), 2145–2161. https:// doi.org/10.1098/rsta.2007.2074

Uchida, T., Jamet, Q., Dewar, W. K., Le Sommer, J., Penduff, T., & Balwada, D. (2022). Diagnosing the thickness-weighted averaged eddy–mean flow interaction from an eddying North Atlantic ensemble: The Eliassen-Palm flux. *Journal of Advances in Modeling Earth Systems*, 14(5), e2021MS002866. https://doi.org/10.1029/2021ms002866

Visbeck, M., Marshall, J., Haine, T., & Spall, M. (1997). Specification of eddy transfer coefficients in coarse-resolution ocean circulation models. *Journal of physical oceanography*, 27(3), 381–402. https://doi.org/10.1175/1520-0485(1997)027<0381:soetci>2.0.co;2

Waldman, R. (2016). *Etude multi-échelle de la convection océanique profonde en mer méditerranée: De l'observation à la modélisation climatique* (Unpublished doctoral dissertation). Université Paul Sabatier-Toulouse III.

Webster, F. (1961). The effect of meanders on the kinetic energy balance of the Gulf Stream. *Tellus*, 13(3), 392–401. https://doi.org/10.1111/j.2153- 3490.1961.tb00100.x

Webster, F. (1965). Measurements of eddy fluxes of momentum in the surface layer of the Gulf Stream. *Tellus*, 17(2), 240–245. https://doi.org/10.1111/j.2153-3490.1965.tb01415.x

Wunsch, C., & Ferrari, R. (2004). Vertical mixing, energy, and the general circulation of the oceans. *Annual Review of Fluid Mechanics*, 36(1), 281–314. https://doi.org/10.1146/annurev.fluid.36.050802.122121