Balance in non-hydrostatic rotating shallow-water ﬂows

Unsteady nonlinear shallow-water ﬂows typically emit inertia–gravity waves through a process called ‘spontaneous adjustment-emission’. This process has been studied extensively within the rotating shallow-water model, the simplest geophysical model having the required capability. Here we consider what happens when the hydrostatic assumption underpinning the shallow-water model is dropped. This assumption is in fact not necessary for the derivation of a two-dimensional or single-layer ﬂow model. All one needs is that the horizontal ﬂow ﬁeld be independent of height in the ﬂuid layer. Then, vertical averaging yields a single-layer ﬂow model, with the full range of expected conservation laws, similar to the shallow-water model yet allowing for non-hydrostatic effects. These effects become important for horizontal scales comparable to or less than the depth of the ﬂuid layer. In a rotating ﬂow, such scales may be activated if the Rossby deformation length (the ratio of the characteristic gravity-wave speed to the Coriolis frequency) is comparable to the the depth of the ﬂuid layer. Then, the range of frequencies supporting inertia–gravity waves is compressed, and the group velocity of these waves is reduced. We ﬁnd that this change in wave properties has the effect of strongly suppressing spontaneous adjustment-emission and trapping inertia–gravity waves near regions of relatively strong circulation.


I. INTRODUCTION
In the atmosphere and oceans, dynamical and thermodynamical fields often closely satisfy certain relations, called balance relations, the simplest of which are hydrostatic and geostrophic.The latter apply when the fluid acceleration is small compared to the remaining terms in the momentum equations.Such balance relations (which can be more complicated [1][2][3][4][5][6][7][8] ) are useful for separating relatively high frequency waves, such as inertia-gravity waves, from the remaining relatively low frequency balanced flow 3,7,[9][10][11][12][13][14] .This balanced flow can be thought of as arising from the potential vorticity (PV) field (the master variable), and the balanced relations provide a means to determine all of the (balanced) dynamical and thermodynamical fields by 'PV inversion' 15 .
Balance has been an important theoretical concept, underlying the derivation of reduced sets of equations like quasigeostrophic [16][17][18] which have proved immensely fruitful.This concept has also enabled researchers to better understand the various ways nearly-balanced flows spontaneously emit inertia-gravity-waves [19][20][21][22][23][24][25][26][27][28][29][30][31][32][33][34] , or how imbalanced flows adjust to a nearly balanced state [35][36][37][38][39][40] .Balance has also been important practically, for example in estimating oceanic flow fields from limited data, [41][42][43][44] or in weather forecasting during data assimilation. 45ere we will not attempt to review the vast literature on the subject of balance, as excellent reviews can be found in the two collections 22,23 , and more recently in a special issue on spontaneous imbalance 46 .Here we focus on single-layer rotating shallow-water (SW) flows, possibly the most widelystudied model in this context due to its relative simplicity: there are just three scalar evolution equations, one of which can be taken to express material conservation of PV, while the other two allow for inertia-gravity waves (IGWs), the imbalance 3 .Replacing the latter by balance relations results in a balanced model with no IGWs.The novelty in this paper is to relax the hydrostatic approximation which forms the basis of the traditional SW model.We shall see that this can have a profound impact on IGWs, generally weakening their emission.
The SW model has a long history going back to Saint-Venant 47 in 1871.The model reduces the parent three-dimensional Euler equations to a single-layer twodimensional set of equations by assuming the horizontal flow is independent of depth and by imposing the hydrostatic approximation.These are valid so long as horizontal scales L are large compared to the mean fluid depth H.In a rotating flow with Coriolis frequency f , there is an additional length scale L D = c/ f called the 'Rossby deformation length' where c = √ gH is the short-scale gravity wave speed and g is the acceleration due to gravity (or reduced gravity). 18The tacit assumption is that L D H, or at least that L H even when L D ∼ H.However, commonly, rotating SW flows develop small scales, especially in PV, as a result of nonlinear flow interactions 48,49 .Horizontal scales with L < L D inevitably form, and so the validity of the rotating SW model requires H L D , i.e. a mean fluid depth much smaller than the intrinsic length scale imposed by rotation.This is trivially satisfied by a non-rotating flow since then L D → ∞.
Regardless of rotation, the hydrostatic approximation underpinning the SW model breaks down when L ∼ H. Nevertheless, one can still construct a single-layer two-dimensional flow model that conserves PV and all the integral invariants associated with symmetries.This model was first derived by Serre 50 in 1953, then Su & Gardner 51 in 1969, but is often credited to Green & Naghdi 52 in their 1976 paper.4][55][56][57][58][59][60][61][62][63][64][65][66][67] .Given the many people who could be credited for this model, it might be fairer to call it the 'non-hydrostatic shallow-water model', especially since the only difference from the SW model is the relaxation of the hydrostatic approximation.Notably, this model and the original SW model both still assume that the horizontal flow is independent of depth.Moreover, one can show that the non-hydrostatic SW model is simply the vertical average of the three-dimensional Euler equations. 49ll conserved quantities, including PV, are just the vertical averages of their three-dimensional counterparts.For this reason, we call this model the 'Vertically-Averaged (VA) model' below.
The inclusion of rotation, or the Coriolis acceleration (relevant for applications to geophysical fluid dynamics), is relatively recent 55,[68][69][70] .Pearce and Esler (2010) 69 derived the VA equations in their vorticity-divergence form, widely used in SW studies of atmosphere/ocean dynamics and convenient for a pseudo-spectral numerical treatment.They verified their numerical model for an analytic solution of a propagating uni-directional Cnoidal wave, and studied the PV evolution in an unstable jet (we have also done similar tests for the method used here 49 ).Recently, Alemi Ardakani (2021) 70 derived generalised variational SW and VA models for a shallow fluid sloshing inside a container subject to arbitrary (threedimensional) translations and rotations.They proved that there exists a materially-conserved PV that is a combination of the PV found by Miles and Salmon (1985)  53 for the non-rotating VA equations, and the PV found by Dellar and  Salmon (2005)  55 for the SW equations under arbitrary rotation (complete Coriolis force).
Most of these works have been theoretical.Prior to our own work 49,71,72 , only Pearce and Esler 69 studied the evolution of nonlinear rotating flows (and they restricted attention to the PV dynamics in a single example).Little, therefore, is known about the actual differences between the SW and VA models, or about the general properties of rotating VA flows.Significantly, we have shown that the VA model is substantially more accurate than the SW model by comparing simulations of horizontal shear instability directly with simulations of the three-dimensional Euler equations with a free surface. 49The greater accuracy of the VA model is well known in studies of unidirectional flows without rotation, and this is due in part to the better representation of wave dispersion, a feature entirely absent in the SW model (without rotation).
In the present paper, we focus on the IGWs generated from initially balanced, turbulent, rotating shallow-water flows, comparing and contrasting the SW and VA models.These flows develop small scale features, not only in PV, but also in vorticity and (horizontal) divergence.We find that the key parameter is the frequency ratio f /N, where f is the Coriolis frequency and N = 3g/H may be called the 'buoyancy frequency' since, like in the parent three-dimensional model, linear wave frequencies lie between f and N. When f /N 1, the SW and VA models agree closely (and agree perfectly in the limit f /N → 0).However, when f /N ∼ 1, there is a stark difference in the form and amplitude of the IGWs emitted in the SW and VA model simulations, to the extent that IGWs in the VA model are almost entirely suppressed.Notably f /N ∼ 1 implies L D ∼ H, i.e. that the Rossby deformation length is comparable to the mean fluid depth.This parameter regime is not typical in the Earth's atmosphere and oceans, where f /N ∼ 10 −2 -10 −1 .But this is an average estimate, [16][17][18] and there are significant regional vari-ations.Regions of weak stratification do occur, such as in the Mediterranean Sea 73,74 and in polar oceans 75 .Furthermore, even when f /N 1, strong dispersion occurs for horizontal scales comparable to or smaller than the fluid depth, an effect not captured by the SW model.More generally, studying the impact of this dispersion contributes to our fundamental understanding of geophysical fluid dynamics, both on Earth and on other planets.
The plan of the paper is as follows.In the next section, we briefly review the SW and VA models and the numerical method used to perform simulations.In section 3, we present our results, comparing a series of simulations differing only in the mean depth H, with the limit H → 0 corresponding to the SW model.Here we compare the evolution of various fields, balanced and imbalanced, and also examine wavenumber and frequency spectra.Conclusions are offered in section 4, where we discuss the significance of our results for geophysical rotating, weakly-stratified flows.

II. THE FLOW MODELS AND THEIR NUMERICAL TREATMENT
In the momentum-mass formulation, both the SW and VA models can be written in the form where u = (u, v) is the two-dimensional vector velocity, h is the free surface height (here above a flat bottom at z = 0), f is the Coriolis frequency, u ⊥ = (−v, u), and p is the verticallyintegrated pressure (divided by the constant density).In the SW model, where g is the acceleration due to gravity (or reduced gravity).
Then the pressure gradient term in (1) reduces to the familiar −g∇h hydrostatic acceleration.In the VA model, there is an additional non-hydrostatic pressure p n so that where p n is determined from h and u by solving the linear elliptic equation 49 where in which ζ = −∇ • u ⊥ is the (vertical component of the) vorticity, δ = ∇ • u is the (horizontal) divergence, and J(•, •) is the Jacobian operator.
Both models possess a material invariant, potential vorticity (PV) q, a direct consequence of Kelvin's circulation theorem.That is, PV is conserved following fluid 'particles', i.e.
In the VA model, PV is the vertical integral of the 3D Rossby-Ertel PV 49 , first shown by Miles & Salmon (1985). 53.It is similar to the well-known SW form but has an additional Jacobian term: Besides PV, both models conserve energy and (in an infinite domain) linear and angular momentum. 49.They also conserve all integrals of functionals of PV, the 'Casimir invariants', as a consequence of material conservation of PV.
7][78][79][80] .This improves the representation of the balanced part of the flow.To also improve the representation of the imbalanced part, the other two variables should be chosen to represent the departure from balance.However, there is no exact definition of balance, so a practical choice is made that still allows for efficiency.A simple but fruitful choice is to use the horizontal divergence δ = ∇ • u and the quantity γ = f ζ − g∇ 2 h, which is the acceleration divergence in the SW model (and the linearised part of it in the VA model).This choice was also made in the paper on which the present work is based. 49Notably, this choice allows one to recover the original variables h and u by solving linear elliptic equations.Under geostrophic balance, both δ and γ vanish.Full details of the numerical method can be found in a previous paper. 49The only difference here is that we use contour advection 81 to accurately represent the PV evolution.In this method, the PV is represented by a set of contours, here 80 contour levels equally spaced in PV, which are advected by the velocity field interpolated from the underlying regular grid.Each contour is represented by a variable number of nodes connected together by local cubic splines, and 'surgery' 82 is used to limit complexity (surgery operates at a sixteenth of the grid spacing).Besides PV contours, two additional PV fields are used to improve the representation of large scales and energy conservation; these fields are blended with the PV associated with the contours to provide a highly accurate representation of PV at all scales. 81][85] We consider flow in a doubly-periodic domain of side length 2π.The basic grid resolution is 256 2 or 512 2 , but the effective resolution for the PV field is 16 times finer in each direction.We mainly report on the 512 2 simulations, but the 256 2 simulations produce qualitatively, and often quantitatively, similar results.The numerical algorithm settings are otherwise the standard recommended ones. 86

III. RESULTS
In this section, we first explain the set-up of the numerical experiments, specifically the initialisation and the physical parameters used.We then illustrate the flow evolution in one case, comparing various fields in the SW simulation (the limit H → 0) with those in the VA simulation.Next we turn to the imbalanced fields -the IGWs -comparing fields, r.m.s.(root mean square) norms, wavenumber and frequency spectra.Finally, we briefly discuss analogous results for different initial conditions.

A. Initialisation and balance
Eight sets of simulations were performed, composed of two Rossby deformation lengths, two Rossby numbers and two resolutions.In each set, H → 0 (the SW case), H = 0.1, 0.2 or 0.4 (recall that domain width is 2π).We consider two Rossby deformation wavenumbers, k D = L −1 D = 6 and k D = 12, and two Rossby numbers ε = 0.2 and 0.6 (defined below).Without loss of generality, we take the Coriolis frequency f = 4π so that a unit of time is a 'day'.In the SW model, note that only the squared short-scale gravity wave speed c 2 = gH appears in the governing equations, after scaling h by the mean depth H.The definition of L D = c/ f then provides c.In the VA model, H remains in the equations after this scaling (through p n ), but in the limit H → 0, the VA model reduces to the SW model.Since we are free to choose a characteristic time scale (here T = 4π/ f = 1) and the domain width (here L dom = 2π), the only independent parameters in the problem are k D and ε in the SW model, and additionally H in the VA model.
We specify the initial PV field through its power spectrum where k is the wavenumber, k 0 is the peak enstrophy wavenumber, and C is determined by requiring |q| max = ε f where ε is the specified PV-based Rossby number.Note, S q is the squared spectral amplitude | q(k)| 2 summed over all wavevectors k lying in the shells Within each shell q(k) is otherwise chosen randomly.We mainly report on simulations with k 0 = k D , motivated by the fact that oceanic eddies typically have scales comparable to L D and dominate the oceanic enstrophy spectrum, see Venaille and Vallis (2011) 87 and references therein.We have also performed simulations starting with k 0 k D , i.e. large-scale conditions, to verify that non-hydrostatic effects still strongly suppress high frequency IGWs (see below).While atmospheric eddies (cyclones) may also exhibit scales comparable to L D , 88 the variation of the Earth's planetary vorticity (Coriolis frequency) over this scale is much greater than in the oceans and cannot be neglected.
For each value of k D , we use the same random number seed in all simulations so that the initial PV field has the same form (for k 0 = k D ).This allows the closest possible comparison between flows.But the PV field is not the only field required for initialisation: we also need to prescribe the divergence δ and (linearised) acceleration divergence γ.Here, these fields are determined by the balance relations which, when explicitly written out, are nonlinear equations to determine h and u (and hence δ and γ) after additionally using the definition of PV in (8). 49Explicitly, the balance relations are which must be solved together with (8), ∇ • u ⊥ = −ζ and ∇ • u = δ , given only the PV field q.In practice, these equations are solved iteratively by a simple adaption of the full, time-dependent numerical algorithm.Solutions converge strongly, even at high Rossby numbers.The same balance relations are used diagnostically to find the balanced fields at any time in a simulation.Note: the mean value of PV is determined by the requirement that the domain integral of ζ must be zero (a consequence of Stokes' theorem). 3gure 1 shows the initial fields of q, δ and γ for the VA simulation with k D = 6, ε = 0.6 and H = 0.4.The non-zero values of δ and γ imply the flow is ageostrophic.Compared with the vorticity ζ (not shown), both δ and the ageostrophic vorticity γ/ f = ζ − g∇ 2 h/ f are small, respectively 0.12% and 8.0% in an r.m.s.measure.Recall, these are balanced fields, not IGWs.The vorticity-based Rossby number |ζ | max / f starts at 0.404 and peaks at 0.511, while the Froude number (|u|/ √ gh) max starts at 0.241 and peaks at 0.266.Another indicator of the flow ageostrophy is the displacement of the free surface, relative to the mean depth, measured by the dimensionless height anomaly Initially, the min/max values are −0.238/0.367,and these peak at −0.269/0.529.[91]

B. Flow evolution
The evolution of the PV field is shown at a few characteristic times in figure 2 for both the SW model (top row) and the VA model for H = 0.4 (bottom row).Initially the flow grows in complexity as like-signed PV regions, vortices, merge and weak filamentary debris stretch and mix in between.In time, the number of vortices decreases and the PV becomes more well mixed between the vortices.By t = 500, the final time of the simulations, one clearly sees that the anti-cyclones (in blue) are much more compact and circular than the cyclones (in red).There are no qualitative differences between the SW and VA evolution, and the similarity in the PV fields at t = 50 is striking.The differences by t = 500 are expected given the slightly different initial conditions and the different models used.Indeed one might have expected greater differences.
We next examine the divergence evolution.Figure 3 compares δ in the SW model with that in the VA model, now using a depth four times smaller, H = 0.1.At early times, the fields compare well, but by t = 500 the VA divergence field is significantly broader scale and nearly twice as large in amplitude.This difference occurs earlier for larger H (not shown).The key point is that even a small value of H can have a significant impact on the evolution of δ , significantly altering its spatial structure.Notably, comparing the balanced part of δ (not shown), we find much closer agreement, with both fields of larger scale.The small-scale features in the SW simulation are therefore mainly IGWs, examined in more detail below.

C. Diagnosis of imbalance
As explained above in section III A, the balanced fields are defined to be those which satisfy the balance relations (10).Henceforth, these are subscripted by b for clarity, e.g.δ b for the balanced divergence.The imbalance is just the difference from the full field at the time of diagnosis, e.g.δ i = δ − δ b for the imbalanced divergence.While there is no perfect balance in a general flow, for low to moderate Rossby and Froude numbers the balance defined this way provides a good estimate of the balanced flow, certainly much better than hydrostatic-geostrophic balance (p n = δ = γ = 0). 3,13gure 4 compares the imbalance divergence δ i at the final time in the SW simulation (H → 0) and in three VA simulations (H = 0.1, 0.2 and 0.4) for the same case illustrated above.Qualitatively similar results are found at earlier times (not shown).In the SW simulation, δ i is of moderately small scale and appears to be randomly distributed.There is no clear association between δ i and the PV field, whose changes induce (weak) IGW emission.A significant portion of the full divergence field δ (shown in the top right panel of figure 3) is imbalanced, about 67.9% in an r.m.s.measure.For H = 0.1, the structure of δ i is similar but at a larger scale and a smaller amplitude: now only 20.5% is imbalanced.For H = 0.2, there is a qualitative change.Now the imbalance appears to be trapped around intense vortices, predominantly anti-cyclones (see right panels in figure 2 for q at this time and for H → 0 and H = 0.4).The amplitude of this imbalance is even smaller, now only 3.00% of the full divergence.Similar results are found for H = 0.4, with 2.96% of the divergence being imbalanced.Clearly finite H has a major influence on IGWs.The trapping effect seen here can be partly explained by the dispersion relation for IGWs on a basic state at rest.The frequency ω of such waves satisfies 71 where k is the wavenumber as before.Notably, in the SW case (H → 0), there is no upper bound on |ω| and moreover short waves are non-dispersive: all have phase and group velocities equal to ±c.Hence, all IGWs propagate away from their source, which explains the random pattern seen in δ i in the left panel of figure 4. For finite H, by contrast, |ω| is bounded between the Coriolis frequency f (when k → 0) and the buoyancy frequency N ≡ √ 3c/H = 3g/H (when k → ∞).The latter name is given due to the similarity with IGWs in a threedimensional rotating statified Boussinesq flow, 16,18,38 which also have frequencies between f and N.In large parts of the atmosphere and oceans, f /N < 1 (even 1), so N is normally the maximum frequency.However, when f /N > 1 as in a weakly-stratified (or strongly-rotating) flow, the situation is reversed.
For H > 0, the phase velocity c p = |ω|/k is a monotonically decreasing function of k which vanishes as k → ∞.However, the group velocity reaches a finite maximum of c g,max when µ ≡ k 2 H 2 = ( 1 + 3ξ − 1)/ξ , where ξ = N 2 / f 2 (this is most easily shown by expressing c 2 g /c 2 as a function of µ).The result is that c g,max /c depends only on the frequency ratio f /N, and in particular c g,max /c vanishes when f /N = 1.When this occurs, linear disturbances of all wavelengths are trapped.
The dependence of c g,max /c on f /N is shown in figure 5, along with the wavenumber k = k max of maximum group velocity.For small f /N, k max ≈ 1/ √ L D H and c g,max ≈ c: this corresponds to the hydrostatic, SW limit.As f /N increases toward 1, there is a rapid decrease of c g,max /c, then for f /N > 1 an increase, ultimately tending to (2/ √ 27) f /N for large f /N (meanwhile k max H → 3/2).
For the VA simulations conducted, f /N = 0.3464, 0.6928 and 1.3856 approximately for H = 0.1, 0.2 and 0.4 respectively.The two values straddling unity occur for H = 0.2 and 0.4, where we see the most trapping in figure 4. Nonlinear effects, neglected in this analysis, do not appear to play a major role.However, an additional simulation conducted for f /N = 1 (for which IGWs are completely trapped in linear theory) does not exhibit any qualitative differences from the simulations having H = 0.2 and 0.4, suggesting that nonlinearity prevents complete trapping (further results are provided below).
We turn next to the (vertically-integrated) pressure field p, which in the VA model contains a non-hydrostatic part, p n , see ( 4).This field is perhaps better thought of as a potential energy density, but then only the hydrostatic part gh 2 /2 is relevant. 49,71The non-hydrostatic part p n is here assumed to be entirely imbalanced, see (10), but there is also some (often very weak) imbalance in the hydrostatic part gh 2 /2.For example, for a VA simulation with H = 0.4, figure 6 shows how the imbalance present in h largely cancels that of p n , especially at small scales.The result is that the imbalanced pressure p i is both substantially smaller in amplitude and larger in scale than p n .Note, these results are for the most nonlinear flow considered, yet the amount of imbalance is exceedingly small.The variation of p i with H is shown in figure 7. Compared to the imbalanced divergence δ i in figure 4, we see that p i is generally of larger scale but otherwise exhibits a similar variation with H.At small H, p i is broadly distributed indicating weak or non-existent IGW trapping.At larger H, when f /N ∼ 1, p i exhibits the same trapping only with fewer finescale features.To appreciate the degree of balance here, in an r.m.s.measure, p i is only 0.023% of the total (verticallyintegrated) pressure in the SW simulation (at this time), and this decreases to 0.013% for H = 0.1, then to 0.0015% for H = 0.2 and then increases slightly to 0.0046% for H = 0.4.Pressure, like height h, is very close to balance in these flows, despite the moderate values of the Rossby and Froude numbers (about 0.34 and 0.13 at this time).Flows with f /N ∼ 1 are particularly devoid of IGWs.
The time variation of the imbalance, measured by the r.m.s.divergence δ i and pressure p i , is shown in figure 8 (other quantities exhibit similar behaviour).Notably, in the SW simulation (H → 0), the imbalance decreases only slightly after a small initial increase.When H > 0, the decrease is more marked, especially for H = 0.2 and 0.4, which show a continued decrease until late in the simulation and drop by a factor of around 10. Except for the very earliest times, the imbalance when H > 0 is significantly weaker than in the SW simulation.In divergence δ i , both H = 0.2 and 0.4 exhibit a similar decrease in the r.m.s.norm.However, for p i , H = 0.2 clearly exhibits the least imbalance at all times.Notably, the non-hydrostatic pressure p n , which contributes to p i , actually grows with H (see (5) and recall p n = 0 in the SW limit H → 0).So while the flow becomes increasingly nonhydrostatic as H increases, the degree of imbalance appears to depend mainly on f /N, with a minimum occurring around f /N = 1 where linear IGWs are completely trapped.Similar results are found for the lower Rossby number ε = 0.2 examined (not shown), except that δ i and p i are about 10-20 times smaller.
The decay of IGWs in the VA simulations must be due the reduced range of frequencies f < |ω| < N over which IGWs can be excited by spontaneous adjust emission.As the flow evolves, the PV field in particular develops increasingly sharp gradients (absent at t = 0), and as these fronts move, they excite all frequencies.The high-frequency cut-off for IGWs in the VA model means that frequencies |ω| > N do not excite IGWs, unlike in the SW model (this is shown explicitly below).
We turn next to the scale distribution of the balance and imbalance by considering the spatial power spectrum of divergence, S δ (k).The power spectrum is defined as usual as the sum of squared spectral amplitudes in wavenumber shells k = constant in Fourier space (see discussion following ( 9)).Results for k D = 6 and for two Rossby numbers ε = 0.2 and 0.6 at t = 500 are shown in figure 9 (left and right panels).At the smaller Rossby number (left panel), the imbalanced divergence is much smaller than the balanced or total divergence (which are indistinguishable here), especially at large scales.The turn-up in imbalanced divergence at the highest wavenumbers is a hyperdiffusion effect: the full fields are evolved using hyperdiffusion while the balanced fields were obtained without hyperdiffusion.The difference results in a spurious imbalance in the range of wavenumbers affected by hyperdiffusion.
Comparing the SW case (H → 0) with the VA one for H = 0.1, the most significant feature is the strong reduction in power at small to intermediate scales, mainly at wavenumbers k > √ 3/H = 10 √ 3 (or log 10 k > 1.24 approximately).At such wavenumbers, strong dispersive effects take hold, as seen e.g. in the frequency dispersion relation (14) for linear IGWs.Similar results are found at higher Rossby number (right panel), except that the reduction in power occurs over a wider range in wavenumbers, and now the SW flow is significantly less balanced.In fact, the imbalance exceeds the bal- ance over an intermediate range of wavenumbers near k = 10.By contrast, the VA flow, even for this small value of H, remains well balanced across all scales.
The reduction in power is even greater at larger H, as shown in figure 10 which focuses on the imbalanced divergence only, but now includes an additional value of H for which f /N = 1 (i.e.H = 0.288675).This demonstrates the strong effect that wave trapping has on suppressing IGW emission, at both Rossby numbers considered.At the higher Rossby number, the linear theory behind this wave trapping mechanism is less accurate, but even so the IGW suppression is remarkably strong when f /N ∼ 1. Lower resolution simulations (not shown) indicate that the high wavenumber IGWs are partially the result of the numerical discretisation: they exhibit increased power at lower resolution in the range k > k max /4 approximately.Higher resolution simulations, therefore, are likely to exhibit even less power at high k than shown here.
On the other hand, resolution has a negligible effect on the balanced flow.What figure 10 shows is the tiny difference between the full and balanced divergence fields.A complementary view is afforded by considering the frequency power spectrum, measuring the amplitude of imbalance present at each frequency ω.Here this is done for the divergence field but the results are similar for other fields.At each time step, the divergence is recorded at a regular array of 16 grid points.This produces a time series of divergence (at each of these points) which can be Fourier analysed to create a power spectrum as a function of sidereal frequency T −1 = ω/(2π).The 16 spectra thus formed are finally averaged to produce the frequency power spectrum P δ .
The results are shown in figure 11 for Rossby numbers ε = 0.2 (left) and 0.6 (right), for various values of H as indicated.The vertical dashed lines mark the location of the Coriolis frequency f (cyan) and the buoyancy frequencies N (coloured the same way as the spectra).For the SW case (H → 0), there is no corresponding dashed line for N since N → ∞.Recall that (linear) IGWs have frequencies between f and N, and the results in figure 11 are consistent: there is a bulge in power in this range, likely due to IGWs.Note, the power spectrum also includes the balanced flow, as it is practically impossible to generate a time series of δ i (this would require balancing the flow at every time step).Thus, the higher power at lower frequencies mostly corresponds to the balanced flow.
The key finding here is that finite H not only closes the frequency gap between f and N, leading to wave trapping, but also significantly reduces wave amplitudes.Non-hydrostatic effects reduce imbalance, so long as f /N is not much larger than unity.

D. Smaller deformation length
We next consider the effect of halving the deformation length L D or doubling the deformation wavenumber k D from 6 to 12. Likewise, the initial scale of the flow is halved (we take k 0 = k D in ( 9)).Again, two Rossby numbers are considered, ε = 0.2 and 0.6, together with the same four values of H (occasionally adding a fifth corresponding to f /N = 1).The flow evolves in a similar way to that already illustated for k D = 6, albeit at a somewhat slower rate due to the increased value of k D . 89Results are not shown, but the pattern of the PV field is closely similar to that already shown in figure 2 except on half the scale.
The imbalanced fields are also similar, except that the effects of finite H come into play sooner since Figure 12 shows this for the imbalanced pressure p i .Note that f /N = 1 corresponds to H = 0.144338 approximately, a value which lies between H = 0.1 and H = 0.2 (the middle two panels in figure 12).A new feature not seen for k D = 6 is found in the right panel for H = 0.4: here the IGWs are widely distributed across the domain.This is because the trapping effect weakens again as f /N increases above unity (here f /N = 2.771 approximately).Moreover, the IGWs are larger scale than found in the SW case.Their scale appears to be dictated by the wavenumber at which strong dispersion first occurs, k = √ 3/H which is about 4.33 for H = 0.4.This is consistent with the scale of the IGWs seen in the right panel of figure 12.
The time variation of the r.m.s.values of imbalanced divergence δ i and pressure p i are shown in figure 13.In δ i , the two values of H either side of H = 0.144338 (for which f /N = 1) exhibit the least imbalance, as expected, while H = 0.4 exhibits significantly larger imbalance (though not as large as for H → 0).The large value of H has a correspondingly large value of f /N, for which wave-trapping is weak.Similar results are found for p i , but now H = 0.1 exhibits the least imbalance (though only 2-3 times smaller than when H = 0.2).Here H = 0.4 exhibits the greatest imbalance (in p i ), slightly  greater than found for H → 0 at late times.But the key result is that we again find a suppression of IGW generation for f /N ∼ 1.
The distribution across scales of the imbalanced divergence is shown in figure 14, for both Rossby numbers and at the final time.The special value of H corresponding to f /N = 1 (yellow curve) shows a striking reduction in imbalance, particularly at large and intermediate scales.IGW suppression is clearly in evidence when f /N ∼ 1, but large H now stimulates large-scale IGWs, as already seen in p i in figure 12 (right panel).There is a systematic reduction of imbalance as f /N increases towards 1, followed by a growth as f /N increases beyond 1.
This suppression of imbalance is also seen in the frequency power spectrum of divergence shown in figure 15, again for Rossby numbers ε = 0.2 (left) and 0.6 (right), for various values of H as indicated.Recall that the vertical dashed lines mark the location of the Coriolis frequency f (cyan) and the buoyancy frequencies N (coloured the same way as the spectra).The bulges in power between f and N are consistent with IGWs being (largely) confined to these frequencies (this is strictly true only in linear theory).The confinement is poorer at larger Rossby number (right), but a substantial reduction in IGW activity still occurs, especially when H = 0.1 and H = 0.2, values for which f /N ∼ 1.

E. Larger scale initial conditions
We finally briefly consider a flow initialised at scales much larger than the Rossby deformation length.To this end, we discuss one case with k 0 = 3 and k D = 12 in the initial PV power spectrum (9), and a Rossby number ε = 0.6.We compare a VA simulation for H = 0.1 with a SW simulation (H → 0).The initial scale of the flow L is here four times larger than L D .This case has f /N = 1.2/ √ 3 = 0.6928....The flow evolution is slower than when k 0 = k D (considered previously), but inevitably 92 small-scale frontal features develop in PV as a result of nonlinear advection.These frontal features generate high frequencies as they move, and thus provide a source for the generation of IGWs by spontaneous adjustment emission.The slower pace of evolution is reflected in the values of the vorticity-based Rossby number |ζ | max / f and Froude number (|u|/ √ gh) max , which both remain around 0.15 (ranging from 0.13 to 0.18).The depth anomaly h ranges from −0.38 in cyclones to 0.94 in anti-cyclones.
The key results are summarised in figure 16, showing the divergence spectrum at t = 500 on the left and the frequency spectrum on the right, for both the SW and the VA simulation.The divergence spectrum shows that the flow is well balanced across all scales, with the greatest imbalance at small to intermediate scales as expected.As above when k 0 = k D (see figures 10 and 14), the simulation consistently exhibits less imbalance than the SW one.This is not just true at the time shown, but also at all earlier times.Here the r.m.s.δ i is approximately two times smaller in the VA simulation than in the SW one (not shown).Likewise, the frequency power spectrum of divergence in the right panel shows that finite H reduces imbalance considerably, especially at high frequencies.Again this is consistent with the dispersion relation (14), which shows that, for finite H, IGWs are confined to a narrower range between f and N (in linear theory).Thus, even when the initial flow has a scale L much larger than L D , finite H reduces spontaneous adjustment emission.The most important factor is not the scale of the initial flow, but the frequency ratio f /N.

IV. CONCLUSIONS
This paper has considered non-hydrostatic effects which are normally neglected in shallow-water theory on the assumption that horizontal scales L are much larger than the mean fluid depth H. Practically, however, there always exist scales L < H, only these are often too small to model directly in numerical simulations.The hydrostatic approximation no longer holds for these scales, and a generalised model is required.One can resort to a fully three-dimensional model, 49 but this is considerably more costly than the shallow-water model.An alternative is to directly average the three-dimensional model, assuming only that the horizontal flow is independent of depth as in traditional shallow-water theory.By making no other approximation, one arrives at a non-hydrostatic form of the shallow-water model, first derived (without rotation) by Serre 50 in 1953, and re-derived by numerous authors since.We call this model the vertically-averaged (VA) model.In a previous study 49 it was shown to be significantly more accurate than the traditional shallow-water model when compared to solutions of the full three-dimensional Euler equations with a free surface.
Here we have examined how non-hydrostatic effects in the VA model modify inertia-gravity wave (IGW) emission and propagation, an important topic in atmospheric and oceanic fluid dynamics. 22,23,46The key parameter is the Coriolisbuoyancy frequency ratio f /N where N ≡ 3g/H and g is the gravity (or reduced gravity for applications to the upper ocean 18,93 or the lower atmosphere 94 ).Linear waves on a basic state of rest have frequencies between f and N.In the traditional shallow-water model, N → ∞ and so there is no upper limit to the frequency of IGWs; moreover such highfrequency waves are non-dispersive, with constant phase and group velocities.In the VA model, by contrast, waves are dispersive and in particular their group velocities tend to zero with wavelength.As f /N → 1, the maximum group velocity vanishes, implying that waves at all scales are trapped: they cannot propagate away from their source (in linear theory).A wide range of numerical simulations of the nonlinear VA equations demonstrate that this wave-trapping effect can have a profound impact on both the emission and the propagation of IGWs.For values of f /N near unity, IGW emission weakens considerably compared to that occurring in the traditional shallow-water model, and moreover the waves remain largely confined to active regions of circulation (intense vortices, especially anti-cyclones).Their spatial form is also modified: the characteristic wave scale L w increases with H, and this scale is where non-hydrostatic dispersive effects become important (specifically, dispersion becomes strong for wavenumbers k > √ 3/H ≈ π/L w ).Even when f /N 1, as is typically found in the Earth's atmosphere and oceans 18,75 , non-hydrostatic effects should not be neglected at horizontal scales L comparable to the mean or nominal fluid depth.While the non-hydrostatic model is often considered to be more complicated, it has the advantage that it suppresses short-scale high-frequency motions, physically, without the need for numerical damping.Another advantage is having an upper frequency limit, which means larger time steps can be taken in numerical modelling.
The frequency ratio f /N is not uniformly small in geophysical flows.Regions of weak stratification (small N) exist in various parts of the oceans, including the Mediterranean Sea 73,74 , deep sea trenches 95 and in polar oceans 75 .Experimental studies of rotating stratified flows 96 also commonly study regimes where f /N ∼ 1.While a vertically-averaged model of the fluid motion may be too idealised, it at least offers a relatively straightforward way to comprehensively examine non-hydrostatic effects within a simplified framework.Here, we have found that such effects can greatly reduce inertia-gravity wave emission, and can lead to wave trapping for f /N ∼ 1.
In future work, we would like to explore the role of bottom topography, 58,97 specifically how the topographic generation of IGWs is altered by non-hydrostatic effects.Another extension would be to include a background planetary vorticity gradient to model the variation of the Coriolis frequency with latitude, or to consider full spherical geometry.One can also include the full Coriolis force, not just the component aligned with gravity. 55,70A two-layer model would additionally permit the study of baroclinic processes (vertical shear).][100] A different starting assumption is required that ensures continuity of the horizontal velocity at the layer interface, and this requires at least a linear dependence on height. 101In short, there is significant scope for further research.

FIG. 5 .
FIG.5.Maximum group velocity c g,max relative to c (black), and the wavenumber k max of maximum group velocity scaled in two ways (blue and red) as indicated, all as a function of the frequency ratio f /N.

FIG. 9 .
FIG. 9. Power spectra of divergence δ (black), balanced divergence δ b (blue) and imbalanced divergence δ i (red, as labelled) for both H → 0 (SW, dotted lines) and H = 0.1 (VA, solid lines).The panel on the left is for Rossby number ε = 0.2 while that on the right is for ε = 0.6.Both are for flows with k D = 6 at the final time t = 500.Note that the black curves often lie beneath the blue ones.

FIG. 10 .FIG. 11 .
FIG. 10.Power spectra of imbalanced divergence δ i for various values of H (as indicated) for Rossby numbers ε = 0.2 (left) and ε = 0.6 (right).Both are for flows with k D = 6 at the final time t = 500.

FIG. 14 .
FIG. 14.Power spectra of imbalanced δ i for various values of H (as indicated) for Rossby numbers ε = 0.2 (left) and ε = 0.6 (right).Both are for flows with k D = 12 at the final time t = 500.The analogous results for k D = 6 are shown in figure 10.

FIG. 15 .
FIG. 15.Frequency power spectra of divergence δ for various values of H (as indicated) for Rossby numbers ε = 0.2 (left) and ε = 0.6 (right).Both are for flows with k D = 12.The analogous results for k D = 6 are shown in figure 11.

1 FIG. 16 .
FIG.16.Left panel: power spectra of divergence δ (black), balanced divergence δ b (blue) and imbalanced divergence δ i (red, as labelled) at t = 500 for both H → 0 (SW, dotted lines) and H = 0.1 (VA, solid lines).Right panel: frequency power spectra of divergence δ , comparing the SW and VA simulations (the cyan vertical dashed line corresponds to the Coriolis frequency f and the blue one corresponds to the buoyancy frequency N in the VA simulation).Both panels are for Rossby number ε = 0.6 and k D = 12, and for large-scale initial conditions having k 0 = 3.