Abstract

Gravity wave emission by geostrophically balanced flow is diagnosed in numerical simulations of lateral and vertical shear instabilities. The diagnostic method in use allows for a separation of balanced flow and residual wave signal up to fourth order in the Rossby number (Ro). While evidence is found for a small but finite gravity wave emission from balanced flow in a single-layer model with large lateral shear and large Ro, a vertically resolved model with moderate velocity amplitudes appropriate to the interior ocean hardly shows any wave emission. Only when static instabilities generated by the shear instability of the balanced flow are allowed can a gravity wave signal similar to the ones reported in earlier studies be detected in the vertically resolved case. This result suggests a relatively small role of spontaneous wave emission in the classical sense of Lighthill radiation, and emphasizes the role of convective or symmetric instabilities during frontogenesis for the generation of internal gravity waves in the ocean and atmosphere.

1. Introduction

Gravity waves are a ubiquitous feature of the ocean and are generated by a variety of processes. One of them, the generation by spontaneous wave emission (e.g., Vanneste 2013), has recently attracted considerable attention. The process refers to the generation of gravity waves by balanced flow spontaneously in the absence of any external forcing, which appears to happen during baroclinic or barotropic instability (e.g., Plougonven and Snyder 2007; Hien et al. 2018; Chouksey et al. 2018), and is thought to be similar to sound generation in aerodynamics (Lighthill 1978; Ford et al. 2000), a process referred to as Lighthill radiation. Although it is intuitively clear to a physical oceanographer what balanced flow and gravity waves are, the exact definition of slow balanced flow and fast gravity waves is surprisingly difficult, since it is complicated by the presence of the nonlinear terms coupling both motions with each other. It is known that such nonlinearities can generate wavelike structures, which are nevertheless part of the balanced flow, the so-called ageostrophic balanced modes or slaved modes (e.g., Warn et al. 1995; McIntyre and Norton 2000; Kafiabad and Bartello 2018). The wavelike structures often seen in numerical simulations of lateral and vertical shear instabilities of the balanced flow could therefore be misinterpreted as spontaneous wave emission by the balanced flow, where in fact they are part of the balanced flow without the occurrence of spontaneous emission.

The most familiar model for the balanced mode is the result of the quasigeostrophic approximation, a first-order asymptotic expansion in the limit of small Rossby (Ro) and Froude number. Since the zero-order geostrophic velocity is free of horizontal divergence in that expansion, any significant horizontal divergence or vertical velocity seen in model simulations or observations is often interpreted as a gravity wave signal. It is, however, clear that the first-order ageostrophic vertical velocity—which can be calculated, for example, from the so-called omega equation (Hoskins et al. 1978)—although of first order in Ro and thus usually small, could also contribute to the observed vertical velocity. The ageostrophic vertical and horizontal first-order velocities in the quasigeostrophic approximation correspond to the ageostrophic balanced mode or slaved wave mode (to first order in Ro), and it is shown in Warn et al. (1995) how this concept of balanced models and wave modes can be extended to higher orders in Ro. This concept allows in principle to uniquely define the slow balanced flow and to differentiate the “true” fast gravity wave signal from that to arbitrary precision in terms of orders of Ro. The limit for infinitive order is called the slow manifold (and the difference to the full state the fast manifold), but it remains unclear if this limit exists (Lorenz 1980; Leith 1980; Lorenz 1986; Ford et al. 2000). In fact, for simple analogs to the equation of motions of the ocean and atmosphere it can be demonstrated that such an asymptotic limit does not exist (Vanneste 2013). In that case, the expansion in Ro is expected to diverge at large orders of Ro, although the lower orders of the expansion would still be useful for practical analysis.

In practice, however, the differentiation of waves and balanced flow in numerical models turns out to be rather difficult even at low order, such that the issue of the actual existence of the slow manifold appears to be of less concern. McIntyre and Norton (2000) applied a method to construct balanced models valid to higher order in Ro to a numerical model of the shallow water equations on a (half) sphere. They find that the balanced model simulation agrees surprisingly well with a model simulation using the full equations including the waves, such that they conclude that spontaneous wave emission remains unimportant in that simulation. To our knowledge, similar higher-order balance models have not yet been applied to numerical simulations of vertically resolved models. Therefore, we perform here a similar analysis as in McIntyre and Norton (2000) but for simulations of lateral shear instability in an idealized configuration applicable to the ocean for a range of values of Ro using a novel decomposition method at the discrete level of the model valid up to fourth order in Ro, and also extend this to the more realistic case of a vertically resolved model having also vertical shear instability.

In the next section, the theoretical background of the concept to differentiate waves from the balanced flow is outlined, while in section 3 the method is applied to numerical models of lateral and vertical shear instability. The last section provides a summary and discussion of the results.

2. Theoretical background

The theoretical background for the method to diagnose the model simulations following Warn et al. (1995) and Kafiabad and Bartello (2018) is outlined in this section. The focus is here on the analytical formulation, but since it is found that the numerical evaluation of the balanced modes needs exact consistency also on the discrete level of the numerical models, all relations on the discrete level which are actually used for the diagnosis are given in the  appendix.

a. The model

We consider the equations of motions in hydrostatic and Boussinesq approximation for a constant Coriolis parameter f and constant background stability frequency N. The resulting primitive equations are scaled using 1/f, L, H, and U as time scale, horizontal and vertical length scale, and horizontal velocity scale, respectively. The vertical velocity scale follows from the continuity equation, the pressure scale from geostrophic balance, and the buoyancy perturbation scale from hydrostatic balance. This yields1

 
tu+u¬+p=Ro(uu+wzu),tb+w=Ro(ub+wzb)
(1)

together with the diagnostic relations zp=b and u+zw=0. The only parameter in Eq. (1) controlling the nonlinear terms is the Rossby number Ro=U/(Lf) by setting Ro=Fr, where Fr is the Froude number Fr=U/(NH). By this setting, the (unscaled) first baroclinic Rossby radius Lr=NH/f becomes Lr=L. The case Ro1 is appropriate to the stratified interior ocean, while a larger Ro approaching unity is more appropriate to the so-called submesoscale flow which occurs for instance close to or in the mixed layer of the ocean. The scaled linearized kinematic surface and bottom boundary conditions become

 
tp|z=0=gh0udz,wz=h=0
(2)

with the scaled water depth h, the parameter g=L02/L21 and the (unscaled) barotropic Rossby radius L0=H(9.81m s2)/f. Vertically integrating the buoyancy equation and the hydrostatic relation yields

 
tp+M(u)=Roz0dz(ub+wzb)
(3)

with the integral operator M=gh0dz+z0dzhzdz. Equation (3) can be used to replace the hydrostatic relation and the buoyancy equation in Eq. (1).

b. The vertical modes

The operator M yields the vertical or baroclinic modes by its scaled eigenvalues cn2 and eigenfunctions ϕn. Taking zz of the eigenvalue equation M(ϕn)=cn2ϕn yields

 
cn2zzϕn+ϕn=0,
(4)

which is the equivalent eigenvalue equation formulated as differential equation. Considering zM(ϕ) at z=h and M(ϕ) and zM(ϕ) at z=0, the boundary conditions to the eigenvalue problem are derived as

 
ϕn|z=0+gzϕn|z=0=0,zϕn|z=h=0.
(5)

Note that by the surface boundary condition, the eigenvectors ϕn are not the ones for a rigid lid, but modified by the free surface and given by ϕn=Ancos(z+h)/cn with tanh/cn=cn/g. Approximate solutions of cn for g1 correspond to the different zero crossings of the tangent function and are c02gh and cn2h2/(nπ)2 for n>0. The constant amplitude An is chosen to satisfy the orthonormality condition h0ϕmϕndz=δn,m.

The eigenfunctions of the operator M allow to diagonalize the linear part of the system. The variables p and u are therefore expressed in terms of the eigenfunctions or vertical modes as

 
p=npn(x,y,t)ϕn(z),u=nun(x,y,t)ϕn(z)
(6)

and b=npnzϕn and w=ncn2unzϕn, which follows from the hydrostatic relation and the continuity equation. Using the expansion Eq. (6) in Eq. (7) and in the momentum equation from Eq. (1), multiplying with ϕm and integrating over depth yields

 
tun+un¬+pn=Ronun,tpn+cn2un=Ronpn
(7)

with nun=h0dzϕn(uu+wzu) and npn=h0dzϕnz0dz(ub+wzb). Equation (7) has in principle the same algebraically structure as a single-layer model, since the different vertical modes n are only coupled by the nonlinear terms nu and np. The equivalent nonlinear terms for a single-layer model with layer thickness pn and layer velocity un are nun=unun and npn=pnun, which would replace the right hand sides in Eq. (7).

c. The spectral representation

Assuming a double-periodic domain and applying the Fourier ansatz

 
un(x,t)=dku^n(k,t)eikx,pn(x,t)=dkp^n(k,t)eikx
(8)

with wavenumber vector k=(kx,ky), yields in Eq. (7) after multiplication with eikx and integration over x,

 
tz^iAz^=Ron^,A=(0ikxi0kxcn2kxcn2ky0),
(9)

with the state vector z^=(u^,p^)T, the vector of the nonlinearities n=(nu,np)T, and its Fourier transform n^=(2π)2dxneikx. The index n denoting the baroclinic mode number is omitted in Eq. (9) for simplicity. The matrix A has eigenvalues ωs and right and left eigenvectors qs=(qus,qυs,qps)T and ps=(pus,pυs,pps), respectively. The index s denotes here either the slow linear geostrophic mode (for s=0) or the fast linear gravity wave mode (for s=±). For the latter there are two linear waves modes with s=+ and s=. The eigenvalues ωs are given by the geostrophic mode ω0=0 and the fast mode ω±=±1+cn2|k|2. The latter is also the dispersion relation of long gravity waves in the system.

It is useful to project Eq. (9) on the fast or slow modes using the eigenvectors of A. The fast or slow mode amplitude gs=psz^ is then governed by

 
tgsiωsgs=Ropsn^=iRoIs(g0,g±),s=0,±,
(10)

where the case s=0 corresponds to the slow geostrophic mode with amplitude g0, and s=± to the (two) fast gravity wave mode with g+ and g. For Ro=0, that is, in the linear case, Eq. (10) describes the fast oscillation in time of the gravity waves g±=eiω±t, and the stationary solution for the slow mode g0=const. The nonlinear terms in n^ couple the modes gs at different k, but also the slow and fast modes which is expressed here by the nonlinear interaction integral Is(g0,g±).

d. The interaction integral

The exact form of the interaction integral Is is in principle not relevant here since it will not be evaluated directly—which is also possible, see, for example, Eden et al. (2019)—but calculated from numerical models. For the expansion into orders of the small parameter Ro which is performed below, however, we need to know about the general structure of the interaction integral. For the single-layer model with nu=uu and np=pu, the nonlinear terms become

 
n^(k)=idk1dk2Nδ(k1+k2k),N={[u^(k1)k2]u^(k2)[u^(k1)k2]υ^(k2)p^(k2)u^(k1)(k1+k2)},
(11)

and the integral Is

 
Is(g0,g±)=ipsn^=dk1dk2s1,s2=0,±gs1(k1)gs2(k2)Ck,k1,k2s,s1,s2δ(k1+k2k)
(12)

for the modes s=0,±. The so-called interaction coefficients C result from the projection of N on the eigenvector ps and are given in Eden et al. (2019). For the primitive equation model the structure of the integral Is is similar to Eq. (12), but a double sum over the vertical modes will add and the interaction coefficient are different [and also given for the rigid lid case in Eden et al. (2019)]. In any case, the slow and fast mode amplitudes gs will show up quadratically in a double sum over the fast and slow modes s, which is true for all kinds of models with quadratic nonlinearities.

By resolving and rearranging the double sum of the slow and fast modes s1,s2, the integral Is in Eq. (12) can therefore also be written as

 
Is(g0,g±)=Is(g0,0)+Is(0,g±)+Ks(g0,g±)
(13)

with

 
Is(g0,0)=dk1dk2g0(k1)g0(k2)Ck,k1,k2s,0,0δ(k1+k2k),Is(0,g±)=dk1dk2s1,s2=±gs1(k1)gs2(k2)Ck,k1,k2s,s1,s2δ(k1+k2k),Ks(g0,g±)=2dk1dk2s2=±g0(k1)gs2(k2)Ck,k1,k2s,0,s2δ(k1+k2k),
(14)

and similar for the primitive equation model but involving sums over all vertical modes. Equation (13) and Eq. (14) are used below to assess the orders of Is(g0,g±) in terms of Ro expanding the amplitudes gs into a power series of Ro.

e. The balanced models

The spectral representation of the primitive equations given by Eq. (10) describes the evolution of the gravity wave amplitudes g± and the slow geostrophic mode g0, and how they are coupled by the interaction integral Is. Due to this coupling, the definition of a slow balanced flow must involve not only the linear geostrophic mode g0, but also the linear gravity wave amplitudes g±. A method to uniquely describe the balanced part of the flow in this framework to any order in Ro is described by Warn et al. (1995). The method originates from the problem of balanced model initialization in numerical weather prediction, pioneered by Machenhauer (1977) and Baer and Tribbia (1977). The task in that problem is to provide an initial condition for a model for which a subsequent integration shows no or only a weak fast adjustment by gravity wave activity. Such an initial condition is then considered to be a balanced state.

Warn et al. (1995) assume accordingly a state in which the gravity waves are initially vanishing and then weakly growing in amplitude, such that g±=O(Ro) and expand the gravity wave amplitudes accordingly as

 
g±=Rof1±+Ro2f2±+Ro3f3±+
(15)

The time derivative in Eq. (10) includes the fast gravity wave oscillations with frequency ω± and the slow growth and decay of the amplitudes g0,± of both slow and fast mode by the nonlinear interactions. Therefore, a fast and slow time scale is introduced with T=Rot* and t=RoT+t*. The slow mode g0 is a function of the slow time T only since ω0=0, and thus Eq. (10) becomes for the geostrophic mode s=0

 
RoTg0=iRoI0(g0,g±).
(16)

Since the wave mode g± is assumed to be of first order in Ro, from Eq. (13) and the form of the integrals in Eq. (14), it follows that the first order of Eq. (16) simply becomes

 
Tg0=iI0(g0,0).
(17)

This first-order balanced model is identical to the familiar (first order) quasigeostrophic approximation as shown by Leith (1980). Only the slow mode amplitude g0 is involved in this first-order balanced model, and in this way Eq. (17)—which is a spectral representation of the quasigeostrophic potential vorticity equation—is closed. The second order of Eq. (10) for the geostrophic mode s=0 becomes

 
0=K0(g0,f1±)=I0(g0,f1±)I0(g0,0)I0(0,f1±),
(18)

using again Eq. (13) and thus with Eq. (17)

 
Tg0=iI0(g0,f1±)+iI0(0,f1±).
(19)

Together with the diagnostic expression for the auxiliary wave amplitude f1±, which is derived below, Eq. (19) is a second-order balanced model. Since the expression for f1± involves the geostrophic mode amplitudes g0 only, the second-order balanced model is also closed. Third-order yields accordingly a third-order balanced model given by

 
Tg0=iI0(g0,f2±)+iI0(0,f2±)iI0(0,f1±)
(20)

and so on for higher order. Again, the expression for f2± given below involves the geostrophic mode amplitudes g0 only. If the asymptotic expansion in Ro would converge, balanced models to arbitrary order of Ro could in principle be constructed, which are closed with respect to the geostrophic mode g0. The remaining task to find the auxiliary wave modes fn± is demonstrated in the next section. We will define the deviation from an actual state from the state given by g0 together with the auxiliary wave modes fn±, as the nonbalanced or fast gravity wave mode.

f. The ageostrophic balanced modes

The auxiliary wave modes fn± are also called ageostrophic balanced modes or slaved wave modes (Warn et al. 1995; McIntyre and Norton 2000; Kafiabad and Bartello 2018) and are part of the balanced mode since they evolve only slowly as seen below. The lowest order of those modes, f1±, corresponds to the first-order (ageostrophic) variables in the quasigeostrophic approximation (Leith 1980). The ageostrophic balanced modes are not needed to predict the evolution of the geostrophic variables given by g0 to the first order, that is, they are not needed and in general unknown in the quasigeostrophic model, but for the higher-order balanced models they are needed and diagnostic relations for fn± need to be derived. To first order in Ro, Eq. (10) becomes for the gravity wave modes s=±

 
t*f1±iω±f1±=iI±(g0,0),
(21)

using again Eq. (13) and Eq. (14) to infer the first order in Ro of the interaction integral I±. The time derivative in Eq. (21) generates the fast wave oscillation with frequency ω±. To suppress this fast oscillation in the balanced state, it is necessary to balance f1± as

 
f1±=I±(g0,0)/ω±.
(22)

This balancing of f1± yields t*f1±=0 in Eq. (21) and thus no waves are generated. Equation (22) is therefore the diagnostic relation which completes the second-order balance model in Eq. (19). It was shown by Leith (1980), that the ageostrophic balanced mode f1± corresponds to the first iteration step in the balanced initialization technique by Machenhauer (1977), and that f1± also corresponds to the ageostrophic variables in the (first order) quasigeostrophic approximation.

Setting t*fn±=0 to suppress wave activity in general, Eq. (10) for s=± becomes

 
n=1(Ron+1Tiω±Ron)fn±=iRoI±(g0,0)iRoI±(0,n=1Ronfn±)in=1Ron+1K±(g0,fn±),
(23)

using again Eq. (13). From the form of the integrals in Eq. (14) it becomes clear that the second, third, and fourth orders are given by

 
Tf1±iω±f2±=iK±(g0,f1±)Tf2±iω±f3±=iI±(0,f1±)iK±(g0,f2±)Tf3±iω±f4±=iI±(0,f1±+f2±)+iI±(0,f1±)+iI±(0,f2±)iK±(g0,f3±).
(24)

From the first line of Eq. (24), f2± is calculated, from the second line f3±, etc. Since only slow time derivatives show up, the ageostrophic balanced modes fn± are only slowly evolving in time as the geostrophic mode. The combination of geostrophic mode amplitude g0 and fn± given by g0q0+iRoifi±q± defines the balanced mode in spectral space, and inverse Fourier transform yields the balanced flow in physical space

In practice, we evaluate the interaction integral I± from a numerical model to find fn± up to fourth order. The term K±(g0,fn±) is given by Eq. (13) replacing g± with fn± and found by evaluating the integrals I± using the numerical model. To find Tf1±, Eq. (22) is analytically differentiated with respect to T as in Kafiabad and Bartello (2018), which yields

 
Tf1±=TI±(g0,0)/ω±=[I±(g0+Tg0,0)I±(g0,0)I±(Tg0,0)]/ω±
(25)

using the relation I±(g0+Tg0,0)=I±(g0,0)+TI±(g0,0)+I±(Tg0,0), which derives from the definition of the integral Is, and using Tg0=iI0(g0,0). To find Tf2± and Tf3±, the corresponding expressions get too complicated, and f2± and f3± are numerically differentiated in time by integrating the full model a few time steps as suggested by Baer and Tribbia (1977).

g. The effect of damping

In the numerical model it is necessary to include damping to suppress the gridscale noise generated by artificial numerical wave dispersion and the finite representation of the interaction integral Is. The primitive equations are therefore equipped with biharmonic operators A4u and A4b, with a viscosity (or diffusivity) A~O(Ro) to damp the gridscale noise. The same parameter A is used for momentum and buoyancy since then the only effect of the numerical damping is that the matrix A gets identical entries iA|k|4 on the diagonal [instead of zero for the case of vanishing damping as in Eq. (9)], the eigenvectors will stay identical as for the case A=0, and the eigenvalues become complex with identical real part as for A=0, but with imaginary part A|k|4.

There is however a complication: in Eq. (3), a term A4p shows up—which enters the matrix A on the diagonal in the spectral representation—but also a damping term A4p|z=0 involving the surface pressure. By this term, the vertical modes get coupled not only by the nonlinear terms but also by the numerical damping. To avoid this coupling, the surface boundary condition is also equipped with damping and Eq. (2) is replaced by

 
tp|z=0=gh0udzA4p|z=0.
(26)

Vertical friction and mixing are not used here, since we find it difficult to avoid the coupling of the vertical modes in the discrete version of the equations. However, the lateral biharmonic damping is found to be sufficient to control gridscale noise in the numerical integrations.

Introducing the numerical damping terms, the eigenvalues ωs obtain an imaginary part of O(Ro) which needs to be taken into account in the expansion. The first to fourth orders of the fast mode amplitude Eq. (10) then become

 
R(ω±)f1±=I±(g0,0)Tf1±iR(ω±)f2±=iK±(g0,f1±)f1±I(ω±)/RoTf2±iR(ω±)f3±=iI±(0,f1±)iK±(g0,f2±)f2±I(ω±)/RoTf3±iR(ω±)f4±=iI±(0,f1±+f2±)+iI±(0,f1±)+iI±(0,f2±)iK±(g0,f3±)f3±I(ω±)/Ro.
(27)

Since the magnitude of I(ω±)=A|k|4 strongly depends on wavenumber, we find that the damping effect can in practice become zero order in Ro for large k. In that case, Eq. (22) and Eq. (24) are used replacing Eq. (27) to determine fn± for large k, but using the complex ω±.

3. Numerical experiments

Numerical simulation of shear instabilities of balanced flow are diagnosed in terms of ageostrophic balanced modes given by Eq. (27). A single-layer model and a primitive equation model are used. Although less realistic than the primitive equation model, we use the single-layer model here because of its reduced complexity for which the method can be easier applied and tested. Further, it is possible to generate larger lateral shear in the initial conditions in the single-layer model, since in the primitive equation model too strong shear can generate convective or symmetric instabilities during the subsequent integration, even when the initial conditions are convectively stable.

The spatial discretization of the primitive equation model is outlined in the  appendix, and it follows the example of standard C-grid ocean models such as the MITgcm (Marshall et al. 1997). The discretization of the single-layer model is similar and detailed in Eden et al. (2019). The time-stepping scheme in both models is Euler forward for the first time step and afterward an Adam–Bashforth two-time level interpolation with adjusted weights to allow for a stable integration. In the primitive equation model, all vertical modes including the barotropic mode are integrated explicitly with sufficiently small time step, since we find that using an implicit method to solve for the barotropic mode to allow for a larger time step, the diagnostic of the higher-order ageostrophic balanced modes becomes inaccurate. The time step in the model is chosen according to the criteria given in the  appendix.

a. The wave diagnostic

The following method to differentiate gravity waves from balanced flow in the model solution is applied: any given model state (u,p) from the primitive equation model is projected on the vertical eigenfunctions using the discrete version of the eigenvalue problem appropriate to the numerical model given in the  appendix. For the single-layer model, this first step is not necessary. Then an horizontal Fourier transform is applied to (un,pn) and the spectral amplitudes are projected on the left eigenvector p0 of the discrete version of the system matrix A to obtain the slow mode amplitude g0. The inverse (vertical and horizontal) transform of g0q0—where q0 is the right eigenvector of A—is then used to calculate the nonlinear tendency terms in the numerical model, which yields I±(g0,0) after forward transform and projection on p±, and thus finally f1± from Eq. (27). Corresponding steps yield the other ageostrophic balanced modes fi±. The inverse (vertical and horizontal) transform of g0q0+i=1KRoifi±q± yields the balanced part of order K of the model state. The difference to the actual model state is interpreted as the residual gravity wave signal.

b. The single-layer model

The single-layer model on a double periodic domain size of 10 × 5 is initialized with an unstable zonal flow with an initial profile of u(y)=(2exp[(y1.25)2/0.12(y3.75)2/0.12],0) to which a small random perturbation is added. The initial conditions are balanced to first order, higher-order balances do not change the results discussed here. Since we use scaled equations as detailed after Eq. (1), the scaled Coriolis parameter is set to one and we also choose c=1, such that the single layer corresponds to the first baroclinic mode using the scaled depth h=π. The domain is resolved with 400 × 200 grid points and we choose biharmonic friction and mixing operators as detailed in section 2g with diffusivity A=(Ro/2)Δx4. The zonal jet begins to meander and dissolves finally into eddies in the simulations. The length scale of the eddies is close to the scaled (first baroclinic) Rossby radius c=1, which determines our choice of the model domain. Increasing the resolution of the model domain therefore does not change much the results presented here since the lateral shear instability is well resolved.

We present simulations with the single layer using different parameters, from Ro=0.02 to Ro=0.3, that is, ranging from mesoscale conditions with small Ro to the so-called submesoscale conditions with finite Ro. Table 1 presents statistics of the local Rossby number Rol=Ro|xυyu| after the onset of the instability in the simulations, showing that locally much larger Rol can be reached in the simulations than specified by the parameter Ro.

Table 1.

Spatial mean Rol¯, spatial maximum Rom=max(Rol), and mean of the largest 5% Ro5% of the actual local Rossby number Rol=Ro|xυyu| in the simulations with the single-layer model for different parameters Ro after the onset of the instability.

Spatial mean Rol¯, spatial maximum Rom=max⁡(Rol), and mean of the largest 5% Ro5% of the actual local Rossby number Rol=Ro|∂xυ−∂yu| in the simulations with the single-layer model for different parameters Ro after the onset of the instability.
Spatial mean Rol¯, spatial maximum Rom=max⁡(Rol), and mean of the largest 5% Ro5% of the actual local Rossby number Rol=Ro|∂xυ−∂yu| in the simulations with the single-layer model for different parameters Ro after the onset of the instability.

Figure 1 shows a time series of the ageostrophic balanced zonal velocity and the residual wave signal up to fourth order and the full thickness p at x=0 in a simulation with Ro=0.1. With the onset of the instability at t10, the ageostrophic balanced mode increases in amplitude in all orders. Figure 1a shows the zonal velocity of the first-order ageostrophic mode u1 and Fig. 1b shows uu0. We call u0 and ui the zonal velocity corresponding to the inverse transform of g0q0 and Roifi±q±, respectively, while u denotes the actual zonal velocity in the model. Note that u1 is the first-order ageostrophic velocity in quasigeostrophic approximation. Since the difference between u1 and uu0 is very small, there is no residual wave signal in the first order. The same is true for the second order, shown in Figs. 1c and 1d, and also for the third order (Figs. 1e,f).

Fig. 1.

(top) The ageostrophic balanced zonal velocity and (bottom) residual wave signal up to fourth order at x=0 for Ro=0.1 in the single-layer model. The terms u0 and ui denote the zonal velocity corresponding to the inverse transform of g0q0 and Roifi±q±, respectively. Also shown is the full thickness p as contour lines. Note the differences in the color scales.

Fig. 1.

(top) The ageostrophic balanced zonal velocity and (bottom) residual wave signal up to fourth order at x=0 for Ro=0.1 in the single-layer model. The terms u0 and ui denote the zonal velocity corresponding to the inverse transform of g0q0 and Roifi±q±, respectively. Also shown is the full thickness p as contour lines. Note the differences in the color scales.

Only in the fourth order (Figs. 1g,h) do small differences between u4 and uu0u1u2u3 begin to show up, which points to a small residual gravity wave signal in the fourth order. This residual wave signal increases going to larger Ro. Figure 2 shows that the wave signal begins to emerge for Ro=0.15 already in the third order, and that at fourth order a clear large-scale near-inertial wave signal is present.

Fig. 2.

As in Fig. 1, but for Ro=0.15.

Fig. 2.

As in Fig. 1, but for Ro=0.15.

Figure 3a shows the residual wave energy EK=|un=0Kun|2/2+(pn=0Kpn)2/2 as a function of Ro normalized with the total energy |u|2/2+p2/2 (which does not vary much for different Ro) during the peak of the instability in the simulations. The first-order residual energy E1 scales as Ro4 such that all the energy can be attributed to the second-order ageostrophic balanced mode for all Ro, since the amplitude of the second-order ageostrophic balanced mode scales as Ro2. Correspondingly, E2 scales roughly as Ro6 and E4 as Ro8 for small Ro, but this scaling behavior gets distorted for larger Ro due to the emergence of the residual gravity wave signal seen in Fig. 2. An exponential fit ~exp(2/Ro) seems to fit the residual gravity wave signal. For small Ro, however, the exponential fit deviates, and also E4 does not scale as Ro10 for small Ro. This points toward a problem in the diagnosis of the fourth-order ageostrophic balanced mode with our method, most likely due to numerical precision errors. In any case, the diagnostics reveals a finite gravity wave signal for Ro>0.1. This wave signal could be related to spontaneous gravity wave emission generated by shear instability of the balanced flow, but note that local Rossby numbers Rol in Table 1 exceed unity such that symmetric instability due to lateral shear becomes possible. Therefore, the wave signal could also relate to this instability process instead of spontaneous gravity wave emission in the classical Lighthill sense.

Fig. 3.

(a) Total residual wave energy EK=|un=0Kun|2/2+p(n=0Kpn)2/2 for K=1,,4 normalized with the total energy integrated over the model domain at the peak of the instability in the single-layer model. The term E1 corresponds to the uppermost dots and E4 to the lowest dots. (b) As in (a), but for the primitive equation model and with total residual wave energy EK=|un=0Kun|2/2+(bn=0Kbn)2/2 also normalized with the total energy. Dashed lines show different power laws Ro4, Ro6, and Ro8, and the dotted line in (a) an exponential scaling law exp(2/Ro).

Fig. 3.

(a) Total residual wave energy EK=|un=0Kun|2/2+p(n=0Kpn)2/2 for K=1,,4 normalized with the total energy integrated over the model domain at the peak of the instability in the single-layer model. The term E1 corresponds to the uppermost dots and E4 to the lowest dots. (b) As in (a), but for the primitive equation model and with total residual wave energy EK=|un=0Kun|2/2+(bn=0Kbn)2/2 also normalized with the total energy. Dashed lines show different power laws Ro4, Ro6, and Ro8, and the dotted line in (a) an exponential scaling law exp(2/Ro).

c. The primitive equation model

The primitive equation model is configured similar to the single-layer model on a double periodic domain of size 5 × 5 × 1 resolved by 255 × 255 × 40 grid points, in the zonal, meridional, and vertical direction, respectively. The model is initialized with a zonal velocity shaped in the horizontal identical to the single-layer model, but here with the first vertical mode as vertical structure. The initial velocity, pressure, and buoyancy fields are balanced to first order in Ro; higher-order balancing does not change the results. Since we use scaled equations, the (scaled) background stratification (or Brunt–Väisälä frequency) is one and the Coriolis parameter is also one. Compared to the single-layer model, the amplitude of the initial zonal velocity u is reduced to 0.2 to exclude the occurrence of static instabilities with zb1 that arise at higher amplitudes of u and show up locally in the simulations. The scaled first baroclinic Rossby radius is Lr=c1=1/π, which is as in the single-layer model close to the length scale of the instabilities showing up after finite time in the integrations, and also determines the choice of the model domain. Therefore, increasing the horizontal and vertical resolution does not change the integrations significantly and also does not affect the decomposition into balanced and unbalanced flow.

The scaled barotropic Rossby radius L0 is set to 5 by adjusting the parameter g=L02/Lr2 accordingly. Larger values of L0 do not change the simulations much, but require a much smaller time step, therefore this rather small ratio between L0 and Lr is chosen. The lateral damping is identical to the single-layer model. While we time step the model only once in the single-layer model to calculate Tfn± in Eq. (27), we find it necessary to calculate 50 time steps in the primitive equation model since otherwise large gridscale noise shows up in the higher-order fn±. Furthermore, it is necessary to exchange Eq. (27) with Eq. (22) and Eq. (24) but with complex ω± for the 40 largest wavenumbers, as discussed in section 2g.

As for the single-layer model, we present simulations for different Ro with a range from Ro=0.10 to Ro=0.60, that is, from mesoscale conditions with small Ro to the so-called submesoscale conditions with finite Ro. Table 2 presents statistics of the local Rossby number Rol=Ro|xυyu| after the onset of the instability in the simulations. Although smaller in general than in the single-layer model because of the reduced initial zonal velocity amplitude, also here larger Rol are reached locally in the simulations than specified by the parameter Ro.

Table 2.

As in Table 1, but at the surface of the primitive equation model. The last line shows the same quantities but for the simulation with a larger amplitude of 1 in the zonal velocity of the initial condition instead of 0.2.

As in Table 1, but at the surface of the primitive equation model. The last line shows the same quantities but for the simulation with a larger amplitude of 1 in the zonal velocity of the initial condition instead of 0.2.
As in Table 1, but at the surface of the primitive equation model. The last line shows the same quantities but for the simulation with a larger amplitude of 1 in the zonal velocity of the initial condition instead of 0.2.

Figure 4 shows the ageostrophic balanced vertical velocity up to fourth order, that is, from w1 to w4 and the corresponding residual wave signals in the primitive equation model after the onset of the instability for Ro=0.5. Again, the first-order ageostrophic vertical velocity w1 is very similar to the actual vertical velocity, and the second-order balanced vertical velocity w2 is very similar to ww1, which is the second-order residual signal. This similarity between w2 and ww1 indicates that most of the residual signal seen in ww1 is in fact related to the balanced mode (and is the so-called ageostrophic balanced mode or slaved mode), and not to the unbalanced wave signal it appears to be. Note that w1 corresponds to the quasigeostrophic vertical velocity, which can also be diagnosed from the omega equation (not shown), and that w2 is an order of magnitude smaller but still entirely related to the balanced mode. There is also good agreement between w3 and ww1w2 but small differences already show up, while w4 deviates more from the corresponding residual wave signal. This behavior is also seen for smaller Ro, that is, agreement up to w3 but deviation for w4.

Fig. 4.

(top) Ageostrophic balanced vertical velocity (color; from w1 to w4) and (bottom) residual wave signal (from ww1 to wi=13wi) up to fourth order at z=0.2 in the primitive equation model for Ro=0.5. The term wi denotes the vertical velocity corresponding to the inverse transform of Roifi±q±. Shown is a time step after the onset of the instability. Also shown is the full pressure as contour lines. Note the different color scales.

Fig. 4.

(top) Ageostrophic balanced vertical velocity (color; from w1 to w4) and (bottom) residual wave signal (from ww1 to wi=13wi) up to fourth order at z=0.2 in the primitive equation model for Ro=0.5. The term wi denotes the vertical velocity corresponding to the inverse transform of Roifi±q±. Shown is a time step after the onset of the instability. Also shown is the full pressure as contour lines. Note the different color scales.

The apparent wave signal at fourth order (i.e., in ww1w2w3 in Fig. 4) for small Ro is most likely related to numerical precision errors and not to a wave signal. This can be seen in Fig. 3b in which the total residual energy at the fourth order (E4) shows a very similar scaling to that of the third order (E3). Since E3 contains mostly ageostrophic balanced modes (as seen in Fig. 4 for Ro=0.5 which also holds for different Ro), this indicates that the fourth-order residual energy E4 is also largely composed of the ageostrophic balanced modes. Thus, the wavelike signals seen in vertical velocity in Fig. 4h are only an apparent signal, arising most likely from a numerical artifact, which we have not been able to remove2, although we use the discretized equations consistent on the grid level to eliminate numerical inaccuracy (we found that small inconsistencies on the discrete level can lead to first-order inaccuracies in the flow decomposition). Similar to Fig. 3a, however, for Ro>0.5 a deviation from the corresponding power laws is seen in all En in Fig. 3b, but here for much larger Ro. These results leads us to conclude that the residual wave signal is very small in this configuration of the primitive equation model and begins to be detectable with our decomposition method only at very large Ro. Then it appears to follow also an exponential power law as in Fig. 3a. We therefore conclude that there is hardly any spontaneous wave emission by the shear instability of the balanced flow in this model simulation for Ro<0.5.

The situation changes when allowing for amplitudes of the zonal velocity in the initial conditions larger than 0.2. In this case, however, static instabilities with zb1 are generated by the growing shear instability, even when the initial conditions are stably stratified. This is not the case in the previous model experiment, where the amplitude was chosen as large as possible but without static instabilities in the simulations. For the experiment with larger velocity amplitudes, we thus apply a convective adjustment to the model with a vertical diffusivity, which takes a value of K=1 if zb<1 and K=0 elsewhere, but we do not account for this convective adjustment in the diagnostic of the ageostrophic balanced mode.

Figure 5 shows the ageostrophic balanced vertical velocity up to third order, w1, w2, and w3, and the corresponding residual wave signals (i.e., w, ww1, ww1w2) for an initial zonal velocity amplitude of 1 and for Ro=0.1. The actual local Rossby numbers Rol for this simulation are also listed in Table 2. They are larger than in the previous simulations with the primitive equation model avoiding static instabilities. Already in the vertical velocity w a wave signal is present which dominates in the second and higher orders over the balanced mode which is shown in Figs. 5a, 5c, and 5e. The wave crests in Figs. 5b, 5d, and 5f are oriented roughly perpendicular with a wavelength of about 0.2–0.3 with a larger vertical mode than the initial conditions. The emergence of the wave signal shows up for all Ro we have tested (0.02–0.3) and can be related to the regions of static instabilities with (or close to) zb1 in space and time. We therefore conclude that the waves we do detect in our simulations are actually generated by convective or symmetric instabilities occurring in regions where zb~1—or where potential vorticity becomes negative as criterion for symmetric instability (Olbers et al. 2012)—but not by spontaneous wave emission.

Fig. 5.

As in Fig. 4, but with larger amplitude of 1 in the zonal velocity of the initial condition and for Ro=0.1.

Fig. 5.

As in Fig. 4, but with larger amplitude of 1 in the zonal velocity of the initial condition and for Ro=0.1.

Experiment with the primitive equation model without strong lateral but vertical shear, where predominantly baroclinic instability shows up, behaves the same as the above discussed experiments. For amplitudes in vertical shear that do not generate static instabilities, hardly any wave signal is detected, only with increased shear and occasional static instabilities a clear wave signal similar to Fig. 5 is generated.

4. Summary and discussion

This study is an elaborate attempt to find spontaneous gravity wave emission during lateral and vertical shear instability of balanced flow. We have designed a novel numerical tool to uniquely differentiate the slow balanced flow from the fast gravity waves up to fourth order based on an asymptotic expansion in the Rossby number Ro as suggested in, for example, Warn et al. (1995) and Kafiabad and Bartello (2018). Here, we apply the method the first time up to fourth order and in a model of growing instabilities to answer the specific question of spontaneous wave emission and its dependency on Ro, whereas the method was applied to a model of decaying isotropic turbulence by Kafiabad and Bartello (2018) only up to second order in Ro. The concept is applied here on the discrete level of the numerical models since we find that this high accuracy is necessary to allow for orders of the expansion larger than one. In principle, our tool allows for arbitrary precision, in practice this limit is hampered by numerical (lack of precision) and mathematical (lack of asymptotic limit) issues. Unlike the available approximate methods to differentiate between waves and balanced flow, such as the Lagrangian frequency method (Shakespeare and Taylor 2015) or the optimal potential vorticity balance method (Viúdez and Dritschel 2004), a Fourier transform of the model state is necessary in Warn et al.’s approach. The application of our decomposition tool is therefore more straightforward to idealized settings, whereas the approximate methods can be more readily applied to realistic scenarios but at the cost of a less accurate decomposition and perhaps a misinterpretation of the wavelike signals. We plan a comparison of the performance of the different methods in idealized models to transfer the results and to use the methods in the realistic ocean models as well.

In a numerical simulation of lateral shear instability with a single-layer model our novel tool can be successfully applied up to fourth order, which allows us to diagnose the gravity wave signal to this order of accuracy. Only in the case of a strong lateral shear, a near inertial wave signal is detected at higher orders of Ro in the model simulations. This result of weak wave generation is in agreement to McIntyre and Norton (2000) who also reported a minor role of wave generation in a single layer appropriate to the global atmosphere. The small but finite wave generation at large Ro we found could be related to spontaneous gravity wave emission in the sense of Lighthill radiation, but since local Rossby numbers well exceed unity, symmetric instability due to the lateral shear could also be the cause of the wave signal. The amplitude of the wave signal, however, shows an exponential scaling with respect to Ro in agreement to previous suggestions of classical spontaneous wave generation (Vanneste and Yavneh 2004; Vanneste 2013).

As compared to the decomposition up to fourth order in Ro in the single-layer model, our novel tool can robustly identify wave signals at least up to third order in Ro in the vertically resolved model, where the reduction in accuracy is most likely related to the error in numerical precision. In the vertically resolved model with lateral and vertical shear small enough to avoid static instabilities, the diagnosis and the scaling behavior with respect to Ro suggest that there is a negligibly small wave signal. Only if the shear becomes strong enough to generate static instabilities the emission of waves is detected with lateral wavelengths of a fraction of the Rossby radius and a higher vertical mode than the balanced flow. This emission is localized in space and time to the occurrence of static instabilities and the waves closely resemble the wavelike signals often seen before in numerical simulations (e.g., Plougonven and Snyder 2007; Hien et al. 2018; Chouksey et al. 2018).

Relatively larger amplitudes of waves generated by convective or symmetric instabilities compared to the ones generated by spontaneous emission are also found and discussed in Chouksey (2018) in a suite of numerical experiments with different configurations and for a range of different Ro. There, however, the method to differentiate waves from balanced flow follows the one by Machenhauer (1977) or uses the omega equation, and is thus only accurate to first order in Ro, such that the amplitude of spontaneous wave emission is most likely overestimated.

We cannot rule out the possibility that in other model configurations than what we have discussed here (and in Chouksey 2018) spontaneous wave emission in the sense of Lighthill radiation plays a more important role. We have tested several different configurations including also ones with vertical shear as before but only weak lateral shear and always found the same result: hardly any wave generation in general, but enhanced wave generation when static instabilities are generated by the balanced flow. In all such cases, the wave signal is localized in space and time to the sites of static instability. We also cannot rule out the possibility that waves at much larger or smaller wavelengths than allowed by our model grid can be generated. On the other hand, the model and laboratory experiments in the literature refer to spontaneously generated waves at similar wavelength as the dynamical process under investigation, in agreement to the (true) wave signal seen in our experiments with convective and/or symmetric instabilities.

These results suggest a relatively small role of spontaneous wave emission in the classical sense of Lighthill radiation, and emphasize the role of convective or symmetric instabilities during frontogenesis for the generation of internal waves. We speculate that the internal wave signals seen and discussed in previous studies based on numerical simulations and laboratory experiments (e.g., Plougonven and Snyder 2007; Hien et al. 2018; Chouksey et al. 2018) may also be generated by convective or symmetric instabilities instead of the classical spontaneous emission mechanism.

Acknowledgments

This paper is a contribution to the Collaborative Research Centre TRR 181 “Energy Transfer in Atmosphere and Ocean” funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation), Projektnummer 274762653.

APPENDIX

Discrete Relations

The discretization of the primitive equations is given by

 
tuυ¯i+¯j+δi+p=nuA(δiδi++δjδj+)2u,tυ+u¯i¯j++δj+p=nυA(δiδi++δjδj+)2υ,tb+wk¯=nbA(δiδi++δjδj+)2b,
(A1)

and follows a standard C-grid arrangement. The indices i, j, k denote discretization in the x, y, and z directions, but shown only when they are relevant in the individual terms. Finite differencing operators δi+p=(pi+1pi)/Δx=δi+1p, where Δx is the grid spacing in x direction, and averaging operators pi¯+=(pi+pi+1)/2=pi+1¯ and similar for the other directions with grid spacing Δy and Δz are introduced. The vertical grid has grid points k=1,,N with pressure and horizontal velocity at depth level zk=Δz(Nk1/2), vertical velocity at depth level zk+Δz/2, and surface pressure ps at zN+Δz/2=0. The surface boundary condition is then given by

 
tps=gΔzk=1NukA(δiδi++δjδj+)2ps.
(A2)

As discussed above, the damping term is included here since otherwise the vertical modes get coupled by the damping. The nonlinear terms nu, nυ, and nb are given for the primitive equations by

 
nu=δi(ui¯+ui¯+)δj(υi¯+uj¯+)δk(wi¯+uk¯+),nυ=δi(uj¯+υi¯+)δj(υj¯+υj¯+)δk(wj¯+υk¯+),nb=δi(ubi¯+)δj(υbj¯+)δk(wbk¯+).
(A3)

For the single-layer model, the discretization of the nonlinear terms in the momentum equation follows the energy-conserving scheme by Sadourny (1975) and similar to nb in Eq. (A3) for the thickness equation.

Integrating the discrete buoyancy equation vertically and using the hydrostatic relation δk+p=b¯k+, the surface boundary condition, and the integrated continuity equation wk=Δzk=1kuk yields

 
tpk=M(uk)Δz(nkb/2+k=k+1k=Nnkb)A(δiδi++δjδj+)2pk
(A4)

with the discrete vertical mode operator

 
M(ϕk)=Δz(gΔz2)k=1Nϕk+Δz2k=kk=Nk=1kϕkΔz24ϕk.
(A5)

The discrete operator M has eigenvalues cn2 and eigenfunctions ϕkn with M(ϕkn)=cn2ϕkn. Taking two finite differences of the eigenvalue equation yields

 
(cn2+Δz24)δkδk+ϕn=ϕkn,(g+Δz2)δN+ϕn+ϕNn=0,δ1ϕn=0.
(A6)

Both the integral form and the solution to the differential form yield identical eigenfunctions and eigenvalues for the discrete grid. Expressing ui,j,k=nui,jnϕkn and pi,j,k=npi,jnϕkn and applying a horizontal Fourier transform yields the discrete version of the system matrix A of Eq. (9) given in Eden et al. [2019, their Eq. (50)], although here equipped with damping terms on the diagonal and for vertical mode number n. The eigenvectors and eigenvalues of A are therefore also identical to the ones in Eden et al. [2019, their Eqs. (51)–(53)] but with c=cn, and the eigenvalues also get an imaginary part related to the damping. The corresponding vertical mode expressions for the auxiliary variables b and w are wi,j,k=n(cn2+Δz2/4)unδk+ϕn and bi,j,k=npn[δkϕnΔz/(4cn2)(ϕkn+ϕk1n)].

The time stepping scheme solving Eq. (9) (i.e., its inverse Fourier transform) is

 
zn+1zn=Δt[iA(αzn+βzn1)+αn^(zn)+βn^(zn+1)]
(A7)

with time levels n, the time step Δt, and with the interpolation weights α=1.5+ϵ, β=(0.5+ϵ). The parameter ϵ>0 allows for a stable integration if Δt is chosen appropriately small. The value of Δt can be estimated from the eigenvalues of the linear discrete system ω¯s, which differ in general from the physical eigenvalues ωs and are given by

 
ω¯1,2=iln(z1,2)/Δt,z1,2=12iαΔt2ωs±1214iΔtωs(β+α2)(Δtωsα)2.
(A8)

The complex frequency ω¯1s corresponds to a spurious numerical mode, which is always strongly damped and cannot be seen in the integrations. The other frequency ω¯2s converges to ωs for Δt0, but for finite Δt it differs from ωs and becomes complex. For ϵ=0, I(ω¯2s)>0 for all wavenumbers such that the scheme becomes unstable. But the parameter ϵ>0 can be chosen for given Δt, such that I(ω¯2s)<0 for all vertical modes and lateral wavenumbers possible on the model grid. We choose ϵ=0.1 and adjust Δt accordingly. The integrations are then always stable.

REFERENCES

REFERENCES
Baer
,
F.
, and
J. J.
Tribbia
,
1977
:
On complete filtering of gravity modes through nonlinear initialization
.
Mon. Wea. Rev.
,
105
,
1536
1539
, https://doi.org/10.1175/1520-0493(1977)105<1536:OCFOGM>2.0.CO;2.
Chouksey
,
M.
,
2018
: Disentangling gravity waves from balanced flow. Ph.D. thesis, Universität Hamburg, https://10.17617/2.3000108.
Chouksey
,
M.
,
C.
Eden
, and
N.
Brüggemann
,
2018
:
2018: Internal gravity wave emission in different dynamical regimes
.
J. Phys. Oceanogr.
,
48
,
1709
1730
, https://doi.org/10.1175/JPO-D-17-0158.1.
Eden
,
C.
,
M.
Chouksey
, and
D.
Olbers
,
2019
:
Mixed Rossby–gravity wave–wave interactions
.
J. Phys. Oceanogr.
,
49
,
291
308
, https://doi.org/10.1175/JPO-D-18-0074.1.
Ford
,
R.
,
M.
McIntyre
, and
W.
Norton
,
2000
:
Balance and the slow quasimanifold: Some explicit results
.
J. Atmos. Sci.
,
57
,
1236
1254
, https://doi.org/10.1175/1520-0469(2000)057<1236:BATSQS>2.0.CO;2.
Hien
,
S.
,
J.
Rolland
,
S.
Borchert
,
L.
Schoon
,
C.
Zülicke
, and
U.
Achatz
,
2018
:
Spontaneous inertia–gravity wave emission in the differentially heated rotating annulus experiment
.
J. Fluid Mech.
,
838
,
5
41
, https://doi.org/10.1017/jfm.2017.883.
Hoskins
,
B.
,
I.
Draghici
, and
H.
Davies
,
1978
:
A new look at the ω-equation
.
Quart. J. Roy. Meteor. Soc.
,
104
,
31
38
, https://doi.org/10.1002/qj.49710443903.
Kafiabad
,
H. A.
, and
P.
Bartello
,
2018
:
Spontaneous imbalance in the non-hydrostatic Boussinesq equations
.
J. Fluid Mech.
,
847
,
614
643
, https://doi.org/10.1017/jfm.2018.338.
Leith
,
C. E.
,
1980
:
Nonlinear normal mode initialization and quasi-geostrophic theory
.
J. Atmos. Sci.
,
37
,
958
968
, https://doi.org/10.1175/1520-0469(1980)037<0958:NNMIAQ>2.0.CO;2.
Lighthill
,
J.
,
1978
: Waves in Fluids. Cambridge University Press, 516 pp.
Lorenz
,
E. N.
,
1980
:
Attractor sets and quasi-geostrophic equilibrium
.
J. Atmos. Sci.
,
37
,
1685
1699
, https://doi.org/10.1175/1520-0469(1980)037<1685:ASAQGE>2.0.CO;2.
Lorenz
,
E. N.
,
1986
:
On the existence of a slow manifold
.
J. Atmos. Sci.
,
43
,
1547
1558
, https://doi.org/10.1175/1520-0469(1986)043<1547:OTEOAS>2.0.CO;2.
Machenhauer
,
B.
,
1977
:
On the dynamics of gravity oscillations in a shallow water model with applications to normal mode initialization
.
Beitr. Phys. Atmos.
,
50
,
253
271
.
Marshall
,
J.
,
A.
Adcroft
,
C.
Hill
,
L.
Perelman
, and
C.
Heisey
,
1997
:
A finite-volume, incompressible Navier Stokes model for studies of the ocean on parallel computers
.
J. Geophys. Res.
,
102
,
5753
5766
, https://doi.org/10.1029/96JC02775.
McIntyre
,
M. E.
, and
W. A.
Norton
,
2000
:
Potential vorticity inversion on a hemisphere
.
J. Atmos. Sci.
,
57
,
1214
1235
, https://doi.org/10.1175/1520-0469(2000)057<1214:PVIOAH>2.0.CO;2.
Olbers
,
D.
,
J.
Willebrand
, and
C.
Eden
,
2012
: Ocean Dynamics. Springer, 703 pp.
Plougonven
,
R.
, and
C.
Snyder
,
2007
:
Inertia–gravity waves spontaneously generated by jets and fronts. Part I: Different baroclinic life cycles
.
J. Atmos. Sci.
,
64
,
2502
2520
, https://doi.org/10.1175/JAS3953.1.
Sadourny
,
R.
,
1975
:
The dynamics of finite-difference models of the shallow-water equations
.
J. Atmos. Sci.
,
32
,
680
689
, https://doi.org/10.1175/1520-0469(1975)032<0680:TDOFDM>2.0.CO;2.
Shakespeare
,
C. J.
, and
J.
Taylor
,
2015
:
The spontaneous generation of inertia–gravity waves during frontogenesis forced by large strain: Numerical solutions
.
J. Fluid Mech.
,
772
,
508
534
, https://doi.org/10.1017/jfm.2015.197.
Vanneste
,
J.
,
2013
:
Balance and spontaneous wave generation in geophysical flows
.
Annu. Rev. Fluid Mech.
,
45
,
147
172
, https://doi.org/10.1146/annurev-fluid-011212-140730.
Vanneste
,
J.
, and
I.
Yavneh
,
2004
:
Exponentially small inertia–gravity waves and the breakdown of quasigeostrophic balance
.
J. Atmos. Sci.
,
61
,
211
223
, https://doi.org/10.1175/1520-0469(2004)061<0211:ESIWAT>2.0.CO;2.
Viúdez
,
Á.
, and
D. G.
Dritschel
,
2004
:
Optimal potential vorticity balance of geophysical flows
.
J. Fluid Mech.
,
521
,
343
352
, https://doi.org/10.1017/S0022112004002058.
Warn
,
T.
,
O.
Bokhove
,
T.
Shepherd
, and
G.
Vallis
,
1995
:
Rossby number expansions, slaving principles, and balance dynamics
.
Quart. J. Roy. Meteor. Soc.
,
121
,
723
739
, https://doi.org/10.1002/qj.49712152313.

Footnotes

© 2019 American Meteorological Society. For information regarding reuse of this content and general copyright information, consult the AMS Copyright Policy (www.ametsoc.org/PUBSReuseLicenses).

1

All vectors are horizontal. The vector u¬ denotes anticlockwise rotation of the vector u by π/2, that is, u¬=(υ,u) for u=(u,υ).

2

Note that tests with balanced initializations of the primitive equation model without any subsequent instability also show the absence of a fast adjustment by gravity waves up to third order, but a fast wave signal due to an adjustment from the initial conditions in the fourth order, indicating an imperfect balanced initialization on this order. These integrations are, however, not presented here.