Skip to main content Accessibility help
Hostname: page-component-684899dbb8-t7hbd Total loading time: 0.983 Render date: 2022-05-21T20:34:50.300Z Has data issue: true Feature Flags: { "shouldUseShareProductTool": true, "shouldUseHypothesis": true, "isUnsiloEnabled": true, "useRatesEcommerce": false, "useNewApi": true }

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

Published online by Cambridge University Press:  15 February 2021

Jeffrey J. Early*
NorthWest Research Associates, Redmond, WA98052, USA
M.P. Lelong
NorthWest Research Associates, Redmond, WA98052, USA
M.A. Sundermeyer
School for Marine Science and Technology University of Massachusetts Dartmouth, New Bedford, MA02744, USA
Email address for correspondence:
Rights & Permissions[Opens in a new window]


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.

JFM Papers
Creative Commons
Creative Common License - CCCreative Common License - BY
This is an Open Access article, distributed under the terms of the Creative Commons Attribution licence (, which permits unrestricted re-use, distribution, and reproduction in any medium, provided the original work is properly cited.
© The Author(s), 2021. Published by Cambridge University Press

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.

Figure 1. Buoyancy frequency as a function of depth for a location in the Eastern Mediterranean Sea. Black dots indicate regularly spaced grid points, while horizontal lines are the roots (Gauss quadrature points) of the 19th internal mode. Note that the vertical scale changes at 150 m depth and the deep buoyancy frequency decreases to order $10^{-3}$ cycles per hour (c.p.h.) below $\sim$900 m.

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

(2.1a)\begin{gather} \partial_t u - \,f_0 v ={-} \frac{1}{\rho_0} \partial_x p \end{gather}
(2.1b)\begin{gather}\partial_t v + \,f_0 u ={-} \frac{1}{\rho_0} \partial_y p \end{gather}
(2.1c)\begin{gather}\partial_t w ={-}\frac{1}{\rho_0} \partial_z p - g \frac{\rho}{\rho_0} \end{gather}
(2.1d)\begin{gather}\partial_x u + \partial_y v + \partial_z w = 0 \end{gather}
(2.1e)\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.1a)–(2.1e) 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.1e) becomes $w=\partial _t \eta$ 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, $\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

(2.2ac)\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,

(2.3)\begin{equation} {{\rm PV}} \equiv \partial_x v - \partial_y u - \,f_0 \partial_z \eta \end{equation}

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 (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

(2.4)\begin{align} \textrm{Ertel PV} &\equiv \frac{\left( \boldsymbol{\zeta} + \boldsymbol{\hat{k}}\,f_0\right) {\boldsymbol{\cdot}} \boldsymbol{\nabla} \rho_{tot}}{\rho_{tot}} \end{align}
(2.5)\begin{align} & \approx \frac{1}{\rho_0} \left(\zeta^{z} \bar{\rho}_z +\, f_0 \rho_z +\, f_0 \bar{\rho}_z \right), \end{align}

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

(2.6)\begin{equation} \frac{\textrm{d}}{\textrm{d}t} \left( \textrm{Ertel PV} \right) \approx \frac{1}{\rho_0} \left( \bar{\rho}_z \partial_t \zeta^{z} +\, f_0 \partial_t \rho_z +\, f_0 w \bar{\rho}_{zz} \right) \end{equation}

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),

(2.7)\begin{equation} \frac{\textrm{d}}{\textrm{d}t} \left( \textrm{Ertel PV} \right) \approx \frac{\bar{\rho}_z }{\rho_0} \frac{\partial}{\partial t}\left( \zeta^{z} +\, f_0 \frac{\partial}{\partial z} \left( \frac{\rho}{\bar{\rho}_z}\right) \right) \end{equation}

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.1a)–(2.1e) are assumed to take the separable form

(3.1)\begin{equation} f(x,y,z,t) = \sum_{jkl} \tfrac{1}{2} \tilde{f}_{jkl}(t) \exp({\textrm{i} (k x + ly )}) F_{jkl}(z) + \textrm{c.c.}, \end{equation}

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

(3.2)\begin{equation} g(x,y,z,t) =\sum_{jkl} \tfrac{1}{2} \tilde{g}_{jkl}(t) \exp({\textrm{i} (k x + ly )}) G_{jkl}(z) + \textrm{c.c.}, \end{equation}

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.1e) to replace $w$ with $\eta _t$, solutions to (2.1a)–(2.1d) must satisfy

(3.3a)\begin{gather} \textrm{i} \omega \hat{u} -\, f_0 \hat{v} ={-} \textrm{i} k \frac{\hat{p}}{\rho_0} \end{gather}
(3.3b)\begin{gather}\textrm{i} \omega \hat{v} +\, f_0 \hat{u} ={-} \textrm{i} l \frac{\hat{p}}{\rho_0} \end{gather}
(3.3c)\begin{gather}\left(N^{2} - \omega^{2} \right) \hat{\eta} G ={-}\frac{\hat{p}}{\rho_0} \partial_z F \end{gather}
(3.3d)\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.3a)–(3.3d) reduce to

(3.4a)\begin{gather} -\, f_0 \hat{v} ={-} \textrm{i} k \frac{\hat{p}}{\rho_0} , \end{gather}
(3.4b)\begin{gather}f_0 \hat{u} ={-} \textrm{i} l \frac{\hat{p}}{\rho_0} , \end{gather}
(3.4c)\begin{gather}N^{2} \hat{\eta} G_g ={-} \frac{\hat{p}}{\rho_0} \partial_z F_g , \end{gather}
(3.4d)\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.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 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,

(3.5)\begin{equation} \left[\begin{array}{@{}c@{}} u_g \\ v_g \\ \eta_g \end{array}\right] = \frac{\hat{A}_0}{2} \left[\begin{array}{@{}c@{}} -\textrm{i} \dfrac{g}{f_0} l F_g(z) \\ \textrm{i} \dfrac{g}{f_0} k F_g(z) \\ G_g(z) \end{array}\right] \textrm{e}^{\textrm{i}\theta_0} + \textrm{c.c.}, \end{equation}

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.3d) 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,

(3.6)\begin{equation} \frac{\textrm{d}^{2} G^{HS}_j}{\textrm{d}z^{2}} ={-} \frac{N^{2}}{gh_j} G^{HS}_j \end{equation}

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

(3.7)\begin{equation} \frac{\textrm{d}}{\textrm{d}z} \left( \frac{1}{N^{2}} \frac{\textrm{d} F^{HS}_j}{\textrm{d}z} \right) ={-} \frac{1}{gh_j} F^{HS}_j \end{equation}

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

(3.8)\begin{equation} \frac{1}{g} \int_{{-}D}^{0} N^{2}(z) G^{HS}_i G^{HS}_j \, \textrm{d}z = \delta_{ij}, \end{equation}


(3.9)\begin{equation} \int_{{-}D}^{0} F^{HS}_i F^{HS}_j \, \textrm{d}z = h_i \delta_{ij}, \end{equation}

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,

(3.10)\begin{gather} {{\rm HKE}}_g=\frac{\hat{A}_0^{2}}{4} \frac{g^{2}}{f_0^{2}}K^{2} \int_{{-}D}^{0} F_g^{2}(z) \, \textrm{d}z \end{gather}
(3.11)\begin{gather}{{\rm PE}}_g=\frac{\hat{A}_0^{2}}{4} \int_{{-}D}^{0} N^{2}(z) G_g^{2}(z)\, \textrm{d}z, \end{gather}

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,

(3.12)\begin{equation} {{\rm PV}}_g={-}\frac{\hat{A}_0 g}{2\, f_0 } \left( K^{2}F_g(z) - \frac{\textrm{d}}{\textrm{d}z} \left( \frac{\,f_0^{2}}{N^{2}} \frac{\textrm{d} F_g(z)}{\textrm{d}z} \right) \right) \textrm{e}^{\textrm{i}\theta_0} + \textrm{c.c.}, \end{equation}

as is traditionally written, or simply

(3.13)\begin{equation} {{\rm PV}}_g ={-}\frac{\hat{A}_0}{2\, f_0 h} \left( ghK^{2} +\, f_0^{2} \right) F^{HS}(z) \,\textrm{e}^{\textrm{i}\theta_0} + \textrm{c.c.}, \end{equation}

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,

(3.14)\begin{equation} \textrm{Ertel PV}_g= \frac{\bar{\rho}_z}{\rho_0} \left[ {{\rm PV}}_g - \frac{\hat{A}_0\, f_0}{2} \left( \partial_z \ln \bar{\rho}_z \right) G(z) \,\textrm{e}^{\textrm{i} \theta_0} +\, f_0 \right] + \textrm{c.c.}, \end{equation}

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,

(3.15)\begin{equation} \left[\begin{array}{@{}c@{}} u_I \\ v_I \\ \eta_I \end{array}\right] = \left[\begin{array}{@{}c@{}} U_I \cos(\,f_0 t + \phi_0) F_I(z) \\ -U_I \sin(\,f_0 t + \phi_0) F_I(z) \\ 0 \end{array}\right]. \end{equation}

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

(3.16)\begin{equation} \left[\begin{array}{@{}ccc@{}} \textrm{i} \omega & -\,f_0 & {\textrm{i}} g k \\ f_0 & \textrm{i} \omega & \textrm{i} g l \\ k h & l h & \omega \end{array}\right] \left[\begin{array}{@{}c@{}} \hat{u} \\ \hat{v} \\ \hat{\eta} \end{array}\right] = \left[\begin{array}{@{}c@{}} 0 \\ 0 \\ 0 \end{array}\right]. \end{equation}

This system of equations admits the internal wave solutions when

(3.17)\begin{equation} \omega = \sqrt{g h K^{2} +\, f_0^{2}}. \end{equation}

The $\pm$ wave solutions are given by,

(3.18)\begin{equation} \left[\begin{array}{@{}c@{}} u_\pm \\ v_\pm \\ \eta_\pm \end{array}\right] = \frac{\hat{A}_\pm}{2} \left[\begin{array}{@{}c@{}} \dfrac{k\omega \mp \textrm{i}l\, f_0}{\omega K} F(z) \\ \dfrac{l \omega \pm \textrm{i} k\, f_0}{\omega K } F(z) \\ \mp \dfrac{K h}{\omega} G(z) \end{array}\right] \textrm{e}^{\textrm{i} \theta_\pm} + \textrm{c.c.}, \end{equation}

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),

(3.19)\begin{equation} \frac{\textrm{d}^{2} G_j}{\textrm{d}z^{2}} - K^{2} G_j ={-}\frac{N^{2} -\, f_0^{2}}{g h_j }G_j. \end{equation}

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,

(3.20)\begin{gather} {\rm HKE}_\pm{=} \frac{\hat{A}_\pm^{2}}{4} \left( 1+ \frac{f_0^{2}}{\omega_j^{2}} \right) \int_{{-}D}^{0} F_j^{2}(z) \, \textrm{d}z \end{gather}
(3.21)\begin{gather}{\rm VKE}_\pm{=} \frac{\hat{A}_\pm^{2}}{4} K^{2} h_j^{2} \int_{{-}D}^{0} G_j^{2}(z) \, \textrm{d}z \end{gather}
(3.22)\begin{gather}{\rm PE}_\pm{=} \frac{\hat{A}_\pm^{2}}{4} \frac{K^{2} h_j^{2}}{\omega_j^{2}} \int_{{-}D}^{0} N^{2}(z) G_j^{2}(z) \, \textrm{d}z \end{gather}

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),

(3.23)\begin{equation} \textrm{Ertel PV}_{{\pm}} = \frac{\bar{\rho}_z}{\rho_0} \left[{\pm} \frac{\hat{A}_\pm}{2 } \frac{K h_j\, f_0}{\omega_j} ( \partial_z \ln \bar{\rho}_z ) G(z) \,\textrm{e}^{\textrm{i} \theta_\pm} +\, f_0 \right] + \textrm{c.c.}, \end{equation}

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,

(4.1)\begin{equation} \frac{1}{g} \int_{{-}D}^{0} (N^{2}(z)-\,f_0^{2}) G_i G_j \, \textrm{d}z = \delta_{ij}, \end{equation}

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

(4.2)\begin{equation} \partial_z \left( \frac{\partial_z F_j}{N^{2} - g h_j K^{2} -\, f_0^{2}} \right) ={-}\frac{1}{g h_j} F_j, \end{equation}

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

(4.3)\begin{equation} \int_{{-}D}^{0} \left( F_i F_j + h_i h_j K^{2} G_i G_j \right) \textrm{d}z = h_i \delta_{ij}. \end{equation}

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.

(4.4)\begin{equation} \rho(x,y,z,t) = \sum_{jkl}\tfrac{1}{2} \tilde{\rho}_{jkl}(t) G_{jlk}(z) \exp({\textrm{i}(kx+ly)}) + \textrm{c.c.}, \end{equation}

where the coefficients are recovered with

(4.5)\begin{equation} \tilde{\rho}_{jkl}(t) = \frac{1}{g} \int_{{-}D}^{0} (N^{2}(z) -\, f_0^{2}) \left[ \frac{1}{N_x N_y} \sum_{xy} \rho(x,y,z,t) \exp({-\textrm{i}(kx+ly)}) \right] G_{jkl}(z) \, \textrm{d}z. \end{equation}

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,

(4.6)\begin{gather} \delta_j(t) = \int_{{-}D}^{0} \left( \delta(t) F_j(z) - \textrm{i} K^{2} w(t) h_j G_j(z) \right) \textrm{d}z , \end{gather}
(4.7)\begin{gather}\zeta_j(t) = \int_{{-}D}^{0} \left( \zeta(t) F_j(z) - \textrm{i}\, f_0 K^{2} \eta(t) h_j G_j(z) \right) \textrm{d}z, \end{gather}

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

(4.8)\begin{equation} u(x,y,z,t) = \sum_{j,k,l} \frac{\tilde{u}_{jkl}(t)}{2} \exp({\textrm{i} (k x + ly)}) F_{jkl}(z) + \textrm{c.c.}, \end{equation}

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

(4.9)\begin{equation} U(x,y,z,t) = \sum_{j,k,l} \frac{\tilde{u}_{jkl}(t)}{2} \exp({\textrm{i} (k x + ly)}) h_{ijk} G_{jkl}(z) + \textrm{c.c.}, \end{equation}

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

(5.1)\begin{gather} U(x,y,z,t) = \sum_{j,k,l} \frac{\tilde{U}_{jkl}(t)}{2} \exp({\textrm{i} (k x + ly)}) G_{jkl}(z) + \textrm{c.c.} \end{gather}
(5.2)\begin{gather}V(x,y,z,t) = \sum_{j,k,l} \frac{\tilde{V}_{jkl}(t)}{2} \exp({\textrm{i} (k x + ly)}) G_{jkl}(z) + \textrm{c.c.} \end{gather}
(5.3)\begin{gather}\eta(x,y,z,t) = \sum_{j,k,l} \frac{\tilde{\eta}_{jkl}(t)}{2} \exp({\textrm{i} (k x + ly)}) G_{jkl}(z) + \textrm{c.c.}, \end{gather}

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,

(5.4)\begin{equation} \hat{A}_0 ={-} \textrm{i}\frac{f_0}{g K^{2}} \left( k \tilde{v}(t) - l \tilde{u}(t) \right), \end{equation}

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

(5.5)\begin{equation} \left[\begin{array}{@{}c@{}}\tilde{U}(t) \\ \tilde{V}(t) \\ \tilde{\eta}(t)\end{array}\right] = \left[\begin{array}{@{}ccc@{}} \dfrac{k\omega - \textrm{i}l\, f_0}{\omega K}h & \dfrac{k\omega + \textrm{i}l\, f_0}{\omega K}h & -\textrm{i} \dfrac{gh}{f_0} l \\ \dfrac{l \omega + \textrm{i} k\, f_0}{\omega K }h & \dfrac{l\omega - \textrm{i} k\, f_0}{\omega K }h & \textrm{i}\dfrac{gh}{f_0} k \\ - \dfrac{K h}{\omega} & \dfrac{K h}{\omega} & 1\end{array}\right] \left[\begin{array}{@{}c@{}} \textrm{e}^{\textrm{i} \omega t} \hat{A}_+ \\ \textrm{e}^{-\textrm{i} \omega t} \hat{A}_- \\ \hat{A}_0 \end{array}\right] \end{equation}

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

(5.6a)\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}
(5.6b)\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,

(5.7a,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,

(5.8a)\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}
(5.8b)\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}
(5.8c)\begin{gather}\hat{A}_0 ={-} \frac{f}{\omega^{2}} \tilde{\varPi}(t). \end{gather}

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.

(5.9a)\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}
(5.9b)\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}
(5.9c)\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.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 $\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.8a) 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,

(6.1)\begin{equation} \psi ={-}\frac{1}{4} \frac{A}{\alpha} \exp({-2 \alpha (x^{2} + y^{2}) - 2\beta z^{2}}), \end{equation}

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

(6.2)\begin{equation} u_I = U_I \,\textrm{e}^{-\gamma z}, v_I = 0, \end{equation}

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 12b).

Figure 2. Energy time series for the Cyprus eddy simulation inferred via the wave-vortex decomposition, showing that total energy and geostrophic energy are approximately conserved, while inertial energy decreases and internal gravity wave increases by comparable amounts. Residual energy increases slightly (from 0.2 % to 0.3 %) during the same period.

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.

Figure 3. Change in energy of inertial, wave and geostrophic components as a function of mode and horizontal scale on day 6 of the Cyprus eddy simulation. Contours on the wave plot indicate the frequency of oscillation, in units of $f_0$. The change in energy is dominated by a sustained loss of energy in the $j=1$ inertial mode, and a significant gain in wave energy at scales of approximately 250 km. Changes in geostrophic energy are an order of magnitude smaller (note the log colour scale), and rapidly oscillate between signs with no sustained gain or loss.

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,

(6.3)\begin{equation} \partial_t \left[\begin{array}{@{}c@{}} \hat{A}_+ \\ \hat{A}_- \\\hat{A}_0 \end{array}\right] = \left[\begin{array}{@{}c@{}} F_+ \\ F_- \\ F_0 \end{array} \right], \end{equation}

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.

Figure 4. Total change in energy (black) and computed total flux (dashed black) of (a) the internal gravity wave energy and (b) inertial energy. Panel (a) also shows the flux from inertial oscillations advecting geostrophic motions (blue). Panel (b) includes the flux from internal gravity waves advecting geostrophic motions (orange), and self-advection of internal gravity waves (purple).

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,

(6.4)\begin{equation} F_\pm^{{io}\boldsymbol{\nabla}{g}} ={\pm} \frac{f_0 \,\textrm{e}^{{\mp} \textrm{i}\omega t}}{2 \omega K} \left( \mathcal{F} \cdot \mathcal{DFT}\left[u_{io} \partial_x \zeta_{g} + v_{io} \partial_y \zeta_{g} \right] \right), \end{equation}

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).

Figure 5. Same as figure 3, but showing only the two dominant mechanisms from figure 4. The wave panel shows flux from inertial oscillations advecting geostrophic motions and the inertial panel shows flux from internal gravity waves advecting geostrophic motions.

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 (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 $\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


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


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,

(A1a)\begin{gather} G^{N_0}_j(z) = A \sin \left( m_j ( z + D) \right) \end{gather}
(A1b)\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,

(A2a)\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}
(A2b)\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

(A3a)\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}
(A3b)\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.8b)–(5.8c) 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$.

(A4)\begin{gather} \hat{A}_+{=} \frac{\textrm{e}^{-\textrm{i} \omega t}}{2 K} \left[ - \textrm{i}\tilde{\delta}(t) - \frac{1}{\omega} \left( f \widetilde{{{\rm PV}}}(t) + \frac{\omega^{2}}{h} \tilde{\eta}(t) \right) \right], \end{gather}
(A5)\begin{gather}\hat{A}_-{=} \frac{\textrm{e}^{\textrm{i} \omega t}}{2 K} \left[ - \textrm{i}\tilde{\delta}(t) + \frac{1}{\omega} \left( f \widetilde{{{\rm PV}}}(t) + \frac{\omega^{2}}{h} \tilde{\eta}(t) \right) \right], \end{gather}
(A6)\begin{gather}\hat{A}_0 ={-} \frac{\,f_0 h}{\omega^{2}} \widetilde{{{\rm PV}}}(t), \end{gather}


(A7a,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,

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

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

  3. (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

(B1)\begin{equation} \mathcal{F}\left[\, f(x) \right] = \frac{1}{L} \int_{{-}D}^{0} f(x) \exp({-\textrm{i} k_j x}) \, {\textrm{d} x}, \end{equation}

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

(B2)\begin{equation} \hat{f}(k_j) = \mathcal{DFT}\left[\, f(x_n) \right] = \frac{1}{L} \sum_{n=0}^{N-1} f(x_n) \exp({-\textrm{i} k_j x_n}) \varDelta \end{equation}

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

(B3)\begin{equation} \frac{1}{N}\sum_{n=0}^{N-1} |\,f(x_n)|^{2} = \sum_{j=0}^{N-1} S(k_j) \,\textrm{d}k, \end{equation}

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.

Table 1. Table of FFT coefficients for wavenumbers $(k,l)$. We have let $-4 \leqslant k \leqslant 3$ and $-4 \leqslant l \leqslant 3$, consistent with an $8\times 8$ 2D FFT. The grey shaded region shows the redundant coefficients that are determined from Hermitian conjugacy by changing the sign on the $l$ component. The pink shaded components are Hermitian conjugate by changing the sign on the $k$ component. The orange components are Nyquist components, and thus not full resolved. The cyan shaded component (including $(k,l)=(-4,0)$, $(-4,4)$ and $(0,-4)$) are self-symmetric, and therefore strictly real.

B.2. Vertical transformation

The finite-length sine and cosine transforms are given by

(B4)\begin{equation} \mathcal{S}\left[ g(z) \right] = \frac{2}{D} \int_{{-}D}^{0} g(z) \sin( m_j z) \, \textrm{d}z \end{equation}


(B5)\begin{equation} \mathcal{C}\left[ f(z) \right] = \frac{2}{D} \int_{{-}D}^{0} f(z) \cos( m_j z) \, \textrm{d}z, \end{equation}

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

(B6)\begin{equation} \hat{g}(m_j) = \mathcal{DST}\left[ g(z_n) \right] = \frac{2}{D} \sum_{n=0}^{N_z} g(z_n) \sin( m_j z_n ) \varDelta \end{equation}


(B7)\begin{equation} \hat{f}(m_j) = \mathcal{DCT}\left[ f(z_n) \right] = \frac{2}{D}\left( \frac{f(0)}{2} + \sum_{n=1}^{N_z-1} f(z_n) \cos( m_j z_n ) + ({-}1)^{j} \frac{f(D)}{2}\right) \varDelta \end{equation}

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,

(B8)\begin{equation} \frac{1}{N_z}\sum_{n=1}^{N_z-1} |g(x_n)|^{2} = \sum_{j=1}^{N_z-1} S(m_j) \,\textrm{d}m \end{equation}

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

(B9)\begin{equation} \frac{1}{N_z} \left( \frac{|\,f(0)|^{2}}{2} + \sum_{n=1}^{N_z-1} |\,f(z_n)|^{2} + \frac{|\,f(D)|^{2}}{2} \right)= \left(\frac{S(m_0)}{2} + \sum_{j=1}^{N_z-1} S(m_j) + 2S(m_{Nz}) \right) \textrm{d}m \end{equation}

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,

(B10a)\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}
(B10b)\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}
(B10c)\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,

(B11)\begin{equation} \hat{A}_0 ={-} \textrm{i} \frac{f_0}{g K^{2}} \left( k \hat{v}_{kl0} - l \hat{u}_{kl0} \right) \end{equation}

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,

(B12a)\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}
(B12b)\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}
(B12c)\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,

(B13)\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} \bar{u}_{klj} + \frac{l \omega \mp \textrm{i} k\, f_0}{\omega K}\bar{v}_{klj} \mp \frac{gK}{\omega} \bar{\eta}_{klj} \right] \end{gather}
(B14)\begin{gather}\hat{A}_0 = \textrm{i} \frac{l h\, f_0}{\omega^{2}} \bar{u}_{klj} - \textrm{i} \frac{k h\, f_0}{\omega^{2}} \bar{v}_{klj} + \frac{\, f_0^{2}}{ \omega^{2}} \bar{\eta}_{klj}. \end{gather}

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

(C1a,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,

(C2)\begin{equation} T_\omega = \left[\begin{array}{@{}ccc@{}} \textrm{e}^{\textrm{i}\omega t} & 0 & 0 \\ 0 & \textrm{e}^{-\textrm{i}\omega t} & 0 \\ 0 & 0 & 1 \end{array}\right] \end{equation}

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

(C3)\begin{equation} T_\omega^{{-}1} = \left[\begin{array}{@{}ccc@{}} \textrm{e}^{-\textrm{i}\omega t} & 0 & 0 \\ 0 & \textrm{e}^{\textrm{i}\omega t} & 0 \\ 0 & 0 & 1 \end{array}\right]. \end{equation}

The projection onto physical variables is defined by,

(C4)\begin{equation} S = \mathcal{DFT}^{{-}1} \cdot \left[\begin{array}{@{}ccccc@{}} \mathcal{F}^{{-}1} & 0 & 0 & 0 & 0 \\ 0 & \mathcal{F}^{{-}1} & 0 & 0 & 0 \\ 0 & 0 & \mathcal{G}^{{-}1} & 0 & 0 \\ 0 & 0 & 0 & \mathcal{G}^{{-}1} & 0 \\ 0 & 0 & 0 & 0 & \mathcal{F}^{{-}1} \end{array}\right] \left[\begin{array}{@{}ccc@{}} \dfrac{k\omega - \textrm{i}l\, f_0}{\omega K} & \dfrac{k\omega + \textrm{i}l\, f_0}{\omega K} & -\textrm{i} \dfrac{g}{f_0} l \\ \dfrac{l \omega + \textrm{i} k\, f_0}{\omega K } & \dfrac{l\omega - \textrm{i} k\, f_0}{\omega K } & \textrm{i}\dfrac{g}{f_0} k \\ - \dfrac{K h}{\omega} & \dfrac{K h}{\omega} & 1 \\ -\textrm{i}Kh & -\textrm{i}Kh & 0 \\ -\rho_0 g \dfrac{K h}{\omega} & \rho_0 g \dfrac{K h}{\omega} & \rho_0 g \end{array}\right] \end{equation}

with inverse,

(C5)\begin{equation} S^{{-}1} = \left[\begin{array}{@{}ccc@{}} \dfrac{k \omega + \textrm{i} l\, f_0}{2\omega K} & \dfrac{l \omega - \textrm{i} k\, f_0}{2\omega K} & - \dfrac{gK}{2\omega} \\ \dfrac{k \omega - \textrm{i} l\, f_0}{2\omega K} & \dfrac{l \omega + \textrm{i} k\, f_0}{2\omega K} & \dfrac{gK}{2\omega} \\ \textrm{i} \dfrac{l h\, f_0}{\omega^{2}} & - \textrm{i} \dfrac{k h\, f_0}{\omega^{2}} & \dfrac{ f_0^{2}}{ \omega^{2}} \end{array}\right] \cdot \left[\begin{array}{@{}ccc@{}} \mathcal{F} & 0 & 0 \\ 0 & \mathcal{F} & 0 \\ 0 & 0 & \mathcal{G} \end{array}\right] \cdot \mathcal{DFT}. \end{equation}

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,

(C6a,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

(C7a,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,

(C8)\begin{equation} \partial_t \rho + w \partial_z \bar{\rho} + u \partial_x \rho + v \partial_y \rho + w \partial_z \rho = 0 \end{equation}

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

(C9)\begin{equation} \partial_t \eta - w + u \partial_x\eta + v \partial_y\eta + w \left(\partial_z\eta +\eta \partial_z \ln N^{2} \right) = 0. \end{equation}

The nonlinear equations of motion written in matrix form are

(C10)\begin{equation} \partial_t \underbrace{ \left[\begin{array}{@{}c@{}} u \\ v \\ \eta \\ 0 \\ 0 \end{array} \right]}_{\psi} + \underbrace{ \left[\begin{array}{@{}ccccc@{}} 0 & -\,f_0 & 0 & 0 & \dfrac{1}{\rho_0}\partial_x \\ f_0 & 0 & 0 & 0 & \dfrac{1}{\rho_0}\partial_y \\ 0 & 0 & 0 & -1 & 0 \\ 0 & 0 & 1 & 0 & N^{{-}2}\dfrac{1}{\rho_0}\partial_z \\ \partial_x & \partial_y & 0 & \partial_z & 0 \end{array}\right] }_{\varLambda} \underbrace{ \left[\begin{array}{@{}c@{}} u \\ v \\ \eta \\ w \\ p \end{array} \right]}_{\psi}+ \underbrace{ \left[\begin{array}{@{}c@{}} {uNL} \\ {vNL} \\ {nNL} \\ 0 \\ 0 \end{array} \right]}_{\varLambda^{{NL}}}\underbrace{ \left[\begin{array}{@{}c@{}} u \\ v \\ \eta \\ w \\ p \end{array} \right]}_{\psi}= 0, \end{equation}

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

(C11a)\begin{gather} {uNL}\left[u,v,\eta,w \right]=u \partial_x u + v \partial_y u + w \partial_z u, \end{gather}
(C11b)\begin{gather}{vNL}\left[u,v,\eta,w \right]=u \partial_x v + v \partial_y v + w \partial_z v, \end{gather}
(C11c)\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),

(C12)\begin{gather} \partial_t \psi + \varLambda \psi + \varLambda^{{NL}} \psi = 0, \end{gather}
(C13)\begin{gather}T_\omega^{{-}1} S^{{-}1}\partial_t \psi + T_\omega^{{-}1} S^{{-}1}\varLambda \psi + T_\omega^{{-}1} S^{{-}1}\varLambda^{{NL}} \psi = 0 , \end{gather}
(C14)\begin{gather}\left(T_\omega^{{-}1} S^{{-}1}\partial_t S T_\omega\right) \boldsymbol{A} + \left(T_\omega^{{-}1} S^{{-}1}\varLambda S T_\omega\right) \boldsymbol{A} + T_\omega^{{-}1} S^{{-}1}\varLambda^{{NL}} \psi = 0, \end{gather}

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,

(C15)\begin{equation} \partial_t \left[\begin{array}{@{}c@{}} \hat{A}_+ \\ \hat{A}_- \\\hat{A}_0 \end{array}\right] = \left[\begin{array}{@{}c@{}} F_+ \\ F_- \\ F_0 \end{array} \right]\quad \textrm{where}\ \left[\begin{array}{@{}c@{}} F_+ \\ F_- \\ F_0 \end{array} \right] \equiv{-}T_\omega^{{-}1} S^{{-}1}\left[\begin{array}{@{}c@{}} {uNL} \\ {vNL} \\ {nNL} \end{array} \right] \psi. \end{equation}

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,

(C16)\begin{equation} \frac{\textrm{d}}{\textrm{d}t} A^{2} = 2 F \bar{F} {\rm \Delta} t + 2 \mathcal{R} \left( F \bar{A} \right). \end{equation}

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,

(C17a,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

(C18a)\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}
(C18b)\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}
(C18c)\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,

(C19a)\begin{gather} u = u_{io} + u_{igw} + u_{g} \end{gather}
(C19b)\begin{gather}v = v_{io} + v_{igw} + v_{g} \end{gather}
(C19c)\begin{gather}\eta = \eta_{igw} + \eta_{g} \end{gather}
(C19d)\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 (C18b) and (C18c),

(C20)\begin{align} F_\pm^{{io}\boldsymbol{\nabla}{g}} &={\pm} \frac{\textrm{i}\, f_0 \,\textrm{e}^{{\mp} \textrm{i}\omega t}}{2 \omega K} \left( k \cdot \mathcal{F} \cdot \mathcal{DFT} \left[u_{io} \partial_x v_{g} + v_{io} \partial_y v_{g} \right] \right.\nonumber\\ &\quad \left.- l \cdot \mathcal{F} \cdot \mathcal{DFT}\left[u_{io} \partial_x u_{g} + v_{io} \partial_y u_{g} \right] \right), \end{align}

and can be further simplified to

(C21)\begin{equation} F_\pm^{{io}\boldsymbol{\nabla}{g}} ={\pm} \frac{\,f_0 \,\textrm{e}^{{\mp} \textrm{i}\omega t}}{2 \omega K} \left( \mathcal{F} \cdot \mathcal{DFT}\left[u_{io} \partial_x \zeta_{g} + v_{io} \partial_y \zeta_{g} \right] \right), \end{equation}

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 (C18c), or with the transformation,

(C22)\begin{equation} S^{{-}1} = \frac{1}{2} \left[\begin{array}{@{}ccc@{}} 0 & 0 & 0 \\ 1 & i & 0 \\ 0 & 0 & 0 \end{array}\right]. \end{equation}

The total forcing on inertial waves is thus

(C23)\begin{equation} F_{io} ={-}\tfrac{1}{2} \textrm{e}^{\textrm{i}\,f_0 t} \left( \widehat{{uNL}} + \textrm{i} \widehat{{vNL}} \right). \end{equation}

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

(C24)\begin{equation} F_{io}^{{igw}\boldsymbol{\nabla}{g}} ={-}\tfrac{1}{2} \textrm{e}^{\textrm{i}\,f_0 t} \mathcal{F} \left[ \overline{\boldsymbol{u}_{igw} \boldsymbol{\cdot} \boldsymbol{\nabla} \left( u_{g} + \textrm{i} v_{g} \right) } \right] \end{equation}


(C25)\begin{equation} F_{io}^{{igw}\boldsymbol{\nabla}{igw}} ={-}\tfrac{1}{2} \textrm{e}^{\textrm{i}\,f_0 t} \mathcal{F} \left[ \overline{\boldsymbol{u}_{igw} \boldsymbol{\cdot} \boldsymbol{\nabla} \left( u_{igw} + \textrm{i} v_{igw} \right) } \right], \end{equation}

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


Arbic, B.K., Scott, R.B., Flierl, G.R., Morten, A.J., Richman, J.G. & Shriver, J.F. 2012 Nonlinear cascades of surface oceanic geostrophic kinetic energy in the frequency domain. J. Phys. Oceanogr. 42 (9), 15771600.CrossRefGoogle Scholar
Asselin, O. & Young, W.B. 2020 Penetration of wind-generated near-inertial waves into a turbulent ocean. J. Phys. Oceanogr. 50 (6), 16991716.CrossRefGoogle Scholar
Bartello, P. 1995 Geostrophic adjustment and inverse cascades in rotating stratified turbulence. J. Atmos. Sci. 52 (24), 44104428.2.0.CO;2>CrossRefGoogle Scholar
Bühler, O., Callies, J. & Ferrari, R. 2014 Wave vortex decomposition of one-dimensional ship-track data. J. Fluid Mech. 756, 10071026.CrossRefGoogle Scholar
Canuto, C., Hussaini, M.Y., Quarteroni, A. & Zang, T.A. 2006 Spectral Methods: Fundamentals in Single Domains. Springer.CrossRefGoogle Scholar
Cuypers, Y., Bouruet-Aubertot, P., Marec, C. & Fuda, J.-L. 2012 Characterization of turbulence from a fine-scale parameterization and microstructure measurements in the Mediterranean Sea during the BOUM experiment. Biogeosciences 9 (8), 31313149.CrossRefGoogle Scholar
Early, J.J., Lelong, M.P. & Smith, K.S. 2020 Fast and accurate computation of vertical modes. J. Adv. Model. Earth Syst. 12 (2), 126.CrossRefGoogle Scholar
Hernandez-Dueñas, G., Smith, L.M. & Stechmann, S.N. 2014 Investigation of Boussinesq dynamics using intermediate models based on wave–vortical interactions. J. Fluid Mech. 747, 247287.CrossRefGoogle Scholar
Kelly, S.M. 2016 The vertical mode decomposition of surface and internal tides in the presence of a free surface and arbitrary topography. J. Phys. Oceanogr. 46 (12), 37773788.CrossRefGoogle Scholar
Kunze, E. 1985 Near-inertial wave propagation in geostrophic shear. J. Phys. Oceanogr. 15 (5), 544565.2.0.CO;2>CrossRefGoogle Scholar
Lelong, M.-P. & Riley, J.J. 1991 Internal wave-vortical mode interactions in strongly stratified flows. J. Fluid Mech. 232, 119.CrossRefGoogle Scholar
Lelong, M.P., Cuypers, Y. & Bouruet-Aubertot, P. 2020 Near-inertial energy propagation inside a mediterranean anticyclonic eddy. J. Phys. Oceanogr. 50 (8), 22712288.CrossRefGoogle Scholar
Lien, R.-C. & Müller, P. 1992 Normal-mode decomposition of small-scale oceanic motions. J. Phys. Oceanogr. 22 (12), 15831595.2.0.CO;2>CrossRefGoogle Scholar
Lien, R.-C. & Sanford, T.B. 2019 Small-scale potential vorticity in the upper ocean thermocline small-scale potential vorticity in the upper ocean thermocline. J. Phys. Oceanogr. 49 (7), 18451872.CrossRefGoogle Scholar
Oscroft, S., Sykulski, A.M. & Early, J.J. 2021 Separating mesoscale and submesoscale flows from clustered drifter trajectories Fluids 6 (1), 134.Google Scholar
Pedlosky, J. 1987 Geophysical Fluid Dynamics. Springer.CrossRefGoogle Scholar
Remmel, M. 2010 New models for the rotating shallow water and Boussinesq equations by subsets of mode interactions. PhD thesis.Google Scholar
Remmel, M., Smith, L.M. & Sukhatme, J. 2010 Nonlinear inertia-gravity wave-mode interactions in three dimensional rotating stratified flows. Commun. Math. Sci. 8 (2), 357376.CrossRefGoogle Scholar
Savage, A.C., et al. . 2017 Spectral decomposition of internal gravity wave sea surface height in global models. J. Geophys. Res.-Oceans 122 (10), 78037821.CrossRefGoogle Scholar
Smith, K.S. & Vanneste, J. 2013 A surface-aware projection basis for quasigeostrophic flow. J. Phys. Oceanogr. 43 (3), 548562.CrossRefGoogle Scholar
Smith, L.M. & Waleffe, F. 2002 Generation of slow large scales in forced rotating stratified turbulence. J. Fluid Mech. 451, 145168.CrossRefGoogle Scholar
Torres, H.S., Klein, P., Menemenlis, D., Qiu, B., Su, Z., Wang, J., Chen, S. & Fu, L.-L. 2018 Partitioning ocean motions into balanced motions and internal gravity waves: a modeling study in anticipation of future space missions. J. Geophys. Res.-Oceans 123 (11), 80848105.CrossRefGoogle Scholar
Wagner, G.L. & Young, W.R. 2015 Available potential vorticity and wave-averaged quasi-geostrophic flow. J. Fluid Mech. 785, 401424.CrossRefGoogle Scholar
Waite, M.L. & Bartello, P. 2006 The transition from geostrophic to stratified turbulence. J. Fluid Mech. 568, 89108.CrossRefGoogle Scholar
Winters, K.B. & de la Fuente, A. 2012 Modelling rotating stratified flows at laboratory-scale using spectrally-based DNS. Ocean Model. 49–50(C), 4759.CrossRefGoogle Scholar
Figure 0

Figure 1. Buoyancy frequency as a function of depth for a location in the Eastern Mediterranean Sea. Black dots indicate regularly spaced grid points, while horizontal lines are the roots (Gauss quadrature points) of the 19th internal mode. Note that the vertical scale changes at 150 m depth and the deep buoyancy frequency decreases to order $10^{-3}$ cycles per hour (c.p.h.) below $\sim$900 m.

Figure 1

Figure 2. Energy time series for the Cyprus eddy simulation inferred via the wave-vortex decomposition, showing that total energy and geostrophic energy are approximately conserved, while inertial energy decreases and internal gravity wave increases by comparable amounts. Residual energy increases slightly (from 0.2 % to 0.3 %) during the same period.

Figure 2

Figure 3. Change in energy of inertial, wave and geostrophic components as a function of mode and horizontal scale on day 6 of the Cyprus eddy simulation. Contours on the wave plot indicate the frequency of oscillation, in units of $f_0$. The change in energy is dominated by a sustained loss of energy in the $j=1$ inertial mode, and a significant gain in wave energy at scales of approximately 250 km. Changes in geostrophic energy are an order of magnitude smaller (note the log colour scale), and rapidly oscillate between signs with no sustained gain or loss.

Figure 3

Figure 4. Total change in energy (black) and computed total flux (dashed black) of (a) the internal gravity wave energy and (b) inertial energy. Panel (a) also shows the flux from inertial oscillations advecting geostrophic motions (blue). Panel (b) includes the flux from internal gravity waves advecting geostrophic motions (orange), and self-advection of internal gravity waves (purple).

Figure 4

Figure 5. Same as figure 3, but showing only the two dominant mechanisms from figure 4. The wave panel shows flux from inertial oscillations advecting geostrophic motions and the inertial panel shows flux from internal gravity waves advecting geostrophic motions.

Early et al. supplementary movie

Change in energy of inertial, wave and geostrophic components as a function of mode and horizontal scale during the Cyprus eddy simulation. Contours on the wave plot indicate the frequency of oscillation, in units of f. The change in energy is dominated by a sustained loss of energy in the j = 1 inertial mode, and a significant gain in wave energy at scales around 250 km. Changes in geostrophic energy are an order of magnitude smaller (note the log color scale), and rapidly oscillate between signs with no sustained gain or loss.

Download Early et al. supplementary movie(Video)
Video 8 MB
You have Access
Open access
Cited by