A generalized wave-vortex decomposition for rotating Boussinesq flows with arbitrary stratification

Abstract The energetically independent linear wave and geostrophic (vortex) solutions are shown to be a complete basis for velocity and density variables $(u,v,w,\rho )$ in a rotating non-hydrostatic Boussinesq fluid with arbitrary stratification and non-periodic vertical boundaries. This work extends the familiar wave-vortex decomposition for triply periodic domains with constant stratification. As a consequence of the decomposition, the fluid can be unambiguously separated into decoupled linear wave and geostrophic components at each instant in time, without the need for temporal filtering. The fluid can then be diagnosed for temporal changes in wave and geostrophic coefficients at each unique wavenumber and mode, including those that inevitably occur due to nonlinear interactions. We demonstrate that this methodology can be used to determine which physical interactions cause the transfer of energy between modes by projecting the nonlinear equations of motion onto the wave-vortex basis. In the particular example given, we show that an eddy in geostrophic balance superimposed with inertial oscillations at the surface transfers energy from the inertial oscillations to internal gravity wave modes. This approach can be applied more generally to determine which mechanisms are involved in energy transfers between wave and vortices, including their respective scales. Finally, we show that the nonlinear equations of motion expressed in a wave-vortex basis are computationally efficient for certain problems. In cases where stratification profiles vary strongly with depth, this approach may be an attractive alternative to traditional spectral models for rotating Boussinesq flow.


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. 2017;Torres et al. 2018). 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 1992;Bühler, Callies & Ferrari 2014;Lien & Sanford 2019;Oscroft et al. 2021). In idealized Boussinesq models with triply periodic domains and constant stratification, such a decomposition can be made unambiguously at each instant in time (Bartello 1995;Smith & Waleffe 2002;Waite & Bartello 2006). For each resolved wavenumber the decomposition splits the flow into two inertia-gravity waves (A ± ), 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, ρ) 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 (1991) and Bartello (1995) to diagnose transfer in the turbulence cascade, and also, for example, by Arbic et al. (2012) 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 2010;Remmel, Smith & Sukhatme 2010;Hernandez-Dueñas, Smith & Stechmann 2014). 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 2010). 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 2020), 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 (2020), 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.

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 Here, p (x, y, z, t) and ρ(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 ρ tot (x, y, z, t) = ρ 0 +ρ(z) + ρ (x, y, z, t) where ∂ z p 0 (z) = −g(ρ 0 +ρ(z)). All variables in (2.1a)-(2.1e) are functions of x, y, z and t, exceptρ, which is only a function of z. We use the usual definition of buoyancy frequency, N 2 (z) ≡ −(g/ρ 0 )∂ zρ . Throughout this manuscript we use the linear approximation to isopycnal displacement η ≡ −ρ/ρ z rather than density anomaly. With this notation (2.1e) becomes w = ∂ t η and (2.1c) 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, ρ(−D) = 0. Similarly, the upper boundary is taken to be a free-slip rigid lid with w(0) = 0 and also no density anomaly, ρ(0) = 0. The depth-integrated energy densities of the flow are where HKE, VKE and 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.1a)-(2.1e), or found as the linear approximation to the available potential vorticity (APV) as defined by Wagner & Young (2015). 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 ζ is vorticity and ζ z = ∂ x v − ∂ 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 ∂ z η is not proportional to ∂ z ρ/ρ 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.1e) 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 (2015) is the relevant conserved quantity.

Wave-vortex solutions
Solutions to (2.1a)-(2.1e) are assumed to take the separable form for w and η. 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 coefficientsf jkl (t),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 domainf jkl (t), it is useful to express solutions in the frequency domain, in which case we denote the variables with·, e.g.f jkl (t) =f jkl e iω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.1e) to replace w with η t , solutions to (2.1a)-(2.1d) must satisfy 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, ω = 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ρ.
3.2. Geostrophic solutions, ω = 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.3a)-(3.3d) reduce to where F g (z), G g (z) denote the geostrophic vertical structure functions. The only equation of consequence for the vertical structure is (3.4c). 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.4c) such thatp = ρ 0 gη and N 2 G g = −g∂ 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Â 0 is a complex-valued amplitude containing the phase information, and θ 0 = kx + ly.
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 (1987)), have non-zero vertical velocities and therefore still require that three-dimensional continuity be satisfied. To maintain continuity we take (3.3d) and set the separation constant to h, such that F(z) = h∂ z G(z) for all z. This additional requirement, combined with the hydrostatic vertical momentum condition N 2 G = −g∂ z F, results in two Sturm-Liouville eigenvalue problems for hydrostatic (HS) vertical modes, 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 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 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 2015). Under rigid lid conditions, there also exists a barotropic mode ( j = 0) where F HS 0 (z) = const with no associated buoyancy anomaly, G HS 0 (z) = 0. This case will be handled separately in the decomposition.

Inertial oscillation solution
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 (2002) 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, ∂ z F I (0) and ∂ 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 Similar to the geostrophic solution where we assumed thatp = ρ 0 gη, the vertical momentum equation requires that (N 2 − ω 2 )G = −g∂ z F. Combined with continuity F = h∂ 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 ± wave solutions are given by, where the horizontal phase is given by θ ± = kx + ly ± ωt + φ and the amplitude is chosen so that depth-integrated total energy isÂ 2 h/2, as will be shown below.
Combining the vertical constraints from non-hydrostatic vertical momentum (N 2 − ω 2 )G = −g∂ z F and continuity F = h∂ z G with the dispersion relation (3.17) results in the K-constant, non-hydrostatic Sturm-Liouville problem (Early et al. 2020), The eigendepth h j and eigenfrequency ω 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 ω (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 jth internal wave mode at total wavenumber K are, which sum to a depth-integrated total energy ofÂ 2 ± h j /2. The internal wave solutions have zero potential vorticity per (2.3), PV ± = 0; but they do have Ertel PV per (2.5), again suggesting that Ertel PV may not be the appropriate quantity for this model.

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.

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.

Projection
A dynamical variable that expands in G, such as density anomaly, ρ(z), can be written as in (3.2), e.g.
where the coefficients are recovered with (4.5) The projection operation (4.5) first requires taking a Fourier transform of the variable, then invoking the orthogonality condition (4.1) with the jth vertical mode G jkl (z) for wavenumber K = √ 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 hG. For example, the divergence, However, this only works for wave solutions since the geostrophic solution does not have the same relationships between (u, v) and η. 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 ∂ z G j (z), and consider the depth-integrated quantities. That is, if which can then be projected using (4.5) to recoverũ 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.

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, η) 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. (2020) 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 The horizontal Fourier transform followed by the vertical projection then recoversŨ ijk (t),Ṽ ijk (t) andη 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,η(t) = 0, and no divergence, δ(t) = ikũ(t) + ilṽ(t) = 0. This leaves only the amplitude and phase of the vorticitỹ ζ(t) = ikṽ(t) − ilũ(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Â + ,Â − andÂ 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Â + ,Â − , andÂ 0 , There is some insight to be gained by defining depth-integrated versions of horizontal divergence and potential vorticity, Now the solution has the form, Importantly, (5.8b)-(5.8c) 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.

Summary of the decomposition
A key feature of the decomposition is that the recovered coefficients (5.4), (5.8a) 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 A + for a given j, k, l at time t = t 0 is exactly equal toÂ + 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.8a) and (5.9), by definition implies nonlinearity.

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. (2020) using the Boussinesq model described by Winters & de la Fuente (2012). The study by Lelong et al. (2020) 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. 2012), 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, , and the strength and extent of the eddy were set by parameters A, α, and β, 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. (2020), 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. (2020) and Asselin & Young (2020) and references therein. The transfer of energy from inertial oscillations to internal gravity waves was estimated using spatial-temporal averaging by Lelong et al. (2020), 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. (2020) (their figure 12b).
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. (2020), 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 1985). 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. (2020) 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. (2020) 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 ±0 . The transfer of energy is then proportional to R(F ±0Ā±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 R(F ±0Ā±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 ζ g by inertial oscillations (u io , v io ) can be written as, where F is the projection operator onto the F j modes and 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. (2020). 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).

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, ρ) 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. (2020), 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 σ (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. (2020) 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 + 1th 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 where 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. 2020).
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 (C18a). 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 5 2 log 2 N − 3N operations when transforming real variables (Canuto et al. 2006). 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 2 z 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 5 2 N z N x N y log 2 N x N y − 3N 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.

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. (2014), but now also for vertically bounded flows with variable stratification.

B.2. Vertical transformation
The finite-length sine and cosine transforms are given by where m j = jπ/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Δ where n = 0..N z and Δ = D/N z . Note that this differs from the discretization for the DFT by including endpoints. This choice of discretization results in the discrete transformŝ with vertical wavenumbers at m j = jπ/D where j = 0..N z . The sum in the DCT treats the end points separately, as they have only half the width of the other points, Δ/2. For the 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, 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 η results in the following matrices,û 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û kl0 andv kl0 , but zero forη kl0 . The DCT inflates the power by a factor of two, but the 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, phase winding operator, with inverse T −1 ω given by, The projection onto physical variables is defined by, These transformations use hydrostatic versions of (5.5) and (5.6), as well as the DFT as defined in appendix B. The vertical projection operators are defined with the hydrostatic modes, with inverses F −1 f j ≡ jf j F HS j (z), G −1 g j ≡ jg j G HS j (z). (C7a,b) Note that, unlike the non-hydrostatic case, the vertical projection and DFT operator commute.
Changing to wave-vortex space, we apply the transformation T −1 ω S −1 to (C10), ∂ t ψ + Λψ + Λ NL ψ = 0, (C12) where A = T −1 ω S −1 ψ 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 ω in the basis transformation, the equations for the two wave coefficients, A ± , 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 (2016), although without the free surface. The physical variables (u, v, η, 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 nonlinear pathway is part of the second term in parenthesis in (C18b) and (C18c), and can be further simplified to where ζ g ≡ ∂ x v g − ∂ 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 (C18c), 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 F igw∇g io = − 1 2 e i f 0 t F u igw · ∇ u g + iv g and F igw∇igw io = − 1 2 e i f 0 t F u igw · ∇ u igw + iv igw , where the (·) indicates horizontal averaging, equivalent to the K = 0 component of the DFT operator that it replaces.