## 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 yields^{1}

together with the diagnostic relations $\u2202zp=b$ and $\u2207\u22c5u+\u2202zw=0$. The only parameter in Eq. (1) controlling the nonlinear terms is the Rossby number $Ro=U/\u2061(Lf)$ by setting $Ro=Fr$, where Fr is the Froude number $Fr=U/\u2061(NH)$. By this setting, the (unscaled) first baroclinic Rossby radius $Lr=NH/f$ becomes $Lr=L$. The case $Ro\u226a1$ 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

with the scaled water depth *h*, the parameter $g=L02/L2\u226b1$ and the (unscaled) barotropic Rossby radius $L0=H\u2061(9.81\u2009m s\u22122)/f$. Vertically integrating the buoyancy equation and the hydrostatic relation yields

### b. The vertical modes

The operator *M* yields the vertical or baroclinic modes by its scaled eigenvalues $cn2$ and eigenfunctions $\varphi n$. Taking $\u2202zz$ of the eigenvalue equation $M\u2061(\varphi n)=cn2\varphi n$ yields

which is the equivalent eigenvalue equation formulated as differential equation. Considering $\u2202zM\u2061(\varphi )$ at $z=\u2212h$ and $M\u2061(\varphi )$ and $\u2202zM\u2061(\varphi )$ at $z=0$, the boundary conditions to the eigenvalue problem are derived as

Note that by the surface boundary condition, the eigenvectors $\varphi n$ are not the ones for a rigid lid, but modified by the free surface and given by $\varphi n=An\u2061\u2009cos\u2061(z+h)/cn$ with $tanh/cn=cn/g$. Approximate solutions of $cn$ for $g\u226b1$ correspond to the different zero crossings of the tangent function and are $c02\u2248gh$ and $cn2\u2248h2/\u2061(n\pi )2$ for $n>0$. The constant amplitude $An$ is chosen to satisfy the orthonormality condition $\u222b\u2212h0\varphi m\varphi n\u2009dz=\delta 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

and $b=\u2211npn\u2202z\varphi n$ and $w=\u2211ncn2\u2207\u22c5un\u2202z\varphi 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 $\varphi m$ and integrating over depth yields

with $nun=\u2212\u222b\u2212h0dz\varphi n\u2061(u\u22c5\u2207u+w\u2202zu)$ and $npn=\u222b\u2212h0dz\varphi n\u222bz0dz\u2032\u2061(u\u22c5\u2207b+w\u2202zb)$. 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=\u2212un\u22c5\u2207un$ and $npn=\u2212\u2207\u22c5pnun$, which would replace the right hand sides in Eq. (7).

### c. The spectral representation

Assuming a double-periodic domain and applying the Fourier ansatz

with wavenumber vector $k=\u2061(kx,ky)$, yields in Eq. (7) after multiplication with $e\u2212ik\u22c5x$ and integration over **x**,

with the state vector $z^=\u2061(u^,p^)T$, the vector of the nonlinearities $n=\u2061(nu,np)T$, and its Fourier transform $n^=\u2061(2\pi )\u22122\u222b\u2009dxne\u2212ik\u22c5x$. The index *n* denoting the baroclinic mode number is omitted in Eq. (9) for simplicity. The matrix $A$ has eigenvalues $\omega s$ and right and left eigenvectors $qs=\u2061(qus,q\upsilon s,qps)T$ and $ps=\u2061(pus,p\upsilon 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=\xb1$). For the latter there are two linear waves modes with $s=+$ and $s=\u2212$. The eigenvalues $\omega s$ are given by the geostrophic mode $\omega 0=0$ and the fast mode $\omega \xb1=\xb11+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=ps\u22c5z^$ is then governed by

where the case $s=0$ corresponds to the slow geostrophic mode with amplitude $g0$, and $s=\xb1$ to the (two) fast gravity wave mode with $g+$ and $g\u2212$. For $Ro=0$, that is, in the linear case, Eq. (10) describes the fast oscillation in time of the gravity waves $g\xb1=ei\omega \xb1t$, 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\u2061(g0,g\xb1)$.

### 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=\u2212u\u22c5\u2207u$ and $np=\u2212\u2207\u22c5pu$, the nonlinear terms become

and the integral $Is$

for the modes $s=0,\xb1$. 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

with

### e. The balanced models

The spectral representation of the primitive equations given by Eq. (10) describes the evolution of the gravity wave amplitudes $g\xb1$ 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\xb1$. 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\xb1=O\u2061(Ro)$ and expand the gravity wave amplitudes accordingly as

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

Since the wave mode $g\xb1$ 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

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

Together with the diagnostic expression for the auxiliary wave amplitude $f1\xb1$, which is derived below, Eq. (19) is a second-order balanced model. Since the expression for $f1\xb1$ 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

and so on for higher order. Again, the expression for $f2\xb1$ 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\xb1$ 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\xb1$, as the nonbalanced or fast gravity wave mode.

### f. The ageostrophic balanced modes

The auxiliary wave modes $fn\xb1$ 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\xb1$, 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\xb1$ need to be derived. To first order in Ro, Eq. (10) becomes for the gravity wave modes $s=\xb1$

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

This balancing of $f1\xb1$ yields $\u2202t*f1\xb1=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\xb1$ corresponds to the first iteration step in the balanced initialization technique by Machenhauer (1977), and that $f1\xb1$ also corresponds to the ageostrophic variables in the (first order) quasigeostrophic approximation.

Setting $\u2202t*fn\xb1=0$ to suppress wave activity in general, Eq. (10) for $s=\xb1$ becomes

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

From the first line of Eq. (24), $f2\xb1$ is calculated, from the second line $f3\xb1$, etc. Since only slow time derivatives show up, the ageostrophic balanced modes $fn\xb1$ are only slowly evolving in time as the geostrophic mode. The combination of geostrophic mode amplitude $g0$ and $fn\xb1$ given by $g0q0+\u2211iRoifi\xb1q\xb1$ 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\xb1$ from a numerical model to find $fn\xb1$ up to fourth order. The term $K\xb1\u2061(g0,fn\xb1)$ is given by Eq. (13) replacing $g\xb1$ with $fn\xb1$ and found by evaluating the integrals $I\xb1$ using the numerical model. To find $\u2202Tf1\xb1$, Eq. (22) is analytically differentiated with respect to *T* as in Kafiabad and Bartello (2018), which yields

using the relation $I\xb1\u2061(g0+\u2202Tg0,0)=I\xb1\u2061(g0,0)+\u2202TI\xb1\u2061(g0,0)+I\xb1\u2061(\u2202Tg0,0)$, which derives from the definition of the integral $Is$, and using $\u2202Tg0=\u2212iI0\u2061(g0,0)$. To find $\u2202Tf2\xb1$ and $\u2202Tf3\xb1$, the corresponding expressions get too complicated, and $f2\xb1$ and $f3\xb1$ 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 $\u2212A\u22074u$ and $\u2212A\u22074b$, with a viscosity (or diffusivity) $A~O\u2061(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 $\u2212A\u22074p$ shows up—which enters the matrix $A$ on the diagonal in the spectral representation—but also a damping term $A\u22074p|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

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 $\omega s$ obtain an imaginary part of $O\u2061(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

Since the magnitude of $I\u2061(\omega \xb1)=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\xb1$ for large **k**, but using the complex $\omega \xb1$.

## 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\xb1\u2061(g0,0)$ after forward transform and projection on $p\xb1$, and thus finally $f1\xb1$ from Eq. (27). Corresponding steps yield the other ageostrophic balanced modes $fi\xb1$. The inverse (vertical and horizontal) transform of $g0q0+\u2211i=1KRoifi\xb1q\xb1$ 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\u2061(y)=\u2061(2\u2061\u2009exp\u2061[\u2212\u2061(y\u22121.25)2/0.12\u2212\u2061(y\u22123.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=\pi $. 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=\u2061(Ro/2)\Delta 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|\u2202x\upsilon \u2212\u2202yu|$ 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.

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 $t\u224810$, 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 $u\u2212u0$. We call $u0$ and $ui$ the zonal velocity corresponding to the inverse transform of $g0q0$ and $Roifi\xb1q\xb1$, 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 $u\u2212u0$ 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).

Only in the fourth order (Figs. 1g,h) do small differences between $u4$ and $u\u2212u0\u2212u1\u2212u2\u2212u3$ 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.

Figure 3a shows the residual wave energy $EK=|u\u2212\u2211n=0Kun|2/2+\u2061(p\u2212\u2211n=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 $~\u2009exp\u2061(\u22122/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.

### 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 $\u2202zb\u2264\u22121$ that arise at higher amplitudes of *u* and show up locally in the simulations. The scaled first baroclinic Rossby radius is $Lr=c1=1/\pi $, 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 $\u2202Tfn\xb1$ 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\xb1$. Furthermore, it is necessary to exchange Eq. (27) with Eq. (22) and Eq. (24) but with complex $\omega \xb1$ 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|\u2202x\upsilon \u2212\u2202yu|$ 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.

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 $w\u2212w1$, which is the second-order residual signal. This similarity between $w2$ and $w\u2212w1$ indicates that most of the residual signal seen in $w\u2212w1$ 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 $w\u2212w1\u2212w2$ 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$.

The apparent wave signal at fourth order (i.e., in $w\u2212w1\u2212w2\u2212w3$ 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 remove^{2}, 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 $\u2202zb\u2264\u22121$ 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 $\u2202zb<\u22121$ 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*, $w\u2212w1$, $w\u2212w1\u2212w2$) 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) $\u2202zb\u2264\u22121$ 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 $\u2202zb~\u22121$—or where potential vorticity becomes negative as criterion for symmetric instability (Olbers et al. 2012)—but not by spontaneous wave emission.

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

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 $\delta i+p=\u2061(pi+1\u2212pi)/\Delta x=\delta i+1\u2212p$, where $\Delta x$ is the grid spacing in *x* direction, and averaging operators $pi\xaf+=(pi+pi+1)/2=pi+1\xaf\u2212$ and similar for the other directions with grid spacing $\Delta y$ and $\Delta z$ are introduced. The vertical grid has grid points $k=1,\u2026,N$ with pressure and horizontal velocity at depth level $zk=\Delta z\u2061(N\u2212k\u22121/2)$, vertical velocity at depth level $zk+\Delta z/2$, and surface pressure $ps$ at $zN+\Delta z/2=0$. The surface boundary condition is then given by

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

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 $\delta k+p=b\xafk+$, the surface boundary condition, and the integrated continuity equation $wk=\u2212\Delta z\u2211k\u2033=1k\u2207\u22c5uk\u2033$ yields

with the discrete vertical mode operator

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

Both the integral form and the solution to the differential form yield identical eigenfunctions and eigenvalues for the discrete grid. Expressing $ui,j,k=\u2211nui,jn\varphi kn$ and $pi,j,k=\u2211npi,jn\varphi 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=\u2211n\u2061(cn2+\Delta z2/4)\u2207\u22c5un\delta k+\varphi n$ and $bi,j,k=\u2211npn\u2061[\delta k\u2212\varphi n\u2212\Delta z/\u2061(4cn2)\u2061(\varphi kn+\varphi k\u22121n)]$.

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

with time levels *n*, the time step $\Delta t$, and with the interpolation weights $\alpha =1.5+\u03f5$, $\beta =\u2212\u2061(0.5+\u03f5)$. The parameter $\u03f5>0$ allows for a stable integration if $\Delta t$ is chosen appropriately small. The value of $\Delta t$ can be estimated from the eigenvalues of the linear discrete system $\omega \xafs$, which differ in general from the physical eigenvalues $\omega s$ and are given by

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

## REFERENCES

*Waves in Fluids*. Cambridge University Press, 516 pp.

*Ocean Dynamics*. Springer, 703 pp.

## 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\xac$ denotes anticlockwise rotation of the vector $u$ by $\pi /2$, that is, $u\xac=\u2061(\u2212\upsilon ,u)$ for $u=\u2061(u,\upsilon )$.

^{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.