Effect of Prandtl’s ratio on balance in geophysical turbulence

Abstract

The fluid dynamics of the atmosphere and oceans is to a large extent controlled by the slow evolution of a scalar field called ‘potential vorticity’ (PV), with relatively fast motions such as inertia-gravity waves playing only a minor role. This state of affairs is commonly referred to as ‘balance’. Potential vorticity is a special scalar field which is materially conserved in the absence of diabatic effects and dissipation, effects that are generally weak in the atmosphere and oceans. Moreover, in a balanced flow, PV induces the entire fluid motion and its thermodynamic structure (Hoskins et al., Q. J. R. Meteorol. Soc., vol. 111, 1985, pp. 877–946). While exact balance is generally not achievable, it is now well established that balance holds to a high degree of accuracy in rapidly rotating and strongly stratified flows. Such flows are characterised by both a small Rossby number, $\mathit{Ro}\equiv |{\it\zeta}|_{max}/f$, and a small Froude number, $\mathit{Fr}\equiv |{\bf\omega}_{h}|_{max}/N$, where ${\it\zeta}$ and ${\bf\omega}_{h}$ are the relative vertical and horizontal vorticity components, while $f$ and $N$ are the Coriolis and buoyancy frequencies. In fact, balance can even be a good approximation when $\mathit{Fr}\lesssim \mathit{Ro}\sim \mathit{O}(1)$. In this study, we examine how balance depends specifically on Prandtl’s ratio, $f/N$, in unforced freely evolving turbulence. We examine a wide variety of turbulent flows, at a mature and complex stage of their evolution, making use of the fully non-hydrostatic equations under the Boussinesq and incompressible approximations. We perform numerical simulations at exceptionally high resolution in order to carefully assess the degree to which balance holds, and to determine when it breaks down. For this purpose, it proves most useful to employ an invariant PV-based Rossby number ${\it\varepsilon}$, together with $f/N$. For a given ${\it\varepsilon}$, our key finding is that – for at least tens of characteristic vortex rotation periods – the flow is insensitive to $f/N$ for all values for which the flow remains statically stable (typically $f/N\lesssim 1$). Only the vertical velocity varies in proportion to $f/N$, in line with quasi-geostrophic (QG) scaling for which $\mathit{Fr}^{2}\ll \mathit{Ro}\ll 1$. We also find that as ${\it\varepsilon}$ increases towards unity, the maximum $f/N$ attainable decreases towards 0. No statically stable flows occur for ${\it\varepsilon}\gtrsim 1$. For all stable flows, balance is found to hold to a remarkably high degree: as measured by an energy norm, imbalance never exceeds more than a few per cent of the balance, even in flows where $\mathit{Ro}>1$. The vertical velocity $w$ remains a tiny fraction of the horizontal velocity $\boldsymbol{u}_{h}$, even when $w$ is dominantly balanced. Finally, typical vertical to horizontal scale ratios $H/L$ remain close to $f/N$, as found previously in QG turbulence for which $\mathit{Fr}\sim \mathit{Ro}\ll 1$.

1. Introduction

The dynamics of the Earth’s atmosphere and oceans is significantly influenced by both the planetary rotation and density stratification. These two features play vital but contrary roles in the large-scale fluid dynamics on Earth (and very likely in other planetary atmospheres, Read (2011)). Rotation alone tends to form deep two-dimensional flows which vary weakly in the vertical (known as ‘Taylor columns’, Taylor (1923)), whereas stratification induces shallow or layered flows having strong variations across stratification surfaces and motion parallel to these surfaces (Riley & Lelong 2000). Their combination leads to the formation of structures that are coupled over a vertical scale approximately given by $H\sim fL/N$ , where $H$ and $L$ are typical vertical and horizontal length scales respectively, $f$ is the Coriolis frequency associated with the Earth’s rotation rate and $N$ is the buoyancy frequency (Charney 1971; Herring 1980; Dritschel, de la Torre Juárez & Ambaum 1999; Reinaud & Dritschel 2002; Reinaud, Dritschel & Koudella 2003; Praud, Sommeria & Fincham 2006). Over much of the atmosphere and oceans, $f$ tends to be small compared with $N$ , and therefore the ratio $f/N$ , known as ‘Prandtl’s ratio’, tends to be small (typically around $10^{-2}$ in the atmosphere and $10^{-1}$ in the oceans, see Gill (1982) and Vallis (2008)). However, this does not imply that the effects of rotation are less important than those of stratification. Indeed, when the scale ratio $H/L\sim f/N$ , both effects contribute comparably.

One of the mathematical consequences of rotation and stratification is that certain terms in the equations of motion tend to dominate other terms. For fast rotation, the Coriolis acceleration and horizontal pressure gradient dominate the horizontal acceleration. For strong stratification, the buoyancy and vertical pressure gradient dominate the vertical acceleration. If one omits the acceleration entirely, the resultant balances are called ‘geostrophic’ and ‘hydrostatic’ respectively. Together, they are called ‘thermal-wind’ balance (see e.g. Vallis 2008). This underlying balance forms the basis for the ‘quasi-geostrophic’ (QG) model of geophysical fluid dynamics (Charney 1948), valid when $\mathit{Fr}^{2}\ll \mathit{Ro}\ll 1$ , where $\mathit{Fr}\equiv |{\bf\omega}_{h}|_{max}/N$ is the Froude number and $\mathit{Ro}\equiv |{\it\zeta}|_{max}/f$ is the Rossby number (ignoring here effects of variable planetary vorticity). In these expressions, ${\bf\omega}_{h}$ and ${\it\zeta}$ are the horizontal and vertical parts of the relative vorticity (relative to the planetary vorticity $f$ ). Notably, when $\mathit{Fr}\sim \mathit{Ro}$ rotation and stratification contribute comparably to the potential vorticity (PV), ${\bf\omega}_{a}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\rho}$ , where ${\bf\omega}_{a}$ is the absolute vorticity (including the planetary vorticity) and ${\it\rho}$ is the density. Moreover, when $\mathit{Fr}\sim \mathit{Ro}$ typical vertical to horizontal scale ratios $H/L$ are comparable to Prandtl’s ratio $f/N$ .

The QG model in fact only approximates a more fundamental higher-order (in $\mathit{Ro}$ ) balance widely exhibited by geophysical flows (see Baer & Tribbia 1977; Leith 1980; McWilliams & Gent 1980; Vallis 1996; Bokhove 1997; Olsson & Cotton 1997; Muraki, Snyder & Rotunno 1999; Dritschel & Viúdez 2007; McKiver & Dritschel 2008; Vanneste 2013 for a sample of the vast literature on the subject). The QG model makes use of only the leading-order thermal-wind balance, and closes the asymptotic system of equations at $\mathit{O}(\mathit{Ro}^{2})$ , where one obtains the dynamical evolution equation for an approximation of PV. The horizontal velocity, buoyancy and pressure fields are obtained by linear ‘PV inversion’, and only at first order in $\mathit{Ro}$ . Moreover, at this order, the velocity is entirely horizontal, resulting in layerwise-2D motion. The QG model thus uses the simplest balance possible. This is an attractive feature which has led to the model’s widespread use in studying basic aspects of atmospheric and oceanic fluid dynamics. However, many previous studies have demonstrated that balance may run much deeper than what is explicitly used in the QG model (see Ford, McIntyre & Norton 2000; Mohebalhojeh & Dritschel 2001; Mohebalhojeh 2002; Dritschel & Viúdez 2003; Viúdez & Dritschel 2004; McKiver & Dritschel 2008 and references therein). This is important as it enables one to better separate ‘balanced’ motions directly arising from PV from all other ‘imbalanced’ motions associated with inertia-gravity waves (IGWs). Thereby, one can more accurately quantify the generation of IGWs in a general flow and assess the degree of balance.

The QG approximation itself provides a balanced estimate of the vertical velocity at $\mathit{O}(\mathit{Ro}^{2})$ , even though this is not used in the material transport of PV. A variety of more complete balance models, including all terms at $\mathit{O}(\mathit{Ro}^{2})$ , can be found in the above cited works. Here, we quantify the degree of balance using nonlinear quasi-geostrophic (NQG) balance (McKiver & Dritschel 2008) and optimal potential vorticity (OPV) balance (Viúdez & Dritschel 2004), both developed for the non-hydrostatic rotating Boussinesq equations. These methods are unique in using the unapproximated form of the PV to improve the estimate of balance. An accurate representation of PV is crucial, as previously shown by Mohebalhojeh & Dritschel (2000) in the shallow-water context.

In this paper, we examine freely evolving rotating stratified turbulence and investigate how balance and, generally, the flow evolution depend on both the Rossby number $\mathit{Ro}$ and Prandtl’s ratio $f/N$ . The closest previous work is that of Praud et al. (2006), who carried out a comprehensive experimental examination of freely decaying rotating stratified turbulence in the large rotating tank in Grenoble. They were able to study the regime $f/N<1$ , which is typically difficult experimentally. They quantified a wide range of flow properties but were not able to directly infer the impact of IGWs. Here, we use a highly accurate numerical approach, designed for stably stratified flows (Dritschel & Viúdez 2003), which enables us to quantify the impact of IGWs as well as to study in detail the balanced vortical part of the flow. In particular, our approach allows us to control the initial IGW activity, and indeed minimise it in order to study spontaneous adjustment emission over the course of the flow evolution. In this way, our work is complementary to that of Praud et al. (2006), although our key findings are in agreement.

Much previous research has focused on forced turbulence (see e.g. Métais et al. 1996; Smith & Waleffe 2002; Waite & Bartello 2006; Bartello 2010; Molemaker, McWilliams & Capet 2010; Deusebio, Vallgren & Lindborg 2013). This is convenient for reaching a statistically steady state when numerical or molecular-like dissipation is draining energy away. However, forcing is problematic as regards balance. It is difficult to apply forcing that does not strongly excite IGWs, even when the forcing is, say, in geostrophic balance. Rotating stratified flows at small Rossby number may exhibit much higher-order balance (Dritschel & Viúdez 2007; McKiver & Dritschel 2008; Tsang & Dritschel 2015), which can be completely masked by the forcing. For this reason, we focus on the unforced case.

A pivotal early work in this context was carried out by Bartello (1995), who developed a way of decomposing a flow into ‘geostrophic’ and ‘ageostrophic’ modes, thereby enabling one to study their interactions and the resulting spectral energy cascades. Because the decomposition is based on the linearised equations about a state of rest, the ‘ageostrophic’ modes may in fact be dominantly balanced (as indeed we find below). We therefore prefer here to decompose a flow into a higher-order balanced state and the residual imbalance. This is especially important when studying carefully balanced initial flow states, in which geostrophic adjustment is minimal. When initial states are not well balanced, Bartello (1995) demonstrates that a significant forward cascade of ‘ageostrophic’ energy occurs, presumably in the form of steepening and breaking IGWs. No such forward cascade is observed here.

The plan of the paper is as follows. In § 2 we describe how we initialise the flow – in a state close to balance – from an existing complex turbulent QG solution at a mature stage of evolution. We briefly describe the numerical model, together with the NQG and OPV balance procedures used for diagnosing balance. Then in § 3 we present key characteristics of the flow evolution, including spectra, vertical velocity, balance and imbalance. We round off the paper in § 4 with a summary and discussion of the main results.

2. Problem formulation

2.1. Non-hydrostatic equations

Here we consider the non-hydrostatic (NH) equations under the Oberbeck–Boussinesq approximation, where the density ${\it\rho}$ varies weakly from a mean background value, ${\it\rho}_{0}$ , i.e.

(2.1) $$\begin{eqnarray}{\it\rho}(\boldsymbol{x},t)={\it\rho}_{0}+{\it\varrho}_{z}z+{\it\rho}^{\prime }(\boldsymbol{x},t),\end{eqnarray}$$

where ${\it\varrho}_{z}z$ is the mean linear density ( ${\it\varrho}_{z}<0$ is a constant) and ${\it\rho}^{\prime }(\boldsymbol{x},t)$ is the anomalous density. The NH equations are obtained by neglecting terms of $\mathit{O}(({\it\rho}-{\it\rho}_{0})^{2}/{\it\rho}_{0}^{2})$ , giving

(2.2a ) $$\begin{eqnarray}\displaystyle & \dot{\boldsymbol{u}}+f\boldsymbol{k}\times \boldsymbol{u}=-{\it\rho}_{0}^{-1}\boldsymbol{{\rm\nabla}}{\it\Phi}+b\boldsymbol{k}, & \displaystyle\end{eqnarray}$$

(2.2b ) $$\begin{eqnarray}\displaystyle & {\dot{b}}+N^{2}w=0, & \displaystyle\end{eqnarray}$$

(2.2c ) $$\begin{eqnarray}\displaystyle & \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}=0, & \displaystyle\end{eqnarray}$$

where $\boldsymbol{u}\equiv (u,v,w)$ is the three-dimensional velocity field, $(\dot{~})=\text{D}(~)/\text{D}t=(~)_{t}+\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}(~)$ denotes the material time derivative (in the rotating frame), ${\it\Phi}$ is the geopotential, $b\equiv -g{\it\rho}^{\prime }/{\it\rho}_{0}$ is the buoyancy, $g$ is gravity and $\boldsymbol{k}$ denotes the vertical unit vector; $N$ is the mean buoyancy frequency defined by $N^{2}\equiv -g{\it\varrho}_{z}/{\it\rho}_{0}$ .

Here, we consider the NH equations in the form introduced in Dritschel & Viúdez (2003). There, the equations of motion are recast to explicitly use material conservation of PV, denoted ${\it\Pi}$ , as well as a pair of variables – the two components $\mathscr{A}$ and $\mathscr{B}$ of the ageostrophic horizontal vorticity $\boldsymbol{A}_{h}$ – representing the leading-order departure from QG balance. The dimensionless Rossby–Ertel PV is given by

(2.3) $$\begin{eqnarray}{\it\Pi}=(\boldsymbol{k}+{\bf\omega}/f)\boldsymbol{\cdot }(\boldsymbol{k}-\boldsymbol{{\rm\nabla}}\mathscr{D}),\end{eqnarray}$$

where ${\bf\omega}$ is the vorticity, $\boldsymbol{k}$ is a unit vector in the vertical direction and $\mathscr{D}$ is the isopycnal displacement defined by $\mathscr{D}=-b/N^{2}$ , where $b$ is the buoyancy anomaly (the departure from the mean buoyancy, $N^{2}z$ ). The dimensionless ageostrophic horizontal vorticity is given by

(2.4) $$\begin{eqnarray}\boldsymbol{A}_{h}={\bf\omega}_{h}/f+\boldsymbol{{\rm\nabla}}_{h}b/f^{2},\end{eqnarray}$$

where $\text{}_{h}$ denotes ‘the horizontal components of’. Thermal-wind balance corresponds to $\boldsymbol{A}_{h}=0$ , to a high degree of accuracy (Dritschel & Viúdez 2003). These variables, ${\it\Pi}$  and $\boldsymbol{A}_{h}$ , can be inverted using a vector potential ${\bf\varphi}$ , in terms of which the velocity and displacement fields are given by $\boldsymbol{u}=-f\boldsymbol{{\rm\nabla}}\times {\bf\varphi}$ and $\mathscr{D}=-f^{2}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }{\bf\varphi}/N^{2}$ (details may be found in Dritschel & Viúdez 2003).

The evolution equations for ${\it\Pi}$ and $\boldsymbol{A}_{h}$ are

(2.5) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{D}{\it\Pi}}{\text{D}t}=0, & \displaystyle\end{eqnarray}$$

(2.6) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{D}\boldsymbol{A}_{h}}{\text{D}t}+f\boldsymbol{k}\times \boldsymbol{A}_{h}=f^{-1}({\bf\omega}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}})\boldsymbol{u}_{h}+(1-N^{2}/f^{2})\boldsymbol{{\rm\nabla}}_{h}w+(\boldsymbol{{\rm\nabla}}_{h}\boldsymbol{u})\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\mathscr{D}, & \displaystyle\end{eqnarray}$$

where $\text{D}/\text{D}t=\partial /\partial t+\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}$ is the material derivative, $\boldsymbol{u}_{h}=(u,v)$ is the horizontal velocity and $\boldsymbol{{\rm\nabla}}_{h}$ is the horizontal gradient operator. Notably, the QG equations are recovered by setting $\boldsymbol{A}_{h}=0$ and approximating the PV by ${\it\Pi}=1+{\it\zeta}/f-\partial \mathscr{D}/\partial z$ . Then the vertical velocity $w=0$ , enabling one to recover $\boldsymbol{u}_{h}$ and $\mathscr{D}$ from a scalar potential (streamfunction) ${\it\phi}$ determined by inverting Poisson’s equation ${\rm\nabla}_{h}^{2}{\it\phi}+(f^{2}/N^{2})\partial ^{2}{\it\phi}/\partial z^{2}={\it\Pi}-1$ .

In the numerical method, ${\it\Pi}$ is represented by contours on isopycnal surfaces ( $z-\mathscr{D}=$ constant) and evolved using the contour-advective semi-Lagrangian (CASL) algorithm (Dritschel & Ambaum 1997; Dritschel & Viúdez 2003). Using the PV as a prognostic variable has the advantage of exploiting the underlying balance, since PV is inextricably linked with the balanced dynamics (see Dritschel & Viúdez 2007; McKiver & Dritschel 2008). The remaining variables, the components of $\boldsymbol{A}_{h}$ , are evolved conventionally, on a grid. Note, in typical geophysical flows close to a (thermal-wind) balanced state, these components make only a minor contribution to the overall dynamics (and they may contain, in part, imbalanced motions). Full details of the numerical method may be found in appendix 1 of Dritschel & Viúdez (2003).

2.2. Initialisation and parameter settings

The turbulent initial condition used here is the same as previously used in McKiver & Dritschel (2008). Nearly 200 simulations were conducted to examine the degree of balance and other characteristics of rotating stratified turbulence. We widely varied both Prandtl’s ratio $f/N$ and the PV-based Rossby number ${\it\varepsilon}\equiv |{\it\Pi}-1|_{max}$ , which, unlike $|{\it\zeta}|_{max}/f$ , is constant in freely evolving flows.

This non-standard choice for the Rossby number enables us here to clearly delineate the domain of dominantly balanced flows in the ${\it\varepsilon}$ – $f/N$ parameter space. Using $U/fL$ for chosen velocity and length scales $U$ and $L$ is generally ambiguous, and the choice ${\it\zeta}/f$ (extreme or r.m.s. value) is undesirable since it can vary strongly in time. It does not provide a single number characterising the amplitude of an evolving flow. Further justification for ${\it\varepsilon}$ is given below. Note, we are not claiming that ${\it\zeta}/f$ (or the Froude number) is irrelevant as regards balance – we are simply advocating an alternative view of parameter space using parameters that are time invariant.

To begin, a single QG simulation was performed starting from an isotropic PV field (after stretching the vertical coordinates by $N/f$ , as is natural in the QG equations). In this case the PV in QG flow is given by $q(\boldsymbol{x},t)={\rm\nabla}^{2}{\it\psi}$ , where the streamfunction ${\it\psi}$ generates the QG layerwise-2D velocity field $\boldsymbol{u}=(-\partial {\it\psi}/\partial y,\partial {\it\psi}/\partial x,0)$ . The initial QG PV field consists of 500 positive (cyclonic) and 500 negative (anticyclonic) spherical vortices of uniform QG PV, of magnitude $q=\pm 4{\rm\pi}$ , without loss of generality. They are placed randomly, without overlapping, in a triply periodic cube (in stretched coordinates $x$ , $y$ and $Nz/f$ ), and their sizes are chosen from a frequently observed power-law number density distribution (Reinaud et al. 2003), with the volume of the largest vortex being 20 times that of the smallest vortex. Figure 1(a) shows the initial configuration.

Figure 1. Distribution of the QG PV anomaly $q$ at times (a) $t=0$ and (b) $t=40$ . The view is orthographic, $60^{\circ }$ from the vertical, and from the $y{-}z$ plane. From this view we can see the top and front faces of the domain, and a faint white line indicates where these faces intersect. Cyclonic vortices are lightly shaded while anticyclonic vortices are darkly shaded.

In both the QG and NH simulations, a basic grid resolution of $128\times 128\times 128$ was used, in a triply periodic domain of dimensions $2{\rm\pi}N/f\times 2{\rm\pi}N/f\times 2{\rm\pi}$ , anticipating the characteristic scaling $H/L\sim f/N$ (see below). In the QG simulation, $f/N$ is absorbed into the definition of the vertical coordinate $z$ and does not need to be specified. In the CASL algorithm, the PV field is represented by contours on isopycnal surfaces (uniformly spaced in density), using four times as many layers as grid points. For consistency a grid four times finer in each horizontal direction is used to convert PV contours to gridded PV values – this is standard in the CASL algorithm and results in a more accurate representation of the flow associated with PV (Dritschel & Ambaum 1997). Contours are retained down to a 20th of the horizontal grid resolution, below which thin filaments are removed by ‘contour surgery’ (Dritschel 1988) to limit the otherwise near exponential growth in contour length. All other aspects of the algorithm rely on fast Fourier transforms, particularly in the dealiased pseudospectral evolution of the ageostrophic horizontal vorticity $\boldsymbol{A}_{h}$ . While the grid resolution may seem modest, the effective resolution of the CASL algorithm is much higher, more than 10 times in each direction, as demonstrated in Dritschel & Viúdez (2003); see also Dritschel & Scott (2009) and Dritschel & Tobias (2012).

The QG simulation was run for 40 time units, well into the decaying stage characterised by a significant forward cascade of the PV spectral density (enstrophy). The complex QG PV field at this time, $q(\boldsymbol{x},40)$ (see figure 1 b), was then used to generate the initial states of NH simulations at finite PV-based Rossby number ${\it\varepsilon}$ and prescribed $f/N$ . Starting from a state of rest in the NH model, the PV anomaly ${\it\varpi}={\it\Pi}-1$ was slowly ramped from 0 to ${\it\varpi}={\it\varepsilon}q(\boldsymbol{x},40)/4{\rm\pi}$ , for a given PV-based Rossby number ${\it\varepsilon}$ , while holding the distribution of PV fixed (that is, the PV contours were held fixed). Over this ramping period, the full dynamical equations for $\boldsymbol{A}_{h}$ were integrated starting from $\boldsymbol{A}_{h}=0$ . This generates an initial state nearly void of IGWs (see Dritschel & Viúdez 2007; McKiver & Dritschel 2008 and references therein), as long as the ramping period is greater than a few inertial periods, $T_{ip}=2{\rm\pi}/f$ . Here, consistent with previous works, we integrate for $5T_{ip}$ . The initial state thus generated is hereafter referred to as ‘ $t=0$ ’.

Non-hydrostatic simulations for many values of ${\it\varepsilon}$ and $f/N$ were then integrated forwards over the equivalent of 20 QG time units (or more in some cases, see below). To ensure that the IGWs, having frequencies between $f$ and $N$ , were well resolved in time, we used an explicit third-order Adams–Bashforth time stepping procedure with a time step ${\rm\Delta}t=0.025$ to $0.1T_{\mathit{bp}}$ , where $T_{\mathit{bp}}=2{\rm\pi}/N$ is the buoyancy period (smaller ${\rm\Delta}t$ values are required for larger ${\it\varepsilon}$ and $f/N$ ). Moreover, to control the generation of grid-scale noise during the time integration, a weak bi-harmonic hyperdiffusion was added to the $\boldsymbol{A}_{h}$ tendencies. The maximum damping rate on the highest wavenumber in spectral space was taken to be $1+10{\it\varepsilon}^{4}$ per inertial period $T_{\mathit{ip}}=2{\rm\pi}/f$ . As demonstrated in Dritschel & Viúdez (2003), this damping rate is much less than that required in a conventional (contour-free) pseudospectral numerical method.

2.3. Balance diagnosis procedures

To quantify the importance of IGWs in the dynamics of rotating stratified turbulence, we decompose the flow at any time into a balanced component entirely due to PV and a residual imbalanced component, which we attribute to IGWs. In fact, such a decomposition is never exact in a general nonlinear flow (Ford et al. 2000; Vanneste & Yavneh 2004; Vanneste 2013), and as a result it is not possible to define balance precisely. Instead, we must settle for estimates of balance. Such estimates nevertheless can be highly accurate, as the residual imbalance can be shown to exhibit the dispersion characteristics of IGWs (Dritschel & Viúdez 2003, 2007). Here, we estimate balance in two ways, using optimal PV (OPV) balance (Viúdez & Dritschel 2004) and nonlinear QG (NQG) balance (McKiver & Dritschel 2008). Both procedures are distinct from all other procedures in one major respect: the use of the unapproximated form of the Rossby–Ertel PV. That is, the PV is retained at all orders in Rossby number. This enables one to define balance more accurately, as explicit comparisons with procedures using approximated PV have demonstrated (McKiver & Dritschel 2008).

Nonlinear QG balance is in other ways conventional in that it makes use of a specified pair of ‘balance relations’, obtained by eliminating a pair of time derivatives from the governing equations (a concise summary is provided in appendix C of Tsang & Dritschel 2015). In NQG balance, we set $\partial \boldsymbol{A}_{h}/\partial t=0$ since this is $\mathit{O}({\it\varepsilon}^{3})$ for small PV-based Rossby number ${\it\varepsilon}$ . This turns the prognostic equations (2.6) for $\boldsymbol{A}_{h}$ into diagnostic ones, giving $\boldsymbol{u}$ and $\mathscr{D}$ in terms of PV (using incompressibility). In fact, because the exact definition of PV is used, a further condition is required (Mohebalhojeh 2002), namely that the linear part of the PV at $\mathit{O}({\it\varepsilon}^{2})$ vanishes. The resulting system of equations can be solved iteratively, and convergence is obtained even when ${\it\varepsilon}=\mathit{O}(1)$ and $f/N>1$ (see below).

Optimal PV balance (Viúdez & Dritschel 2004) differs significantly from all other balance procedures. Instead of specifying a pair of balance relations, which is arguably ad hoc (regarding which pair is ‘best’), an estimate of balance is obtained by integrating the full equations of motion, apart from the PV (again, a concise summary of the method is provided in appendix D of Tsang & Dritschel 2015). The PV is modified, artificially, by slowly varying it on fluid particles so as to minimise the generation of IGWs. Essentially, OPV balance seeks a ‘rest state’ configuration of fluid particles at some time $t-{\it\Delta}$ in the past that, after evolution through a time ${\it\Delta}$ , arrive at the actual configuration of fluid particles at the diagnostic time $t$ . Over this period ${\it\Delta}$ , the PV anomaly ${\it\varpi}={\it\Pi}-1$ on each fluid particle is ramped from 0 (a state of rest ${\bf\varphi}=0$ ) to its actual value at the diagnostic time $t$ . In practice, finding the rest state configuration at $t-{\it\Delta}$ must be done by a sequence of forwards and backwards integrations, until convergence. The result, if ${\it\Delta}$ is sufficiently long, is a flow with very weak IGW activity (more details and tests can be found in Dritschel & Viúdez (2007)). In this sense, OPV balance attempts to minimise imbalance but does not entirely eliminate it. The procedure depends only on the ramp period ${\it\Delta}$ , which is constrained to be short enough to avoid severe distortion of the PV field (otherwise, the procedure does not converge). In practice, ${\it\Delta}=5T_{ip}$ normally works best in that the resulting imbalance is minimal.

We denote the balanced fields obtained from either NQG or OPV balance by $\boldsymbol{u}_{\mathit{b}}$ and $\mathscr{D}_{\mathit{b}}$ . Then, the imbalanced fields are just the differences from the actual fields at the given diagnostic time $t$ , i.e.  $\boldsymbol{u}_{\mathit{i}}=\boldsymbol{u}-\boldsymbol{u}_{\mathit{b}}$ (note that in what follows all balanced and imbalanced fields are denoted by subscripts $\mathit{b}$ and $\mathit{i}$ respectively). A useful measure of imbalance is the energy norm

(2.7) $$\begin{eqnarray}E_{\mathit{i}}=\langle |\boldsymbol{u}_{\mathit{i}}^{2}|+N^{2}\mathscr{D}_{\mathit{i}}^{2}\rangle \text{}^{1/2},\end{eqnarray}$$

where the angled brackets denote a domain average. The percentage of imbalance is defined by $\%E_{\mathit{i}}=100E_{\mathit{i}}/E$ , where $E$ is defined in the same way as $E_{\mathit{i}}$ but using the full fields.

In previous studies, OPV balance has been generally found to attribute a smaller fraction of a flow to imbalance than other procedures, including NQG balance. In this sense, OPV balance appears to give a more accurate estimate of balance (and by subtraction, imbalance). However, in those studies, only small values of Prandtl’s ratio, in fact $f/N\leqslant 0.1$ , were examined. Here, at larger $f/N$ , the situation reverses, as shown by $E_{\mathit{i}}(t)$ in figure 2, for a flow at a moderate PV-based Rossby number, ${\it\varepsilon}=0.5$ , and a moderate Prandtl ratio, $f/N=0.5$ . Results are presented for various choices of the ramp period ${\it\Delta}$ , but none give a better estimate of the balance than NQG. Evidently, NQG balance is much less sensitive to $f/N$ than OPV balance. Since the focus of this paper is on the effect of moderate to large $f/N$ , in the results presented below we use NQG balance exclusively.

3. Results

3.1. Potential vorticity evolution

We first examine qualitatively how the PV evolution varies with $f/N$ , for a fixed PV-based Rossby number ${\it\varepsilon}$ . Figure 3 compares the end state (after 20 equivalent QG time units) of four simulations, all having ${\it\varepsilon}=0.25$ but widely different Prandtl ratios, $f/N=0.1$ , 0.5, 1.0 and 1.5. Note, only a small portion of the domain is shown to ease comparison. This figure demonstrates clearly that the PV evolution is virtually unaffected by the value of $f/N$ , up to the maximum value of $f/N$ for which we have static stability (see § 3.3 below). Many of the finest details in the PV field are seen at the same places in all four cases.

Figure 2. Evolution of the imbalance energy norm $E_{\mathit{i}}$ for ${\it\varepsilon}=0.5$ and $f/N=0.5$ as determined from NQG balance (bold solid line with triangles), and OPV balance using the ramp periods ${\it\Delta}=1.25T_{\mathit{ip}}$ (thin solid line), ${\it\Delta}=2.5T_{\mathit{ip}}$ (short dashed line), ${\it\Delta}=3.75T_{\mathit{ip}}$ (long dashed line) and ${\it\Delta}=5T_{\mathit{ip}}$ (bold line). The time $t$ in this and subsequent figures is given in equivalent QG time units.

Figure 3. Comparison of the PV anomaly fields ${\it\varpi}$ for various Prandtl ratios $f/N$ ((a) 0.1, (b) 0.5, (c) 1.0, (d) 1.5), for ${\it\varepsilon}=0.25$ and at 20 QG time units. The view and shading are as in figure 1, but only the inner eighth of the domain is shown. Note, the QG solution at this time differs substantially – see figure 4 of McKiver & Dritschel (2008).

A complementary view is afforded by the spectral density of horizontal kinetic energy, $\mathscr{E}(k_{h},k_{z})=\overline{\hat{u} ^{2}+\hat{v}^{2}}$ , where the overline denotes a circular average over the wavevector phase ${\it\theta}$ , defined through $k_{x}=k_{h}\cos {\it\theta}$ and $k_{y}=k_{h}\sin {\it\theta}$ . Here, $k_{h}$ is the magnitude of the horizontal wavevector. Figure 4 plots $\mathscr{E}(k_{h},k_{z})$ , time averaged over the last half of each simulation, for the same four cases as illustrated in figure 3. In these plots, $k_{h}$ is scaled by $f/N$ , the natural QG scaling. Notably, there are scarcely any differences at all in the results. Even fine details of the fluid motion are insensitive to the value of $f/N$ . These results also indicate that the large-to-intermediate scales are approximately isotropic after the $f/N$ scaling, as previously found in QG turbulence (Dritschel et al. 1999; Reinaud et al. 2003). Anisotropy is evident only at small scales (high $k_{h}$ and $k_{z}$ ), where the increased power to the left of the $k_{h}$ – $k_{z}$ diagonal indicates that small-scale structures are relatively flat.

Figure 4. Comparison of the time-averaged spectral density of the horizontal kinetic energy $\mathscr{E}(k_{h},k_{z})$ for the simulations with ${\it\varepsilon}=0.25$ and for various values of $f/N$ ((a) 0.1, (b) 0.5, (c) 1.0, (d) 1.5). The horizontal wavenumber $k_{h}$ is scaled by $f/N$ . The time average is taken between $t=10$ and $t=20$ QG time units. Note, there is a 10-fold difference in $\mathscr{E}$ between adjacent contours, with $\mathscr{E}=1$ on the first dashed contour and $\mathscr{E}=10$ on the adjacent solid contour.

3.2. Vertical velocity

In large-scale geophysical flows it has been widely observed that the vertical velocity component $w$ tends to be much weaker than the horizontal velocity components, as the stable density stratification inhibits vertical motion. The weak vertical motion that remains, while ageostrophic, may not contain significant IGW activity. In fact, there can be a dominant balanced component $w_{\mathit{b}}\gg w_{\mathit{i}}$ , especially in carefully initialised flows at small $\mathit{Ro}$ and $\mathit{Fr}$ (see Dritschel & Viúdez 2003, 2007; Viúdez & Dritschel 2003; McKiver & Dritschel 2008).

Figure 5. Comparison of the full (a–d), balanced (e–h) and imbalanced (i–l) components of the vertical velocity field (in a $y=0$ cross section) at 20 QG time units ((a,e,i)  $f/N=0.1$ ; (b,f,j)  $f/N=0.5$ ; (c,g,k)  $f/N=1.0$ ; (d,h,l)  $f/N=1.5$ ). The plotted contours have values $\pm {\it\Delta}/2$ , $\pm 3{\it\Delta}/2,\dots \,$ , where ${\it\Delta}$ is the contour interval (negative contours have dashed lines, positive contours have solid lines; the zero value is omitted). For the full and balanced cases the contour interval is ${\it\Delta}=0.0002f/N$ . The imbalanced contour interval is one fifth of the balanced contour intervals.

Since in QG theory, $w_{\mathit{b}}$ is proportional to $f/N$ (see § 3.4 below), vertical motion increases with Prandtl’s ratio and might also enhance IGW activity. To explore this possibility, we use NQG balance to diagnose $w_{\mathit{b}}$ and $w_{\mathit{i}}$ in the same four cases as illustrated in figures 3 and 4. The results, at the final time $t=20$ , are given in figure 5, which compares $w$ , $w_{\mathit{b}}$ and $w_{\mathit{i}}$ (a–d, e–h and i–l) for increasing $f/N$ (left to right). Note, the contour intervals are chosen to be proportional to $f/N$ to facilitate comparison. Here, we show only a $y=0$ (vertical) cross-section, but this is typical of other cross-sections. The contour interval for $w_{i}$ is here five times smaller than that used for $w$ and $w_{b}$ , demonstrating just how well balanced the dynamics is over the entire range of $f/N$ values considered. In fact, $w_{i}/(f/N)$ decreases with $f/N$ and becomes more localised near strong PV anomalies. On the other hand, the structure of the full and balanced fields varies little with $f/N$ . This is because $w_{b}$ is controlled entirely by the PV, and the PV does not vary significantly with $f/N$ . The full field, $w$ , varies only slightly due to the differences in $w_{i}$ .

3.3. Imbalance

We next quantify the level of imbalance occurring in the same four simulations as illustrated above, principally to understand the dependence on $f/N$ . To this end, we consider the percentage of imbalance in various fields, defined by $\%{\it\xi}_{\mathit{i}}\equiv 100\langle {\it\xi}_{\mathit{i}}\rangle /\langle {\it\xi}\rangle$ for any field ${\it\xi}(\boldsymbol{x},t)$ . Figure 6 shows the time evolution of $\%u_{\mathit{i}}$ , $\%w_{\mathit{i}}$ and $\%\mathscr{D}_{\mathit{i}}$ together with the energy estimate $\%E_{\mathit{i}}$ (see (2.7)), for four widely different Prandtl numbers, all for ${\it\varepsilon}=0.25$ as before. Note, $\%v_{\mathit{i}}$ is virtually identical to $\%u_{\mathit{i}}$ and is therefore not shown. Each field exhibits a different behaviour. For $u$ and $w$ , the percentage of imbalance generally decreases with $f/N$ , although the levels of imbalance differ by almost a factor of 100: imbalance is exceptionally weak in the horizontal velocity field. Even in $w$ , the imbalance is only around 10 % or less, demonstrating that only a small portion of this field contains IGWs (and probably smaller than indicated, as NQG balance is not perfect – it is an overestimate of the true imbalance). For $\mathscr{D}$ , by contrast, the percentage of imbalance generally increases with $f/N$ , although it is again very weak. Finally, the percentage of imbalance in the energy norm, $\%E_{\mathit{i}}$ , first decreases slightly with $f/N$ then increases, although it remains ${\sim}0.1\,\%$ . The conclusion is that, for this PV-based Rossby number ${\it\varepsilon}=0.25$ , balanced motions due entirely to PV account for nearly the entire dynamical evolution, across a very wide range of Prandtl numbers $f/N$ . In § 3.5 below, we show that this remains true even for ${\it\varepsilon}=\mathit{O}(1)$ .

Figure 6. The time evolution of the percentage of imbalance in $u$ , $w$ , $\mathscr{D}$ and energy (a–d) for four Prandtl ratios $f/N=0.1$ (short dashed line), 0.5 (bold line), 1.0 (long dashed line) and 1.5 (thin line). Time in QG time units is shown along the horizontal axis. Here, as in figures 3–5, the PV-based Rossby number ${\it\varepsilon}=0.25$ .

The low level of imbalance seen here, we argue, is not principally due to the inevitable numerical diffusion required for numerical stability. Undoubtably, such diffusion removes some of the IGWs, yet the diffusion is also highly scale selective (and weak), removing mainly the smallest-scale features. Energetically, these features account for a very small proportion of the IGWs present, since IGWs are excited by the vortices and their interactions at scales comparable to the individual vortices, as seen for example in figure 5 (and in many previous works, e.g. Dritschel & Viúdez 2003, 2007; Viúdez & Dritschel 2003; McKiver & Dritschel 2008; Viúdez 2008; Tsang & Dritschel 2015). In a different norm (e.g. enstrophy) magnifying the importance of small-scale features, our conclusions may well be different. Yet in terms of energy, we find that flows close to a state of balance remain there, within limits imposed by the parameters ${\it\varepsilon}$ and $f/N$ .

3.4. Validity of quasi-geostrophic scaling

The results above indicate that the fluid motion remains dominantly balanced for values of $f/N$ up to unity or greater, at least for small to moderate PV-based Rossby numbers ${\it\varepsilon}$ . Here, we examine how well QG scale analysis applies to the fully NH dynamics. In QG theory, flow variables are assumed to take the following characteristic values:

(3.1a−c ) $$\begin{eqnarray}x,y\sim L,\quad z\sim H,\quad t\sim L/U,\end{eqnarray}$$

(3.2a−d ) $$\begin{eqnarray}u,v\sim U,\quad w\sim W,\quad b\sim B,\quad {\it\Phi}/{\it\rho}_{0}\sim P\end{eqnarray}$$

(cf. (2.2)). Assuming leading-order geostrophic and hydrostatic balance, we have that $U=P/fL$ and $B=P/H$ . However, geostrophic balance requires ${\it\varepsilon}\sim U/fL\ll 1$ , implying therefore $U={\it\varepsilon}fL$ and $B={\it\varepsilon}(fL)^{2}/H$ . The scale of the vertical velocity $W$ is found from the buoyancy equation $\partial b/\partial t+\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}b+N^{2}w=0$ , leading to

(3.3) $$\begin{eqnarray}W=\frac{BU}{N^{2}L}=\frac{{\it\varepsilon}^{2}f^{3}L^{2}}{N^{2}H}.\end{eqnarray}$$

In particular, the characteristic ratio of the vertical velocity to the horizontal velocity is given by

(3.4) $$\begin{eqnarray}\frac{W}{U}={\it\varepsilon}\frac{f}{N}\left(\frac{fL}{NH}\right).\end{eqnarray}$$

Thus, $W/U$ depends on the PV-based Rossby number, ${\it\varepsilon}$ , Prandtl’s ratio, $f/N$ , and the quantity $fL/NH$ . The latter is the inverse of the vertical-to-horizontal scale ratio $H/L$ divided by $f/N$ , namely $NH/fL$ , also known as (the square root of) the Burger number. This ratio is also equal to the Rossby number $U/fL$ divided by the Froude number, $U/NH$ . The validity of QG theory requires $(U/NH)^{2}\ll U/fL$ , which is satisfied when $U/NH\sim U/fL\ll 1$ , i.e. when $NH/fL\sim 1$ . In fact, $NH/fL\gg {\it\varepsilon}\text{}^{1/2}$ is sufficient.

We next estimate this ratio in the NH simulations using the domain-averaged horizontal kinetic energy $\mathit{K}=\langle |\boldsymbol{u}_{h}|^{2}\rangle /2$ and potential energy $\mathit{P}=\langle b^{2}/N^{2}\rangle /2$ . The QG scaling above implies $\mathit{K}\sim U^{2}=({\it\varepsilon}fL)^{2}$ and $\mathit{P}\sim B^{2}/N^{2}=({\it\varepsilon}(fL)^{2})^{2}/(NH)^{2}$ , so that we can estimate $NH/fL$ from

(3.5) $$\begin{eqnarray}\frac{NH}{fL}=\left(\frac{\mathit{K}}{2\mathit{P}}\right)^{1/2}.\end{eqnarray}$$

A factor of 2 is included here since, in QG theory, $|\boldsymbol{u}_{h}|^{2}/f^{2}=(\partial {\it\phi}/\partial x)^{2}+(\partial {\it\phi}/\partial y)^{2}$ while $b/N^{2}=\partial {\it\phi}/\partial \tilde{z}$ , where $\tilde{z}=Nz/f$ is the natural ‘stretched’ vertical coordinate. (For a spherical QG vortex in the coordinates $x$ , $y$ and $\tilde{z}$ , one can show that $\mathit{K}=2\mathit{P}$ , so that $NH/fL=1$ in (3.5).)

Figure 7(a) plots the time evolution of $NH/fL$ for the same four values of $f/N$ as illustrated in the previous figures, again for ${\it\varepsilon}=0.25$ , while figure 7(b) plots $w_{\mathit{rms}}/u_{\mathit{rms}}$ scaled by the QG estimate $W/U$ in (3.4). There is here very little difference across the wide range of Prandtl ratios $f/N$ considered, indicating that QG scaling works well even for $f/N=\mathit{O}(1)$ . Additionally, $NH/fL$ remains $\mathit{O}(1)$ , in fact a little less than 1, as found previously both in typical QG vortex interactions and in QG turbulence (Reinaud & Dritschel 2002; Reinaud et al. 2003). This is due to the fact that vertical shear is more destabilising than horizontal shear, and so a vortex must compensate by adopting an oblate shape (in the vertically stretched coordinates; see Reinaud & Dritschel (2002)). A wide selection of much longer simulations, extending to 100 QG time units or more, indicate that this scaling persists. The same is found for the scaled ratio of $w_{\mathit{rms}}/u_{\mathit{rms}}$ . This is consistently found to be around 0.04, even for much larger PV-based Rossby numbers (see below) and for much longer simulations. While QG scaling correctly predicts that $w_{\mathit{rms}}/u_{\mathit{rms}}$ is proportional to $W/U$ , it greatly overestimates the magnitude of this quantity. The vertical velocity is much smaller than a simple estimate would indicate, perhaps explaining why balance dominates the flow evolution. In the next subsection, we return to this issue in a wider context.

Figure 7. Evolution of (a) $NH/fL$ as estimated from the ratio of the kinetic to potential energy and (b) the ratio of the root-mean-square (r.m.s.) vertical velocity to the r.m.s. horizontal velocity scaled by the QG estimate for $W/U$ given in (3.4). Results are shown for $f/N=0.1$ (short dashed line), 0.5 (bold line), 1 (long dashed line) and 1.5 (thin line). All cases have ${\it\varepsilon}=0.25$ .

3.5. Limits of balance

We next explore the ${\it\varepsilon}$ – $f/N$ parameter space more widely to determine the limits of balance. More than 120 simulations were carried out, most extending to 100 QG time units or beyond, well into the decaying stage of the turbulent evolution. For a range of PV-based Rossby numbers ${\it\varepsilon}$ between 0 and 1, we progressively increased $f/N$ until the numerical code failed due to static instability, $S=-\partial \mathscr{D}/\partial z<-1$ (the code requires monotonically decreasing density with $z$ ). Such instability is often preceded by a period when the vorticity-based Rossby number $\mathit{Ro}={\it\zeta}/f$ is less than $-1$ somewhere in the flow (a necessary condition for inertial instability, cf. Knox (1997), Lazar, Stegner & Heifetz (2013)) or by a period when the Richardson number $\mathit{Ri}=N^{2}(1+S)/|\partial \boldsymbol{u}_{h}/\partial z|^{2}$ is less than $1/4$ (a necessary condition for shear or Kelvin–Helmholtz instability, cf. Howard (1961), Miles (1961), Hazel (1972)).

These conditions are generally not sufficient, and in particular we observe cases where $\mathit{Ro}$ is substantially below $-1$ over extended periods with no development of static instability. This is demonstrated in figure 8(a,b) for the case ${\it\varepsilon}=0.8$ and $f/N=0.8$ (the largest value of $f/N$ for this value of ${\it\varepsilon}$ ). The Rossby number $\mathit{Ro}$ drops to as low as $-1.717$ around $t=3.264$ QG time units, but $S_{\mathit{min}}>-0.85$ throughout the entire simulation. Moreover, the Froude number $\mathit{Fr}=\mathit{Ri}^{-1/2}$ remains less than 1.35, well less than the critical value of 2 necessary for shear instability (see panel c).

Figure 8. Time evolution (in QG time units) of the extreme values of (a) the vorticity-based Rossby number $\mathit{Ro}={\it\zeta}/f$ , (b) the static stability parameter $S=-\partial \mathscr{D}/\partial z$ , (c) the Froude number $\mathit{Fr}=|\partial \boldsymbol{u}_{h}/\partial z|/(N\sqrt{1+S})$ and (d) Rellich’s parameter $R={\it\Pi}-(f/N)^{2}|\boldsymbol{A}_{h}|^{2}/4$ (see text) for the case with ${\it\varepsilon}=0.8$ and $f/N=0.8$ .

Hence, despite having ${\it\zeta}/f<-1$ , the anticyclones in our simulations do not appear to be unstable. However, arguably, inertial instability is not relevant to stratified three-dimensional vortices: one cannot ignore the stabilising effects of stratification. Instead, instability requires that the (total) PV be negative: ${\it\Pi}<0$ (Sawyer 1947; Ooyama 1966; Charney 1973) (this is called ‘symmetric instability’ – see Lazar et al. (2013) for a detailed discussion). Notably, this never occurs in our simulations since PV is conserved and ${\it\Pi}>0$ in all cases.

The Rellich parameter $R$ shown in figure 8(d) is associated with the elliptic–hyperbolic character of the PV inversion equation (a double Monge–Amp\`ere equation for the vertical component of ${\bf\varphi}$ , see Dritschel & Viúdez (2003), in particular appendix A.2). This equation is considered to be elliptic if $R>0$ and hyperbolic if $R<0$ . Rellich’s parameter $R={\it\Pi}-(f/N)^{2}|\boldsymbol{A}_{h}|^{2}/4$ varies throughout the domain, and both elliptic and hyperbolic subdomains can exist side-by-side, as in the case illustrated (and in another originally studied in Dritschel & Viúdez (2003)). The elliptic/hyperbolic character of the equation is determined by linearising the Monge–Amp\`ere equation about a presumed solution; the resulting linear equation is then classified in the usual way (Bakelman 1994). Note, $R>0$ is a sufficient condition for both inertial and static stability, but $R<0$ does not necessarily imply instability (Dritschel & Viúdez 2003).

Notably, in this extreme case close to the limits of balance, the percentage of imbalance $\%E_{\mathit{i}}$ as measured by the energy norm (2.7) remains around 2 % throughout the entire simulation. This is shown in figure 9, which also plots the total energy versus time. The total energy here is seen to decrease by approximately 16.2 %, which may be the result of a forward cascade of the imbalance spectral energy density, as suggested in Bartello (1995). Such energy would be dissipated in the numerical code by the hyperdiffusion acting on $\boldsymbol{A}_{h}$ at small scales. However, other cases at smaller ${\it\varepsilon}$ and $f/N$ exhibit a comparable loss in total energy; e.g. when ${\it\varepsilon}=0.25$ and $f/N=0.1$ , the loss of total energy is 14.6 % whereas the time average $\%E_{\mathit{i}}=0.16$ , which is 12 times smaller than that found for ${\it\varepsilon}=0.8$ and $f/N=0.8$ . The conclusion is that, in these highly turbulent simulations, the balanced vortical motions contribute dominantly to the total energy dissipation; the imbalance contributes also, but to a much weaker degree.

Figure 9. Time evolution (in QG time units) of the total energy $E$ (solid curve) and the imbalanced part $E_{i}$ (multiplied by 50, dashed curve), for the case with ${\it\varepsilon}=0.8$ and $f/N=0.8$ .

The observed loss in energy mainly comes from the intense filamentation of PV occurring as vortices merge and strongly deform one other, cf. figure 1(b). This results in a strong forward spectral enstrophy cascade, ultimately leading to significant enstrophy dissipation. This also leads to energy dissipation, since the PV filaments not only contain enstrophy but also some energy, even if a much smaller proportion of the total.

A summary of the time-averaged minimum and maximum Rossby numbers ${\it\zeta}/f$ for all values of $f/N$ is given in figure 10(a). There is essentially no variation with $f/N$ . The same has been found for the other diagnostics shown in figure 8. Hence, the main dependence is on the PV-based Rossby number, and this dependence is shown in figure 10(b) (averaged over $f/N$ ) for the minimum and maximum Rossby numbers and the Froude number. For small ${\it\varepsilon}$ , $\mathit{Ro}_{\mathit{min}}\approx -{\it\varepsilon}$ (time averaged) while $\mathit{Ro}_{max}\approx {\it\varepsilon}$ (as indicated by the dashed lines). For larger ${\it\varepsilon}$ , the dependence becomes nonlinear, with anticyclones strengthening relative to cyclones. The Froude number shows a nonlinear dependence on ${\it\varepsilon}$ throughout, and rises steeply as ${\it\varepsilon}\rightarrow 1$ .

Figure 10. (a) Time-mean values of the minimum and maximum Rossby numbers, $\mathit{Ro}_{\mathit{min}}$ and $\mathit{Ro}_{max}$ , as a function of $f/N$ for all simulations conducted, namely ${\it\varepsilon}=0.125$ ( $+$ ), 0.25 ( $\times$ ), 0.375 (▵), 0.5 (▿), 0.625 (▫), 0.75 (◃), 0.8 (▹) and 0.9 (*). (b) Average across $f/N$ of the time-mean Rossby numbers (thin lines), with reference lines $\mathit{Ro}=\pm {\it\varepsilon}$ shown as dashed lines, together with the average time-mean Froude number (bold line), all as a function of the PV-based Rossby number ${\it\varepsilon}$ .

Cross-sections of various fields are shown in figure 11. The $y{-}z$ cross-section chosen passes through a strong anticyclone (in the lower left of each panel) having $\mathit{Ro}_{\mathit{min}}=-1.1427$ (located just above the diagonal $y=z$ , above and a little to the right of a similar-sized but slightly flatter vortex). The striking feature here is the conical shape of the anticyclones – this is in fact seen in all cross-sections and appears to be a generic feature. Cyclones are inverted but generally less well defined. The fact that all strong anticyclones adopt this shape suggests that it is the most stable shape for this Rossby number. Note, $S<0$ in anticyclones, implying less stability through decreased stratification ( $N\sqrt{1+S}$ ), while the opposite is true for cyclones (see Tsang & Dritschel 2015 for details). Rellich’s parameter $R$ most clearly displays the conical vortex shapes; the PV contribution dominates $R$ , so in this image we are mainly seeing the PV structure of the vortices. Finally, the vertical velocity $w$ is weak everywhere apart from in the vicinity of the strongly interacting cyclone and anticyclone in the upper right part of the domain ( $u$ and $v$ are $\mathit{O}(1)$ ). The anticyclone is in the process of splitting the cyclone vertically in two.

Figure 11. Cross-sections at $x=-43{\rm\pi}/64$ of the Rossby number $\mathit{Ro}$ (a), the static stability parameter $S=N^{-2}\partial b/\partial z$ (b), Rellich’s parameter $R$ (c) and the vertical velocity $w$  (d) for the flow with ${\it\varepsilon}=0.8$ and $f/N=0.8$ at $t=17$ QG units. The contour intervals are 0.1, 0.1, 0.1 and 0.005 respectively. The contour levelling is identical to that in figure 5.

We next consider general properties of the simulations that remained statically stable over the full time evolution. The highest PV-based Rossby number attainable was ${\it\varepsilon}=0.9$ , and then only for $f/N\leqslant 0.25$ . The domain of stability in the ${\it\varepsilon}$ – $f/N$ parameter space is shown in figure 12 alongside the time-averaged percentage of imbalance $\%E_{\mathit{i}}$ as measured by the energy norm (see (2.7) and the text below it). The percentage of imbalance is remarkably small, ${<}3\,\%$ , throughout the parameter space. It increases with ${\it\varepsilon}$ , approximately as ${\sim}{\it\varepsilon}\text{}^{3/2}$ when ${\it\varepsilon}\ll 1$ , and rises sharply near the highest PV-based Rossby number. Notably, this differs from the exponentially small scaling ( ${\sim}{\it\varepsilon}\text{}^{-1/2}\exp (-{\it\alpha}/{\it\varepsilon})$ ) found in Vanneste & Yavneh (2004) for weak disturbances to an idealised uniform horizontal shear flow. The differences may arise from the fact that in turbulence there is no scale separation between the wave scale and the underlying balanced vortical flow. We also find that as a function of $f/N$ , $\%E_{\mathit{i}}$ exhibits a shallow minimum near $f/N=0.5$ and tends to rise steeply near the maximum $f/N$ attainable for a given ${\it\varepsilon}$ . This suggests that the spontaneous emission of IGWs is weakest, albeit only marginally, when $f\approx N/2$ . Notably, there is no abrupt transition across the line $f/N=1$ .

Figure 12. (a) Summary of parameter space, with statically stable flows marked by a square and unstable flows marked by a cross (the bold curve is an estimate of the stability boundary, see text); (b) the time-averaged percentage of imbalance $\%E_{\mathit{i}}$ , as measured by the energy norm, plotted as a function of $f/N$ for ${\it\varepsilon}=0.125$ ( $+$ ), 0.25 ( $\times$ ), 0.375 (▵), 0.5 (▿), 0.625 (▫), 0.75 (◃), 0.8 (▹) and 0.9 (*). Only the statically stable flows are considered in (b).

As shown in figure 12(a), the maximum attainable $f/N$ is well approximated by an elliptic curve found by a least-squares fit of $(f/N)_{max}^{2}$ to a linear function of ${\it\varepsilon}^{2}$ . This gives

(3.6) $$\begin{eqnarray}\frac{{\it\varepsilon}^{2}}{0.920^{2}}+\frac{(f/N)^{2}}{1.567^{2}}=1.\end{eqnarray}$$

This curve defines the approximate limits of balance for the fully developed turbulent flow investigated.

Two other diagnostics are shown in figure 13, namely the late-time-averaged values of $w_{max}/u_{max}$ , scaled by the factor $W/U$ given in (3.4), and the ratio $NH/fL=(\mathit{K}/2\mathit{P})\text{}^{1/2}$ , computed from the domain-averaged kinetic and potential energies $\mathit{K}$ and $\mathit{P}$ . Throughout parameter space, the vertical velocity is much weaker than the horizontal velocity, and the ratio $w_{max}/u_{max}$ is consistently overestimated (by a factor of 20–40) by the QG estimate $W/U$ . This, however, may be the main reason why balance retains such a tight control on the dynamics. The QG estimate $W/U$ in (3.4) simply gives the scaling $w_{max}/u_{max}$ with $f/N$ and ${\it\varepsilon}$ . The prefactor in the scaling varies by only 25 % throughout parameter space. Furthermore, the ratio $NH/fL\approx 0.86$ with less than a 3 % variation. This is fully consistent with QG scaling, in which $H/L\sim f/N$ (Charney 1971; Herring 1980; Dritschel et al. 1999; Reinaud & Dritschel 2002; Reinaud et al. 2003), even though the flows considered here are strongly ageostrophic when ${\it\varepsilon}=\mathit{O}(1)$ .

Figure 13. Late-time averages (over $15\leqslant t\leqslant 20$ QG time units) of (a) $w_{max}/u_{max}$ scaled by the QG estimate $W/U$ and (b) $NH/fL$ , versus $f/N$ for ${\it\varepsilon}=0.125$ ( $+$ ), 0.25 ( $\times$ ), 0.375 (▵), 0.5 (▿), 0.625 (▫), 0.75 (◃), 0.8 (▹) and 0.9 (*).

4. Discussion and conclusions

In this work we have shown that the dynamics of a rotating stably stratified freely evolving turbulent flow has little dependence on the Prandtl ratio $f/N$ . This statement holds throughout a broad parameter space spanned by the PV-based Rossby number ${\it\varepsilon}$ and $f/N$ , which may take $\mathit{O}(1)$ values. Outside of this parameter space, roughly beyond the elliptical curve $({\it\varepsilon}/0.9)^{2}+((f/N)/1.6)^{2}=1$ , all flows are found to be statically unstable, leading to overturning density surfaces and a breakdown of the numerical method employed. While this conclusion is based on simulations of uniform PV patches, arguably such patches are less regular than continuous PV vortices and are consequently a more demanding test of balance.

Where flows remain statically stable, they are also found to be well balanced, in the sense that one can deduce nearly all information on the flow field from knowledge of the PV field alone. The residual imbalance, characterised by IGWs, plays a very minor role, especially in the horizontal velocity field and in the displacement of density surfaces. Even in the vertical velocity field, balance dominates across the parameter space. As measured by an energy norm, imbalance never contributes more than 3 % to the entire flow evolution, and often much less. The key point is that, for carefully initialised flows where IGWs are weak, they remain weak. This is not the result of numerical damping, but rather a physical property of the governing equations when the PV-based Rossby number is less than unity.

Higher PV-based Rossby numbers lead to static instability, as do large values of $f/N$ . Under such conditions, balance cannot be expected to hold well everywhere. Local breakdown is inevitable, leading to strong emission of IGWs. Yet, even under these circumstances, the breakdown may be highly localised and short-lived before stable stratification is re-established, resulting in only a small contribution to the overall imbalance. This regime, however, is beyond the scope of the present analysis – an entirely different numerical approach is required.

As the effects of rotation and stratification weaken further, one expects to find highly active three-dimensional turbulence characterised by density overturning and mixing, strong vortex stretching, direct energy cascades to small scales and significant energy dissipation. Such behaviour has been frequently reported in past studies of geophysical turbulence (Molemaker, McWilliams & Yavneh 2005; Waite & Bartello 2006; Molemaker et al. 2010; Deusebio et al. 2013). It is claimed to occur even for small Rossby numbers, i.e. close to a state of QG balance. However, crucially, those studies have not examined the PV-based Rossby number, the only Rossby number that is invariant in a freely decaying flow. These past studies have used large-scale forcing (with little attention to balance, or employing only leading-order geostrophic balance) to sustain the turbulence, and so PV is not conserved. Understandably, a PV-based Rossby number in these circumstances appears to have no advantage. However, Rossby numbers based on measured velocity and length scales, we argue, are much less relevant to understanding the conditions under which geophysical turbulence remains close to a state of balance – certainly for freely decaying flows and likely also for forced flows. We have found that the key parameter controlling balance is the PV-based Rossby number: this must remain less than unity and Prandtl’s ratio $f/N$ must not be much larger than unity.

For imbalanced turbulent flows, the strong forward cascade of energy can only be maintained by a significant source of energy, without which the turbulence will ultimately decay, restratify and approach a balanced state by dispersing IGWs (Bartello 1995; Praud et al. 2006). Thus, for freely evolving rotating stratified turbulence, a natural attracting state appears to be balance.