A Geometric Perspective on the Modulation of Potential Energy Release by a Lateral Potential Vorticity Gradient

: The release of available potential energy by growing baroclinic instability requires the slope of the eddy ﬂuxes to be shallower than that of mean density surfaces, where the amount of energy released depends on both the ﬂux angle and the distance of ﬂuid parcel excursions against the background density gradient. The presence of a lateral potential vorticity (PV) gradient is known to affect the growth rate and energy release by baroclinic instability, but often makes the mathematics of formal linear stability analysis intractable. Here the effects of a lateral PV gradient on baroclinic growth are examined by considering its effects on the slope of the eddy ﬂuxes. It is shown that the PV gradient systematically shifts the unstable modes toward higher wavenumbers and creates a cutoff to the instability at large scales, both of which steepen the eddy ﬂux angle and limit the amount of energy released. This effect may contribute to the severe inhibition of baroclinic turbulence in systems dominated by barotropic jets, making them less likely to transition to turbulence-dominated ﬂow regimes.


Introduction
Lateral gradients of potential vorticity (PV) are responsible for several of the most remarkable features studied in the field of geophysical fluid dynamics. This "β-effect" (so named because it often refers to the meridional gradient of planetary vorticity via a linear approximation to the Coriolis parameter) is in many ways paradoxical, in that it implicates many fundamental concepts from different branches of fluid dynamics in a somewhat opaque fashion. The effect is associated with the development of laterally sheared jets which simultaneously suppress cross-stream mixing (e.g., [1][2][3]) while enhancing the mixing in the streamwise direction (e.g., [4]), establishes fundamental length scales where both waves and turbulence are of leading-order dynamical importance (e.g., [5]), imparts anisotropy into the inverse cascade of kinetic energy (e.g., [6]), and reinforces its own flow structures through upgradient momentum fluxes (e.g., [7]). These dynamics are most readily observed on gas planets (e.g., [8,9]), but also have terrestrial analogs (e.g., [10,11]) which can play important roles in weather and climate [12,13].
The β-effect is known to exert a strong influence on the energetics of chaotic flows, in both the transition from linear instabilities to nonlinear turbulence (e.g., [14]) and in the "arrest" phase where the upscale energy cascade is either steered into predominantly zonal motions (e.g., [15]) or dissipated by large-scale friction (e.g., [16]). The zonation caused by β results in the appearance of persistent coherent structures that can emerge even without an upscale cascade [17][18][19] and can be studied based on their statistical equilibria (e.g., [20][21][22]). The most well-recognized of these coherent structures are laterally sheared jets that can suppress or completely eliminate baroclinic instabilities (hereafter BCI; [23][24][25][26]). Recent studies of BCI suppression (e.g., [14,27]) noted that the suppression process tends to occur in distinct stages, where linear instability first begets the onset of baroclinic turbulence, followed by evolution of the flow's shear and stratification that stabilizes the flow and prevent further significant changes. References [24,28] attribute the suppression partly to meridional confinement of normal modes of instability, wherein the PV gradient restricts the horizontal displacement of fluid parcels and thus limits their ability to tap available potential energy. The suppression is also partly due to a positive feedback loop in which upgradient momentum fluxes reinforce the jet, both enhancing its PV gradient and distorting the normal modes so that they become less coherent. The role of β in this process is to generate a horizontally convergent momentum flux that begins the process of jet reinforcement [28].
Though much interest in β-plane dynamics has justifiably tended toward these later-stage, nonlinear effects, PV gradients were also shown to have an effect on the linear stability of flows (e.g., [14,[29][30][31]). Linear stability has been shown to be key for understanding the tendency of β-plane jets to be organized into flow states that are marginally stable to baroclinic instability ("baroclinically adjusted", e.g., [14,32]), and may also influence how flows transition from turbulent-to jet-dominated states (e.g., [33]). It was recognized in [34] that the linear instability of baroclinic flows can be attributed mechanistically to the release of available potential energy when fluid parcels follow trajectories that are shallower than density surfaces, where the amount of potential energy released depends on the angle the trajectory follows [29,35]. This simple geometric description of baroclinic instability has become a commonly used tool for teaching the basics of instabilities and turbulence (cf. [31]).
The purpose of this study is to examine the effect of β on baroclinic instability using the above geometric interpretation; that is, how β affects the angle of fluid parcel trajectories for growing baroclinic modes. A combination of theory and numerical modeling will be used to demonstrate that β suppresses baroclinic growth by forcing parcels to follow trajectories that are suboptimal for potential energy release. The article concludes with perspectives on these results with regard to the behavior of more realistic flows, and opportunities for future work.

Theory
Consider a Boussinesq fluid under the quasigeostrophic (QG) approximation, and for simplicity assume a doubly periodic domain in the horizontal bounded by rigid surfaces in the vertical. The set of state variables that will be sufficient for this study include the horizontal geostrophic velocity, u = (u, v), the ageostrophic vertical velocity, w, buoyancy perturbation, b, and QG potential vorticity, q (QGPV). The Coriolis parameter, f = f 0 + βy, will consist of a fixed reference value, f 0 , and meridionally varying component with β = ∂ f /∂y > 0. The buoyancy frequency associated with the Boussinesq background state will be denoted N 2 , and will be treated as a positive constant.
Following Taylor and Ferrari [36], the buoyancy and velocity fields will include a background state which is assumed to be in thermal wind balance, such that the total buoyancy and velocity are given by and are related via a constant lateral buoyancy gradient, M 2 . The prognostic equations will solve for the perturbations away from this background state. Of primary interest in this study are the budgets for the buoyancy and QGPV, which for simplicity will use no diffusivity or viscosity, and where the QGPV is related to the state variables via the QG streamfunction, ψ: Here D/Dt = ∂ t + u T · ∇ h is the material derivative using the horizontal gradient operator ∇ h = ∂ x , ∂ y , which will be distinct from the three dimensional gradient operator, ∇ = ∂ x , ∂ y , ∂z , used further on. With the doubly periodic domain a Reynolds averaging will be employed, indicated by (·), such that the overbar indicates a horizontal average in both directions and primes indicate deviations from this average. Please note that with this flow and domain configuration the case where β = 0 essentially constitutes the classic Eady problem [34]. The emergent dynamics when β > 0 will be the focus of this study.
Consider a scenario of growing baroclinic instabilities, which are marked by a temporal increase in the magnitude of the eddy available potential energy, b 2 /2N 2 (EPE), and eddy potential enstrophy, q 2 /2. The budgets for these variables can be obtained through multiplication and averaging of Equations (3) and (4), and are given by 1 2 Please note that the periodicity of the domain and horizontal averaging has allowed the advection and triple correlation terms to be eliminated from these budgets. We first observe that in order for the eddy potential enstrophy to grow in time we must have a QGPV flux that is directed down the background QGPV gradient, which for β > 0 implies that v q < 0. Furthermore, assuming w b is positive (indicating an increase in the eddy kinetic energy, whose budget is not shown here), in order for the EPE to grow in time we must have v b M 2 < 0, so that the horizontal eddy buoyancy flux is also downgradient. Finally, we must have The parameter C b describes the slope of the eddy buoyancy fluxes, which for growing instabilities must fall between the horizontal (C b = 0) and the mean isopycnal slope, These bounds on C b establish the concept of the "wedge of instability" (e.g., [29,34,35,37]), inside which the buoyancy force performs positive work on the parcel and tends to vertically accelerate the parcel away from its initial position, increasing the growth of the instability. The seminal work of Eady [34] established the well-known result that the EPE release is maximized when C b = 1/2, though this value must be interpreted with caution because it is based on a heuristic argument concerning the displacement of individual fluid parcels rather than the fluid as a whole. Note also that C b must be zero at the vertical boundaries in order to satisfy the zero flux boundary conditions, implying that C b varies with z.
The concept of the eddy flux slope optimizing (or limiting) potential energy release is a useful tool for teaching the mechanics of baroclinic instability and its energetics. Here the intent is to take this concept a step further to explore how β modulates this energy release through its influence on C b . To this end, note that one can solve for C b directly by substituting (9) into (7), which leads to This expression is consistent with the assumption laid out in (9)-assuming a growing mode the numerator of the second term on the right is positive and (as mentioned previously) its denominator is negative, so that their ratio is negative. Thus, the role of this term is to decrease C b below one, so that the buoyancy flux slope remains within the wedge of instability.
A lesser-known budget can be obtained by multiplying (3) by q and (4) by b, and after application of the Reynolds averaging one arrives at a prognostic equation for the eddy covariance of the buoyancy and QGPV, In similar fashion to (9), one may define the slope of the eddy QGPV flux in terms of the mean isopycnal slope, which can in turn be substituted back into (11) to obtain Please note that there is no a priori expectation that the tendency of q b must be positive (e.g., q and b could be anticorrelated), or that C q would have the same bounds as C b .
The significance of C q can be seen by using thermal wind balance and leveraging the double periodicity and Reynolds averaging to derive a relationship between the slope of the eddy QGPV and buoyancy fluxes, where the term on the right involves the dot product of the three-dimensional gradient of w with the QGPV induction vector, Assuming v b varies rapidly in the vertical compared to C b , one may substitute (9) and (12) into this expression to obtain This solution for C q may then be substituted into (13), yielding The terms of this equation are labeled with circled numbers for easy reference here and in Section 3. Within this expression for C b consider specifically term 3 . By the discussion following (7)-(9), both the numerator and denominator of this term are negative for growing instabilities, implying that the full term is guaranteed to be positive as long as β > 0. The effect of this term is thus to always increase the slope of the eddy buoyancy flux. Assuming that C b tends toward the optimal value of 1/2 for the β = 0 case, the hypothesis pursued here is that the effect of β is thus to drive the slope of the buoyancy fluxes steeper than their optimal value, limiting the generation of EPE. This hypothesis is tested, and the role of each term in (16) evaluated, using the numerical experiments discussed next.

Numerical Simulations and Results
The spectral flow solver Dedalus [39] was used to create a 20-layer quasigeostrophic (QG) model with which to test the behavior of C b and the growth of eddy energy due to changes in β. Since the focus here is on the growth phase of baroclinic instabilities, the linearized QG potential vorticity (QGPV) equations are solved within each layer, which prevents both nonlinear energy transfers between wavenumbers and the formation of jets. To be consistent with (7) and (8), the equations are inviscid and employ no bottom drag. The model is doubly periodic in the horizontal directions, with each layer containing 512 × 512 gridpoints and featuring a deformation radius of L d = 4.5 km. The primary advantage of this domain geometry is that the Dedalus solver allows the QGPV equation for each layer to be defined independently from one another, so that no explicit discretization or spectral basis for the vertical direction is required (the layers are coupled via the vortex stretching terms by defining algebraic equations for q corresponding to each layer-see the Supplemental Materials). The grid spacing is set to be ∆x = 2.5 L d , and the background velocity in each layer n is U n = 0.025 − (n − 1)/19 m s −1 . The strength of the background PV gradient is measured by defining the nondimensional variable β * = βL 2 d /U 1 . Five simulations were run, spanning choices of β * = {0, 0.5, 1.0, 1.5, 2.0}. The simulations were initially seeded with random noise at all wavenumbers in the QGPV field in each layer. Only a relatively small subset of the wavenumbers resolved by the grid were linearly unstable and exhibited exponential growth (e.g., [31]), and these modes quickly became distinguishable from the ambient noise. The models were run for 20,000 timesteps, with an initial timestep ∆t = 3600 s. Since Dedalus uses an adaptive timestep, the actual simulated time for each run varied between 1600 and 1800 days, or approximately 40 to 45 times the Eady growth timescale. Snapshots of the model output were taken every 100 timesteps (for 200 snapshots in total). At each output time t the linear growth rate for each zonal wavenumber k was measured as which essentially measures the ratio of the eddy kinetic energies at successive snapshots (for details see [14], Equations (11) and (12)). Because there is no dissipation or bottom drag the fields grow exponentially until the simulation is terminated, and would continue to grow in this way if the simulations were not stopped. Thus, the simulations were stopped after the diagnostics exhibited a steady growth rate, which comprised approximately the final 150 snapshots. The model output was averaged over the final 100 snapshots to produce the results shown here. To be consistent with the averaging operation defined previously, the values of each variable at each gridpoint are treated as "prime" quantities, which are multiplied and averaged over the full horizontal domain to produce the mean eddy covariances. The first and last layers are ignored in this analysis to avoid the QGPV gradient associated with the temperature gradient at the boundaries (e.g., [31]), and one more layer is sacrificed from the output in order to calculate b via the vertical derivative of ψ . Other layerwise quantities such as q and u are averaged to the layer interfaces so that they are colocated with b . Thus, only output from the remaining seventeen interior layers is considered. Figure 1 shows the time-averaged values of C b (panel a), C q (panel b), and the terms comprising the solution for C b in (16) (panel c) as a function of the model layer. Please note that values for the β * = 0 case (blue line) are only shown in panel (a), since there is negligible potential enstrophy generation and the theory for C q becomes invalid with no background QGPV gradient. The average of each curve over all layers is shown using the colored dots on the left side of each panel. It is immediately clear that the slope of the eddy buoyancy fluxes varies significantly across layers, changing by around a factor of two for C b and up to a factor of three for C q at larger values of β * . This supports the aforementioned cautionary note about using parcel excursion theory to interpret the behavior of the full fluid continuum, as the layerwise values of C b in the β * = 0 case deviate significantly from the predicted optimal value of 1/2 even though the domain average is relatively close (0.4, blue dot). Both C b and C q decrease toward the vertical boundaries because the vertical eddy fluxes must tend to zero to satisfy the no-flux boundary conditions, while the horizontal fluxes do not. It is also evident that increasing β * serves to increase the mean value of both C b and C q , consistent with the hypothesis following (16). Of more significant interest is the behavior of the terms contributing to C b in (16), namely how each of them is affected by changes in β * (Figure 1c). In this panel both terms 2 and 3 are calculated directly from the model output, and term 1 is shown to be the residual from subtracting the other terms from the calculation of C b . The gross behavior of these terms is consistent with the predictions from theory; term 3 is positive definite and has mean values which increase with larger β * , and term 2 must be negative to offset the positive values from 3 to keep the eddy flux slope within the wedge of instability. Term 1 also provides a negative-definite contribution to this budget, whose magnitude increases with larger β * . Perhaps the most interesting feature of these terms is that increases in C b with increasing β * are not driven only by the explicit dependence of β * in 3 , which increases relatively little as β * grows. Rather, C b is also affected by a decrease in the magnitude of 2 , which is partially offset by an increase in the magnitude of 1 . Thus, all three labeled terms in (16) play a role in the steepening of the buoyancy fluxes as β * becomes larger.
To investigate further, each output snapshot was Fourier transformed to obtain the values of each term in as a function of zonal (k) and meridional (l) wavenumber. Since the fastest growing mode in the linear instability problem corresponds to the gravest meridional mode, only the time-averaged results for l = 0 are shown here. Figure 2a shows the diagnosed growth rates from layer 10, following the method described in (17) for calculating σ t . The curve for the β * = 0 case, corresponding to the Eady model, matches the theory shown by [31], where the short-wave cutoff occurs at the normalized (nondimensional) wavenumber µ = kL d = 2.4 and maximal nondimensional growth rate σ = 0.31 occurs at µ = 1.61. The curves for β > 0 are shifted and tail off toward higher wavenumbers, akin to solutions of the Charney problem [40,41]; however, modes smaller than µ ≈ 3 failed to demonstrate steady growth rates and featured erratic behavior of C b , and are thus excluded from this analysis.  Figure 2b shows the behavior of the terms in (16) as a function of µ. As predicted, the values of 3 (solid lines) are always positive and become greater as β * increases. These curves bend upward toward smaller wavenumbers, unlike those of 1 and 2 , which appear to be constant and to trend linearly at small wavenumbers, respectively. Conversely, at large wavenumbers the effect of 3 tends to zero, and the tendency for the slope to remain in the wedge of instability becomes a competition between 1 (negative values) and 2 (positive values). As shown in (14) the presence of the gradient operator in the numerator of 1 suggests that this term should trend linearly with µ, which appears to be approximately the case for larger wavenumbers. Interestingly, 2 transitions from negative to positive as one moves toward higher wavenumbers, with the transition point occurring near the short-wave cutoff for the Eady problem. Investigation of these behaviors at large wavenumbers are beyond the scope of this study, but are certainly worthy of further consideration.
Panel (c) shows values of C b as a function of µ, where the slope of the unstable modes essentially overlap each other regardless of the value of β * except at large wavenumbers. This suggests that near the deformation radius the slope of the modes becomes set by their physical size (i.e., distance of their meridional oscillations) rather than by β * . In the Eady case (blue line) where β * imposes no constraint, the slope of the modes transitions smoothly from the isopycnal slope (C b = 1) at small scales to horizontal (C b = 0) at the domain scale. The gray vertical dashed lines indicate the most unstable Eady mode at µ = 1.6, at which conversion of EPE to eddy kinetic energy is maximized. This corresponds to the optimal buoyancy flux slope of C b = 0.5 (horizontal dashed line), affirming that energy conversion is maximized along the direction of the most unstable linear mode (e.g., [35]). This behavior stands in contrast to the cases β * > 0, where because β * imposes a large-scale cutoff (panel a) the curves terminate at large scales with C b > 0. Because the growth curves are shifted toward larger wavenumbers, the most unstable linear modes no longer correspond to the optimal slope C b = 0.5; rather, they occur at steeper slopes where potential energy conversion is less efficient.
Finally, these curves can be interpreted in the context of Figure 1a, which showed the domain-averaged values of C b calculated in physical space. The domain averaging essentially represents an average over all wavenumbers, but is weighted toward the more energetic, faster-growing modes. To this end, Figure 2d shows values of C b in panel (c) weighted by the growth rates in panel (a), and normalized so that the β * = 2.0 curve has a peak value of one. Here it is evident that the tendency for the growth curves to shift toward higher wavenumbers at larger β * means that the domain-averaged C b becomes weighted toward steeper modes.
In summary, the tendency for C b to increase with β * occurs because a nonzero β imposes a threshold wavenumber above which the flow remains linearly stable. This prevents larger, shallower modes from contributing to the domain-averaged value of C b , and restricts instability to smaller waves which are less efficient at generating EPE (e.g., Appendix B in [42]). These small waves also are oriented at steeper slopes, further inhibiting EPE generation. A schematic depiction of this behavior is shown in Figure 2e. For a given vertical fluid parcel excursion δz, increasing β restricts the horizontal distance δy the parcel can travel. The overall parcel trajectory thus becomes steeper, occurring along paths that are suboptimal for EPE generation. Even if the growth rates are not reduced significantly (panel a), limiting δy means that less EPE is generated overall irrespective of the trajectory. Together these effects can severely inhibit EPE generation (and thus eddy formation), possibly helping to explain observed "threshold behavior" (e.g., [33]) where the flow rapidly switches from eddy-dominated to jet-dominated states despite only small increases in the effective β.

Conclusions
A useful mnemonic for teaching the concept of slope is "rise over run", with the "rise" referring to the change in vertical position and the "run" referring to the change in horizontal position. The same mnemonic can be used to understand the effect of β in inhibiting potential energy release by baroclinic instability. Two principle restoring forces exist in geophysical fluid dynamics: a vertical restoring force set by the density stratification that controls the "rise", and a horizontal restoring force set by the potential vorticity gradient (a "half-cycle Coriolis force", e.g., [43]) that controls the "run". Thus, when β is increased and the "run" is reduced, the slope of the eddy fluxes is increased.
The slope of the eddy fluxes depends on the horizontal wavelength, with larger waves naturally having shallower slopes and shorter waves steeper slopes. When β = 0 the slopes within the range of unstable wavenumbers smoothly transition from the isopycnal slope at the short-wave cutoff to zero at the domain scale (Figure 2a, blue line). Thus, the set of unstable modes occupies the full extent of the wedge of instability. Nonzero β restricts the unstable modes to a progressively smaller portion of this wedge; by imposing a long-wave cutoff, the eddy transport is forced to become steeper, inhibiting its ability to extract potential energy both by forcing it along suboptimal slopes and by limiting the distance of its horizontal excursions. This manifests in reduced linear growth rates and would likely imply a reduced ability to overcome the barotropic governor (e.g., [24]) and manifest baroclinic turbulence from a marginally stable barotropic jet.
This study focused on a relatively narrow portion of the baroclinic life cycle, namely the linear growth phase in the limit of constant β. An interesting extension of these results might be found by exploring a flow configuration similar to that of [14], who found that linear growth rates increased with the strength of a spatially varying β. Much of the theoretical framework introduced here would still apply when β is nonconstant or even negative, albeit any expectations about its effects would have to be made depending on the specific flow configuration. It would be interesting to explore whether this theory has applications for flows where sloping bathymetry introduces an "effective β" and affects the instability of the baroclinic modes (e.g., [44]). It would likely be possible to apply some of the theory here to nonlinear or baroclinically adjusted flows as well (i.e., including jets or lateral shears) with appropriate modifications to the averaging, if necessary. These topics would make for interesting new perspectives on longstanding knowledge about β-plane dynamics.