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

^{*}

Next Article in Journal

Previous Article in Journal

Previous Article in Special Issue

Previous Article in Special Issue

Université Grenoble Alpes, CNRS, IRD, Grenoble-INP, Institut des Géosciences de l’Environnement, CS 40 700, CEDEX 9, 38058 Grenoble, France

Author to whom correspondence should be addressed.

Academic Editor: Pavel S. Berloff

Received: 9 April 2021
/
Revised: 25 May 2021
/
Accepted: 27 May 2021
/
Published: 2 June 2021

(This article belongs to the Collection Geophysical Fluid Dynamics)

With the advent of submesoscale $O(1\phantom{\rule{3.33333pt}{0ex}}\mathrm{km})$ permitting basin-scale ocean simulations, the seasonality of mesoscale $O(50\phantom{\rule{3.33333pt}{0ex}}\mathrm{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(Ro)\ll 1,O(Fr)\ll 1$ ) while submesoscale dynamics are associated with $O(Ro)\sim 1,O(Fr)\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 destablizes and results in overall lower mean and eddy kinetic energies despite the domain being more susceptible to baroclinic instability from weaker stratification.

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,4,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\phantom{\rule{3.33333pt}{0ex}}\mathrm{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\phantom{\rule{3.33333pt}{0ex}}\mathrm{km})$[10]. In terms of the Rossby number ($Ro\phantom{\rule{3.33333pt}{0ex}}(=U/{f}_{0}L)$) and Froude number ($Fr\phantom{\rule{3.33333pt}{0ex}}(=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,12,13,14,15,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,18,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,22,23,24,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,27,28,29,30,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,33,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\phantom{\rule{3.33333pt}{0ex}}\mathrm{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.

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\phantom{\rule{4pt}{0ex}}(=50\phantom{\rule{3.33333pt}{0ex}}\mathrm{km})$ and horizontal resolution is ∼$4\phantom{\rule{3.33333pt}{0ex}}\mathrm{km}\phantom{\rule{4pt}{0ex}}(={\delta}_{\widehat{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,44,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(Ro)\ll 1,O(Fr)\ll 1$).

MSQG solves prognostically for the non-dimensionalized QG potential vorticity (PV):
where $q={\zeta}_{g}+\beta y-\frac{{f}_{0}}{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 $\widehat{(\xb7)}$ 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:
where ${\nabla}_{\mathrm{h}}$ is the horizontal gradient operator, and $\widehat{y}\phantom{\rule{4pt}{0ex}}(\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,47,48,49]. The dissipation term is implemented via a biharmonic viscosity:

$$\frac{D\widehat{q}}{D\widehat{t}}=\widehat{\mathcal{F}}+\widehat{\mathcal{D}},$$

$$\widehat{\mathcal{F}}=\frac{{\widehat{\nabla}}_{\mathrm{h}}\times \widehat{\tau}(\widehat{y})}{{\widehat{H}}_{1}}=-\frac{{\widehat{\tau}}_{0}}{R{o}^{m}{\widehat{H}}_{1}}sin\left(\frac{2\pi}{N}\widehat{y}\right)sin\left(\frac{\pi}{N}\widehat{y}\right),$$

$$\widehat{\mathcal{D}}=-{R{e}_{4}}^{-1}{\widehat{\nabla}}_{\mathrm{h}}^{4}\widehat{q}.$$

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:
where ${H}_{i}^{\u2020}=({H}_{i}+{H}_{i+1})/2$, ${g}^{\prime}$ is the reduced gravity and subscript i is the layer index (Figure 1). We vary $F{r}_{1}$ in time but keep $F{r}_{2}$ stationary (${\widehat{A}}_{F{r}_{2}}=0$), which is consistent with the seasonal variability of stratification being confined in the upper few hundred meters in the real ocean [50].

$$F{r}_{i}=\frac{U}{\sqrt{{g}_{i}^{\prime}{H}_{i}^{\u2020}}}=F{r}_{i}^{m}{\left[1+{\widehat{A}}_{F{r}_{i}}sin(2\pi {\widehat{f}}_{F{r}_{i}}\widehat{t})\right]}^{-1/2},$$

We spin up the model for 10 years from a spun-up run with lower resolution ($N=256$, equivalently ${\delta}_{\widehat{x}}L\sim 15\phantom{\rule{3.33333pt}{0ex}}\mathrm{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 ($\mathcal{T}=L/U=5\times {10}^{5}\phantom{\rule{3.33333pt}{0ex}}\mathrm{sec}\mathrm{onds}\sim $5.8 days).

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\phantom{\rule{3.33333pt}{0ex}}(\in [1,3]$; Figure 1) neglecting viscous and external forcing terms:
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}_{(\xb7)}$ with respect to $(t,z,y,x)$. The stream function is defined as ${\psi}_{i}={\varphi}_{g;i}/{f}_{0}$ where ${\varphi}_{g;i}$ is the geostrophic pressure anomaly and relative vorticity can be written as ${\zeta}_{g;i}={\nabla}_{\mathrm{h}}^{2}{\psi}_{i}$. The layer-thickness equation on the other hand is [7]:
We leave the derivation of the layered QGPV and its relation to the continuously stratified framework to Appendix A.

$$\begin{array}{cc}\hfill {\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}\phantom{\rule{1.em}{0ex}}& =-{f}_{0}({\partial}_{x}{u}_{a;i}+{\partial}_{y}{v}_{a;i})\hfill \\ \hfill \phantom{\rule{1.em}{0ex}}& ={f}_{0}{\partial}_{z}{w}_{a;i},\hfill \end{array}$$

$$\begin{array}{cc}\hfill {\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})\hfill \\ \hfill \phantom{\rule{1.em}{0ex}}& ={H}_{i}{\partial}_{z}{w}_{a;i}\hfill \end{array}$$

The ageostrophic vertical velocity can be diagnosed via the QG omega equation (Appendix B; [52,53]):
where ${N}_{i}^{2}={g}_{i}^{\prime}/{H}_{i}^{\u2020}$, and ${b}_{i}={f}_{0}\frac{{\psi}_{i}-{\psi}_{i+1}}{{H}_{i}^{\u2020}}$ is the buoyancy. The $\mathbf{Q}$ tensor is:
where ${\mathbf{u}}_{g;i}^{\u2020}=-{\partial}_{y}{\psi}_{i}^{\u2020}\mathbf{i}+{\partial}_{x}{\psi}_{i}^{\u2020}\mathbf{j}$ is the geostrophic velocity derived from the inter-facial stream function $({\psi}_{i}^{\u2020}=\frac{{H}_{i}{\psi}_{i+1}+{H}_{i+1}{\psi}_{i}}{{H}_{i}+{H}_{i+1}}$; [3]). $\mathbf{i}$ and $\mathbf{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}$):
where ${\delta}_{E}=E{k}^{b}{H}_{3}$ is the bottom Ekman-layer thickness.

$${N}_{i}^{2}{\nabla}_{\mathrm{h}}^{2}{w}_{a;i}+{f}_{0}^{2}{\partial}_{zz}{w}_{a;i}=\beta {\partial}_{x}{b}_{i}-2{\nabla}_{\mathrm{h}}\xb7{\mathbf{Q}}_{i}-{\nabla}_{\mathrm{h}}^{2}{b}_{i}{N}_{i}^{2}{\partial}_{t}\frac{1}{{N}_{i}^{2}},$$

$${\mathbf{Q}}_{i}={Q}_{i}^{1}\mathbf{i}+{Q}_{i}^{2}\mathbf{j}=\left({\partial}_{x}{\mathbf{u}}_{g;i}^{\u2020}\xb7{\nabla}_{\mathrm{h}}{b}_{i}\right)\mathbf{i}+\left({\partial}_{y}{\mathbf{u}}_{g;i}^{\u2020}\xb7{\nabla}_{\mathrm{h}}{b}_{i}\right)\mathbf{j},$$

$${w}_{E;0}=-\frac{1}{{f}_{0}}{\nabla}_{\mathrm{h}}\times \tau =-\frac{UH}{L}{\widehat{\tau}}_{0}{sin}^{2}\left[\frac{2\pi y}{{L}_{0}}\right]sin\left[\frac{\pi y}{{L}_{0}}\right],$$

$${w}_{E;3}=\frac{{\delta}_{E}}{2}{\zeta}_{g;3},$$

Now, multiplying Equation (5) by $-{\psi}_{i}$ and integrating over the depth of each layer gives the kinetic energy (KE) budget:
Dropping the divergence terms as they vanish upon area integration, for each layer, we get:

$$\begin{array}{cc}\hfill {H}_{i}[\frac{{D}_{i}}{Dt}& \frac{|{\nabla}_{\mathrm{h}}{\psi}_{i}{|}^{2}}{2}-{\nabla}_{\mathrm{h}}\xb7({\mathbf{u}}_{g;i}{\psi}_{i}{\nabla}_{\mathrm{h}}^{2}{\psi}_{i})-\frac{\beta}{2}{\partial}_{x}{\psi}_{i}^{2}]\hfill \\ \hfill \phantom{\rule{1.em}{0ex}}& =-{f}_{0}\int {\psi}_{i}{\partial}_{z}{w}_{a;i}dz\hfill \\ \hfill \phantom{\rule{1.em}{0ex}}& ={f}_{0}\left[-({w}_{a;i-1}{\psi}_{i-1}^{\u2020}-{w}_{a;i}{\psi}_{i}^{\u2020})+\int {w}_{a;i}{\partial}_{z}{\psi}_{i}dz\right]\hfill \\ \hfill \phantom{\rule{1.em}{0ex}}& ={f}_{0}\left[-({w}_{a;i-1}{\psi}_{i-1}^{\u2020}-{w}_{a;i}{\psi}_{i}^{\u2020})+{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].\hfill \end{array}$$

$$\begin{array}{cc}\hfill \frac{{H}_{1}}{2}{\partial}_{t}{\left|{\nabla}_{\mathrm{h}}{\psi}_{1}\right|}^{2}\phantom{\rule{1.em}{0ex}}& ={f}_{0}\left[{w}_{a;1}{\psi}_{1}^{\u2020}+{w}_{a;1}{H}_{1}\frac{{\psi}_{1}-{\psi}_{2}}{{H}_{2}+{H}_{1}}\right],\hfill \end{array}$$

$$\begin{array}{cc}\hfill \frac{{H}_{2}}{2}{\partial}_{t}{\left|{\nabla}_{\mathrm{h}}{\psi}_{2}\right|}^{2}\phantom{\rule{1.em}{0ex}}& ={f}_{0}\left[-({w}_{a;1}{\psi}_{1}^{\u2020}-{w}_{a;2}{\psi}_{2}^{\u2020})+{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],\hfill \end{array}$$

$$\begin{array}{cc}\hfill \frac{{H}_{3}}{2}{\partial}_{t}{\left|{\nabla}_{\mathrm{h}}{\psi}_{3}\right|}^{2}\phantom{\rule{1.em}{0ex}}& ={f}_{0}\left[-{w}_{a;2}{\psi}_{2}^{\u2020}+{w}_{a;2}{H}_{3}\frac{{\psi}_{2}-{\psi}_{3}}{{H}_{3}+{H}_{2}}\right].\hfill \end{array}$$

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}$:
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}$:
where $\frac{{D}_{i}^{\u2020}}{Dt}={\partial}_{t}+{\mathbf{u}}_{g;i}^{\u2020}\xb7{\nabla}_{\mathrm{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:
and Equation (16) with ${f}_{0}({\psi}_{2}-{\psi}_{3})$:

$$\begin{array}{cc}\hfill \frac{{D}_{1}^{\u2020}}{Dt}\left[\frac{{f}_{0}}{{g}_{1}^{\prime}}({\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}^{\prime}}({\psi}_{3}-{\psi}_{2})\right]\hfill \\ \hfill \phantom{\rule{1.em}{0ex}}& =-{w}_{a;1}+\frac{{f}_{0}{H}_{1}}{{g}_{2}^{\prime}({H}_{1}+{H}_{2})}({\mathbf{u}}_{3}-{\mathbf{u}}_{2})\xb7{\nabla}_{\mathrm{h}}({\psi}_{3}-{\psi}_{2})\hfill \\ \hfill \phantom{\rule{1.em}{0ex}}& =-{w}_{a;1},\hfill \end{array}$$

$$\begin{array}{cc}\hfill \frac{{D}_{2}^{\u2020}}{Dt}\left[\frac{{f}_{0}}{{g}_{2}^{\prime}}({\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}^{\prime}}({\psi}_{2}-{\psi}_{1})\right]\hfill \\ \hfill \phantom{\rule{1.em}{0ex}}& =-{w}_{a;2},\hfill \end{array}$$

$$\begin{array}{c}\hfill {\partial}_{t}\left[\frac{{f}_{0}^{2}}{2{g}_{1}^{\prime}}{({\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}^{\prime}}^{-1},\end{array}$$

$$\begin{array}{c}\hfill {\partial}_{t}\left[\frac{{f}_{0}^{2}}{2{g}_{2}^{\prime}}{({\psi}_{2}-{\psi}_{3})}^{2}\right]=-{f}_{0}({\psi}_{2}-{\psi}_{3}){w}_{a;2}.\end{array}$$

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}^{\prime}$ is inside the derivative.

Now, the mean KE (MKE; ${K}^{\#}$), eddy KE (EKE; $\mathcal{K}$), mean APE (MAPE; ${P}^{\#}$) and eddy APE (EAPE; $\mathcal{P}$) can be defined as:
where $\overline{(\xb7)}$ is the ensemble mean and the eddy is defined as fluctuations about the ensemble mean, viz. ${(\xb7)}^{\prime}=(\xb7)-\overline{(\xb7)}$. We note that the ensemble mean of the fluctuations vanish ($\overline{{(\xb7)}^{\prime}}=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,57,58]. The ensemble means of total KE and APE each satisfy $\overline{{K}_{i}}=\frac{{H}_{i}}{2}\overline{|{\nabla}_{\mathrm{h}}{\psi}_{i}{|}^{2}}={K}_{i}^{\#}+{\mathcal{K}}_{i}$, $\overline{{P}_{i}}=\frac{{f}_{0}^{2}}{2{g}_{i}^{\prime}}\overline{{({\psi}_{i}-{\psi}_{i+1})}^{2}}={P}_{i}^{\#}+{\mathcal{P}}_{i}$. Hence, the exchanges ($\Pi $) of KE and APE within and between layers are:
where $\langle \xb7\rangle =\phantom{\rule{0.277778em}{0ex}}\int \int \phantom{\rule{0.277778em}{0ex}}(\xb7)dxdy$ is the area integration. Further details regarding the sign convention and forcing/dissipation terms are given in Appendix C and Appendix D. Summing up each layer gives the volume integrated energy exchanges:

$${K}_{i}^{\#}=\frac{{H}_{i}}{2}{\left|{\nabla}_{\mathrm{h}}\overline{{\psi}_{i}}\right|}^{2},\phantom{\rule{4pt}{0ex}}{\mathcal{K}}_{i}=\frac{{H}_{i}}{2}\overline{|{\nabla}_{\mathrm{h}}{\psi}_{i}^{\prime}{|}^{2}},$$

$${P}_{i}^{\#}=\frac{{f}_{0}^{2}}{2{g}_{i}^{\prime}}{(\overline{{\psi}_{i}}-\overline{{\psi}_{i+1}})}^{2},\phantom{\rule{4pt}{0ex}}{\mathcal{P}}_{i}=\frac{{f}_{0}^{2}}{2{g}_{i}^{\prime}}\overline{{({\psi}_{i}^{\prime}-{\psi}_{i+1}^{\prime})}^{2}},$$

$${\Pi}_{{K}_{1}^{\#}\to {\mathcal{K}}_{1}}=-{H}_{1}\langle \overline{{\psi}_{1}}{\nabla}_{\mathrm{h}}\xb7\overline{{\mathbf{u}}_{g;1}^{\prime}{\nabla}_{\mathrm{h}}^{2}{\psi}_{1}^{\prime}}\rangle ,$$

$${\Pi}_{{K}_{1}^{\#}\to {K}_{2}^{\#}}=-{f}_{0}\langle \overline{{w}_{a;1}}\overline{{\psi}_{1}^{\u2020}}\rangle ,\phantom{\rule{4pt}{0ex}}\phantom{\rule{4pt}{0ex}}{\Pi}_{{\mathcal{K}}_{1}\to {\mathcal{K}}_{2}}=-{f}_{0}\langle \overline{{w}_{a;1}^{\prime}{{\psi}_{1}^{\u2020}}^{\prime}}\rangle ,$$

$${\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 ,\phantom{\rule{4pt}{0ex}}\phantom{\rule{4pt}{0ex}}{\Pi}_{{\mathcal{P}}_{1}\to {\mathcal{K}}_{1}}=\frac{{f}_{0}{H}_{1}}{{H}_{2}+{H}_{1}}\langle \overline{{w}_{a;1}^{\prime}({\psi}_{1}^{\prime}-{\psi}_{2}^{\prime})}\rangle ,$$

$${\Pi}_{{P}_{1}^{\#}\to {\mathcal{P}}_{1}}=\frac{{f}_{0}^{2}}{{g}_{1}^{\prime}}\langle (\overline{{\psi}_{1}}-\overline{{\psi}_{2}}){\nabla}_{\mathrm{h}}\xb7\overline{{{\mathbf{u}}_{g;1}^{\u2020}}^{\prime}({\psi}_{1}^{\prime}-{\psi}_{2}^{\prime})}\rangle ,$$

$${\Pi}_{{B}_{1}^{\#}\to {P}_{1}^{\#}}=-\frac{{f}_{0}^{2}}{2}\langle {(\overline{{\psi}_{1}}-\overline{{\psi}_{2}})}^{2}{\partial}_{t}{{g}_{1}^{\prime}}^{-1}\rangle ,\phantom{\rule{4pt}{0ex}}\phantom{\rule{4pt}{0ex}}{\Pi}_{{B}_{1}^{\#}\to {\mathcal{P}}_{1}}=-\frac{{f}_{0}^{2}}{2}\langle \overline{{({\psi}_{1}^{\prime}-{\psi}_{2}^{\prime})}^{2}}{\partial}_{t}{{g}_{1}^{\prime}}^{-1}\rangle ,$$

$${\Pi}_{{K}_{2}^{\#}\to {\mathcal{K}}_{2}}=-{H}_{2}\langle \overline{{\psi}_{2}}{\nabla}_{\mathrm{h}}\xb7\overline{{\mathbf{u}}_{g;2}^{\prime}{\nabla}_{\mathrm{h}}^{2}{\psi}_{2}^{\prime}}\rangle ,$$

$${\Pi}_{{K}_{2}^{\#}\to {K}_{3}^{\#}}=-{f}_{0}\langle \overline{{w}_{a;2}}\overline{{\psi}_{2}^{\u2020}}\rangle ,\phantom{\rule{4pt}{0ex}}\phantom{\rule{4pt}{0ex}}{\Pi}_{{\mathcal{K}}_{2}\to {\mathcal{K}}_{3}}=-{f}_{0}\langle \overline{{w}_{a;2}^{\prime}{{\psi}_{2}^{\u2020}}^{\prime}}\rangle ,$$

$${\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 ,\phantom{\rule{4pt}{0ex}}\phantom{\rule{4pt}{0ex}}{\Pi}_{{\mathcal{P}}_{1}\to {\mathcal{K}}_{2}}=\frac{{f}_{0}{H}_{2}}{{H}_{2}+{H}_{1}}\langle \overline{{w}_{a;1}^{\prime}({\psi}_{1}^{\prime}-{\psi}_{2}^{\prime})}\rangle ,$$

$${\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 ,\phantom{\rule{4pt}{0ex}}\phantom{\rule{4pt}{0ex}}{\Pi}_{{\mathcal{P}}_{2}\to {\mathcal{K}}_{2}}=\frac{{f}_{0}{H}_{2}}{{H}_{3}+{H}_{2}}\langle \overline{{w}_{a;2}^{\prime}({\psi}_{2}^{\prime}-{\psi}_{3}^{\prime})}\rangle ,$$

$${\Pi}_{{P}_{2}^{\#}\to {\mathcal{P}}_{2}}=\frac{{f}_{0}^{2}}{{g}_{2}^{\prime}}\langle (\overline{{\psi}_{2}}-\overline{{\psi}_{3}}){\nabla}_{\mathrm{h}}\xb7\overline{{{\mathbf{u}}_{g;2}^{\u2020}}^{\prime}({\psi}_{2}^{\prime}-{\psi}_{3}^{\prime})}\rangle ,$$

$${\Pi}_{{K}_{3}^{\#}\to {\mathcal{K}}_{3}}=-{H}_{3}\langle \overline{{\psi}_{3}}{\nabla}_{\mathrm{h}}\xb7\overline{{\mathbf{u}}_{g;3}^{\prime}{\nabla}_{\mathrm{h}}^{2}{\psi}_{3}^{\prime}}\rangle ,$$

$${\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 ,\phantom{\rule{4pt}{0ex}}\phantom{\rule{4pt}{0ex}}{\Pi}_{{\mathcal{P}}_{2}\to {\mathcal{K}}_{3}}=\frac{{f}_{0}{H}_{3}}{{H}_{3}+{H}_{2}}\langle \overline{{w}_{a;2}^{\prime}({\psi}_{2}^{\prime}-{\psi}_{3}^{\prime})}\rangle ,$$

$${\Pi}_{{P}^{\#}\to {K}^{\#}}=\sum _{i=1}^{2}{f}_{0}\u2329\overline{{w}_{a;i}}(\overline{{\psi}_{i}}-\overline{{\psi}_{i+1}})\u232a,$$

$${\Pi}_{\mathcal{P}\to \mathcal{K}}=\sum _{i=1}^{2}{f}_{0}\u2329\overline{{w}_{a;i}^{\prime}({\psi}_{i}^{\prime}-{\psi}_{i+1}^{\prime})}\u232a,$$

$${\Pi}_{{P}^{\#}\to \mathcal{P}}=\sum _{i=1}^{2}\frac{{f}_{0}^{2}}{{g}_{i}^{\prime}}\u2329(\overline{{\psi}_{i}}-\overline{{\psi}_{i+1}}){\nabla}_{\mathrm{h}}\xb7\overline{{{\mathbf{u}}_{g;i}^{\u2020}}^{\prime}({\psi}_{i}^{\prime}-{\psi}_{i+1}^{\prime})}\u232a,$$

$${\Pi}_{{K}^{\#}\to \mathcal{K}}=-\sum _{i=1}^{3}{H}_{i}\u2329\overline{{\psi}_{i}}{\nabla}_{\mathrm{h}}\xb7\overline{{\mathbf{u}}_{g;i}^{\prime}{\nabla}_{\mathrm{h}}^{2}{\psi}_{i}^{\prime}}\u232a,$$

$${\Pi}_{{B}^{\#}\to {P}^{\#}}=-\frac{{f}_{0}^{2}}{2}\u2329{(\overline{{\psi}_{1}}-\overline{{\psi}_{2}})}^{2}{\partial}_{t}{{g}_{1}^{\prime}}^{-1}\u232a,$$

$${\Pi}_{{B}^{\#}\to \mathcal{P}}=-\frac{{f}_{0}^{2}}{2}\u2329\overline{{({\psi}_{1}^{\prime}-{\psi}_{2}^{\prime})}^{2}}{\partial}_{t}{{g}_{1}^{\prime}}^{-1}\u232a.$$

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 ($\langle {K}_{1}\rangle $), most relevant for studies interested in surface seasonal dynamics, is in sync with the background stratification (${g}_{1}^{\prime}$), viz. higher $\langle {K}_{1}\rangle $ 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 ($\langle |{\nabla}_{\mathrm{h}}\Psi {|}^{2}\rangle $ where $\Psi ={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.

In Figure 3, we show the mean and eddy KE in the first layer (${K}_{1}^{\#},{\mathcal{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}^{\prime}$) from the CTRL run from which we see coherent features of mesoscale eddies (Figure 3g,h).

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}^{\prime}$ in its denominator while MKE lags ${g}_{1}^{\prime}$ 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 ($\mathcal{P},\mathcal{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}_{\mathcal{P}\to \mathcal{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}_{\mathcal{P}\to \mathcal{K}}<0$), whereas baroclinic instability would predict EAPE→EKE (${\Pi}_{\mathcal{P}\to \mathcal{K}}>0$). Regarding the dissipation terms, only the bottom drag for EKE (${D}_{b}^{\mathcal{K}}$) shows a notable seasonality and has a similar magnitude to the energy flux from MKE to EKE (${\Pi}_{{K}^{\#}\to \mathcal{K}}$). The amplitude of bottom drag ($|{D}_{b}^{\mathcal{K}}|$) lags EKE by ∼41 days and aligns well with the ensemble-mean barotropic TKE (Figure 2a).

To provide a climatological view of the energy fluxes ($\Pi $), 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 \mathcal{K}}$) is highest during the summer (Figure 4b and Figure 5a). We attribute the larger energy conversion from EAPE to EKE during the winter (${\Pi}_{\mathcal{P}\to \mathcal{K}}$) to the flow being more susceptible to baroclinic instability with reduced stratification.

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}^{\prime}$) implied from Figure 2. It is perhaps interesting to note that although the ensemble-mean barotropic TKE lags ${g}_{1}^{\prime}$ 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}^{\prime}$, 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 ${\mathcal{K}}_{2}$ and ∼68 days for ${\mathcal{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 {\mathcal{K}}_{2,3}}$). The vertical transfer of EKE (${\Pi}_{{\mathcal{K}}_{1}\to {\mathcal{K}}_{2}},{\Pi}_{{\mathcal{K}}_{2}\to {\mathcal{K}}_{3}}$) and conversion from EAPE (${\Pi}_{\mathcal{P}\to \mathcal{K}}$), on the other hand, show a coherent seasonal pattern with the maxima of ${\mathcal{K}}_{2}$ and ${\mathcal{K}}_{3}$ falling in between the maxima of the two fluxes. We, therefore, attribute the time lag in ${\mathcal{K}}_{2}$ and ${\mathcal{K}}_{3}$ to the balance between baroclinic instability and vertical transfers of EKE.

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,28,29,30,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}_{\mathrm{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 (Figure 4b and Figure 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.

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.

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.

Not applicable.

Not applicable.

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.

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.

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.

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:
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}^{\prime}=\frac{{U}^{2}}{{H}_{1}^{\u2020}}F{r}_{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}^{\prime}}({\psi}_{i+1}-{\psi}_{i})$. The layer thickness can, therefore, be written using the stream function as [7]:
where $\frac{{D}_{1}}{Dt}{\eta}_{0}=\frac{{D}_{3}}{Dt}{\eta}_{3}=0$ due to rigid-lid and flat-bottom boundary conditions.

$$\frac{{D}_{i}}{Dt}{\zeta}_{g;i}+\beta {v}_{g;i}=\frac{{f}_{0}}{{H}_{i}}\frac{{D}_{i}}{Dt}{h}_{i},$$

$$\begin{array}{cc}\hfill {h}_{i}& ={H}_{i}+{\eta}_{i-1}-{\eta}_{i}\hfill \\ \hfill \phantom{\rule{1.em}{0ex}}& ={H}_{i}+\frac{{f}_{0}}{{g}_{i-1}^{\prime}}({\psi}_{i}-{\psi}_{i-1})-\frac{{f}_{0}}{{g}_{i}^{\prime}}({\psi}_{i+1}-{\psi}_{i}),\hfill \end{array}$$

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]:
where $b=B-{B}_{0}$ is the QG fluctuation about the background buoyancy (${B}_{0}$). Hence, we get:
and taking the material derivative gives the buoyancy equation in the continuously stratified framework:

$$\begin{array}{cc}\hfill \eta & =\frac{\partial z}{\partial B}{|}_{z=H}[{B}_{0}(t,z=H)-B(t,z=H+\eta )]\hfill \\ \hfill \phantom{\rule{1.em}{0ex}}& =-\frac{\partial z}{\partial B}{|}_{z=H}b,\hfill \end{array}$$

$$\frac{b}{{N}^{2}}=-\eta ,$$

$$\frac{D}{Dt}\frac{b}{{N}^{2}}=-w.$$

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

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:
and multiplying them by ${f}_{0}$ gives:
and the horizontal gradients of the buoyancy Equation (A5) with the diffusive term yields:

$${\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}{\varphi}_{a;i},$$

$${\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}{\varphi}_{a;i},$$

$$\frac{D}{Dt}({f}_{0}{\partial}_{z}{u}_{g})+{\partial}_{y}{\mathbf{u}}_{g}\xb7{\nabla}_{\mathrm{h}}b-{f}_{0}^{2}{\partial}_{z}{v}_{a}-\beta y{f}_{0}{\partial}_{z}{v}_{g}=-{f}_{0}{\nu}_{4}{\partial}_{z}{\nabla}_{\mathrm{h}}^{4}{u}_{g},$$

$$\frac{D}{Dt}({f}_{0}{\partial}_{z}{v}_{g})-{\partial}_{x}{\mathbf{u}}_{g}\xb7{\nabla}_{\mathrm{h}}b+{f}_{0}^{2}{\partial}_{z}{u}_{a}+\beta y{f}_{0}{\partial}_{z}{u}_{g}=-{f}_{0}{\nu}_{4}{\partial}_{z}{\nabla}_{\mathrm{h}}^{4}{v}_{g},$$

$$\frac{1}{{N}^{2}}\frac{D}{Dt}{\partial}_{x}b+{\partial}_{x}b{\partial}_{t}\frac{1}{{N}^{2}}+\frac{{\partial}_{x}{\mathbf{u}}_{g}}{{N}^{2}}\xb7{\nabla}_{\mathrm{h}}b+{\partial}_{x}{w}_{a}=-{\kappa}_{4}{\partial}_{x}{\nabla}_{\mathrm{h}}^{4}b,$$

$$\frac{1}{{N}^{2}}\frac{D}{Dt}{\partial}_{y}b+{\partial}_{y}b{\partial}_{t}\frac{1}{{N}^{2}}+\frac{{\partial}_{y}{\mathbf{u}}_{g}}{{N}^{2}}\xb7{\nabla}_{\mathrm{h}}b+{\partial}_{y}{w}_{a}=-{\kappa}_{4}{\partial}_{y}{\nabla}_{\mathrm{h}}^{4}b.$$

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

$$2{\partial}_{y}{\mathbf{u}}_{g}\xb7{\nabla}_{\mathrm{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.$$

$$2{\partial}_{x}{\mathbf{u}}_{g}\xb7{\nabla}_{\mathrm{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.$$

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}_{\mathrm{h}}^{2}b$) and their parameters being set identical (viz. ${\nu}_{4}={\kappa}_{4}\phantom{\rule{4pt}{0ex}}(=R{e}_{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}_{\mathrm{h}}^{2}{w}_{a}+{f}_{0}^{2}{\partial}_{zz}{w}_{a}=\beta {\partial}_{x}b-2{\nabla}_{\mathrm{h}}\xb7\mathbf{Q}-{\nabla}_{\mathrm{h}}^{2}b{N}^{2}{\partial}_{t}\frac{1}{{N}^{2}}.$$

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).

In this section, we derive the mean and eddy KE equations. Equation (5) can be split into its mean and eddy component:
where $\frac{{D}^{\#}}{Dt}={\partial}_{t}+\overline{{\mathbf{u}}_{g}}\xb7{\nabla}_{\mathrm{h}}$. Multiplying this by $-\overline{\psi}$ gives:
and taking its ensemble mean yields the mean KE equation:

$$\frac{{D}^{\#}}{Dt}{\nabla}_{\mathrm{h}}^{2}\overline{\psi}+\frac{{D}^{\#}}{Dt}{\nabla}_{\mathrm{h}}^{2}{\psi}^{\prime}+{\mathbf{u}}_{g}^{\prime}\xb7{\nabla}_{\mathrm{h}}[{\nabla}_{\mathrm{h}}^{2}(\overline{\psi}+{\psi}^{\prime})]+\beta {\partial}_{x}(\overline{\psi}+{\psi}^{\prime})={f}_{0}{\partial}_{z}(\overline{w}+{w}^{\prime}),$$

$$\begin{array}{cc}\hfill \frac{{D}^{\#}}{Dt}\frac{|{\nabla}_{\mathrm{h}}\overline{\psi}{|}^{2}}{2}-{\nabla}_{\mathrm{h}}\xb7\overline{{\mathbf{u}}_{g}}\overline{\psi}{\nabla}_{\mathrm{h}}^{2}\overline{\psi}-\overline{\psi}\frac{{D}^{\#}}{Dt}{\nabla}_{\mathrm{h}}^{2}{\psi}^{\prime}& -\overline{\psi}{\mathbf{u}}_{g}^{\prime}\xb7{\nabla}_{\mathrm{h}}[{\nabla}_{\mathrm{h}}^{2}(\overline{\psi}+{\psi}^{\prime})]\hfill \\ \hfill \phantom{\rule{1.em}{0ex}}& -\beta {\partial}_{x}\frac{{\overline{\psi}}^{2}}{2}-\overline{\psi}\beta {\partial}_{x}{\psi}^{\prime}=\overline{w}\overline{b}-\overline{\psi}{f}_{0}{\partial}_{z}{w}^{\prime},\hfill \end{array}$$

$$\frac{{D}^{\#}}{Dt}\frac{|{\nabla}_{\mathrm{h}}\overline{\psi}{|}^{2}}{2}-{\nabla}_{\mathrm{h}}\xb7\overline{{\mathbf{u}}_{g}}\overline{\psi}{\nabla}_{\mathrm{h}}^{2}\overline{\psi}-\beta {\partial}_{x}\frac{{\overline{\psi}}^{2}}{2}=\overline{w}\overline{b}+\overline{\psi}{\nabla}_{\mathrm{h}}\xb7\overline{{\mathbf{u}}_{g}^{\prime}{\nabla}_{\mathrm{h}}^{2}{\psi}^{\prime}}.$$

On the other hand, the ensemble mean of the total KE Equation (11) is:
which can be expanded as:

$$\overline{\frac{D}{Dt}\frac{|{\nabla}_{\mathrm{h}}{\psi |}^{2}}{2}}-{\nabla}_{\mathrm{h}}\xb7\overline{{\mathbf{u}}_{g}\psi {\nabla}_{\mathrm{h}}^{2}\psi}-\beta {\partial}_{x}\frac{\overline{{\psi}^{2}}}{2}=\overline{wb},$$

$$\begin{array}{cc}\hfill \frac{{D}^{\#}}{Dt}\frac{|{\nabla}_{\mathrm{h}}\overline{\psi}{|}^{2}}{2}+\frac{{D}^{\#}}{Dt}\frac{\overline{|{\nabla}_{\mathrm{h}}{\psi}^{\prime}{|}^{2}}}{2}& +{\nabla}_{\mathrm{h}}\xb7\overline{{\mathbf{u}}_{g}^{\prime}\frac{|{\nabla}_{\mathrm{h}}{\psi}^{\prime}{|}^{2}}{2}}+\overline{{\mathbf{u}}_{g}^{\prime}\xb7{\nabla}_{\mathrm{h}}({\partial}_{x}\overline{\psi}{\partial}_{x}{\psi}^{\prime}+{\partial}_{y}\overline{\psi}{\partial}_{y}{\psi}^{\prime})}\hfill \\ \hfill \phantom{\rule{1.em}{0ex}}& -{\nabla}_{\mathrm{h}}\xb7\overline{{\mathbf{u}}_{g}\psi {\nabla}_{\mathrm{h}}^{2}\psi}-\beta {\partial}_{x}\frac{\overline{{\psi}^{2}}}{2}=\overline{wb}.\hfill \end{array}$$

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

$$\begin{array}{cc}\hfill \frac{{D}^{\#}}{Dt}\frac{\overline{|{\nabla}_{\mathrm{h}}{\psi}^{\prime}{|}^{2}}}{2}& +{\nabla}_{\mathrm{h}}\xb7\overline{{\mathbf{u}}_{g}^{\prime}\frac{|{\nabla}_{\mathrm{h}}{\psi}^{\prime}{|}^{2}}{2}}+{\nabla}_{\mathrm{h}}\xb7\overline{{\mathbf{u}}_{g}^{\prime}({\partial}_{x}\overline{\psi}{\partial}_{x}{\psi}^{\prime}+{\partial}_{y}\overline{\psi}{\partial}_{y}{\psi}^{\prime})}\hfill \\ \hfill \phantom{\rule{1.em}{0ex}}& -{\nabla}_{\mathrm{h}}\xb7\left(\overline{{\mathbf{u}}_{g}\psi {\nabla}_{\mathrm{h}}^{2}\psi}-\overline{{\mathbf{u}}_{g}}\overline{\psi}{\nabla}_{\mathrm{h}}^{2}\overline{\psi}\right)-\beta {\partial}_{x}\frac{\overline{{{\psi}^{\prime}}^{2}}}{2}=\overline{{w}^{\prime}{b}^{\prime}}-\overline{\psi}{\nabla}_{\mathrm{h}}\xb7\overline{{\mathbf{u}}_{g}^{\prime}{\nabla}_{\mathrm{h}}^{2}{\psi}^{\prime}}.\hfill \end{array}$$

Since the divergence terms vanish upon area integration, we can see the mean and eddy KE exchanging the term $-\overline{\psi}{\nabla}_{\mathrm{h}}\xb7\overline{{\mathbf{u}}_{g}^{\prime}{\nabla}_{\mathrm{h}}^{2}{\psi}^{\prime}}$ (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.

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:

$$\begin{array}{cc}\hfill {\partial}_{t}{K}_{1}^{\#}={f}_{0}\left[\overline{{w}_{a;1}}\overline{{\psi}_{1}^{\u2020}}+{H}_{1}\overline{{w}_{a;1}}\frac{\overline{{\psi}_{1}}-\overline{{\psi}_{2}}}{{H}_{2}+{H}_{1}}\right]& +{H}_{1}\overline{{\psi}_{1}}{\nabla}_{\mathrm{h}}\xb7\overline{{\mathbf{u}}_{g;1}^{\prime}{\nabla}_{\mathrm{h}}^{2}{\psi}_{1}^{\prime}}\hfill \\ \hfill \phantom{\rule{1.em}{0ex}}& -\overline{{\psi}_{1}}{\nabla}_{\mathrm{h}}\times \tau +{H}_{1}\overline{{\psi}_{1}}{\nu}_{4}{\nabla}_{\mathrm{h}}^{4}({\nabla}_{\mathrm{h}}^{2}\overline{{\psi}_{1}}),\hfill \end{array}$$

$${\partial}_{t}{\mathcal{K}}_{1}={f}_{0}\left[\overline{{w}_{a;1}^{\prime}{{\psi}_{1}^{\u2020}}^{\prime}}+{H}_{1}\overline{{w}_{a;1}^{\prime}\frac{{\psi}_{1}^{\prime}-{\psi}_{2}^{\prime}}{{H}_{2}+{H}_{1}}}\right]-{H}_{1}\overline{{\psi}_{1}}{\nabla}_{\mathrm{h}}\xb7\overline{{\mathbf{u}}_{g;1}^{\prime}{\nabla}_{\mathrm{h}}^{2}{\psi}_{1}^{\prime}}+{H}_{1}\overline{{\psi}_{1}^{\prime}{\nu}_{4}{\nabla}_{\mathrm{h}}^{4}({\nabla}_{\mathrm{h}}^{2}{\psi}_{1}^{\prime})},$$

$$\begin{array}{cc}\hfill {\partial}_{t}{P}_{1}^{\#}=-{f}_{0}\overline{{w}_{a;1}}(\overline{{\psi}_{1}}-\overline{{\psi}_{2}})& -\frac{{f}_{0}^{2}}{{g}_{1}^{\prime}}(\overline{{\psi}_{1}}-\overline{{\psi}_{2}}){\nabla}_{\mathrm{h}}\xb7\overline{{\mathbf{u}}_{{g;1}^{\prime}}^{\u2020}({\psi}_{1}^{\prime}-{\psi}_{2}^{\prime})}\hfill \\ \hfill \phantom{\rule{1.em}{0ex}}& -\frac{{f}_{0}^{2}}{{g}_{1}^{\prime}}(\overline{{\psi}_{1}}-\overline{{\psi}_{2}}){\kappa}_{4}{\nabla}_{\mathrm{h}}^{4}(\overline{{\psi}_{1}}-\overline{{\psi}_{2}})-\frac{{f}_{0}^{2}{(\overline{{\psi}_{1}}-\overline{{\psi}_{2}})}^{2}}{2}{\partial}_{t}{{g}_{1}^{\prime}}^{-1},\hfill \end{array}$$

$$\begin{array}{cc}\hfill {\partial}_{t}{\mathcal{P}}_{1}=-{f}_{0}\overline{{w}_{a;1}^{\prime}({\psi}_{1}^{\prime}-{\psi}_{2}^{\prime})}& +\frac{{f}_{0}^{2}}{{g}_{1}^{\prime}}(\overline{{\psi}_{1}}-\overline{{\psi}_{2}}){\nabla}_{\mathrm{h}}\xb7\overline{{\mathbf{u}}_{{g;1}^{\prime}}^{\u2020}({\psi}_{1}^{\prime}-{\psi}_{2}^{\prime})}\hfill \\ \hfill \phantom{\rule{1.em}{0ex}}& -\frac{{f}_{0}^{2}}{{g}_{1}^{\prime}}\overline{({\psi}_{1}^{\prime}-{\psi}_{2}^{\prime}){\kappa}_{4}{\nabla}_{\mathrm{h}}^{4}({\psi}_{1}^{\prime}-{\psi}_{2}^{\prime})}-\frac{{f}_{0}^{2}\overline{{({\psi}_{1}^{\prime}-{\psi}_{2}^{\prime})}^{2}}}{2}{\partial}_{t}{{g}_{1}^{\prime}}^{-1}.\hfill \end{array}$$

For the second layer:

$$\begin{array}{cc}\hfill {\partial}_{t}{K}_{2}^{\#}={f}_{0}[-(\overline{{w}_{a;1}}\overline{{\psi}_{1}^{\u2020}}-\overline{{w}_{a;2}}\overline{{\psi}_{2}^{\u2020}})\phantom{\rule{1.em}{0ex}}& +{H}_{2}\left(\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)]\hfill \\ \hfill \phantom{\rule{1.em}{0ex}}& +{H}_{2}\overline{{\psi}_{2}}{\nabla}_{\mathrm{h}}\xb7\overline{{\mathbf{u}}_{g;2}^{\prime}{\nabla}_{\mathrm{h}}^{2}{\psi}_{2}^{\prime}}+{H}_{2}\overline{{\psi}_{2}}{\nu}_{4}{\nabla}_{\mathrm{h}}^{4}({\nabla}_{\mathrm{h}}^{2}\overline{{\psi}_{2}}),\hfill \end{array}$$

$$\begin{array}{cc}\hfill {\partial}_{t}{\mathcal{K}}_{2}={f}_{0}[-(\overline{{w}_{a;1}^{\prime}{{\psi}_{1}^{\u2020}}^{\prime}}-\overline{{w}_{a;2}^{\prime}{{\psi}_{2}^{\u2020}}^{\prime}})\phantom{\rule{1.em}{0ex}}& +{H}_{2}\left(\overline{{w}_{a;1}^{\prime}\frac{{\psi}_{1}^{\prime}-{\psi}_{2}^{\prime}}{{H}_{2}+{H}_{1}}}+\overline{{w}_{a;2}^{\prime}\frac{{\psi}_{2}^{\prime}-{\psi}_{3}^{\prime}}{{H}_{3}+{H}_{2}}}\right)]\hfill \\ \hfill \phantom{\rule{1.em}{0ex}}& -{H}_{2}\overline{{\psi}_{2}}{\nabla}_{\mathrm{h}}\xb7\overline{{\mathbf{u}}_{g;2}^{\prime}{\nabla}_{\mathrm{h}}^{2}{\psi}_{2}^{\prime}}+{H}_{2}\overline{{\psi}_{2}^{\prime}{\nu}_{4}{\nabla}_{\mathrm{h}}^{4}({\nabla}_{\mathrm{h}}^{2}{\psi}_{2}^{\prime})},\hfill \end{array}$$

$${\partial}_{t}{P}_{2}^{\#}=-{f}_{0}\overline{{w}_{a;2}}(\overline{{\psi}_{2}}-\overline{{\psi}_{3}})-\frac{{f}_{0}^{2}}{{g}_{2}^{\prime}}(\overline{{\psi}_{2}}-\overline{{\psi}_{3}}){\nabla}_{\mathrm{h}}\xb7\overline{{\mathbf{u}}_{{g;2}^{\prime}}^{\u2020}({\psi}_{2}^{\prime}-{\psi}_{3}^{\prime})}-\frac{{f}_{0}^{2}}{{g}_{2}^{\prime}}(\overline{{\psi}_{2}}-\overline{{\psi}_{3}}){\kappa}_{4}{\nabla}_{\mathrm{h}}^{4}(\overline{{\psi}_{2}}-\overline{{\psi}_{3}}),$$

$${\partial}_{t}{\mathcal{P}}_{2}=-{f}_{0}\overline{{w}_{a;2}^{\prime}({\psi}_{2}^{\prime}-{\psi}_{3}^{\prime})}+\frac{{f}_{0}^{2}}{{g}_{2}^{\prime}}(\overline{{\psi}_{2}}-\overline{{\psi}_{3}}){\nabla}_{\mathrm{h}}\xb7\overline{{\mathbf{u}}_{{g;2}^{\prime}}^{\u2020}({\psi}_{2}^{\prime}-{\psi}_{3}^{\prime})}-\frac{{f}_{0}^{2}}{{g}_{2}^{\prime}}\overline{({\psi}_{2}^{\prime}-{\psi}_{3}^{\prime}){\kappa}_{4}{\nabla}_{\mathrm{h}}^{4}({\psi}_{2}^{\prime}-{\psi}_{3}^{\prime})}.$$

For the third layer:

$$\begin{array}{cc}\hfill {\partial}_{t}{K}_{3}^{\#}={f}_{0}\left[-\overline{{w}_{a;2}}\overline{{\psi}_{2}^{\u2020}}+{H}_{3}\overline{{w}_{a;2}}\frac{\overline{{\psi}_{2}}-\overline{{\psi}_{3}}}{{H}_{3}+{H}_{2}}\right]& +{H}_{3}\overline{{\psi}_{3}}{\nabla}_{\mathrm{h}}\xb7\overline{{\mathbf{u}}_{g;3}^{\prime}{\nabla}_{\mathrm{h}}^{2}{\psi}_{3}^{\prime}}\hfill \\ \hfill \phantom{\rule{1.em}{0ex}}& +{H}_{3}\overline{{\psi}_{3}}[{\nu}_{4}{\nabla}_{\mathrm{h}}^{4}({\nabla}_{\mathrm{h}}^{2}\overline{{\psi}_{3}})+\u03f5{\nabla}_{\mathrm{h}}^{2}\overline{{\psi}_{3}}],\hfill \end{array}$$

$$\begin{array}{cc}\hfill {\partial}_{t}{\mathcal{K}}_{3}={f}_{0}\left[-\overline{{w}_{a;2}^{\prime}{{\psi}_{2}^{\u2020}}^{\prime}}+{H}_{3}\overline{{w}_{a;2}^{\prime}\frac{{\psi}_{2}^{\prime}-{\psi}_{3}^{\prime}}{{H}_{3}+{H}_{2}}}\right]& -{H}_{3}\overline{{\psi}_{3}}{\nabla}_{\mathrm{h}}\xb7\overline{{\mathbf{u}}_{g;3}^{\prime}{\nabla}_{\mathrm{h}}^{2}{\psi}_{3}^{\prime}}\hfill \\ \hfill \phantom{\rule{1.em}{0ex}}& +{H}_{3}\overline{{\psi}_{3}^{\prime}[{\nu}_{4}{\nabla}_{\mathrm{h}}^{4}({\nabla}_{\mathrm{h}}^{2}{\psi}_{3}^{\prime})+\u03f5{\nabla}_{\mathrm{h}}^{2}{\psi}_{3}^{\prime}]}.\hfill \end{array}$$

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].

- Phillips, N.A. The general circulation of the atmosphere: A numerical experiment. Q. J. R. Meteorol. Soc.
**1956**, 82, 123–164. [Google Scholar] [CrossRef] - Lorenz, E.N. The Nature and Theory of the General Circulation of the Atmosphere; World Meteorological Organization: Geneva, Switzerland, 1967; Volume 218. [Google Scholar]
- 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. [Google Scholar] [CrossRef] - Holland, W.R.; Schmitz, W.J. Zonal penetration scale of model midlatitude jets. J. Phys. Oceanogr.
**1985**, 15, 1859–1875. [Google Scholar] [CrossRef] - Kang, D.; Curchitser, E.N. Energetics of Eddy–Mean Flow Interactions in the Gulf Stream Region. J. Phys. Oceanogr.
**2015**, 45, 1103–1120. [Google Scholar] [CrossRef] - McWilliams, J.C. A perspective on the legacy of Edward Lorenz. Earth Space Sci.
**2019**, 6, 336–350. [Google Scholar] [CrossRef] - Vallis, G.K. Atmospheric and Oceanic Fluid Dynamics, 2nd ed.; Cambridge University Press: Cambridge, UK, 2017. [Google Scholar]
- 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. [Google Scholar] [CrossRef] - 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. [Google Scholar] [CrossRef] - Lévy, M.; Franks, P.J.; Smith, K.S. The role of submesoscale currents in structuring marine ecosystems. Nat. Commun.
**2018**, 9, 1–16. [Google Scholar] [CrossRef] - Thomas, L.N.; Tandon, A.; Mahadevan, A. Submesoscale processes and dynamics. Ocean Model. Eddying Regime
**2008**, 177, 17–38. [Google Scholar] [CrossRef] - McWilliams, J.C. Submesoscale currents in the ocean. Proc. R. Soc. A Math. Phys. Eng. Sci.
**2016**, 472, 20160117. [Google Scholar] [CrossRef] [PubMed] - 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. [Google Scholar] [CrossRef] - 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. [Google Scholar] [CrossRef] - 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. [Google Scholar] [CrossRef] - Tagklis, F.; Bracco, A.; Ito, T.; Castelao, R. Submesoscale modulation of deep water formation in the Labrador Sea. Sci. Rep.
**2020**, 10, 1–13. [Google Scholar] [CrossRef] - Charney, J.G. Geostrophic turbulence. J. Atmos. Sci.
**1971**, 28, 1087–1095. [Google Scholar] [CrossRef] - 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. [Google Scholar] [CrossRef] - 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. [Google Scholar] [CrossRef] - 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**. [Google Scholar] [CrossRef] - 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. [Google Scholar] - 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. [Google Scholar] [CrossRef] - Eden, C.; Böning, C. Sources of eddy kinetic energy in the Labrador Sea. J. Phys. Oceanogr.
**2002**, 32, 3346–3363. [Google Scholar] [CrossRef] - Pujol, M.I.; Larnicol, G. Mediterranean sea eddy kinetic energy variability from 11 years of altimetric data. J. Mar. Syst.
**2005**, 58, 121–142. [Google Scholar] [CrossRef] - 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. [Google Scholar] [CrossRef] - 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. [Google Scholar] [CrossRef] - 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. [Google Scholar] [CrossRef] - Uchida, T.; Abernathey, R.; Smith, S. Seasonality of eddy kinetic energy in an eddy permitting global climate model. Ocean Model.
**2017**, 118, 41–58. [Google Scholar] [CrossRef] - Dong, J.; Fox-Kemper, B.; Zhang, H.; Dong, C. The seasonality of submesoscale energy production, content, and cascade. Geophys. Res. Lett.
**2020**. [Google Scholar] [CrossRef] - 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**. [Google Scholar] [CrossRef] - 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**. [Google Scholar] [CrossRef] - Fox-Kemper, B.; Ferrari, R.; Hallberg, R. Parameterization of mixed layer eddies. Part I: Theory and diagnosis. J. Phys. Oceanogr.
**2008**, 38, 1145–1165. [Google Scholar] [CrossRef] - Johnson, L.; Lee, C.M.; D’Asaro, E.A. Global Estimates of Lateral Springtime Restratification. J. Phys. Oceanogr.
**2016**, 46, 1555–1573. [Google Scholar] [CrossRef] - 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. [Google Scholar] [CrossRef] - 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. [Google Scholar] [CrossRef] - Pedlosky, J. Geophysical Fluid Dynamics, 2nd ed.; Springer: New York, NY, USA, 1987. [Google Scholar]
- 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. [Google Scholar] [CrossRef] - Young, W.R. An Exact Thickness-Weighted Average Formulation of the Boussinesq Equations. J. Phys. Oceanogr.
**2012**, 42, 692–707. [Google Scholar] [CrossRef] - Deremble, B.; Martinez, E.M. MSOM: Multiple Scale Ocean Model; MEOM Research Group: Grenoble, France, 2020. [Google Scholar] [CrossRef]
- Popinet, S. A quadtree-adaptive multigrid solver for the Serre-Green-Naghdi equations. J. Comput. Phys.
**2015**, 302, 336–358. [Google Scholar] [CrossRef] - Verron, J. Nudging satellite altimeter data into quasi-geostrophic ocean models. J. Geophys. Res. Ocean.
**1992**, 97, 7479–7491. [Google Scholar] [CrossRef] - 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. [Google Scholar] [CrossRef] - 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. [Google Scholar] [CrossRef] - Karabasov, S.; Berloff, P.S.; Goloviznin, V. CABARET in the ocean gyres. Ocean Model.
**2009**, 30, 155–168. [Google Scholar] [CrossRef] - Hallberg, R. Using a resolution function to regulate parameterizations of oceanic mesoscale eddy effects. Ocean Model.
**2013**, 72, 92–103. [Google Scholar] [CrossRef] - Rhines, P.B.; Schopp, R. The wind-driven circulation: Quasi-geostrophic simulations and theory for nonsymmetric winds. J. Phys. Oceanogr.
**1991**, 21, 1438–1469. [Google Scholar] [CrossRef] - 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. [Google Scholar] [CrossRef] - Simonnet, E. Quantization of the low-frequency variability of the double-gyre circulation. J. Phys. Oceanogr.
**2005**, 35, 2268–2290. [Google Scholar] [CrossRef] - 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. [Google Scholar] [CrossRef] - 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. [Google Scholar] [CrossRef] - 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. [Google Scholar] [CrossRef] - Hoskins, B.J.; Draghici, I.; Davies, H.C. A new look at the ω-equation. Q. J. R. Meteorol. Soc.
**1978**, 104, 31–38. [Google Scholar] [CrossRef] - Gill, A.E. Atmosphere-Ocean Dynamics; Academic Press: Cambridge, MA, USA, 1982; Volume 30, 662p. [Google Scholar]
- Grooms, I.; Julien, K.; Fox-Kemper, B. On the interactions between planetary geostrophy and mesoscale eddies. Dyn. Atmos. Ocean.
**2011**, 51, 109–136. [Google Scholar] [CrossRef] - 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. [Google Scholar] [CrossRef] - 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. [Google Scholar] [CrossRef] - 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. [Google Scholar] [CrossRef] - 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. [Google Scholar] [CrossRef] - Chassignet, E.P.; Marshall, D.P. Gulf Stream separation in numerical ocean models. Geophys. Monogr. Ser.
**2008**, 177. [Google Scholar] [CrossRef] - 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). - Griffies, S.M. The Gent–McWilliams skew flux. J. Phys. Oceanogr.
**1998**, 28, 831–841. [Google Scholar] [CrossRef]

Parameter | Notation | Value | Unit |
---|---|---|---|

Number of horizontal grids | N | 1024 | - |

Number of vertical layers | ${n}_{l}$ | 3 | - |

Non-dim. horizontal domain size | ${\widehat{L}}_{0}$ | 80 | - |

Non-dim. horizontal resolution | ${\delta}_{\widehat{x}}$ | ${N}^{-1}{\widehat{L}}_{0}$ | - |

Background Rossby number | $R{o}^{m}$ | $0.025$ | - |

Non-dim. Coriolis parameter | ${\widehat{f}}_{0}$ | ${R{o}^{m}}^{-1}$ | - |

Bottom Ekman number | $E{k}^{b}$ | $0.004$ | - |

Non-dim. surface Ekman pumping | ${\widehat{\tau}}_{0}$ | $0.0001$ | - |

Biharmonic Reynolds number | $R{e}_{4}$ | 4000 | - |

Non-dim. beta | $\widehat{\beta}$ | $0.5$ | - |

Background Froude number | $F{r}_{1}^{m};F{r}_{2}^{m}$ | $0.00409959;0.01319355$ | - |

Amplitude of $F{r}_{i}$ | ${\widehat{A}}_{F{r}_{1}};{\widehat{A}}_{F{r}_{2}}$ | $0.1;0$ | - |

Non-dim. frequency of $F{r}_{i}$ | ${\widehat{f}}_{F{r}_{1}};{\widehat{f}}_{F{r}_{2}}$ | $62.{2}^{-1};62.{2}^{-1}$ | - |

Non-dim. layer thickness | ${\widehat{H}}_{1};{\widehat{H}}_{2};{\widehat{H}}_{3}$ | $0.06;0.14;0.8$ | - |

Non-dim. reduced gravity | $\widehat{{g}_{i}^{\prime}}$ | ${F{r}_{i}}^{-2}\widehat{{H}_{i}^{\u2020}}$ | - |

Non-dim. maximum time stepping | ${\delta}_{\widehat{t}}^{\mathrm{max}}$ | $5\times {10}^{-2}$ | - |

CFL condition | - | $0.4$ | - |

Horizontal velocity | U | $0.1$ | [m s${}^{-1}$] |

Length scale | L | 50 | [km] |

Total layer thickness | H | 5000 | [m] |

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

© 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/).