The buoyancy staircase limit in surface quasigeostrophic turbulence

Surface buoyancy gradients over a quasigeostrophic fluid permit the existence of surface-trapped Rossby waves. The interplay of these Rossby waves with surface quasigeostrophic turbulence results in latitudinally inhomogeneous mixing that, under certain conditions, culminates in a surface buoyancy staircase: a meridional buoyancy profile consisting of mixed-zones punctuated by sharp buoyancy gradients, with eastward jets centred at the sharp gradients and weaker westward flows in between. In this article, we investigate the emergence of this buoyancy staircase limit in surface quasigeostrophic turbulence and we examine the dependence of the resulting dynamics on the vertical stratification. Over decreasing stratification [$\mathrm{d}N(z)/\mathrm{d}z\leq0$, where $N(z)$ is the buoyancy frequency], we obtain flows with a longer interaction range (than in uniform stratification) and highly dispersive Rossby waves. In the staircase limit, we find straight jets that are perturbed by eastward propagating along jet waves, similar to two-dimensional barotropic $\beta$-plane turbulence. In contrast, over increasing stratification [$\mathrm{d}N(z)/\mathrm{d}z\geq0$], we obtain flows with shorter interaction range and weakly dispersive Rossby waves. In the staircase limit, we find sinuous jets with large latitudinal meanders whose shape evolves in time due to the westward propagation of weakly dispersive along jet waves. These along jet waves have larger amplitudes over increasing stratification than over decreasing stratification, and, as a result, the ratio of domain-averaged zonal to meridional speeds is two to three times smaller over increasing stratification than over decreasing stratification. Finally, we find that, for a given Rhines wavenumber, jets over increasing stratification are closer together than jets over decreasing stratification.


2
McIntyre 2008). The ultimate limit of such inhomogeneous mixing, which can be achieved for sufficiently large values of , is a potential vorticity staircase: a piecewise constant potential vorticity profile consisting well-mixed regions separated by isolated discontinuities, with eastward jets centred at the discontinuities and westward flows in between (Danilov & Gurarie 2004;Dunkerton & Scott 2008;Scott & Dritschel 2012, 2019. Analogously, a buoyancy gradient at the surface of a quasigeostrophic fluid supports the existence of surface-trapped Rossby waves that are less dispersive than their barotropic counterparts (Held et al. 1995;Lapeyre 2017). The purpose of this article is to investigate the formation of zonal jets in the presence of a background surface buoyancy gradient and to examine the realizability of surface buoyancy staircases in the surface quasigeostrophic model. Although the present study is the first to systematically investigate the emergence of surface quasigeostrophic jets, there are previous studies which make use of the uniformly stratified surface quasigeostrophic model with a background buoyancy gradient. These include Smith et al. (2002), who derive the dependence of the diffusion coefficient of a passive tracer in the presence a background buoyancy gradient. Another is Sukhatme & Smith (2009), who, in their investigation of -turbulence models with a background gradient, note that, because of the decreased interaction range, surface quasigeostrophic jets in uniform stratification should be narrower than their counterparts in the barotropic model. Finally, Lapeyre (2017) demonstrates that jets can indeed form in the uniformly stratified surface quasigeostrophic model.
We also investigate how surface quasigeostrophic jets depend on the underlying vertical stratification. Yassin & Griffies (2022) show that the vertical stratification modifies the interaction range of vortices in the surface quasigeostrophic model. Suppose we have an infinitely deep fluid governed by the time-evolution of geostrophic buoyancy anomalies at its upper boundary. Then if the stratification is decreasing [ ( ) 0, where ( ) is buoyancy frequency] towards the fluid's surface (that is, the upper boundary), then the interaction range is longer than in the uniformly stratified model and the resulting turbulence is characterized by thin buoyancy filaments -analogous to the thin vorticity filaments in two-dimensional barotropic turbulence. Conversely, if the stratification is increasing [ ( ) 0] towards the surface, then the interaction range is shorter than in uniform stratification, and the buoyancy field appears spatially diffuse and lacks thin filamentary structures. In this article, we find that the interaction range is related to Rossby wave dispersion: flows with a longer interaction range have more dispersive Rossby waves whereas flows with a shorter interaction range have less dispersive Rossby waves. One of our aims is to characterize the dependence of surface quasigeostrophic jets on the functional form of the vertical stratification.
There are two motivations behind the present work. The first is its potential relevance to the upper ocean. Buoyancy anomalies at the ocean's surface are governed by the surface quasigeostrophic model (Lapeyre & Klein 2006;LaCasce & Mahadevan 2006;Isern-Fontanet et al. 2006). Both numerical (Isern-Fontanet et al. 2008;Lapeyre 2009;Qiu et al. 2016Qiu et al. , 2020Miracca-Lage et al. 2022) as well as observational (González-Haro & Isern-Fontanet 2014) studies indicate that a significant fraction of the surface geostrophic velocity is induced by sea surface buoyancy anomalies, especially over wintertime extratropical currents. Moreover, upper ocean turbulence has been found to be anisotropic (Maximenko et al. 2005;Scott et al. 2008), with significant differences in anisotropy between major extratropical currents and other regions in the ocean (Wang et al. 2019). However, our neglect of the planetary effect, as well our assumption of vanishing interior potential vorticity, may limit the direct relevance of this study to the upper ocean.
The second motivation is that the variable stratification surface quasigeostrophic model is a simple two-dimensional model in which we can investigate how jet dynamics depend on the stratification's vertical structure. Another such model is the equivalent barotropic 3 model for which the deformation radius represents the rigidity of the free surface. Small values of the deformation radius lead to a pliable free surface allowing a significant degree of horizontal divergence. The resulting flow then has an exponentially short interaction range, with a horizontal attenuation on the order of the deformation radius (Polvani et al. 1989), and with approximately non-dispersive Rossby waves. Consequently, for a finite deformation radius, we obtain jets whose width is on the order of the deformation radius with a fixed meandering shape (Scott et al. 2022). In contrast, for the variable stratification surface quasigeostrophic model, rather than just specifying a constant (i.e., the deformation wavenumber), one instead has to specify the stratification's functional form, ( ). Over decreasing stratification [ ( ) < 0], because of the longer interaction range and the more dispersive waves, we obtain jets similar to the two-dimensional barotropic model. Conversely, over increasing stratification [ ( ) > 0], the shorter interaction range along with the weakly dispersive waves lead to sinuous jets whose shape evolves in time through the propagation of weakly dispersive along jet waves. Moreover, because of these along jet waves, a smaller fraction of the total energy is contained in the zonal mode over increasing stratification (with a shorter interaction range) than over decreasing stratification (with a longer interaction range).
The remainder of this article is organized as follows. Section 2 introduces the variable stratification surface quasigeostrophic model and shows how the stratification's vertical structure controls both the interaction range of point vortices as well as the dispersion of surfacetrapped Rossby waves. Then, in section 3, we introduce two wavenumbers, and , whose ratio, / , forms the key non-dimensional parameter of this study; here, is a wavenumber depending on the energy injection rate whereas is a wavenumber depending on surface damping rate. This non-dimensional number is a generalization of the non-dimensional number used in previous studies (Danilov & Gurarie 2002;Sukoriansky et al. 2007;Scott & Dritschel 2012). By considering an idealized buoyancy staircase, we also investigate how the Rhines wavenumber relates to the jet spacing under decreasing, increasing, and uniform stratification. Section 4 then presents numerical experiments detailing the emergence of the staircase limit as / is increased for various stratification profiles. In addition, we also present experiments where we fix the external parameters and vary the vertical stratification alone. Finally, we conclude in section 5.

Equations of motion
Consider an infinitely deep fluid with zero interior potential vorticity. The geostrophic streamfunction, , then satisfies in the fluid interior, ∈ (−∞, 0). The horizontal Laplacian is denoted by ∇ 2 = 2 + 2 and the non-dimensional stratification is given by where ( ) is the buoyancy frequency and is the constant local value of the Coriolis parameter. Time-evolution is determined by the material conservation of surface potential vorticity (Bretherton 1966), at the upper boundary, = 0, where 0 = (0). Explicitly, the time-evolution equation is represents the advection of by the geostrophic velocity, =ˆ× ∇ . The frequency, Λ, is given by where ( ) is a background zonal geostrophic flow. Without loss of generality, we have assumed that (0) = 0 in the time-evolution equation (2.4) to eliminate a constant advective term. The dissipation, , consists of linear damping and small-scale dissipation, where is the damping rate. The forcing, , and the small-scale dissipation, ssd, are described in section 4. The surface buoyancy anomaly, | =0 , is related to the surface potential vorticity, , through Therefore, the time-evolution equation (2.4) equivalently states that surface buoyancy anomalies are materially conserved in the absence of forcing and dissipation. In addition, the frequency, Λ, corresponds to a meridional buoyancy gradient, where ( , ) is the buoyancy field that is in geostrophic balance with background zonal velocity, ( ).
If we further assume a doubly periodic domain in the horizontal, then we can expand the streamfunction as and lower boundary condition (2.12) The upper boundary condition (2.11) is a normalization for the vertical structure, Ψ ( ), chosen so that (2.13) Focus on Fluids articles must not exceed this page length The corresponding Fourier expansion of the surface potential vorticity is given by (2.14) whereˆ= − ( )ˆ, (2.15) and the function ( ) is given by (2.16) The function ( ) relatesˆtoˆin the Fourier space inversion relation (2.15) and so we call ( ) the inversion function.

The inversion function and spatial locality
The inversion function ( ), which is determined by the stratification's vertical structure, controls the spatial locality of the resulting turbulence. We illustrate this point with the 6 following piecewise stratification profile, (2.18) where Δ = 0 − pyc . At small horizontal scales, where , and then ( ) ≈ / 0 , as in the uniformly stratified model of Held et al. (1995). Likewise, in the large-scale limit, where pyc , and (2.20) then ( ) ≈ / pyc . However, for wavenumbers between pyc , the inversion function takes an approximate power law form where 0 > 0 and 0. The power depends on the ratio pyc / 0 between the deep and surface stratification. If the stratification decreases towards the surface [ ( ) 0, or pyc / 0 > 1] then > 1, with pyc / 0 → ∞ sending → 2. In contrast, if the stratification increases towards the surface [ ( ) 0, or pyc / 0 < 1] then < 1, with pyc / 0 → 0 sending → 0. Thus, for wavenumbers pyc , the inversion relation (2.15) has the approximate formˆ= −ˆ, (2.22) whereˆ=ˆ/ 0 , which is the inversion relation for -turbulence (Pierrehumbert et al. 1994;Smith et al. 2002;Sukhatme & Smith 2009). Figure 1 provides two examples, one with decreasing stratification (with ≈ 1.50) and another with increasing stratification (with ≈ 0.49). To see how the parameter modifies the resulting dynamics, consider a point vortex at the origin, given by = (| |), where | | is the horizontal distance from the vortex centre, and (| |) is the Dirac delta. If = 2, then the streamfunction induced by the point vortex is logarithmic, . Smaller leads to vortices with velocities decaying more quickly with the horizontal distance | |, and hence a shorter interaction range. Thus, the vertical stratification modifies the relationship between a surface buoyancy anomaly and its induced velocity field: a surface buoyancy anomaly over decreasing stratification [ ( ) 0] generates a longer range velocity field than an identical buoyancy anomaly over increasing stratification [ ( ) 0].

Wave dispersion in variable stratification
The background gradient term, Λ, in the time-evolution equation (2.4) allows for the propagation of surface-trapped Rossby waves. Substituting a wave solution of the form , where the vertical structure Ψ ( ) satisfies the boundary value problem (2.10)-(2.12), into the time-evolution equation (2.4) yields the angular frequency Given the relationship (2.8) between the meridional surface buoyancy gradient d /d | =0 and the frequency Λ, a poleward decreasing buoyancy gradient ( d /d < 0) implies westward propagating ( < 0) Rossby waves. The dispersion relation (2.23) shows that Rossby wave dispersion is coupled to the flow's interaction range and hence the stratification's vertical structure. If we approximate the inversion function as a power law (2.21) between pyc , then the zonal phase speed, = / , becomes ∼ 1/ . Therefore, at these horizontal scales, Rossby waves are more dispersive over decreasing stratification (with > 1) than over increasing stratification (with < 1). In the limit that 0 pyc in which → 0, then ≈ constant, and so Rossby waves become non-dispersive.

From edge waves to surface-trapped jets
The emergence of jets in barotropic -plane turbulence is due to two properties of the potential vorticity ( , which consists of uniform regions of potential vorticity punctuated by sharp potential vorticity gradients. The second property is that, through potential vorticity inversion, strong (positive) latitudinal gradients in potential vorticity correspond to eastward jets. Therefore, inverting a potential vorticity staircase produces a flow with eastward zonal jets centred at the sharp frontal zones, with weaker westward flows in between (Scott & Dritschel 2019).
However, the limit of a potential vorticity staircase is only achieved for sufficiently large values of the non-dimensional number / Rh (Scott & Dritschel 2012), which is a ratio of the forcing intensity wavenumber, , to the Rhines wavenumber, Rh . The forcing intensity wavenumber is given by (Maltrud & Vallis 1991) where K is the kinetic energy injection rate in the barotropic model, and is obtained by setting the turbulent strain rate equal to the Rossby wave frequency (Vallis & Maltrud 1993). The Rhines wavenumber is given by (Rhines 1975) where rms is the rms velocity. Scott & Dritschel (2012) found that the ratio / Rh controls the structure of zonal jets in barotropic -plane turbulence; as / Rh is increased, the zonal jet strength increases and the potential vorticity gradient at the jet core becomes larger, with the staircase limit approached as / Rh ∼ (10). Jet formation in surface quasigeostrophic turbulence proceeds similarly, with the surface buoyancy (which is proportional to ) taking the role of the potential vorticity and the frequency, Λ, taking the role of the potential vorticity gradient, . In this section, we first derive a non-dimensional number analogous to / Rh for surface quasigeostrophy. Then we consider how vertical stratification (and the non-locality parameter ) modifies jet structure in the buoyancy staircase limit, as well as how it modifies the relationship between the Rhines wavenumber and the jet spacing.
Before proceeding, we comment on two differences between two-dimensional barotropic turbulence and its surface quasigeostrophic counterpart. First, in the absence of forcing and dissipation, the kinetic energy, is a conserved constant in two-dimensional barotropic turbulence (the overline denotes an area average). With a constant kinetic energy injection rate, K , and a linear damping rate, , the equilibrium kinetic energy is K = K /2 . By definition, the rms velocity is given by Combining this expression with the definition of the kinetic energy (3.3) and substituting into the definition of the Rhines wavenumber (3.2) yields a Rhines wavenumber expressed in terms of external parameters alone, In contrast, in surface quasigeostrophy, the total energy, is a conserved constant in the absence of forcing and dissipation and there is no general relationship between the rms velocity, rms , and the equilibrium total energy, E = /2 , where is the total energy injection rate in the surface quasigeostrophic model. Therefore, we are not generally able to express the Rhines wavenumber in terms of the external parameters , Λ, and . Second, because E and K have different dimensions, the kinetic energy injection in the barotropic model, K , has different dimensions than the total energy injection rate in the surface quasigeostrophic model, . In particular, has dimensions of 2 / 3 .

The forcing intensity wavenumber
To obtain the forcing intensity wavenumber, , we compare the Rossby wave frequency (2.23) to the turbulent strain rate, ( ). If the inversion function is not approximately constant (i.e., ≠ 0) then the strain rate is (Yassin & Griffies 2022) In particular, if ( ) = 0 , then ( ) ∼ 1/3 0 1/3 (4− )/3 . Setting the absolute value of the Rossby wave frequency for waves with = equal to the turbulent strain rate (3.6) yields the condition (3.9) 9 However, unlike in two-dimensional barotropic turbulence where rms = √ 2K = √︁ K / , we do not have a general relationship between rms and the external parameters and in surface quasigeostrophic turbulence. To obtain a second wavenumber that depends on the damping rate, , we follow Smith et al. (2002). From dimensional considerations, the energy spectrum at small wavenumbers is (3.10) Then, defining as the wavenumber at which the inverse cascade halts, we obtain where the second equality follows because the integral is dominated by its peak at low wavenumbers. Solving for and neglecting any non-dimensional coefficients, we obtain . (3.12) Note that the damping rate wavenumber, , has the same dependence on Λ, , and as the Rhines wavenumber, Rh , only if = 2.
3.3. Surface potential vorticity inversion A perfect surface potential vorticity staircase consists of mixed zones of halfwidth , where d /d = −Λ, separated by jump discontinuities at which d /d = ∞. We find it more conveniant to work with the relative surface potential vorticity, , rather than the total surface potential vorticity, + Λ . In this case, if the total surface potential vorticity, + Λ , is a perfect staircase with step width 2 , then the relative surface potential vorticity, , is a 2 -periodic sawtooth wave.
Our first question is whether such a staircase is possible for general ( ). To answer this question, we consider the velocity field induced by a jump discontinuity in . For a jump discontinuity in an infinite domain, the zonal velocity is given by (3.14) If ( ) = 0 , then this expression is proportional to | | −1 if ≠ 1 and logarithmic otherwise, and so the zonal velocity diverges at = 0 if 1. Consequently, we expect that a perfect staircase should not be possible over constant or increasing stratification due to the divergence of the zonal velocity at a jump discontinuity.
We therefore consider the more general case of a sloping staircase, where there is a finite frontal zone of width 2 between the mixed zones. In this case, is a 2( + )-periodic sloping sawtooth wave (see figure 2), and is given by the periodic extension of for − ( + ) < < − . The meridional gradient d /d is then a piecewise constant 2( + )-periodic function (3.16) Therefore the gradient in the frontal zones exceeds the gradient in the mixed zones by a factor of / , which approaches infinity as / → ∞ in the sawtooth wave limit. The zonal velocity, = − , is obtained by using the inversion relation (2.15) to solve for the streamfunction. Alternatively, taking the meridional derivative of surface potential vorticity (2.3) gives which shows that the induced zonal velocity is obtained by smoothing d /d by the function ( ). An immediate consequence is that the east-west asymmetry in the zonal velocity is fundamentally due to the east-west asymmetry in the gradient d /d . Figure 2 shows an example of sloping sawtooth profile along with the induced zonal velocities. For a power law inversion function, ( ) = 0 , the parameter modifies the zonal velocity in two ways. First, in more local flows (with smaller ), the zonal velocity decays more rapidly away from the jet centre, as expected. Second, the degree of smoothing increases with , and so more local regimes (with smaller ) are more east-west asymmetric, with the ratio | min | / max taking smaller values for smaller . Figure 3(b) shows | min | / max as a function of / for ∈ {1/2, 1, 3/2, 2}. For = 2, we obtain | min | / max → 1/2 in the limit / → 0 so that eastward jets are only twice as strong as westward flows in the perfect staircase limit (Danilov & Gurarie 2004;Dritschel & McIntyre 2008). At = 3/2, we find | min | / max ≈ 0.29 in the / → 0 limit so that eastward jets are now more than three time as strong as westward flows. Once 1, then the maximum jet velocity diverges as → 0 [figure 3(a)] and so | min | / max → 0 as / → 0.
If ( ) is not a power law, then the results are similar so long as ( ) can be approximated by a power law at small wavenumbers. Figure 2 shows the induced velocity for the inversion functions computed from idealized stratifications profiles (shown in figure 1). Because these inversion functions can be approximated by power laws ( ) ≈ 0.49 and ( ) ≈ 1.50 at small wavenumbers, the induced velocity fields nearly coincide with the velocity fields computed from power law inversion functions with = 0.5 and = 1.5.
Finally, we examine how the Rhines wavenumber, Rh , relates to jet spacing. Let be the half-separation between the jets, i.e., the half distance between consecutive zonal velocity maxima. For two-dimensional barotropic turbulence (i.e., the = 2 case), we have = 45 1/4 / Rh ≈ 2.59/ Rh in the staircase limit (i.e, for / → 0, Dritschel & McIntyre 2008; Scott & Dritschel 2012). This result is found by solving for the zonal velocity induced by a staircase with halfwidth = , taking the rms of the zonal velocity, and then substituting into the definition of the generalized Rhines wavenumber (3.9). As figure 3(d) shows, because the velocity field induced by a perfect staircase depends on the inversion function, ( ), the relationship between and Rh also depends on the inversion function. For ( ) = 3/2 , an analogous calculation gives ≈ 2.35/ Rh in the staircase limit. For = 1, even though the maximum velocity diverges at / → 0, the rms velocity asymptotes to a constant value, and so we obtain a half jet-separation of ≈ 1.73/ Rh (figure 3). Finally in the = 1/2 case, although the rms speed has not converged by / = 10 −6 , the product Rh is approaching values close to zero. with = 23.6 and 0 = 0.65 Nyq where Nyq = is the Nyquist wavenumber. The forcing is isotropic, centred at wavenumber = 80, and normalized so that the energy injection rate is = 1 (see appendix B in Smith et al. 2002). However, the effective energy injection rate, eff , is smaller than due to dissipation. To determine eff from numerical simulations, we use eff = 2 E where E is the equilibrated total energy diagnosed from the model. In what follows, we report values of / using eff instead of . The model is integrated forward in time until at least = 5/ to allow the fluid to reach equilibrium. All model runs use 1024 2 horizontal grid points.

For what values of / do jets form?
For our first set of simulations, we vary / over the values shown in figure 4. We do so by fixing = 8 and varying . For a given value of , we choose Λ and so as to maintain = 8 (the energy injection rate, , is fixed at unity for all model runs). Given and , we rearrange the definition of (3.12) to solve for = Λ 2 , and finally use the definition = Λ 2 to solve for Λ.

Power law inversion functions
We first describe the results from three series of simulations with power law inversion functions, ( ) = , with ∈ { 1 /2, 1, 3 /2}. Summary diagnostics from these simulations are shown in figure 4. In panel (a), we observe that the ratio of energy in the zonal mode to total energy, E zonal /E, increases with / , and that the majority of the total energy is in the zonal mode for sufficiently large / . For a fixed / , more of the total energy is zonal in more non-local flows (with larger ) than in more local flows (with smaller ); for = 3 /2, we have E zonal /E ≈ 1 by / ≈ 6 as compared to / ≈ 12 for = 1. Moreover, for = 1 /2, we find that E zonal /E asymptotes to approximately 0.9 once / ≈ 18 with little subsequent change for larger values of / . In panel (b), we observe a striking contrast in the ratio | |/| | between different values of (the overline denotes a domain average). For = 3 /2, the domain averaged zonal speed, | |, is approximately eight times larger than the domain averaged meridional speed, | |, for large / . In contrast, for = 1 /2, | | only exceeds | | by a multiple of two for large / .
Next, we examine the jet structure for different as a function of / . Figure 5 shows -snapshots from model runs with ( ) = . For each value of , two model runs are shown: one where jets have just become visible in the -snapshot and another with the largest value of / , which we expect to be closest to the staircase limit. The jets are visible in these snapshots as the regions with strong gradients. Because these are -snapshots rather than ( + Λ )-snapshots, the ( + Λ )-staircase is instead a -sawtooth, and the mixed zones between the jets are approximately linear in . We confirm this to be the case in figure 6, where the zonal averages of the total surface potential vorticity, +Λ , and the zonal velocity are shown. For the = 3 /2 and = 1 cases, we observe an approximate staircase structure with nearly uniform mixed zones separated by frontal zones, and with jets centred at sharp gradients. As expected from the idealized staircases of section 3, close to the staircase limit, the = 1 jets are narrower than the = 3 /2 jets, and the ratio of maximum westward speed to maximum eastward speed, | min |/| max |, is smaller at = 1 than at = 3 /2.
In contrast to the = 3 /2 and the = 1 series, the = 1 /2 series approaches the staircase limit slowly with / . The = 1 /2 staircase remains smooth even at / = 42 [figure 6(c)]. The ratio of frontal zone width to mixed zone width, / , is between 0.5 and 0.65 for = 1 /2 jets. In contrast, this ratio is between 0.15 and 0.2 for the = 3 /2 and = 1 jets. In part, the broadness of the = 1 /2 frontal zones is a consequence of zonal averaging in Figure 5: Snapshots of the relative surface potential vorticity, , for simulations with power law inversion functions, ( ) = . In each snapshot, the field is normalized by its maximum value in the snapshot. Only one quarter of the domain is shown (i.e., 512 2 grid points).
the presence of large amplitude undulations. However, it is evident from the -snapshots of figure 5 that the = 1 /2 frontal zones are indeed broader than the = 3 /2 and = 1 frontal zones [e.g., compare panels (a) and (d) with (f) in figure 5], even without zonal averaging. We now examine how the generalized Rhines wavenumber, Rh , relates to the jet spacing. From figure 3(d), a ratio of / ≈ 0.2 leads to a Rh ≈ 2.2 for = 3 /2 and Rh ≈ 2.0 for = 1. But as figure 4(d) shows, we find values closer to Rh ≈ 3 for both of these cases. In contrast, for the = 1 /2 jets, figure 3(d) predicts 1.98 Rh 2.5 for the observed range of 0.5 / 0.65, but we find Rh ≈ 1.5 for / 18, which is smaller than predicted.
Returning to figure 5, we observe that there are undulations along the jets, with smaller values of corresponding to larger amplitude undulations. These undulations propagate as waves and are less dispersive for smaller , propagating eastward for = 3 2 , westward for = 1 /2, and are nearly stationary for = 1. Moreover, the waves in the = 1 /2 case maintain their shape as they propagate for a significant fraction of the domain, although they eventually disperse or merge with other along jet waves. That we obtain larger amplitude along jet undulations for smaller is a consequence of the more local inversion operator (2.15) at smaller . A jet in a highly local flow (with small ) is "a coherent structure that hangs together strongly while being easy to push sideways" (McIntyre 2008, in the context of equivalent barotropic jets). However, although both an equivalent barotropic jet and an = 1 /2 jet exhibit large meridional undulations, the undulations in the equivalent barotropic case are frozen in place (because of a vanishing group velocity at large scales, McIntyre 2008) and so the equivalent barotropic jet behaves like a meandering river with a fixed shape. In contrast, the = 1 /2 jet behaves like a flexible string whose shape evolves in time with the propagation of weakly dispersive waves. Another difference between the two cases is that an equivalent barotropic jet has a width given by the deformation radius. In contrast, there is no analogous characteristic scale for = 1 /2 jets and, in principle, the jets should become infinitely thin as / → ∞.
Energy spectra for the three power law simulations are shown in figure 7. The energy spectrum obtained from dimensional analysis (3.10) gives a − −3 wavenumber dependence, which leads to the familiar −5 spectrum for beta-plane barotropic turbulence ( = 2). Although early investigations (Chekhlov et al. 1996;Huang et al. 2000;Danilov & Gryanik 2004) found a −5 spectrum in barotropic -plane turbulence, Scott & Dritschel (2012) instead found a shallower −4 spectrum in the staircase limit (suggested earlier by Danilov & Gryanik 2004;Danilov & Gurarie 2004), which they explained as a consequence of the sharp discontinuities of the staircase. Generalizing their argument to the present case, a one dimensional ( ) series with discontinuities implies a Fourier series with coefficients  If ( ) ∼ , then we obtain a spectrum ( ) ∼ − −2 , which yields the −4 spectrum observed in Scott & Dritschel (2012), where = 2. For = 3 /2, = 1, and = 1 /2, the predicted spectrum is proportional to −3.5 , −3 , and −2.5 , respectively. The diagnosed spectra shown in figure 7 are consistent with these shallow spectra, instead of energy spectrum (3.10) obtained from dimensional considerations. 60. As seen in figure 4, the ( ) 0 case is similar to the = 3 /2 case, with the various diagnostics close to the = 3 /2 counterpart. In contrast, there are significant differences between the ( ) 0 simulations and the = 1 /2 simulations. In the ( ) 0 series, the ratio of energy in the zonal mode to total energy continues to increase as / is increased, whereas it asymptotes to a constant in the = 1 /2 series. Moreover, the ratio of domain average zonal speed to domain averaged meridional speed, | |/| |, is generally larger in the 0 series than in the = 1 /2 series. Finally, for the largest values of / , the product Rh reaches smaller values in the 0 simulations than in the = 1 /2 simulations. These differences can be explained by the snapshots of figure 9 as well as the zonal averages of figure 6. As expected from the model diagnostics, both the snapshots and the zonal average from the 0 simulation are qualitatively similar to the = 3 /2 simulation. In contrast, the 0 snapshot is evidently closer to the staircase limit than the = 1 /2 snapshot: the mixed zones are more homogeneous and the frontal zones are sharper. The zonal average of the 0 simulation in figure 6 also shows how the 0 simulation is closer to the staircase limit than the = 1 /2 simulation, although, again, zonal averaging in the presence of large amplitude undulations is artificially smoothing the jets. Therefore, the differences in the diagnostics between the 0 series and the = 1 /2 series stem from the more rapid approach (i.e., at smaller / ) of the 0 series to the staircase limit.

Simulations with fixed parameters
The dependence of the non-dimensional number / on the external parameters , Λ, and depends on the functional form of ( ). For example, if ( ) ∼ , then (4.5) Because the forcing intensity wavenumber, , is obtained by solving the implicit equation for (3.7), an analogous expression for / is not possible for general ( ). However, at sufficiently large , the inversion function asymptotes to ( ) ≈ / 0 and so, using -turbulence expression for (3.8) with = 1, we obtain for large , where is the approximate power law dependence of ( ) near . Therefore, simulations with identical / but distinct inversion functions cannot be directly compared because they have different values of Λ and . Here, we investigate how the stratification modifies jet structure as all other parameters are held fixed. We therefore run two additional series of simulations with the stratification profiles and inversion functions shown in figure 1. The stratification profiles were chosen so that they both have identical stratification at the upper boundary. One case corresponds to an increasing stratification profile, 0, with an approximate power law dependence of ( ) ∼ 0.49 at small wavenumbers. The second case consists of a decreasing stratification profile, 0, with a ( ) ∼ 1.50 at small wavenumbers. Aside from the different stratification profiles, these two series of simulations are run under the same conditions as the constant stratification ( = 1) simulations of section 4.2, with identical values of Λ, , and .
Summary diagnostics are shown in figure 10. We see that, at a fixed value of Λ and , more of the total energy is in the zonal mode in the ( ) 0 simulation than in the constant stratification simulation, which in turn is larger than the ( ) 0 simulation (and similarly for the ratio of area averaged zonal to meridional speeds, | |/| |). Therefore, increased non-locality (larger ) promotes anisotropy in the velocity field and leads to larger zonal velocities relative to meridional velocities. Indeed, figure 11 shows snapshots from these simulations; the more local, 0, simulations have larger meridional undulations along the jets. Moreover, compared to the / = 15 constant stratification simulation in

Conclusion
We have examined the emergence of staircase-like buoyancy structures in surface quasigeostrophic turbulence with a mean background buoyancy gradient. We found that the stratification's vertical structure controls the locality of the inversion operator and the dispersion of surface-trapped Rossby waves. As we go from decreasing stratification profiles [ ( ) 0] to increasing stratification profiles [ ( ) 0], the inversion operator becomes more local and Rossby wave less dispersive. In all cases, we find that the nondimensional ratio, / , controls the extent of inhomogeneous buoyancy mixing. Larger / correspond to sharper buoyancy gradients at jet centres with larger peak jet velocities that are separated by more homogeneous mixed-zones. Moreover, we found that the staircase limit is reached at smaller / in more non-local flows; the staircase limit is reached by / ≈ 15 for our 0 simulations compared to / ≈ 25 for our 0 simulations.

20
In addition, once the staircase limit is reached, the dynamics of the jets depends on the locality of the inversion operator and, hence, on the stratification's vertical structure. In flows with a more non-local inversion operator [or decreasing stratification, ( ) 0], we obtain straight jets that are perturbed by dispersive, eastward propagating, along jet waves. In contrast, for more local flows [or over increasing stratification, ( ) 0], we obtain jets with latitudinal meanders on the order of the jet spacing. The shape of these jets evolves in time as these meanders propagate westwards as weakly dispersive waves.
The inversion operator's locality is also reflected in two more aspects of the dynamics. First, the domain-averaged zonal speed exceeds the domain-averaged meridional speed by approximately a factor of eight in our most non-local simulations, whereas this ratio is merely two in our most local simulations. This observation is consistent with the fact that jets are narrower and exhibit larger latitudinal meanders in more local flows. Second, for a given Rhines wavenumber, jets in more local flows are closer together. Indeed, we found Rh ≈ 3 − 4 in our most non-local simulations, where is the jet half spacing, as compared to Rh ≈ 0.5 − 1.5 in our most local simulations. Several open questions remain. First, we have not examined the dynamics of the along jet waves. As we observed, these waves propagate eastwards in our most non-local simulations [with ( ) 0] but westwards for our most local simulations [with ( ) 0]. These waves are not described by the dispersion relation (2.23); rather, the relevant model is that of freely propagating edge waves along a buoyancy discontinuity (McIntyre 2008). However, the difficulty here is that a jump discontinuity in the buoyancy field results in infinite velocities over constant or increasing stratification. In addition, the relationship of the along jet waves in the staircase limit to the non-linear zonons found by Sukoriansky et al. (2008) remains unclear.
The divergence of the velocity at a buoyancy discontinuities raises a second question. Is there a limit to how close the staircase limit can be approached? In barotropic dynamics, the velocity remains finite at a jump continuity in the vorticity, and, in this case, Scott & Dritschel (2012) report that a vorticity staircase case can be approached arbitrarily. Whether this result continues to hold for arbitrarily sharp buoyancy gradients and arbitrarily large zonal velocities is not clear. Because the rms velocity seems to converge for arbitrarily sharp staircases, even for the most local inversion relations we considered, there may not be any energetic reason precluding arbitrarily sharp buoyancy gradients.
Finally, there remains the question of how relevant these results are for the upper ocean, which, in addition to surface buoyancy gradients, has interior potential vorticity gradients as well. In particular, our neglect of the -effect limits the direct relevance of this model to the upper ocean. Whether surface buoyancy staircases can emerge under more realistic oceanic conditions requires further investigation.