## 1. Introduction

The linearized rotating Boussinesq equations admit two solution types – inertia–gravity waves and geostrophic (vortex) motions. The wave and geostrophic solutions form the foundation for how we understand and interpret ocean and atmosphere observations in a wide variety of contexts. These two types of linear motion are predictive in the sense that knowledge of the solution at one time enables knowledge of the solution for all time. Such solutions thus guide our intuition for how the ocean evolves, while deviations from those predictions also serve as a direct measurement of nonlinearity. For these reasons, a great deal of effort goes into separating fluid motions into these two types of solutions.

Wave and vortex solutions can be separated in the frequency–wavenumber domain by utilizing the dispersion relation of the linear wave solutions, a method well suited to model output (e.g. Savage *et al.* Reference Savage2017; Torres *et al.* Reference Torres, Klein, Menemenlis, Qiu, Su, Wang, Chen and Fu2018). With more sparse *in situ* observations, other methods have been developed to make this same separation; however, these typically require additional statistical assumptions to overcome the limitations of sparse sampling (Lien & Müller Reference Lien and Müller1992; Bühler, Callies & Ferrari Reference Bühler, Callies and Ferrari2014; Lien & Sanford Reference Lien and Sanford2019; Oscroft *et al.* Reference Oscroft, Sykulski and Early2021). In idealized Boussinesq models with triply periodic domains and constant stratification, such a decomposition can be made unambiguously at each instant in time (Bartello Reference Bartello1995; Smith & Waleffe Reference Smith and Waleffe2002; Waite & Bartello Reference Waite and Bartello2006). For each resolved wavenumber the decomposition splits the flow into two inertia–gravity waves ($A_\pm$), with frequencies that lie between the Coriolis, $f_0$, and the buoyancy frequency, $N$; and a zero-frequency, geostrophic solution ($A_0$), which accounts for all linear potential vorticity (PV). Thus, the wave-vortex decomposition is a linear transformation that projects the variables $(u,v,\rho )$ onto an equivalent representation $(A_+, A_-, A_0)$ of two wave and one vortex mode without loss of information. Note, because the transformation uses vertical eigenmodes that guarantee the continuity equation is satisfied, the vertical velocity, $w$, is redundant and not needed in the transformation. Furthermore, the inverse transformation can recover both $w$ and pressure from the wave-vortex components $(A_+, A_-, A_0)$.

Aside from being an alternative and compact representation of the dynamical variables, the wave-vortex decomposition has a number of applications. Unlike other spectral representations of fluid flow, the wave-vortex projection is useful because it projects directly onto solutions of the equations of motion, including, as shown in this manuscript, in the case of arbitrary stratification. Each wave and vortex solution is energetically independent and coefficients of the projection are thus physically meaningful, directly encoding the amplitude and phase of the unique wave and geostrophic solutions. For example, applying the wave-vortex projection to output from a perfectly linear wave model will show no changes in amplitude and phase over time, while in contrast, a nonlinear wavelike process will have amplitude and phase that become decorrelated with time. Alternatively, for flows that bear no resemblance to the linear solutions, the projection may not be meaningful – thus the interpretation of the components as representations of wave and geostrophic components is ultimately problem specific.

One of the simplest diagnostics utilizing the wave-vortex decomposition is assessing how total energy shifts between inertial, internal gravity waves and geostrophic solutions. The physical mechanisms that transfer energy between wave-vortex modes can further be diagnosed by projecting the nonlinear equations of motion into wave-vortex space. This approach was used by Lelong & Riley (Reference Lelong and Riley1991) and Bartello (Reference Bartello1995) to diagnose transfer in the turbulence cascade, and also, for example, by Arbic *et al.* (Reference Arbic, Scott, Flierl, Morten, Richman and Shriver2012) to diagnose energy transfers across modes and wavenumbers in quasi-geostrophic turbulence. With the nonlinear equations of motion projected onto the wave-vortex modes, it is also possible to create a series of reduced-interaction models, as has been done for triply periodic domains (Remmel Reference Remmel2010; Remmel, Smith & Sukhatme Reference Remmel, Smith and Sukhatme2010; Hernandez-Dueñas, Smith & Stechmann Reference Hernandez-Dueñas, Smith and Stechmann2014). These models are reduced versions of the equations of motions that restrict interactions between certain modes. For example, restricting interactions between only PV modes results in the quasi-geostrophic equations, while restricting interactions between only wave modes results in an extension of the weak wave turbulence model (Remmel Reference Remmel2010).

While the aforementioned studies have addressed and utilized the wave-vortex decomposition for the case of constant stratification, typical stratification profiles in the ocean often resemble an exponential-like function, as shown in figure 1. A computational challenge that arises when solving the equations of motion on a regular grid with such stratification is that, while the grid resolution (black dots in the figure) is more than adequate at depth, rapid variations near the surface are not resolved, even with large numbers of grid points (257 in this case). To address this limitation, in this manuscript we extend the wave-vortex decomposition for Boussinesq flows to arbitrary non-constant stratification. This involves solving an eigenvalue problem (EVP) to obtain the vertical dependence (Early, Lelong & Smith Reference Early, Lelong and Smith2020), rather than Fourier mode expansions in the three spatial directions as in the case of constant stratification. We begin in §§ 2 and 3 with the linearized equations of motion and their solutions. Section 4 then details the projection onto the vertical modes, while § 5 shows the decomposition itself.

In § 6 we provide an example application in which we diagnose the results of the Cyprus eddy studied by Lelong, Cuypers & Bouruet-Aubertot (Reference Lelong, Cuypers and Bouruet-Aubertot2020), in which an eddy in geostrophic balance superimposed with inertial oscillations at the surface transfers energy from the inertial oscillations to generate internal gravity waves. The present decomposition is performed at each instant in time using the same model output, and is shown to agree with the temporal filtering based method used in the original study. Results of the present analysis show that advection of geostrophic vorticity by the inertial oscillations accounts for all the energy transfer from the inertial oscillations to internal gravity waves, confirming the hypothesis in the original study.

Section 7 discusses the implications of these results and addresses general challenges encountered in numerical modelling, including the important implications that, in cases of variable stratification, regularly spaced grids may only resolve a small fraction of the physically relevant vertical modes, and as such, it may be more computationally efficient to integrate the nonlinear equations of motion in wave-vortex space than in physical space. Finally, § 8 offers some concluding remarks. Appendix A shows how the results are simplified for constant stratification, and appendix B details the numerical implementation. The projection of the nonlinear equations of motion onto the wave-vortex modes is documented in appendix C.

## 2. Background

The linearized, unforced, inviscid equations of motion for fluid velocity $u(x,y,z,t)$, $v(x,y,z,t)$, $w(x,y,z,t)$, on an $f$-plane are

*a*)\begin{gather} \partial_t u - \,f_0 v ={-} \frac{1}{\rho_0} \partial_x p \end{gather}

*b*)\begin{gather}\partial_t v + \,f_0 u ={-} \frac{1}{\rho_0} \partial_y p \end{gather}

*c*)\begin{gather}\partial_t w ={-}\frac{1}{\rho_0} \partial_z p - g \frac{\rho}{\rho_0} \end{gather}

*d*)\begin{gather}\partial_x u + \partial_y v + \partial_z w = 0 \end{gather}

*e*)\begin{gather}\partial_t \rho + w \partial_z \bar{\rho} = 0. \end{gather}

Here, $p(x,y,z,t)$ and $\rho (x,y,z,t)$ are perturbation pressure and density, respectively, defined such that total pressure $p_{{tot}}(x,y,z,t) = p_0(z) + p(x,y,z,t)$ and total density $\rho _{{tot}}(x,y,z,t) = \rho _0 + \bar {\rho }(z) + \rho (x,y,z,t)$ where $\partial _z p_0(z) = -g( \rho _0 + \bar {\rho }(z) )$. All variables in (2.1*a*)–(2.1*e*) are functions of $x$, $y$, $z$ and $t$, except $\bar {\rho }$, which is only a function of $z$. We use the usual definition of buoyancy frequency, $N^{2}(z) \equiv -({g}/{\rho _0} )\partial _z \bar {\rho }$. Throughout this manuscript we use the linear approximation to isopycnal displacement $\eta \equiv -\rho /\bar {\rho }_z$ rather than density anomaly. With this notation (2.1*e*) becomes $w=\partial _t \eta$ and (2.1*c*) can be similarly rewritten.

We assume boundaries are periodic in the horizontal, $(x,y)$, and bounded in the vertical, $z$. The lower boundary is assumed flat at $z=-D$ with free slip and $w(-D)=0$, and no density anomaly, $\rho (-D)=0$. Similarly, the upper boundary is taken to be a free-slip rigid lid with $w(0)=0$ and also no density anomaly, $\rho (0)=0$.

The depth-integrated energy densities of the flow are

*a*–

*c*)\begin{equation} {{\rm HKE}}\equiv \tfrac{1}{2}\int_{{-}D}^{0} (u^{2}+v^{2})\, \textrm{d}z,\quad {{\rm VKE}}\equiv \tfrac{1}{2}\int_{{-}D}^{0} w^{2}\, \textrm{d}z, \quad {{\rm PE}} \equiv \tfrac{1}{2} \int_{{-}D}^{0} N^{2} \eta^{2}\, \textrm{d}z, \end{equation}

where ${{\rm HKE}}$, ${{\rm VKE}}$ and ${{\rm PE}}$ are the horizontal kinetic energy, vertical kinetic energy and potential energy per unit mass, respectively. The other conserved quantity of interest is the quantity typically identified as quasi-geostrophic PV,

which can be directly derived from the linear equations (2.1*a*)–(2.1*e*), or found as the linear approximation to the available potential vorticity (APV) as defined by Wagner & Young (Reference Wagner and Young2015).

It is noteworthy that linearized Ertel PV does not correspond to a useful quantity in this model – it is neither conserved nor time independent for the internal gravity wave solutions (see (3.23)). Linearized Ertel PV is

where $\boldsymbol {\zeta }$ is vorticity and $\zeta ^{z}=\partial _x v - \partial _y u$ is its vertical component. Notable is that linear Ertel PV per (2.5) does not equal the conserved PV quantity (2.3) – the primary difficulty being that $\partial _z \eta$ is not proportional to $\partial _z \rho / \bar {\rho }_z$ for non-constant stratification. Applying the total derivative to (2.5) results in

which is not a conservation equation, but a balance between three terms: local changes in the vertical component of vorticity, local changes in the vertical gradient of the density anomaly and the vertical advection of the background density gradient. The connection between (2.6) and the conserved PV (2.3) is found using the thermodynamic equation (2.1*e*) and re-arranging, which reproduces (2.3),

up to a scaling factor. The key difference between quasi-geostrophic PV (2.3) and linear Ertel PV (2.5) is that the latter neglects vertical advection of the background density gradient. In the present context then, APV as defined in Wagner & Young (Reference Wagner and Young2015) is the relevant conserved quantity.

## 3. Wave-vortex solutions

Solutions to (2.1*a*)–(2.1*e*) are assumed to take the separable form

for $u$, $v$, $p$ and

for $w$ and $\eta$. This presumes a Fourier basis satisfying the periodic boundary conditions in $x$, $y$ and real-valued functions $F_{jkl}(z)$, $G_{jkl}(z)$ satisfying the vertical boundary problem. The summation is over all wavenumbers $k$, $l$, but also over eigenmodes $j$ from the bases $\{ F_{jkl}(z) \}$ and $\{ G_{jkl}(z) \}$, indicated with subscripts to emphasize their dependence on wavenumbers $k$ and $l$. The coefficients $\tilde {f}_{jkl}(t)$, $\tilde {g}_{jkl}(t)$ are complex, encoding both amplitude and phase. The ‘c.c.’ refers to the complex conjugate, which contains half the power of the real-valued solutions, but no new information. Although the wave-vortex decomposition is performed at fixed time in the time domain $\tilde {f}_{jkl}(t)$, it is useful to express solutions in the frequency domain, in which case we denote the variables with $\hat {\cdot }$, e.g. $\tilde {f}_{jkl}(t)= \hat {f}_{jkl}\,\textrm {e}^{\textrm {i}\omega t}$. Finally, we will often drop the subscripts ${jkl}$ entirely, and simply work with the coefficients at a fixed $j$, $k$ and $l$.

Using the thermodynamic equation (2.1*e*) to replace $w$ with $\eta _t$, solutions to (2.1*a*)–(2.1*d*) must satisfy

*a*)\begin{gather} \textrm{i} \omega \hat{u} -\, f_0 \hat{v} ={-} \textrm{i} k \frac{\hat{p}}{\rho_0} \end{gather}

*b*)\begin{gather}\textrm{i} \omega \hat{v} +\, f_0 \hat{u} ={-} \textrm{i} l \frac{\hat{p}}{\rho_0} \end{gather}

*c*)\begin{gather}\left(N^{2} - \omega^{2} \right) \hat{\eta} G ={-}\frac{\hat{p}}{\rho_0} \partial_z F \end{gather}

*d*)\begin{gather}\left( \textrm{i} k \hat{u} + \textrm{i} l \hat{v}\right) F ={-}\textrm{i} \omega \hat{\eta} \partial_z G. \end{gather}

We now examine all possible solutions to these equations of motion by considering zero and non-zero frequency, as well as zero and non-zero horizontal wavenumbers. The vertical structure is treated in detail for each solution.

### 3.1. Geostrophic solutions, $\omega = 0, k^{2} + l^{2} = 0$

Geostrophic solutions require a horizontal density anomaly, and because there is no mean vertical or horizontal density anomaly by definition, there are no mean geostrophic currents in this model. Any mean density anomaly should be subsumed into the definition of the mean density $\bar {\rho }$.

### 3.2. Geostrophic solutions, $\omega = 0, k^{2} + l^{2} > 0$

Geostrophic solutions have no time variation, and the thermodynamic equation therefore implies that $w=0$. Assuming non-zero horizontal wavenumber $k^{2} + l^{2} > 0$, the equations of motion (3.3*a*)–(3.3*d*) reduce to

*a*)\begin{gather} -\, f_0 \hat{v} ={-} \textrm{i} k \frac{\hat{p}}{\rho_0} , \end{gather}

*b*)\begin{gather}f_0 \hat{u} ={-} \textrm{i} l \frac{\hat{p}}{\rho_0} , \end{gather}

*c*)\begin{gather}N^{2} \hat{\eta} G_g ={-} \frac{\hat{p}}{\rho_0} \partial_z F_g , \end{gather}

*d*)\begin{gather}\left( \textrm{i} k \hat{u} + \textrm{i} l \hat{v}\right) F_g = 0 , \end{gather}

where $F_g(z)$, $G_g(z)$ denote the geostrophic vertical structure functions. The only equation of consequence for the vertical structure is (3.4*c*). With no vertical velocity, the rigid lid boundary conditions place no constraint on $F_g(z)$, $G_g(z)$. The decision to disallow density anomalies at the boundaries implies that $G_g(z)$ is an odd function, and therefore $F_g(z)$ is an even function. Although gravity $g$ does not enter into (3.4) without a free surface, it is still convenient to set the separation constant in (3.4*c*) such that $\hat {p}=\rho _0 g \hat {\eta }$ and $N^{2} G_g = -g \partial _z F_g$, with $G_g(0)=G_g(-D)=0$ at the boundaries. This allows the amplitude of the solution to be expressed in terms of sea-surface height, analogous to typical notation for geostrophic motions.

The geostrophic solution, or vortex solution, is given by,

where $\hat {A}_0$ is a complex-valued amplitude containing the phase information, and $\theta _0=k x + l y$.

As a consequence of only having one constraint connecting $F_g(z)$ and $G_g(z)$, there is no preferred set of vertical basis functions for the geostrophic solution. Any complete basis can be used to represent the geostrophic solution. However, near-geostrophic theories with a different choice of scalings, such as quasi-geostrophy (QG; e.g. see Pedlosky (Reference Pedlosky1987)), have non-zero vertical velocities and therefore still require that three-dimensional continuity be satisfied. To maintain continuity we take (3.3*d*) and set the separation constant to $h$, such that $F(z) = h \partial _z G(z)$ for all $z$. This additional requirement, combined with the hydrostatic vertical momentum condition $N^{2} G = -g \partial _z F$, results in two Sturm–Liouville eigenvalue problems for hydrostatic (HS) vertical modes,

with boundary conditions $G^{HS}(0)=G^{HS}(-D)=0$ or,

with $\partial _z F^{HS}(0)=\partial _z F^{HS}(-D)=0$ where $j$ is the mode number and eigenvalue $h_j$ is the equivalent depth. It follows directly from Sturm–Liouville theory that the vertical modes resulting from the HS EVPs satisfy the orthogonality conditions

and

where we have implicitly normalized the amplitude of the modes. The ${1}/{g}$ normalization in (3.8) arises naturally when using a free-surface boundary condition, and is kept here for consistency.

The importance of the $\{ G^{HS}_j(z) \}$ and $\{ F^{HS}_j(z) \}$ bases are twofold. First, Sturm–Liouville theory guarantees that they are complete, and therefore capable of representing any function. Second, the specific relationship between these modes is such that both continuity and the linearized vertical momentum equation are satisfied. In practice, this means that they often reflect the vertical structure of various linear solutions. It is in this sense that $\{ G^{HS}_j(z) \}$ and $\{ F^{HS}_j(z) \}$ are ‘preferred’ bases for representing certain flows, including quasi-geostrophy and hydrostatic linear internal waves.

The horizontal kinetic energy and potential energy of the geostrophic solution (3.5) as a function of depth are found by averaging over time and horizontally, including the energy from the complex conjugate,

where $K^{2}=k^{2}+l^{2}$. Vertical kinetic energy is identically zero. If we use the hydrostatic normal modes $F^{HS}_j$, $G^{HS}_j$ then depth-integrated horizontal kinetic energy reduces to ${{\rm HKE}}_g=({\hat {A}_0^{2}}/{4})({g^{2} h_j}/{\,f_0^{2}})K^{2}$ and depth-integrated potential energy reduces to ${{\rm PE}}_g = ({\hat {A}_0^{2}}/{4})g$.

The linearized potential vorticity is,

as is traditionally written, or simply

after using the hydrostatic modes and (3.7) to rewrite $F_g$. These expressions are exactly the potential vorticity identified in the quasi-geostrophic potential vorticity equation. In contrast, the Ertel PV is,

which does not correctly account for changes in the density gradient (see also Wagner & Young Reference Wagner and Young2015).

Under rigid lid conditions, there also exists a barotropic mode ($j=0$) where $F^{HS}_0(z)=\textrm {const}$ with no associated buoyancy anomaly, $G^{HS}_0(z)=0$. This case will be handled separately in the decomposition.

### 3.3. Inertial oscillation solution, $\omega \neq 0, k^{2} + l^{2} = 0$

This solution has no vertical velocity, density anomaly or pressure gradients. It is simply a horizontally uniform oscillating horizontal velocity field, with no constraints on vertical structure other than the boundary conditions. In the triply periodic model used in Smith & Waleffe (Reference Smith and Waleffe2002) this solution is referred to as the vertically sheared horizontal mode, while in the bounded domain it is identified as the inertial oscillation solution,

Here, since there is no conjugate to $k^{2}+l^{2}=0$, the amplitude is purely real. $F_I(z)$ is an arbitrary function, and can be expanded in any complete basis. This is noteworthy because it essentially leaves the boundary conditions for $F_I(z)$ unspecified, and unlike other solutions considered here, $\partial _z F_I(0)$ and $\partial _z F_I(-D)$ are not necessarily zero. Therefore, one must be careful not to expand $F_I(z)$ in a basis with unnecessarily restrictive boundary conditions. That said, there is not necessarily any physical insight to be gained from this additional freedom at the boundaries, and it would certainly be reasonable to restrict the model to solutions where $\partial _z F_I(0)=\partial _z F_I(-D)=0$.

### 3.4. Wave solutions, $\omega \neq 0, k^{2} + l^{2} > 0$

Similar to the geostrophic solution where we assumed that $\hat {p}=\rho _0 g \hat {\eta }$, the vertical momentum equation requires that $(N^{2}-\omega ^{2})G = -g \partial _z F$. Combined with continuity $F = h \partial _z G$, the vertical dependence vanishes from the problem and we are left with

This system of equations admits the internal wave solutions when

The $\pm$ wave solutions are given by,

where the horizontal phase is given by $\theta _\pm =k x + l y \pm \omega t + \phi$ and the amplitude is chosen so that depth-integrated total energy is $\hat {A}^{2} h/2$, as will be shown below.

Combining the vertical constraints from non-hydrostatic vertical momentum $(N^{2}-\omega ^{2})G = -g \partial _z F$ and continuity $F = h \partial _z G$ with the dispersion relation (3.17) results in the $K$-constant, non-hydrostatic Sturm–Liouville problem (Early *et al.* Reference Early, Lelong and Smith2020),

The eigendepth $h_j$ and eigenfrequency $\omega _j$ are interchangeable using the dispersion relation (3.17) with fixed $K$. Note that the EVP could have been written in terms of a fixed frequency $\omega$ (with no subscript $j$), with eigendepth $h_j$ and eigenwavenumber $K_j$ (with subscript $j$), but the constant frequency formulation is not relevant for the decomposition problem at fixed time.

The depth-integrated energies for the $j$th internal wave mode at total wavenumber $K$ are,

which sum to a depth-integrated total energy of ${\hat {A}_\pm ^{2} h_j}/{2}$. The internal wave solutions have zero potential vorticity per (2.3), ${{\rm PV}}_\pm =0$; but they do have Ertel PV per (2.5),

again suggesting that Ertel PV may not be the appropriate quantity for this model.

## 4. Orthogonality and projection

The primary challenge that separates this wave-vortex decomposition from previous ones is dealing with the vertical modes resulting from the $K$-constant EVP in (3.19). In a vertically periodic domain with constant stratification in $z$, Fourier series are an appropriate basis. For a vertically bounded domain with arbitrary stratification in $z$ and no buoyancy anomaly at the boundaries, the appropriate basis are the eigenmodes $G_j$ of (3.19) with $G(0)=G(-D)=0$.

### 4.1. Orthogonality

The non-hydrostatic Sturm–Liouville problem given by (3.19) implies that for a given wavenumber $K$, two modes $G_i(z)$, $G_j(z)$ satisfy the orthogonality condition,

where we have normalized the modes. Unlike the hydrostatic case, there does not appear to be an equivalent Sturm–Liouville problem for the non-hydrostatic $F_j$ modes (with constant $K$) and therefore no associated orthogonality condition. The expression

as far as we know, cannot be coerced to Sturm-Liouville form. The closest relationship we are able to find is

The difference between (3.9) and (4.3) is significant – the former can be used on any function, while the latter requires a specific relationship between the dynamical variables to project on the $F_j$ modes.

### 4.2. Projection

A dynamical variable that expands in $G$, such as density anomaly, $\rho (z)$, can be written as in (3.2), e.g.

where the coefficients are recovered with

The projection operation (4.5) first requires taking a Fourier transform of the variable, then invoking the orthogonality condition (4.1) with the $j$th vertical mode $G_{jkl}(z)$ for wavenumber $K=\sqrt {k^{2}+l^{2}}$. However, in order to use orthogonality condition (4.3) as a projection operator, dynamical variables expanded in $F$ must be added to a related dynamical variable that scales like $h G$. For example, the divergence, $\delta =\partial _x u + \partial _y v$, and vertical vorticity, $\zeta =\partial _x v - \partial _y u$, can be recovered from the wave solution (3.18) with,

where $\delta (z,t) = \sum \delta _j(t) F_j(z)$ and $\zeta (z,t) = \sum \zeta _j(t) F_j(z)$. However, this only works for wave solutions since the geostrophic solution does not have the same relationships between $(u,v)$ and $\eta$. It thus appears that (4.3) is not particularly useful in recovering solutions.

To project variables $u$ and $v$ (and also $p$) that are expanded in $F$, we instead use the relationship derived from continuity, $F_j(z)=h_j \partial _z G_j(z)$, and consider the depth-integrated quantities. That is, if

then we compute $U=\int _{-D}^{z} u \, \textrm {d}z^{\prime }$ so that,

which can then be projected using (4.5) to recover $\tilde {u}_{jkl}(t)$. Notable here is that the depth-integrated quantities represented by (4.9) are themselves depth dependent.

As discussed in the next section, the only part of the solution that must be handled in a special manner is the barotropic $j=0$ mode $F_0(z)$, which as previously discussed has no projection on the $G$ modes in the rigid lid case. In practice, the integration linking (4.8) to (4.9) can be performed by projecting $u$ onto either the $\{F^{HS}_j(z) \}$ basis or a cosine basis (either of which satisfy the correct boundary conditions and have a constant/barotropic mode), integrating spectrally, and then transforming back to the spatial domain.

## 5. Wave-vortex decomposition

Per the previous discussion, the wave-vortex decomposition requires integrating $(u,v)$ to get $(U,V)$, taking the Fourier transform in the horizontal of $(U,V,\eta )$ and then projecting the vertical structure at each horizontal wavenumber $k$ and $l$ onto the vertical eigenmodes found via the $K$-constant EVP, (3.19). Early *et al.* (Reference Early, Lelong and Smith2020) developed a methodology for the computation and projection onto these modes. Written as a sum of individual linear solutions, and explicitly including the dependence on $j,k,l$, the three required variables are expressed as

where $\tilde {U}_{ijk}(t)= \tilde {u}_{ijk}(t) h_{ijk}$ and $\tilde {V}_{ijk}(t)= \tilde {v}_{ijk}(t) h_{ijk}$. The horizontal Fourier transform followed by the vertical projection then recovers $\tilde {U}_{ijk}(t)$, $\tilde {V}_{ijk}(t)$ and $\tilde {\eta }_{ijk}(t)$.

### 5.1. Non-zero wavenumber solutions, $k^{2} + l^{2} > 0$, $j=0$

When vertically integrating the horizontal velocities $u$, $v$ to project onto the vertical modes, the amplitude of the $j=0$ mode must be handled separately. The $j=0$ mode for the rigid lid boundary condition has no density anomaly, $\tilde {\eta }(t)=0$, and no divergence, $\tilde {\delta }(t) = \textrm {i} k \tilde {u}(t) + \textrm {i} l \tilde {v}(t)=0$. This leaves only the amplitude and phase of the vorticity $\tilde {\zeta }(t) = \textrm {i} k \tilde {v}(t) - \textrm {i} l \tilde {u}(t)$. The only valid solution is therefore the vortex solution,

valid for all $k^{2} + l^{2} > 0$.

### 5.2. Non-zero wavenumber solutions, $k^{2} + l^{2} > 0$, $j>0$

For each wavenumber $(k,l)$ and mode $j$ there are six unknowns: the amplitudes and phases of the three different solutions. We denote the complex amplitudes as $\hat {A}_+$, $\hat {A}_-$ and $\hat {A}_0$, for the positive and negative wave, and geostrophic solutions, respectively. In matrix form the three linearly independent solutions from (3.5) and (3.18) at wavenumbers $k$, $l$ and mode $j$ are given by

which can be inverted to solve for $\hat {A}_+$, $\hat {A}_-$, and $\hat {A}_0$,

*a*)\begin{gather} \hat{A}_\pm{=} \frac{\textrm{e}^{{\mp} \textrm{i}\omega t}}{2} \left[ \frac{k \omega \pm \textrm{i} l\, f_0}{\omega K h} \tilde{U}(t) + \frac{l \omega \mp \textrm{i} k\, f_0}{\omega Kh}\tilde{V}(t) \mp \frac{gK}{\omega} \tilde{\eta}(t) \right] \end{gather}

*b*)\begin{gather}\hat{A}_0 =\textrm{i} \frac{l\, f_0}{\omega^{2}} \tilde{U}(t) - \textrm{i} \frac{k\, f_0}{\omega^{2}}\tilde{V}(t) + \frac{\, f_0^{2}}{ \omega^{2}} \tilde{\eta}(t). \end{gather}

There is some insight to be gained by defining depth-integrated versions of horizontal divergence and potential vorticity,

*a*,

*b*)\begin{equation} \tilde{\varDelta}(t) = \textrm{i} \left( k \tilde{U}(t) + l \tilde{V}(t) \right), \quad \tilde{\varPi}(t) \equiv \textrm{i} \left( k \tilde{V}(t) - l \tilde{U}(t) \right) -\, f_0 \tilde{\eta}(t). \end{equation}

Now the solution has the form,

*a*)\begin{gather} \hat{A}_+{=} \frac{\textrm{e}^{-\textrm{i} \omega t}}{2 K h} \left[ - \textrm{i}\tilde{\varDelta}(t) - \omega \left( \frac{f}{\omega^{2}} \tilde{\varPi}(t) + \tilde{\eta}(t) \right) \right] \end{gather}

*b*)\begin{gather}\hat{A}_-{=} \frac{\textrm{e}^{\textrm{i} \omega t}}{2Kh} \left[ - \textrm{i} \tilde{\varDelta}(t) + \omega \left( \frac{f}{\omega^{2}} \tilde{\varPi}(t) + \tilde{\eta}(t) \right) \right] \end{gather}

*c*)\begin{gather}\hat{A}_0 ={-} \frac{f}{\omega^{2}} \tilde{\varPi}(t). \end{gather}

Importantly, (5.8*b*)–(5.8*c*) show that the vortex solution is recovered directly from potential vorticity and the sum of the two wave solutions is recovered from the divergence of the transport. Extracting the phase information and energetics of individual wave solutions still requires additional information from vorticity and isopycnal displacement.

### 5.3. Zero wavenumber solutions, $k^{2} + l^{2} = 0$

The only $k^{2} + l^{2} = 0$ solution is still inertial oscillations, per (3.15), with simple rotation and zero isopycnal displacement, i.e.

*a*)\begin{gather} u_I(0) = \left[\tilde{u}(t) \cos(\,f_0 t) - \tilde{v}(t) \sin(\,f_0 t)\right] F_I(z) \end{gather}

*b*)\begin{gather}v_I(0) = \left[\tilde{u}(t) \sin(\,f_0 t) + \tilde{v}(t) \cos(\,f_0 t)\right] F_I(z) \end{gather}

*c*)\begin{gather}\eta_I(0) = 0. \end{gather}

### 5.4. Summary of the decomposition

A key feature of the decomposition is that the recovered coefficients (5.4), (5.8*a*) and (5.9) are strictly independent of time when applied to time-dependent linear solutions. That is, the left-hand sides of these equations are time independent, while the right-hand sides contain terms that are time dependent. This is not a contradiction; it simply reflects the fact that for unforced inviscid flow, the amplitude and phase of the linear solutions will remain fixed for all time. Applying the linear decomposition to nonlinear flows, an important implication of the latter result is that the actual linearity or nonlinearity of a flow can be made precise by assessing time variation in the recovered coefficients. For example, if $\hat {A}_+$ for a given $j,k,l$ at time $t=t_0$ is exactly equal to $\hat {A}_+$ computed at time $t=t_1$, then that component of the flow was perfectly linear in the sense that the wave solution (3.18) exactly described its evolution. The key takeaway when applying the above decomposition is that any time variation in the recovered coefficients, (5.4), (5.8*a*) and (5.9), by definition implies nonlinearity.

## 6. Nonlinear transfers between wave and vortex solutions

Having derived a generalized wave-vortex decomposition for arbitrary stratification, we next demonstrate its utility by analysing output from a nonlinear numerical simulation performed by Lelong *et al.* (Reference Lelong, Cuypers and Bouruet-Aubertot2020) using the Boussinesq model described by Winters & de la Fuente (Reference Winters and de la Fuente2012). The study by Lelong *et al.* (Reference Lelong, Cuypers and Bouruet-Aubertot2020) was motivated by evidence of intense near-inertial wave activity at the base of the quasi-permanent Cyprus eddy in the Eastern Mediterranean (Cuypers *et al.* Reference Cuypers, Bouruet-Aubertot, Marec and Fuda2012), and was designed to explain the origin and dynamics of the observed waves.

Background stratification in the region surrounding the Cyprus eddy changes rapidly near the surface, following an approximately exponential-like profile as shown in figure 1. Within this stratification, the Cyprus eddy was modelled as an axisymmetric vortex in geostrophic equilibrium via the streamfunction,

where $(u_{g},v_{g},\rho _{g}) = (-\partial _y \psi ,\partial _x \psi ,-\rho _0\, f_0 \partial _z \psi /g)$, and the strength and extent of the eddy were set by parameters $A$, $\alpha$, and $\beta$, chosen to closely match the observations. This eddy initial condition projects exactly onto to the geostrophic solutions in § 3, and remains stable in the nonlinear model. To model the effects of an impulsive wind stress at the surface, superimposed on the geostrophic vortex initial condition is an inertial oscillation of the form

which itself projects exactly onto the inertial solution in § 3. In the absence of the eddy, the inertial oscillation also remains stable in a nonlinear $f$-plane model. However, as shown by Lelong *et al.* (Reference Lelong, Cuypers and Bouruet-Aubertot2020), the combined presence of the anticyclonic eddy and the inertial oscillation causes the inertial oscillation to lose energy while generating slightly subinertial internal gravity waves that propagate downward into the eddy core. This classic problem has a rich history, as discussed, for example, in recent papers by Lelong *et al.* (Reference Lelong, Cuypers and Bouruet-Aubertot2020) and Asselin & Young (Reference Asselin and Young2020) and references therein.

The transfer of energy from inertial oscillations to internal gravity waves was estimated using spatial-temporal averaging by Lelong *et al.* (Reference Lelong, Cuypers and Bouruet-Aubertot2020), but can also be computed diagnostically at each instant in time using the wave-vortex decomposition described in the previous sections. Figure 2 shows energy time series for the inertial, geostrophic and wave mode solutions computed via the wave-vortex decomposition, which are consistent with the transfer of energy observed in Lelong *et al.* (Reference Lelong, Cuypers and Bouruet-Aubertot2020) (their figure 12*b*).

In addition to total energies, the present methodology also enables us to examine more precisely which scales are involved in the energy transfers, and further identify exactly the dynamical mechanisms that are responsible. Figure 3 shows the change in energy among the different vertical modes and horizontal scales for geostrophic, inertial and wave components on day 6 of the simulation. The dominant energy transfer for all three components is in the lowest modes. Only the geostrophic flow shows a weak signature of energy transfer in higher modes, although at later times energy transfers occur at higher internal gravity wave modes as well (not shown). As anticipated by Lelong *et al.* (Reference Lelong, Cuypers and Bouruet-Aubertot2020), the peak energy transfer into the waves occurs at horizontal wavelengths similar to the length scales of the geostrophic flow, with corresponding near-inertial frequencies deviating from $f_0$ by only a few per cent. Note, however, that the original study found wave frequencies to be slightly subinertial, an effect caused by the eddy's anticyclonic geostrophic vorticity (Kunze Reference Kunze1985). This shift to subinertial frequencies is not directly captured by the linear decomposition, and instead the subinertial waves alias into other superinertial frequencies.

Lelong *et al.* (Reference Lelong, Cuypers and Bouruet-Aubertot2020) identified the vertical gradient of advection of geostrophic vorticity by inertial oscillations as the most likely dynamical mechanism for transferring energy from inertial oscillations to internal gravity waves. This is consistent with figure 2, which suggests that the waves draw their energy entirely from the inertial flow, with the geostrophic flow acting as a catalyst in facilitating energy transfer. Although the evidence presented in Lelong *et al.* (Reference Lelong, Cuypers and Bouruet-Aubertot2020) is entirely consistent with this hypothesis, the authors were not able to show conclusively whether this mechanism can account for the total energy transferred from inertial oscillations to internal gravity waves.

To compute the energy transfers between inertial, wave and vortex modes, we rewrite the nonlinear equations of motion by projecting them onto wave-vortex space in appendix C. For the problem considered here, the internal wave frequencies do not exceed $3\,f_0$ and we are therefore able to make the hydrostatic approximation, simplifying the mathematics and reducing numerical complexity. Summarizing the results from the appendix, the equations of motion take the form,

where the nonlinear terms are encapsulated in the three terms $F_{\pm 0}$. The transfer of energy is then proportional to $\mathcal {R}(F_{\pm 0} \bar {A}_{\pm 0})$ according to (C17). To confirm that the energy flux term is computed correctly, figure 4 compares the total change in energy of the constituent parts determined via the wave-vortex decomposition at each time step, to the energy flux terms computed from $\mathcal {R}(F_{\pm 0} \bar {A}_{\pm 0})$ also computed at each time step. The two lines are nearly indiscernible, indicating that the wave-vortex projection of the nonlinear equations of motion is correctly reproducing the dynamics of the Boussinesq model.

Appendix C further shows that the nonlinear flux of energy into internal gravity waves via advection of geostrophic vorticity $\zeta _{g}$ by inertial oscillations $(u_{io},v_{io})$ can be written as,

where $\mathcal {F}$ is the projection operator onto the $F_j$ modes and $\mathcal {DFT}$ is the discrete Fourier transform. Figure 4(*a*) shows that this mechanism accounts for all of the transfer of energy into internal gravity waves, directly confirming the hypothesis of Lelong *et al.* (Reference Lelong, Cuypers and Bouruet-Aubertot2020). Additionally, figure 4(*b*) shows that the primary mechanism draining energy from inertial oscillations is the advection of the geostrophic flow by internal gravity waves, with a small contribution from self-advection by internal gravity waves.

Having identified the two dominant physical mechanisms responsible for the nonlinear transfer of energy from inertial oscillations to internal gravity waves, we can now examine exactly which modes and scales are involved in the transfer. Figure 5 shows the nonlinear fluxes from only the two transfer mechanisms identified above, revealing the same dominant modes of transfer as in figure 3. Here, it is clear that the advection of geostrophic flow by internal gravity waves primarily drains energy from the $j=1$ inertial mode, but also shifts some energy to the $j=2$ mode. The advection of geostrophic flow by inertial oscillations transfers energy into the mode $j=1$ internal gravity waves at scale between 50 and 500 km. The results of figure 5 further imply that broader range of modes and scales seen to be exchanging energy in figure 3 must be explained by other mechanisms that transfer energy within the internal gravity wave modes. Indeed, at later times we find that internal gravity wave energy cascades to higher modes (not shown).

## 7. Modelling implications

In addition to the physical insights gained from applying the wave-vortex decomposition, there are also several implications for numerically modelling fluid flows. The first is the rather startling recognition that for the exponential stratification profile in figure 1, 257 evenly spaced vertical grid points in a pseudospectral model only resolves approximately 19 vertical modes. Conversely, this suggests that in cases of challenging stratification, modelling the equations of motion in wave-vortex space, as in appendix C, may actually be more computationally efficient than traditional spectral approaches.

One of the central claims of this manuscript is that the above wave-vortex decomposition can account for all variance of $(u,v,w,\rho )$ at any instant in time. In practice, however, ocean data and numerical models have a finite number of grid points, $N$, and it is not true in general that $N$ vertical modes will be resolved with $N$ grid points. As noted in Early *et al.* (Reference Early, Lelong and Smith2020), only if the grid points are at (or near) the Gaussian quadrature points of the vertical modes, $F(z)$ and $G(z)$, for all resolved $K$, will all variance project onto the modes. For the case of constant stratification, the Gaussian quadrature points are evenly spaced and the vertical modes coincide with cosine and sine bases. In this special case, the $N$ vertical grid points in the Boussinesq model will coincide with $N-2$ resolved internal modes, leaving only the Nyquist frequency unresolved plus a barotropic mode. However, as figure 1 demonstrates, for variable stratification, regions of rapidly changing stratification will lack resolution if an evenly spaced grid is used.

Because of the above issue, many numerical models use alternative vertical coordinates such as $\sigma$ (pressure) coordinates and density coordinates, or other more complicated hybrids, in order to better resolve the solutions. To resolve internal wave modes, Early *et al.* (Reference Early, Lelong and Smith2020) showed that a Wentzel–Kramers–Brillouin (WKB)-scaled coordinate more efficiently positions grid points than a density coordinate when capturing vertical modes. Once the vertical modes are computed, the Gauss quadrature points for the first $N$ modes can computed from the roots of the $N+1$th vertical mode. When creating a numerical model, these points are the optimal choice for resolving the vertical modes.

So what happens when grid points in a numerical model are not able to fully resolve the physics? From a diagnostic point of view, any variance not captured by the $M$ resolvable modes, results in a residual. The residual of projecting a function $f$ onto the $F$ modes is defined as $f_R \equiv f-\mathcal {F}^{-1}[ \mathcal {F}[f] ]$ where $\mathcal {F}$ projects using $M$ modes. For the Cyprus eddy example, the residual energy is shown in green in figure 2, and represents at most $0.3\,\%$ of the total energy, a $48\,\%$ increase from its initial value. Thus the initial conditions start with variance unresolved by the modal decomposition, and nonlinear processes further shift some of the resolved variance into unresolved variance over the course of the simulation. Because the evenly spaced $z$ grid in the model fully resolves the higher modes at lower depths, but only 19 modes near the surface, we expect an accumulation of residual energy at depth. Indeed, we find two peaks in residual energy between 450 and 500 m depth on day 37 of the simulation. How exactly this affects the resulting physics is not entirely clear. All we can say is that, given higher resolution, this residual energy would be associated with higher-mode internal waves which are presently not being represented correctly.

In the case of ocean observations with limited vertical sampling resolution, the issue is completely different than with a numerical model. In this case the physics is certainly evolving correctly, but the unresolved wave modes alias into the lower modes. Depending on the spectral slope of the unresolved modes, this aliasing may or may not have a significant impact on the coefficients of the resolved modes (Early *et al.* Reference Early, Lelong and Smith2020).

A second but related implication for modelling is that it may be advantageous to model the nonlinear equations of motions directly in wave-vortex space. This has the advantage of establishing which wave and vortex solutions are valid *a priori*, using the methods discussed above. Additionally, damping and/or small scale variance removal is then performed directly on the wave and vortex coefficients, which correspond to wave energy and potential enstrophy damping.

In light of the effective vertical resolution issues discussed above, the computational efficiency of modelling the equations of motion in wave-vortex space is relatively favourable for cases of nonlinear stratification. The rate limiting steps are the horizontal and vertical transformations required to compute the terms in (C18*a*). The two basic numerical operations to be performed on a vector of length $N$ are a matrix multiplication, which requires $2N^{2}$ operations, and a fast Fourier transformation (FFT), which requires $\frac {5}{2} \log _2 N - 3 N$ operations when transforming real variables (Canuto *et al.* Reference Canuto, Hussaini, Quarteroni and Zang2006). The vertical transformation must be computed as a matrix multiplication applied to each of the $N_x N_y/2$ wavenumber vectors of length $N_z$, for total computational cost of $N_z^{2} N_x N_y$. The horizontal transformation can be computed using an FFT algorithm applied to each depth $N_z$ for a total computational cost of $\frac {5}{2} N_z N_x N_y \log _2 N_x N_y - 3 N_z N_x N_y$. The transformation from wave-vortex space to physical space requires applying the vertical transformation 10 times (7 for the hydrostatic case) and the horizontal transformation 10 times. To finish the pseudospectral multiplication, the results must be projected back onto wave-vortex space for a grand total of 13 horizontal transformations and 10 vertical transformations. Assuming that $N_x = N_y$, the total computational cost of the horizontal and vertical transforms are approximately equal when $10 \log _2 N_x = N_z$, or $13 \log _2 N_x = N_z$ for the hydrostatic case. This means that for a horizontal resolution of $256^{2}$ the horizontal transformations dominate the computation time until approximately 80–100 vertical modes are used.

## 8. Conclusion

The wave-vortex decomposition presented in this paper unambiguously separates linear wave and geostrophic motions under arbitrary stratification into decoupled modes at any given instant in time. The present decomposition has been fully implemented for arbitrary stratification, as well as the special case of constant stratification (see appendices A and B). The methodology has been validated against output from a linear simulation of a Boussinesq model by confirming that the initial conditions can be exactly recovered at all output times. We have further shown that this method successfully reproduces the results of more traditional methods that rely on spatial-temporal filtering.

In addition to these basic validations, the hydrostatic nonlinear equations of motion projected in wave-vortex space (see appendix C) are shown to successfully reproduce changes in wave-vortex amplitude from a nonlinear Boussinesq model. This suggests that the nonlinear equations of motion can be integrated in wave-vortex space with little modification from the work presented here. Estimates of the computational complexity of the method, presented in § 7, show that numerical integration of the wave-vortex modes may actually perform better than integration in physical space when the stratification profile varies strongly with depth.

One of the more useful aspects of the decomposition is the ability to determine the exact nonlinear pathways that move energy between the wave and geostrophic solutions at different scales. This can be done diagnostically (as in § 6) or while directly integrating the equations of motion. By selectively turning off transfer mechanisms, one can also derive reduced-interaction models, such as in Hernandez-Dueñas *et al.* (Reference Hernandez-Dueñas, Smith and Stechmann2014), but now also for vertically bounded flows with variable stratification.

While the present work is a step towards separating wave and vortex motions in generalized flows, there are still limitations to this methodology that prevent its application to, for example, output from global circulation models. In particular, the flows considered here currently lack (i) surface and bottom buoyancy anomalies (ii) a free surface, (iii) horizontal dependence on stratification or background mean flow and (iv) bottom topography. Two recent studies offer significant progress towards addressing the first two issues. First, Smith & Vanneste (Reference Smith and Vanneste2013) showed that geostrophic motions can be decomposed into uncoupled structures that include a surface buoyancy anomaly, and are also orthogonal (or ‘diagonalize energy’ to use the terminology therein), albeit under rigid lid conditions. Second, Kelly (Reference Kelly2016) showed that in the presence of a free surface, there exists a complete set of modes as well as an orthogonality condition that decouples the surface mode from the interior modes, allowing for an unambiguous partitioning of wave energy. Taken together, these results suggest that it should be possible to create a complete framework that includes both a free surface and a surface buoyancy anomaly.

## Supplementary movie

Supplementary movie is available at https://doi.org/10.1017/jfm.2020.995.

## Acknowledgements

We thank P. Bartello and 2 anonymous reviewers for helpful suggestions that improved the manuscript.

## Funding

We would like to acknowledge National Science Foundation awards OCE-1658564, OCE-1536747 and OCE-1536439 for funding this research.

## Declaration of interests

The authors report no conflict of interest.

## Appendix A. Constant stratification

Assuming constant stratification, $N^{2}(z) = N_0^{2}$, significantly simplifies the problem because the $F$-modes become orthogonal. The vertical modes take the form,

*a*)\begin{gather} G^{N_0}_j(z) = A \sin \left( m_j ( z + D) \right) \end{gather}

*b*)\begin{gather}F^{N_0}_j(z) = A h_j m_j \cos \left( m_j ( z + D) \right), \end{gather}

with eigendepth $h_j=({1}/{g})(({N_0^{2} -\, f_0^{2}})/({k^{2}+l^{2}+m_j^{2}}))$ and vertical wavenumber $m_j = {j{\rm \pi} }/{D}$. Using the normalization $A^{2} = ({1}/{D})({2g}/({N^{2} - f^{2}}))$ results in the following orthogonality conditions,

*a*)\begin{gather} \frac{1}{g}\int_{{-}D}^{0} (N_0^{2} -\, f_0^{2}) G^{N_0}_i(z) G^{N_0}_j(z) \, \textrm{d}z = \delta_{ij} \end{gather}

*b*)\begin{gather}\int_{{-}D}^{0} F^{N_0}_i(z) F^{N_0}_j(z) \, \textrm{d}z = \frac{g h_j^{2} m_j^{2}}{N_0^{2} -\, f_0^{2}} \delta_{ij}. \end{gather}

The orthogonality conditions imply that if $\eta (z) = \sum _n \hat {\eta }_n G^{N_0}(z)$ or $u(z) = \sum _n \hat {u}_n F^{N_0}(z)$, then

*a*)\begin{gather} \hat{\eta}_n = \frac{N_0^{2} -\, f_0^{2}}{g} \int_{{-}D}^{0} G^{N_0}_n(z) \eta(z) \, \textrm{d}z \end{gather}

*b*)\begin{gather}\hat{u}_n = \frac{N_0^{2} -\, f_0^{2}}{g h_n^{2} m_n^{2}} \int_{{-}D}^{0} F^{N_0}_n(z) u(z) \, \textrm{d}z. \end{gather}

The consequence is that $(\hat {u},\hat {v})$ can be recovered directly, without integrating to get transport quantities, that is (5.8*b*)–(5.8*c*) with $\tilde {\varDelta }(t)$ replaced by $\tilde {\delta }(t)$, $\tilde {\varPi }(t)$ replaced by $\widetilde {{\rm PV}}(t)$, and $\widetilde {A_\pm }$ no longer normalized by $h$.

where,

*a*,

*b*)\begin{equation} \tilde{\delta}(t) = \textrm{i} \left( k \tilde{u}(t) + l \tilde{v}(t) \right), \quad \widetilde{PV}(t) \equiv \textrm{i} \left( k \tilde{v}(t) - l \tilde{u}(t) \right) -\, f_0 \tilde{\eta}(t)/h. \end{equation}

## Appendix B. Numerical implementation

The decomposition was tested with output from a linear simulation with a rotating spectral Boussinesq model (Winters & de la Fuente Reference Winters and de la Fuente2012) with constant stratification. Implementation of the methodology requires the following steps,

(i) Discrete Fourier transforms of $u$, $v$ and $\eta$ in $x$ and $y$.

(ii) A discrete cosine transform of $u$ and $v$ and a discrete sine transform of $\eta$ in $z$.

(iii) Computation of the wave-vortex coefficients from the transformed variables.

The last step requires careful bookkeeping to ensure that all terms are properly accounted for and not double counted.

The domain is assumed to have $[N_x\times N_y \times (N_z+1)]$ points, where $N_x$ and $N_y$ are typically in powers of 2 to take advantage of the FFT, while $N_z+1$ has $2^{n}+1$ points to accommodate the type-I discrete cosine transforms (DCT-I) and type-I discrete sine transforms (DST-I) used by the Winters model.

#### B.1. Horizontal transformation

The finite-length Fourier transform in a periodic domain is given by

where $k_j\!=\!{2{\rm \pi} j}/{L}$. For a discretized domain with points at $x_n\!=\!n \varDelta$ where $n\!=\![0 \, \ldots \, N-1]$ and $\varDelta =L/N$, the discrete Fourier transform is

with wavenumbers $k_j={2{\rm \pi} j}/{L}$ now limited to $j=[0 \,\ldots \,N-1]$. Variance is preserved following Plancherel's theorem,

where $\textrm {d}k={2{\rm \pi} }/{L}$ and $S(k_j)={L}/{2{\rm \pi} } |\,\hat {f}(k_j)|^{2}$ is defined as the spectrum.

Applying the $\mathcal {DFT}$ in both $x$ and $y$ using the usual numerical algorithms on a real value function results in a two-dimensional transformed matrix as shown in table 1. For a real-valued function the power is split between the two conjugate pairs, and therefore $\hat {f}(k_j)$ has to be doubled to be compared to $\tilde {f}_{jkl}$ in (3.1). The grey and pink regions in table 1 are Hermitian conjugates of other values in the table. The Nyquist frequency $j=N/2$ is unresolved since the sine at the Nyquist is zero, and thus the orange regions are also ignored. Only the white regions and the cyan component at $k=l=0$ contain the information for the inversion.

#### B.2. Vertical transformation

The finite-length sine and cosine transforms are given by

and

where $m_j = {j {\rm \pi}}/{D}$. The discretized versions of these transforms, the DST-I and DCT-I used by the Winters model, are defined with points at $z_n=n\varDelta$ where $n=0..N_z$ and $\varDelta =D/N_z$. Note that this differs from the discretization for the $\mathcal {DFT}$ by including endpoints. This choice of discretization results in the discrete transforms

and

with vertical wavenumbers at $m_j = {j {\rm \pi}}/{D}$ where $j=0..N_z$. The sum in the $\mathcal {DCT}$ treats the end points separately, as they have only half the width of the other points, $\varDelta /2$. For the $\mathcal {DST}$ the function is zero at the endpoints, $g(z_0)=g(z_n)=0$, and the $m_0$ and $m_{N_z}$ wavenumber components have zero power.

With these definitions of the transform, Plancherel's theorem states that,

with $S(m_j)={D}/{2{\rm \pi} } |\hat {g}(m_j)|^{2}$ and

with $S(m_j)={D}/{2{\rm \pi} } |\hat {f}(m_j)|^{2}$ where $\textrm {d}m={{\rm \pi} }/{D}$. The variance of the $m_0$ wavenumber is notably a factor of 2 larger than the variance of a constant function.

#### B.3. Wave-vortex coefficients

Applying the discrete transformations exactly as defined above to $u$, $v$ and $\eta$ results in the following matrices,

*a*)\begin{gather} \hat{u}_{klj} = \mathcal{DCT}_z \left[\mathcal{DFT}_y \left[\mathcal{DFT}_x \left[ u(x,y,z) \right] \right] \right] \end{gather}

*b*)\begin{gather}\hat{v}_{klj} = \mathcal{DCT}_z \left[\mathcal{DFT}_y \left[\mathcal{DFT}_x \left[ v(x,y,z) \right] \right] \right] \end{gather}

*c*)\begin{gather}\hat{\eta}_{klj} = \mathcal{DST}_z \left[\mathcal{DFT}_y \left[\mathcal{DFT}_x \left[ \eta(x,y,z) \right] \right] \right]. \end{gather}

##### B.3.1 Coefficients, $k^{2} + l^{2} > 0$, $j=0$

Starting with the $j=0$ mode, the discrete transforms give non-zero coefficient functions for $\hat {u}_{kl0}$ and $\hat {v}_{kl0}$, but zero for $\hat {\eta }_{kl0}$. The $\mathcal {DCT}$ inflates the power by a factor of two, but the $\mathcal {DFT}$ returns only half the power of the real-valued function. The result is that,

exactly as written before.

##### B.3.2 Coefficients, $k^{2} + l^{2} > 0$, $j>0$

The projection operations as defined in (A3) can be related to the sine and cosine transformations in (B4) and (B5) by a scaling,

*a*)\begin{gather} \bar{u}_{klj} = \frac{\sqrt{D(N_0^{2} -\, f_0^{2})}}{h_{klj} m_{j} \sqrt{2g} } \hat{u}_{klj}, \end{gather}

*b*)\begin{gather}\bar{v}_{klj} = \frac{\sqrt{D(N_0^{2} -\, f_0^{2})}}{h_{klj} m_{j} \sqrt{2g}} \hat{v}_{klj}, \end{gather}

*c*)\begin{gather}\bar{\eta}_{klj} = \frac{\sqrt{D(N_0^{2} -\, f_0^{2})}}{\sqrt{2g}} \hat{\eta}_{klj}. \end{gather}

The wave-vortex coefficients are then recovered with,

## Appendix C. Nonlinear equations of motion

Here, we rewrite the hydrostatic nonlinear equations of motion for $(u,v,\eta ,w,p)$ as nonlinear equations of motion for $(A_+,A_-,A_0)$.

#### C.1. Linear transforms

We separate the linear transformations into two parts $S \cdot T_\omega$ where $S : \,(A_+,A_-,A_0) \,\,\to\! (u,v,\eta ,w,p)$ maps from wave-vortex space to physical variables and $T_\omega : (A_+,A_-,A_0) \to (A_+,A_-,A_0)$ winds the wave phases from the initial conditions to the current time in wave-vortex space. The total transformations are therefore defined so that

*a*,

*b*)\begin{equation} \left[\begin{array}{@{}c@{}} u \\ v \\ \eta \\ w \\ p \end{array} \right] = S \cdot T_\omega \left[\begin{array}{@{}c@{}} \hat{A}_+ \\ \hat{A}_- \\\hat{A}_0 \end{array}\right] \quad\textrm{and}\quad \left[\begin{array}{@{}c@{}} \hat{A}_+ \\ \hat{A}_- \\\hat{A}_0 \end{array}\right] = T_\omega^{{-}1} S^{{-}1} \left[\begin{array}{c}u \\v \\ \eta \end{array}\right], \end{equation}

where it is understood that the physical variables are functions of $(x,y,z,t)$ and wave-vortex coefficients are functions of $t$ and indexed with $jkl$. The operator $T_\omega$ is the phase winding operator,

with inverse $T_\omega ^{-1}$ given by,

The projection onto physical variables is defined by,

with inverse,

These transformations use hydrostatic versions of (5.5) and (5.6), as well as the $\mathcal {DFT}$as defined in appendix B. The vertical projection operators are defined with the hydrostatic modes,

*a*,

*b*)\begin{equation} \mathcal{F} \left[ f(z) \right] \equiv \frac{1}{h_j} \int_{{-}D}^{0} f(z) F^{HS}_j \, \textrm{d}z, \quad \mathcal{G} \left[ g(z) \right] \equiv \frac{1}{g} \int_{{-}D}^{0} g(z) N^{2}(z) G^{HS}_j \, \textrm{d}z \end{equation}

with inverses

*a*,

*b*)\begin{equation} \mathcal{F}^{{-}1} \left[ \tilde{f}_{j} \right] \equiv \sum_{j} \tilde{f}_{j} F^{HS}_{j}(z), \quad \mathcal{G}^{{-}1} \left[ \tilde{g}_{j} \right] \equiv \sum_{j} \tilde{g}_{j} G^{HS}_{j}(z). \end{equation}

Note that, unlike the non-hydrostatic case, the vertical projection and $\mathcal {DFT}$ operator commute.

#### C.2. Nonlinear equations of motion

As a preliminary step before re-writing the nonlinear equations of motion, the density equation,

can be expressed using $\eta =-\rho /\partial _z\bar {\rho }$ and $N^{2}=-({g}/{\rho _0})\partial _z\bar {\rho }$ as

The nonlinear equations of motion written in matrix form are

where the nonlinear operator $\varLambda ^{{NL}}$ is defined as

*a*)\begin{gather} {uNL}\left[u,v,\eta,w \right]=u \partial_x u + v \partial_y u + w \partial_z u, \end{gather}

*b*)\begin{gather}{vNL}\left[u,v,\eta,w \right]=u \partial_x v + v \partial_y v + w \partial_z v, \end{gather}

*c*)\begin{gather}{nNL}\left[u,v,\eta,w \right]=u \partial_x \eta + v \partial_y \eta + w \left(\partial_z \eta +\eta \partial_z \ln N^{2} \right). \end{gather}

Changing to wave-vortex space, we apply the transformation $T_\omega ^{-1} S^{-1}$ to (C10),

where $\boldsymbol {A}=T_\omega ^{-1} S^{-1} \psi$ as in (C1) above. The result of this transformation are the nonlinear equations of motion in wave-vortex space,

If we had not included the time-winding operator $T_\omega$ in the basis transformation, the equations for the two wave coefficients, $A_\pm$, would include linear terms to evolve the phases, and thus the nonlinear equations would resemble forced wave equations. It is also worth noting that the first part of this transformation, the projection onto vertical modes, is identical to Kelly (Reference Kelly2016), although without the free surface.

The physical variables $(u,v,\eta ,w)$ in (C15) can be expressed in terms of wave-vortex coefficients $(A_+,A_-,A_0)$, in which case the multiplication of physical variables becomes a convolution of wave-vortex coefficients. However, we have instead left this in pseudospectral form, where the multiplication occurs in physical space and the result is transformed into wave-vortex space.

#### C.3. Energy flux

The nonlinear equations of motion (C15) each have the form $\partial _t A = F$, where $F$ is a complicated function of the other dependent variables. To compute the flux diagnostically, we take an Euler time step, for which the solution is $A(t+{\rm \Delta} t) = F(t){\rm \Delta} t + A(t)$. The change in energy is computed by multiplying by the complex conjugate and taking the time derivative so that,

The first term, dependent on the Euler time step ${\rm \Delta} t$, is only of consequence when $A$ is zero and can be neglected for these diagnostics. The squared amplitude is converted into total energy using the expression in § 3, resulting in expressions for the change in total wave and geostrophic energies,

*a*,

*b*)\begin{equation} \frac{\textrm{d}}{\textrm{d}t} E_\pm = \frac{1}{2} h \mathcal{R} \left( F_\pm \bar{A}_\pm \right) \quad\textrm{and} \quad \frac{\textrm{d}}{\textrm{d}t} E_0 = \frac{1}{2} g\left( 1+ \frac{ghK^{2}}{f_0^{2}}\right) \mathcal{R} \left( F_0 \bar{A}_0 \right). \end{equation}

#### C.4. Nonlinear pathways

While the total energy flux is computed using (C17) with (C15), this does not immediately tell us which nonlinear interactions are important. Writing out the nonlinear terms from (C15) more explicitly, we have

*a*)\begin{gather} \textrm{e}^{\textrm{i}\omega t}F_+{=} -\frac{1}{2 K} \left( k \cdot \widehat{{uNL}} + l \cdot \widehat{{vNL}} \right) - \textrm{i}\frac{\,f_0}{2\omega K} \left( l \cdot \widehat{{uNL}} - k \cdot \widehat{{vNL}} \right) + \frac{g K}{2 \omega}\widehat{{nNL}}, \end{gather}

*b*)\begin{gather}\textrm{e}^{-\textrm{i}\omega t}F_-{=} -\frac{1}{2 K} \left( k \cdot \widehat{{uNL}} + l \cdot \widehat{{vNL}} \right) + \textrm{i}\frac{\,f_0}{2\omega K} \left( l \cdot \widehat{{uNL}} - k \cdot \widehat{{vNL}} \right) - \frac{g K}{2 \omega}\widehat{{nNL}}, \end{gather}

*c*)\begin{gather}F_0 ={-}\textrm{i} h \frac{\, f_0}{\omega^{2}} \left( l \cdot \widehat{{uNL}} - k \cdot \widehat{{vNL}} \right) - \frac{\,f_0^{2}}{ \omega^{2}}\widehat{{nNL}} , \end{gather}

where we have simplified the notation so that $\widehat {{uNL}}=\mathcal {F} \cdot \mathcal {DFT}[{uNL}]$ and $\widehat {{nNL}}=\mathcal {G} \cdot \mathcal {DFT}[{nNL}]$. The idea now is to single out the specific nonlinear pathways that are move energy between modes.

Using the wave-vortex decomposition, the physical variables in $uNL$, $vNL$ and $nNL$ can be expressed in terms of their constituent parts. In this case we separate the inertial oscillations ($io$) from gravity waves ($igw$), and combine the positive and negative waves into a single internal gravity wave part. With the addition of the geostrophic component ($g$), the physical variables are unambiguously separated as,

*a*)\begin{gather} u = u_{io} + u_{igw} + u_{g} \end{gather}

*b*)\begin{gather}v = v_{io} + v_{igw} + v_{g} \end{gather}

*c*)\begin{gather}\eta = \eta_{igw} + \eta_{g} \end{gather}

*d*)\begin{gather}w = w_{igw}. \end{gather}

The primary hypothesis described in § 6 is that the advection of geostrophic vorticity by inertial oscillations is the energy source for internal gravity waves. This particular nonlinear pathway is part of the second term in parenthesis in (C18*b*) and (C18*c*),

and can be further simplified to

where $\zeta _g \equiv \partial _x v_g - \partial _y u_g$.

It is useful to treat the forcing of inertial oscillations separately from the rest of the waves because so many terms cancel. The forcing term can be found by considering the negative wave forcing (C18*c*), or with the transformation,

The total forcing on inertial waves is thus

Individual nonlinear pathways can be computed by considering the constituent parts. The two most relevant pathways for this problem are

and

where the $\overline {(\cdot )}$ indicates horizontal averaging, equivalent to the $K=0$ component of the $\mathcal {DFT}$ operator that it replaces.