# The Seasonal Variability of the Ocean Energy Cycle from a Quasi-Geostrophic Double Gyre Ensemble

Takaya Uchida, Bruno Deremble, Thierry Penduff

► To cite this version:

Takaya Uchida, Bruno Deremble, Thierry Penduff. The Seasonal Variability of the Ocean Energy Cycle from a Quasi-Geostrophic Double Gyre Ensemble. Fluids, 2021, 6 (6), pp.206. 10.3390/fluids6060206. insu-03275460

HAL Id: insu-03275460

https://insu.hal.science/insu-03275460

Submitted on 5 Aug 2021

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
Article

# The Seasonal Variability of the Ocean Energy Cycle from a Quasi-Geostrophic Double Gyre Ensemble

Takaya Uchida *ID*, Bruno Deremble ID and Thierry Penduff ID

Université Grenoble Alpes, CNRS, IRD, Grenoble-INP, Institut des Géosciences de l'Environnement, CS 40 700, CEDEX 9, 38058 Grenoble, France; bruno.deremble@univ-grenoble-alpes.fr (B.D.); Thierry.PENDUFF@cnrs.fr (T.P.)

* Correspondence: takaya.uchida@univ-grenoble-alpes.fr

**Abstract:** With the advent of submesoscale $O(1 \text{ km})$ permitting basin-scale ocean simulations, the seasonality of mesoscale $O(50 \text{ km})$ eddies with kinetic energies peaking in summer has been commonly attributed to submesoscale eddies feeding back onto the mesoscale via an inverse energy cascade under the constraint of stratification and Earth's rotation. In contrast, by running a 101-member, seasonally forced, three-layer quasi-geostrophic (QG) ensemble configured to represent an idealized double-gyre system of the subtropical and subpolar basin, we find that the mesoscale kinetic energy shows a seasonality consistent with the summer peak without resolving the submesoscales; by definition, a QG model only resolves small Rossby and Froude number dynamics ($O(R_o) \ll 1, O(F_r) \ll 1$) while submesoscale dynamics are associated with $O(R_o) \sim 1, O(F_r) \gtrsim 1$. Here, by quantifying the Lorenz cycle of the mean and eddy energy, defined as the ensemble mean and fluctuations about the mean, respectively, we propose a different mechanism from the inverse energy cascade. During summer, when the Western Boundary Current is stabilized and strengthened due to increased stratification, stronger mesoscale eddies are shed from the separated jet. Conversely, the opposite occurs during the winter; the separated jet destabilizes and results in overall lower mean and eddy kinetic energies despite the domain being more susceptible to baroclinic instability from weaker stratification.

**Keywords:** ocean circulation; Lorenz energy cycle; quasi-geostrophic flows; ensemble modelling

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

check for
updates

**Citation:** Uchida, T.; Deremble, B.; Penduff, T. The Seasonal Variability of the Ocean Energy Cycle from a Quasi-Geostrophic Double Gyre Ensemble. *Fluids* 2021, *6*, 206. https://doi.org/10.3390/fluids6060206

Academic Editor: Pavel S. Berloff

Received: 9 April 2021

Accepted: 27 May 2021

Published: 2 June 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

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

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

## 1. Introduction

The energy cycle of the atmospheric system, namely the energy exchange between the mean flow and fluctuations about the mean, have long been of interest due to the fluctuating flow being attributed to what is commonly known as the “weather” [1,2]. Similarly, the oceanographic community has had a long-standing interest in eddies, the weather system of the oceans [3–5]. In a seminal paper, Lorenz [2] provided a framework in understanding the eddy–mean flow interaction, a framework often referred to as the Lorenz energy cycle (herein LEC; [6]).

LEC generally decomposes the flow into four energy reservoirs: the mean and eddy available potential energy (APE) and kinetic energy (KE), respectively. The concept of APE is perhaps unique to the field of geophysical fluid dynamics where the gravitational force plays a dominant role in the governing equations. Although all geophysical fluids store gravitational potential energy, only a small fraction of it is available to generate fluid motion, hence the prefix “available”. The energy exchanges between each reservoir elucidate the balance of physical processes responsible for causing the eddy flow [5], e.g., exchanges between the mean and eddy KE are associated with barotropic instability while exchanges between the eddy APE and eddy KE are associated with baroclinic instability. Barotropic instability is generated via horizontal shear in the mean flow while baroclinic instability occurs when the effect of gravity, due to weak vertical stratification, has a similar order of magnitude as the effects of Earth’s rotation [7]. The balance between the two instabilities
results in the weather and eddies we commonly observe in the atmosphere and ocean. With the recent increase in computational power and advent of eddy resolving simulations of the ocean, there has been a growing interest in the interlinkage between the energy exchanges and temporal variability, namely seasonality, in the eddy flow [8,9].

In the context of Physical Oceanography, the eddies can be further separated into meso- and submesoscale eddies. Mesoscale eddies are roughly on the spatial scales of the first baroclinic Rossby radius of deformation ($NH/f \sim O(50 \text{ km})$) where *N* and *H* are the vertical stratification and ocean depth, respectively, and *f* is the Coriolis parameter) while submesoscale eddies are on the scale of *O*(1 km) [10]. In terms of the Rossby number ($Ro (=U/f_0L)$) and Froude number ($Fr (=U/NH)$) where *U* and *L* are the characteristic scales of velocity and length, the spatial scales translate as mesoscale dynamics being on the order of $O(Ro) \ll 1$, $O(Fr) \ll 1$, and submesoscale flows being associated with $O(Ro) \sim 1$, $O(Fr) \gtrsim 1$ [11–16]. In other words, mesoscale dynamics are more constrained by Earth's rotation and stratification, leading to the well-known phenomenon of inverse energy cascade where KE is transferred from scales about the Rossby radius to larger scales [17–19]. To what extent the framework of inverse energy cascade is applicable for scales smaller than the Rossby radius remains an open question [11,12,20].

Although there is some geographical variability [21–25], many studies using meso- and submesoscale permitting ocean simulations have attributed the seasonality in mesoscale KE to energy being transferred upscale from the submesoscales where the seasonal modulation of the mixed-layer depth leads to a strong signal [14,26–31]. Instabilities within the mixed layer are inherently submesoscale due to the reduced stratification and shallow depth scale, and are most active during late winter/early spring when the mixed-layer is the deepest [32–34]. The summertime peak in mesoscale KE has consequently been explained by the time required for the submesoscale energy to cascade upscale. Other mechanisms, such as air–sea interaction, have also been argued for the cause of mesoscale seasonality [35]. While we agree that submesoscale instabilities and air–sea interaction affect mesoscale variability, here, we examine another mechanism on the other end of the spectrum in modulating the mesoscale seasonality: the basin-scale ($O(1000 \text{ km})$) affecting the mesoscale.

In order to quantify the exchanges between the energy reservoirs, we run a seasonally forced, three-layer quasi-geostrophic (QG) ensemble with a double-gyre configuration and examine the LEC. By definition, a QG model only resolves small Rossby number dynamics based on asymptotic expansion of the governing equations [36], i.e., the simulated eddy field only consists of mesoscale variability. The background state in quasi-geostrophy can be considered as the basin-scale state. In particular, we define the mean via the ensemble mean and eddies as the fluctuations about the mean. The ensemble mean: (i) negates the ergodic assumption where one treats the temporal mean equivalent to an ensemble mean, which is questionable for a temporally varying system; (ii) removes the arbitrary temporal and/or spatial scale in defining the mean [37]; (iii) is consistent with the Reynold's definition of eddy–mean decomposition [38]; and (iv) retains the temporal, namely seasonal, variability of the LEC.

The paper is organized as follows: We describe the model configuration in Section 2 and re-derive the layered QG equations and LEC in Section 3, which will aid our discussion later on. We present our results in Section 4 and conclude in Section 5.

## 2. Model Description

We use the quasi-geostrophic (QG) configuration of the Multiple Scale Ocean Model (MSOM; [39], herein referred to as MSQG), based on the Basilisk language [40], to simulate a three-layer double-gyre flow with a rigid lid and flat bottom. No-flux conditions are applied at the lateral boundaries. The parameters used are similar to prior QG studies which examine the dynamics of a double-gyre system [3,41,42] and are summarized in Table 1. The characteristic length scale of the Rossby radius (viz. radii of mesoscale eddies) is prescribed as $L (=50 \text{ km})$ and horizontal resolution is $\sim 4 \text{ km} (=\delta_x L)$ and therefore we have
roughly 12 grid points per radius; our simulation can be considered mesoscale resolving under the numerics of a second-order Arakawa advection scheme [43–45] (we note that our kilometric resolution does not allow for the submesoscales to be permitted in our model due to the QG constraint: $O(\text{Ro}) \ll 1, O(\text{Fr}) \ll 1$).

<table><caption>**Table 1.** Parameters used to configure the three-layer QG simulation and dimensionalized characteristic scales. The bottom Ekman number is the ratio between the bottom Ekman-layer thickness and $\hat{H}_3$ and bottom friction is $\epsilon = Ek^b/(2Ro^m\hat{H}_3)$. Beta is dimensionalized as $\beta = \hat{\beta}U/L^2$ and the dimensionalized domain size is 4000 km ($=\hat{L}_0L$). The frequency of Fr translates approximately to a 360-day year ($=f_{\text{Fr}}^{-1}L/U$). The prognostic time stepping is determined via the CFL condition within values smaller than $\delta_{\hat{t}}^{\text{max}}$.</caption><thead><tr><th>Parameter</th><th>Notation</th><th>Value</th><th>Unit</th></tr></thead><tbody><tr><td>Number of horizontal grids</td><td>N</td><td>1024</td><td>-</td></tr><tr><td>Number of vertical layers</td><td>n<sub>l</sub></td><td>3</td><td>-</td></tr><tr><td>Non-dim. horizontal domain size</td><td>L̂<sub>0</sub></td><td>80</td><td>-</td></tr><tr><td>Non-dim. horizontal resolution</td><td>δ̂<sub>x</sub></td><td>N<sup>-1</sup>L̂<sub>0</sub></td><td>-</td></tr><tr><td>Background Rossby number</td><td>Ro<sup>m</sup></td><td>0.025</td><td>-</td></tr><tr><td>Non-dim. Coriolis parameter</td><td>f̂<sub>0</sub></td><td>Ro<sup>m-1</sup></td><td>-</td></tr><tr><td>Bottom Ekman number</td><td>Ek<sup>b</sup></td><td>0.004</td><td>-</td></tr><tr><td>Non-dim. surface Ekman pumping</td><td>τ̂<sub>0</sub></td><td>0.0001</td><td>-</td></tr><tr><td>Biharmonic Reynolds number</td><td>Re<sub>4</sub></td><td>4000</td><td>-</td></tr><tr><td>Non-dim. beta</td><td>β̂</td><td>0.5</td><td>-</td></tr><tr><td>Background Froude number</td><td>Fr<sub>1</sub><sup>m</sup>; Fr<sub>2</sub><sup>m</sup></td><td>0.00409959; 0.01319355</td><td>-</td></tr><tr><td>Amplitude of Fr<sub>i</sub></td><td>Â<sub>Fr1</sub>; Â<sub>Fr2</sub></td><td>0.1; 0</td><td>-</td></tr><tr><td>Non-dim. frequency of Fr<sub>i</sub></td><td>f̂<sub>Fr1</sub>; f̂<sub>Fr2</sub></td><td>62.2<sup>-1</sup>; 62.2<sup>-1</sup></td><td>-</td></tr><tr><td>Non-dim. layer thickness</td><td>Î̂<sub>1</sub>; Î̂<sub>2</sub>; Î̂<sub>3</sub></td><td>0.06; 0.14; 0.8</td><td>-</td></tr><tr><td>Non-dim. reduced gravity</td><td>ĝ'<sub>i</sub></td><td>Fr<sub>i</sub><sup>-2</sup>Ĥ<sub>i</sub><sup>†</sup></td><td>-</td></tr><tr><td>Non-dim. maximum time stepping</td><td>δ<sup>max</sup><sub>l</sub></td><td>5 × 10<sup>-2</sup></td><td>-</td></tr><tr><td>CFL condition</td><td>-</td><td>0.4</td><td>-</td></tr><tr><td>Horizontal velocity</td><td>U</td><td>0.1</td><td>[m s<sup>-1</sup>]</td></tr><tr><td>Length scale</td><td>L</td><td>50</td><td>[km]</td></tr><tr><td>Total layer thickness</td><td>H</td><td>5000</td><td>[m]</td></tr></tbody></table>

MSQG solves prognostically for the non-dimensionalized QG potential vorticity (PV):

$$ \frac{D\hat{q}}{D\hat{t}} = \hat{\mathcal{F}} + \hat{\mathcal{D}}, \qquad (1) $$

where $q = \zeta_g + \beta y - \frac{f_0}{\hat{H}}h$ is the QGPV (details are given in Appendix A; [7]) and the $\beta$-plane approximation is applied ($f = f_0 + \beta y$). $\mathcal{F}$ and $\mathcal{D}$ are the forcing and dissipative terms, and $\hat{\cdot}$ are non-dimensionalized variables. The forcing term is the wind stress curl without any buoyancy forcing at the surface, and is kept stationary with the formulation:

$$ \hat{\mathcal{F}} = \frac{\hat{\nabla}_h \times \hat{\tau}(\hat{y})}{\hat{H}_1} = -\frac{\hat{\tau}_0}{Ro^m \hat{H}_1} \sin\left(\frac{2\pi}{N}\hat{y}\right) \sin\left(\frac{\pi}{N}\hat{y}\right), \qquad (2) $$

where $\nabla_h$ is the horizontal gradient operator, and $\hat{y}$ ($\in [0.5, N - 0.5]$) is the non-dimensionalized meridional extent of the domain. Only the wind stress curl is prescribed in the model and not the wind stress itself ($\tau$) but we denote it for clarity in notation. We have kept the wind stress curl axisymmetric as low-frequency variability is not the focus of this study [46–49]. The dissipation term is implemented via a biharmonic viscosity:

$$ \hat{\mathcal{D}} = -Re_4^{-1} \hat{\nabla}_h^4 \hat{q}. \qquad (3) $$
The background stratification is defined at each layer interface via the Froude number where we enforce the seasonality by varying it in time according to:

$$ Fr_i = \frac{U}{\sqrt{g'_i H_i^\dagger}} = Fr_i^m [1 + \hat{A}_{Fr_i} \sin(2\pi \hat{f}_{Fr_i} \hat{t})]^{-1/2}, \quad (4) $$

where $H_i^\dagger = (H_i + H_{i+1})/2$, $g'_i$ is the reduced gravity and subscript $i$ is the layer index (Figure 1). We vary $Fr_1$ in time but keep $Fr_2$ stationary ($\hat{A}_{Fr_2} = 0$), which is consistent with the seasonal variability of stratification being confined in the upper few hundred meters in the real ocean [50].

<figure><img src="image_4.png" /><figcaption>Figure 1. Vertical structure of the three-layer QG model with a rigid lid and flat bottom. The layer interface displacement ($\eta_i$) is shown in the thin curvy lines and net layer thickness is $h_i = H_i + \eta_{i-1} - \eta_i$. The stream functions ($\psi_i$) are defined within each layer.</figcaption></figure>

We spin up the model for 10 years from a spun-up run with lower resolution ($N = 256$, equivalently $\delta_{\hat{x}}L \sim 15$ km) and then perturb the first-layer stream function at a single, random grid point with a perturbation on the order of ($O(10^{-5})$) to generate 100 slightly perturbed stream function fields. We use the perturbed fields as the initial conditions to generate 100 ensemble members. The surface wind stress and temporally varying background stratification are kept identical during the spin up and amongst ensemble members after the spin up. We run each ensemble member for another 10 years and for reference, we also have a control (CTRL) run without any perturbations to the initial condition; in total, we have 101 ensemble members and the CTRL run is there to show that the perturbations do not lead to a bifurcation in the dynamical regime within the 10 years of our simulation [51]. The model outputs were saved as instantaneous snapshots at every characteristic time scale ($T = L/U = 5 \times 10^5$ seconds $\sim 5.8$ days).

### 3. Derivation of the Lorenz Energy Cycle

Although the layered QG equations have been derived countless times [1,3,7,36], here, we re-derive the energy equations for a rigid-lid and flat-bottom three-layer QG model with a seasonally varying background stratification, in which the latter leads to some subtleties. In the remainder of the study, we only discuss dimensionalized variables. We start off with the order Rossby number relative vorticity equation for a given layer $i$ ($\in [1,3]$; Figure 1) neglecting viscous and external forcing terms:

$$ \begin{aligned} \partial_t \zeta_{g,i} + u_{g,i} \partial_x \zeta_{g,i} + v_{g,i} \partial_y \zeta_{g,i} + \beta v_{g,i} &= -f_0 (\partial_x u_{a,i} + \partial_y v_{a,i}) \\ &= f_0 \partial_z w_{a,i}, \end{aligned} \quad (5) $$
which are derived by taking the cross product of the momentum Equations (A6) and (A7). The subscripts g and a denote the geostrophic and ageostrophic components, respectively (e.g., $\zeta = \zeta_g + \zeta_a$). We denote the partial derivatives as $\partial_\cdot$ with respect to $(t, z, y, x)$. The stream function is defined as $\psi_i = \phi_{g;i}/f_0$ where $\phi_{g;i}$ is the geostrophic pressure anomaly and relative vorticity can be written as $\zeta_{g;i} = \nabla_h^2 \psi_i$. The layer-thickness equation on the other hand is [7]:

$$ 
\begin{aligned}
\partial_t h_i + u_{g;i} \partial_x h_i + v_{g;i} \partial_y h_i &= -H_i (\partial_x u_{a;i} + \partial_y v_{a;i}) \\
&= H_i \partial_z w_{a;i}
\end{aligned} 
\qquad (6) $$

We leave the derivation of the layered QGPV and its relation to the continuously stratified framework to Appendix A.

The ageostrophic vertical velocity can be diagnosed via the QG omega equation (Appendix B; [52,53]):

$$ N_i^2 \nabla_h^2 w_{a;i} + f_0^2 \partial_{zz} w_{a;i} = \beta \partial_x b_i - 2 \nabla_h \cdot Q_i - \nabla_h^2 b_i N_i^2 \partial_t \frac{1}{N_i^2}, \qquad (7) $$

where $N_i^2 = g_i'/H_i^\dagger$, and $b_i = f_0 \frac{\psi_i - \psi_{i+1}}{H_i^\dagger}$ is the buoyancy. The Q tensor is:

$$ Q_i = Q_i^1 i + Q_i^2 j = (\partial_x u_{g;i}^\dagger \cdot \nabla_h b_i) i + (\partial_y u_{g;i}^\dagger \cdot \nabla_h b_i) j, \qquad (8) $$

where $u_{g;i}^\dagger = -\partial_y \psi_i^\dagger i + \partial_x \psi_i^\dagger j$ is the geostrophic velocity derived from the inter-facial stream function ($\psi_i^\dagger = \frac{H_i \psi_{i+1} + H_{i+1} \psi_i}{H_i + H_{i+1}}$; [3]). i and j are the horizontal Cartesian unit vectors. The last term on the right-hand side of (7) is due to the temporally varying background stratification (Appendix B). We solved Equation (7) iteratively for $w_a$ via a two-dimensional geometric multigrid solver with the boundary conditions of Ekman pumping ($w_E$):

$$ w_{E;0} = -\frac{1}{f_0} \nabla_h \times \tau = -\frac{UH}{L} \hat{\tau}_0 \sin^2 \left[ \frac{2\pi y}{L_0} \right] \sin \left[ \frac{\pi y}{L_0} \right], \qquad (9) $$

$$ w_{E;3} = \frac{\delta_E}{2} \zeta_{g;3}, \qquad (10) $$

where $\delta_E = Ek^b H_3$ is the bottom Ekman-layer thickness.

Now, multiplying Equation (5) by $-\psi_i$ and integrating over the depth of each layer gives the kinetic energy (KE) budget:

$$ 
\begin{aligned}
H_i \left[ \frac{D_i}{Dt} \frac{|\nabla_h \psi_i|^2}{2} - \nabla_h \cdot (u_{g;i} \psi_i \nabla_h^2 \psi_i) - \frac{\beta}{2} \partial_x \psi_i^2 \right] \\
&= -f_0 \int \psi_i \partial_z w_{a;i} dz \\
&= f_0 \left[ - (w_{a;i-1} \psi_{i-1}^\dagger - w_{a;i} \psi_i^\dagger) + \int w_{a;i} \partial_z \psi_i dz \right] \\
&= f_0 \left[ - (w_{a;i-1} \psi_{i-1}^\dagger - w_{a;i} \psi_i^\dagger) + H_i \left( w_{a;i-1} \frac{\psi_{i-1} - \psi_i}{H_i + H_{i-1}} + w_{a;i} \frac{\psi_i - \psi_{i+1}}{H_{i+1} + H_i} \right) \right]. \quad (11)
\end{aligned}
$$

Dropping the divergence terms as they vanish upon area integration, for each layer, we get:

$$ \frac{H_1}{2} \partial_t |\nabla_h \psi_1|^2 = f_0 \left[ w_{a;1} \psi_1^\dagger + w_{a;1} H_1 \frac{\psi_1 - \psi_2}{H_2 + H_1} \right], \qquad (12) $$

$$ \frac{H_2}{2} \partial_t |\nabla_h \psi_2|^2 = f_0 \left[ - (w_{a;1} \psi_1^\dagger - w_{a;2} \psi_2^\dagger) + H_2 \left( w_{a;1} \frac{\psi_1 - \psi_2}{H_2 + H_1} + w_{a;2} \frac{\psi_2 - \psi_3}{H_3 + H_2} \right) \right], \qquad (13) $$
$$ \frac{H_3}{2} \partial_t |\nabla_h \psi_3|^2 = f_0 \left[ -w_{a;2} \psi_2^\dagger + w_{a;2} H_3 \frac{\psi_2 - \psi_3}{H_3 + H_2} \right]. \quad (14) $$

On the other hand, using relation (A2), the layer-thickness equations can be manipulated as $\frac{H_2}{H_1+H_2}(6)|_{i=1} - \frac{H_1}{H_1+H_2}(6)|_{i=2}$:

$$ \begin{aligned} \frac{D_1^\dagger}{Dt} \left[ \frac{f_0}{g_1'} (\psi_1 - \psi_2) \right] &= -w_{a;1} + \frac{H_1}{H_1+H_2} \left[ w_{a;2} - \frac{D_2}{Dt} \frac{f_0}{g_2'} (\psi_3 - \psi_2) \right] \\ &= -w_{a;1} + \frac{f_0 H_1}{g_2'(H_1+H_2)} (\mathbf{u}_3 - \mathbf{u}_2) \cdot \nabla_h (\psi_3 - \psi_2) \\ &= -w_{a;1}, \end{aligned} \quad (15) $$

where the second term on the right-hand side above (15) vanishes due to thermal wind. Similarly, $\frac{H_3}{H_2+H_3}(6)|_{i=2} - \frac{H_2}{H_2+H_3}(6)|_{i=3}$:

$$ \begin{aligned} \frac{D_2^\dagger}{Dt} \left[ \frac{f_0}{g_2'} (\psi_2 - \psi_3) \right] &= -w_{a;2} + \frac{H_3}{H_2+H_3} \left[ w_{a;1} - \frac{D_2}{Dt} \frac{f_0}{g_1'} (\psi_2 - \psi_1) \right] \\ &= -w_{a;2}, \end{aligned} \quad (16) $$

where $\frac{D_1^\dagger}{Dt} = \partial_t + \mathbf{u}_{g;1}^\dagger \cdot \nabla_h$. The available potential energy (APE) equations can, therefore, be derived by multiplying Equation (15) with $f_0(\psi_1 - \psi_2)$ and again dropping the divergence terms:

$$ \partial_t \left[ \frac{f_0^2}{2g_1'} (\psi_1 - \psi_2)^2 \right] = -f_0(\psi_1 - \psi_2)w_{a;1} - \frac{f_0^2 (\psi_1 - \psi_2)^2}{2} \partial_t g_1'^{-1}, \quad (17) $$

and Equation (16) with $f_0(\psi_2 - \psi_3)$:

$$ \partial_t \left[ \frac{f_0^2}{2g_2'} (\psi_2 - \psi_3)^2 \right] = -f_0(\psi_2 - \psi_3)w_{a;2}. \quad (18) $$

We see from Equation (17) that there is an additional source of APE due to the temporally varying background potential energy (BPE; $B^\#$), which then feeds back onto the KE via Equations (12) and (13) through baroclinic instability. BPE takes the same form as APE except that only $g'$ is inside the derivative.

Now, the mean KE (MKE; $K^\#$), eddy KE (EKE; $K^\circ$), mean APE (MAPE; $P^\#$) and eddy APE (EAPE; $P^\circ$) can be defined as:

$$ K_i^\# = \frac{H_i}{2} |\nabla_h \overline{\psi_i}|^2, \quad K_i = \frac{H_i}{2} |\overline{\nabla_h \psi_i'}|^2, \quad (19) $$

$$ P_i^\# = \frac{f_0^2}{2g_i'} (\overline{\psi_i} - \overline{\psi_{i+1}})^2, \quad P_i = \frac{f_0^2}{2g_i'} (\overline{\psi_i'} - \overline{\psi_{i+1}'})^2, \quad (20) $$

where $\overline{\cdot}$ is the ensemble mean and the eddy is defined as fluctuations about the ensemble mean, viz. $(\cdot)' = (\cdot) - \overline{\cdot}$. We note that the ensemble mean of the fluctuations vanish $(\overline{\cdot})' = 0$. The strength of defining the mean as such is that in addition to the ensemble-mean operator commuting with the derivatives with respect to $(t, z, y, x)$ [38], it provides a unique decomposition between the mean and eddy. In other words, the mean does not depend on an arbitrary temporal or spatial scale, which is beneficial in our case as the separated jet is on QG scaling in the cross-jet direction while on planetary-geostrophic scaling in the along-jet direction [54,55]. The ensemble mean can be interpreted as the QG response to external forcing while the eddies are a result of intrinsic variability [56–58]. The ensemble means of total KE and APE each satisfy $\overline{K_i} = \frac{H_i}{2} |\overline{\nabla_h \psi_i}|^2 = K_i^\# + K_i$, $\overline{P_i} =$
$\frac{f_0^2}{2g_i'}(\overline{\psi_i - \psi_{i+1}})^2 = P_i^\# + P_i$. Hence, the exchanges (Π) of KE and APE within and between layers are:

$$ \Pi_{K_1^\# \to K_1} = -H_1 \langle \overline{\psi_1} \nabla_h \cdot \overline{u_{g;1}'} \nabla_h^2 \psi_1' \rangle, \quad (21) $$

$$ \Pi_{K_1^\# \to K_2^\#} = -f_0 \langle \overline{w_{a;1}} \overline{\psi_1^\dagger} \rangle, \quad \Pi_{K_1 \to K_2} = -f_0 \langle \overline{w'_{a;1}} \overline{\psi_1^\dagger} \rangle, \quad (22) $$

$$ \Pi_{P_1^\# \to K_1^\#} = \frac{f_0 H_1}{H_2 + H_1} \langle \overline{w_{a;1}}(\overline{\psi_1} - \overline{\psi_2}) \rangle, \quad \Pi_{P_1 \to K_1} = \frac{f_0 H_1}{H_2 + H_1} \langle \overline{w'_{a;1}}(\psi_1' - \psi_2') \rangle, \quad (23) $$

$$ \Pi_{P_1^\# \to P_1} = \frac{f_0^2}{g_1'} \langle \overline{\psi_1} - \overline{\psi_2} \rangle \nabla_h \cdot \overline{u_{g;1}^\dagger} (\psi_1' - \psi_2') \rangle, \quad (24) $$

$$ \Pi_{B_1^\# \to P_1^\#} = -\frac{f_0^2}{2} \langle (\overline{\psi_1} - \overline{\psi_2})^2 \partial_t g_1'^{-1} \rangle, \quad \Pi_{B_1^\# \to P_1} = -\frac{f_0^2}{2} \langle (\overline{\psi_1} - \overline{\psi_2})^2 \partial_t g_1'^{-1} \rangle, \quad (25) $$

$$ \Pi_{K_2^\# \to K_2} = -H_2 \langle \overline{\psi_2} \nabla_h \cdot \overline{u_{g;2}'} \nabla_h^2 \psi_2' \rangle, \quad (26) $$

$$ \Pi_{K_2^\# \to K_3^\#} = -f_0 \langle \overline{w_{a;2}} \overline{\psi_2^\dagger} \rangle, \quad \Pi_{K_2 \to K_3} = -f_0 \langle \overline{w'_{a;2}} \overline{\psi_2^\dagger} \rangle, \quad (27) $$

$$ \Pi_{P_1^\# \to K_2^\#} = \frac{f_0 H_2}{H_2 + H_1} \langle \overline{w_{a;1}}(\overline{\psi_1} - \overline{\psi_2}) \rangle, \quad \Pi_{P_1 \to K_2} = \frac{f_0 H_2}{H_2 + H_1} \langle \overline{w'_{a;1}}(\psi_1' - \psi_2') \rangle, \quad (28) $$

$$ \Pi_{P_2^\# \to K_2^\#} = \frac{f_0 H_2}{H_3 + H_2} \langle \overline{w_{a;2}}(\overline{\psi_2} - \overline{\psi_3}) \rangle, \quad \Pi_{P_2 \to K_2} = \frac{f_0 H_2}{H_3 + H_2} \langle \overline{w'_{a;2}}(\psi_2' - \psi_3') \rangle, \quad (29) $$

$$ \Pi_{P_2^\# \to P_2} = \frac{f_0^2}{g_2'} \langle (\overline{\psi_2} - \overline{\psi_3}) \nabla_h \cdot \overline{u_{g;2}^\dagger} (\psi_2' - \psi_3') \rangle, \quad (30) $$

$$ \Pi_{K_3^\# \to K_3} = -H_3 \langle \overline{\psi_3} \nabla_h \cdot \overline{u_{g;3}'} \nabla_h^2 \psi_3' \rangle, \quad (31) $$

$$ \Pi_{P_2^\# \to K_3^\#} = \frac{f_0 H_3}{H_2 + H_3} \langle \overline{w_{a;2}}(\overline{\psi_2} - \overline{\psi_3}) \rangle, \quad \Pi_{P_2 \to K_3} = \frac{f_0 H_3}{H_3 + H_2} \langle \overline{w'_{a;2}}(\psi_2' - \psi_3') \rangle, \quad (32) $$

where $\langle \cdot \rangle = \iint (\cdot) dxdy$ is the area integration. Further details regarding the sign convention and forcing/dissipation terms are given in Appendices C and D. Summing up each layer gives the volume integrated energy exchanges:

$$ \Pi_{P^\# \to K^\#} = \sum_{i=1}^{2} f_0 \langle \overline{w_{a;1}}(\overline{\psi_i} - \overline{\psi_{i+1}}) \rangle, \quad (33) $$

$$ \Pi_{P \to K} = \sum_{i=1}^{2} f_0 \langle \overline{w'_{a;1}}(\psi_i' - \psi_{i+1}') \rangle, \quad (34) $$

$$ \Pi_{P^\# \to P} = \sum_{i=1}^{2} \frac{f_0^2}{g_i'} \langle (\overline{\psi_i} - \overline{\psi_{i+1}}) \nabla_h \cdot \overline{u_{g;1}^\dagger} (\psi_i' - \psi_{i+1}') \rangle, \quad (35) $$

$$ \Pi_{K^\# \to K} = - \sum_{i=1}^{3} H_i \langle \overline{\psi_i} \nabla_h \cdot \overline{u_{g;1}'} \nabla_h^2 \psi_i' \rangle, \quad (36) $$

$$ \Pi_{B^\# \to P^\#} = -\frac{f_0^2}{2} \langle (\overline{\psi_1} - \overline{\psi_2})^2 \partial_t g_1'^{-1} \rangle, \quad (37) $$

$$ \Pi_{B^\# \to P} = -\frac{f_0^2}{2} \langle (\overline{\psi_1'} - \overline{\psi_2'})^2 \partial_t g_1'^{-1} \rangle. \quad (38) $$

## 4. Results

We start by showing the total kinetic energy (TKE) during the spin-up phase and for the 10 years of output we have (viz. 20 years in total; Figure 2). The ensemble spread starts to grow after 1.5 years of integration from the perturbed initial conditions and plateaus
roughly for the latter seven years. The area-integrated TKE in the first layer (⟨⟨$K_1$⟩⟩), most relevant for studies interested in surface seasonal dynamics, is in sync with the background stratification ($g'_1$), viz. higher ⟨$K_1$⟩ during summer when stratification is stronger and visa versa (Figure 2b). For the lower layers, there is a temporal lag evident by the barotropic TKE (⟨|$\nabla_h$Ψ|²⟩ where $Ψ = H^{-1}\sum_i H_i \psi_i$ is the barotropic stream function; Figure 2a). Although it is difficult to detect a clear seasonal signal for the barotropic TKE from an individual ensemble member such as in the CTRL run, their ensemble mean shows a robust seasonality. For the remainder of the study, we use the last five years of output in order to maximize the signal of intrinsic variability amongst members.

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Time series of the horizontally averaged barotropic (a) and first-layer TKE (b). Each ensemble member is shown with thin gray lines and standard deviation of the ensemble mean in black shading. The CTRL run is shown with a red dashed line and ensemble mean with a cyan dot-dashed line. The reduced gravity is shown in blue plotted against the right y axis.</figcaption></figure>

In Figure 3, we show the mean and eddy KE in the first layer ($K_1^\#$, $K_1$) during summer and winter for the last year of output and their difference. The seasons were defined at the time step when the reduced gravity was at its maximum and minimum, respectively. We see the characteristic feature of a robust separated Western Boundary Current in a double-gyre system with very little meandering while the EKE is more meridionally spread out. Consistent with Figure 2, summertime has a stronger mean jet and EKE than winter (Figure 3e,f). We also show snapshots of eddy PV ($q_{g';1}'$) from the CTRL run from which we see coherent features of mesoscale eddies (Figure 3g,h).
<figure><img src="image_6.png" /><figcaption>**Figure 3.** The summer and wintertime mean and eddy KE and their difference during the last year of output (a–f). Note the differences are plotted on a logarithmic scale. (g,h) Snapshot of eddy PV for summer and winter during the last year of output from the CTRL run. All panels show the variable in the first layer.</figcaption></figure>

4.1. *The Domain Integrated Lorenz Energy Cycle*

We now move on to quantifying the LEC in order to examine the processes responsible
for generating higher KE during summertime. As we define the mean as the ensemble
mean (as opposed to a temporal mean which has commonly been applied), we are able
to examine the temporal variability of LEC. We compute the terms in Equations (33)–(38)
for the last five years of output and show them in Figure 4. The time series of MAPE is in
sync with the background stratification dominated by *g*' in its denominator while MKE
lags $g'_1$ by ~11 days upon taking the lag correlation ($P^\#, K^\#$; Figure 4a). MAPE has the largest magnitude amongst the reservoirs by an order of magnitude and for KE, the eddies are more energetic than the mean. The energy flux from MAPE to MKE is negative year round ($\Pi_{P^\# \to K^\#} < 0$; black solid line in Figure 4b) due to Ekman pumping steepening the isopycnals. The energy input due to wind stress ($F_s^{K^\#}$) is in sync with MKE with energetic surface currents resulting in a stronger surface stress. The eddy energy reservoirs ($P, K$), on the other hand, lag the stratification by ~17 days but their peaks precede winter when the domain is most susceptible to baroclinic instability and energy conversion from EAPE to EKE takes its yearly maximum ($\Pi_{P \to K}$; Figure 4a). It is perhaps interesting to note that the sign of flux between EAPE and EKE occasionally reverses during summer with barotropic instability over compensating for baroclinic instability; the energy pathway becomes MKE→EKE→EAPE ($\Pi_{P \to K} < 0$), whereas baroclinic instability would predict EAPE→EKE ($\Pi_{P \to K} > 0$). Regarding the dissipation terms, only the bottom drag for EKE ($D_b^K$) shows a notable seasonality and has a similar magnitude to the energy flux from MKE to EKE ($\Pi_{K^\# \to K}$). The amplitude of bottom drag ($|D_b^K|$) lags EKE by ~41 days and aligns well with the ensemble-mean barotropic TKE (Figure 2a).

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Time series of each term in the domain-integrated LEC. (a) The mean and eddy KE (black) and APE (red) reservoirs in the units of ×10¹³ [(J/kg) m³] and stratification of the first layer interface (blue; $g'_1$). MAPE is multiplied by 0.1 to have it fit on the same y axis. (b) The energy fluxes between each energy reservoir and forcing terms due to surface wind stress ($F_s^{K^\#}$) and temporally varying BPE. (c) Dissipation terms due to horizontal viscosity ($D_h$) and bottom friction ($D_b$). The mean and eddy horizontal dissipation terms are lumped together and fluxes are in the units of ×10⁶ [(W/kg) m³]. The forcing and dissipation terms are detailed in Appendix D.</figcaption></figure>
To provide a climatological view of the energy fluxes (II), we take the yearly average of the last five years and show the LEC diagram for a climatological summer and winter (Figure 5). Each season per year is defined as four time steps; summer is when the reduced gravity takes its maxima and four time steps about it, and four time steps about the minima in reduced gravity for winter. The seasonal climatology is then taken as the average of the five years. Again, we see that all reservoirs are more energetic during the summer. Focusing on MKE, except for the surface wind stress, the reservoir has loss terms year round and yet stores more energy during the summer. We attribute the summertime maxima in MKE to the separated jet stabilizing due to increased stratification, which results in the jet shedding stronger eddies. Indeed, the energy flux from MKE to EKE ($\Pi_{K^{#} \to K}$) is highest during the summer (Figures 4b and 5a). We attribute the larger energy conversion from EAPE to EKE during the winter ($\Pi_{P \to K}$) to the flow being more susceptible to baroclinic instability with reduced stratification.

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Time series of the seasonal climatology of energy fluxes between the energy reservoirs (a). (b,c) The LEC diagram for the climatological summer and winter averaged over the last five years of output. The energies are in the units of ×10¹³ [(J/kg) m³] and fluxes are in ×10⁶ [(W/kg) m³]. The energy exchanges do not exactly cancel out due to each reservoir having temporal variability.</figcaption></figure>

## 4.2. Time Lag in Lower-Layer Energetics

In this section, we investigate the mechanism for the lag in KE in the lower layers ($K_2$, $K_3$) from KE in the first layer ($K_1$) and stratification ($g'_1$) implied from Figure 2. It is perhaps interesting to note that although the ensemble-mean barotropic TKE lags $g'_1$ by ~41 days, neither the domain-integrated MKE nor EKE show such a long lag (Figure 4a). This has to do with MKE and EKE being volume-weighted variables of quadratic terms, while the barotropic TKE being a quadratic term of a volume-weighted variable; MKE and EKE have a larger weighting on the surface stream function, which is in sync with $g'_1$, than the barotropic TKE. The lag within lower layers becomes apparent for the time series of area integrated EKE within each layer (~128 days for $K_2$ and ~68 days for $K_3$; Figure 6a). We also focus on EKE for the remainder of this section as EKE is always larger than MKE
by a factor of three (Figure 4). Examining the energy fluxes, Figure 6b shows that the
contribution from barotropic instability becomes negligible within the lower two layers
with the relative significance of the separated jet diminishing with depth, and shows
no clear seasonality ($\Pi_{K_{2,3}^{#} \to K_{2,3}}$). The vertical transfer of EKE ($\Pi_{K_1 \to K_2}$, $\Pi_{K_2 \to K_3}$) and
conversion from EAPE ($\Pi_{P \to K}$), on the other hand, show a coherent seasonal pattern with
the maxima of $K_2$ and $K_3$ falling in between the maxima of the two fluxes. We, therefore,
attribute the time lag in $K_2$ and $K_3$ to the balance between baroclinic instability and vertical
transfers of EKE.

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Time series of volume-integrated EKE over each layer, and fluxes within and between layers ($\Pi_{K_1 \to K_2}$, $\Pi_{K_2 \to K_3}$) plotted along with the reduced gravity ($g_1'$). (a) The EKEs have their temporal mean removed so as to plot against the same y axis. (b) A rolling mean by five time steps (~29 days) is applied to the time series of the energy fluxes. The energy flux from MKE to EKE within the two bottom layers is summed up ($\Pi_{K_{2,3}^{#} \to K_{2,3}}$) and conversion from APE is shown as the conversion rate volume integrated over the three layers as the amount that goes into each layer is simply the total conversion weighted by layer thickness (cf. Equations (23), (28), (29) and (32)). The conversion from $P_1$ and $P_2$ were in sync with each other (not shown). For further details regarding each term, see Appendix D.</figcaption></figure>

**5. Discussion and Conclusions**

By running a seasonally forced 101-member ensemble of a three-layer quasi-geostrophic (QG) model in an idealized double-gyre configuration, we have shown that the kinetic energy (KE) peaks during summer when the (basin-scale) stratification is strongest during the year (Figure 2). Such seasonality in mesoscale eddy KE (EKE) has been observed in other studies using realistic simulations of the ocean [14,27–31]. Due to the air–sea interaction, the seasonal modulation of the mixed-layer depth leads to a strong seasonal signal in submesoscale instabilities. The submesoscale EKE takes its maximum during late winter/early spring and previous studies have commonly explained the summer peak in the mesoscale range as the time lag for the submesoscale EKE to cascade upscale. The mechanism of inverse energy cascade fails, however, to explain the mesoscale seasonality in our model, as a QG model by definition cannot resolve any submesoscale instabilities.

Using the framework of the Lorenz energy cycle (LEC; [2]), we have quantified the
reservoirs of mean and eddy available potential energy (APE) and KE, and energy fluxes
amongst them. We note that our ensemble framework has allowed us to examine the
seasonal variability of LEC. Our results show that all four reservoirs store more energy during the summer than winter (Figure 4a). For the mean KE (MKE), we attribute the summertime maximum to increased stratification leading to a more baroclinically stable and stronger jet. Conceptually, this can be understood based on a mass-flux balance. Since the wind stress is kept stationary, the Sverdrup transport ($\beta^{-1}\nabla_h \times \tau$) remains constant throughout the simulation. Based on mass balance, the accumulating transport towards the north/south boundaries must be fluxed out via the Western Boundary Current. Figure 3c shows an intensification of MKE during summer along the Western Boundary resulting from less energy lost to the eddies within the gyre interior. Hence, a more stable jet results in a stronger mean flow. Our results of jet stabilization and its zonal penetration into the gyre are complementary to earlier studies where they attributed the penetration scale to parameters of lateral friction, vertical resolution and topography [4,59]. Here, we have investigated the effect of a seasonally varying background stratification.

Shifting our focus to EKE, based on baroclinic instability, one might expect the opposite to be true, namely, wintertime having more EKE than summertime due to weaker stratification. The LEC shows that year round, energy fluxes from MKE to EKE associated with barotropic instability over compensate for the fluxes from eddy APE (EAPE) to EKE, a pathway associated with baroclinic instability. Since MKE is higher during summer, the larger flux of energy from MKE to EKE results in EKE peaking in summer (Figures 4b and 5). Although our simulation is highly idealized, we argue that barotropic processes dominating in the separated jet region are consistent with a recent study on energetics using a realistic simulation of the North Atlantic Ocean [55]. We note that the balance between barotropic and baroclinic instability in our LEC is in the domain integrated sense. In a domain without a jet, we would expect baroclinic instability to be the dominant mechanism in generating eddies so long as the background state is baroclinically unstable.

To our knowledge, Qiu et al. [26] is the only study using a realistic ocean simulation showing how the seasonality in background state can modulate the mesoscale variability. Their results differ from ours, however, in that they attribute the mesoscale seasonality to the classical Phillips-like baroclinic instability arising from the interior background stratification and vertical shear in horizontal velocity [1]. In addition to the submesoscale variability modulating mesoscale seasonality, our results suggest that, in reality, it is possible that the basin-scale variability does so as well. We note that since our QG model does not permit submesoscales, the baroclinic energy flux from EAPE to EKE is likely underestimated compared to the real ocean. It would be interesting to revisit the LEC for realistic ocean ensembles [57,58] to see whether we would see a stabilization of the separated Gulf Stream during summer and consequently larger energy fluxes from MKE to EKE.

**Author Contributions:** Conceptualization, T.U.; methodology, T.U. and B.D.; software, B.D.; validation, T.U., B.D.; formal analysis, T.U.; investigation, T.U. and B.D.; computational resources, T.P.; data curation, T.U.; writing, T.U.; visualization, T.U.; project administration, T.P.; funding acquisition, T.P. All authors have read and agreed to the published version of the manuscript.

**Funding:** 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. We acknowledge high-performance computing resources for generating the ensemble and analyzing our model outputs on Occigen maintained by CINES with the reference A0090112020.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The last time step from the low-resolution ensemble stream function and example analysis code in Python is available on Github (doi:10.5281/zenodo.4667532). For more output files, please contact the corresponding author.
**Acknowledgments:** We thank William K. Dewar for his always insightful and fruitful discussion. The authors would like to acknowledge the developers of the Basilisk language (http://basilisk.fr/, accessed on 9 April 2021) upon which MSOM is based, and the developers of the xarray Python package [60], which we used to post process the model outputs.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

## Appendix A. Derivation of the Layered Quasi-Geostrophic Potential Vorticity

As the relative vorticity Equation (5) and layer-thickness Equation (6) have a common term on the right-hand side, they can be combined as:

$$ \frac{D_i}{Dt} \zeta_{g; i} + \beta v_{g; i} = \frac{f_0}{H_i} \frac{D_i}{Dt} h_i, \qquad (A1) $$

and we get the governing equation for QGPV $q_i = \zeta_{g; i} + \beta y - \frac{f_0}{H_i} h_i$ [7]. It is perhaps interesting to note that the QGPV remains identical for a stationary and temporally varying background stratification (viz. $g'_1 = \frac{U_1^2}{H_1^2} Fr_1^{-2}(t)$) although we have shown that this is not the case for the energy budget. The stream function is related to the layer displacement via $\eta_i = \frac{f_0}{g'_i}(\psi_{i+1} - \psi_i)$. The layer thickness can, therefore, be written using the stream function as [7]:

$$ \begin{aligned} h_i &= H_i + \eta_{i-1} - \eta_i \\ &= H_i + \frac{f_0}{g'_{i-1}}(\psi_i - \psi_{i-1}) - \frac{f_0}{g'_i}(\psi_{i+1} - \psi_i), \end{aligned} \qquad (A2) $$

where $\frac{D_1}{Dt} \eta_0 = \frac{D_3}{Dt} \eta_3 = 0$ due to rigid-lid and flat-bottom boundary conditions.

Now, suppose at any given time, we have total buoyancy (B) defined on a layer interface (Figure A1). Based on Taylor expansion, the layer interface displacement can be expanded as [7]:

$$ \begin{aligned} \eta &= \left. \frac{\partial z}{\partial B} \right|_{z=H} [B_0(t, z = H) - B(t, z = H + \eta)] \\ &= - \left. \frac{\partial z}{\partial B} \right|_{z=H} b, \end{aligned} \qquad (A3) $$

where $b = B - B_0$ is the QG fluctuation about the background buoyancy ($B_0$). Hence, we get:

$$ \frac{b}{N^2} = -\eta, \qquad (A4) $$

and taking the material derivative gives the buoyancy equation in the continuously stratified framework:

$$ \frac{D}{Dt} \frac{b}{N^2} = -w. \qquad (A5) $$

Equation (A4) gives the physical intuition that the material derivative of $b/N^2$ leads to vortex stretching.

<figure><img src="image_10.png" /><figcaption>**Figure A1.** Schematic of a relation between buoyancy (B) and layer interface displacement ($\eta$). The background buoyancy is $B_0$ defined at $z = H$.</figcaption></figure>
## Appendix B. The Omega Equation with a Temporally Varying Background Stratification

We derive the QG omega equation using the continuously stratified framework. Taking the vertical derivative of the order-Rossby number momentum equations with the viscous term:

$$ \partial_t u_{g;i} + u_{g;i} \partial_x u_{g;i} + v_{g;i} \partial_y u_{g;i} - f_0 v_{a;i} - \beta y v_{g;i} = - \partial_x \phi_{a;i}, \quad (A6) $$

$$ \partial_t v_{g;i} + u_{g;i} \partial_x v_{g;i} + v_{g;i} \partial_y v_{g;i} + f_0 u_{a;i} + \beta y u_{g;i} = - \partial_y \phi_{a;i}, \quad (A7) $$

and multiplying them by $f_0$ gives:

$$ \frac{D}{Dt} (f_0 \partial_z u_g) + \partial_y u_g \cdot \nabla_h b - f_0^2 \partial_z v_a - \beta y f_0 \partial_z v_g = -f_0 v_4 \partial_z \nabla_h^4 u_g, \quad (A8) $$

$$ \frac{D}{Dt} (f_0 \partial_z v_g) - \partial_x u_g \cdot \nabla_h b + f_0^2 \partial_z u_a + \beta y f_0 \partial_z u_g = -f_0 v_4 \partial_z \nabla_h^4 v_g, \quad (A9) $$

and the horizontal gradients of the buoyancy Equation (A5) with the diffusive term yields:

$$ \frac{1}{N^2} \frac{D}{Dt} \partial_x b + \partial_x b \partial_t \frac{1}{N^2} + \frac{\partial_x u_g}{N^2} \cdot \nabla_h b + \partial_x w_a = - \kappa_4 \partial_x \nabla_h^4 b, \quad (A10) $$

$$ \frac{1}{N^2} \frac{D}{Dt} \partial_y b + \partial_y b \partial_t \frac{1}{N^2} + \frac{\partial_y u_g}{N^2} \cdot \nabla_h b + \partial_y w_a = - \kappa_4 \partial_y \nabla_h^4 b. \quad (A11) $$

Summing Equation (A8) with (A11), and -(A9) with (A10), and using the thermal wind relation, we get:

$$ 2 \partial_y u_g \cdot \nabla_h b + N^2 \partial_y w_a - \beta y \partial_x b - f_0^2 \partial_z v_a + \partial_y b N^2 \partial_t \frac{1}{N^2} = 0. \quad (A12) $$

$$ 2 \partial_x u_g \cdot \nabla_h b + N^2 \partial_x w_a + \beta y \partial_y b - f_0^2 \partial_z u_a + \partial_x b N^2 \partial_t \frac{1}{N^2} = 0. \quad (A13) $$

The viscous and diffusive terms do not appear as they cancel out due to the thermal-wind relation ($f_0 \partial_z \zeta_g = \nabla_h^2 b$) and their parameters being set identical (viz. $v_4 = \kappa_4 (= Re_4^{-1} L^3 U)$). Taking $\partial_y$ (A12) + $\partial_x$ (A13) gives the omega equation for a temporally varying background stratification:

$$ N^2 \nabla_h^2 w_a + f_0^2 \partial_{zz} w_a = \beta \partial_x b - 2 \nabla_h \cdot Q - \nabla_h^2 b N^2 \partial_t \frac{1}{N^2}. \quad (A14) $$

Although the last term on the right-hand side involves a time derivative, there is no time dependency in our case as we know the analytical form of the background stratification (Equation (4)). Its contribution to the omega equation turned out to be negligible (not shown).

## Appendix C. Decomposing the Mean and Eddy Energetics

In this section, we derive the mean and eddy KE equations. Equation (5) can be split into its mean and eddy component:

$$ \frac{D^\#}{Dt} \nabla_h^2 \bar{\psi} + \frac{D^\#}{Dt} \nabla_h^2 \psi' + u'_g \cdot \nabla_h [\nabla_h^2 (\bar{\psi} + \psi')] + \beta \partial_x (\bar{\psi} + \psi') = f_0 \partial_z (\bar{w} + w'), \quad (A15) $$
where $\frac{D^#}{Dt} = \partial_t + \overline{u_g} \cdot \nabla_h$. Multiplying this by $-\bar{\psi}$ gives:

$$ \frac{D^#}{Dt} \frac{|\nabla_h \bar{\psi}|^2}{2} - \nabla_h \cdot \overline{u_g \bar{\psi} \nabla_h^2 \bar{\psi}} - \bar{\psi} \frac{D^#}{Dt} \nabla_h^2 \psi' - \bar{\psi} \overline{u_g' \cdot \nabla_h [\nabla_h^2 (\bar{\psi} + \psi')]} \\ - \beta \partial_x \frac{\bar{\psi}^2}{2} - \bar{\psi} \beta \partial_x \psi' = \overline{w} \bar{b} - \bar{\psi} f_0 \partial_z w', \quad (A16) $$

and taking its ensemble mean yields the mean KE equation:

$$ \frac{D^#}{Dt} \frac{|\nabla_h \bar{\psi}|^2}{2} - \nabla_h \cdot \overline{u_g \bar{\psi} \nabla_h^2 \bar{\psi}} - \beta \partial_x \frac{\bar{\psi}^2}{2} = \overline{w} \bar{b} + \bar{\psi} \nabla_h \cdot \overline{u_g' \nabla_h^2 \psi'}. \quad (A17) $$

On the other hand, the ensemble mean of the total KE Equation (11) is:

$$ \frac{D}{Dt} \frac{|\nabla_h \psi|^2}{2} - \nabla_h \cdot \overline{u_g \psi \nabla_h^2 \psi} - \beta \partial_x \frac{\bar{\psi}^2}{2} = \overline{wb}, \quad (A18) $$

which can be expanded as:

$$ \frac{D^#}{Dt} \frac{|\nabla_h \bar{\psi}|^2}{2} + \frac{D^#}{Dt} \frac{|\nabla_h \psi'|^2}{2} + \nabla_h \cdot \overline{u_g' \frac{|\nabla_h \psi'|^2}{2}} + \overline{u_g' \cdot \nabla_h (\partial_x \bar{\psi} \partial_x \psi' + \partial_y \bar{\psi} \partial_y \psi')} \\ - \nabla_h \cdot \overline{u_g \psi \nabla_h^2 \psi} - \beta \partial_x \frac{\psi^2}{2} = \overline{wb}. \quad (A19) $$

Taking the difference between Equations (A17) and (A19) gives the eddy KE equation:

$$ \frac{D^#}{Dt} \frac{|\nabla_h \psi'|^2}{2} + \nabla_h \cdot \overline{u_g' \frac{|\nabla_h \psi'|^2}{2}} + \nabla_h \cdot \overline{u_g' (\partial_x \bar{\psi} \partial_x \psi' + \partial_y \bar{\psi} \partial_y \psi')} \\ - \nabla_h \cdot (\overline{u_g \psi \nabla_h^2 \psi} - \overline{u_g \bar{\psi} \nabla_h^2 \psi}) - \beta \partial_x \frac{\psi'^2}{2} = \overline{w' b'} - \bar{\psi} \nabla_h \cdot \overline{u_g' \nabla_h^2 \psi'}. \quad (A20) $$

Since the divergence terms vanish upon area integration, we can see the mean and eddy KE exchanging the term $-\bar{\psi} \nabla_h \cdot \overline{u_g' \nabla_h^2 \psi'}$ (Equations (A17) and (A20)). The same procedure can be done for Equation (6) or the buoyancy equation to derive the mean and eddy APE equations.

## Appendix D. The Three-Layer QG Lorenz Energy Cycle

The Lorenz energy cycle [2] for the first layer, dropping the divergence terms in Equations (A17) and (A20), while bringing back the viscous and diffusive terms becomes:

$$ \partial_t K_1^# = f_0 \left[ \overline{w_{a;1} \psi_1^+} + H_1 \overline{w_{a;1} \frac{\psi_1^- - \psi_2^-}{H_2 + H_1}} \right] + H_1 \overline{\psi_1^-} \nabla_h \cdot \overline{u_{g;1}' \nabla_h^2 \psi_1'} \\ - \overline{\psi_1^-} \nabla_h \times \tau + H_1 \overline{\psi_1^-} \nu_4 \nabla_h^4 (\nabla_h^2 \psi_1^-), \quad (A21) $$

$$ \partial_t K_1 = f_0 \left[ w'_{a;1} \psi_1^{+\prime} + H_1 w'_{a;1} \frac{\psi_1' - \psi_2'}{H_2 + H_1} \right] - H_1 \overline{\psi_1^-} \nabla_h \cdot \overline{u_{g;1}' \nabla_h^2 \psi_1'} + H_1 \psi_1' \nu_4 \nabla_h^4 (\nabla_h^2 \psi_1'), \quad (A22) $$

$$ \partial_t P_1^# = -f_0 \overline{w_{a;1}} (\overline{\psi_1} - \overline{\psi_2}) - \frac{f_0^2}{g_1'} (\overline{\psi_1} - \overline{\psi_2}) \nabla_h \cdot \overline{u_{g;1}^+ (\psi_1' - \psi_2')} \\ - \frac{f_0^2}{g_1'} (\overline{\psi_1} - \overline{\psi_2}) \kappa_4 \nabla_h^4 (\overline{\psi_1} - \overline{\psi_2}) - \frac{f_0^2 (\overline{\psi_1} - \overline{\psi_2})^2}{2} \partial_t g_1'^{-1}, \quad (A23) $$
$$ \partial_t \mathcal{P}_1 = -f_0 \overline{w'_{a;1}(\psi'_1 - \psi'_2)} + \frac{f_0^2}{g'_1} (\overline{\psi_1} - \overline{\psi_2}) \nabla_h \cdot \overline{u'_{g;1}(\psi'_1 - \psi'_2)} \\ - \frac{f_0^2}{g'_1} (\overline{\psi'_1 - \psi'_2}) \kappa_4 \nabla_h^4 (\psi'_1 - \psi'_2) - \frac{f_0^2 (\psi'_1 - \psi'_2)^2}{2} \partial_t g'_1^{-1}. \quad (A24) $$

For the second layer:

$$ \partial_t K_2^\# = f_0 \left[ - (\overline{w_{a;1}} \overline{\psi_1^\dagger} - \overline{w_{a;2}} \overline{\psi_2^\dagger}) + H_2 (\overline{w_{a;1}} \frac{\overline{\psi_1} - \overline{\psi_2}}{H_2 + H_1} + \overline{w_{a;2}} \frac{\overline{\psi_2} - \overline{\psi_3}}{H_3 + H_2}) \right] \\ + H_2 \overline{\psi_2} \nabla_h \cdot \overline{u'_{g;2}} \nabla_h^2 \psi'_2 + H_2 \overline{\psi_2} \kappa_4 \nabla_h^4 (\nabla_h^2 \overline{\psi_2}), \quad (A25) $$

$$ \partial_t K_2 = f_0 \left[ - (\overline{w'_{a;1}} \overline{\psi_1^\dagger} - \overline{w'_{a;2}} \overline{\psi_2^\dagger}) + H_2 (\overline{w'_{a;1}} \frac{\overline{\psi'_1} - \overline{\psi'_2}}{H_2 + H_1} + \overline{w'_{a;2}} \frac{\overline{\psi'_2} - \overline{\psi'_3}}{H_3 + H_2}) \right] \\ - H_2 \overline{\psi_2} \nabla_h \cdot \overline{u'_{g;2}} \nabla_h^2 \psi'_2 + H_2 \overline{\psi'_2} \kappa_4 \nabla_h^4 (\nabla_h^2 \overline{\psi'_2}), \quad (A26) $$

$$ \partial_t P_2^\# = -f_0 \overline{w_{a;2}} (\overline{\psi_2} - \overline{\psi_3}) - \frac{f_0^2}{g'_2} (\overline{\psi_2} - \overline{\psi_3}) \nabla_h \cdot \overline{u'_{g;2}(\psi'_2 - \psi'_3)} - \frac{f_0^2}{g'_2} (\overline{\psi_2} - \overline{\psi_3}) \kappa_4 \nabla_h^4 (\overline{\psi_2} - \overline{\psi_3}), \quad (A27) $$

$$ \partial_t P_2 = -f_0 \overline{w'_{a;2}(\psi'_2 - \psi'_3)} + \frac{f_0^2}{g'_2} (\overline{\psi_2} - \overline{\psi_3}) \nabla_h \cdot \overline{u'_{g;2}(\psi'_2 - \psi'_3)} - \frac{f_0^2}{g'_2} (\overline{\psi'_2} - \overline{\psi'_3}) \kappa_4 \nabla_h^4 (\overline{\psi'_2} - \overline{\psi'_3}). \quad (A28) $$

For the third layer:

$$ \partial_t K_3^\# = f_0 \left[ - \overline{w_{a;2}} \overline{\psi_2^\dagger} + H_3 \overline{w_{a;2}} \frac{\overline{\psi_2} - \overline{\psi_3}}{H_3 + H_2} \right] + H_3 \overline{\psi_3} \nabla_h \cdot \overline{u'_{g;3}} \nabla_h^2 \psi'_3 \\ + H_3 \overline{\psi_3} [\kappa_4 \nabla_h^4 (\nabla_h^2 \overline{\psi_3}) + \epsilon \nabla_h^2 \overline{\psi_3}], \quad (A29) $$

$$ \partial_t K_3 = f_0 \left[ - \overline{w'_{a;2}} \overline{\psi_2^\dagger} + H_3 \overline{w'_{a;2}} \frac{\overline{\psi'_2} - \overline{\psi'_3}}{H_3 + H_2} \right] - H_3 \overline{\psi_3} \nabla_h \cdot \overline{u'_{g;3}} \nabla_h^2 \psi'_3 \\ + H_3 \overline{\psi'_3} [\kappa_4 \nabla_h^4 (\nabla_h^2 \psi'_3) + \epsilon \nabla_h^2 \psi'_3]. \quad (A30) $$

Although the biharmonic diffusive terms in the APE Equations (A23), (A24), (A27) and (A28), which originate from diffusive terms in the layer-thickness Equation (6), are applied solely for numerical stability and their similarity with buoyancy in primitive equations, their formulation is conceptually similar to the Gent-McWilliams' skew diffusivity (GM; [61]). GM represents the process of baroclinic instability upon which isopycnal displacements are smoothed out adiabatically within the isopycnal layer. Considering the quasi two-dimensional and adiabatic nature of the QG system, the interpretation of layer-thickness diffusivity becomes similar to the GM skew diffusivity. A major difference here is that the diffusivity is set as the bihamonic diffusivity and as such, should be negligible in damping the resolved eddies [14,45].

## References

1. Phillips, N.A. The general circulation of the atmosphere: A numerical experiment. Q. J. R. Meteorol. Soc. **1956**, *82*, 123–164. [CrossRef]

2. Lorenz, E.N. *The Nature and Theory of the General Circulation of the Atmosphere*; World Meteorological Organization: Geneva, Switzerland, 1967; Volume 218.

3. Holland, W.R. The role of mesoscale eddies in the general circulation of the ocean—Numerical experiments using a wind-driven quasi-geostrophic model. *J. Phys. Oceanogr.* **1978**, *8*, 363–392. [CrossRef]

4. Holland, W.R.; Schmitz, W.J. Zonal penetration scale of model midlatitude jets. *J. Phys. Oceanogr.* **1985**, *15*, 1859–1875. [CrossRef]

5. Kang, D.; Curchitser, E.N. Energetics of Eddy–Mean Flow Interactions in the Gulf Stream Region. *J. Phys. Oceanogr.* **2015**, *45*, 1103–1120. [CrossRef]

6. McWilliams, J.C. A perspective on the legacy of Edward Lorenz. *Earth Space Sci.* **2019**, *6*, 336–350. [CrossRef]

7. Vallis, G.K. *Atmospheric and Oceanic Fluid Dynamics*, 2nd ed.; Cambridge University Press: Cambridge, UK, 2017.
8. Su, Z.; Wang, J.; Klein, P.; Thompson, A.F.; Menemenlis, D. Ocean submesoscales as a key component of the global heat budget. *Nat. Commun.* **2018**, *9*. [CrossRef]

9. Ajayi, A.; Le Sommer, J.; Chassignet, E.; Molines, J.M.; Xu, X.; Albert, A.; Cosme, E. Spatial and Temporal Variability of the North Atlantic Eddy Field From Two Kilometric-Resolution Ocean Models. *J. Geophys. Res. Ocean.* **2020**, *125*, e2019JC015827. [CrossRef]

10. Lévy, M.; Franks, P.J.; Smith, K.S. The role of submesoscale currents in structuring marine ecosystems. *Nat. Commun.* **2018**, *9*, 1–16. [CrossRef]

11. Thomas, L.N.; Tandon, A.; Mahadevan, A. Submesoscale processes and dynamics. *Ocean Model. Eddying Regime* **2008**, *177*, 17–38. [CrossRef]

12. McWilliams, J.C. Submesoscale currents in the ocean. *Proc. R. Soc. A Math. Phys. Eng. Sci.* **2016**, *472*, 20160117. [CrossRef] [PubMed]

13. Buckingham, C.E.; Naveira Garabato, A.C.; Thompson, A.F.; Brannigan, L.; Lazar, A.; Marshall, D.P.; George Nurser, A.; Damerell, G.; Heywood, K.J.; Belcher, S.E. Seasonality of submesoscale flows in the ocean surface boundary layer. *Geophys. Res. Lett.* **2016**, *43*, 2118–2126. [CrossRef]

14. Uchida, T.; Balwada, D.; Abernathey, R.; McKinley, G.; Smith, S.; Levy, M. The contribution of submesoscale over mesoscale eddy iron transport in the open Southern Ocean. *J. Adv. Model. Earth Syst.* **2019**, *11*, 3934–3958. [CrossRef]

15. Siegelman, L.; Klein, P.; Rivière, P.; Thompson, A.F.; Torres, H.S.; Flexas, M.; Menemenlis, D. Enhanced upward heat transport at deep submesoscale ocean fronts. *Nat. Geosci.* **2020**, *13*, 50–55. [CrossRef]

16. Tagklis, F.; Bracco, A.; Ito, T.; Castelao, R. Submesoscale modulation of deep water formation in the Labrador Sea. *Sci. Rep.* **2020**, *10*, 1–13. [CrossRef]

17. Charney, J.G. Geostrophic turbulence. *J. Atmos. Sci.* **1971**, *28*, 1087–1095. [CrossRef]

18. Arbic, B.K.; Polzin, K.L.; Scott, R.B.; Richman, J.G.; Shriver, J.F. On eddy viscosity, energy cascades, and the horizontal resolution of gridded satellite altimeter products. *J. Phys. Oceanogr.* **2013**, *43*, 283–300. [CrossRef]

19. Aluie, H.; Hecht, M.; Vallis, G.K. Mapping the Energy Cascade in the North Atlantic Ocean: The Coarse-Graining Approach. *J. Phys. Oceanogr.* **2018**, *48*, 225–244. [CrossRef]

20. Yang, Y.; McWilliams, J.C.; Liang, X.S.; Zhang, H.; Weisberg, R.H.; Liu, Y.; Menemenlis, D. Spatial and Temporal Characteristics of the Submesoscale Energetics in the Gulf of Mexico. *J. Phys. Oceanogr.* **2020**. [CrossRef]

21. Font, J.; Garcialadona, E.; Gorriz, E. The seasonality of mesoscale motion in the northern current of the western mediterranean-several years of evidence. *Oceanol. Acta* **1995**, *18*, 207–219.

22. Strub, P.T.; James, C. Altimeter-derived variability of surface velocities in the California Current System: 2. Seasonal circulation and eddy statistics. *Deep Sea Res. Part II Top. Stud. Oceanogr.* **2000**, *47*, 831–870. [CrossRef]

23. Eden, C.; Böning, C. Sources of eddy kinetic energy in the Labrador Sea. *J. Phys. Oceanogr.* **2002**, *32*, 3346–3363. [CrossRef]

24. Pujol, M.I.; Larnicol, G. Mediterranean sea eddy kinetic energy variability from 11 years of altimetric data. *J. Mar. Syst.* **2005**, *58*, 121–142. [CrossRef]

25. Jouanno, J.; Sheinbaum, J.; Barnier, B.; Molines, J.M.; Candela, J. Seasonal and interannual modulation of the eddy kinetic energy in the Caribbean Sea. *J. Phys. Oceanogr.* **2012**, *42*, 2041–2055. [CrossRef]

26. Qiu, B.; Chen, S.; Klein, P.; Sasaki, H.; Sasai, Y. Seasonal mesoscale and submesoscale eddy variability along the North Pacific Subtropical Countercurrent. *J. Phys. Oceanogr.* **2014**, *44*, 3079–3098. [CrossRef]

27. Sasaki, H.; Klein, P.; Qiu, B.; Sasai, Y. Impact of oceanic-scale interactions on the seasonal modulation of ocean dynamics by the atmosphere. *Nat. Commun.* **2014**, *5*, 1–8. [CrossRef]

28. Uchida, T.; Abernathey, R.; Smith, S. Seasonality of eddy kinetic energy in an eddy permitting global climate model. *Ocean Model.* **2017**, *118*, 41–58. [CrossRef]

29. Dong, J.; Fox-Kemper, B.; Zhang, H.; Dong, C. The seasonality of submesoscale energy production, content, and cascade. *Geophys. Res. Lett.* **2020**. [CrossRef]

30. Schubert, R.; Jonathan, G.; Greatbatch, R.J.; Baschek, B.; Biastoch, A. The Submesoscale Kinetic Energy Cascade: Mesoscale Absorption of Submesoscale Mixed-Layer Eddies and Frontal Downscale Fluxes. *J. Phys. Oceanogr.* **2020**. [CrossRef]

31. Ajayi, A.; Le Sommer, J.; Chassignet, E.; Molines, J.M.; Xu, X.; Albert, A.; Dewar, W.K. Diagnosing cross-scale kinetic energy exchanges from two submesoscale permitting ocean models. *J. Adv. Model. Earth Syst.* **2021**. [CrossRef]

32. Fox-Kemper, B.; Ferrari, R.; Hallberg, R. Parameterization of mixed layer eddies. Part I: Theory and diagnosis. *J. Phys. Oceanogr.* **2008**, *38*, 1145–1165. [CrossRef]

33. Johnson, L.; Lee, C.M.; D'Asaro, E.A. Global Estimates of Lateral Springtime Restratification. *J. Phys. Oceanogr.* **2016**, *46*, 1555–1573. [CrossRef]

34. Thompson, A.F.; Lazar, A.; Buckingham, C.; Naveira Garabato, A.C.; Damerell, G.M.; Heywood, K.J. Open-Ocean Submesoscale Motions: A Full Seasonal Cycle of Mixed Layer Instabilities from Gliders. *J. Phys. Oceanogr.* **2016**, *46*, 1285–1307. [CrossRef]

35. Rieck, J.K.; Böning, C.W.; Greatbatch, R.J.; Scheinert, M. Seasonal variability of eddy kinetic energy in a global high-resolution ocean model. *Geophys. Res. Lett.* **2015**, *42*, 9379–9386. [CrossRef]

36. Pedlosky, J. *Geophysical Fluid Dynamics*, 2nd ed.; Springer: New York, NY, USA, 1987.

37. Bachman, S.; Fox-Kemper, B.; Bryan, F. A tracer-based inversion method for diagnosing eddy-induced diffusivity and advection. *Ocean Model.* **2015**, *86*, 1–14. [CrossRef]
38. Young, W.R. An Exact Thickness-Weighted Average Formulation of the Boussinesq Equations. *J. Phys. Oceanogr.* **2012**, *42*, 692–707.
[CrossRef]

39. Deremble, B.; Martinez, E.M. MSOM: Multiple Scale Ocean Model; MEOM Research Group: Grenoble, France, 2020. [CrossRef]

40. Popinet, S. A quadtree-adaptive multigrid solver for the Serre-Green-Naghdi equations. *J. Comput. Phys.* **2015**, *302*, 336–358.
[CrossRef]

41. Verron, J. Nudging satellite altimeter data into quasi-geostrophic ocean models. *J. Geophys. Res. Ocean.* **1992**, *97*, 7479–7491.
[CrossRef]

42. Barnier, B.; Le Provost, C. Influence of bottom topography roughness on the jet and inertial recirculation of a mid-latitude gyre. *Dyn. Atmos. Ocean.* **1993**, *18*, 29–65. [CrossRef]

43. Arakawa, A. Computational design for long-term numerical integration of the equations of fluid motion: Two-dimensional incompressible flow. Part I. *J. Comput. Phys.* **1997**, *135*, 103–114. [CrossRef]

44. Karabasov, S.; Berloff, P.S.; Goloviznin, V. CABARET in the ocean gyres. *Ocean Model.* **2009**, *30*, 155–168. [CrossRef]

45. Hallberg, R. Using a resolution function to regulate parameterizations of oceanic mesoscale eddy effects. *Ocean Model.* **2013**, *72*, 92–103. [CrossRef]

46. Rhines, P.B.; Schopp, R. The wind-driven circulation: Quasi-geostrophic simulations and theory for nonsymmetric winds. *J. Phys. Oceanogr.* **1991**, *21*, 1438–1469. [CrossRef]

47. Hogg, A.M.C.; Killworth, P.D.; Blundell, J.R.; Dewar, W.K. Mechanisms of decadal variability of the wind-driven ocean circulation. *J. Phys. Oceanogr.* **2005**, *35*, 512–531. [CrossRef]

48. Simonnet, E. Quantization of the low-frequency variability of the double-gyre circulation. *J. Phys. Oceanogr.* **2005**, *35*, 2268–2290.
[CrossRef]

49. Berloff, P.; Hogg, A.M.C.; Dewar, W. The turbulent oscillator: A mechanism of low-frequency variability of the wind-driven ocean gyres. *J. Phys. Oceanogr.* **2007**, *37*, 2363–2386. [CrossRef]

50. Chelton, D.B.; DeSzoeke, R.A.; Schlax, M.G.; El Naggar, K.; Siwertz, N. Geographical variability of the first baroclinic Rossby radius of deformation. *J. Phys. Oceanogr.* **1998**, *28*, 433–460. [CrossRef]

51. Dijkstra, H.A.; Wubs, F.W.; Cliffe, A.K.; Doedel, E.; Dragomirescu, I.F.; Eckhardt, B.; Gelfgat, A.Y.; Hazel, A.L.; Lucarini, V.; Salinger, A.G.; et al. Numerical bifurcation methods and their application to fluid dynamics: Analysis beyond simulation. *Commun. Comput. Phys.* **2014**, *15*, 1–45. [CrossRef]

52. Hoskins, B.J.; Draghici, I.; Davies, H.C. A new look at the ω-equation. *Q. J. R. Meteorol. Soc.* **1978**, *104*, 31–38. [CrossRef]

53. Gill, A.E. *Atmosphere-Ocean Dynamics*; Academic Press: Cambridge, MA, USA, 1982; Volume 30, 662p.

54. Grooms, I.; Julien, K.; Fox-Kemper, B. On the interactions between planetary geostrophy and mesoscale eddies. *Dyn. Atmos. Ocean.* **2011**, *51*, 109–136. [CrossRef]

55. Jamet, Q.; Deremble, B.; Wienders, N.; Uchida, T.; Dewar, W.K. On Wind-driven Energetics of Subtropical Gyres. *J. Adv. Model. Earth Syst.* **2021**, e2020MS002329. [CrossRef]

56. Sérazin, G.; Jaymond, A.; Leroux, S.; Penduff, T.; Bessières, L.; Llovel, W.; Barnier, B.; Molines, J.M.; Terray, L. A global probabilistic study of the ocean heat content low-frequency variability: Atmospheric forcing versus oceanic chaos. *Geophys. Res. Lett.* **2017**, *44*, 5580–5589. [CrossRef]

57. Leroux, S.; Penduff, T.; Bessières, L.; Molines, J.M.; Brankart, J.M.; Sérazin, G.; Barnier, B.; Terray, L. Intrinsic and Atmospherically Forced Variability of the AMOC: Insights from a Large-Ensemble Ocean Hindcast. *J. Clim.* **2018**, *31*, 1183–1203. [CrossRef]

58. Jamet, Q.; Dewar, W.K.; Wienders, N.; Deremble, B. Spatiotemporal Patterns of Chaos in the Atlantic Overturning Circulation. *Geophys. Res. Lett.* **2019**, *46*, 7509–7517. [CrossRef]

59. Chassignet, E.P.; Marshall, D.P. Gulf Stream separation in numerical ocean models. *Geophys. Monogr. Ser.* **2008**, *177*. [CrossRef]

60. Hoyer, S.; Hamman, J.; Roos, M.; Cherian, D.; Fitzgerald, C.; Fujii, K.; Maussion, F.; Crusaderky, Kleeman, A.; Clark, S.; et al. *xarray*: N-D Labeled Arrays and Datasets; USA, 2020. Available online: http://doi.org/10.5281/zenodo.4299126 (accessed on 9 April 2021)

61. Griffies, S.M. The Gent–McWilliams skew flux. *J. Phys. Oceanogr.* **1998**, *28*, 831–841. [CrossRef]