1. Introduction
Potential vorticity (PV) mixing and the formation of fronts, concentrated bands of strong PV gradients at the edges of mixed regions, are important processes in large-scale geophysical flows, with examples ranging from the Earth's ocean and winter stratosphere to Jupiter's banded structure and Great Red Spot (see, e.g. Galperin & Read Reference Galperin and Read2019). Because of the PV invertibility principle (Dritschel & McIntyre Reference Dritschel and McIntyre2008), along these fronts run strong, fast-moving jets containing large amounts of kinetic energy. Fronts and jets are important, for example, in mesoscale oceanic flows, where fronts act as barriers to transport of heat, momentum and chemical constituents across the front, while jets may transport quantities long distances in the along-front direction.
In this paper the long-time dynamics of a canonical geophysical model, shallow water quasi-geostrophic flow, are studied for a range of initial conditions and model parameters leading to flow regimes dominated by PV fronts and their collocated jets. Shallow water quasi-geostrophic flow describes a rapidly rotating shallow fluid layer with a deformable free surface. This model simplifies the dynamics while retaining key processes active in real geophysical flows, and has been used to model dynamical features occurring at oceanic fronts (see, e.g. Pratt & Stern Reference Pratt and Stern1986). The system is governed by the Charney–Hasegawa–Mima equation (see, e.g. Pedlosky Reference Pedlosky1987)
\begin{equation} \frac{\partial q}{\partial t} + J(\psi,q) = 0, \end{equation}
which describes material advection of PV 
$q = (\varDelta - L_D^{-2})\psi$, where 
$\varDelta$ is the two-dimensional (2-D) Laplacian, 
$\psi$ is the streamfunction, 
$\Delta \psi$ is the vorticity and 
$J(\cdot,\cdot )$ is the 2-D Jacobian. The deformation wavenumber 
$k_D = L_D^{-1}$ is the inverse of the Rossby deformation length 
$L_D = \sqrt {gH}/f$, where 
$g$ is the gravitational acceleration, 
$H$ is the mean layer depth and 
$f$ is the Coriolis parameter; 
$L_D$ measures the relative tendencies of gravity to relax the free surface and of background planetary rotation to maintain free surface height anomalies. Equation (1.1) also governs quasi-2-D fluctuations of the electrostatic potential 
$\psi$ for a plasma in a uniform strong magnetic field, where then 
$L_D$ is the ion Larmor radius (Hasegawa & Mima Reference Hasegawa and Mima1978).
The inviscid quadratic invariants of (1.1) are the total energy
\begin{equation} E \equiv \tfrac{1}{2}\langle -\psi q \rangle = \tfrac{1}{2} [ \langle |\nabla \psi|^2 \rangle + k_D^2 \langle \psi^2 \rangle ] = K + P, \end{equation}
where 
$K \equiv \langle |\nabla \psi |^2 \rangle /2$ is the kinetic energy and 
$P \equiv k_D^2 \langle \psi ^2 \rangle /2$ is the potential energy, and the potential enstrophy
\begin{equation} Q \equiv \tfrac{1}{2}\langle \Delta \psi q \rangle = \tfrac{1}{2} [ \langle |\Delta \psi|^2 \rangle + k_D^2 \langle |\nabla\psi|^2 \rangle ], \end{equation}
which is the sum of the barotropic enstrophy 
$Z \equiv \langle |\Delta \psi |^2 \rangle /2$ and the rescaled kinetic energy 
$k_D^2 K \equiv k_D^2\langle |\nabla \psi |^2 \rangle /2$. Here the angle brackets 
$\langle \cdot \rangle$ denote an integral over the domain. The inviscid dynamics also conserve an infinite hierarchy of PV norms (Casimirs), including the 
$L^2$ norm 
$C \equiv \langle q^2 \rangle /2$.
The dynamics of (1.1) have been extensively studied in both the limits 
$L \ll L_D$ and 
$L \gg L_D$, where 
$L$ is the characteristic length scale of the flow. In the first limit, 
$k_D L \ll 1$, deformations of the free surface are negligible and the turbulence is effectively governed by the 2-D Euler equations, or material advection of vorticity 
$\Delta \psi$. In this case the quadratic invariants are kinetic energy 
$K$, which cascades to large scales, and enstrophy 
$Z$, which cascades to small scales.
Conversely, when 
$L \gg L_D$, i.e. the characteristic length scale of the flow is much larger than 
$L_D$, the streamfunction dominates the PV and the turbulent evolution slows down (Larichev & McWilliams Reference Larichev and McWilliams1991; Iwayama, Shepherd & Watanabe Reference Iwayama, Shepherd and Watanabe2002). The limit 
$k_D L \rightarrow \infty$ yields the asymptotic model (AM),
\begin{equation} \frac{\partial \psi}{\partial \tau} + J(\Delta \psi, \psi) = 0, \end{equation}
which describes material advection of the streamfunction 
$\psi$ by the vorticity 
$\Delta \psi$ on a rescaled slow time 
$\tau = t/k_D^2$ (Larichev & McWilliams Reference Larichev and McWilliams1991). The AM is a singular limit because 
$\tau \rightarrow 0$ as 
$L_D \rightarrow 0$, and the evolution slows to a halt; see, e.g. Cushman-Roisin (Reference Cushman-Roisin1987). This model assumes that small scales never develop. The quadratic invariants of the AM are the kinetic energy 
$K$ and the rescaled potential energy 
$P/k_D^2$, which are expected to undergo forward and inverse cascades, respectively. As 
$L/L_D \rightarrow \infty$, the slow time 
$\tau \rightarrow 0$, consistent with the frozen vortical quasi-crystals observed in simulations of smooth PV fields with length scale 
$L \gg L_D$ (Larichev & McWilliams Reference Larichev and McWilliams1991; Boffetta, De Lillo & Musacchio Reference Boffetta, De Lillo and Musacchio2002; Iwayama et al. Reference Iwayama, Shepherd and Watanabe2002), and with the analytical results of Tran & Dritschel (Reference Tran and Dritschel2006), which indicate that smooth vortices much larger than 
$L_D$ are inactive, cannot merge and are thus prevented from transferring potential energy to larger scales. The AM regime may not be self-consistent, however, since forward cascades of kinetic energy may generate flow features with scales 
${O}(L_D)$ (already anticipated by Larichev & McWilliams Reference Larichev and McWilliams1991).
Large-scale geophysical flows generically form steep PV gradients, or fronts, separating regions in which PV is well mixed (McIntyre Reference McIntyre1982; Dritschel & McIntyre Reference Dritschel and McIntyre2008; Dunkerton & Scott Reference Dunkerton and Scott2008; Dritschel & Scott Reference Dritschel and Scott2011; Scott & Dritschel Reference Scott and Dritschel2018), so theories for perfectly smooth PV distributions are of limited relevance to the dynamics observed in real oceanic and atmospheric flows. For equivalent-barotropic turbulence in particular, the above numerical and analytical results pertaining to smooth PV distributions cannot be generalized to flows containing steep gradients, even when 
$L \gg L_D$ and the flow is dominated by potential energy. In the thin jet limit, in which the jets are much thinner than the length scale of the mixed regions, the dynamics of the waves and nonlinear disturbances that propagate on the fronts are described to first order by a modified KdV equation, which links frequency to the length scale of the mixed regions (Nycander, Dritschel & Sutyrin Reference Nycander, Dritschel and Sutyrin1993; Burgess & Dritschel Reference Burgess and Dritschel2019; Burgess Reference Burgess2020).
An outstanding question, addressed here, is whether sharp fronts inevitably form even for initial conditions with a characteristic length scale 
$L$ much larger than the Rossby deformation length 
$L_D$, which sets the width of the fronts and their collocated jets. Answering this question is crucially important to understanding the long-time dynamics of geophysical flows. Another question is how flow properties such as the number of mixed levels, typical area of mixed regions and jet spacing depend on the initial conditions and model parameters.
To begin answering these questions, this paper studies quasi-geostrophic shallow water flows with initial characteristic length scales 
$L_0$ larger than or equal to the deformation length 
$L_D$, 
$L_0/L_D \geq 1$. (Note, in geophysical fluid dynamics, the ratio 
$(L_D/L_0)^2$ is known as the Burger number; see, e.g. Pedlosky Reference Pedlosky1987.) Two sets of simulations are analysed, one set with 
$k_D = L_D^{-1} = 25$ and the other with 
$k_D = 40$. In each set, the peak wavenumber 
$k_0 = L_0^{-1}$ of the initial energy spectrum, which provides an initial length scale 
$L_0 = 1/k_0$ for the flow, is varied from 
$k_0 = 5$ to 
$k_0 = k_D$.
This parameter regime in which the deformation length 
$L_D$ is small compared with the characteristic length scale of the flow is most relevant to the ocean. Like oceanic jets, the jets that emerge are very ‘wiggly’, exhibiting large meanders, and enclose regions of well-mixed PV, forming rings reminiscent of those seen in the Gulf Stream. We stress that oceanic flows are much more complicated and the present idealised model omits important effects active in the ocean, such as baroclinicity. Nonetheless, the model can be used to derive important insights about meandering jets in the small 
$L_D$ limit, and has been used previously to investigate oceanic jets (see, e.g. Pratt & Stern Reference Pratt and Stern1986).
The paper is structured as follows. Section 2 describes the numerical simulations, including the choice of initial conditions and governing parameters. We then examine the time evolution of energy and characteristic scales for these simulations in § 3, and the structure of the PV and kinetic energy fields in § 4. The PV field in all cases is found to form a distinct staircase structure at sufficiently late times, with kinetic energy concentrated on sharp PV gradients. A simple theory accounting for the observed variation of this staircase structure with 
$L_0/L_D$ is developed in § 5. This theory uses the initial energy to predict the number of steps that develop and the area between them. Our conclusions are offered in § 6.
2. Numerical simulations
Equation (1.1) is solved by contour advection using the combined Lagrangian advection method (Dritschel & Fontane Reference Dritschel and Fontane2010) on a 
$1024^2$ basic inversion grid with effective resolution 
$(16 \times 1024)^2 = 16\,384^2$ within a domain of side length of 
$2 {\rm \pi}$. Contour surgery removes PV filaments at 
$1/16\,384$ the domain width, but preserves sharp gradients indefinitely, so is ideally suited to investigating dynamics associated with PV fronts. In contrast, pseudospectral simulations with Laplacian or hyperviscosity inevitably degrade PV fronts, so long-time dynamics cannot be investigated as accurately as with contour advection.
The initial spectrum takes the form
\begin{equation} \mathcal{C}(k) = c(k_D^2 + k^2)k^{2p-3}\,{\rm e}^{-(p-1)(k/k_0)^2}, \end{equation}
where 
$\mathcal {C}(k)$ is the power spectrum of 
$C \equiv \langle q^2 \rangle /2$ and 
$p = 3$ for all simulations. Without loss of generality, we take the maximum PV 
$|{q_{max}}| = 4{\rm \pi}$, and this only controls the time scale of the flow evolution (the ‘eddy turn-around time’ is then a unit time in the limit 
$k_0/k_D \rightarrow \infty$). The maximum PV determines the constant 
$c$ in (2.1). Two sets of simulations are conducted, one with deformation wavenumber 
$k_D = 25$ and the other with 
$k_D = 40$. These are sufficiently large to keep domain effects insignificant over the duration of the simulations. In each set the spectral peak wavenumber is varied from 
$k_0 = 5$ to 
$k_0 = k_D$ in increments of 
$\Delta k_0 = 5$.
The simulations start with spatially random smooth PV fields, examples of which are shown in figure 1 for simulations with 
$k_D = 1/25$ and 
$L_0 = 1/5$, 
$L_0 = 1/15$, and 
$L_0 = 1/25$. Potential vorticity concentrations in the initial field merge and evolve into small intense vortices, around which the mixed PV regions subsequently organise.
Figure 1. Initial PV for simulations with 
$L_D=1/25$ and, from (a) to (c), 
$L_0/L_D = 5$, 
$L_0/L_D = 1.67$ and 
$L_0/L_D = 1$.
3. Scale evolution and decay rates
Total energy is very nearly conserved in all simulations, as shown in figure 2. As 
$L_0/L_D$ increases, potential energy makes up an increasingly large fraction of the initial energy. No matter the value of 
$L_0/L_D$ the kinetic energy decreases at sufficiently late times. This occurs even for the simulation with 
$L_0/L_D = 8$ (b, red), though the decline in kinetic energy (red dotted line) is somewhat delayed, commencing at about 
$t = 3000$, whereas simulations with smaller values of 
$L_0/L_D$ show an immediate decrease of kinetic energy. The delay in kinetic energy decay is related to the fact that the dynamics slow down and front formation is delayed as 
$L_0/L_D$ increases. However, the amount of kinetic energy at later times once the decay has commenced is very similar between simulations with the same value of 
$L_D$, and shows less dependence on the value of the initial characteristic length scale 
$L_0$.
Figure 2. Time evolution of total energy 
$E$ (solid line), potential energy 
${P}$ (dash-dot line) and kinetic energy 
${K}$ (dotted line) for 
$L_D = 1/25$ (a) with 
$L_0/L_D=5$ (red), 
$L_0/L_D=2.5$ (blue) and 
$L_0/L_D=1.25$ (black), and 
$L_D = 1/40$ (b) with 
$L_0/L_D=8$ (red), 
$L_0/L_D=2.67$ (blue) and 
$L_0/L_D=1.33$ (black).
This is especially true for the simulations with 
$L_D = 1/40$. Past a certain time, all the kinetic energy curves roughly collapse and follow the same decay law 
${K} \approx t^{-0.4}$. This steeper decay than the 
$t^{-1/3}$ found for the frontal jets alone (see Burgess Reference Burgess2020) reflects the kinetic energy contribution of the vortex cores within the mixed regions of non-zero PV (as opposed to the free cores located in the zero-PV regions). The kinetic energy in these cores decays more quickly than the kinetic energy in the frontal jets. The flow structure is discussed further in § 4.
To study how the characteristic length scales of total, potential and kinetic energy evolve, and in particular to look for evidence of inverse transfer of potential energy, we compute the energy centroid length scales. The centroid length scale of total energy is
\begin{equation} \ell_E = \frac{ \langle (-\varDelta)^{{-}1/2} \psi q \rangle }{ \langle \psi q \rangle } = \frac{ \langle | (-\varDelta)^{1/4} \psi |^2 \rangle + k_D^2 \langle | (-\varDelta)^{{-}1/4} \psi |^2 \rangle }{ \langle | \nabla \psi |^2 \rangle + k_D^2 \langle \psi^2 \rangle }, \end{equation}
where the operator 
$(-\varDelta )^{\theta }$ is defined by 
$(-\varDelta )^{\theta } \hat {\psi }(\boldsymbol {k}) = k^{2\theta } \hat {\psi }(\boldsymbol {k})$ (Tran & Dritschel 2006). Here, 
$k = | \boldsymbol {k} |$ is the wavenumber and 
$\hat {\psi }(\boldsymbol {k})$ is the Fourier transform of 
$\psi$. The characteristic length scales of the kinetic and potential energy are
\begin{equation} \ell_{K} = \frac{\langle | (-\varDelta)^{1/4} \psi |^2 \rangle}{ \langle | \nabla \psi |^2 \rangle}, \quad \ell_{P} = \frac{\langle | (-\varDelta)^{{-}1/4} \psi |^2 \rangle }{\langle \psi^2 \rangle }. \end{equation}
In practice, 
$\ell _E$, 
$\ell _P$ and 
$\ell _K$ are computed by summing over spectra and then dividing by the relevant energy.
At sufficiently late times in all simulations, potential energy is transferred to larger scales, as evidenced by the increase in the centroid length scale 
$\ell _P$ of the potential energy as shown in figure 3 for simulations with 
$L_D = 1/40$. For the simulation with 
$L_0/L_D = 8$ shown in panel (a), the potential energy length scale is approximately constant until 
$t \approx 3000$, after which it begins to increase, demonstrating that inverse transfer of potential energy occurs even in this case, in which the flow starts deep in the AM regime, and theories assuming smooth PV fields would predict impeded inverse potential energy transfer and an effectively infinite time scale for flow evolution (Tran & Dritschel Reference Tran and Dritschel2006). Note that the inverse transfer of potential energy and the decay of kinetic energy begin at roughly the same time (
$t \approx 3000$) for the simulation with 
$L_0/L_D = 8$.
Figure 3. Time evolution of centroid length scales 
$\ell _E$ of the total energy (solid black line), 
$\ell _P$ of the potential energy (dashed magenta line) and 
$\ell _K$ of the kinetic energy (dotted blue line) for simulations with 
$L_D=1/40$ and 
$L_0/L_D=8$ (a), 
$L_0/L_D=4$ (b), 
$L_0/L_D=2.67$ (c), 
$L_0/L_D=2$ (d), 
$L_0/L_D=1.6$ (e) and 
$L_0/L_D=1.33$ (f). The ratio 
$\ell _E/\ell _P$ (black dash-dot line) is also shown.
For the simulations with smaller 
$L_0/L_D$, i.e. 
$L_0/L_D = 2.67$ to 
$L_0/L_D = 1.33$ shown in figure 3(c)–(f), the length scale of the kinetic energy first increases and then begins to flatten out. This indicates the emergence of PV fronts and the concentration of kinetic energy in the jets collocated with fronts, which have a characteristic width 
${O}(L_D)$.
4. Flow structure
The structure of the flow varies significantly with 
$L_0/L_D$ as evident in figure 4, which shows PV fields (panels a–c) and kinetic energy density (panels d–f) for simulations with L D = 1/25. For 
$L_0/L_D = 5$ (figure 4a), the domain is almost completely occupied by mixed regions of non-zero PV, and there are about nine mixed levels (including both positive and negative PV regions). The mixed regions are arranged into what we call ‘multi-level vortices’, which have a ‘wedding cake’ structure, with ascending terraces on which mixed PV takes increasingly high values. Levels with low and intermediate values of 
$|q|$ occupy the most area, and levels with the highest values of 
$|q|$, which are close to the vortex cores where mixing is inhibited, occupy the least area. As 
$L_0/L_D$ decreases to 1.67 and then 1 (b and c), the multi-level vortices are separated by larger and larger regions on which 
$q \approx 0$, and the non-zero mixed regions become island-like rather than forming continuously connected regions spanning the domain. For smaller values of 
$L_0$, there are also fewer non-zero mixed levels, with about four mixed levels for the case 
$L_0/L_D = 1$ (c).
Figure 4. Potential vorticity (a–c) and kinetic energy (d–f) for simulations with 
$L_D=1/25$ and 
$L_0/L_D = 5$ (a,d), 
$L_0/L_D = 1.67$ (b,e) and 
$L_0/L_D = 1$ (c,f), all at 
$t=200 \, 000$. The domain length 
$L = 2 {\rm \pi}$ is 25 deformation wavelengths 
$2{\rm \pi} L_D$.
The kinetic energy fields and jet intensity also vary depending on 
$L_0/L_D$, as can be seen in plots (d)–(f) of figure 4. Note that the same colour range is used for all three plots to facilitate comparison of the jet intensities. Jets are on average weaker and more densely bunched together when 
$L_0/L_D$ is larger (figure 4d) and the domain is filled with more non-zero mixed levels, while smaller 
$L_0/L_D$ yields fewer mixed levels and stronger jets (figure 4f). The total length of the jets also varies depending on 
$L_0/L_D$: for larger 
$L_0/L_D$ (figure 4d), the total jet length is greater, while for smaller 
$L_0/L_D$ (figure 4f), there is less total jet length in the domain but the jets are stronger. This accounts for the very similar values of kinetic energy, independent of 
$L_0$, as seen in figure 2. The difference in jet strength is also easy to understand because jet speed is proportional to the PV jump across the front (see, e.g. Nycander et al. Reference Nycander, Dritschel and Sutyrin1993). Since the maximum and minimum PV are 
$\approx \pm 4{\rm \pi}$ in the flows shown in figures 4 and 5, the PV jumps across the fronts are larger on average when there are fewer fronts, yielding stronger jets.
Figure 5. Potential vorticity (a–c) and kinetic energy (d–f) for simulations with 
$L_D=1/40$ and 
$L_0/L_D = 8$ (a,d, 
$t=999 \, 900$), 
$L_0/L_D = 2.67$ (b,e, 
$t = 874 \, 200$) and 
$L_0/L_D = 1.33$ (c,f, 
$t = 433 \, 600$). The domain length 
$L = 2 {\rm \pi}$ is 40 deformation wavelengths 
$2{\rm \pi} L_D$.
The trends seen in the simulations with 
$L_D = 1/40$, whose fields are shown in figure 5, are similar: more mixed levels fill more of the domain and frontal jets are weaker for the larger scale initial conditions (
$L_0/L_D=2.67$, b,e). For the smaller scale initial conditions (
$L_0/L_D = 1.33$, c,f), the mixed regions are separated into island-like vortices surrounded by a sea of near-zero PV and bounded by strong jets.
For the simulation with 
$L_D = 1/40$ and 
$L_0/L_D = 8$ shown in plots (a,d) in figure 5, the mixed regions and jets are not yet well defined, though the PV contours are gathering together in bundles foreshadowing the emergence of fronts. This simulation has the most slowly evolving flow, and despite the fact that there is an inverse cascade of potential energy, a PV staircase is only just starting to emerge by the end of the simulation. Hence, well-mixed regions of PV and frontal jets coinciding with sharp fronts have not yet had time to emerge.
5. Potential vorticity staircases
5.1. Jet spacing and Rhines scale
Theoretical arguments link the jet spacing, which is equivalent to the staircase structure, in zonal flows on the 
$\beta$ plane, to the Rhines scale 
$L_{Rh} = \sqrt {U/\beta }$, where 
$U$ is a characteristic flow speed, which may be taken as either the root-mean-square (r.m.s.) speed or the jet speed (Rhines Reference Rhines1975; Dritschel & McIntyre Reference Dritschel and McIntyre2008). There are important differences in the present system as compared with the extensively studied system of zonal jets on the 
$\beta$ plane. First, in the latter system, merger between mixed regions must be associated with the loss of a step from the staircase, increasing the jet spacing. By contrast, when 
$\beta = 0$, the staircase structure and jet spacing are not equivalent: the staircases studied here are fixed once they emerge, even though the jet spacing continues to grow as mixed PV regions merge. The staircase structure does not change (other than sharpening) because the area occupied by each mixed PV level is fixed after the staircase emerges.
Since 
$\beta = 0$ in the present system, the Rhines scale 
$L_{Rh}$ also cannot be defined in the same way. If 
$L_{Rh}$ is associated with a fixed staircase structure (not with an evolving jet spacing) then we should expect it to be fixed once the staircase structure is established. On the other hand, if 
$L_{Rh}$ is associated with the evolving jet spacing, which increases as the mixed regions merge, then we should expect it to grow as well. Since the jet speed emerges along with the staircase structure and is not expected to change significantly after the staircase is established, changes in 
$L_{Rh}$ cannot be due to changes in the jet speed. Rather, they must be associated with reductions in an effective vorticity gradient, 
$\beta _{eff}$, which we can define as
\begin{equation} \beta_{eff} = \frac{q_{max} - q_{min}}{L(t)}, \end{equation}
where 
$L(t)$ is the distance between PV maxima and minima, which grows as the mixed regions grow through merger; 
$\beta _{eff}$ decreases in time as the length scale of the flow increases, leading to increases in 
$L_{Rh}$.
Turning to the simulation results, the evolution towards a PV staircase is illustrated for the case 
$L_D=1/25$, 
$L_0/L_D = 5$ in figure 6, which shows the area 
$A$ occupied by PV values greater than a given threshold 
$q_{thr}$, which is allowed to vary from 
$-$12 to 12 in increments of 
$dq$, where 
$dq$ is the fixed increment in PV from one contour to the next. This case starts in the AM regime (large 
$L_0/L_D$) yet nonetheless establishes a distinct staircase structure at late times. In this diagnostic, there is little evolution apart from sharpening of the staircase. The number of steps does not change. The flow however remains highly dynamic (see (a,d) of figure 4), and the kinetic energy continues to decay (see figure 2). Vortices continue to merge and PV fronts reduce in curvature, albeit at a diminishing rate, in a self-similar manner (Burgess & Dritschel Reference Burgess and Dritschel2019; Burgess Reference Burgess2020).
Figure 6. Formation of a PV staircase in the simulation with 
$L_D = 1/25$ and 
$L_0/L_D = 5$ at the times given. Here 
$A(q) \ge q_{thr}$ is the area occupied by PV exceeding a threshold value 
$q_{thr}$.
As 
$L_0/L_D$ decreases, the area of the central step (or steps) in the PV staircase increases, as shown in figure 7. This implies that the mixing of low-level PV is more vigorous and extensive as 
$L_0/L_D$ decreases. The mixed area, and indeed the entire staircase structure, appears to depend only on the ratio 
$L_0/L_D$, as can be seen in figure 8, which shows staircases for 
$L_0/L_D = 1$ with 
$L_D = 1/25$ (a) and 
$L_D = 1/40$ (b). Notably, the PV levels that emerge depend only on 
$L_0/L_D$ (for 
$L_D$ small enough to be able to ignore finite domain effects). The staircases in figure 7 are better established for 
$L_D=1/25$ (a–c) due to the significantly faster pace of evolution; longer integrations with 
$L_D=1/40$ as in (d–f) and 
$L_0/L_D$ large are likely to produce similarly sharp staircases.
Figure 7. Potential vorticity staircases for simulations with 
$L_D = 1/25$ and 
$L_0/L_D = 5$ (a, 
$t = 200 \, 000$), 
$L_0/L_D = 2.5$ (b, 
$t = 200 \, 000$) and 
$L_0/L_D = 1.67$ (c, 
$t = 200 \, 000$); and 
$L_D = 1/40$ and 
$L_0/L_D = 8$ (d, 
$t = 999 \, 900$), 
$L_0/L_D = 4$ (e, 
$t = 605 \, 200$) and 
$L_0/L_D = 2.67$ (f, 
$t = 447 \, 400$), respectively.
Figure 8. Potential vorticity staircases for simulations with 
$L_D = 1/25$, 
$L_0/L_D = 1$ (a, 
$t = 200 \, 000$) and 
$L_D = 1/40$, 
$L_0/L_D = 1$ (b, 
$t = 604 \, 500$). The structures are the same, suggesting the staircase structure depends only on the ratio 
$L_0/L_D$.
5.2. A model for the staircase structure
In developing its staircase structure, the system must interpolate between the maximum and minimum PV values, which are protected from dissipation by their strength and fixed due to material advection of PV, while conserving total energy 
$E = K +P$ and keeping the domain integral of PV zero. As the flow evolves, the energy 
$E$ is increasingly dominated by the potential energy, 
$E \approx P$. If at late times we have 
$n$ mixed levels with areas 
$A_1$, 
$A_2$,…, 
$A_n$ and corresponding values of the streamfunction 
$\psi _1,\psi _2,\ldots, \psi _n$, then
\begin{equation} E = E_0 \approx \tfrac{1}{2}k_D^2 ( \psi_1^2 A_1 + \psi_2^2 A_2 + \cdots+ \psi_n^2 A_n ), \end{equation}
where 
$E_0$ is the fixed initial energy and we have neglected the unmixed PV in the vortex cores.
In the mixed regions, which cover most of the domain, the PV 
$q = \nabla ^2 \psi - k_D^2 \psi$ is dominated by the second term, with relative vorticity 
$\nabla ^2 \psi$ being concentrated on the flanks of the jets bounding the mixed regions. We may therefore approximate 
$q \approx -k_D^2 \psi$, and the domain integral is
\begin{equation} \int_{\mathcal{D}} q \,{\rm d}A = 0 \approx{-}k_D^2(\psi_1 A_1 +\psi_2 A_2 + \cdots+ \psi_n A_n). \end{equation} In predicting the staircase structure we are trying to solve for 
$n$, and for the 
$A_i$ and 
$\psi _i$ given the initial conditions and these constraints. At this point, to progress further, we must make the following simplifying assumptions.
• Firstly, we assume a staircase with
$n$ equally sized steps, i.e. that the mixed regions on which $q \ne 0$ all have the same area $A$.• We also assume that the PV jumps
$\Delta q$ between mixed regions take on a fixed constant value
(5.4)\begin{equation} \Delta q = 2 \frac{q_{m+}}{n}, \end{equation}
where 
$q_{m+}$ is the maximum value taken on by mixed PV, with 
$q_{m-}$ being the corresponding minimum mixed PV value. Note that 
$q_{m+} < q_{max}$ and 
$q_{m-} > q_{min}$ because the extremal values of PV occur within small strong vortices, which are protected from mixing.
Substituting 
$A_i = A$ and expressing the 
$\psi_i$ in terms of 
$\Delta q$ in (5.3) just yields 
$q_{m+} = -q_{m-}$, which is true for a staircase that is symmetric about 
$q=0$. Making the same substitutions into (5.2) yields, after some algebra,
\begin{equation} \frac{6 k_D^2 E_0}{A q_{m+}^2} = n + 3 + \frac{2}{n}. \end{equation} Note that examining the staircases that emerge in the present system shows that as 
$L_0/L_D$ decreases, a large step emerges at 
$q \approx 0$, and the area of this step increases as 
$L_0/L_D$ decreases, while the sizes of the steps on which 
$q \ne 0$ remain roughly the same. This large step at 
$q \approx 0$ does not contribute to the energy or to the domain integral of PV, so we can omit it from the step count, and take the area 
$A$ of the 
$n$ steps with non-zero PV values to be a constant.
Despite the approximations made, the predicted relationship between 
$E_0$ and the measured quantities 
$n$, 
$q_{m+}$ and 
$A_{av}$, where 
$A_{av}$ is the typical step area defined in equation (5.7) below, holds very well, and in fact holds almost exactly for the simulations in which the non-zero mixed PV levels completely fill the domain, as shown in figure 9. These simulations are with 
$L_D = 1/25$ and 
$L_0/L_D = 5$ (open black circle falling on the red line), 
$L_D = 1/40$ and 
$L_0/L_D = 8$ (upper open blue square falling on the red line), and 
$L_D = 1/40$ with 
$L_0/L_D = 4$ (lower open blue square falling on the red line). The agreement becomes progressively worse as 
$L_0/L_D$ decreases and the 
$q \approx 0$ mixed level grows in area. The values of 
$q_{m+}$ used are 
$q_{m+} = 5.75$ for 
$L_D=1/25$ and 
$q_{m+} = 5.25$ for 
$L_D=1/40$, and these are estimated from the staircase structures. The number of steps was counted using an automated procedure, which finds the area 
$A_2$ of the second-largest step (or third, in the special case of 
$L_0 = 1/10$ – see figure 7b, where there are two large steps at low values of PV), and then calculates 
$n$ as
\begin{equation} n = \sum_{A_i \leq A_2} \frac{A_i}{A_2} + \sum_{A_i > A} 1, \end{equation}
and finally subtracts one from the step count if the largest step is located at 
$q \approx 0$. This procedure ensures that tiny steps at higher absolute values of PV are neither completely ignored nor counted as full steps, and that larger steps at low values of PV (such as appear for the case 
$L_0 = 1/10$) are counted as single steps. It also omits the energetically irrelevant step at 
$q \approx 0$ from the step count. The typical area 
$A_{av}$ is calculated as
\begin{equation} A_{av} = \frac{A_{tot}}{n}, \end{equation}
where 
$A_{tot}$ is the total area occupied by steps with non-zero values of mixed PV (except in the case 
$L_D=1/25$, 
$L_0=1/10$, where both large steps with weak values of PV are omitted from both the step count n and the calculation of 
$A_{tot}$). Figure 10 shows that 
$A$ is nearly independent of 
$k_0/k_D$, motivating use of the average step area 
$A_{av}$ (across 
$k_0/k_D$) in figure 9. Hence, the simple assumptions of equal areas between PV steps, and equal PV jumps at each step, provide a good estimate of the number of steps as a function of the initial energy 
$E_0$.
Figure 9. Total energy 
$E_0$ as a function of the number of steps 
$n$ with 
$q \ne 0$, typical non-zero step size 
$A_{\mbox {av}}$, and maximum mixed PV value 
$q_{m+}$ for simulations with 
$L_D = 1/25$ (black open circles) and 
$L_D=1/40$ (blue open squares). The red line corresponds to the theoretical prediction, (5.5).
Figure 10. Measured values of 
$A$ together with 
$A_{\mbox {av}}$ (dashed lines) for 
$k_D=25$ (black) and 
$k_D=40$ (blue).
The area of the central step or steps, whose size reflects the degree of low-level PV mixing, decreases approximately linearly as 
$L_0/L_D$ increases, at least up to 
$L_0/L_D \approx 2.4$; see figure 11(b) which shows the largest step size (jump in area) 
$\Delta A_{max}$ (or 
$\Delta A_0$ in the case that the largest step falls at 
$q \approx 0$) normalized by the domain area 
$\mathcal{D}$ as a function of 
$L_0/L_D$. The flow is evidently more vigorous, especially at early times, at smaller values of 
$L_0/L_D$. Indeed, the initial shear, as measured, for example, by the r.m.s. vorticity, increases as 
$L_0/L_D$ decreases; see figure 11(a). There is a striking correlation between this shear and the area of the largest step. Another way to see this is provided in figure 12, showing how the largest step in PV, 
$\Delta q_1$, varies with the initial shear. There is an almost linear relationship, consistent with simple models of the effects of shear on vortices, where shear on the order of a tenth of the vortex PV is sufficient to overwhelm the vortex and elongate it indefinitely (Dritschel Reference Dritschel1990; Dritschel & McIntyre Reference Dritschel and McIntyre2008). Here, at early times, the shear mixes the weak PV levels until the resulting step in PV is sufficient to withstand the shear, which decays as the flow evolves.
Figure 11. Initial shear, i.e. r.m.s. vorticity as a function of 
$L_0/L_D$ (a) and largest step size (jump in area) 
$\Delta A_{max}$ (or 
$\Delta A_0$ in the case that 
$q \approx 0$ on the largest mixed area) normalized by the domain area 
$\mathcal{D}$ as a function of 
$L_0/L_D$ (b) for simulations with 
$L_D = 1/25$ (open circles) and 
$L_D = 1/40$ (open squares).
Figure 12. Potential vorticity jump 
$\Delta q_1$ bounding region of zero PV plotted against shear (initial r.m.s. vorticity) for 
$L_D=1/25$ (black) and 
$L_D=1/40$ (blue).
6. Conclusions
This paper has examined the late-time evolution of freely evolving quasi-geostrophic shallow water turbulence in the inviscid limit. We have focused on the less studied regime in which the Rossby deformation length, 
$L_D$, the fundamental scale above which flow interactions are strongly suppressed, is much smaller than the domain scale. Under these conditions, the flow evolution proceeds much more slowly than in classical 2-D turbulence (
$L_D\to \infty$, see Dritschel et al. Reference Dritschel, Scott, Macaskill, Gottwald and Tran2008). We have considered how the initial characteristic scale of the flow 
$L_0$, in relation to 
$L_D$, affects the flow structures which emerge at sufficiently late times.
For large 
$L_0/L_D$, it was theorised that the flow remains large scale (of characteristic scale 
$L(t) \gg L_D$) for all time (Larichev & McWilliams Reference Larichev and McWilliams1991; Iwayama et al. Reference Iwayama, Shepherd and Watanabe2002). Our simulations show, however, that small-scale frontal structures eventually emerge – sharp PV fronts co-located with jets of width 
$O(L_D)$ – and these subsequently dominate the evolution. Despite the slow initial evolution when 
$L_0/L_D\gg 1$, smaller scales inevitably develop as PV is mixed by Rossby wave breaking. While this is weak, when it is sustained for long enough, the mixing becomes so complete that near staircases form in the PV: these take the form of plateaus of near uniform PV separated by virtual discontinuities. These staircases form in separate islands, which slowly move through a sea of near-zero PV. Near the centre of each island there is a small vortex, also of width 
${O}(L_D)$. These vortices contain the highest-magnitude PV values in the flow.
The formation of PV staircases occurs for all 
$L_0/L_D$ studied, down to 
$L_0/L_D=1$, but likely for all smaller values as well since the initial mixing is even more vigorous in these cases. The key difference between different values of 
$L_0/L_D$ is the form of the final staircase, specifically the area occupied by PV values below a given threshold. Large 
$L_0/L_D$ results in a staircase with nearly equal steps in PV and in area. Smaller 
$L_0/L_D$ results in greater mixing of low-level PV, resulting in a larger area of mixed PV close to the zero value. This mixing was associated with the initial shear: the region of mixed PV is bounded by PV jumps proportional to the initial shear, which grows as 
$L_0/L_D$ reduces.
Finally, we developed a simple theory to describe the PV staircase which emerges, and how it depends on 
$L_0/L_D$. The theory predicts the average number of steps, the average area of the PV levels (excluding the zero PV level) and the associated average PV jump, purely in terms of the initial energy of the flow.
In this work we have ignored the effects of a background planetary vorticity gradient 
$\beta$, which has historically been associated with zonal jet formation and inhomogeneous PV mixing (Rhines Reference Rhines1975; Dritschel & McIntyre Reference Dritschel and McIntyre2008; Galperin & Read Reference Galperin and Read2019). Here, even with 
$\beta =0$, jets emerge and come to dominate the flow, and the PV mixing can be so complete as to form a near PV staircase – without forcing. Small 
$L_D$ facilitates this process by sustaining large-amplitude waves on scales 
$L>L_D$ which break and mix, or merge during collisions with other vortices. This wave breaking and vortex merger continues indefinitely and self-similarly (in an unbounded domain), progressively sharpening the PV staircase. It appears inevitable that the staircase becomes ‘perfect’ as 
$t\to \infty$, i.e. the PV becomes piecewise uniform, except perhaps in the vortices occupying the island centres.
While the model examined in this paper is a great simplification of naturally occurring geophysical flows, it nonetheless provides new insight into the way tracers like PV can mix in turbulent flows. Most previous studies have focused on the simpler barotropic limit 
$L/L_D\to 0$, in the presence of a background planetary vorticity gradient 
$\beta$, or their multi-layer (predominantly two-layer) analogues. Here, we have focused on the less well-examined situation in which 
$L/L_D\gg 1$, and when 
$\beta =0$. One might have expected no jet formation in this case, but in fact it is widespread. Instead of a background planetary vorticity gradient, local variations in PV provide the necessary ingredient for inhomogeneous PV mixing and consequent jet formation, seen here as PV staircases within slowly drifting and weakly interacting islands of anomalous PV. In the barotropic limit, mixing weakens and all but ceases without some forcing mechanism to maintain mixing: freely decaying flows tend to form weak jets without sharp PV variations. Here, by contrast, when 
$L/L_D\gg 1$, freely decaying flows continue to mix so long as the PV anomalies are permitted to grow through occasional mergers. This alone is adequate for steepening PV gradients, and for maintaining their steepness. No forcing or damping is required. In this sense, the present model is the simplest model capable of jet formation. It isolates the essential mechanism: inhomogeneous PV mixing.
Acknowledgements
BHB acknowledges support for this research from a Leverhulme Trust Early Career Fellowship.
Declaration of interests
The authors report no conflict of interest.
Open access
