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

With the advent of submesoscale O(1 km) permitting basin-scale ocean simulations, the seasonality of mesoscale O(50 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) 1, O(Fr) 1) while submesoscale dynamics are associated with O(Ro) ∼ 1, O(Fr) & 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.


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][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 ∼ O (50 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 0 L)) 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) 1, O(Fr) 1, and submesoscale flows being associated with O(Ro) ∼ 1, O(Fr) 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 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.

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 km) and horizontal resolution is ∼4 km (=δ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) 1, O(Fr) 1). MSQG solves prognostically for the non-dimensionalized QG potential vorticity (PV): where q = ζ g + βy − f 0 H h is the QGPV (details are given in Appendix A; [7]) and the β-plane approximation is applied ( f = f 0 + βy). F and D are the forcing and dissipative terms, and (·) 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 ∇ h is the horizontal gradient operator, andŷ (∈ [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 (τ) 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: 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 = (H i + H i+1 )/2, g is the reduced gravity and subscript i is the layer index ( Figure 1). We vary Fr 1 in time but keep Fr 2 stationary (Â 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]. We spin up the model for 10 years from a spun-up run with lower resolution (N = 256, equivalently δx L ∼ 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 × 10 5 seconds ∼5.8 days).

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 (∈ [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., ζ = ζ g + ζ a ). We denote the partial derivatives as ∂ (·) with respect to (t, z, y, x). The stream function is defined as ψ i = φ g;i / f 0 where φ g;i is the geostrophic pressure anomaly and relative vorticity can be written as ζ g;i = ∇ 2 h ψ 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. The ageostrophic vertical velocity can be diagnosed via the QG omega equation (Appendix B; [52,53]): where is the buoyancy. The Q tensor is: ; [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 ): where δ E = Ek b H 3 is the bottom Ekman-layer thickness. Now, multiplying Equation (5) by −ψ 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: On the other hand, using relation (A2), the layer-thickness equations can be manipu- where the second term on the right-hand side above (15) vanishes due to thermal wind. where The available potential energy (APE) equations can, therefore, be derived by multiplying Equation (15) with f 0 (ψ 1 − ψ 2 ) and again dropping the divergence terms: and Equation (16) 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), mean APE (MAPE; P # ) and eddy APE (EAPE; P) can be defined as: where (·) is the ensemble mean and the eddy is defined as fluctuations about the ensemble mean, viz. (·) = (·) − (·). We note that the ensemble mean of the fluctuations vanish ((·) = 0). The strength of defining the mean as such is that in addition to the ensemblemean 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 Hence, the exchanges (Π) of KE and APE within and between layers are: where · = (·)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:

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 ( |∇ h Ψ| 2 where Ψ = H −1 ∑ i H i ψ 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 , 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).

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 (Π P # →K # < 0; black solid line in Figure 4b) due to Ekman pumping steepening the isopycnals. The energy input due to wind stress (F K # s ) 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 (Π P →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 (Π P →K < 0), whereas baroclinic instability would predict EAPE→EKE (Π P →K > 0). Regarding the dissipation terms, only the bottom drag for EKE (D K b ) shows a notable seasonality and has a similar magnitude to the energy flux from MKE to EKE (Π K # →K ). The amplitude of bottom drag (|D K b |) 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 (Π), 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 (Π K # →K ) is highest during the summer (Figures 4b and 5a). We attribute the larger energy conversion from EAPE to EKE during the winter (Π P →K ) to the flow being more susceptible to baroclinic instability with reduced stratification.

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 (Π K # 2,3 →K 2,3 ). The vertical transfer of EKE (Π K 1 →K 2 , Π K 2 →K 3 ) and conversion from EAPE (Π P →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. a) b) Figure 6. Time series of volume-integrated EKE over each layer, and fluxes within and between layers (Π K 1 →K 2 , Π K 2 →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 (Π K # 2,3 →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.

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][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 (β −1 ∇ h × τ) 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.
and we get the governing equation for QGPV q i = ζ g;i + βy − 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 = U 2 H † 1 Fr −2 1 (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 The layer thickness can, therefore, be written using the stream function as [7]: where D 1 Dt η 0 = D 3 Dt η 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]: where b = B − B 0 is the QG fluctuation about the background buoyancy (B 0 ). Hence, and taking the material derivative gives the buoyancy equation in the continuously stratified framework:  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: ∂ t v g;i + u g;i ∂ x v g;i + v g;i ∂ y v g;i + f 0 u a;i + βyu g;i = −∂ y φ a;i , and multiplying them by f 0 gives: and the horizontal gradients of the buoyancy Equation (A5) with the diffusive term yields: Summing Equation (A8) with (A11), and −(A9) with (A10), and using the thermal wind relation, we get: The viscous and diffusive terms do not appear as they cancel out due to the thermalwind relation ( f 0 ∂ z ζ g = ∇ 2 h b) and their parameters being set identical (viz. ν 4 = κ 4 (= Re −1 4 L 3 U)). Taking ∂ y (A12) + ∂ x (A13) gives the omega equation for a temporally varying background stratification: 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: where D # Dt = ∂ t + u g · ∇ h . Multiplying this by −ψ gives: and taking its ensemble mean yields the mean KE equation: On the other hand, the ensemble mean of the total KE Equation (11) is: which can be expanded as: Taking the difference between Equations (A17) and (A19) gives the eddy KE equation: Since the divergence terms vanish upon area integration, we can see the mean and eddy KE exchanging the term −ψ∇ h · u g ∇ 2 h ψ (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.