Fast and Accurate Computation of Vertical Modes

Abstract

The vertical modes of linearized equations of motion are widely used by the oceanographic community in numerous theoretical and observational contexts. However, the standard approach for solving the generalized eigenvalue problem using second‐order finite difference matrices produces $O(1)$ errors for all but the few lowest modes, and increasing resolution quickly becomes too slow as the computational complexity of eigenvalue algorithms increases as $O(n^3)$. Existing methods are therefore inadequate for computing a full spectrum of internal waves, such as needed for initializing a numerical model with a full internal wave spectrum. Here we show that rewriting the eigenvalue problem in stretched coordinates and projecting onto Chebyshev polynomials results in substantially more accurate modes than finite differencing at a fraction of the computational cost. We also compute the surface quasigeostrophic modes using the same methods. All spectral and finite difference algorithms are made available in a suite of Matlab classes that have been validated against known analytical solutions in constant and exponential stratification.

1 Introduction

Vertical modes arise as part of the separable solution to both the internal wave problem and quasigeostrophic theory. The eigenvalue problem (EVP) is treated in many introductory physical oceanography textbooks, (e.g., Cushman‐Roisin & Beckers, ; Gill, ), and the resulting vertical modes describe the vertical structure of the linear solutions for a given density profile. While there are an infinite number of bases that can be used to represent ocean currents and density anomalies satisfying certain boundary conditions, the vertical modes correspond to $O(1)$ dynamical solutions of the equations of motion and are therefore both diagnostic and prognostic. For this reason, vertical modes are the standard basis with which to represent the vertical structure of ocean currents, and it would be hard to overstate their usefulness for describing and modeling the ocean.

There are two primary uses for the vertical modes: (1) projecting a given flow field onto the vertical modes to determine its spectrum (a forward transformation) or (2) creating a dynamically consistent linear flow field from a given spectrum (an inverse transformation). For example, projection onto the vertical modes (the forward transformation) was used to construct the Garrett‐Munk internal wave spectrum (Munk, ; Polzin & Lvov, ), analyze regional and global ocean simulations (Buijsman et al., ; Kelly et al., ; Zilberman et al., ), interpret satellite altimetry observations (LaCasce & Wang, ; Zhao, ), and to describe the vertical structure of balanced flow in the ocean (Wortham & Wunsch, ; Wunsch, ). Examples of using the inverse transform are less common but typically involve initializing a full spectrum of internal waves (Sun & Kunze, ; Winters & D'Asaro, ).

This study was motivated by two situations where current techniques for computing vertical modes were found to have significant errors for reasonable computation times. In the first, perhaps less common situation, a numerical model needed to be initialized with a full spectrum of internal waves—a task which requires solving an EVP at each resolved wavenumber in the model. In the second situation, we sought to compare an observed horizontal velocity spectrum of internal waves to the Garrett‐Munk spectrum near a very strong pycnocline. This requires computing vertical modes at high frequencies—which requires appropriately resolving the mode variability near the pycnocline.

There are a number of sources of error that arise in performing either the forward or inverse transformation. These include the following:

1. A poorly defined mean density function.
2. Measurement noise and uncertainty in the density function and in the case of a forward transform the dynamical variables.
3. Aliasing error, due to the location of the grid points of the dynamical variables (relevant only for forward transform).
4. Interpolation error, due to the location of (and lack of) data points specifying the density function.
5. Numerical truncation error of the modes generated from the vertical EVP.

Whether performing a forward or inverse transformation, these sources of error will compound in some fashion to create error in either the resulting flow field or the resulting spectrum.

The first source of error, a poorly defined mean density function, can result from lack of data or often just uncertainty in what qualifies as a mean (e.g., questions of what time and length scales to average over or even whether averaging is the correct approach for a nonlinear system). Measurement noise and uncertainty in data is a topic unto itself, so the methods here proceed without concern for measurement noise when possible. The aliasing error arises when the grid points are placed such that higher modes project onto the lower modes. This source of error is relatively easily controlled by computing the condition number of the projection matrix, which provides a fairly precise cutoff for the number of resolvable modes. On the rare occasion that a density function can be specified analytically, the interpolation error can be minimized to numerical precision. However, in the usual case where the density is given on some grid with uneven spacing, the density function must be interpolated in between those grid points. Our work suggests that density interpolation error does not dominate the error for most cases. This manuscript is therefore largely concerned with the final source of errors—arising from the numerical representation of the modes in the vertical EVP.

The standard method for solving the vertical EVP is to discretize the problem and construct the matrices using second‐order finite difference matrices (e.g., Cushman‐Roisin & Beckers, ). However, this approach produces unacceptable errors for all but the very lowest modes in the simplest stratifications. This is problematic because numerical algorithms for solving EVPs scale as $O(n^3)$, where $n$ is the number of discretization points in the vertical, so small increases in resolution come at a large computational cost. Instead of using finite difference methods, Kelly () solves the hydrostatic form of the EVP spectrally using Galerkin's method, at a fraction of the computational cost and with much higher accuracy. Dunphy () uses a Chebyshev collocation method to solve the nonhydrostatic case with fixed frequency. Here, we extend the ideas of Kelly () and Dunphy () to the more general nonhydrostatic cases where the EVP must be solved for each frequency in the spectrum or each resolved wavenumber in a numerical model. Our approach solves the EVP spectrally with Chebyshev polynomials to produce high‐quality vertical modes, even with relatively small $n$. The same techniques are then applied to solve for surface quasigeostrophic modes, used to describe the effect of density anomalies at the ocean boundaries.

To illustrate the differences in performance of these two methods, we consider the 10th vertical mode at two different frequencies in a realistic stratification profile. The left panel in Figure shows an example stratification profile taken from the Sargasso Sea during the Latmix 2011 summer field experiments (Shcherbina et al., ). The 10th vertical mode is computed with 512 grid points using both second‐order finite difference methods (blue) and spectral methods (purple). The center panel shows that, at hydrostatic frequencies, finite difference methods perform well. However, the right panel shows the 10th mode at a nonhydrostatic frequency (corresponding to the dashed line in the left panel) for which the finite difference method fails. Generally, we find that finite difference methods with $N$ grid points will return about $N/10$ good modes, although this is stratification and grid point location dependent. For many low‐mode low‐frequency studies then, finite difference methods are fully up to the task. However, for nonhydrostatic frequencies or high modes, the methods described below are more suitable.

The left panel shows a stratification profile taken from the Sargasso Sea in the summertime (Shcherbina et al., ). The two vertical lines indicate a nearly hydrostatic frequency (dashed) and nonhydrostatic frequency (dotted). The 10th mode computed using finite difference methods (blue) and spectral methods (purple) is compared at the hydrostatic frequency (center panel) and nonhydrostatic frequency (right panel). Both methods use 512 evenly spaced grid points as input then interpolate to a high‐resolution grid as output.

This paper begins with a derivation of the two most relevant forms of the vertical EVP that arise from the linearized equations of motion in section . This provides the necessary context for orthogonality relations that form the basis of the normalization of the vertical modes, which in turn shows limitations of using certain vertical modes as a basis. Section demonstrates some of the problems associated with finite differencing, while section shows how these can be remedied using spectral methods. Details of the numerical implementation are described in section , and some of the other sources of error are examined in section . Section discusses some best practices and potential pitfalls. Finally, the appendices include the exact analytical solutions for constant and exponential stratification that are employed for unit testing, the WKB (Wentzel‐Kramers‐Brillouin) approximated solution used for the adaptive grid method, and the class inheritance tree of the publicly available Matlab implementation of these methods.

2 Background

The linearized equations of motion for the fluid velocity $u(x,y,z,t)$, $v(x,y,z,t)$, $w(x,y,z,t)$, on the $f$‐plane are [Image Omitted. See PDF] [Image Omitted. See PDF] [Image Omitted. See PDF] [Image Omitted. See PDF] [Image Omitted. See PDF] where $p(x,y,z,t)$ and $\rho (x,y,z,t)$ are the perturbation pressure and density, respectively. These are defined such that the total pressure $p_{\text{tot}}(x,y,z,t)=p(x,y,z,t)+p_0(z)$ and the total density ${\rho }_{\text{tot}}(x,y,z,t)={\rho }_0+\overline{\rho }(z)+\rho (x,y,z,t)$, where ${\partial }_z{p}_0(z)=-g\overline{\rho }(z)$. All variables in the equations of motion are functions of $x$, $y$, $z$ and $t$, except for $\overline{\rho }$ which is strictly a function of $z$. The operator ${\partial }_z$ is understood to reduce to $\frac{d}{dz}$ when applied to univariate functions. We use the usual definition of the buoyancy frequency $N^2(z)\equiv -\frac{g}{{\rho }_0}{\partial }_z\overline{\rho }$.

There are three linearly independent solutions to equations , assuming periodic horizontal boundary conditions and a flat bottom: two wave solutions and the geostrophic solution.

2.1 Wave Solutions

The positive frequency wave solution is given by [Image Omitted. See PDF] where the functions $F(z)$ and $G(z)$ are the eigenfunctions with corresponding eigenvalue $h$, to be discussed in detail below. The frequency is determined through the dispersion relation, [Image Omitted. See PDF] and the negative rotating wave solution is found by flipping the sign on the frequency, $\omega \mapsto -\omega$. In this notation the phase angle of the wave is given by $\alpha ={\tan }^{-1}\left(\frac{l}{k}\right)$ and $K=\sqrt{k^2+l^2}$ is the total horizontal wavenumber. The horizontal phase is given by ${\theta }_{\pm }=kx+ly\pm \omega t$. The eigenvalue $h$ is referred to as the equivalent depth and is related to the wave group velocity, $c_g=\sqrt{gh}$. The value $h$ can be replaced in favor of eigenfrequency $\omega$ using the dispersion relation, but equation  uses both in order to avoid singularities at $K=0$ and for notational compactness. Applying this solution to equations  leads to two coupled equations for the vertical structure functions, [Image Omitted. See PDF] which can be combined into various second‐order EVPs (see section ).

2.2 Geostrophic Solution

The geostrophic solution is given by [Image Omitted. See PDF] where ${\theta }_0=kx+ly$. The solution can also be written entirely in terms of stream function $\psi (x,y,z)=\frac{g}{f_0}\cos {\theta }_{0}F(z)$, where $(u_g,v_g,{\rho }_g)=(-{\partial }_y\psi ,{\partial }_x\psi ,-{\rho }_0{f}_0{\partial }_z\psi /g)$. Satisfying the vertical momentum equation requires that $N^{2}G=-g{\partial }_{z}F$, but unlike the wave solutions, geostrophic solutions already satisfy continuity. For a given wavenumber $(k,l)$ the geostrophic solution is entirely specified by a vertical profile of any one of the variables, from which all others are immediately deduced. For example, $\widehat{\rho }(k,l,z)$ determines $G(z)$, from which $F(z)$ is determined by integration—the thermal wind balance. There is no preferred basis for the geostrophic solution.

Although the scalings that lead to equations  result in the linear geostrophic solution where $w_g=0$, near‐geostrophic theories with a different choice of scalings, such as quasi‐geostrophy (Pedlosky, ), have nonzero vertical velocities $w_g\ne 0$ and therefore require full continuity, that is, that $F=h{\partial }_{z}G$. As with the wave solution, this requirement combined with $N^{2}G=-g{\partial }_{z}F$ results in an EVP, detailed in section .

An eigenbasis constructed with the rigid lid boundary condition ( $w(0)=G(0)=0$) precludes nonzero density anomalies at the surface. One work‐around to this limitation is to further decompose the geostrophic solution into three parts: two parts resulting from the density anomaly at the boundaries and one part from the remaining density anomaly in the interior. Following Lapeyre and Klein (), the idea is then to let [Image Omitted. See PDF] where both the surface and bottom components of the flow are required to have no potential vorticity in the interior, [Image Omitted. See PDF] but account for the density anomaly at the boundaries, for example, [Image Omitted. See PDF]

The resulting modes can be solved for a given wavenumber and are referred to as the surface quasi‐geostrophic (SQG) modes. This methodology has been used to construct the interior velocity field from sea surface height (SSH) and temperature data Wang et al. ().

Smith and Vanneste () formulate a new EVP that results in modes capable of capturing surface density anomalies for quasigeostrophic flows. Taking equation  from Smith and Vanneste () and writing it in the present notation, we have [Image Omitted. See PDF] with surface boundary condition $G(0)=\frac{D}{\alpha }{\partial }_{z}G$, where $\alpha$ is a tunable parameter. Importantly, these modes remain orthogonal, unlike the combined set of SQG modes and interior modes described above.

An alternative to both the SQG and the Smith and Vanneste () approaches is to use the internal wave eigenbasis constructed with the free surface boundary condition (detailed below) which also results in orthogonal modes capable of capturing nonzero density anomaly at the ocean surface.

2.3 Vertical EVP

The vertical EVP is formed using the two coupled equations from section for the vertical structure functions, $(N^2-{\omega }^2)G=-g{\partial }_{z}F$ (vertical momentum) and $F=h{\partial }_{z}G$ (continuity). In combination with the dispersion relation, one of the eigenfunctions can be eliminated, resulting in various single second‐order EVPs for the vertical structure functions $F$ or $G$. The two most practical EVPs to solve are for $G(z)$ with constant $\omega$, [Image Omitted. See PDF] and for $G(z)$ with constant $K$, [Image Omitted. See PDF] The near‐geostrophic EVP is found by combining the vertical momentum condition $N^{2}G=-g{\partial }_{z}F$ with the continuity condition $F=h{\partial }_{z}G$, which can be treated as the $(\omega \text{‐constant EVP})$ with $\omega =0$. Note that this is equivalent to making the hydrostatic approximation (Kelly, ) and removes all dependence on frequency $\omega$ and wavenumber $K$.

For the cases considered here we take the lower boundary condition at $z=-D$ to be either free slip, where $w(-D)=0$, or no slip, where $u(-D)=0$. These correspond to $G(-D)=0$ and $F(-D)=0$, respectively. These conditions can be seen as the limiting cases of sloped bottom topography (LaCasce, ). The surface boundary condition at $z=0$ is taken to be either a rigid lid with $w(0)=0$, $G(0)=0$ or a free surface approximated as $p(x,y,0)={\rho }_{0}g\eta (x,y,0)$, where $\eta \equiv -{({\partial }_z\overline{\rho })}^{-1}\rho$ is the linear approximation of the isopycnal displacement. In terms of the vertical modes, the free surface boundary condition becomes $h{\partial }_{z}G(0)=G(0)$. Finally, there are many cases where the density profile does not extend to the full depth of the ocean and no boundary conditions (beyond the EVP itself) should be added.

Solving either EVP yields a set of eigenvalues $h_j$ that can be ordered such that $h_1>h_2>h_3>..>h_n$, each with corresponding eigenfunction $G_j$. This means that solving the $(\omega \text{‐constant EVP})$ results in wave solutions with the same frequency $\omega$, but different wavenumbers $K_j$, and similarly solving the $(K\text{‐constant EVP})$ results in wave solutions with the same wavenumber $K$ and different frequencies ${\omega }_j$. Although we do not implement this numerically, note that rearranging the $(K\text{‐constant EVP})$ poses the EVP for fixed group velocity, $gh$, with eigenvalue $K_j^2$.

The equations for the $(\omega \text{‐constant EVP})$ and $(K\text{‐constant EVP})$ are both Sturm‐Liouville problems and share the property that their eigenmodes are orthogonal. Following the procedure in Kelly (), the eigenmodes found with the $(\omega \text{‐constant EVP})$ satisfy [Image Omitted. See PDF] and [Image Omitted. See PDF] while the eigenmodes found with $(K\text{‐constant EVP})$ satisfy [Image Omitted. See PDF] and [Image Omitted. See PDF] where $\beta$ and $\gamma$ are unspecified constants that depend on the chosen normalization, as discussed below. It is important to note that these orthogonality conditions only apply for a particular choice of $\omega$ or $K$. For example, an eigenmode $G_j(z,k_1)$ found using $K=k_1$ is not orthogonal to an eigenmode $G_j(z,k_2)$ found using $K=k_2$ if $k_1\ne k_2$.

The most significant difference between the two EVPs is that eigenmodes from the $(K\text{‐constant EVP})$ often form a complete basis set for typical ocean stratification profiles, while the eigenmodes from the $(\omega \text{‐constant EVP})$ do not. The $(K\text{‐constant EVP})$ is a regular Sturm‐Liouville problem when the weighting function $w_K(z)\equiv N^2-f_0^2>0$ for all $z$, a condition typically met in the ocean. We note that although it is fair to say that stratification with $N>f_0$ is typical of the world oceans, after examining 30,000 CTD profiles, Kunze () finds that 10% of the data have $N<2f$ and a full 30% of the data suggest $N<2f$ within 380 m of the bottom. In contrast, the weighting function $w_{\omega }(z)\equiv N^2-{\omega }^2$ in the $(\omega \text{‐constant EVP})$ switches sign at turning points $z_T$, where $N(z_T)=\omega$. Consequently, the norm of an arbitrary function defined on the domain $[-D,0]$ and satisfying the boundary conditions is not guaranteed to be positive using the norm implied by equation , a necessary condition for completeness. Intuitively, this can be seen in Figure , which shows that the high‐frequency modes have no variance beyond the turning points and are therefore incapable of representing arbitrary functions on the domain.

2.4 Normalization

The amplitude of each vertical mode can be scaled by an arbitrary constant, so one must choose a normalization appropriate for the problem being considered. The four most common scenarios are setting the total energy, a horizontal velocity ( $U$), a vertical velocity ( $W$), and the SSH of a given wave.

To set the total energy of the internal wave solution in equation , we use the modes from the $(K\text{‐constant EVP})$ and therefore use the norm implied by equation , [Image Omitted. See PDF] where we have chosen $\gamma =1$ as the normalization constant in order to keep the vertical modes unitless. Another reasonable choice is to take $\gamma =N_0^{2}D/g$, but using $\gamma =1$ keeps the norm universal, rather than problem specific. Taking the total energy of the wave, [Image Omitted. See PDF] and then depth integrating and averaging over space and time produces a wave with energy $P^2/2$ if we set the coefficient $A=P/\sqrt{h_j-({\omega }^2-f_0^2)G^2(0)/{\omega }^2}$ in equation .

Setting the maximum initial eastward velocity to $U$ can be accomplished by imposing $\max F_j=1$ and $A=U$. The maximum vertical velocity $W$ is set using the norm $\max G_j=1$ where $A=W/(K{h}_j)$, but is clearly singular for inertial waves which have no vertical velocity. The SSH is set using the pressure at the surface with $F(0)=1$ and $A=\text{SSH}·\frac{\omega }{Kh}$.

To set the total energy of the interior geostrophic solution in equation , we assume the solution uses the typical geostrophic modes from equation $(\omega \text{‐constant EVP})$ with $\omega =0$ and therefore use the norm implied by equation , [Image Omitted. See PDF] where we have taken $\beta =D/h_i$. This produces a mode with energy $P^2$ if we let, [Image Omitted. See PDF]

Setting the maximum eastward velocity requires $B=U\frac{f_0}{gl}$ using $\max F_j=1$. The SSH is set using the pressure at the surface by setting $F(0)=1$ and $B=\text{SSH}$.

3 The Problem With Finite Differencing

Computing the lowest vertical modes with finite differencing methods can be relatively fast and accurate when considering a single wavenumber or frequency. Although one can encounter problems with the higher modes, this can usually be ameliorated by increasing resolution. The primary issues with finite differencing arise when needing to compute many modes across a broad range of frequencies and wavenumbers—the two scenarios that motivated the present study. To compute a complete internal wave frequency spectrum requires solving the $(\omega \text{‐constant EVP})$ at each resolved frequency between the Coriolis frequency and the maximum buoyancy frequency, roughly $O(10^2)$ EVPs. This is especially challenging near the buoyancy frequency, where all oscillations occur in the narrow region where $N(z)<\omega$. On the other hand, initializing a numerical model with an internal wave spectrum involves solving the $(K\text{‐constant EVP})$ for each resolved wavenumber in the model, which easily requires $O(10^4)$ computations or more.

A prerequisite to initializing a numerical model with a precise internal wave spectrum is that the modes must be computed for each unique horizontal wavenumber $K$ resolved by the model using the $(K\text{‐constant EVP})$. If the numerical model has $(N_x,N_y)$ horizontal grid points, approximately $N_x{N}_y/2$ unique EVPs must be solved (up to another factor of 2 can be eliminated with isotropic horizontal resolution). Unfortunately, eigenvalue algorithms scale as $O(n^3)$ for $n$ by $n$ matrices. This means that initialization of an internal wave spectrum scales as $O(N_x{N}_y{N}_z^3),$ and thus, with any reasonable vertical resolution, this will quickly become a rate limiting step to a model run. In practical terms, the computation time of these EVPs takes approximately 1, 10, and 100 s for $N_z$ of 512, 1,024, and 2,048, respectively, on consumer hardware from 2015.

The problem is further exacerbated by the poor performance of finite difference methods. To demonstrate, we compare different numerical methods against an analytical solution. Consider an exponential density profile—the canonical deep ocean stratification profile which has known analytical solutions for both the nonhydrostatic internal modes (Garrett & Munk, ), as well as the SQG modes (LaCasce, ). Using a numerical method to solve the $(K\text{‐constant EVP})$, we can compute the relative error of the numerical solution to the analytical solution. We define the relative error as [Image Omitted. See PDF] where $f(z_i)$ is the true solution evaluated at the grid points $z_i$ and $f_i$ is the numerical approximation. Figure shows the maximum relative error of the eigenmodes $F$, $G$ and eigenvalue $h$ found by solving the $(K\text{‐constant EVP})$ using second‐order finite difference methods for standard exponential stratification with 64 vertical grid points (blue). Details of the numerical implementation of the analytical solutions are given in Appendix .

Reasonable error magnitudes are $O(10^{-2})$ ; however, the top panel of Figure shows that no modes computed with finite differencing satisfy this condition. The situation is even worse for the $K=\frac{2\pi }{500\mathrm{m}}$ case shown in the bottom panel, where even the lowest mode has an $O(0.1)$ error. The blue dotted line in Figure shows that increasing the resolution tenfold for second‐order finite differencing decreases the error by a factor of 100, as would be expected. However, this comes at 1,000 times the computational cost and still barely produces any usable modes.

Relative error as a function of vertical mode number using 64 evenly spaced grid points in exponential stratification N(z)=N0ez/b, where N0=3 cph and b=1,300 m at latitude 33° in a 5,000 m deep ocean for K=0.0 (top panel) and K=2π500m (bottom panel). Shown are two second‐order finite differencing methods with (1) 64 grid points (blue), (2) 631 grid points (blue dotted), and 3 spectral methods using Chebyshev polynomials with coordinates in (3) depth (red), (4) WKB scaled (purple), and (5) density (orange).

The accuracy of finite differencing can be increased by going to higher orders (Fornberg, ), since the truncation error at order $s$ scales as ${(\mathrm{\Delta }z)}^s$. The truncation error of the 10th mode in exponential stratification is shown in the top panel of Figure where the blue (solid, dotted) line shows the (second and sixth)‐order finite difference method converging at its predicted rate. The bottom panel shows that even with 1,024 grid points, there are only eight usable modes for the second‐order finite differencing method, while sixth‐order gives up to 50 modes. However, while increasing the order of the method does provide some gains in accuracy, the most efficient way to proceed is simply to use spectral methods, which promise exponentially decreasing truncation error, rather than the polynomial truncation errors offered by finite differencing. When using an analytical density function (dashed lines, Figure ) rather than gridded data, there is no interpolation error and the spectral methods truncation errors reach a noise floor somewhere between 64 and 128 grid points. Furthermore, the number of usable modes is an order of magnitude higher than even the sixth‐order finite difference method. In practical terms, the second‐order finite difference method is producing about 10 good modes in 100 s, while the spectral methods are producing about 100 good modes in 1 s. The increase in truncation error at higher resolution is likely due to increasingly compounded errors of the eigenvalue solvers.

The top panel shows relative error as a function of resolution for the 10th mode in exponential stratification with K=2π500m. The density function is specified on an evenly spaced grid (solid lines) or passed directly as an analytical function (dashed lines). The bottom panel shows the number of usable modes as a function of resolution, defined as the number of modes with truncation errors less than 10−2. The convergence rates of the second‐order and sixth‐order finite difference methods are found to be (Δz)2.0 and (Δz)5.8, respectively.

4 Chebyshev Polynomials

Written in matrix form the $(K\text{‐constant EVP})$ is [Image Omitted. See PDF] where $\mathbf{v}$ is the vector representation of the normal mode $G$ at grid points in $z$, $\mathsf{A}={\partial }_{zz}-K^2$ and $\mathsf{B}=(f_0^2-N^2)/g$. For finite differencing, $\mathsf{A}={\mathsf{D}}_{zz}-K^2\mathsf{I}$, where ${\mathsf{D}}_{zz}$ is the $N_z\times N_z$ differentiation matrix and $\mathsf{I}$ is the $N_z$ dimensional identity matrix. To use Chebyshev polynomials, we project vector $\mathbf{v}$ onto a Chebyshev basis using $\widehat{\mathbf{v}}={\mathsf{T}}^{-1}\mathbf{v}$, where $\mathsf{T}$ is the matrix that transforms a vector from a Chebyshev basis to the coordinate basis. In a practical sense, the columns of $\mathsf{T}$ are the Chebyshev polynomials. Then the EVP becomes [Image Omitted. See PDF] or simply [Image Omitted. See PDF]

The vector $\widehat{\mathbf{v}}$ contains the coefficients needed to reconstruct eigenfunctions, and ${\partial }_{zz}\mathsf{T}$ are the second derivatives of the Chebyshev poynomials.

The optimal choice of grid for Chebyshev polynomials is a Gauss‐Lobatto grid (Boyd, ; Canuto et al., ); for example, equation  below and thus the eigenmatrices and eigenfunctions are always created on a Gauss‐Lobatto grid for any chosen coordinate. Because the basis functions are continuous functions of $z$, the resulting vertical modes can be interpolated onto any grid at any resolution by evaluating the functions at the points of interest.

It is, however, rarely the case that density is given as an analytical function or that observations are made on a Gauss‐Lobatto grid, which means that typically the density needs to be interpolated on the appropriate grid. Interpolation is performed using B‐splines implemented with the numerical framework described in Early and Sykulski (). The advantage to using B‐splines to represent gridded density data is that it is easy to accommodate arbitrary grids. Despite being a low‐order method, this is generally not a limitation (see section ). In the cases shown in Figure , the algorithms are given the density ( $\rho$) on a uniform grid in $z$ of 64 points and the resulting modes are returned on the same grid (except where noted for the high‐resolution finite differencing case). This is, of course, suboptimal for the spectral cases which use a Gauss‐Lobatto grid on various coordinates to compute the EVP. When given an analytical function for density, these methods perform even better, as can be seen in Figure .

Despite the potential limitations imposed by interpolating the density with B‐splines onto a Gauss‐Lobatto grid, the red line in Figure shows that the Chebyshev method performs extremely well, even while interpolating from an evenly spaced grid and outputting to the same grid. The first 20 and 14 modes have error less than $O(10^{-2})$ for the $K=0$ and $K=2\pi /500{\mathrm{m}}^{-1}$ cases, respectively. However, at higher horizontal resolution (larger wavenumbers $K$), even the spectral method's errors grow large, because the points at which the functions are evaluated do not sufficiently capture the oscillations of the modes. This can be remedied by using a stretched coordinate, $s$.

4.1 Stretched Coordinates

In order to find an independent coordinate better suited to capturing the structure of the eigenmodes, we rewrite the EVPs in terms of a generic coordinate $s$ and then consider two concrete examples. Taking $z$ to be a function of $s$ and applying the chain rule leads to [Image Omitted. See PDF] and [Image Omitted. See PDF]

For example, the $(K\text{‐constant EVP})$ becomes [Image Omitted. See PDF] where now $F=h{({\partial }_{s}z)}^{-1}{\partial }_{s}G$. The free surface boundary condition in these coordinates becomes ${({\partial }_{s}z)}^{-1}{\partial }_{s}G=G/h$, while the normalization conditions are now [Image Omitted. See PDF] and [Image Omitted. See PDF]

A necessary condition for using stretched coordinates is that the function $s(z)$ must be strictly monotonic.

4.2 Density Coordinates

For density coordinates $s=-g\overline{\rho }/{\rho }_0$, equation  can be written as [Image Omitted. See PDF] where $F=h{N}^2{\partial }_{s}G$. Note that the derivatives of the density are computed on the $z$ coordinate then projected onto the $s$ coordinate in the EVP. This avoids using inverses of functions that tend towards zero and therefore has greater numerical stability. While the method does well for the high wavenumber case (Figure , lower panel), it performs somewhat poorly with a uniform relative error of $O(0.03)$ for all modes in the low wavenumber (upper panel), as shown by the orange line. Evidently, density coordinates cluster points inefficiently in this case.

4.3 WKB Stretched Coordinates

A compromise between depth ( $z$) and density ( $\overline{\rho }$) coordinates is the WKB stretched coordinate, $s={\int }_D^{z}N(z^{\prime })d{z}^{\prime }$. In this case the EVP becomes [Image Omitted. See PDF] where $F=hN{\partial }_{s}G$.

The purple line in Figure shows that the vertical modes computed on WKB coordinates have uniform accuracy of $O(10^{-3})$ for $K=0$, outperforming the density coordinate case and also performing nearly as well as density coordinates in the high wavenumber case.

4.4 Adaptive Grid for $\omega$‐Constant EVP

Solving the $(\omega \text{‐constant EVP})$ suffers from the additional challenge that as the frequency increases and the distance between turning points decreases, the grid spacing necessary to capture the mode structure becomes ever smaller. As noted in section , this issue arises when considering internal waves near a pycnocline. An example stratification profile and vertical mode found at a frequency approaching the maximum frequency in the pycnocline is shown in Figure . The relative error as a function of mode for this example is shown in Figure , from which it is clear that even the WKB stretched coordinate that performed so well for the $(K\text{‐constant EVP})$ does relatively poorly in this scenario. Sturm‐Liouville theory tells us that the $n$th $F$ mode will have $n$ zero crossings in the oscillatory region where $N^2(z)>{\omega }^2$ (Arfken, ). Thus, in order to resolve these oscillations, one would require at least $2n$ optimally placed points in that region, as well as additional points to capture the variance in the decay regions. Simply increasing resolution of the Chebyshev grid cannot efficiently solve the problem, as grid points will continue to be poorly placed.

The left panel shows a stratification profile with pycnocline taken from Cushman‐Roisin and Beckers (). The vertical dashed line represents a frequency with two turning points in the pycnocline. The right panel shows the fourth vertical F mode at that frequency.

Relative error as a function of vertical mode number using 64 evenly spaced grid points for the frequency and stratification shown in Figure . Shown is the WKB scaled spectral method (purple) as in Figure and also the adaptive grid method (green) that clusters points near the regions of oscillation. Only the first 16 modes are shown.

To resolve this issue, we devise an ad hoc method for clustering points in regions of interest. Our approach is to partition the domain into regions where the modes are hypothesized to be nonzero, formulate the EVP for each region (using WKB stretched coordinates), then couple the equations at the region boundaries. This enables us to assign most of the points to the regions where the solution is assumed to be nonzero, and allocate a few remaining points in the other regions. A comparison of this adaptive method and the standard WKB stretched coordinates is shown in Figure , where the adaptive method is able to capture a few more usable modes than the standard WKB stretched coordinate method with one EVP. The value of this method becomes more pronounced as the maximum frequency is reached. The number of usable modes (error $<10^{-2}$) drops to zero as the maximum buoyancy frequency is approached when using the single EVP, as shown in Figure . However, using the adaptive grid algorithm, we are able to guarantee a minimum number of usable modes as points cluster around the turning frequencies.

The equation boundaries are established by using the WKB approximated solution to identify the regions where the modes are expected to be nonzero. Specifically, the equation boundaries are the points where WKB solution decays to $10^{-5}$ of its value from the turning point. This is an adjustable tolerance, chosen to be small enough that only a few points are needed to capture the nearly zero function but large enough that the nonzero regions aren't unnecessarily large. The gray vertical lines in Figure show the number of coupled equations being used to solve the EVP. At the lowest frequencies only one equation is used and the method is identical to the WKB stretched coordinate method described in section . As the frequency increases, the algorithm eventually separates into two coupled equations: one for the top boundary and pycnocline and another for the deep region where no mode variance is expected (refer to Figure ). At high enough frequency the region above the pycnocline is decoupled as well, and three coupled EVP problems are solved.

The EVPs are coupled by requiring that the function and its first derivative are continuous at the equation boundaries, following the procedure described in section 22.3 of Boyd (). The “eigenvalue rule‐of‐thumb” as discussed in Boyd () states that roughly $n/2$ modes will be accurate when using $n+1$ Chebyshev polynomials away from boundary layers or critical levels. Solving the EVP with turning points near the maximum buoyancy certainly does not satisfy this criterion, but the rule‐of‐thumb can be modified to use half of the modes with eigenvalues greater than zero. Although we make no attempt at proving the general validity of this modification, the dashed line in Figure indicates that the rule‐of‐thumb generally does well at predicting how many modes are good quality.

The number of usable modes (error <10−2) versus frequency for the stratification in Figure using 64 points. The adaptive grid method and standard WKB method are shown in green and purple, respectively. The dashed green line is the rule‐of‐thumb number of good modes by the adaptive grid method. The two vertical gray lines separate the regions where the adaptive algorithm used 1, 2, or 3 coupled equations.

5 Numerical Implementation

One of the primary products of this paper is the implementation of these methods as classes in Matlab (see Appendix  for more details). Figure shows the flowchart followed by the initialization algorithm for the InternalModes class, described in this section.

Algorithm flowchart for the initialization of InternalModes class.

The two methods for initializing the classes are both called using

• im = InternalModes(rho,z,zOut,latitude);
•  

where the arguments rho,z are either a gridded density field at locations z or function handle valid in the domain spanned by [min(z) max(z)]. When the function handle is given, the density function is projected onto Chebyshev polynomials. If gridded data are provided, then the density is interpolated using B‐splines. The argument zOut specifies the grid on which all output is given, which need not span the full depth.

After initialization, all classes support setting the upper and lower boundary conditions as well as setting the normalization to any of the choices discussed in section .

The two primary functions for computing internal modes are

• [F,G,h,omega] = im.ModesAtWavenumber(k);
•  

for the $(K\text{‐constant EVP})$ and

• [F,G,h,k] = im.ModesAtFrequency(omega);
•  

for the $(\omega \text{‐constant EVP})$.

The implementation of these methods for finite differencing is straightforward—the EVP is solved either on the gridded input data as given or on a grid that matches the output grid if specified as a function. The differentiation matrices are created using the algorithms described in Fornberg (). However, the spectral implementations require additional choices.

The EVP being solved is [Image Omitted. See PDF] where $s$ is a generic stretched coordinate, $a(s)$, $b(s)$, $c(s)$, and $d(s)$ are referred to here as coefficient functions.

The algorithm can be separated into the three parts. First, we compute the coefficient functions for each EVP, for example, $N^2$ and ${\partial }_{z}N$ for equation . Second, the EVP is solved on the appropriate coordinate with $n_{\text{evp}}$ points. Finally, the resulting modes are normalized and projected onto the output grid with an arbitrary number of points.

5.1 Initialization With an Analytical Function

If the InternalModes class is initialized with a function handle for the density, it is projected onto Chebyshev polynomials which are then used to compute the coefficient functions and, if necessary, the stretched coordinate.

To project onto Chebyshev polynomials, we define a grid with $n_{\text{grid}}$ points using [Image Omitted. See PDF] where $i$ is an integer index ranging from 0 to $n_{\text{grid}}-1$. We evaluate the density function on that grid, $\overline{\rho }(z_{\text{lobatto}}^i)$. The density function is then expanded in a Chebyshev polynomial basis such that [Image Omitted. See PDF] where ${\widehat{\rho }}^k$ indicates the $k$th coefficient for Chebyshev polynomial defined on $z$ coordinates. The coefficients for the derivative of the function, denoted ${\widehat{\rho }}_z^k$, are then computed using a recursion formula [Image Omitted. See PDF] where $c_k=2$ for $k=0$ and $c_k=1$ otherwise. Because a Gauss‐Lobatto grid was used for the $z$‐coordinate, the Chebyshev transformation is performed with a rescaled fast Fourier transformation in $O(n_{\text{grid}}\log n_{\text{grid}})$ operations. The differentiation requires only $O(n_{\text{grid}})$ operations, which means that all of the coefficient functions for the EVP can be computed on a relatively fine grid. For example, $n_{\text{grid}}=O(2^{14})$ takes a fraction of a second on consumer hardware from 2015.

The stretched coordinates implemented here are either $s=z$, $s=-g\overline{\rho }/{\rho }_0$, or $s={\int }_D^{z}N(z^{\prime })\mathrm{d}{z}^{\prime }$, where the latter two cases require density to be strictly monotonic. For those two cases if ${\overline{\rho }}_z\ge 0$ anywhere in the domain, then an error is thrown. For the WKB coordinate, $s={\int }_D^{z}N(z^{\prime })\mathrm{d}{z}^{\prime }$, the integral is computed spectrally using equation .

5.2 Initialization With Gridded Data

In the more typical scenario where a user initializes the InternalModes class with gridded data from observations or a numerical model, B‐splines are used to interpolate the data and compute the coefficient functions. The primary advantage to using B‐splines in this scenario is that B‐splines can be created for arbitrary grids without suffering from Runge's phenomena at lower orders. We fit the data to sixth‐order interpolating spline (with five nonzero derivatives) using the methodology and numerical implementation described in Early and Sykulski ().

If the method requires that density remain monotonic (e.g., for WKB and density coordinates), then the B‐spline fits are constrained to be monotonic following Pya and Wood (). If the data are not monotonic, then this implicitly smooths to find the nearest monotonic fit in a least squares sense.

Computing the stretched coordinate and the coefficient functions from the spline interpolant requires derivatives and integrals of the B‐splines, which are relatively straightforward to compute because they are just piecewise polynomials (De Boor, ). The WKB method requires computing the square root of the B‐spline interpolant. The approach taken here is to build a new interpolating spline of the same order that interpolates between the square root of the data points.

5.3 EVP

All three Chebyshev methods solve their respective EVPs on a Gauss‐Lobatto grid of their respective coordinate, that is, $s=z$, $s=-g\overline{\rho }/{\rho }_0$, or $s={\int }_D^{z}N(z^{\prime })\mathrm{d}{z}^{\prime }$. The Gauss‐Lobatto grid in $s$ is defined in the usual way as [Image Omitted. See PDF] where the number of points $n_{\text{evp}}$ reflects the size of the EVP and is therefore also an absolute upper bound to the number of modes that may be computed.

Once the Gauss‐Lobatto grid $s_{\text{lobatto}}^i$ is created, the corresponding value $z(s_{\text{lobatto}}^i)$ is computed using the bisection method as implemented in Driscoll et al. (). The method is set to terminate with a relative error of $O(10^{-12})$.

The coefficient functions in equation  are now simply evaluated onto the $s_{\text{lobatto}}^i$ grid using the interpolant (either a B‐spline or Chebyshev). The Chebyshev polynomials and their derivatives ( $T,T_s$, and $T_{ss}$) are computed using the standard recursion formulas and then multiplied by the coefficient functions to create eigenvalue matrices. The boundary conditions are implemented by replacing the first and last rows of the matrices $\mathsf{A}$ and $\mathsf{B}$ in .

The EVP is then solved using the standard generalized EVP solver. This is typically the rate limiting step in the process, taking $O(n_{\text{evp}}^3)$ operations. The resulting eigenvectors now contain coefficients to the Chebyshev polynomials defined on the stretched coordinate.

5.4 Adaptive Grid

The adaptive grid is created by locating regions in the domain where the solution is expected to be small and then allocating fewer points (and therefore fewer Chebyshev polynomials) to those regions. After identifying the turning points $z_T$ where $N^2(z_T)={\omega }^2$, the WKB solution is used to identify the equation boundaries $z_{\text{bnd}}$, the points where $F_j^{\text{WKB}}(z_{\text{bnd}},\omega )/F_j^{\text{WKB}}(z_T,\omega )=10^{-5}$. The WKB solution is assumed to work locally in the stratification and is therefore applied at the turning point, $z_T$, in the direction of decaying variance. The eigenvalue is assumed to be set globally by integration over all oscillatory regions. If $m$ equation boundaries $z_{\text{bnd}}$ are identified, they delineate $m+1$ regions: the “decay” regions, where the solution is anticipated to be small and governed by the decaying exponential and “oscillatory” regions where the solution is expected to dominate and include the oscillatory solution. These decaying and oscillating regions are necessarily alternating.

The $m+1$ regions are coupled using the same technique described in section 22.3 of Boyd () by requiring continuity across boundaries for $G$ and ${\partial }_{z}G$. The key benefit to this algorithm is that the decay regions are allocated as few as 6 points each, while oscillatory regions are apportioned the remaining points relative to their WKB length, $L^m=\int N(z^{\prime })\mathrm{d}{z}^{\prime }$. While we find that 6 points appear to be sufficient for the decay regions, in practice, we do in fact apportion 1/16th of the total points evenly between the decay regions as a hedge that this ad hoc method will fail for some unforeseen cases.

The adaptive grid algorithms use low‐order interpolation to identify $z_T$ and $z_{\text{bnd}}$ because high accuracy is not required, and this keeps the number of computations $O(n_{\text{grid}})$.

5.5 Normalization

The final step is to normalize the resulting eigenvectors and project them onto $z_{\text{out}}$. Normalization using the two integral conditions can be performed exactly invoking the fact that the integral of each Chebyshev polynomial $T^k(z)$ is exactly known [Image Omitted. See PDF]

We have defined $w^k$ such that it can be summed with a vector of Chebyshev coefficients to produce the integral. In other words, if ${\widehat{v}}^k$ is a Chebyshev coefficient vector, then $I=\sum {\widehat{v}}^k{w}^k$ is the definite integral. The integrands in equations  and are computed pseudospectrally before integrating (by transforming to the spatial domain, multiplying, then transforming back Chebyshev coefficients).

Computing the max U and and max W norm is more problematic because the function extrema do not necessarily lie on the grid points. For the implementation, we simply take the maximum at the resolved grid points, but if higher accuracy is required, one could locate the extrema using the methods in Boyd ().

Finally, the normalized eigenmodes are projected onto the output grid using the slow Chebyshev transforms, [Image Omitted. See PDF]

If a large number of output points are requested, this operation could dominate the total computation time.

The algorithm flowchart for the mode computation is shown in Figure

Algorithm flowchart for the mode computation.

5.6 SQG Modes

The two functions for computing the SQG modes are

• psi = im.SurfaceModesAtWavenumber(k);
•  

and

• psi = im.BottomModesAtWavenumber(k);
•  

where k is an array of wavenumbers.

The SQG modes are found from equation  using a linear solver after replacing the top and bottom points of the matrix with the boundary conditions. As noted in Tulloch and Smith (), the SQG modes require a high density of grid points near the boundaries, a task well suited to the Gauss‐Lobatto grid in equation . The number of Chebyshev polynomials is chosen so that the Gauss‐Lobatto grid captures at least 10 points over the $e$‐folding scale. The $e$‐fold scale for constant stratification (see Appendix ) is $\mathrm{\Delta }{z}_{\text{efold}}=\frac{f_0}{K{N}_0}$, and the distance between the first two points in a Lobatto grid, equation , is [Image Omitted. See PDF] where we have defined $D=z_{\max }-z_{\min }$ as the depth of the domain. Setting $\mathrm{\Delta }{z}_{\text{boundary}}=\frac{1}{10}\mathrm{\Delta }{z}_{\text{efold}}$, we find that we need [Image Omitted. See PDF] points (and therefore also $n_{\text{grid}}$ polynomials) to sufficiently capture the SQG mode. The resulting SQG modes are projected onto the output grid using equation .

5.7 Unit Testing

In order to ensure that each of the algorithm implementations is correct, the numerically generated eigenmodes and eigenvalues are compared against the analytical solutions for constant stratification and exponential stratification shown in Appendix . The comparison is performed across a range of wavenumbers and frequencies for both surface boundary conditions and all four norms.

The computed SQG modes are also compared against analytical solutions for constant and exponential stratifications, shown in Appendix .

6 Other Sources of Error

In section we noted five sources of errors that contribute to the total error when computing the forward or inverse transformation, but this manuscript has primarily focused on one source of error: the numerics of accurately representing the vertical modes in the EVP. We now discuss these other sources of error and how they are dealt with in the numerical implementation.

6.1 Aliasing Error

When performing a forward transformation, where a given dynamical field is projected onto the vertical modes, the data grid will determine how many modes are resolvable. As with a Fourier transformation, higher‐frequency modes alias into the lower‐frequency modes. Unlike the Fourier transformation, however, the optimal grid for performing a transformation is not an evenly spaced grid but depends on the stratification profile, and therefore, the EVP being solved. Here we show that there is a relatively easy method for determining the number of resolvable modes for a given grid using the condition number of the resulting matrix.

To show the effect of different grid choices on the forward transformation, we use the analytical solution of the vertical modes in exponential stratification in combination with an imposed spectrum to generate stochastic isopycnal displacement profiles typical of the world oceans. In particular, we use the Garrett‐Munk spectrum, [Image Omitted. See PDF] where $j_{*}$ is the roll‐off mode, usually set to 3 but possibly as high as 20, $p$ is the slope which is usually very nearly $5/2$, and $H_0$ normalizes sum over $j=1..\infty$ to unity. For each stochastically generated set of coefficients, $m^j$, a profile $\eta (z)={\sum }_{j=1}^N{G}^j(z)m^j$ is created. The profile is then evaluated on three different grids: $z_{\text{even}}^i$, $z_{\text{lobatto}}^i$, and $z_{\text{quadrature}}^i$, where $z_{\text{quadrature}}^i$ is the grid of Gaussian quadrature points, determined by the roots of a mode one higher than is trying to be used (Boyd, ; Press et al., ). Using the first $n$ modes, we then attempt to recover the coefficients using least squares—in practice this is Matlab's mldivide (\) operator. For example, the first $n$ coefficients of the evenly spaced grid are determined by [Image Omitted. See PDF] where $j=1..n$. The root‐mean‐square error (rmse) is defined as the error in the sum of recovered and missing coefficients [Image Omitted. See PDF]

The top panel shows the root‐mean‐quare error as a function of total modes used to recover the coefficients with an inverse transformation. Vertical dashed lines are the predicted cutoffs for the different grids, based on the matrix condition number at which the modes are no longer resolvable. The condition number as a function of total modes is shown in the bottom panel.

Figure shows the result of trying to recover the mode coefficients, ${\tilde{m}}^j$ using successively more modes for the three different grids. In all three cases the rmse decreases as modes are added until a dramatic increase occurs, correlating with a similarly dramatic increase in the condition number of the matrix. The quadrature grid, defined as the roots of the $G^{N-1}$ mode plus the boundaries, performs best, as expected. Including the boundaries in this definition means that the top and bottom boundary points provide no useful information, and therefore, only $N-2$ modes are recoverable for $G$ and $N-1$ for $F$ with a rigid lid. This definition is chosen so that the forward and inverse transform for constant stratification coincides with the discrete sine and cosine transform (and their associated grids).

While the condition number of the matrix is clearly controlled by the grid being used, the choice of norm also affects the condition number. Generally speaking, the $K\text{‐constant norm}$ performs well for transformations with the $G$ modes and the $\omega \text{‐constant norm}$ with the $F$ modes. The exception to this is when the free surface boundary condition is used, the barotropic mode has a substantially different $L^2$ norm and should be rescaled.

6.2 Mode Error

Here, we test whether or not the truncation error in the vertical modes is a limitation in recovering the coefficients of the spectrum, $m_j$. To this end, we created profiles on a quadrature grid with variance distributed using $H(j)$ for a range of parameters including white noise (large $j_{\star }$), and very smooth (large $p$), as described above, and recorded how many mode coefficients were recoverable for some error tolerance across all the different numerical methods in this paper. We found that discarding modes with truncation errors exceeding the requested error tolerance of the modes guarantee coefficients recovered with the requested error tolerance. In other words, if one wants relative errors in coefficients of less than $10^{-2}$, one needs modes with errors less than $10^{-2}$. The only exception is for very steep spectra (e.g., $p=-10$), where the coefficients of the higher modes are indistinguishable from zero. The reverse of this is not true—in fact, including modes with truncation errors exceeding the error tolerance can often return coefficients within the bounds of the error tolerance. Evidently, the errors in these ordered, orthogonal bases work systematically in our favor.

6.3 Interpolation Error

Another source of error may arise from interpolation of the density function. This issue is treated separately from a poorly defined or noisy mean density function and therefore assumes that the data given are gridded and without error. For gridded data, our method uses a relatively low order B‐spline to interpolate between grid points which is then evaluated for the coefficient functions where needed. Does this low‐order interpolation method limit the mode recoverability described in the previous section?

To address this question, we compute the vertical modes for exponential stratification from an evenly spaced density function with variable number of grid points. These modes are then used to recover the coefficients ${\tilde{m}}^j$, as above, on a high‐resolution quadrature grid. We find that with as few as 16 grid points for the density function, we are able to recover 90 mode coefficients with less than 1% error.

6.4 Poorly Defined or Noisy Mean Density

It is often the case that the mean density function, $\overline{\rho }(z)$, is not easily defined. Averaging over a mooring time series of density will not necessarily result in a monotonic density function—and averaging often removes the sharp gradients that exist in individual profiles, which may not be desirable. Even output from numerical models can suffer from these same issues, depending on the boundary conditions. Noisy data, where errors in the observed value of $\overline{\rho }(z)$ stem from instrument errors, also effectively constitute a poorly defined mean. One way to frame these issues is to ask how a misspecified mean density affects our ability to infer the vertical spectrum of a given flow.

First, we note that all methods, as implemented, will proceed without error for noisy data. However, the most notable difference between methods is that the WKB and density coordinate methods use a density function constrained to the nearest monotonic spline fit, as previously described. It is not clear that this implicit smoothing is necessarily “better” than using the unaltered density function with the $z$‐coordinate method for noisy data. This decision, and how to deal with measurement noise in general, is beyond the scope of this manuscript. Additional smoothing of the density data can be done using many techniques, including the constrained smoothing splines described in Early and Sykulski ().

In order to test the effects of misspecifying a mean density profile, we generate density profiles in exponential stratification that follow the GM spectrum, as above, but we attempt to recover the coefficients using vertical modes computed from a noisy mean density profile. The results are consistent with section . For example, as long as the modes from the noisy profile have errors less than $10^{-2}$ relative to the modes generated from the true profile, the recovered coefficients ${\tilde{m}}_j$ will have errors less than $10^{-2}$ relative to $m_j$. The relative error of the vertical modes increases as a function of mode number until eventually the mode errors and therefore coefficient errors reach $O(1)$. Exactly where mode errors reach $O(1)$ depends on the details of the noise, or misspecification, of the mean density profile.

While accurately recovering the coefficients $m_j$ of the spectrum becomes impossible with a noisy mean density profile, we are able to accurately infer the spectrum from which $m_j$ was generated by either ensemble averaging over additional synthetic profiles or bin averaging nearby modes. One way to see why this might be true is to note that the coefficient errors never exceed $O(1)$—so although the exact coefficient is incorrect, the magnitude is correct on average. In practice, using modes from the misspecified mean density profile causes variance that should be associated with mode $j=33$, for example, to be assigned to the variance of nearby modes.

That the spectrum is recoverable despite a noisy mean is consistent with previous analysis methods. In Polzin and Lvov (), internal wave spectra are found using WKB approximated modes and a WKB stretched grid. It is also important to note that computing the spectrum with a misspecified mean density function still requires an orthogonal set of modes, and therefore, fast and accurate mode computation is still helpful.

7 Discussion

The spectral methods presented here solve the most relevant forms of the vertical eigenvalue and surface quasigeostrophic mode problems efficiently and accurately. The methods also include an algorithm for computing modes in challenging stratification profiles at high frequencies near turning points. The algorithms are implemented in a publicly available Matlab suite (https://github.com/JeffreyEarly/GLOceanKit/) using the class hierarchy described in Appendix , and the implementations are validated against known analytical solutions, under a wide range of conditions. However, the methods do not always perform well under all conditions.

Poorly resolved features and discontinuities in the density profile will produce the Gibbs phenomenon, where “ringing” occurs in the vicinity of the discontinuity. In one example, we found that a narrow 5‐m‐wide pycnocline in a 5,000‐m‐deep ocean produced strong spurious oscillations unless the pycnocline was sufficiently well resolved with enough grid points. In another example, we defined an analytical profile with a discontinuity in ${\overline{\rho }}_{zz}$ (the highest derivative used in the EVP) and this also produced the Gibbs phenomenon. Interestingly, in these cases lower‐order finite differencing produced better modes for the same vertical resolution, because these methods implicitly smooth the derivatives. A logical extension of this work would be to apply the “splitting” algorithms in chebfun (Driscoll et al., ) to handle such discontinuities.

Our recommendations for projecting observed or modeled fields onto the modes are as follows:

• When possible, use a quadrature grid for the fields and modes, for example, the function GaussQuadraturePointsForModesAtWavenumber in the InternalModesSpectral class computes the Gauss quadrature points for $G$. This is a relatively expensive operation (it involves solving the EVP) but provides near‐optimal point placement.
• In the usual case where there is no freedom to choose grid points, compute the condition number as a function of mode as described in section and limit the number of modes to a low condition number. Alternatively, if computation time is not a limitation, one can perform the least squares fit of the fields to successively more modes until the coefficients are no longer stable.

Many of the errors described in section may still be a concern but may be quantified with some of the techniques in section , by using the density function and spectrum specific to the problem.

A AppendixInternal Mode Solutions

We present analytical solutions for the internal mode eigenvalue problem in three different scenarios: constant stratification, exponential stratification, and the WKB approximated solution for arbitrary stratification. These solutions are used to validate the numerical implementations.

A.1 Constant Stratification

The internal baroclinic modes in constant stratification are given as [Image Omitted. See PDF] [Image Omitted. See PDF] with eigendepth $h_j$ and vertical wavenumber $m_j$ given by [Image Omitted. See PDF] where we have assumed that $w=0$ at the lower boundary $z=-D$. In the case of a rigid lid ( $w=0$ at $z=0$), the correction to the vertical wavenumber is $\xi =0$. However, if the linearly approximated free surface boundary condition is used, $h_j{G}_j(0)=G_j(0)$, then $\xi$ is nonzero. The equations for $\xi$ are transcendental and are therefore solved with a numerical root finding algorithm. The equations for $\xi$ are written in a form conducive for finding the desired root.

• For fixed wavenumber, $k$, the vertical wavenumber correction $\xi$ is found by solving [Image Omitted. See PDF] near $\xi =0$ and the eigendepth $h_j$ is given by [Image Omitted. See PDF]
• For fixed frequency, $\omega$, the vertical wavenumber correction $\xi$ is found by solving [Image Omitted. See PDF] near $\xi =0$ and the eigendepth $h_j$ is given by [Image Omitted. See PDF]

A1 Inverse of the Normalization Constants ( $A^{-1}$) for All Combinations of Solutions and Norms Considered Here

	Trigonometric	Linear	Hyperbolic
$K$‐constant	${(-1)}^j\sqrt{{\sin }^2\left(m_{j}D\right)+\frac{(N_0^2-f_0^2)D}{2g}\left(1-\frac{\sin \left(2m_{j}D\right)}{2m_{j}D}\right)}$	$D\sqrt{1+\frac{(N_0^2-f_0^2)D}{3g}}$	$\sqrt{{\sinh }^2\left(m_{0}D\right)+\frac{(N_0^2-f_0^2)D}{2g}\left(\frac{\sinh (2m_{0}D)}{2m_{0}D}-1\right)}$
$\omega$‐constant	${(-1)}^j{h}_j{m}_j\sqrt{\frac{1}{2}+\frac{\sin \left(2m_{j}D\right)}{4m_{j}D}}$	$D$	$h_0{m}_0\sqrt{\frac{\sinh \left(2m_{0}D\right)}{4m_{0}D}+\frac{1}{2}}$
$U_{\max }$	${(-1)}^j{h}_j{m}_j$	$D$	$h_0{m}_0\cosh (m_{0}D)$
$W_{\max }$	${(-1)}^j,\sin (m_{0}D)\text{for}j=0$	$D$	$\sinh (m_{0}D)$

The normalization for these modes is given in the first column of Table .

In the case of the free surface boundary condition, there also exists a barotropic mode ( $j=0$), the solution of which changes from trigonometric to hyperbolic at $\omega =N_0$, or $k=k_{*}$ where [Image Omitted. See PDF]

In these cases then, the vertical wavenumber reduces to $m_0=\xi /D$.

• Trigonometric case, $k<k_{*}$ or $\omega <N_0$

  The solution is exactly the same as the baroclinic solutions given above in equations  and , but now equation  is solved when $j=0$. As a practical matter, it is numerically more stable to solve [Image Omitted. See PDF] for the positive root near $\xi =\sqrt{(N_0^2-f_0^2)D/g-k^2{D}^2}$. For fixed frequency, equation  can be used to find the root near $\xi =\sqrt{D(N_0^2-{\omega }^2)/g}$.

• Linear case, $k=k_{*}$ or $\omega =N_0$

  The solution is given by [Image Omitted. See PDF] [Image Omitted. See PDF] where $h_0=D$.

• Hyperbolic case, $k>k_{*}$ or $\omega >N_0$ The solution is given by [Image Omitted. See PDF] [Image Omitted. See PDF]

  • $-$
  • For fixed wavenumber, $k$, the vertical wavenumber correction $\xi$ is found by solving [Image Omitted. See PDF] for the root near $\xi =\sqrt{k^2{D}^2-(N_0^2-f_0^2)D/g}$ and the eigendepth $h_0$ is given by [Image Omitted. See PDF]
  • $-$
  • For fixed frequency, $\omega$, the vertical wavenumber correction $\xi$ is found by solving [Image Omitted. See PDF] for the root near $\xi =\sqrt{D({\omega }^2-N_0^2)/g}$ and the eigendepth $h_0$ is given by [Image Omitted. See PDF]

A.2 Exponential Stratification

Exponential stratification serves as the $O(1)$ representation of the ocean stratification away from the poles and continental boundaries. As formulated and first solved in Garrett and Munk (), we take the stratification to be $N^2=N_0^2{e}^{2z/b}$, where $N_0$ is buoyancy frequency and $b$ is the thermocline depth (with canonical values of three cycles per hour and 1,300 m).

Letting the stretched coordinate $s=N_0{e}^{z/b}$ as in section , the $(\omega \text{‐constant EVP})$ becomes [Image Omitted. See PDF] which has solution [Image Omitted. See PDF] where [Image Omitted. See PDF] is chosen to satisfy the lower boundary condition $G(-D)=0$. The Bessel function $Y_{\nu }(s)$ has a singularity at $s=0$, so for many choices of $\omega$ and $c_j$, ${\alpha }_j\ll 1$ and the second term needs to be neglected for stable numerical evaluation. The order of the Bessel function is set by the frequency $\nu =\frac{b\omega }{c}$, or, using the dispersion relation wavenumber $\nu =\sqrt{\frac{b^2{f}_0^2}{c^2+b^2{k}^2}}$. The $F$ modes are found by taking the derivative [Image Omitted. See PDF]

The discretization into modes is a result of applying the second boundary condition, in this case either a rigid lid $G(0)=0$ or a free surface $G(0)=F(0)$ and then finding the eigenmode speeds $c_i=\sqrt{g{h}_i}$ that satisfy those conditions. The $c_i$s are therefore determined by the roots of Bessel functions, for which there is no general closed form solution. The challenge then becomes finding bounds for the roots and writing the equation in a form suitable for a root finding algorithm.

In practice, it is easiest to find the inverse of the $c_i$s, so we let $x=\frac{b{N}_0}{c}$ and then write the root equation, $f(x)=0$, in terms of parameter function $\nu (x)$ and $s(x)$. Again, stable numerical evaluation requires that we work around the singularity of $Y_{\nu }$, so if $\nu (x)<s(x)e^{-D/b}$, then we take [Image Omitted. See PDF] and when $\nu (x)>s(x)e^{-D/b}$ we take [Image Omitted. See PDF] where in either case the rigid‐lid condition requires [Image Omitted. See PDF] [Image Omitted. See PDF] and the free surface requires [Image Omitted. See PDF] [Image Omitted. See PDF]

• $\omega$‐constant solution

  The parameter functions are defined as [Image Omitted. See PDF] [Image Omitted. See PDF] where the scaling factor $\eta$ is chosen from the WKB approximated solution (Desaubies, ),

  • $-$

  • $\omega >N_0{e}^{-D/b}$

    The root equation must be taken to be $f_{{\nu }_{\text{big}}}(x)$ and [Image Omitted. See PDF](from Desaubies, , equation 2.18); the first $n$ modes are found within $[0.5,n+1]$.

  • $-$

  • otherwise

    The root equation must be taken to be $f_{{\nu }_{\text{small}}}(x)$ and [Image Omitted. See PDF](from Desaubies, equation 2.19); the first $n$ modes are found within $[0.5,n+1]$.

• $K$‐constant solution

  Unlike the $\omega$‐constant solution, the order of the Bessel function changes for each root, and therefore, the root equation being used may have to change during evaluation.

  The parameter functions are defined as [Image Omitted. See PDF] [Image Omitted. See PDF] where $ϵ=\frac{f_0}{N_0}$ and $\lambda =bk$. Desaubies () provides low‐frequency (lf) and high‐frequency (hf) bounds for the eigenvalues, which we can be written in terms of wavenumber [Image Omitted. See PDF] [Image Omitted. See PDF]

  We transition at [Image Omitted. See PDF]

  To set the search bounds for the root algorithm, we want [Image Omitted. See PDF] and [Image Omitted. See PDF] so the first $n$ modes are found within $[x_{\text{lower}},x_{\text{upper}}]$.

Normalization is performed by direct evaluation of the mode functions for the $U_{\max }$ and $W_{\max }$ and by numerical integration for the $K\text{‐constant norm}$ and $\omega \text{‐constant norm}$ norms.

A.3 WKB Solution

Desaubies () found the WKB solution for stratification with at most one turning point $z_T$. Simplified and using the notation of this manuscript, the WKB solution with single turning point $z_T$ is [Image Omitted. See PDF] [Image Omitted. See PDF] where $\xi =sgn({\omega }^2-N^2){\left(\frac{3}{2}q\right)}^{\frac{2}{3}}$ and we have defined [Image Omitted. See PDF] so that it goes to zero at the turning point but is positive everywhere else. The eigenvalue is proportional to the integral of the stratification over the oscillatory region [Image Omitted. See PDF] while normalization coefficient for the $K\text{‐constant norm}$ is [Image Omitted. See PDF]

It is useful to note that away from the turning point, the Airy function $Ai(\xi )$ can be approximated as sinusoidal above the turning point and exponentially decaying below the turning point [Image Omitted. See PDF]

In the case of no turning point ( $\omega <N(z),\forall z\in [-D,0]$) the WKB solution is [Image Omitted. See PDF] [Image Omitted. See PDF] [Image Omitted. See PDF] where we have taken $z_T=-D$ in . The eigenvalue is now [Image Omitted. See PDF]

Note that this solution does reduce to the exact solution for constant stratification.

B AppendixSQG Mode Solutions

The SQG modes are computed by taking the Fourier transform of $\psi$ in equation  and scaling by the Fourier transformed boundary condition $\widehat{\rho }(k,l,z=0)$ such that $\psi (x,y,z)=\int \int \left({\left.-\frac{g}{{\rho }_0}\widehat{\rho }\right|}_{z=0}\right)\varphi (k,l,z)e^{i(kx+ly)}\mathrm{d}k\mathrm{d}l$. Equation  then becomes [Image Omitted. See PDF] with surface boundary condition $f_0{\partial }_z\varphi =1$ and bottom boundary boundary condition ${\partial }_z\varphi =0$ (and vice versa for the bottom boundary modes).

To solve this problem numerically, we project $\varphi$ onto a Chebyshev basis where $\mathbf{v}=\mathsf{T}\widehat{\mathbf{v}}$ is the vector representation of $\varphi$. The equation for the SQG modes becomes [Image Omitted. See PDF] in the interior with boundary conditions $f_0{\partial }_z\mathsf{T}\widehat{\mathbf{v}}=1$ and $f_0{\partial }_z\mathsf{T}\widehat{\mathbf{v}}=0$ for the surface mode.

The numerical implementation of equation  is validated against the known analytical solutions.

B.1 Constant Stratification

The surface quasigeostrophic modes for constant stratification can be found in Tulloch and Smith (). For the upper boundary this is [Image Omitted. See PDF] and the lower boundary [Image Omitted. See PDF] where $\lambda =\frac{K{N}_0}{f_0}$. These solutions cannot be evaluated numerically for all wavenumbers and depths because the $\sinh$ function may overflow. Instead, we use [Image Omitted. See PDF] and [Image Omitted. See PDF] for reliable numerical evaluation.

B.2 Exponential Stratification

The surface SQG modes for exponential stratification are first solved by LaCasce (). Defining the scale $\eta \equiv \frac{N_{0}Kb}{2f_0}$, the surface mode is given by [Image Omitted. See PDF] while the bottom mode can be shown to be [Image Omitted. See PDF]

Bessel functions $K_0$ and $K_1$ should not be confused with the wavenumber $K$.

B.3 WKB Solution

The WKB solution for the SQG mode does not appear to have been previously derived. Assuming a solution of the form $\varphi =e^{q(z)}$ and inserting this into equation , we have that [Image Omitted. See PDF]

We then assume a series expansion $q=q_0+ϵq_1+..$. combined with the assumption that both $q^{\prime \prime }$ and $N^{\prime }$ are $O(ϵ)$. The two lowest order equations are therefore [Image Omitted. See PDF] [Image Omitted. See PDF] which has solution [Image Omitted. See PDF]

Applying boundary condition $f_0{\varphi }_z(-D)=0$ and $f_0{\varphi }_z(0)=1$, the surface mode is [Image Omitted. See PDF] where $s(z)\equiv \frac{K}{f_0}{\int }_0^{z}N(z)dz$ in the integral from the surface to depth, while $a=\frac{f_0{N}_z(-D)}{2K{N}^2(-D)}$, $b=\frac{f_0{N}_z(0)}{2K{N}_0^2}$, and $\alpha =\frac{1-a}{1+a}$.

For constant stratification this reduces to the result in Appendix .

C AppendixClass Hierarchy

The algorithms in this manuscript are implemented as classes in Matlab in order to take advantage of class‐based inheritance. The class hierarchy is shown in Figure , where InternalModeBase is the abstract superclass which defines the primary interface. A class cluster InternalModes (not part of the hierarchy) is included to provide a single interface from which to initialize all the concrete subclasses.

Fig. C1. Class hierarchy for the Matlab implementation of the algorithms in this manuscript. InternalModesBase is an abstract class.

The primary Spectral class uses depth ( $z$) coordinates from section to compute the eigenvalue problem, while the DensitySpectral and WKBSpectral subclasses use the stretched coordinates described in sections  and , respectively. The AdaptiveSpectral class overrides its superclass when the frequency is high enough, as described in section . The ConstantStratification class implements the analytical solution from Appendices and , while the ExponentialStratification class implements the analytical solution from Appendices and . The WKB class implements the WKB approximated solution from Appendices  and and inherits from the Spectral class in order to use the spectrally computed stratification function, $N^2$. The FiniteDifference class uses finite difference matrices of arbitrary order, on arbitrary grids using the algorithm described in Fornberg ().

Acknowledgments

All classes available online (https://github.com/JeffreyEarly/GLOceanKit/). Thanks to Jonathan Lilly for his careful comments on an early draft of this manuscript. This work was supported by the Office of Navy Research Grant N00014‐15‐1‐2465, NASA grant 80NSSC18K0070, and the National Science Foundation Awards 1658564 and 1536747.