# Diagnosing the Thickness-Weighted Averaged Eddy-Mean Flow Interaction From an Eddying North Atlantic Ensemble: The Eliassen-Palm Flux

Takaya Uchida, Quentin Jamet, William Dewar, Julien Le Sommer, Thierry Penduff, Dhruv Balwada

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

## ► To cite this version:

Takaya Uchida, Quentin Jamet, William Dewar, Julien Le Sommer, Thierry Penduff, et al.. 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, 2022, 14 (5), 10.1029/2021MS002866. [https://doi.org/10.1029/2021MS002866](https://doi.org/10.1029/2021MS002866)

HAL Id: hal-03669662

https://hal.science/hal-03669662

Submitted on 17 May 2022

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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.

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

10.1029/2021MS002866

**Key Points:**

* Eddying ensemble runs of the North Atlantic Ocean are used to diagnose the thickness-weighted averaged eddy-mean flow interaction

* A dynamically-consistent approximately neutral surface is implemented to define the buoyancy coordinate for a realistic equation of state

* The Eliassen-Palm flux convergence implies a tendency to force a poleward migration of the Gulf Stream

**Correspondence to:**

T. Uchida,
takaya.uchida@univ-grenoble-alpes.fr

**Citation:**

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*, e2021MS002866.
https://doi.org/10.1029/2021MS002866

Received 12 OCT 2021
Accepted 20 MAR 2022

© 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 License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

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

# Diagnosing the Thickness-Weighted Averaged Eddy-Mean Flow Interaction From an Eddying North Atlantic Ensemble: The Eliassen-Palm Flux

Takaya Uchida<fnref n="1" /> Ⓓ, Quentin Jamet<fnref n="1" /> Ⓓ, William K. Dewar<fnref n="1" />,<fnref n="2" />, Julien Le Sommer<fnref n="1" /> Ⓓ, Thierry Penduff<fnref n="1" /> Ⓓ, and Dhruv Balwada³ Ⓓ

**Abstract** The thickness-weighted average (TWA) framework, which treats the residual-mean flow as the prognostic variable, provides a clear theoretical formulation of the eddy feedback onto the residual-mean flow. The averaging operator involved in the TWA framework, although in theory being an ensemble mean, in practice has often been approximated by a temporal mean. Here, we analyze an ensemble of North Atlantic simulations at mesoscale-permitting resolution (1/12°). We therefore recognize means and eddies in terms of ensemble means and fluctuations about those means. The ensemble dimension being orthogonal to the temporal and spatial dimensions negates the necessity for an arbitrary temporal or spatial scale in defining the eddies. Eddy-mean flow feedbacks are encapsulated in the Eliassen-Palm (E-P) flux tensor and its convergence indicates that eddy momentum fluxes dominate in the separated Gulf Stream. The eddies can be interpreted to contribute to the zonal meandering of the Gulf Stream and a northward migration of it in the meridional direction. Downstream of the separated Gulf Stream in the North Atlantic Current region, the interfacial form stress convergence becomes leading order in the E-P flux convergence.

**Plain Language Summary** We have greatly benefited from global climate simulations in gaining insight into what the climate would look like in an ever warming future. Due to computational constraints, however, the oceanic component of such simulations have been poorly constrained. The storm systems of the ocean, often referred to as eddies, defined as fluctuations about jets such as the Gulf Stream and meandering of the jet itself, have remained challenging to accurately simulate on a global scale. Although relatively small in scale compared to the global ocean, eddies have been known to modulate the climate by transporting heat from the equator to the poles. By running a regional simulation of the North Atlantic Ocean and taking advantage of recent theoretical developments, we implement a new framework to evaluate such simulations in representing the Gulf Stream.

## 1. Introduction

Eddy-mean flow interaction has been a key framework in understanding jet formation in geophysical flows such as in the atmosphere and ocean (Bühler, 2014; Vallis, 2017). A prominent example of such a jet in the North Atlantic ocean is the Gulf Stream. Previous studies have shown how eddies fluxing buoyancy and momentum back into the mean flow energize the western boundary currents including the Gulf Stream (Aluie et al., 2018; Chassignet & Xu, 2017; Lévy et al., 2010; Waterman & Lilly, 2015). Basin-scale simulations, however, often lack sufficient spatial resolution to accurately resolve the eddies and hence, result in underestimating the eddy fluxes of momentum and tracers (Arbic et al., 2013; Balwada et al., 2018; Capet et al., 2008b; Kjellsson & Zanna, 2017; Schubert et al., 2020; Uchida et al., 2019). Due to computational constraints, we will continue to rely on models which only partially resolve the mesoscale, a scale roughly on the order of O (20–200 km) at which the ocean currents are most energetic (Ajayi et al., 2020; Stammer, 1997; Xu & Fu, 2011, 2012), for global ocean and climate simulations. As a result, there has been an on-going effort to develop energy-backscattering eddy parametrizations which incorporate the dynamical effects of eddy momentum fluxes due to otherwise unresolved mesoscale turbulence (e.g., Bachman, 2019; Bachman et al., 2018; Berloff, 2018; Guillaumin & Zanna, 2021; Jansen et al., 2019; Juricke et al., 2020; Kitsios et al., 2013; Perezhogin, 2019; Uchida et al., 2022; Zanna & Bolton, 2020; Zanna et al., 2017).

<footnote n="1">¹Université Grenoble Alpes, CNRS, IRD, Grenoble-INP, Institut des Géosciences de l'Environnement, Grenoble, France,</footnote>

<footnote n="2">²Department of Earth, Ocean and Atmospheric Science, Florida State University, Tallahassee, FL, USA, ³Lamont-Doherty</footnote>

<footnote>Earth Observatory, Columbia University, New York, NY, USA</footnote>
There has been less emphasis, however, on quantifying the spatial and temporal characteristics of the eddy buoyancy and momentum fluxes themselves, which the parametrizations are deemed to represent. The focus of this study is, therefore, to examine the dynamical effects of mesoscale turbulence on the mean flow in realistic, partially air-sea coupled, eddying ensemble runs of the North Atlantic. The thickness-weighted average (TWA) framework, which treats the residual-mean velocity as a prognostic variable, allows for a straightforward theoretical expression of the eddy feedback onto the residual-mean flow (e.g., Andrews, 1983; Aoki, 2014; de Szoeke & Bennett, 1993; Gallimore & Johnson, 1981; Maddison & Marshall, 2013; McDougall & McIntosh, 2001; Young, 2012). It is well known in the atmospheric and Southern Ocean literature that it is the residual-mean flow, which is the residual that emerges upon the partial cancellation between the Eulerian mean flow and eddies, that captures the “mean” flow for heat and tracer transport (Bühler, 2014; Vallis, 2017). The TWA framework has been fruitful in examining eddy-mean flow interaction in idealized modeling studies (e.g., Bire & Wolfe, 2018; Cessi & Wolfe, 2013; D. P. Marshall et al., 2012; Ringler et al., 2017). Here, we extend these studies to a realistic simulation of the North Atlantic. We will examine the TWA eddy diffusivities and mode water formation in subsequent papers.

To our knowledge, Aiki and Richards (2008), Aoki et al. (2016), Stanley (2018) and Zhao and Marshall (2020) are the only studies that diagnose the TWA framework in realistic ocean simulations. Aiki and Richards (2008), however, recompute the hydrostatic pressure using potential density for their off-line diagnosis in defining their buoyancy coordinate, which can result in significant discrepancies from the pressure field used in their on-line calculation and consequently errors in the diagnosed geostrophic shear. Although Aoki et al. (2016) negate this complication between the buoyancy coordinate and mean pressure field by analyzing their outputs in geopotential coordinates, they compute the eddy component of the pressure term ($F^+$ in their paper) using potential density, resulting in errors in the interfacial form stress (viz. this violates Equation 10 described below for $\phi'$ and $m'$). Their truncation in Taylor expansion about the mean position of buoyancy surfaces for the sake of convenience in diagnosing the residual-mean flow in geopotential coordinates limits the accuracy of the eddy terms. Lastly, all four studies assume ergodicity. The ergodic assumption of treating a temporal mean equivalent to an ensemble mean, although a pragmatic one and has its place for examining the climate where the time scales are of interest, prevents examining the temporal evolution of the residual-mean fields and conflates temporal variability with the eddies. The conflation can have leading-order consequences in quantifying the energy cycle; by adjusting the temporal mean from monthly to annual, Aiki and Richards (2008, cf. Table 2 in their paper) show that the amount of kinetic and potential energy stored in the mean and eddy reservoirs can change by up to a factor of four. Eddy-mean flow interaction in the TWA framework, hence, warrants further investigation, and we believe our study is the first to strictly implement an ensemble mean in this context. In Sections 4.1 and 4.2, we show that the ensemble framework provides new insights into turbulence studies.

When discussing eddy versus mean flow, one of the ambiguities lies in how the two are decomposed and interpreted (Bachman et al., 2015). As noted above, often, the eddies are defined from a practical standpoint as the deviation from a temporally and/or spatially coarse-grained field regardless of the coordinate system (e.g., Aiki & Richards, 2008; Aoki et al., 2016; Griffies et al., 2015; Lévy et al., 2012; Sasaki et al., 2014; Uchida et al., 2017; Zhao & Marshall, 2020), which leaves open the question of how the filtering affects the decomposition. Due to the ensemble averaging nature of the TWA framework, we are able to uniquely define the two; the mean flow (ensemble mean) is the oceanic response to the surface boundary state and lateral boundary conditions, and the eddy (fluctuations about the ensemble mean) is the field due to intrinsic variability including mesoscale turbulence (Leroux et al., 2018; Sérazin et al., 2017).

The paper is organized as follows: We describe the model configuration in Section 2 and briefly provide an overview of the TWA framework in Section 3. The results are given in Section 4. In particular, our data set provides a unique opportunity to examine the validity of the often assumed ergodicity when decomposing the flow into its eddy and mean flow components, which we give in Section 4.2. Discussion and conclusions are given in Section 5.
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Bathymetry of the modeled domain. The domain was configured to wrap around zonally in order to save computation and memory allocation when generating the ensemble. The hatches indicate the northern and southern regions excluded from our analysis.</figcaption></figure>

## 2. Model Description

We use the model outputs from the realistic runs described in Jamet et al. (2019b), Jamet et al. (2020) and Uchida, Jamet, et al. (2021), which are an air-sea partially coupled, 48-member ensemble of the North Atlantic ocean at mesoscale-permitting resolution (1/12°; or sometimes referred to as “eddy rich”) using the hydrostatic configuration of the Massachusetts Institute of Technology general circulation model (MITgcm; J. Marshall et al., 1997). We have 46 vertical levels increasing from 6 m near the surface to 250 m at depth. Harmonic, biharmonic horizontal and vertical viscosity values of $A_{h2} = 20 \text{ m}^2\text{s}^{-1}$, $A_{h4} = 10^{10} \text{ m}^4\text{s}^{-1}$ and $A_v = 10^{-5} \text{ m}^2\text{s}^{-1}$ were used respectively. For completeness, we provide a brief summary of the configuration below.

Figure 1 shows the bathymetry of the modeled domain extending from 20°S to 55°N. In order to save computational time and memory allocation, the North Atlantic basin was configured to zonally wrap around periodically. Open boundary conditions are applied at the north and south boundaries of our domain and Strait of Gibraltar, such that oceanic velocities (u) and potential temperature and practical salinity (Θ, S) are restored with a 36 min relaxation time scale toward a state derived by an ocean-only global Nucleus for European Modeling of the Ocean (NEMO) simulation (Molines et al., 2014, ORCA12. L46-MJM88 run in their paper, hereon referred to as ORCA12). The open boundary conditions are prescribed every 5 days from the ORCA12 run and linearly interpolated in between. A sponge layer is further applied to two adjacent grid points from the open boundaries where model variables are restored toward boundary conditions with a 1-day relaxation time scale. In total, relaxation is applied along three grid points from the boundaries with it being the strongest at the boundary along with radiation conditions at the northern/southern most boundary. Although relatively short, no adverse effects were apparent upon inspection in response to these relaxation time scales; for example, changes in the open boundary conditions were seen to induce a physically consistent Atlantic Meridional Overturning Circulation response inside the domain (Jamet et al., 2020).

The 48-member ensemble was constructed as follows: 48 oceanic states separated by 48 hr each were taken during an initial 96-day-long integration beginning 14 November 1962. Simulations initialized with these states were then run under yearly *repeating* 1963 atmospheric and boundary conditions for a year, that is, the atmospheric state and boundary conditions are cyclic for this year. After the 1 year of integration from the 48 states,
the last time step from each simulation was taken as the initial condition for the ensuing ensemble members; each spun-up initial oceanic state is physically consistent with the atmospheric and boundary conditions of 1 January 1963 (details are given in Jamet et al., 2020). At the surface, the ocean is partially coupled to an atmospheric boundary layer model (CheapAML; Deremble et al., 2013). In CheapAML, atmospheric surface temperature and relative humidity respond to ocean surface structures by exchanges of heat and humidity computed according to the Coupled Ocean-Atmosphere Response Experiment (COARE3; Fairall et al., 2003) flux formula, but are strongly restored toward prescribed values over land; there are no zonally propagating signals of climate teleconnection. The prescribed atmospheric state is taken from the Drakkar forcing set and boundary forcing from the ORCA12 run (details are given in Jamet et al., 2019a). The ensemble members are integrated forward in time for 5 years (1963–1967), and exposed to the same prescribed atmospheric state above the boundary layer and relaxation at the north/south boundaries across all ensemble members. (Note that the forcing and relaxation are no longer cyclic after the 1-year spin-up phase.) During this interval, the oceanic state and the atmospheric boundary layer temperature and humidity evolve in time. In the following, we interpret the ensemble mean as the ocean response to the atmospheric state prescribed above the atmospheric boundary layer as well as the oceanic conditions imposed at the open boundaries of the regional domain, while the ensemble spread is attributed to intrinsic ocean dynamics that develop at mesoscale-permitting resolution (Jamt et al., 2019b; Leroux et al., 2018; Sérazin et al., 2017).

The model outputs were saved as 5-day averages. In the context of mesoscale dynamics, which is the focus of this study, some temporal averaging is appropriate in order to filter out temporal scales shorter than the mesoscale eddies themselves. From a probabilistic perspective, the 5-day averaging results in more Gaussian-like eddy statistics (based on the central-limit theorem). From a dynamical point of view, this does not allow us to close the residual-mean and eddy budgets (cf. Stanley, 2018, Section 4.4). Nevertheless, the ensemble dimension of our data set provides an unique opportunity to examine the TWA eddy-mean flow interaction. In the following analysis, we exclude the northern and southern extent of 5° from our analysis to avoid effects from the open boundary conditions and sponge layer (Figure 1) and to maximize the signal of intrinsic variability amongst the ensemble members. We also use the last year of output (1967) for the same reasons.

### 3. Theory and Implementation of Thickness-Weighted Averaging

The ocean is a stratified fluid, and the circulation and advection of tracers tend to align themselves along the stratified density surfaces. Hence, a natural way to understand the circulation is to consider the variables in a buoyancy framework and the residual-mean flow rather than the Eulerian mean flow. We leave the detailed derivation of the TWA framework to Young (2012, and references therein) and here, only provide a brief summary; the primitive equations in geopotential coordinates are first transformed to buoyancy coordinates upon which a thickness weighting and ensemble averaging along constant buoyancy surfaces are applied to obtain the TWA governing equations. Following the notation by Young (2012) and Ringler et al. (2017), the TWA horizontal momentum equations in the buoyancy coordinate system ($\tilde{t}$, $\tilde{x}$, $\tilde{y}$, $\tilde{b}$) are:

$$ \hat{u}_{\tilde{t}} + \hat{u}\hat{u}_{\tilde{x}} + \hat{v}\hat{u}_{\tilde{y}} + \hat{\sigma}\hat{u}_{\tilde{b}} - f\hat{v} + \bar{m}_{\tilde{x}} = -\bar{e}_{1} \cdot (\tilde{\nabla} \cdot \mathbf{E}) + \hat{\chi} \quad (1) $$

$$ \hat{v}_{\tilde{t}} + \hat{u}\hat{v}_{\tilde{x}} + \hat{v}\hat{v}_{\tilde{y}} + \hat{\sigma}\hat{v}_{\tilde{b}} + f\hat{u} + \bar{m}_{\tilde{y}} = -\bar{e}_{2} \cdot (\tilde{\nabla} \cdot \mathbf{E}) + \hat{\mathcal{Y}} \quad (2) $$

where $\bar{\cdot}(\cdot)$ and $\hat{\cdot}(\cdot)$ are the ensemble averaged and TWA variables respectively, $\sigma(=\zeta_b)$ the specific thickness and $\zeta$ the depth of an iso-surface of buoyancy. The subscripts denote partial derivatives. The Montgomery potential is $m = \dot{\phi} - \tilde{b}\zeta$ where $\dot{\phi}$ is the dynamically active part of hydrostatic pressure (the meaning of $\tilde{\cdot}(\cdot)$ will become clearer later). $\sigma$ is the dia-surface velocity across buoyancy contours, which we detail below for a realistic equation of state (EOS) for density. The vectors $\bar{e}_1 = i + \bar{\zeta}_x k$ and $\bar{e}_2 = j + \bar{\zeta}_y k$ form the basis vectors spanning the buoyancy horizontal space where $i, j$ and $k$ are the Cartesian geopotential unit vectors, and $\mathbf{E}$ is the E-P flux tensor described in detail in Section 4.1. Although each ensemble member has an individual basis ($e_1, e_2$), the E-P flux divergence yields no cross terms upon averaging as the TWA operator commutes with the divergence of $\mathbf{E}$ (for mathematical details, see Section 3.4 in Maddison & Marshall, 2013); this allows for the tensor expression in Equations 1 and 2. $\mathcal{X}$ and $\mathcal{Y}$ are the viscous and forcing terms.
One subtle yet important point involves the buoyancy coordinate ($\tilde{b}$) for a realistic, non-linear EOS (Jackett & McDougall, 1995). The analysis in Young (2012) implicitly assumes a linear EOS. With a realistic EOS, defining the vertical coordinate using potential density introduces errors. However, what constitutes a better buoyancy variable is the subject of some debate (e.g., de Szoeke & Springer, 2009; Jackett & McDougall, 1997; Klocker et al., 2009; Lang et al., 2020; McDougall & Jackett, 2005; Tailleux, 2016). Although other choices are possible, we argue for the use of in-situ density anomaly ($\delta$ def = $\rho$ - $\tilde{\rho}(z)$ where $\rho$ is the in-situ density and $\tilde{\rho}$ is a function of only depth; Montgomery, 1937; Stanley, 2018, 2019). With in-situ density anomaly, buoyancy can be defined as:

$$ \tilde{b}(\Theta, S, z) \stackrel{\text{def}}{=} -\frac{g}{\rho_0} \delta \stackrel{\text{def}}{=} \tilde{b}(t, x, y, z) \quad (3) $$

where $\rho_0 = 999.8 \text{ kg m}^{-3}$ the Boussinesq reference density prescribed in MITgcm. $b$ is used to denote a thermodynamic function and $\tilde{b}$ denotes the buoyancy at a point in space-time. The question becomes how to choose $\tilde{\rho}(z)$ so that monotonicity is maintained ($[b]_{\Theta,S} > 0$; the vertical partial derivative is taken in respect to constant potential temperature and practical salinity ($[\cdot|_{\Theta,S}] \cdot b(\Theta, S, z)]_{\Theta,S} > 0$ implies $\frac{\partial}{\partial z}\tilde{b}(t, x, y, z) > 0$ if the stratification is statically stable). The vertical derivative of the in-situ density anomaly can be decomposed as:

$$ [\delta_z]_{\Theta,S} = [\rho_z]_{\Theta,S} - \frac{d}{dz} \tilde{\rho} = [\rho_{\Phi}]_{\Theta,S} \frac{d\Phi}{dz} - \frac{d}{dz} \tilde{\rho} = -\frac{\rho_0 g}{c_s^2} - \frac{d}{dz} \tilde{\rho}, \quad (4) $$

where $\Phi = -\rho_0 g z$ is the dynamically non-active part of hydrostatic pressure, and $c_s$ is the sound speed. We remind the reader that a Boussinesq fluid is not strictly incompressible and a finite sound speed can be diagnosed (Olbers et al., 2012; Vallis, 2017). For simplicity, we can write $\frac{d}{dz}\tilde{\rho} \stackrel{\text{def}}{=} -\rho_0 g C_s^{-2}$ where $C_s = C_s(z)$ is a function of only depth, which yields:

$$ [b]_{\Theta,S} = -\frac{g}{\rho_0} [\delta_z]_{\Theta,S} = g^2 \frac{C_s^2 - c_s^2}{c_s^2 C_s^2}. \quad (5) $$

Denoting $C_s = c_s + \Delta_c$ where $c_s^{-1} \Delta_c \ll 1$, the right-hand side (RHS) of Equation 5 becomes:

$$ g^2 \frac{(c_s + \Delta_c)^2 - c_s^2}{c_s^2 C_s^2} \approx \frac{g^2}{c_s^2} \left[ \left( 1 + \frac{2\Delta_c}{c_s} \right) - 1 \right] = \frac{2g^2 \Delta_c}{c_s C_s^2} \sim O(10^{-6}). \quad (6) $$

Hence, so long as $C_s \gtrsim c_s$, monotonicity is assured while removing a large portion of compressibility, that is, the iso-surfaces of $b$ become close to neutral surfaces. In practice, we chose $C_s$ to be larger than the maximum sound speed at each depth by $10^{-5} \text{ m s}^{-1}$ over the entire ensemble in order to avoid a singularity (viz. $[b]_{\sim z, \Theta,S} = 0$). With $C_s$ determined, integrating for $\tilde{\rho}$ gives:

$$ \tilde{\rho} = - \int_z^0 \frac{\rho_0 g}{C_s} dz + \rho_0, \quad (7) $$

which reduces to $\tilde{\rho}|_{z=0} = \rho_0$. The buoyancy equation using the in-situ density anomaly becomes:

$$ \frac{D}{Dt} b = b \dot{\Theta} + b \dot{S} + b \frac{Dz}{Dt} \quad (8) $$

$$ = B + w g^2 \frac{C_s^2 - c_s^2}{c_s^2 C_s^2}, \quad (9) $$

where $B \stackrel{\text{def}}{=} b \dot{\Theta} + b \dot{S}$, and $\dot{\Theta}$ and $\dot{S}$ are the net diabatic contributions on potential temperature and practical salinity respectively, which we approximate by diagnosing off-line the sum of harmonic and biharmonic diffusion below the mixed layer using the 5-day averaged outputs of $\Theta$ and $S$. We summarize the RHS of (9) as the dia-surface velocity $\sigma \stackrel{\text{def}}{=} B + w g^2 \frac{c_s^2 - c_s^2}{c_s^2 C_s^2}$.

A further requirement of the TWA framework is that the pressure anomaly defined by such buoyancy coordinate translates into a body force in the buoyancy coordinate
$$ \nabla_h \tilde{\phi}(z) \mapsto \nabla_h \tilde{\phi}(\tilde{b}) = \tilde{\nabla}_h m, \qquad (10) $$

where the subscript (·)$_h$ represents the horizontal gradient and $\tilde{\nabla}_h = (\partial_{\tilde{x}}, \partial_{\tilde{y}})$. Using in-situ buoyancy anomaly, the pressure anomaly becomes:

$$ \tilde{\phi}(z) = \int_{\sim} b \, dz. \qquad (11) $$

The (·) is used to denote that the pressure anomaly is defined by the in-situ buoyancy anomaly. The pressure anomaly for a Boussinesq hydrostatic fluid, on the other hand, is:

$$ \phi(z) = \int \frac{g}{\rho_0} (\rho - \rho_0) \, dz. \qquad (12) $$

Since $\tilde{\rho}$ is only a function of depth, the horizontal gradient of the two remain identical ($\nabla_h \tilde{\phi} = \nabla_h \phi$) and Equation 10 holds (We note that Equation 10 does not hold for pressure anomaly defined by potential density when the EOS is non-linear, and while more elaborate techniques may improve the neutrality of $\delta$, the relation to the dynamics is non-trivial for other density variables such as neutral and orthobaric densities.). The use of in-situ density anomaly to define the buoyancy coordinate maintains the desirable properties of a unique, statically stable vertical coordinate and a simple hydrostatic balance ($\sigma = \zeta_b = -m_{bb}$) while removing roughly 99% of the effect of compressibility basin wide at each depth ($\frac{g^2(c_x^2-c_x^2)}{g^2c_x^2} \approx \frac{2c_x\Delta c}{c_x^2} \sim O(10^{-2})$). For a non-linear EOS, a material conservation of potential vorticity (PV) and non-acceleration conditions do not exist (cf. Vallis, 2017, Chapter 4). Discussion regarding the energetics are given in Appendix A.

The raw simulation outputs were in geopotential coordinates so we first remapped all of the variables in Equations 1 and 2 onto 55 buoyancy levels spread across the range of $\tilde{b} \in (-0.196, -0.287)$ m s⁻² (with the mathematical formulation of $\delta = \delta_0 + A_\delta \frac{\tanh(\tau)-\tanh(0)}{\tanh(\tau_{max})-\tanh(0)}$ where $\delta_0 = 20$ kg m⁻³, $A_\delta = 9.2$ kg m⁻³, and $\tau \in [0, 2)$ in order to account for the abyssal weak stratification):

$$ (u, v, b, \nabla_h \tilde{\phi}, \Theta, S, \omega) (t, x, y, z) \mapsto (u, v, \zeta, \tilde{\nabla}_h m, \Theta, S, \omega) (\tilde{t}, \tilde{x}, \tilde{y}, \tilde{b}) \qquad (13) $$

using the fastjmd95 Python package to compute the in-situ density and its partial derivatives (Abernathey, 2020), and the xgcm Python package (Abernathey et al., 2021; Busecke & Abernathey, 2020; Jones et al., 2020) which allows for coordinate remapping consistent with the finite-volume discretization of MITgcm. The horizontal velocity vector becomes $\mathbf{u}\mathbf{i} + \mathbf{v}\mathbf{j} \mapsto \mathbf{u}\mathbf{e}_1 + \mathbf{v}\mathbf{e}_2$. For the horizontal pressure anomaly gradient, we re-computed the pressure anomaly using the 5-day averaged outputs and have invoked the identity (10). In the case where the buoyancy contour outcrops for some members, we treat it by making the layer thickness vanish ($\Delta\zeta = 0$) and carry on with our TWA analysis. This is consistent with the boundary treatment of Young (2012) where he notes that buoyancy contours intersecting the boundary to be continued just beneath the surface.

## 4. Results

We start by showing the time series of domain-averaged horizontal kinetic energy (KE) and potential temperature, and an arbitrary buoyancy iso-surface (Figure 2). Figure 2a shows the simulation has a prominent seasonal cycle with the KE and temperature both peaking in summer. In Figure 2, we also show the residual-mean fields on 3 January 1967, the first day of the year of output we analyze. The depth of the buoyancy level shown in Figure 2c is below the ensemble-mean mixed-layer depth (MLD; Figure 2b) basin wide where diabatic effects are small, but is shallow enough to capture the imprint of the Gulf Stream; the iso-surface shoals drastically across the latitude of ~38°N where the separated Gulf Stream is situated (Figure 2d). The ensemble-mean MLD was computed as the depth at which the potential density computed from ensemble-mean temperature and salinity fields increased by 0.03 kg m⁻³ from the density at 10 m depth ($\overline{\text{MLD}} \stackrel{\text{def}}{=} \text{MLD}(\Theta, \overline{\text{S}})$; de Boyer Montégut et al., 2004). The residual-mean KE field (MKE, $K^\# \stackrel{\text{def}}{=} |\hat{\mathbf{u}}|^2/2$; Figure 2d) shows the characteristic features of the Gulf Stream, North Brazil Current and equatorial undercurrent. The North Brazil Current, although having large values in KE, shows no imprint on the buoyancy depth (Figure 2c). The residual-mean Rossby number ($\text{Ro}^\# \stackrel{\text{def}}{=} f^{-1}(\hat{v}_{\tilde{x}} - \hat{u}_{\tilde{y}})$)
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Time series of the domain-averaged total KE (black) and potential temperature (red) for the 48 ensemble members between 15°S and 50°N. The thick lines show the ensemble mean and the thin lines each ensemble member (a). The ensemble-mean MLD on 3 January 1967 and depth of the iso-surface of buoyancy $\bar{b} = -0.26 \, \text{m s}^{-2}$ (b, c). The residual-mean KE ($K^*$) and Rossby number ($Ro^*$) on the same buoyancy surface (d, e).</figcaption></figure>

is smaller than unity over most of the Atlantic basin (Figure 2e), indicating that the residual-mean flow in the interior is balanced in our model with the exception of regions with energetic currents, for example, the Gulf Stream, loop current in the Gulf of Mexico and the North Brazil Current. Near the equator, the Coriolis parameter becomes small leading to large Rossby numbers. The kinematics of discretizing the gradients in buoyancy coordinates are given in Appendix B. We now move on to examine the eddy feedback onto the (residual) mean flow. Hereon, we drop the prefix “residual” unless required for clarity.
## 4.1. The Eliassen-Palm Flux

The E-P flux tensor (E) in the TWA framework (Equations 1 and 2) is:

$$ \mathbf{E} = \begin{pmatrix} \widehat{u''u''} + \frac{1}{2\bar{\sigma}}\bar{\zeta}^{\prime 2} & \widehat{u''v''} & 0 \\ \widehat{v''u''} & \widehat{v''v''} + \frac{1}{2\bar{\sigma}}\bar{\zeta}^{\prime 2} & 0 \\ \widehat{\sigma''u''} + \frac{1}{\bar{\sigma}}\bar{\zeta}'m_{\bar{x}}' & \widehat{\sigma''v''} + \frac{1}{\bar{\sigma}}\bar{\zeta}'m_{\bar{y}}' & 0 \end{pmatrix} \quad (14) $$

where $(\cdot)'' = (\cdot) - (\hat{\cdot})$ and $(\cdot)' = (\cdot) - (\bar{\cdot})$ are the residual from the thickness-weighted and ensemble averages respectively (Aoki, 2014; Maddison & Marshall, 2013; Ringler et al., 2017). The two are related via the (eddy-induced) bolus velocity (Greatbatch, 1998; McDougall & McIntosh, 2001):

$$ \mathbf{u}'' = \mathbf{u} - \frac{\bar{\sigma}\bar{u}}{\bar{\sigma}} = \bar{u} + \mathbf{u}' - \frac{(\bar{\sigma} + \sigma')(\bar{u} + \mathbf{u}')}{\bar{\sigma}} \quad (15) $$

$$ = \mathbf{u}' + \frac{\bar{\sigma}'\mathbf{u}'}{\bar{\sigma}}. \quad (16) $$

We show each term in Equation 14 in Figure 3. The eddy momentum flux $\widehat{u''v''}$ is often associated with barotropic processes in analogy to atmospheric jets (Figure 3a; Aoki et al., 2016; Chan et al., 2007; Jamet et al., 2021; Vallis, 2017, Chapter 15). The zonal and meridional eddy momentum flux $(\widehat{u''}^2, \widehat{v''}^2)$ exchange momentum between the eddies and mean flow, that is, to accelerate or decelerate the Gulf Stream as they affect the horizontal shear upon taking their gradients. The term due to the vertical displacement of buoyancy layer $(\frac{1}{2\bar{\sigma}}\bar{\zeta}'^2)$ is related to the eddy potential energy (EPE; cf. Equations A15–A17). The interfacial form stress $(\bar{\zeta}'\bar{V}_h m')$; Figures 3e and 3f) often associated with baroclinic instability is “deceivingly” orders of magnitude smaller than the other terms. However, it is the divergence of the E-P flux and not the flux itself that goes into the momentum equations, and the horizontal $(\bar{V}_h)$ and vertical gradient $(\delta_b)$ differ by roughly O (10⁶). The contribution from the diabatic and compressibility effects (i.e., the terms with $\bar{\sigma}$) were smaller than the interfacial form stress by another order of magnitude or more in the subtropics (not shown). It is quite surprising that the signals in the equatorial undercurrent region, although having relatively high KE (Figure 2d), are significantly smaller than in the Gulf Stream and North Brazil Current regions, virtually not visible in Figure 3. This implies that the mean flow dominates over the eddies in the equatorial region.

Writing out the E-P flux divergence in Equations 1 and 2 gives:

$$ -\bar{e}_1 \cdot (\bar{V} \cdot \mathbf{E}) = -\bar{\sigma}^{-1} \left( \left[ \bar{\sigma} \left( \widehat{u''u''} + \frac{1}{2\bar{\sigma}}\bar{\zeta}^{\prime 2} \right) \right]_{\bar{x}} + \left[ \bar{\sigma}\widehat{v''u''} \right]_{\bar{y}} + \left[ \bar{\sigma} \left( \widehat{\sigma''u''} + \frac{1}{\bar{\sigma}}\widehat{\zeta}'m_{\bar{x}}' \right) \right]_{\bar{b}} \right) \quad (17) $$

$$ = -\bar{\sigma}^{-1} \left( \left[ \bar{\sigma}u''u'' + \frac{\bar{\zeta}^{\prime 2}}{2} \right]_{\bar{x}} + \left[ \bar{\sigma}v''u'' \right]_{\bar{y}} + \left[ \bar{\sigma}\bar{\sigma}''u'' + \bar{\zeta}'m_{\bar{x}}' \right]_{\bar{b}} \right), \quad (18) $$

$$ \stackrel{\text{def}}{=} - (E_{\bar{x}}^{00} + E_{\bar{y}}^{10} + E_{\bar{b}}^{20}) \quad (19) $$

$$ -\bar{e}_2 \cdot (\bar{V} \cdot \mathbf{E}) = -\bar{\sigma}^{-1} \left( \left[ \bar{\sigma} \widehat{u''v''} \right]_{\bar{x}} + \left[ \bar{\sigma} \left( \widehat{v''v''} + \frac{1}{2\bar{\sigma}} \bar{\zeta}^{\prime 2} \right) \right]_{\bar{y}} + \left[ \bar{\sigma} \left( \widehat{\sigma''v''} + \frac{1}{\bar{\sigma}} \widehat{\zeta}'m_{\bar{x}}' \right) \right]_{\bar{b}} \right) \quad (20) $$

$$ = -\bar{\sigma}^{-1} \left( \left[ \bar{\sigma}u''v'' \right]_{\bar{x}} + \left[ \bar{\sigma}v''v'' + \frac{\bar{\zeta}^{\prime 2}}{2} \right]_{\bar{y}} + \left[ \bar{\sigma}\bar{\sigma}''v'' + \bar{\zeta}'m_{\bar{y}}' \right]_{\bar{b}} \right), \quad (21) $$

$$ \stackrel{\text{def}}{=} - (E_{\bar{x}}^{01} + E_{\bar{y}}^{11} + E_{\bar{b}}^{21}). \quad (22) $$

As the signal in the North Atlantic basin is the largest in the separated Gulf Stream region (Figure 3), we show each term in the E-P flux divergence north of 25°N (Figure 4). The large signal is consistent with Jamet et al. (2021) where they found the subtropical gyre to be a Fofonoff-like inertial circulation (Fofonoff, 1981), and that the separated jet was where the energy input to the gyre from surface winds was predominantly lost to
<figure><img src="image_6.png" /><figcaption>**Figure 3.** The residual-mean Ertel potential vorticity normalized by the local Coriolis parameter ($\Pi^#/f \stackrel{\text{def}}{=} \sigma^{-1}(1+\text{Ro}^#)$) (a) and terms in the E-P flux tensor (b–f) on 3 January 1967 on the iso-surface of buoyancy as in Figure 2. Note the scaling factors on panels (a, e, and f).</figcaption></figure>

eddies. The convergence of interfacial form stress ($E_b^{20}, E_b^{21}$) becomes larger than the convergence of the eddy momentum flux terms due to cross correlation in the zonal and meridional momentum ($E_y^{10}, E_x^{01}$), which are the smallest amongst the three terms in the E-P flux convergence (Figures 4b and 4c). The contribution from the terms with dia-surface velocity ($\omega''$) was roughly two-orders of magnitude smaller than the other terms in
<figure><img src="image_7.png" /><figcaption>**Figure 4.** The terms in the convergence of E-P flux tensor on 3 January 1967 on the iso-surface of $\bar{b} = -0.26 \, \text{m s}^{-2}$ (a–f). Positive values (red shadings) indicate the eddies fluxing momentum to the mean flow and vice versa. The panels are laid out so that summing up the top three rows per column yields the total zonal $(-\bar{e}_1 \cdot (\bar{\nabla} \cdot \mathbf{E}))$ (g) and meridional E-P flux divergence $(-\bar{e}_2 \cdot (\bar{\nabla} \cdot \mathbf{E}))$ (h) respectively. The contours in gray shading east of 285°E indicate the 400, 300 and 200 m depth of the buoyancy surface.</figcaption></figure>

the E-P flux convergence in the adiabatic interior (not shown), which supports the neutrality of $\delta$ to define the buoyancy surfaces. Right at the separation of the Gulf Stream west of 290°E and around 36°N, the convergence of eddy momentum flux and potential energy ($E_{\bar{x}}^{00}$, $E_{\bar{y}}^{11}$), and interfacial form stress ($E_{\bar{b}}^{20}$, $E_{\bar{b}}^{21}$) tend to counteract each other; in the zonal direction, the eddy momentum flux and potential energy convergence tends to decelerate.
ate the Gulf Stream while the interfacial form stress convergence tends to accelerate it (Figures 4a and 4e). The repeating positive and negative features further downstream are roughly on the scales of the Rossby deformation radius, consistent with Uchida, Deremble, et al. (2021) where they diagnosed the E-P flux convergence from a 101-member quasi-geostrophic (QG) double-gyre ensemble. In the meridional direction, the eddy momentum flux and potential energy convergence also tend to smooth out the Gulf Stream (decelerate the jet in the subpolar gyre by injecting northward momentum, and southward momentum in the subtropical gyre) while the interfacial form stress convergence tends to sharpen it (Figures 4d and 4f). The similar order of magnitude between $E_x^{00}$, $E_y^{11}$ and $E_b^{20}$, $E_b^{21}$ is in contrast, however, from a fully developed QG jet within a wind-driven double-gyre circulation where the interfacial form stress convergence dominated the E-P flux convergence (Uchida, Deremble, et al., 2021). While this does not provide as proof that the Gulf Stream in primitive equation models deviates from quasi geostrophy, the disagreement is consistent with previous studies arguing that western boundary currents, which are on the order of O (100 km) in the across-jet direction but O (1000 km) in the along-jet direction, may not be well approximated by QG dynamics, which is isotropic in its formulation (Grooms et al., 2011; Jamet et al., 2021). Further examinations, however, are required to quantify the level of deviation.

We now examine further details in the separated Gulf Stream, a region where eddies have been shown to modulate the mean flow structure (e.g. Chassignet & Xu, 2021; Cronin, 1996), as seasonal means in order to capture representative features. Winter is defined as the months of January, February, March, and summer as July, August, September. Upon separation, the zonal E-P flux convergence tends to decelerate the Gulf Stream. The repeating features of positive and negative values for the zonal component of the E-P flux convergence persist and are likely associated to the jet meandering (Figures 5a and 5c). In the meridional direction, we again see positive values on the northern flank of the separated Gulf Stream and negative on its southern flank (Figures 5b and 5d). This north-south dipole feature is likely associated with the gradient of the eddy energy, and may be trivial as the energy naturally maximizes near the center of the jet. The zonal and meridional component of the E-P flux convergence can jointly be interpreted to force the Gulf Stream to migrate northwards (decelerate the jet northwards in the subtropical gyre on the North flank of the separated Gulf Stream and southwards in the subpolar gyre; Figures 4b and 4d) although this largely being contained west of 310°E. The interpretation of poleward jet migration is consistent with the zonal E-P flux convergence where the overall structure of the forcing of the zonal equation is a deceleration on one side of the Gulf Stream and an acceleration on the other; the eddy momentum flux in the zonal momentum equation decelerate both the core and the flanks immediately downstream of Cape Hatteras (Figure 4a) and alternate further downstream (a signature of meandering) while the form drag term partially cancels this (Figure 4e). East of 310°E, the E-P flux convergence tends to shift the North Atlantic Current east and southwards in the open ocean, while northwards closer to the continental rise (Figure 4h). Examining the meridional transect averaged over the zonal extent of 290°E–305°E where the separated Gulf Stream is roughly zonal (Figure 2d), the separated Gulf Stream can be identified with the steep shoaling of the iso-surfaces of buoyancy between 36°N and 40°N (Figures 5e–5h). The overall magnitude and reversal in sign at the core of the jet (around 37.5°N) with diminishing amplitude with depth for the zonal E-P flux convergence during winter (–*e*₁ ⋅ (*V* ⋅ *E*); Figures 5g, 6a, and 6b) is roughly in agreement with Ringler et al. (2017, their Figure 6 where the sign convention in Equation 17 is reversed from ours for the eddy forcing term and their units are in [m s⁻¹ day⁻¹]) where they diagnosed an idealized zonally re-entrant jet. It is interesting to note, however, that the vertical structure of the E-P flux convergence is much smoother and barotropic during the summer with a consistent deceleration of the jet on its northern flank and acceleration on its southern flank (Figures 5g, 6e, and 6f). We note that such seasonal features may be specific to the year of 1967, and the temporal evolution of the E-P flux convergence should be addressed in a dedicated study. We leave this for further work, focusing here on the TWA implementation for a realistic model.

In Figure 6, we show the vertical profile of the seasonal E-P flux convergence along with each component in Equations 17 and 20 area averaged over the zonal extent of 290°E–305°E. The E-P flux convergence closely follows that of the interfacial form stress convergence (i.e., baroclinic instability) with the Reynolds stress due to cross correlation between the zonal and meridional eddy momentum ($E_y^{10}$, $E_x^{01}$; orange lines) taking the smallest magnitude. The amplitude of interfacial form stress convergence is larger near the surface (viz. larger buoyancy values), which is expected from the seasonal surface forcing affecting the isopycnal tilt and hence baroclinicity of the surface flow. The meridional smoothing of the separated Gulf Stream is also apparent from the vertical profiles with the meridional E-P flux convergence taking negative values on the southern flank of the jet and positive values on the northern flank. The convergence of eddy momentum flux and potential energy tends to
<figure><img src="image_8.png" /><figcaption>Figure 5.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Vertical profile of the area-averaged, seasonal zonal and meridional E-P flux convergence north and south of the separated Gulf Stream over the zonal extent 290°E–305°E. The area averaging is separated between 35°N–37.5°N and 37.5°N–40°N. The top panels show the seasonal mean for winter and bottom for summer.</figcaption></figure>

mirror that of interfacial form stress (blue and green lines in Figure 6). This counteracting balance is consistent with what Aiki et al. (2016, the terms $\partial_x R^x$ and $\partial_z (R^z + F_a^+)$ in their Figures 5a and 6) found in the Kuroshio extension region.

## 4.2. The Ergodic Assumption

In this section, we replace the averaging operator with the temporal mean of the 50 years of output $(\hat{\cdot})^t, (\cdot)^{t' \text{ def}} = (\cdot) - (\cdot)^t$ from a single arbitrary realization (realization 00 to be specific) to examine the ergodic assumption and compare with our TWA results. Realization area averaging is separated area averaging is separated 00 was taken from a 24-member ensemble originally designed for a different study (Jamet et al., 2019b). The 48 members discussed above were constructed by adding 24 members to the first 5 years of this data set. The TWA operator now becomes $\hat{\sigma}^t \hat{\sigma}^{-1} \overline{\sigma(\cdot)}$ and eddies $(\cdot)^{t'} \text{ def} = (\cdot) - (\cdot)^t$. The maximum sound speed per depth ($C_s$) was recomputed for the 50 years of realization 00 in remapping the coordinate system. Although the averaging operator is now along the time dimension, we note that this is different from the Temporal-Residual Mean (TRM) framework developed by McDougall and McIntosh (2001) in the sense that we proceed with our analysis

**Figure 5.** The seasonal mean of the zonal and meridional E-P flux convergence for winter and summer of 1967 (a–d). The contours in gray shading indicate the 400, 300 and 200 m depth of the buoyancy surface. The zonal-mean transect between 290°E and 305°E of the E-P flux convergence is shown in colored shading and ensemble-mean depth in black contours (e–h). The iso-surface of buoyancy used through Figures 2–4 is shown as the gray dashed line. The masked out region north of 30°N near the surface during winter is where the iso-surfaces of buoyancy outcrop across all ensemble members. We see that more buoyancy surfaces outcrop during winter.
<figure><img src="image_10.png" /><figcaption>**Figure 7.** The terms in the climatological convergence of E-P flux tensor on the iso-surface of $\tilde{b} = -0.26 \, \text{m s}^{-2}$ from realization Figures 4 and 5a–5d 00 (a–d). We do not show the terms due to the Reynolds stress ($\overline{\left(\frac{u''v''}{t}\right)}$) as they were negligible compared to the other terms, and omit the superscript $t$ on variables with primes to avoid the clutter. Climatology of the total zonal ($-\bar{e}_1 \cdot (\tilde{\nabla} \cdot \mathbf{E})$) and meridional E-P flux divergence ($-\bar{e}_2 \cdot (\tilde{\nabla} \cdot \mathbf{E})$) respectively (e, f). The contours in gray shading east of 285°E indicate the 400, 300, and 200 m depth of the buoyancy surface.</figcaption></figure>

in buoyancy coordinate. The hope of applying the ergodic assumption to a temporally varying system, as we have shown in previous sections, is that for a sufficiently long time series, such sub- and inter-annual variability will cancel out with only the stationary feature being extracted in the “mean” flow.

In Figure 7, we show the climatological E-P flux convergence from realization area averaging is separated area averaging is separated 00. In other words, all time scales shorter than 50 years are now relegated to the eddies. While having similar spatial structures to Figures 4 and 5a–5d, they are more spread out with less detail. In particular, the seasonality is obscured by the climatological mean of 50 years and becomes similar to the summer-time of the 48-member ensemble (Figures 5c and 5d). In other words, the wintertime signal seen with the ensemble diagnostics (Figures 5a and 5b) are not well captured by the climatological E-P fluxes convergence. This could either suggest that such signal are peculiar to the year 1967 we analyzed with our 48-member ensemble,
<figure><img src="image_11.png" /><figcaption>**Figure 8.** The climatological zonal-mean transect between 290°E and 305°E of the E-P flux convergence is shown in colored shading and ensemble-mean depth in black contours from realization 00 (a, b). The iso-surface of buoyancy used in Figure 7 is shown as the gray dashed line.</figcaption></figure>

or that summertime signals may have a stronger imprint on the residual time mean. Considering the 50-year time scale of averaging, the signals that emerge in the climatological E-P flux convergence are likely due to transient eddies while the standing eddies would be included in the mean flow. The climatological zonal-mean transect also resemble the ensemble summertime albeit with weaker amplitude (Figures 5e–5h, and 8) where the eddies tend to zonally decelerate the separated Gulf Stream on its northern flank and accelerate it on its southern flank (Figure 8a). In the meridional direction, the eddies tend to decelerate the subpolar gyre on the northern flank of the separated Gulf Stream and the subtropical gyre on its southern flank (Figure 8b).

Taking the climatological time mean of 50 years of output is perhaps the most conservative definition of the mean flow under ergodicity. We, therefore, now loosen the temporal averaging to a climatological annual cycle in defining the residual mean flow. In doing so, we chunk the 50 years into 50 annual segments and take their average to produce a single segment of ~365 days. Namely, we treat each year as an individual realization of the ocean, generating a pseudo 50-member year-long ensemble (hereon pseudo-ensemble for short). The eddies are now defined as fluctuations about this climatological annual cycle. In Figure 9, we show the MKE on a buoyancy level on January 3 with similar depths diagnosed from the ensemble and pseudo-ensemble. While the maximum MKE amplitudes are similar, the mean flow is more spread out in the pseudo-ensemble. This likely comes from the different paths the Gulf Stream takes resulting as a response to different yearly atmospheric states, which get averaged all together. In other words, while the degrees of freedom are similar between the ensemble (48 members) and pseudo-ensemble (50 members assuming a decorrelation time scale of a year), the ensemble mean captures the oceanic response to the atmospheric state specific to 1967. The pseudo-ensemble, on the other hand, implies that 50 years are not sufficient for the “eddies” to emerge as a coherent signal upon averaging for a climatological annual cycle and the mean flow incorporates the signal of atmospheric interannual, decadal and low-frequency variability.

The imprint of fluctuations from each year onto the MKE domain averaged over the depths of ~50–500 m ($\bar{b} \in (-0.25, -0.26)$) result in its seasonality to differ from the ensemble mean; the pseudo-ensemble takes its maximum around March while the ensemble around August (black solid and dashed lines in Figure 9c respectively). However, the seasonality in the area averaged MKE from the pseudo-ensemble on $\bar{b} = -0.26 \text{ m s}^{-2}$ shows a summertime maximum (black dot-dashed line in Figure 9c). This implies that the discrepancy between $K^H$ and $K^{Ht}$ results from the surface ocean being sensitive to the atmospheric state while being less so in the interior. Indeed, the domain averaged eddy KE (EKE; see Appendix A for definition) diagnosed from the ensemble shows a maximum during winter when the surface ocean is more susceptible to baroclinic instability due to atmospheric cooling (red line in Figure 9c; Uchida et al., 2017). We conclude that in the process of creating a climatological annual cycle, we convolute the oceanic response to different atmospheric states (i.e., interannual variability) and contaminate the eddy-mean flow decomposition. The oceanic mean flow conflated with atmospheric variability also imprints itself onto the E-P flux convergence for the climatological winter and summer as we show in
<figure><img src="image_12.png" /><figcaption>**Figure 9.** The (residual) mean KE on 3 January from the ensemble ($K^{\#}$) and pseudo-ensemble ($K^{\#t}$) on buoyancy levels with similar depth (a, b). The regions with outcropping buoyancy surface are masked out. The colors indicate the MKE and contours in gray scaling show the depths for 50, 100, 200, and 300 m. Time series of domain averaged MKE ($K^{\#}$ and $K^{\#t}$) in black plotted against the left y axis and EKE ($\hat{K}$) in red plotted against the right y axis (c) The domain was taken over the horizontal extent shown in panels (a, b). Note the difference in magnitudes of order on the y axes.</figcaption></figure>

Figure 10, which arguably looks noisier than Figures 5a–5d particularly north of the 300 m depth contour in the subpolar gyre.

**5. Discussion and Summary**

By running a 48-member ensemble run of the North Atlantic Ocean at mesoscale-permitting resolution (1/12°) partially coupled to the atmosphere, we have shown that the TWA framework can be employed successfully in diagnosing eddy-mean flow interactions in a realistic ocean simulation. In doing so, we have introduced a new buoyancy variable for a realistic EOS, which is approximately neutral and dynamically consistent; both characteristics are necessary for the TWA analysis (Stanley, 2018). The ensemble approach negates the necessity for any temporal averaging in defining the residual-mean flow; we are able to exclude any temporal variability, such as seasonal and interannual fluctuations, from the eddy term and extract the intrinsic variability of the ocean. We show that the Eliassen-Palm (E-P) flux convergence (i.e., negative divergence), which encapsulates the eddy feedback onto the mean flow (Maddison & Marshall, 2013), tends to accelerate the Gulf Stream northwards on its northern flank ($-\bar{e}_2 \cdot (\tilde{V} \cdot E) > 0$) and decelerate it on its southern flank ($-\bar{e}_2 \cdot (\tilde{V} \cdot E) < 0$; Figures 5b, 5d, 5f, and 5h); that is, the eddies can be interpreted to force the Gulf Stream to migrate northwards on 3 January 1967. However, a more detailed examination of the mechanism of poleward jet migration will likely necessitate studies using idealized simulations where each dynamical mechanism is easier to parse out (cf. Chemke & Kaspi, 2015).
<figure><img src="image_13.png" /><figcaption>**Figure 10.** The E-P flux convergence from the pseudo-ensemble for the climatological winter and summer on the iso-surface of $\bar{b} = -0.26 \, \text{m s}^{-2}$. The contours in gray shading east of 285°E indicate the 400, 300, and 200 m depth of the buoyancy surface.</figcaption></figure>

Here, we have documented a dynamically-consistent implementation of the TWA framework for a realistic ocean simulation and the E-P flux convergence diagnosed in the context of oceanic ensemble simulations.

Modeling studies with varying spatial resolution have shown that the Gulf Stream tends to overshoot northwards and the North Atlantic Current (NAC) flows too zonally in coarse resolution models (e.g., Chassignet & Xu, 2017, 2021; Lévy et al., 2010). The overshooting may partially be attributable to eddy feedback being insufficiently resolved at mesoscale-permitting resolutions, in addition to unresolved submesoscale boundary layer processes (e.g., Renault et al., 2016). In particular, it would be interesting to see whether further increasing the model resolution would increase the amplitude of baroclinic instability near the surface ($E_b^{20}, E_b^{21}$) and convergence of eddy momentum flux and potential energy in the interior ($E_x^{00}, E_y^{11}$), which tend to accelerate the jet southward in the subpolar gyre and decelerate it southward in the subtropical gyre upon the Gulf Stream separation west of 290°E (i.e., shift the jet southwards) as we see from their annual means (Figure 11). The same could be said for a better representation of the NAC path where the eddies in our model tend to flux northward momentum into the mean flow and hence allow for its north-eastward turn near the continental rise of the Grand Banks (Figures 4 and 5). Although it is beyond the scope of this study, the significance of baroclinic processes will likely increase with resolution as mixed-layer instability becomes better resolved (Boccaletti et al., 2007; Capet et al., 2008a, 2008b; Su et al., 2018; Uchida et al., 2019; Yang et al., 2021).

We have also examined the often assumed ergodicity in decomposing the eddy and mean flow by replacing the averaging operator with a 50-year time mean for a single realization within the ensemble. To some extent, the agreement between Figures 4, 5, 7, and 11 implies that the ensemble size of 48 is able to extract the eddy signals that emerge at mesoscale-permitting resolution. The difference between the ensemble and 50-year climatology of an arbitrary realization amongst the ensemble (realization 00), on the other hand, likely comes from seasonal, interannual and decadal variability, and transient eddies, which are obscured in the climatological view. Loosening the time mean to a climatological annual cycle for the mean flow, on the other hand, convolutes the oceanic
<figure><img src="image_14.png" /><figcaption>**Figure 11.** The annual mean of the coverage of eddy momentum flux and potential energy, and interfacial form stress for $\tilde{b} = -0.26 \text{ m s}^{-2}$ (a–d) The contours in gray shading east of 285°E indicate the 400, 300, and 200 m depth of the buoyancy surface. The annual and zonal mean transect between 290°E and 305°E of the E-P flux convergence is shown in colored shading and ensemble-mean depth in black contours (e–h). The iso-surface of buoyancy used through panels (a–d) is shown as the gray dashed line.</figcaption></figure>
response to interannual variability in the atmospheric forcing and contaminates the eddy-mean flow decomposition (Figure 9). This is consistent with Aiki and Richards (2008) where they found the energy stored in the mean and eddy flow to change depending on the duration of the temporal averaging applied. While it is not our intention to claim whether defining the mean flow via a time mean is appropriate or not for realistic simulations, our results imply that one should be mindful of what goes into defining the mean flow and consequently the eddies.

Lastly, ensemble modeling has shown us that a small perturbation such as eddies to the non-linear system can lead to very different states of the ocean and climate (e.g., Bessières et al., 2017; Fedele et al., 2021; Jamet et al., 2019b; Lorenz, 1963; Maher et al., 2019; Uchida, Deremble, & Penduff, 2021). In light of this, we argue that it is important to consider the full spatiotemporal variability of the ocean. The ensemble framework allows one to capture the space-time varying eddy-mean flow interaction and not just its climatological state.

## Appendix A: Energetics Under a Non-Linear Equation of State

In this Appendix, we derive the energetics in a similar manner to Aiki et al. (2016) but in a framework consistent with the ensemble formalism and a realistic equation of state (EOS). The thickness-weighted average (TWA) residual-mean horizontal momentum equation in geopotential coordinates neglecting dissipation is (Ringler et al., 2017; Young, 2012):

$$ \hat{\mathbf{u}}_t + \mathbf{v}^\# \cdot \nabla \hat{\mathbf{u}} + f \mathbf{k} \times \hat{\mathbf{u}} = -\nabla_h \phi^\# - \bar{\mathbf{e}} \cdot (\nabla \cdot \mathbf{E}), \quad (A1) $$

where $\mathbf{v}^\# \stackrel{\text{def}}{=} \hat{\mathbf{u}}\hat{\mathbf{i}} + \hat{\mathbf{v}}\hat{\mathbf{j}} + \omega^\#\mathbf{k}$ and $\phi^\# \stackrel{\text{def}}{=} \bar{m}(\tilde{t}, \tilde{x}, \tilde{y}, b^\#(t, x, y, z)) + b^\#z$ are the residual-mean velocity and hydrostatic pressure anomaly. It is important to keep in mind that the “z” here is the ensemble averaged depth of an iso-surface of buoyancy, viz. $z = \zeta(\tilde{t}, \tilde{x}, \tilde{y}, b^\#(t, x, y, z))$. The residual-mean kinetic energy (MKE; $K^\# = |\hat{\mathbf{u}}|^2/2$) budget becomes:

$$ \begin{aligned} K_t^\# + \mathbf{v}^\# \cdot \nabla K_t^\# &= -\hat{\mathbf{u}} \cdot \nabla_h \phi^\# - \hat{\mathbf{u}} \cdot [\bar{\mathbf{e}} \cdot (\nabla \cdot \mathbf{E})] \\ &= -\hat{\mathbf{u}} \cdot \nabla_h \phi^\# - \omega^\# \phi_z^\# + \omega^\# \phi_z^\# - \hat{\mathbf{u}} \cdot [\bar{\mathbf{e}} \cdot (\nabla \cdot \mathbf{E})] \\ &= -\mathbf{v}^\# \cdot \nabla \phi^\# + \omega^\# b^\# - \hat{\mathbf{u}} \cdot [\bar{\mathbf{e}} \cdot (\nabla \cdot \mathbf{E})]. \end{aligned} \quad (A2) $$

We can now define the dynamic enthalpy for the mean state in a similar manner to McDougall (2003) and Young (2010):

$$ h^\# \stackrel{\text{def}}{=} \int_{\Phi_0}^{\Phi^\#} \frac{b^\#}{g} d\Phi^{\#\prime} = \int_{z}^{0} b^\# dz^\prime, \quad (A3) $$

where $\Phi^\# = \Phi_0 - gz$ is the dynamically non-active part of the hydrostatic pressure to be consistent with the Boussinesq approximation. Note that $h^\#$ is not a function of the TWA temperature and salinity ($\hat{\Theta}, \hat{S}$) due to non-linearities in the EOS, that is, $b(\hat{\Theta}, \hat{S}, z) \neq b(\Theta, S, z) = b = \tilde{b} = b^\#$. While there exist a temperature and salinity variable to evaluate the material derivative of $h^\#$ since an EOS exists for $b^\#$, it is unclear whether they can be analytically expressed for a non-linear EOS. We, therefore, express the material derivative of $h^\#$ as:

$$ \begin{aligned} \frac{D^\#}{Dt} h^\# &= h_z^\# \frac{D^\# z}{Dt} + \mathcal{H}^\# \\ &= -\omega^\# b^\# + \mathcal{H}^\#, \end{aligned} \quad (A4) $$

where $\mathcal{H}^\#$ carries the net sum of the diabatic and non-linear effects. Thus, the residual-mean total energy equation becomes:

$$ \frac{D^\#}{Dt} (K^\# + h^\#) = -\nabla \cdot \mathbf{v}^\# \phi^\# + \mathcal{H}^\# - \hat{\mathbf{u}} \cdot [\bar{\mathbf{e}} \cdot (\nabla \cdot \mathbf{E})], \quad (A5) $$

where we have invoked $\nabla \cdot \mathbf{v}^\# = 0$.
On the other hand, the total KE budget remapped onto buoyancy coordinate is:

$$ \frac{DK}{Dt} = -\tilde{\nabla} \cdot \mathbf{v}\phi + w\tilde{b}, \qquad (A6) $$

where $\mathbf{v} \stackrel{\text{def}}{=} v^1 \mathbf{e}_1 + v^2 \mathbf{e}_2 + v^3 \mathbf{e}_3 = u\mathbf{e}_1 + v\mathbf{e}_2 + \left(\sigma + \frac{\zeta_i}{\sigma}\right) \mathbf{e}_3$ and $\tilde{\nabla} \cdot \mathbf{v} = \sigma^{-1} \left[ (\sigma v^1)_{\bar{x}} + (\sigma v^2)_{\bar{y}} + (\sigma v^3)_{\bar{b}} \right] (=0)$ is the three-dimensional divergence. Unlike the residual-mean dynamic enthalpy, the definition of the total dynamic enthalpy is straight forward (Young, 2010):

$$ h = \int_{\zeta}^{0} b(\Theta, S, \zeta') d\zeta', \qquad (A7) $$

yielding:

$$ \frac{D}{Dt}(K+h) = -\tilde{\nabla} \cdot \mathbf{v}\phi + \mathcal{H}, \qquad (A8) $$

where $\mathcal{H} \stackrel{\text{def}}{=} h_{\Theta} \frac{D\Theta}{Dt} + h_S \frac{DS}{Dt}$. Terms due to non-linearity in the EOS do not emerge in the definition of $\mathcal{H}$ as Equation A8 is not averaged. Ensemble averaging after thickness weighting Equation A8 gives:

$$ \begin{aligned} \overline{\sigma \frac{D}{Dt}(K+h)} &= -\overline{\sigma \tilde{\nabla} \cdot \mathbf{v}\phi} + \overline{\sigma \mathcal{H}} \\ &= -\overline{\sigma \tilde{\nabla} \cdot \mathbf{v}\phi} + \overline{\sigma \hat{\mathcal{H}}}, \end{aligned} \qquad (A9) $$

The total KE can be expanded as:

$$ \begin{aligned} K &= \frac{1}{2} |\hat{\mathbf{u}} + \mathbf{u}''|^2 \\ &= \frac{|\hat{\mathbf{u}}|^2}{2} + \frac{|\mathbf{u}''|^2}{2} + \hat{u}u'' + \hat{v}v'' \end{aligned} \qquad (A10) \\ \stackrel{\text{def}}{=} K^{\#} + \mathcal{K} + \hat{u}u'' + \hat{v}v'', \end{aligned} $$

so plugging in Equation A10, and keeping in mind that $\overline{(\hat{\cdot})} = (\hat{\cdot})$ and $\overline{\sigma(\cdot)''} = 0$, each term on the left-hand side (LHS) of Equation A9 can be written as:

$$ \begin{aligned} \overline{\sigma \frac{DK}{Dt}} &= \overline{\sigma (K_{\bar{i}} + uK_{\bar{x}} + vK_{\bar{y}} + \sigma K_{\bar{b}})} \\ &= (\overline{\sigma K})_{\bar{i}} + (\overline{\sigma uK})_{\bar{x}} + (\overline{\sigma vK})_{\bar{y}} + (\overline{\sigma \sigma K})_{\bar{b}} \end{aligned} \qquad (A11) \\ = \overline{\sigma} \left[ \frac{D^{\#}}{Dt} (K^{\#} + \hat{\mathcal{K}}) + \tilde{\nabla} \cdot (\mathbf{J}^{\mathcal{K}} + \hat{u}\mathbf{J}^{\mathcal{u}} + \hat{v}\mathbf{J}^{\mathcal{v}}) \right], \end{aligned} $$

where $\hat{\mathcal{K}}$ is the eddy kinetic energy (EKE), and $\mathbf{J}^{\mathcal{K}} \stackrel{\text{def}}{=} \overline{u''\mathcal{K}\mathbf{e}_1 + v''\mathcal{K}\mathbf{e}_2 + \overline{\sigma''\mathcal{K}\mathbf{e}_3}$, $\mathbf{J}^{\mathcal{u}} \stackrel{\text{def}}{=} \overline{u''\overline{v''}\mathbf{e}_1 + v''\overline{u''}\mathbf{e}_2 + \overline{\sigma''\overline{v''}\mathbf{e}_3}$ are the eddy fluxes of kinetic energy, eddy zonal and meridional velocities respectively, and

$$ \begin{aligned} \overline{\sigma \frac{D\mathbf{h}}{Dt}} &= \overline{\sigma (h_{\bar{i}} + u h_{\bar{x}} + v h_{\bar{y}} + \sigma h_{\bar{b}})} \\ &= (\overline{\sigma h})_{\bar{i}} + (\overline{\sigma u h})_{\bar{x}} + (\overline{\sigma v h})_{\bar{y}} + (\overline{\sigma \sigma h})_{\bar{b}} \\ &= (\overline{\sigma \hat{h}})_{\bar{i}} + [\overline{\sigma} (\hat{u} \hat{h} + u'' \overline{h''})]_{\bar{x}} + [\overline{\sigma} (\hat{v} \hat{h} + v'' \overline{h''})]_{\bar{y}} + [\overline{\sigma} (\hat{\sigma} \hat{h} + \overline{\sigma'' h''})]_{\bar{b}} \\ &= \overline{\sigma} \left( \frac{D^{\#}}{Dt} \hat{h} + \tilde{\nabla} \cdot \mathbf{J}^{\mathcal{h}} \right), \end{aligned} \qquad (A12) $$

where $\mathbf{J}^{\mathcal{h}} \stackrel{\text{def}}{=} \overline{u''\overline{h''}\mathbf{e}_1 + v''\overline{h''}\mathbf{e}_2 + \overline{\sigma''\overline{h''}\mathbf{e}_3}$ is the eddy flux of fluctuations in dynamic enthalpy, and we have used the relation $\overline{\sigma\phi\theta} = \overline{\sigma} (\overline{\hat{\phi}\theta} + \overline{\hat{\phi}''\theta''})$ (Equation 72 in Young, 2012). Hence, combining Equations A11 and A12, Equation A9 becomes:
$$ \frac{D^{\#}}{Dt} (K^{\#} + \hat{K} + \hat{h}) = -\tilde{\nabla} \cdot (\mathbf{J}^{K} + \mathbf{J}^{h} + \hat{u}\mathbf{J}^{u} + \hat{v}\mathbf{J}^{v}) - \tilde{\nabla} \cdot \widehat{\mathbf{v}\phi} + \hat{H}. \quad (A13) $$

Subtracting Equations A5 from A13 yields the eddy energy budget:

$$ \frac{D^{\#}}{Dt} (\hat{K} + \hat{h} - h^{\#}) = - (\tilde{\nabla} \cdot \widehat{\mathbf{v}\phi} - \nabla \cdot \mathbf{v}^{\#} \phi^{\#}) - \tilde{\nabla} \cdot (\mathbf{J}^{K} + \mathbf{J}^{h} + \hat{u}\mathbf{J}^{u} + \hat{v}\mathbf{J}^{v}) + \hat{H} - \mathcal{H}^{\#} + \hat{u} \cdot [\hat{\mathbf{e}} \cdot (\nabla \cdot \mathbf{E})]. \quad (A14) $$

Equations A5 and A14 are the relations derived by Aoki (2014) but for a non-linear EOS and non-zero dia-surface velocity where the residual-mean flow and eddies exchange energy via the E-P flux divergence and residual vertical buoyancy flux due to non-linearities in the EOS. It is perhaps interesting to note that $h''$ is not the eddy potential energy (EPE; $\hat{H} \stackrel{\text{def}}{=} \hat{h} - h^{\#}$ in Equation A14) and they are related to one another as $h'' = h - (h^{\#} + \hat{H})$.

For a linear EOS, the EPE can be rewritten as:

$$ \hat{H} = -b^{\#} (\hat{\zeta} - \bar{\zeta}) = -b^{\#} \frac{\overline{\sigma' \zeta'}}{\overline{\sigma}}, \quad (A15) $$

by taking advantage of $\hat{h} = -\tilde{b}\hat{\zeta}$, $h^{\#} = -b^{\#}\bar{\zeta}$ and $b = \tilde{b} = b^{\#}(t, x, y, \bar{\zeta}(\tilde{t}, \tilde{x}, \tilde{y}, \tilde{b}))$. Equation A15 provides the physical intuition of EPE being defined as the difference between potential energy at the TWA depth ($\hat{\zeta}$) and ensemble-mean depth ($\bar{\zeta}$). In a similar manner, we can also derive:

$$ h'' = -\tilde{b}(\zeta - \hat{\zeta}) = -\tilde{b}\zeta'', \quad (A16) $$

and hence, $\overline{h''} = -\hat{H}$. Assuming the background buoyancy frequency can be defined as the inverse of ensemble-mean thickness (viz. $\overline{\sigma}^{-1} \sim N^2$) leads to further manipulation of EPE:

$$ \begin{aligned} \hat{H} &\sim -b^{\#} N^2 \overline{\zeta' \zeta'}_{\tilde{b}}' = -b^{\#} N^2 \left( \frac{\overline{\zeta'^2}}{2} \right)^{\frac{1}{2}} \\ &= -N^2 \left[ \left( b^{\#} \frac{\overline{\zeta'^2}}{2} \right)^{\frac{1}{2}} \right]_{\tilde{b}}^{\frac{1}{2}} \end{aligned} \quad (A17) $$

where the last term in Equation A17 further reduces to the available potential energy under quasi-geostrophic approximation ($b't \sim N^2\zeta t$). The first-term on the RHS of Equation A17 vanishes upon volume integration pending on boundary conditions (i.e., rigid lid and a flat bottom).

## Appendix B: Kinematics of Discretization

As in Figure B1, imagine $u_1$ and $u_2$ are on the same buoyancy contour. The relation between the two is:

$$ u_2 \approx u_1 + u_\infty \Delta x + u_\zeta \Delta \zeta. \quad (B1) $$

<figure><img src="image_15.png" /><figcaption>Figure B1. Schematic of discretized gradients.</figcaption></figure>
Now,

$$ 
\begin{align} 
u_{\tilde{x}} & \stackrel{\text{def}}{=} u_x + \frac{\Delta\zeta}{\Delta x} \sigma^{-1} u_{\tilde{b}} \nonumber \\ 
& = u_x + \frac{\Delta\zeta}{\Delta x} u_{\zeta} \tag{B2} \\ 
& = \frac{u_2 - u_1}{\Delta x} \quad (\because \text{equation(B1))}, \nonumber 
\end{align} 
$$

so once all of the variables are remapped onto the buoyancy coordinate from geopotential, the discretized horizontal gradients can be taken along the original Cartesian grid. The gradients on the model outputs were taken using the xgcm Python package (Abernathey et al., 2021; Busecke & Abernathey, 2020). In order to minimize the computational cost, we took the ensemble mean first whenever possible, for example, $\bar{\sigma} = \partial_{\tilde{b}}\zeta = \partial_{\tilde{b}}\bar{\zeta}$, $\tilde{V}_h\bar{\sigma} = \partial_{\tilde{b}}\tilde{V}_h\bar{\zeta}$ etc. The gradient operators commuting with the ensemble mean is also the case for the perturbations, that is,

$$ \tilde{V}_h(\bar{m} + m') = \tilde{V}_h m = \overline{\tilde{V}_h m} + (\tilde{V}_h m)' \quad (B3) $$

Hence, $\tilde{V}_h m' = (\tilde{V}_h m)'$ (cf. Maddison & Marshall, 2013, Section 2.3 in their paper).

## Data Availability Statement

The simulation outputs are available on the Florida State University cluster (http://ocean.fsu.edu/~qjamet/share/data/Uchida2021/). Python scripts used for the off-line diagnosis are available on Github (doi:10.5281/zenodo.5823419). Uchida acknowledges the xarray (Hoyer et al., 2021) and dask (Rocklin et al., 2021) Python package developers for parallelizing the analysis.

## References

* Abernathey, R. P. (2020). fastjmd95: Numba implementation of Jackett & McDougall (1995) ocean equation of state. Retrieved from https://github.com/xgcm/fastjmd95

* Abernathey, R. P., Busecke, J., Smith, T., Banihirwe, A., Fernandes, F., Bourbeau, J., et al. (2021). xgcm: General circulation model postprocessing with xarray. https://doi.org/10.5281/zenodo.3634752

* Aiki, H., & Richards, K. J. (2008). Energetics of the global ocean: The role of layer-thickness form drag. *Journal of Physical Oceanography*, 38(9), 1845–1869. https://doi.org/10.1175/2008JPO3820.1

* Aiki, H., Zhai, X., & Greatbatch, R. J. (2016). Energetics of the global ocean: The role of mesoscale eddies. *World Scientific*, 109–134. https://doi.org/10.1142/9789814696623_0004

* Ajayi, A., Le Sommer, J., Chassignet, E., Molines, J.-M., Xu, X., Albert, A., & Cosme, E. (2020). Spatial and temporal variability of the North Atlantic eddy field from two kilometric-resolution ocean models. *Journal of Geophysical Research: Oceans*, 125, https://doi.org/10.1029/2019JC015827

* 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

* Andrews, D. G. (1983). A finite-amplitude Eliassen-Palm theorem in isentropic coordinates. *Journal of the Atmospheric Sciences*, 40(8), 1877–1883. https://doi.org/10.1175/1520-0469(1983)040<1877:afaep>2.0.co;2

* Aoki, K. (2014). A constraint on the thickness-weighted average equation of motion deduced from energetics. *Journal of Marine Research*, 72(5), 355–382. https://doi.org/10.1357/002224014815469886

* Aoki, K., Kubokawa, A., Furue, R., & Sasaki, H. (2016). Influence of eddy momentum fluxes on the mean flow of the Kuroshio extension in a 1/10 ocean general circulation model. *Journal of Physical Oceanography*, 46(9), 2769–2784. https://doi.org/10.1175/JPO-D-16-0021.1

* Arbic, B. K., Polzin, K. L., Scott, R. B., Richman, J. G., & Shriver, J. F. (2013). On eddy viscosity, energy cascades, and the horizontal resolution of gridded satellite altimeter products. *Journal of Physical Oceanography*, 43(2), 283–300. https://doi.org/10.1175/JPO-D-11-0240.1

* Bachman, S. D. (2019). The GM+E closure: A framework for coupling backscatter with the Gent and McWilliams parameterization. *Ocean Modelling*, 136, 85–106. https://doi.org/10.1016/j.ocemod.2019.02.006

* Bachman, S. D., Anstey, J. A., & Zanna, L. (2018). The relationship between a deformation-based eddy parameterization and the LANS-α turbulence model. *Ocean Modelling*, 126, 56–62. https://doi.org/10.1016/j.ocemod.2018.04.007

* Bachman, S. D., Fox-Kemper, B., & Bryan, F. O. (2015). A tracer-based inversion method for diagnosing eddy-induced diffusivity and advection. *Ocean Modelling*, 86, 1–14. https://doi.org/10.1016/j.ocemod.2014.11.006

* Balwada, D., Smith, S. K., & Abernathey, R. P. (2018). Submesoscale vertical velocities enhance tracer subduction in an idealized Antarctic Circumpolar Current. *Geophysical Research Letters*, 45(18), 9790–9802. https://doi.org/10.1029/2018GL079244

* Berloff, P. (2018). Dynamically consistent parameterization of mesoscale eddies. Part III: Deterministic approach. *Ocean Modelling*, 127, 1–15. https://doi.org/10.1016/j.ocemod.2018.04.009

* Bessières, L., Leroux, S., Brankart, J.-M., Molines, J.-M., Moine, M.-P., Bouttier, P.-A., et al. (2017). Development of a probabilistic ocean modelling system based on NEMO 3.5: Application at eddying resolution. *Geoscientific Model Development*, 10(3), 1091–1106. https://doi.org/10.5194/gmd-10-1091-2017

## Acknowledgments

This research was funded by the French “Make Our Planet Great Again” (MOPGA) initiative managed by the Agence Nationale de la Recherche under the Programme d’Investissement d’Avenir, with the reference ANR-18-MPGA-0002. High-performance computing resources on Cheyenne (doi:10.5065/D6RX99HX) used for running the ensembles were provided by NCAR’s Computational and Information Systems Laboratory, sponsored by the National Science Foundation, under the university large allocation UFSU0011. This study is a contribution to the “Assessing the Role of forced and internal Variability for the Ocean and climate Response in a changing climate” (ARVOR) project supported by the French “Les Enveloppes Fluides et l’Environnement” (LEFE) program. Dewar acknowledges support from the NSF grants OCE-1829856 and OCE-1941963. Balwada acknowledges support from the NSF grant OCE-1756882. We thank the editor Stephen Griffies, along with Laure Zanna and four other anonymous reviewers for their constructive and careful reviews, which lead to significant improvements in the manuscript. We would like to thank Geoffrey Stanley and Trevor McDougall for insightful discussions regarding the buoyancy coordinate for a realistic EOS. We also acknowledge Edward Peirce and Kelly Hirai for maintaining the cluster upon which the data analyses were conducted. Uchida would also like to acknowledge the Eddy Parametrization Reading Group at Columbia University during his Ph.D. where he was first introduced to the literature of mesoscale eddy parametrizations and has kept him intrigued ever since.
Bire, S., & Wolfe, C. L. (2018). The role of eddies in buoyancy-driven eastern boundary currents. *Journal of Physical Oceanography*, 48(12), 2829–2850. https://doi.org/10.1175/JPO-D-18-0040.1

Boccaletti, G., Ferrari, R., & Fox-Kemper, B. (2007). Mixed layer instabilities and restratification. *Journal of Physical Oceanography*, 37(9), 2228–2250. https://doi.org/10.1175/jpo3101.1

Bühler, O. (2014). *Waves and mean flows*. Cambridge University Press. https://doi.org/10.1017/CBO9781107478701

Busecke, J., & Abernathey, R. P. (2020). CMIP6 without the interpolation: Grid-native analysis with Pangeo in the cloud. In *Earthcube annual meeting*. Retrieved from https://github.com/earthcube2020/ec20_busecke_etal

Capet, X., McWilliams, J. C., Molemaker, M. J., & Shchepetkin, A. F. (2008a). Mesoscale to submesoscale transition in the California current system. Part I: Flow structure, eddy flux, and observational tests. *Journal of Physical Oceanography*, 38(1), 29–43. https://doi.org/10.1175/2007JPO3671.1

Capet, X., McWilliams, J. C., Molemaker, M. J., & Shchepetkin, A. F. (2008b). Mesoscale to submesoscale transition in the California Current system. Part III: Energy balance and flux. *Journal of Physical Oceanography*, 38(10), 2256–2269. https://doi.org/10.1175/2008JPO3810.1

Cessi, P., & Wolfe, C. L. (2013). Adiabatic eastern boundary currents. *Journal of Physical Oceanography*, 43(6), 1127–1149. https://doi.org/10.1175/JPO-D-12-0211.1

Chan, C. J., Plumb, R. A., & Cerovecki, I. (2007). Annular modes in a multiple migrating zonal jet regime. *Journal of the Atmospheric Sciences*, 64(11), 4053–4068. https://doi.org/10.1175/2007JAS2156.1

Chassignet, E. P., & Xu, X. (2017). Impact of horizontal resolution (1/12° to 1/50°) on Gulf Stream separation, penetration, and variability. *Journal of Physical Oceanography*, 47(8), 1999–2021. https://doi.org/10.1175/JPO-D-17-0031.1

Chassignet, E. P., & Xu, X. (2021). On the importance of high-resolution in large-scale ocean models. *Advances in Atmospheric Sciences*, 38, 1–1634. https://doi.org/10.1007/s00376-021-0385-7

Chemke, R., & Kaspi, Y. (2015). Poleward migration of eddy-driven jets. *Journal of Advances in Modeling Earth Systems*, 7(3), 1457–1471. https://doi.org/10.1002/2015MS000481

Cronin, M. (1996). Eddy-mean flow interaction in the gulf stream at 68 W. Part II: Eddy forcing on the time-mean flow. *Journal of Physical Oceanography*, 26(10), 2132–2151. https://doi.org/10.1175/1520-0485(1996)026<2132:emfiit>2.0.co;2

de Boyer Montégut, C., Madec, G., Fischer, A. S., Lazar, A., & Iudicone, D. (2004). Mixed layer depth over the global ocean: An examination of profile data and a profile-based climatology. *Journal of Geophysical Research*, 109(C12). https://doi.org/10.1029/2004JC002378

Deremble, B., Wienders, N., & Dewar, W. (2013). CheapAML: A simple, atmospheric boundary layer model for use in ocean-only model calculations. *Monthly Weather Review*, 141(2), 809–821. https://doi.org/10.1175/MWR-D-11-00254.1

de Szoeke, R. A., & Bennett, A. F. (1993). Microstructure fluxes across density surfaces. *Journal of Physical Oceanography*, 23(10), 2254–2264. https://doi.org/10.1175/1520-0485(1993)023<2254:mfads>2.0.co;2

de Szoeke, R. A., & Springer, S. R. (2009). The materiality and neutrality of neutral density and orthobaric density. *Journal of Physical Oceanography*, 39(8), 1779–1799. https://doi.org/10.1175/2009JPO4042.1

Fairall, C. W., Bradley, E. F., Hare, J. E., Grachev, A. A., & Edson, J. B. (2003). Bulk parameterization of air–sea fluxes: Updates and verification for the COARE algorithm. *Journal of Climate*, 16(4), 571–591. https://doi.org/10.1175/1520-0442(2003)016<0571:bpoasf>2.0.co;2

Fedele, G., Penduff, T., Pierini, S., Alvarez-Castro, M. C., Bellucci, A., & Masina, S. (2021). Interannual to decadal variability of the Kuroshio extension: Analyzing an ensemble of global hindcasts from a dynamical system viewpoint. *Climate Dynamics*, 1, 975–992. https://doi.org/10.1007/s00382-021-05751-7

Fofonoff, N. (1981). *The Gulf Stream system. Evolution of physical oceanography*. MIT Press. Retrieved from https://ocw.mit.edu/ans7870/textbooks/Wunsch/Edited/Chapter4.pdf

Gallimore, R. G., & Johnson, D. R. (1981). The forcing of the meridional circulation of the isentropic zonally averaged circumpolar vortex. *Journal of the Atmospheric Sciences*, 38(3), 583–599. https://doi.org/10.1175/1520-0469(1981)038<0583:tfotmc>2.0.co;2

Greatbatch, R. J. (1998). Exploring the relationship between eddy-induced transport velocity, vertical momentum transfer, and the isopycnal flux of potential vorticity. *Journal of Physical Oceanography*, 28(3), 422–432. https://doi.org/10.1175/1520-0485(1998)028<0422:etrb ei>2.0.co;2

Griffies, S. M., Winton, M., Anderson, W. G., Benson, R., Delworth, T. L., Dufour, C. O., et al. (2015). Impacts on ocean heat from transient mesoscale eddies in a hierarchy of climate models. *Journal of Climate*, 28(3), 952–977. https://doi.org/10.1175/JCLI-D-14-00353.1

Grooms, I., Julien, K., & Fox-Kemper, B. (2011). On the interactions between planetary geostrophy and mesoscale eddies. *Dynamics of Atmospheres and Oceans*, 51(3), 109–136. https://doi.org/10.1016/j.dynatmoce.2011.02.002

Guillaumin, A. P., & Zanna, L. (2021). Stochastic-deep learning parameterization of ocean momentum forcing. *Journal of Advances in Modeling Earth Systems*, 13(9), e2021MS002534. https://doi.org/10.1029/2021MS002534

Hoyer, S., Abernathey, R. P., Hamman, J., Bovy, B., Maussion, F., & Fujii, K. (2021). *xarray: N-d labeled arrays and datasets*. https://doi.org/10.5281/zenodo.598201

Jackett, D. R., & McDougall, T. J. (1995). Minimal adjustment of hydrographic profiles to achieve static stability. *Journal of Atmospheric and Oceanic Technology*, 12(2), 381–389. https://doi.org/10.1175/1520-0426(1995)012<0381:maohpt>2.0.co;2

Jackett, D. R., & McDougall, T. J. (1997). A neutral density variable for the world's oceans. *Journal of Physical Oceanography*, 27(2), 237–263. https://doi.org/10.1175/1520-0485(1997)027<0237:andvft>2.0.co;2

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

Jamt, Q., Dewar, W. K., Wienders, N., & Deremble, B. (2019a). Fast warming of the surface ocean under a climatological scenario. *Geophysical Research Letters*, 46(7), 3871–3879. https://doi.org/10.1029/2019GL082336

Jamt, Q., Dewar, W. K., Wienders, N., & Dersemble, B. (2019b). Spatiotemporal patterns of chaos in the Atlantic Overturning circulation. *Geophysical Research Letters*, 46(13), 7509–7517. https://doi.org/10.1029/2019GL082552

Jamt, Q., Dewar, W. K., Wienders, N., Dersemble, B., Close, S., & Penduff, T. (2020). Locally and remotely forced subtropical AMOC variability: A matter of time scales. *Journal of Climate*, 33, 5155–5172. https://doi.org/10.1175/JCLI-D-19-0844.1

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

Jones, S. C., Busecke, J., Uchida, T., & Abernathey, R. P. (2020). Vertical regridding and remapping of CMIP6 ocean data in the cloud. In *Earth-cube annual meeting*. Retrieved from https://github.com/earthcube2020/ec20_jones_etal

Juricke, S., Danilov, S., Koldunov, N., Oliver, M., & Sidorenko, D. (2020). Ocean kinetic energy backscatter parametrization on unstructured grids: Impact on global eddy-permitting simulations. *Journal of Advances in Modeling Earth Systems*, 12(1). https://doi.org/10.1029/2019MS001855

Kitsios, V., Frederiksen, J. S., & Zidikheri, M. J. (2013). Scaling laws for parameterisations of subgrid eddy–eddy interactions in simulations of oceanic circulations. *Ocean Modelling*, 68, 88–105. https://doi.org/10.1016/j.oocemod.2013.05.001
Kjellsson, J., & Zanna, L. (2017). The impact of horizontal resolution on energy transfers in global ocean models. *Fluids*, 2(3), 45. https://doi.org/10.3390/fluids2030045

Klocker, A., McDougall, T. J., & Jackett, D. R. (2009). A new method for forming approximately neutral surfaces. *Ocean Science*, 5(2), 155–172. https://doi.org/10.5194/os-5-155-2009

Lang, Y., Stanley, G. J., McDougall, T. J., & Barker, P. M. (2020). A pressure-invariant neutral density variable for the World's Oceans. *Journal of Physical Oceanography*, 50, 1–3604. https://doi.org/10.1175/JPO-D-19-0321.1

Leroux, S., Penduff, T., Bessières, L., Molines, J.-M., Brankart, J.-M., Sérazin, G., et al. (2018). Intrinsic and atmospherically forced variability of the AMOC: Insights from a large-ensemble ocean hindcast. *Journal of Climate*, 31(3), 1183–1203. https://doi.org/10.1175/JCLI-D-17-0168.1

Lévy, M., Klein, P., Tréguier, A.-M., Iovino, D., Madec, G., Masson, S., & Takahashi, K. (2010). Modifications of gyre circulation by sub-mesoscale physics. *Ocean Modelling*, 34(1–2), 1–15. https://doi.org/10.1016/j.ocemod.2010.04.001

Lévy, M., Resplandy, L., Klein, P., Capet, X., Iovino, D., & Éthé, C. (2012). Grid degradation of submesoscale resolving ocean models: Benefits for offline passive tracer transport. *Ocean Modelling*, 48, 1–9. https://doi.org/10.1016/j.ocemod.2012.02.004

Lorenz, E. N. (1963). Deterministic nonperiodic flow. *Journal of the Atmospheric Sciences*, 20(2), 130–141. https://doi.org/10.1175/1520-0469(1963)020<0130:dnf>2.0.co;2

Maddison, J. R., & Marshall, D. P. (2013). The Eliassen-Palm flux tensor. *Journal of Fluid Mechanics*, 729, 69–102. https://doi.org/10.1017/jfm.2013.259

Maher, N., Milinski, S., Suarez-Gutierrez, L., Botzet, M., Dobrynin, M., Kornblueh, L., et al. (2019). The Max Planck Institute grand ensemble: Enabling the exploration of climate system variability. *Journal of Advances in Modeling Earth Systems*, 11(7), 2050–2069. https://doi.org/10.1029/2019MS001639

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

Marshall, J., Hill, C., Perelman, L., & Adcroft, A. (1997). Hydrostatic, quasi-hydrostatic, and nonhydrostatic ocean modeling. *Journal of Geophysical Research*, 102(C3), 5733–5752. https://doi.org/10.1029/96JC02776

McDougall, T. J. (2003). Potential enthalpy: A conservative oceanic variable for evaluating heat content and heat fluxes. *Journal of Physical Oceanography*, 33(5), 945–963. https://doi.org/10.1175/1520-0485(2003)033<0945:peacov>2.0.co;2

McDougall, T. J., & Jackett, D. R. (2005). An assessment of orthobaric density in the global ocean. *Journal of Physical Oceanography*, 35(11), 2054–2075. https://doi.org/10.1175/JPO2796.1

McDougall, T. J., & McIntosh, P. C. (2001). The temporal-residual-mean velocity. Part II: Isopycnal interpretation and the tracer and momentum equations. *Journal of Physical Oceanography*, 31(5), 1222–1246. https://doi.org/10.1175/1520-0485(2001)031<1222:trmvp>2.0.co;2

Molines, J.-M., Barnier, B., Penduff, T., Treguier, A., & Le Sommer, J. (2014). ORCA12. L46 climatological and interannual simulations forced with DFS4. 4: GJM02 and MJM88. drakkar group experiment rep (Technical Report). In Gdri-drakkar-2014-03-19 (p. 50). Retrieved from http://www.drakkar-ocean.eu/publications/reports/orca12referenceexperiments2014

Montgomery, R. B. (1937). A suggested method for representing gradient flow in isentropic surfaces. *Bulletin of the American Meteorological Society*, 18(6–7), 210–212. https://doi.org/10.1175/1520-0477-18.6-7.210

Olbers, D., Willebrand, J., & Eden, C. (2012). *Ocean dynamics*. Springer Science & Business Media.

Perezhogin, P. (2019). Deterministic and stochastic parameterizations of kinetic energy backscatter in the NEMO ocean model in double-gyre configuration. *IOP Conference Series: Earth and Environmental Science*, 386, 012025. https://doi.org/10.1088/1755-1315/386/1/012025

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

Ringler, T., Saenz, J. A., Wolfram, P. J., & Van Roekel, L. (2017). A thickness-weighted average perspective of force balance in an idealized circumpolar current. *Journal of Physical Oceanography*, 47(2), 285–302. https://doi.org/10.1175/JPO-D-16-0096.1

Rocklin, M., Crist-Harif, J., Bourbeau, J., Augspurger, T., Durant, M., & Signell, J., et al. (2021). *dask: Parallel computing with task scheduling*. Retrieved from https://dask.org/

Sasaki, H., Klein, P., Qiu, B., & Sasai, Y. (2014). Impact of oceanic-scale interactions on the seasonal modulation of ocean dynamics by the atmosphere. *Nature Communications*, 5(1), 1–8. https://doi.org/10.1038/ncomms6636

Schubert, R., Jonathan, G., 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, 2573–2589. https://doi.org/10.1175/JPO-D-19-0311.1

Sérazin, G., Jaymond, A., Leroux, S., Penduff, T., Bessières, L., Llovel, W., et al. (2017). A global probabilistic study of the ocean heat content low-frequency variability: Atmospheric forcing versus oceanic chaos. *Geophysical Research Letters*, 44(11), 5580–5589. https://doi.org/10.1002/2017GL073026

Stammer, D. (1997). Global characteristics of ocean variability estimated from regional TOPEX/POSEIDON altimeter measurements. *Journal of Physical Oceanography*, 27(8), 1743–1769. https://doi.org/10.1175/1520-0485(1997)027<1743:gcoove>2.0.co;2

Stanley, G. J. (2018). *Tales from topological oceans*. Doctoral dissertation, University of Oxford. Retrieved from https://ora.ox.ac.uk/objects/uuid:414355e6-ba26-4004-bc71-51e4fa5fb1bb

Stanley, G. J. (2019). Neutral surface topology. *Ocean Modelling*, 138, 88–106. https://doi.org/10.1016/j.ocemod.2019.01.008

Su, Z., Wang, J., Klein, P., Thompson, A. F., & Menemenlis, D. (2018). Ocean submesoscales as a key component of the global heat budget. *Nature Communications*, 9(1), 775. https://doi.org/10.1038/s41467-018-02983-w

Tailleux, R. (2016). Generalized patched potential density and thermodynamic neutral density: Two new physically based quasi-neutral density variables for ocean water masses analyses and circulation studies. *Journal of Physical Oceanography*, 46(12), 3571–3584. https://doi.org/10.1175/JPO-D-16-0072.1

Uchida, T., Abernathey, R. P., & Smith, S. K. (2017). Seasonality of eddy kinetic energy in an eddy permitting global climate model. *Ocean Modelling*, 118, 41–58. https://doi.org/10.1016/j.ocemod.2017.08.006

Uchida, T., Balwada, D., Abernathey, R. P., McKinley, G. A., Smith, S. K., & Lévy, M. (2019). The contribution of submesoscale over mesoscale eddy iron transport in the open Southern Ocean. *Journal of Advances in Modeling Earth Systems*, 11, 3934–3958. https://doi.org/10.1029/2019MS001805

Uchida, T., Deremble, B., Dewar, W. K., & Penduff, T. (2021). Diagnosing the Eliassen-Palm flux from a quasi-geostrophic double gyre ensemble. In *Earthcube annual meeting*. https://doi.org/10.5281/zenodo.5496375

Uchida, T., Deremble, B., & Penduff, T. (2021). The seasonal variability of the ocean energy cycle from a quasi-geostrophic double gyre ensemble. *Fluids*, 6(6), 206. https://doi.org/10.3390/fluids6060206

Uchida, T., Deremble, B., & Popinet, S. (2022). Deterministic model of the eddy dynamics for a midlatitude ocean model. *Journal of Physical Oceanography*. https://doi.org/10.1175/JPO-D-21-0217.1
Uchida, T., Jamet, Q., Poje, A. C., & Dewar, W. K. (2021). An ensemble-based eddy and spectral analysis, with application to the Gulf Stream. *Journal of Advances in Modeling Earth Systems, 13*, e2021MS002692. https://doi.org/10.1029/2021MS002692

Vallis, G. K. (2017). *Atmospheric and oceanic fluid dynamics* (2nd ed.). Cambridge University Press.

Waterman, S., & Lilly, J. M. (2015). Geometric decomposition of eddy feedbacks in barotropic systems. *Journal of Physical Oceanography, 45*(4), 1009–1024. https://doi.org/10.1175/JPO-D-14-0177.1

Xu, Y., & Fu, L.-L. (2011). Global variability of the wavenumber spectrum of oceanic mesoscale turbulence. *Journal of Physical Oceanography, 41*(4), 802–809. https://doi.org/10.1175/2010JPO4558.1

Xu, Y., & Fu, L.-L. (2012). The effects of altimeter instrument noise on the estimation of the wavenumber spectrum of sea surface height. *Journal of Physical Oceanography, 42*(12), 2229–2233. https://doi.org/10.1175/JPO-D-12-0106.1

Yang, P., Jing, Z., Sun, B., Wu, L., Qiu, B., Chang, P., & Ramachandran, S. (2021). On the upper-ocean vertical eddy heat transport in the Kuroshio extension. Part I: Variability and dynamics. *Journal of Physical Oceanography, 51*(1), 229–246. https://doi.org/10.1175/JPO-D-20-0068.1

Young, W. R. (2010). Dynamic enthalpy, conservative temperature, and the seawater Boussinesq approximation. *Journal of Physical Oceanography, 40*(2), 394–400. https://doi.org/10.1175/2009JPO4294.1

Young, W. R. (2012). An exact thickness-weighted average formulation of the Boussinesq equations. *Journal of Physical Oceanography, 42*(5), 692–707. https://doi.org/10.1175/JPO-D-11-0102.1

Zanna, L., & Bolton, T. (2020). Data-driven equation discovery of ocean mesoscale closures. *Geophysical Research Letters, 47*, e2020GL088376. https://doi.org/10.1029/2020GL088376

Zanna, L., Mana, P. P., Anstey, J., David, T., & Bolton, T. (2017). Scale-aware deterministic and stochastic parametrizations of eddy-mean flow interaction. *Ocean Modelling, 111*, 66–80. https://doi.org/10.1016/j.ocemod.2017.01.004

Zhao, K., & Marshall, D. P. (2020). Ocean eddy energy budget and parameterization in the North Atlantic. In *Ocean sciences meeting*. Retrieved from https://agu.confex.com/agu/osm20/meetingapp.cgi/Paper/646129