# Propagation of internal gravity waves: optimal vertical discretisation and grid design based on numerical error analysis

Gabriel Derrida, Laurent Debreu, Florian Lemarié, Jérome Chanut

► To cite this version:

Gabriel Derrida, Laurent Debreu, Florian Lemarié, Jérome Chanut. Propagation of internal gravity waves: optimal vertical discretisation and grid design based on numerical error analysis. 2026. hal-05514655

HAL Id: hal-05514655

https://inria.hal.science/hal-05514655v1

Preprint submitted on 17 Feb 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
RESEARCH ARTICLE

# Propagation of internal gravity waves: optimal vertical discretisation and grid design based on numerical error analysis

Gabriel Derrida<fnref n="1" /> | Laurent Debreu<fnref n="1" /> | Florian Lemarié<fnref n="1" />

| Jérômenut<fnref n="2" />

Correspondence

Gabriel Derrida, Univ. Grenoble Alpes, Inria,
CNRS, Grenoble INP, LJK, 38000, Grenoble,
France

Email: gabriel.derrida@inria.fr

Funding information

French National Research Agency (ANR),
Grant/Award Number:
ANR-23-CE56-0006-01,
ANR-22-POCE-0003

Focusing on the internal gravity wave dynamic, we introduce a general framework for the vertical discretisation of stratified linearized primitive equations. We derive a discrete formulation that preserves energy conservation and leads to a consistent discrete Sturm-Liouville problem. A $\beta$-indexed family of discretisations is introduced, and a perturbation analysis provides an explicit expression for the leading-order eigenvalue errors. We show that an appropriate choice of parameter $\beta$ significantly reduces the leading-order eigenvalue error and, for uniform grids, makes it vanish, yielding fourth-order instead of second-order convergence. An explicit optimal vertical grid is also derived for a single eigenmode, and the extension to grids optimized for multiple modes is discussed.

KEYWORDS

internal gravity waves, vertical discretisation, grid design,
eigenvalue error, optimal grid, Charney-Phillips, Lorenz

<footnote n="1">¹Univ. Grenoble Alpes, Inria, CNRS,
Grenoble INP, LJK, Grenoble, France</footnote>

<footnote n="2">²Mercator Océan International, Toulouse,
France</footnote>
# 1 | INTRODUCTION

Internal gravity waves (IGWs) are a key component of ocean dynamics, transporting energy and momentum across a broad range of spatial and temporal scales. They contribute to vertical mixing, influence large-scale stratification, and play an essential role in redistributing energy from tides (Smith and Young, 2002; Nycander, 2005) and mesoscale motions (Arbic et al., 2012; Barkan et al., 2021). Because of their impact on water-mass transformations and climate-relevant processes (Saint-Laurent et al., 2012; Nikurashin and Ferrari, 2013), accurately representing IGWs has become a central challenge for modern ocean modelling systems. Moreover, the joint increase in horizontal and vertical resolution in ocean models is now enabling the resolution of an ever larger fraction of the key dynamical processes associated with internal-wave dynamics. The accuracy of the numerical simulation of IGWs is particularly sensitive to several aspects of ocean model formulation: (i) the vertical grid structure including grid spacing, staggering, and coordinate systems (e.g. Thuburn and Woollings, 2005; Holdaway et al., 2013; Stewart et al., 2017), (ii) the discretisation of the horizontal pressure gradient and continuity equation (e.g. Thuburn, 2006; Ullrich, 2014; Demange, 2014); and (iii) the horizontal resolution (e.g. Hallberg, 2013; Buijsman et al., 2020).

Theoretical studies of IGWs in a continuously stratified fluid typically reduce to a Sturm-Liouville (SL) problem (e.g. Gill, 1982; Kelly, 2016; Early et al., 2020). In a bounded domain such as the ocean, this problem yields a discrete set of eigenfunctions, the baroclinic modes, which form a complete basis for the vertical structure of internal waves. Because each mode has well-defined dispersion properties, analysing internal-wave representation through the baroclinic-mode framework provides a natural and rigorous way to assess the accuracy of their numerical simulation. In this study, we develop a quantitative, model-independent framework that links vertical grid design and discretisation choices to the expected numerical errors in the dispersion properties of baroclinic modes.

The construction of vertical coordinate systems has long been recognized as a critical element of ocean model design (Griffies et al., 2000). Over the years, numerous studies have examined the influence of different vertical staggering arrangements (Tokioka, 1978; Arakawa, 1988; Leslie and Purser, 1992; Thuburn, 2006) and have proposed various strategies for optimising the placement of vertical grid levels Baer and Ji (1989);?. Despite improvements in model resolution and numerical schemes, the design of an "optimal" vertical grid tailored to IGW dynamics remains non-trivial. Recently, Robey and Grooms (2024) investigated the behavior of baroclinic modes in a continuously stratified fluid and their representation in discretized-layer models, showing how the number and placement of vertical layers affect the fidelity of baroclinic-mode representation.

Building on these previous studies, we examine how vertical grid design affects the accurate representation of internal wave modes. We aim to identify configurations, both in terms of staggering, discretisation and vertical-level placement, that minimize truncation errors in baroclinic-mode discretisation. Using classical perturbation theory for eigenvalue problems, we derive a leading-order estimate of the discretisation error, providing a quantitative criterion to compare grids and identify designs that are optimal in minimizing this error. Even though the application context of this work is ocean-oriented, it is equally relevant to the atmospheric context.

## Outline of the paper

The paper is organized as follows. In Section 2, we present the continuous theory of normal-mode decomposition, derive the classical Sturm-Liouville problem, and provide the expression for the conserved energy of the continuous system. In Section 3, we formulate a general vertical discretisation framework equipped with energy-conserving constraints at the discrete level, and subsequently formulate several classical discretisation schemes as particular instances of this framework (namely the Lorenz and Charney-Phillips grids with second-order discretizations). In Section 4, we propose a generalisation of these standard discretisations that, for an appropriate choice of parameters,
achieves fourth-order accuracy while retaining the same stencil as the conventional second-order schemes, in the constant-stratification case. In Sections 5 and 6, we extend this numerical error analysis to the cases of an arbitrary stratification and a generalised vertical coordinate, using perturbation theory for eigenvalue problems. To complement this analysis of errors associated with a given discretisation scheme, Section 7 examines the influence of the vertical grid on the numerical error. For constant stratification, we can analytically derive the optimal grid for a given discretisation scheme that minimises the second-order error in the eigenvalues. For slowly varying stratification, using Wentzel-Kramers-Brillouin (WKB) theory, we show that the Charney coordinate proposed by Robey and Grooms (2024) is optimal in the sense that it minimises the sum of the second-order error terms when a large number of baroclinic modes are considered. This result, however, applies only to a specific discretisation that lies between the Lorenz and Charney-Phillips grids. Finally, we propose a numerical grid-optimization protocol that can be applied regardless of the discretisation scheme or the stratification. Section 8 summarises the main findings and conclusions.

## 2 | CONTINUOUS THEORY OF NORMAL-MODE DECOMPOSITION

We recall the classical normal-mode decomposition for internal gravity waves within the 2D linearised framework used here. This formulation naturally reduces to a Sturm-Liouville (SL) problem, whose eigenvalues and eigenfunctions provide, via mode decomposition, a representation of the solution of the original 2D system.

### 2.1 | 2D linearised system

We consider the 2D $(x-z)$ linearised, inviscid primitive equations over a flat bottom, under the hydrostatic and Boussinesq assumptions. For more details about the vertical-mode decomposition see Gill (1982) p.167 and Blayo and Debreu (2005). $u(x, z, t)$ and $w(x, z, t)$ denote the perturbation components of the horizontal and vertical velocities, while $p(x, z, t)$ and $b(x, z, t)$ denote dynamical pressure and buoyancy perturbations around a background state $(\bar{\rho}(z), \bar{b}(z))$ satisfying the hydrostatic balance $(\frac{d\bar{\rho}(z)}{dz} = \bar{b}(z))$. Using these notations, the system of interest is

$$ \partial_x u + \partial_z w = 0 \qquad (1a) $$

$$ \partial_t u + \partial_x p = 0 \qquad (1b) $$

$$ \partial_z p - b = 0 \qquad (1c) $$

$$ \partial_t b + w N^2(z) = 0 \qquad (1d) $$

where $N^2(z)$ is the Brunt-Väisälä frequency defined by $N^2(z) = \frac{d\bar{b}}{dz}$. We assume a stably stratified fluid, so $N^2(z) > 0$. Horizontally, the domain extends from $x = 0$ to $x = L$, while vertically the model extends from the flat bottom at $z = -H$ to the ocean surface at $z = 0$, where a rigid-lid condition is imposed. Under these assumptions, the surface and bottom boundary conditions read:

$$ w(z = 0) = w(z = -H) = 0 \qquad (2) $$
**Remark:**

The rigid-lid approximation removes surface gravity waves and simplifies the analysis by isolating internal modes. A free-surface boundary condition, commonly used in ocean models, could also be considered (see, for example, Kelly, 2016) but this would introduce unnecessary complications without enhancing the relevance of our study.

## 2.2 | Energy conservation

Since the discretization schemes presented later in Sec. 3 will be required to preserve energy, we recall the corresponding conservation property at the continuous level. The linear system (1) is energy-conserving. On $\Omega = [0, L] \times [-H, 0]$, with the mean density normalized to unity, the total kinetic and potential energy are defined by:

$$E_C(t) = \frac{1}{2} \int_{\Omega} u^2 d\Omega \quad \text{and} \quad E_P(t) = \frac{1}{2} \int_{\Omega} \frac{b^2}{N^2(z)} d\Omega \qquad (3)$$

Multiplying equation (1b) by $u$ and (1c) by $w$ gives

$$\partial_t E_C = - \int_{\Omega} \nabla \cdot (v \cdot p) \, d\Omega + \int_{\Omega} w \, b \, d\Omega$$

where $v = (u, w)^T$. Using equation (1d) yields

$$\partial_t E_P = - \int_{\Omega} w \, b \, d\Omega$$

Hence

$$\partial_t (E_C + E_P) = - \int_{\partial \Omega} \rho v \cdot n \, d(\partial \Omega) \qquad (4)$$

where $\partial\Omega$ denotes the boundary of the domain. Thus, the mechanical energy evolves solely through fluxes across the domain boundaries.

## 2.3 | Eigenvalue problem and normal mode decomposition

Assuming that all fields $u$, $p$, $w$ and $b$ are sufficiently smooth, the 2D linearised system (1) can be reduced to a single equation for the vertical velocity $w$ by eliminating $u$, $p$, and $b$:

$$\partial_{zztt} w + N^2(z) \partial_{xx} w = 0 \qquad (5)$$

This equation is obtained by differentiating the continuity equation with respect to $t$ and eliminating $u$ using the horizontal momentum equation (1b). The resulting equation is then differentiated with respect to $z$ to eliminate $p$ via the hydrostatic balance (1c). Finally, differentiating again with respect to $t$ removes $b$ using the buoyancy equation (1d).

The compact representation (5) highlights the coupled propagation of internal waves in the horizontal and vertical directions and serves as the starting point for a normal-mode decomposition. Equation (5) is invariant under translations in $x$ and $t$ because its coefficients do not depend on these variables. We therefore seek plane-wave solutions of the
form

$$w(x, z, t) = W(z) e^{i(kx - \omega t)} \qquad (6)$$

where the vertical structure $W(z)$ is determined by the following SL problem:

$$\frac{d^2}{dz^2} W(z) + \lambda N^2(z) W(z) = 0 \qquad (7a)$$

$$W(0) = W(-H) = 0 \qquad (7b)$$

with $\lambda = \frac{k^2}{\omega^2}$. Because the vertical domain is bounded, system (7) admits a countable set of solutions $W_n(z)$ (the vertical normal modes) with corresponding eigenvalues $\lambda_n = \frac{k_n^2}{\omega_n^2}$. Each vertical mode is associated with a specific horizontal wavenumber $k_n$ for a given frequency $\omega$. The vertical normal modes are orthogonal with respect to the weighted scalar product $\langle f, g \rangle = \frac{1}{H} \int_{-H}^{0} N^2(z) f(z) g(z) dz$.

The general solution of the original linearised system (1) can be expressed as a superposition of these vertical normal modes:

$$\begin{aligned} w(x, z, t) &= \sum_{n=1}^{\infty} w_n(x, t) W_n(z) \quad ; \quad b(x, z, t) = N^2(z) \sum_{n=1}^{\infty} b_n(x, t) W_n(z) \\ \rho(x, z, t) &= \sum_{n=1}^{\infty} \rho_n(x, t) \frac{d}{dz} W_n(z) \quad ; \quad u(x, z, t) = \sum_{n=1}^{\infty} u_n(x, t) \frac{d}{dz} W_n(z) \end{aligned}$$

Using the hydrostatic relation (1c) and the continuity equation (1a), one obtains (here $b_n$ is the projection of $b(z)N^{-2}(z)$ onto the vertical modes $W_n$):

$$\begin{aligned} \rho_n(x, t) &= -\frac{1}{\lambda_n} b_n(x, t) \\ w_n(x, t) &= -\partial_x u_n(x, t). \end{aligned}$$

The momentum (1b) and buoyancy (1d) equations then lead to the following set of 1D equations for each mode $n$:

$$\partial_t u_n - c_n^2 \partial_x b_n = 0, \qquad (8a)$$

$$\partial_t b_n - \partial_x u_n = 0, \qquad (8b)$$

where $c_n = \frac{1}{\sqrt{\lambda_n}}$ is the phase speed of mode $n$. One goal of this paper is to study the discrete counterparts of the SL problem (7).

## 3 | VERTICAL DISCRETISATION AND GENERALISED EIGENVALUE PROBLEM

In this section, we introduce a general framework for analysing vertical discretisations of system (1), focusing on energy conservation and the discrete counterpart of the Sturm-Liouville problem (7). This framework is then applied to classical examples, including second-order schemes on Lorenz (1960) and Charney and Phillips (1953) grids.
## 3.1 | 2D discrete linearised system

We consider a vertical discretisation of system (1) using a generic discretisation scheme. Each continuous variable is represented by a finite set of values along the vertical direction. As a result, the four variables in (1) are represented by four vectors of discrete values:

$$ 
\begin{aligned}
\mathbf{W}(x, t) &= \{W_j(x, t)\}_{N_W} ; \quad \mathbf{U}(x, t) = \{U_j(x, t)\}_{N_U} \\
\mathbf{P}(x, t) &= \{P_j(x, t)\}_{N_P} ; \quad \mathbf{B}(x, t) = \{B_j(x, t)\}_{N_B}
\end{aligned}
$$

where, for a given vector $\mathbf{X}(x, t) = \{X_j(x, t)\}_{N_X}$ ($\mathbf{X} = \mathbf{W}, \mathbf{U}, \mathbf{P}, \mathbf{B}$), $N_X$ denotes the dimension of the vector $\mathbf{X}$. Each variable can be discretized on a different grid, resulting in vectors of potentially different dimensions $N_W, N_U, N_P$, and $N_B$. The discrete equivalent of system (1) reads

$$ \mathcal{A}_U \partial_x \mathbf{U}(x, t) + \mathcal{D}_W \mathbf{W}(x, t) = 0 \qquad (9a) $$

$$ \mathcal{A}_U \partial_t \mathbf{U}(x, t) + \mathcal{A}_P \partial_x \mathbf{P}(x, t) = 0 \qquad (9b) $$

$$ \mathcal{D}_P \mathbf{P}(x, t) - \mathcal{A}_B \mathbf{B}(x, t) = 0 \qquad (9c) $$

$$ \mathcal{A}_B \partial_t \mathbf{B}(x, t) + \mathcal{A}_W \mathbf{W}(x, t) = 0 \qquad (9d) $$

where the matrices $\mathcal{A}_U, \mathcal{D}_W, \mathcal{A}_P, \mathcal{D}_P, \mathcal{A}_B, \mathcal{A}_W$ are determined by the vertical grid and the discretisation method. $\mathcal{A}_U$ and $\mathcal{A}_B$ are symmetric positive-definite matrices that define the discrete inner products for horizontal velocity and buoyancy, respectively. For consistency, the dimensions of these matrices must satisfy

* $\mathcal{A}_U$ is a square matrix of size $N_U \times N_U$,

* $\mathcal{A}_B$ is a square matrix of size $N_B \times N_B$,

* $\mathcal{D}_W$ is a matrix of size $N_U \times N_W$,

* $\mathcal{D}_P$ is a matrix of size $N_B \times N_P$,

* $\mathcal{A}_P$ is a matrix of size $N_U \times N_P$,

* $\mathcal{A}_W$ is a matrix of size $N_B \times N_W$.

The Lorenz and Charney-Phillips grids and associated second-order schemes discussed in Section 3.4 are specific instances of this general framework.

**Remark** The Lorenz and Charney-Phillips discretisations satisfy the following property: there exists two matrices $M_{W \to B}$ and $M_{P \to U}$ such that

$$ \mathcal{A}_W = \mathcal{A}_B M_{W \to B}, \quad \mathcal{A}_P = \mathcal{A}_U M_{P \to U} \qquad (10) $$

The role of $M_{W \to B}$ and $M_{P \to U}$ is to map the grids where B and U are defined to the grids where W and P are defined.
When (10) holds, the discrete system (9) simplifies to

$$ \mathcal{A}_U \partial_x \mathbf{U}(x, t) + \mathcal{D}_W \mathbf{W}(x, t) = 0 \qquad (11a) $$

$$ \partial_t \mathbf{U}(x, t) + \mathcal{M}_{P \to U} \partial_x \mathbf{P}(x, t) = 0 \qquad (11b) $$

$$ \mathcal{D}_P \mathbf{P}(x, t) - \mathcal{A}_B \mathbf{B}(x, t) = 0 \qquad (11c) $$

$$ \partial_t \mathbf{B}(x, t) + \mathcal{M}_{W \to B} \mathbf{W}(x, t) = 0 \qquad (11d) $$

## 3.2 | Discrete conservation of energy

We now identify the constraints on the generic discretisation (9) or (11) that ensure energy conservation at the discrete level. The discrete analogues $\tilde{\mathcal{E}}_C$ and $\tilde{\mathcal{E}}_P$ of the continuous kinetic and potential energies (3) are defined as

$$ \tilde{\mathcal{E}}_C(t) = \int_0^L \frac{1}{2} \mathbf{U}(x,t)^T \mathcal{A}_U \mathbf{U}(x,t) = \int_0^L \frac{1}{2} \langle \mathbf{U}, \mathbf{U} \rangle_{\mathcal{A}_U} \quad \text{and} \quad \tilde{\mathcal{E}}_P(t) = \int_0^L \frac{1}{2} \mathbf{B}(x,t)^T \mathcal{A}_B \mathbf{B}(x,t) = \int_0^L \frac{1}{2} \langle \mathbf{B}, \mathbf{B} \rangle_{\mathcal{A}_B}. \qquad (12) $$

One can show (see details in Appendix A) that the time evolution of the discrete energy can be written as

$$ \partial_t (\tilde{\mathcal{E}}_C + \tilde{\mathcal{E}}_P) = - \int_0^L \mathbf{W}^T (\mathcal{D}_W^T \mathcal{A}_U^{-1} \mathcal{A}_P + \mathcal{A}_W^T \mathcal{A}_B^{-1} \mathcal{D}_P) \mathbf{P} $$

In all configurations used in this work, the vertical velocity $\mathbf{W}$ is defined at cells interfaces and vanishes at boundaries (rigid-lid assumption): $W_0 = W_{N_W-1} = 0$. Therefore the discrete energy is conserved (i.e. $\partial_t (\tilde{\mathcal{E}}_C + \tilde{\mathcal{E}}_P) = 0$ in the absence of external forcing) if the matrix $\mathcal{D}_W^T \mathcal{A}_U^{-1} \mathcal{A}_P + \mathcal{A}_W^T \mathcal{A}_B^{-1} \mathcal{D}_P$ has nonzero elements only in its first and last rows. This can be expressed using a projection matrix $\mathcal{P}_0$ which removes the first and last rows of a vector of size $N_W$, $\mathcal{P}_0 = I_{N_W} - (e_0 e_0^T + e_{N_W-1} e_{N_W-1}^T)$. Then the discrete energy is conserved if

$$ \mathcal{P}_0 (\mathcal{D}_W^T \mathcal{A}_U^{-1} \mathcal{A}_P + \mathcal{A}_W^T \mathcal{A}_B^{-1} \mathcal{D}_P) = 0. \qquad (13) $$

When (10) holds, this condition remarkably simplifies to

$$ \mathcal{P}_0 (\mathcal{D}_W^T \mathcal{M}_{P \to U} + \mathcal{M}_{W \to B}^T \mathcal{D}_P) = 0 \qquad (14) $$

## 3.3 | Discrete eigenvalue problem

As in Sec. 2.3 for the continuous case, we can derive the discrete counterpart of the eigenvalue problem (7). As a first step, we apply the same methodology to obtain the discrete analogue of equation (5). First, the continuity equation (9a) is differentiated with respect to time to obtain

$$ \mathcal{A}_U \partial_t \mathbf{U} = - \mathcal{D}_W \partial_t \mathbf{W} $$

where $\mathbf{U}$ can be eliminated using the horizontal momentum equation (9b), to obtain

$$ \mathcal{D}_W \partial_t \mathbf{W} = \mathcal{A}_P \partial_x^2 \mathbf{P}. $$
Multiplying by $P_0 D_W^T A_U^{-1}$ and using the constraint to conserve the discrete total energy (13), we have

$$P_0 D_W^T A_U^{-1} D_W \partial_t \mathbf{W} = -P_0 A_W^T A_B^{-1} D_P \partial_{x_2} \mathbf{P}$$

From the hydrostatic equation (9c), we can eliminate $\mathbf{P}$ and obtain

$$P_0 D_W^T A_U^{-1} D_W \partial_t \mathbf{W} = -P_0 A_W^T A_B^{-1} A_B \partial_{x_2} \mathbf{B}.$$

Finally, taking the time derivative and employing the buoyancy equation (9d) eliminates $\mathbf{B}$, leading to the discrete analogue of equation (5)

$$-P_0 D_W^T A_U^{-1} D_W \partial_{t_2} \mathbf{W} + P_0 A_W^T A_B^{-1} A_W \partial_{x_2} \mathbf{W} = 0. \quad (15)$$

As in the continuous case (6), looking for solutions in the form $\tilde{W}_j(x, t) = W_j e^{i(kx-\omega t)}$, the discrete version of equation (7) is

$$\Delta \mathbf{W} + \lambda^D \mathbf{I} \mathbf{W} = 0 \quad (16)$$

where $\Delta = -P_0 D_W^T A_U^{-1} D_W$ denotes the matrix of the discrete Laplacian operator<fnref n="1" /> corresponding to the second-order vertical derivative in (7). The operator $\mathbf{I} = P_0 A_W^T A_B^{-1} A_W$ represents the vertical interpolation operator for modes on the computational grid. The quantity $\lambda^D = \frac{k^2}{\omega^2}$ is an eigenvalue of the discrete problem. When (10) holds, matrix $\mathbf{I}$ writes

$$\mathbf{I} = P_0 M_{W \to B}^T A_B M_{W \to B} \quad (17)$$

Under the rigid-lid assumption, system (16) is complemented by homogeneous Dirichlet boundary conditions for the vertical velocity $\mathbf{W}$: $W_0 = W_{N_W-1} = 0$.

**Remark** Including a constant Coriolis parameter $f$ does not change the vertical eigenvalue problem derived above. One obtains, after elimination of the horizontal velocity components,

$$-P_0 D_W^T A_U^{-1} D_W (\partial_t^2 \mathbf{W} + f^2 \mathbf{W}) + P_0 A_W^T A_B^{-1} A_W (\partial_x^2 + \partial_y^2) \mathbf{W} = 0.$$

Looking for separable solutions $\mathbf{W} \propto e^{i(kx+l y-\omega t)}$ again leads to the same discrete Sturm-Liouville form

$$\Delta \mathbf{W} + \lambda^D \mathbf{I} \mathbf{W} = 0,$$

but with the modified eigenparameter

$$\lambda^D = \frac{k^2 + l^2}{\omega^2 - f^2},$$

<footnote n="1">¹Strictly speaking, the proper discrete Laplacian operator should be written as $Q_W^{-1} \Delta$ where the matrix $Q_W$ is defined as the discrete inner product $W^T Q_W W = \int W(z)^2 dz$. A similar correction applies for the operator $\mathbf{I}$.</footnote>
instead of $\lambda^D = k^2/\omega^2$. Thus, as is known, the Coriolis force only shifts the horizontal dispersion relation (via $\omega^2 \mapsto \omega^2 - f^2$) and does not affect the vertical operators $\Delta$ and $I$ nor the vertical structure of the eigenmodes. This invariance is relevant for the error-compensation arguments developed below in Section 4 and thereafter.

## 3.4 | Classical examples of vertical grid and associated discretization

This section demonstrates that common discretization schemes used in oceanic and atmospheric modeling fit within the general framework introduced in the previous subsection.

Let us consider a vertical grid with $M + 1$ points ($M$ layers) bounded by $z_0 = -H$ and $z_M = 0$, such as

$$-H = z_0 < z_1 < z_2 < \dots < z_{M-1} < z_M = 0. \qquad (18)$$

The thickness of the $j^{th}$ cell is denoted by $\Delta z_j = z_j - z_{j-1}$, for $j = 1, \dots, M$.

A common approach in oceanic and atmospheric models is to use staggered grids, such as the Lorenz (1960) and the Charney and Phillips (1953) vertical grids. The primary distinction between these two grids lies in the placement of density (or equivalently buoyancy): in the Lorenz grid, it is stored at the cell centers, whereas in the Charney-Phillips grid, it is stored at cell interfaces. Additionally, the modified Lorenz grid Arakawa (1988), shown on the left side of Figure 2, stores pressure at the cell interfaces. For both Lorenz-type and Charney-Phillips grids, the intermediate “half-levels” are defined as midpoints between adjacent full levels, i.e.

$$z_{j-\frac{1}{2}} = \frac{z_j + z_{j-1}}{2}$$

so that tracer points are always centered between two consecutive interfaces. Extensive literature has compared

<figure><img src="image_2.png" /><figcaption>FIGURE 1 Vertical grid structure with interfaces $z_j$ and cell centers $z_{j-1/2}$.</figcaption></figure>

these grids; see Tokioka (1978); Arakawa and Konor (1996); Thuburn (2006) and references therein. Although both the original and modified Lorenz grids can exhibit a computational mode Arakawa (1988); Bell and White (2017), they are often preferred because they naturally accommodate a conservative discretisation around tracer cells. In the following subsections (3.4.1 and 3.4.2), we describe the corresponding discretisations and their associated discrete SL systems, which naturally arise from second-order approximations of the continuous SL problem (7). Section 4 will demonstrate that a straightforward modification of these schemes leads to a generalized formulation that maintains consistency
<figure><img src="image_3.png" /><figcaption>FIGURE 2 Illustrations of modified Lorenz grid (L grid) and Charney-Phillips grid (CP grid)</figcaption></figure>

with the continuous SL problem while achieving potentially higher-order accuracy, all within the framework of simple second-order discretisations.

### 3.4.1 | Lorenz grid with second-order approximations

We begin by considering a modified Lorenz grid, where buoyancy $B$ and horizontal velocity $U$ are staggered (at mid-levels), while pressure $P$ and vertical velocity $W$ are defined at full levels (i.e. $N_U = N_B = M$ and $N_W = N_P = M + 1$). The size and structure of the matrices arising from a second-order discretization on this grid arrangement are (only non-zero elements are given)

$$ 
\begin{align*}
\mathcal{A}_U & \in \mathbb{R}^{M \times M} & \text{with } (\mathcal{A}_U)_{j,j} &= \Delta z_j, \\
\mathcal{D}_W & \in \mathbb{R}^{M \times (M+1)} & \text{with } (\mathcal{D}_W)_{j,j} &= -1, (\mathcal{D}_W)_{j,j+1} = 1, \\
\mathcal{D}_P & \in \mathbb{R}^{M \times (M+1)} & \text{with } (\mathcal{D}_P)_{j,j} &= -1, (\mathcal{D}_P)_{j,j+1} = 1, \\
\mathcal{A}_P & \in \mathbb{R}^{M \times (M+1)} & \text{with } (\mathcal{A}_P)_{j,j} &= \Delta z_j/2, (\mathcal{A}_P)_{j,j+1} = \Delta z_j/2, \\
\mathcal{A}_B & \in \mathbb{R}^{M \times M} & \text{with } (\mathcal{A}_B)_{j,j} &= N_{j-\frac{1}{2}}^2 \Delta z_j, \\
\mathcal{A}_W & \in \mathbb{R}^{M \times (M+1)} & \text{with } (\mathcal{A}_W)_{j,j} &= N_{j-\frac{1}{2}}^2 \Delta z_j/2, (\mathcal{A}_W)_{j,j+1} = N_{j-\frac{1}{2}}^2 \Delta z_j/2.
\end{align*}
$$

where $N_{j-\frac{1}{2}}$ is the discrete value of the Brunt-Väisälä frequency at cell centers $z = z_{j-1/2}$.

This is a simple example of matrices of the system (9) where $N_U = N_B = M$ and $N_W = N_P = M + 1$.

We also note that matrices satisfy property (10) $(\mathcal{A}_W = \mathcal{A}_B M_{W \to B}, \mathcal{A}_P = \mathcal{A}_U M_{P \to U})$ where $M_{W \to B} = M_{P \to U} = M$ performs averaging and are given by

$$ M = \begin{pmatrix} \frac{1}{2} & \frac{1}{2} & 0 & 0 & \dots & 0 & 0 \\ 0 & \frac{1}{2} & \frac{1}{2} & 0 & \dots & 0 & 0 \\ 0 & 0 & \frac{1}{2} & \frac{1}{2} & \dots & 0 & 0 \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & \dots & 0 & 0 & \frac{1}{2} & \frac{1}{2} & 0 \\ 0 & \dots & 0 & 0 & 0 & \frac{1}{2} & \frac{1}{2} \end{pmatrix}_{M \times (M+1)} $$
One can notice that $D_P = D_W$, so the matrix $D_W^T M + M^T D_P$ has nonzero elements only at the boundaries (i.e. in the first and last rows). Then,

$$P_0(D_W^T M + M^T D_P) = 0$$

so the discrete energy is conserved. The discrete counterpart of system (1), obtained by substituting these matrices into (9) or (11), is

$$\Delta z_j \partial_x U_{j-\frac{1}{2}} + W_j - W_{j-1} = 0 \qquad (19a)$$

$$\partial_t U_{j-\frac{1}{2}} + \partial_x \frac{P_j + P_{j-1}}{2} = 0 \qquad (19b)$$

$$P_j - P_{j-1} - N_{j-\frac{1}{2}}^2 \Delta z_j B_{j-\frac{1}{2}} = 0 \qquad (19c)$$

$$\partial_t B_{j-\frac{1}{2}} + \frac{W_j + W_{j-1}}{2} = 0 \qquad (19d)$$

Using the derivation presented in Sec. 3.3, the discrete version of the SL problem in the form of (16) obtained from (19) is

$$\underbrace{\frac{W_{j+1} - W_j}{\Delta z_{j+1}} - \frac{W_j - W_{j-1}}{\Delta z_j}}_{\Delta W} + \lambda^D \underbrace{\frac{1}{4} \left( N_{j+\frac{1}{2}}^2 \Delta z_{j+1} (W_{j+1} + W_j) + N_{j-\frac{1}{2}}^2 \Delta z_j (W_j + W_{j-1}) \right)}_{I W} = 0, \quad j = 1 \dots M-1 \qquad (20)$$

with boundary conditions $W_0 = W_M = 0$.

The discretisation (20) yields a matrix $\Delta$ that is positive definite, while the matrix $I$ is only positive semidefinite. As a result, the spectral structure of the continuous problem is not fully preserved, and a spurious numerical mode inevitably appears, as was expected for the Lorenz grid.

### 3.4.2 | Charney-Phillips grid with second-order approximations

We now consider a Charney-Phillips grid in which pressure $P$ and horizontal velocity $U$ are staggered (at mid-levels), while buoyancy $B$ and vertical velocity $W$ are defined at full levels (i.e. $N_U = N_P = M$ and $N_W = N_B = M + 1$). Again, assuming a second-order discretization, the size and the structure of the corresponding matrices are:

$$
\begin{align*}
\mathcal{A}_U & \in \mathbb{R}^{M \times M} & \text{with } (\mathcal{A}_U)_{j,j} &= \Delta z_j, \\
\mathcal{D}_W & \in \mathbb{R}^{M \times (M+1)} & \text{with } (\mathcal{D}_W)_{j,j} &= -1, \quad (\mathcal{D}_W)_{j,j+1} = 1, \\
\mathcal{D}_P & \in \mathbb{R}^{(M+1) \times (M)} & \text{with } (\mathcal{D}_P)_{j,j} &= 1, \quad (\mathcal{D}_P)_{j+1,j} = -1, \\
\mathcal{A}_P & \in \mathbb{R}^{M \times M} & \text{with } (\mathcal{A}_P)_{j,j} &= \Delta z_j, \\
\mathcal{A}_B & \in \mathbb{R}^{(M+1) \times (M+1)} & \text{with } (\mathcal{A}_B)_{j,j} &= N_j^2 \frac{\Delta z_j + \Delta z_{j+1}}{2} \quad \text{for } 1 \le j \le M-1, \quad (\mathcal{A}_B)_{0,0} = N_0^2 \frac{\Delta z_1}{2}, \quad (\mathcal{A}_B)_{M,M} = N_M^2 \frac{\Delta z_M}{2}, \\
\mathcal{A}_W & \in \mathbb{R}^{(M+1) \times (M+1)} & \text{with } (\mathcal{A}_W)_{j,j} &= N_j^2 \frac{\Delta z_j + \Delta z_{j+1}}{2} \quad \text{for } 1 \le j \le M-1, \quad (\mathcal{A}_W)_{0,0} = N_0^2 \frac{\Delta z_1}{2}, \quad (\mathcal{A}_W)_{M,M} = N_M^2 \frac{\Delta z_M}{2}
\end{align*}
$$

where $N_j$ is the discrete value of the Brunt-Väisälä frequency at cell interfaces $z = z_j$.

It can be verified that the energy conservation condition (13) is satisfied. Indeed, since the matrices $\mathcal{A}_W$ and $\mathcal{A}_B$, as well as $\mathcal{A}_U$ and $\mathcal{A}_P$, are identical, the property (10) holds with $M_W \to B$ and $M_P \to U$ two identity matrices (of size $M+1$)
and M) so that the constraint (14) reduces to

$$P_0(\mathcal{D}_W^T + \mathcal{D}_P) = 0, \qquad (21)$$

which is fulfilled by the above definitions of $\mathcal{D}_W$ and $\mathcal{D}_P$. The discrete counterpart of system (1) on a Charney-Phillips grid is

$$\Delta z_j \partial_x U_{j-\frac{1}{2}} + W_j - W_{j-1} = 0 \qquad (22a)$$

$$\partial_t U_{j-\frac{1}{2}} + \partial_x P_{j-\frac{1}{2}} = 0 \qquad (22b)$$

$$P_{j+\frac{1}{2}} - P_{j-\frac{1}{2}} - N_j^2 \frac{\Delta z_j + \Delta z_{j+1}}{2} B_j = 0 \qquad (22c)$$

$$\partial_t B_j + W_j = 0 \qquad (22d)$$

The discrete version of the SL problem in the form of (16) obtained from (22) is

$$\frac{W_{j+1} - W_j}{\Delta z_{j+1}} - \frac{W_j - W_{j-1}}{\Delta z_j} + \frac{\lambda^D}{2} N_j^2 W_j (\Delta z_j + \Delta z_{j+1}) = 0, \quad j = 1 \dots M-1 \qquad (23)$$

with boundary conditions $W_0 = W_M = 0$.

This discretisation (23), based on the Charney-Phillips grid, has the advantage that no interpolation is required to evaluate the discrete modes, in contrast with the Lorenz grid. Moreover, the associated matrices $\Delta$ and $I$ are both positive definite, which ensures that the discrete eigenvalues $\lambda^D$ remain real and strictly positive (no spurious modes appear).

As for the Lorenz grid case, (23) is again a second-order approximation of the continuous system (7) at $z = z_j$.

## 4 | GENERALIZATION OF SECOND-ORDER APPROXIMATIONS

In the previous section, we derived the vertical discretisation of the continuous problem and introduced usual discretisations corresponding to Lorenz and Charney-Phillips grids. These discretisations lead to second-order accuracy (see Section 5 for details) and conserve energy. This section proposes a generalisation of these discretisations aimed at improving the order of accuracy while preserving energy conservation at the discrete level, without altering the fundamental structure of the second-order centred schemes.

A detailed modified equation analysis (see Appendix B.1), in the general case of varying layer thicknesses and Brunt-Väisälä frequency, reveals that the discretizations (20) and (23) both approximate the continuous SL problem (7) with second-order accuracy. To illustrate this point and how it can be improved, consider the simplified case of a uniform grid with constant layer thickness $\Delta z$ and constant Brunt-Väisälä frequency $N(z) = N_0$. In this case, the discretizations (20) and (23) reduce to

$$\frac{W_{j+1} - 2W_j + W_{j-1}}{\Delta z^2} + \lambda^D N_0^2 \left( W_j + \frac{\beta}{4} (W_{j+1} - 2W_j + W_{j-1}) \right) = 0 \qquad (24)$$

with $\beta = 1$ for the Lorenz grid case (20), and $\beta = 0$ for the Charney-Phillips grid case (23). System (24) approximates
the second-order derivative (first term of (24)) with a truncation error given by

$$ \frac{W_{j+1} - 2W_j + W_{j-1}}{\Delta z^2} = \left. \frac{d^2 W}{dz^2} \right|_{z=z_j} + \frac{\Delta z^2}{12} \left. \frac{d^4 W}{dz^4} \right|_{z=z_j} + O(\Delta z^4). \quad (25) $$

Using the fact that, for the underlying eigenvalue problem, we have at second order

$$ \left. \frac{d^4 W}{dz^4} \right|_{z=z_j} = -\lambda^D N_0^2 \left. \frac{d^2 W}{dz^2} \right|_{z=z_j} + O(\Delta z^2), \quad (26) $$

this can be rewritten as

$$ \frac{W_{j+1} - 2W_j + W_{j-1}}{\Delta z^2} = \left(1 - \frac{\Delta z^2}{12} \lambda^D N_0^2\right) \left. \frac{d^2 W}{dz^2} \right|_{z=z_j} + O(\Delta z^4). \quad (27) $$

The modified equation associated to (24) can thus be expressed as

$$ \left. \frac{d^2 W}{dz^2} \right|_{z=z_j} + \frac{\lambda^D N_0^2}{1 + \frac{\Delta z^2}{4} \lambda^D N_0^2 (\beta - 1/3)} W(z_j) + O(\Delta z^4) = 0. \quad (28) $$

Comparing with the exact continuous equation $\frac{d^2 W}{dz^2} + \lambda N_0^2 W = 0$, we obtain the relation

$$ \lambda N_0^2 = \frac{\lambda^D N_0^2}{1 + \frac{\Delta z^2}{4} \lambda^D N_0^2 (\beta - 1/3)} + O(\Delta z^4), $$

which yields

$$ \lambda^D = \frac{\lambda}{1 - \frac{\Delta z^2}{4} \lambda N_0^2 (\beta - 1/3)} + O(\Delta z^4) = \lambda - \frac{1 - 3\beta}{12} N_0^2 \lambda^2 \Delta z^2 + O(\Delta z^4). \quad (29) $$

For the Lorenz grid ($\beta = 1$) and for the Charney-Phillips grid ($\beta = 0$), we thus have respectively

$$ \lambda^{D,L} = \frac{\lambda}{1 - \frac{\Delta z^2}{6} \lambda N_0^2} + O(\Delta z^4), \quad \lambda^{D,CP} = \frac{\lambda}{1 + \frac{\Delta z^2}{12} \lambda N_0^2} + O(\Delta z^4). \quad (30) $$

In this simple case of a uniform grid and constant $N$, the discrete eigenvalues can actually be computed analytically to confirm the previous developments. By inserting $W_j = A \sin(j\theta_n)$ with $\theta_n = \frac{n\pi}{M}$ for the bottom and top boundary conditions to be satisfied, they are given by

$$ \lambda_n^D = \frac{4 \sin^2(\frac{n\pi}{2H} \Delta z)}{\Delta z^2 N_0^2 \left(1 - \beta \sin^2(\frac{n\pi}{2H} \Delta z)\right)}, \quad n = 1, \dots, M-1 \quad (31) $$

while the exact continuous eigenvalues are

$$ \lambda_n = \frac{n^2 \pi^2}{N_0^2 H^2}, \quad n = 1, \dots \quad (32) $$
By expanding (31) in Taylor series with respect to $\Delta z$ and using (32), one can verify relations (29). It is also clear from (31) that $\lambda_n^{D,L} > \lambda_n$ and $\lambda_n^{D,CP} < \lambda_n$ for all $n$.

As just shown, both discretizations are second-order accurate. Truncation errors originate from the approximations of both the second derivative and, for the Lorenz grid, the term involving $N^2W$. In order to improve the accuracy of the discretization, the obvious way is to consider higher-order approximations of the second derivative and, for the Lorenz grid, the term involving $N^2W$. However, more simply, we can consider modifying the approximation of the term involving $N^2W$ to counterbalance the leading truncation error arising from the second-derivative approximation. This idea — namely, modifying the discrete representation of the $N^2W$ term to compensate the leading truncation error from the second-derivative approximation — will be developed in detail below.

The discrete form of the second term involving $N^2W$ in (9) corresponds to the $I$ matrix defined by (17) and involves the operator $\mathcal{A}_B$, which defines the scalar product for buoyancy, and the operator $M_{W \to B}$, which performs the interpolation from vertical velocity points to buoyancy points. Modifying the interpolation operator $M_{W \to B}$ is possible but more intricate since it also impacts the energy conservation condition (14). The other option is to modify the operator $\mathcal{A}_B$ only; we consider this latter option below.

Note that the energy conservation (14) is automatically maintained if we modify $\mathcal{A}_B$ while keeping $D_P, D_W, M_{W \to B}, M_{P \to U}$ as for the usual second order discretisation scheme.

The idea is now to derive expressions for $\mathcal{A}_B$ such that we generalize the discretizations leading to (20) and (23) by making $\beta$ a tunable parameter. This generalized discretization will impact the order of accuracy of the scheme according to the value of $\beta$. As we will see later in section 7, it will also affect the design of optimal distributions of vertical levels.

Whether in the Lorenz case or in the Charney-Phillips case, hydrostatic balance does not involve any interpolation because $B$ and the derivative of $P$ with respect to $z$ are collocated. Here, we propose a more general formulation that incorporates neighboring points in the definition of buoyancy in the hydrostatic balance. This generalized discretization originates from the Charney-Phillips grid, in which pressure $P$ and horizontal velocity $U$ are staggered (defined at midpoints), while buoyancy $B$ and vertical velocity $W$ are located at full levels. In matrix notation, this amounts to considering a $\beta$-dependent matrix $\mathcal{A}_B$ in (9). For non constant layer thickness and non constant $N$, the general term of the matrix-vector product $\mathcal{A}_B B$ is defined as

$$ (\mathcal{A}_B B)_j = \frac{\beta \Delta z_{j+1}}{2} N_{j+\frac{1}{2}}^2 \tilde{B}_{j+\frac{1}{2}} + (1-\beta) \frac{\Delta z_j + \Delta z_{j+1}}{2} N_j^2 B_j + \frac{\beta \Delta z_j}{2} N_{j-\frac{1}{2}}^2 \tilde{B}_{j-\frac{1}{2}} \quad (33) $$

with $\tilde{B}_{j+\frac{1}{2}} = \frac{B_{j+1}+B_j}{2}$ and $\beta \in [0, 1]$ such that the matrix $\mathcal{A}_B \in \mathbb{R}^{M \times M}$ is given by

$$ (\mathcal{A}_B)_{jj} = (1-\beta) N_j^2 \frac{\Delta z_j + \Delta z_{j+1}}{2} + \frac{\beta \Delta z_{j+1}}{4} N_{j+\frac{1}{2}}^2 + \frac{\beta \Delta z_j}{4} N_{j-\frac{1}{2}}^2, $$

$$ (\mathcal{A}_B)_{jj-1} = \frac{\beta \Delta z_j}{4} N_{j-\frac{1}{2}}^2, $$

$$ (\mathcal{A}_B)_{jj+1} = \frac{\beta \Delta z_{j+1}}{4} N_{j+\frac{1}{2}}^2. $$

We note that this generalised choice of $\mathcal{A}_B$ is symmetric. It is also easy to check that the matrix is strictly diagonally dominant and thus positive definite for $\beta \in [0, 1]$, which is a necessary condition since $\mathcal{A}_B$ defines the scalar product.
For the constant layer thickness and constant N case, we recover (24) since (33) reduces to

$$ (\mathcal{A}_B \mathbf{B})_j = N_0^2 \Delta z \left[ B_j + \frac{\beta}{4} (B_{j+1} - 2B_j + B_{j-1}) \right] $$

and it is obvious that the choice $\beta = 1/3$ ($\beta/4 = 1/12$) renders the scheme fourth order accurate since (29) becomes

$$ \lambda^D = \lambda + O(\Delta z^4). \qquad (34) $$

In (33), when $\beta = 0$, the discretization coincides with the one arising from the classical Charney-Phillips staggering. The $\beta$-dependent discretization of system (1) reads

$$ \Delta z_j \partial_x U_{j-\frac{1}{2}} + W_j - W_{j-1} = 0 \qquad (35a) $$

$$ \partial_t U_{j-\frac{1}{2}} + \partial_x P_{j-\frac{1}{2}} = 0 \qquad (35b) $$

$$ P_{j+\frac{1}{2}} - P_{j-\frac{1}{2}} - \left( \frac{\beta \Delta z_{j+1}}{2} N_{j+\frac{1}{2}}^2 \tilde{B}_{j+\frac{1}{2}} + (1-\beta) \frac{(\Delta z_{j+1} + \Delta z_j)}{2} N_j^2 B_j + \frac{\beta \Delta z_j}{2} N_{j-\frac{1}{2}}^2 \tilde{B}_{j-\frac{1}{2}} \right) = 0 \qquad (35c) $$

$$ \partial_t B_j + W_j = 0 \qquad (35d) $$

Applying again the derivation from Sec. 3.3, we obtain the following discrete SL eigenvalue problem from (35)

$$ \frac{W_{j+1} - W_j}{\Delta z_{j+1}} - \frac{W_j - W_{j-1}}{\Delta z_j} + \frac{\lambda^D}{4} \left[ \beta W_{j+1} N_{j+\frac{1}{2}}^2 \Delta z_{j+1} + \beta W_{j-1} N_{j-\frac{1}{2}}^2 \Delta z_j + W_j \left\{ (2(1-\beta)N_j^2 + \beta N_{j+\frac{1}{2}}^2) \Delta z_{j+1} + (2(1-\beta)N_j^2 + \beta N_{j-\frac{1}{2}}^2) \Delta z_j \right\} \right] = 0 \qquad (36) $$

with boundary conditions $W_0 = W_M = 0$.

Note that in the limiting case of $\beta = 1$, this discretization of the Sturm-Liouville problem coincides with the one arising from the classical Lorenz staggering. However, on the Charney-Phillips grid, the case $\beta = 1$ is not strictly admissible since it would lead to a non-positive-definite matrix $\mathcal{A}_B$. In the following, the vertical discretization is assumed to be of Charney-Phillips type, and we focus on solving (36), including the case $\beta = 1$. Note, however, that this particular case should be interpreted as the discrete formulation arising from the Lorenz grid (i.e. (19)), rather than as (35).

## 5 | SECOND-ORDER EIGENVALUE ERROR ANALYSIS

In the previous sections, we examined several discretisations of the Sturm-Liouville problem (7). In particular, we introduced a family of discretisations indexed by a parameter $\beta$ (36). We already showed that, in the particular case of a uniform grid and constant Brunt-Väisälä frequency, the choice $\beta = 1/3$ leads to a fourth-order accurate scheme. We now analyse, in the more general case, the influence of $\beta$ on the second-order error in the eigenvalue estimates. In particular, we employ a perturbation analysis approach that interprets the modified equation as a perturbation of the original continuous equation and quantifies its effect on the eigenvalues.
## 5.1 | Second-order eigenvalue error via perturbation study

The vertical computational grid is defined by a generalised vertical coordinate transformation $z(s)$, $0 \le s \le 1$. The locations of the cells interfaces $z_j$ are defined by

$$z_j = z(s_j) \quad \text{with} \quad s_j = j \Delta s \quad \text{and} \quad \Delta s = \frac{1}{M}, \quad j = 0, \dots, M \qquad (37)$$

where $M$ is the number of layers and $s_j$ are regularly spaced points. The transformation $z(s)$ is such that $z(0) = -H$ and $z(1) = 0$ and that the scale factor $h(s) = z'(s)$ is strictly positive. By construction, $h$ satisfies the constraint $\int_0^1 h(s) ds = z(1) - z(0) = H$.

**Remark** For convenience, in the following (and in the appendix), whenever a function (like $W$ or $N^2$) depends only on the variable $z(s)$ (and not on its derivatives), we abbreviate $F(z(s))$ simply as $F$. In this case, the prime notation $F'$ denotes the derivative of $F$ with respect to $z$, evaluated at $z = z(s)$. By contrast, primes on $z$ (such as $z'$, $z''$, etc.) always denote derivatives of $z(s)$ with respect to the parameter $s$.

For the generalised Charney-Phillips discretisations (36), we find (see appendix B.1) the following modified equation

$$W''_{D} + \lambda^{D} N^{2} W_{D} + \frac{1}{24} (\Delta s)^{2} \lambda^{D} \left( z'^{2} \left[ (-2 + 3\beta) (2(N^{2})') W'_{D} + W_{D}(N^{2})'' \right] + 2(-1 + 3\beta) N^{2} W''_{D} \right) + 4(-2 + 3\beta) z''(N^{2} W_{D})' = 0 \qquad (38)$$

Using classical perturbation theory (e.g. Rellich, 1969), the order-(Δs)² differences between the eigenvalue λ^D of the discrete problem (38) and the corresponding exact eigenvalue λ of the continuous problem (7), as well as the differences between their associated eigenvectors, can be computed analytically. The eigenvalue of the discrete problem λ^D and eigenvectors W_D can be expanded as

$$\lambda^D = \lambda + (\Delta s)^2 \mu + o((\Delta s)^2) \quad ; \quad W_D(z(s)) = W(z(s)) + (\Delta s)^2 V(z(s)) + o((\Delta s)^2)$$

where it is shown in the appendix B.2 that

$$\mu = \frac{\lambda}{24H} \int_{0}^{1} h(s)^{3} \left[ (3\beta - 2)(N^{2}W^{2})'' - 2(3\beta - 1)N^{2}WW'' \right] ds \qquad (39)$$

where the eigenfunction $W(z)$ is normalized such that

$$\frac{1}{H} \int_{-H}^{0} N^{2}(z) W^{2}(z) dz = \frac{1}{H} \int_{0}^{1} h(s) N^{2}(z) W^{2}(z) ds = 1.$$

Expression (39) provides an explicit second-order correction to the eigenvalues arising from the vertical discretisation (36). In particular, $\mu$ depends on the choice of the vertical grid $h(s)$ and the parameter $\beta$.

**Remark** The perturbative developments above are asymptotic and valid only in the limit $\Delta s \to 0$. They further rely on sufficient smoothness of $N^2(z)$ and of the eigenfunctions $W$ so that the derivatives appearing in the expansion are well defined and bounded. In practice these hypotheses are best satisfied for the lowest modes, which are the most
regular; for coarse meshes or non-smooth stratification the higher-order and representation errors may dominate and the estimate may lose accuracy. This point is further illustrated in Section 5.2.

To illustrate (39), we consider a simplified framework in which the analytical solution of the eigenvalue problem (7) is known. To do so, we assume a constant Brunt-Väisälä frequency $N(z) = N_0$. Under such an assumption, the solutions to the eigenvalue problem (7) are given by

$$W_n(z) = \frac{\sqrt{2}}{N_0} \sin \frac{n\pi}{H} z \quad \text{and} \quad \lambda_n = \frac{n^2\pi^2}{N_0^2 H^2} \quad \text{with} \quad n \in \mathbb{N}^* \qquad (40)$$

Using (40) in (39), the second-order eigenvalue error $\mu_n$ can be computed for a given function $z(s)$.

$$\mu_n = \frac{n^4 \pi^4}{12 N_0^2 H^5} \int_0^1 h(s)^3 (3\beta - 1 + 3(\beta - 1) \cos \frac{2n\pi z(s)}{H}) ds \qquad (41)$$

This expression can be further simplified when we consider a regular grid (e.g., $z(s) = H(s-1)$, $h(s) = z'(s) = H$),

$$\mu_n = \frac{(3\beta - 1) n^4 \pi^4}{12 N_0^2 H^2} \qquad (42)$$

## 5.2 | Validity of the analytic expression of the second-order error

To illustrate how the analytical error estimate can be used in practice, we can compare the prediction of the perturbation analysis with a numerical counterpart defined as

$$\mu_n^e = \frac{\lambda_n^D - \lambda_n}{\Delta s^2} \qquad (43)$$

where $\lambda_n^D$ is calculated by solving (36) and $\lambda_n$ is the exact eigenvalue. In all this section, we consider a regular grid. Note that when the continuous problem does not have an analytical solution, $\lambda_n$, used in (43) will be computed numerically with a very high-resolution grid to provide a reference value. It is important to emphasize that $\mu_n^e$ corresponds to the total eigenvalue error, while the analytical expression $\mu_n$ only provides the second-order contribution (i.e. the error captured at order $\Delta s^2$ by the perturbation expansion). The comparison between $\mu_n$ and $\mu_n^e$ therefore indicates to which extent the second-order analytical estimate succeeds in representing the total numerical error due to discretisation and under which conditions.

As mentioned above, we expect the perturbation analysis to be more accurate for low-order modes, which are smoother and better resolved on the grid. Higher-order modes, by contrast, exhibit more oscillations and are more sensitive to higher-order effects that are not captured by the second-order expansion. We also examine the accuracy of the error estimate according to how well the computational grid resolves the stratification profile. Indeed, the derivation of the modified equation (38) and the subsequent perturbation analysis assume that $N^2(z)$ and its derivatives are smooth functions that are well represented on the computational grid.

For the sake of simplicity, we only consider here a uniform grid.

We will illustrate these aspects on two profiles of Brunt-Väisälä frequency: a constant profile and a localized profile with a pycnocline.

In the case of a constant Brunt-Väisälä frequency $N(z) = N_0$, the eigenvalues of the continuous and discrete problem
have been obtained (32, 31), which allows for a direct comparison without numerical approximation.

$$ \lambda_n^D = \frac{4 \sin^2(\frac{n\pi}{2H}\Delta z)}{\Delta z^2 N_0^2 \left(1 - \beta \sin^2(\frac{n\pi}{2H}\Delta z)\right)}, \quad \lambda_n = \frac{n^2\pi^2}{N_0^2 H^2} $$

For this constant stratification case, figure 3 compares the analytical second-order error $\mu_n$ given by (42) with the numerical error estimate $\mu_n^e$ defined by (43), for two values of the discretisation parameter $\beta$ ($\beta = 0$ and $\beta = 1$), showing excellent agreement for low- and medium-order modes and an increasing relative difference for higher modes where higher-order effects become significant.

<figure><img src="image_4.png" /><figcaption>FIGURE 3 Validation of the analytical second-order error estimate. Comparison between the analytical expression $\mu_n = (3\beta - 1)n^4\pi^4/(12N_0^2H^2)$ and the numerical estimate $\mu_n^e = (\lambda_n^D - \lambda_n)/\Delta s^2$ for different values of $\beta$ with $M = 40$ grid points. The bottom-right panel shows the relative difference, confirming that the analytical approximation is very accurate for low-order modes ($n < M/4$) and, as expected, slowly degrades for higher modes due to higher-order effects. The second x-axis indicates the number of grid points per wavelength for each mode: $\sigma_n/\Delta z$ with $\sigma_n = 2H/n$.</figcaption></figure>

The localized Brunt-Väisälä frequency profile is of the form

$$ N^2(z) = \frac{A}{\cosh^2\left(\frac{z-z_0}{\sigma/2}\right)} + B. \quad (44) $$
In this function, the parameter $\sigma$ controls the width of the stratification peak centered at $z_0$, while $A$ and $B$ set the approximate amplitudes (at the pycnocline) and background level (at top and bottom) of $N^2(z)$, respectively. For the numerical tests, we set $H = 4000$ m and choose $A = (10^{-3})^2 \text{ s}^{-2}$, $B = (10^{-4})^2 \text{ s}^{-2}$, $\sigma = 160$ m and $z_0 = -100$ m, which are typical values for oceanic stratification. Note that the important parameter to vary here is $\sigma/\Delta z$, which controls how well the stratification profile is resolved on the uniform grid considered here. Figure 4 shows this stratification profile for the chosen parameters.

<figure><img src="image_5.png" /><figcaption>FIGURE 4 Vertical profile of the Brunt-Väisälä frequency $N^2(z)$ used in the numerical tests.</figcaption></figure>

Figure 5 shows the relative difference $\frac{\mu_n^e - \mu_n}{\mu_n^e}$ as a function of the relative width $\sigma/\Delta z$ for the first 19 modes, for the same two values of the discretisation parameter $\beta$ ($\beta = 0$ and $\beta = 1$). When the stratification profile is not resolved ($\sigma/\Delta z = 2$), the relative difference is large, indicating that as expected, the analytical estimate $\mu_n$ does not accurately capture the total numerical error $\mu_n^e$. This is particularly true for higher-order modes, which are more sensitive to the detailed structure of $N^2(z)$. However, as soon as the pycnocline is somewhat better resolved, the agreement improves rapidly: for example, for $\sigma/\Delta z = 4$ the relative difference already becomes moderate (less than 10%). When $\sigma/\Delta z$ reaches 6, the relative difference becomes negligible (below 2%) for all of the first 7 modes considered, indicating that the analytical estimate $\mu_n$ is a very good approximation of the total numerical error $\mu_n^e$ when the stratification profile is at least partially resolved on the grid.

<figure><img src="image_6.png" /><figcaption>FIGURE 5 Effect of resolution of $N^2$ on the accuracy of the eigenvalue error approximation. The relative error is plotted as a function of the relative width $\sigma/\Delta z$ for the first 7 modes. Left: $\beta = 0$, Right: $\beta = 1$.</figcaption></figure>
# 6 | ORDER OF APPROXIMATION

In this section, we analyze the convergence order of the eigenvalue error for the Charney-Phillips family of discretizations (36) as a function of the parameter $\beta$ introduced in Section 4. We consider both the special case of a uniform vertical grid and the more general case of a non-uniform grid.

## 6.1 | Uniform vertical grid

For the uniform grid case, expression (39) can be strongly simplified: the scale factor $h(s)$ is constant and can be taken out of the integrals, while the terms involving total derivatives vanish due to the homogeneous Dirichlet boundary conditions (rigid-lid conditions). We thus obtain

$$ \mu_n = -\frac{\lambda_n}{12} (3\beta - 1) \int_0^1 N^2 W_n W_n'' ds \quad (45) $$

As already observed in (42) in the case of a constant Brunt-Väisälä frequency, we see that for the special value

$$ \beta = \frac{1}{3}, $$

the second-order correction (of order $(\Delta s)^2$) vanishes. We emphasize that relation (45) shows that this holds even in the case of a non-constant Brunt-Väisälä frequency $N^2(z)$.

Figure 6 shows the convergence of the relative error on the first eigenvalue as a function of the number of grid points $M$, for different values of $\beta$, in the case of a uniform grid. The error plotted is the true error between the discrete eigenvalue $\lambda_1^D$ computed using (36) and the reference eigenvalue $\lambda_1$ computed either analytically or numerically at very high resolution. The left panel corresponds to the constant Brunt-Väisälä frequency case, while the right panel corresponds to the variable Brunt-Väisälä frequency profile shown in figure 4.

In both cases, we observe that all second-order schemes exhibit second-order convergence, except for the special case $\beta = \frac{1}{3}$, which shows fourth-order convergence, as expected from the previous analysis. For the variable Brunt-Väisälä frequency case, we also observe that the convergence rate for $\beta = \frac{1}{3}$ becomes even higher than fourth order for sufficiently high resolution. This is probably due to the fact that higher-order error terms (of order $(\Delta s)^4$) also have reduced coefficients for this special value of $\beta$. Note that relation (45) also indicates that the error is symmetric with respect to the value $\beta = \frac{1}{3}$: for example, the cases $\beta = 0$ and $\beta = \frac{2}{3}$ yield the same second-order error but with opposite signs. This is also confirmed by the numerical results shown in figure 6.

## 6.2 | Non-uniform vertical grid

In the more general case of a non-uniform vertical grid, expression (39) shows that the second-order error $\mu_n$ generally does not vanish for any value of $\beta$. This is due to the presence of the scale factor $h(s)$ inside the integral, which prevents the cancellation of all terms for a single value of $\beta$. Figure 8 illustrates this point by showing the convergence of the relative error on the first eigenvalue as a function of the number of grid points $M$, for different values of $\beta$, in the case of a non-uniform grid defined. The non-uniform grid is a very simple one based on a cubic polynomial stretching and where we impose a) a zero second derivative at the surface b) the ratio of the grid spacing at the surface and at the bottom to be a constant $\alpha < 1$. This leads to the following expression for the grid interfaces $z(s)$, where $s \in [0, 1]$
<figure><img src="image_7.png" /><figcaption>FIGURE 6 Convergence of the relative error on the first eigenvalue as a function of the number of grid points M, for different values of β on a uniform grid. Left: constant Brunt-Väisälä frequency N(z) = N₀. Right: variable Brunt-Väisälä frequency as in figure 4.</figcaption></figure>

is the normalized vertical coordinate:

$$z(s) = \frac{H}{1 + 2\alpha} \left[ (-1 + s)^3 - \alpha(2 + (-3 + s)s^2) \right], \quad (46)$$

which satisfies z(0) = -H, z(1) = 0, z''(1) = 0 and $\frac{z'(1)}{z'(0)} = \alpha$. In the following, we set $\alpha = 1/100$. An example of such a grid with M = 60 points is shown in figure 7.

<figure><img src="image_8.png" /><figcaption>FIGURE 7 Example of a non-uniform vertical grid defined by (46) with $\alpha = 1/100$ and M = 60 grid points.</figcaption></figure>

On figure 8, the left panel corresponds to the constant Brunt-Väisälä frequency case, while the right panel corresponds to the variable Brunt-Väisälä frequency profile shown in figure 4. The computations are performed for the first mode (top) and the fourth mode (bottom).

In both cases, we observe that all values of β yield second-order convergence, as expected from the previous analysis. The special case $\beta = \frac{1}{3}$ no longer shows fourth-order convergence, confirming that on a non-uniform grid, no value of β can eliminate the second-order error term in general. However, we note that the case $\beta = \frac{1}{3}$ still yields an order
of magnitude smaller error than other values of $\beta$, indicating that this value remains optimal.

<figure><img src="image_9.png" /><figcaption>FIGURE 8 Convergence of the relative error on the first eigenvalue as a function of the number of grid points $M$, for different values of $\beta$ on a non-uniform grid defined by (37) with $\alpha = 1/100$. Top: First mode. Bottom: Fourth mode. Left: constant Brunt-Väisälä frequency $N(z) = N_0$. Right: variable Brunt-Väisälä frequency as in figure 4.</figcaption></figure>

7 | OPTIMAL VERTICAL GRID

In previous sections, we analysed the impact of the discretisation parameter $\beta$ on the numerical error of internal gravity-wave eigenvalues. We now turn to the complementary question of determining an optimal vertical grid for a given discretisation scheme. The analytical expression of the second-order numerical error $\mu_n$ (39) helps address this question. The second-order eigenvalue error $\mu_n$ (39) is proportional to an integral of the form

$$ \mu_n = \int_0^1 z'(s)^3 f_n(z(s)) ds $$

where $f_n(z)$ depends on the stratification $N^2(z)$ and the eigenfunction $W_n(z)$. The uniqueness of an optimal grid is governed by the sign of $f_n(z)$. Indeed, if $f_n(z)$ changes sign, then $\mu_n$ can be made arbitrarily close to zero by concentrating grid points where $f_n(z)$ is negative (or positive) and coarsening the grid where it is positive (or negative). In that case, there is no unique solution for the optimal-grid problem without additional constraints on $z(s)$. The
minimisation problem is thus ill-posed, and several grid configurations can achieve an arbitrarily small value of $\mu_n$ (see Section 7.2). Obviously, if we want to optimise a grid for several modes simultaneously by minimising $\sum_n \mu_n^2$, the problem is better posed: there are far fewer ways to reduce (or cancel) all the $\mu_n$ at once, which reduces degeneracy and promotes uniqueness of the optimal grid. It will be confirmed in Section 7.2 that this multi-mode optimisation problem is indeed better posed.

We will see below that there are special cases where $f_n(z)$ has a unique sign over the entire domain, leading to a well-posed optimisation problem for a single mode. In that case, as shown in Appendix C, the optimal grid can be derived analytically using the Beltrami identity from the calculus of variations. For $f_n(z) > 0$, it is given by the solution of the ordinary differential equation

$$z'(s) = \left( \frac{K}{f_n(z(s))} \right)^{1/3} \qquad (47)$$

where $K^{1/3} = \int_{-H}^{0} f_n(z)^{1/3} dz$. In addition the minimum value is then

$$\mu_{n,\min} = K = \left( \int_{-H}^{0} f_n(z)^{1/3} dz \right)^3 .$$

## 7.1 | Analytical optimal grids

Below, we study particular cases where $f_n(z)$ has a unique sign, allowing us to derive explicit expressions for the optimal grid.

**The case $N^2 = \text{const}$:**

When $N^2$ is constant, the eigenfunctions and eigenvalues of (7) are explicit and have been given in subsection 5.1. In this case, we have (see (41))

$$\mu_n = \frac{n^4 \pi^4}{12 N_0^2 H^5} \int_0^1 h^3(s) (3\beta - 1 + 3(\beta - 1) \cos(2n\pi z(s))) ds$$

or

$$f_n(z) = \frac{n^4 \pi^4}{12 N_0^2 H^5} (3\beta - 1 + 3(\beta - 1) \cos \frac{2n\pi z}{H}). \qquad (48)$$

It is easy to verify that $f_n(z)$ is positive if $\beta \ge \frac{2}{3}$.

As an example, for $\beta = \frac{2}{3}$, the optimal grid spacing satisfies the differential equation:

$$z'(s) = K \left| \sin \left( \frac{n\pi z(s)}{H} \right) \right|^{-\frac{2}{3}} \quad \text{with} \quad K = \frac{\Gamma\left(\frac{5}{6}\right)}{\sqrt{\pi} \Gamma\left(\frac{4}{3}\right)}. \qquad (49)$$

where $\Gamma$ is the Euler Gamma function.

Also note that, in this case of constant $N^2$, only the value $\beta = 1$ leads to a trivial uniform grid $z'(s) = h(s) = H$ as the optimal solution.
**The case $N^2 \neq \text{const}$:**

For the general non constant case, we have

$$f_n(z) = \frac{\lambda_n}{24H} [(3\beta - 2)(N^2 W_n^2)'' - 2(3\beta - 1)N^2 W_n W_n''] .$$

Using the continuous eigenvalue equation (7), we can rewrite this as

$$f_n(z) = \frac{1}{24H} [(2 - 3\beta)(W_n'' W_n)'' + 2(3\beta - 1)(W_n'')^2]. \quad (50)$$

This expression indicates that enforcing a constant sign is, in general, more challenging than in the constant-$N^2$ case. The only trivial case arises when

$$\beta = \frac{2}{3}, \quad (51)$$

for which the first term vanishes. Then, (50) becomes

$$f_n(z) = \frac{1}{12H} W_n''(z)^2, \quad (52)$$

and $f_n(z)$ is always positive. Using equation (67) from appendix C, the optimal grid spacing is thus given by

$$z'(s) = K (W_n''(z(s)))^{-\frac{2}{3}} \quad \text{with} \quad K = \int_{-H}^{0} (W_n''(z))^{\frac{2}{3}} dz. \quad (53)$$

This result indicates that the optimal grid spacing adapts to the curvature of the eigenfunction $W_n(z)$. Specifically, finer grid resolution is needed where $W_n''(z)$ is large. Obviously, each value of $n$ leads to a different optimal grid. It is therefore interesting to consider the optimal solution when the minimisation is performed over several modes simultaneously. Suppose we want to minimise the sum of the errors of the first $q$ modes, i.e. $\sum_{n=1}^{q} \gamma_n \mu_n$ where $\gamma_n$ are positive weights. The optimal grid spacing is then given by

$$z'(s)^3 \sum_{n=1}^{q} \gamma_n (W_n''(z(s)))^2 = z'(s)^3 \sum_{n=1}^{q} \gamma_n \lambda_n^2 N^4(z(s)) (W_n(z(s)))^2 = \text{constant.} \quad (54)$$

Using the Wentzel-Kramers-Brillouin (WKB) approximation<fnref n="2" />, we have $W_n(z)^2 \sim \frac{1}{N(z)}$. As is known, this approximation is not valid for small $n$, but it may not alter the sum significantly when $q$ is large enough. So that the optimal grid spacing for multiple modes approximately satisfies

$$z'(s)^3 N(z(s))^3 \sum_{n=1}^{q} \gamma_n \lambda_n^2 = \text{constant} \quad \Rightarrow \quad z'(s)^3 N(z(s))^3 = \text{constant}. \quad (55)$$

We remarkably recover the heuristic result that the optimal vertical grid for internal-wave propagation maintains a constant value of $z'(s)N(z(s))$. This coordinate is called the stratification-weighted Charney coordinate and is given

<footnote n="2">²The WKB method provides an analytical framework for approximating eigenvalues and eigenfunctions when exact solutions are not available (e.g. Bender and Orszag, 1978)</footnote>
by

$$s(z) = \frac{1}{C} \int_{-H}^{z} N(z') dz' \quad \text{with} \quad C = \int_{-H}^{0} N(z') dz'. \qquad (56)$$

This particular coordinate transformation has been used in several papers to design vertical grids either for numerical models (?) or to numerically approximate vertical eigenpairs (Early et al., 2020). Here we derive this result rigorously as the optimal grid for minimising the second-order eigenvalue errors, and emphasize that it holds only for the specific discretisation parameter $\beta = \frac{2}{3}$.

## 7.2 | Numerical optimisation of the vertical grid

In the previous section, the admissible range of the parameter $\beta$, to ensure the well-posedness of the optimal grid problem and associated optimal grids, was calculated. While for the general case $N^2 \neq \text{cst}$, only the specific value $\beta = \frac{2}{3}$ ensures a positive sign of the error and gives (53) as the optimal grid expression, in the more specific case of $N^2 = \text{cst}$, a wider range of $\beta (\frac{2}{3} \le \beta \le 1)$ gives admissible optimal grids. We now investigate numerically how these conditions are relevant for different analytical stratification $N^2$.

We know that different parameters can affect the optimal grid, such as the stratification profile $N^2(z)$, the discretisation parameter $\beta$, but also the number of modes considered for the optimisation.

### Spacing parameterisation

To ensure positivity of the spacings and the sum constraint, we first define an auxiliary function $c(s)$ expressed in a B-spline basis as

$$c(s) = \sum_{i=1}^{N_b} c_i B_i(s),$$

where the $c_i$ are the coefficients to be optimized. The $N_b$ B-spline basis functions $B_i(s)$ are defined on a uniform knot vector on $s \in [0, 1]$. The knots (and associated control points) are regularly spaced between $s = 0$ and $s = 1$. The vertical grid spacing is then defined as

$$w_i = \exp(c_i - \max_j c_j), \quad h_i = M H \frac{w_i}{\sum_j w_j}.$$

with M the number of cells. By construction $h_i > 0$ and

$$\sum_i h_i = M H,$$

Varying the number of basis functions $N_b$ directly controls the number of degrees of freedom in the grid parametrisation: fewer basis functions enforce greater smoothness and typically reduce the number of local minima encountered by the optimisation, while more basis functions increase flexibility at the cost of more potential minima. For the numerical experiments presented here, we use $N_b = M/10$ basis functions.

Since we know from the previous analysis that there are likely multiple local minima for most of the values of $\beta$, we perform several independent optimisations, initialised with random guesses for $c_i$, to explore the solution space. For
each optimisation, we use a global optimiser (the Differential Evolution algorithm from the scipy library) to limit the risk of being trapped in local minima.

## Objective function

The objective function to be minimised is the following numerical error:

$$ \mu_{\text{num}}(c_i, N_{\text{eig}}) = \sum_{n=1}^{N_{\text{eig}}} \frac{1}{n^2} \left| \frac{\lambda_n^D - \lambda_n}{\lambda_n} \right| \qquad (57) $$

where $\lambda_n^D$ is the numerically computed eigenvalue using the discretised eigenvalue problem (36) on the grid derived from $c_i, i = 1, \dots, N_b$, and $\lambda_n$ are the eigenvalues of the discrete problem computed at very high resolution. The number of eigenvalues considered in the cost function is denoted by $N_{\text{eig}}$. The scaling factor is $1/n^2$ since we know from previous analysis (see (29)) that, for constant Brunt-Vaisala frequency, the absolute relative error $\left| \frac{\lambda_n^D - \lambda_n}{\lambda_n} \right|$ grows like $\lambda_n \sim n^2$.

## Results

Figures 9 and 10 show the results of the numerical optimisation of the vertical grid for different values of $\beta$, two stratification profiles $N^2(z)$ (constant and variable), and two numbers of eigenvalues considered in the cost function ($N_{\text{eig}} = 1$ (9) and $N_{\text{eig}} = 4$ (10)). For each case, the minimum and maximum values of the cost function obtained from the different independent optimisations are indicated, as well as, for reference, the cost associated with the use of a uniform grid and the cost associated with the WKB grid for variable $N$ (the WKB grid is identical to the uniform grid for constant $N$). For the special value of $\beta = 2/3$, the cost and optimal grids derived from (54) are also shown for reference (green dashed lines).

The variable stratification profile is defined as in Section 5.2 and the number of vertical cells used in the discretisation is $M = 75$. With $H = 4000m$, this leads to a ratio $\sigma/\Delta z = 3$ where $\sigma$ is the width of the stratification perturbation (see Section 5.2). In this case, where the variation of $N$ is only partially resolved, we showed in Section 5.2 that the second-order error approximation remains broadly valid but exhibits non-negligible quantitative discrepancies.

When only one eigenvalue is considered in the cost function (Figure 9), we observe that for $\beta < 2/3$, multiple distinct grids are found that yield very low values (closed to zero) of the cost function. This is consistent with the theoretical analysis presented earlier. For $\beta \ge 2/3$ and in the case of a constant Brunt-Vaisala frequency (top row), a unique optimal grid emerges that minimizes the error. The study done in (7.1) predicts that when $N^2$ is variable, only $\beta = 2/3$ guarantees a well-posed optimisation problem with a unique optimal grid. However, even at $\beta = 2/3$ multiple distinct optimal grids are still found, although they are noticeably less diverse than those obtained for $\beta < 2/3$. This probably stems from the fact that the second-order error approximation is only partially valid for this stratification profile, as discussed above. When four eigenvalues are included in the cost function (Figure 10), the optimisation problem becomes better posed, as anticipated. Even for $\beta < 2/3$, the number of distinct optimal grids is significantly reduced, and a unique optimal grid is obtained for $\beta \ge 2/3$. This highlights that considering multiple modes in the optimisation helps to regularise the problem and promotes uniqueness of the optimal grid. Only for the case $\beta = 1/3$, which leads to the minimum error, do we still observe multiple distinct optimal grids. This gives the possibility of adding additional constraints to the optimisation problem to select a specific optimal grid among this family without altering the minimum error value.

On these figures, the WKB grid defined by (56) is also shown for reference. The WKB grid is independent of the discretization ($\beta$) and is identical to a regular uniform grid when the stratification is constant. We observe that when four modes are considered in the cost function and for $\beta = 2/3$, the numerically optimised grid closely matches the
WKB grid, confirming the theoretical analysis of Section 7.1. For the case $\beta = 1$ which corresponds to the Lorenz discretisation, the optimal grid also approaches the WKB grid when four modes are considered, although the match is not as close as for $\beta = \frac{2}{3}$ (cost WKB: $2.1 \times 10^{-3}$ vs optimal cost: $1.4 \times 10^{-3}$ for variable $N^2$). For other values of $\beta$, the optimal grid deviates significantly from the WKB grid. As an example, for the Charney-Phillips discretisation ($\beta = 0$), the optimal grid has a cost on order of magnitude less than the WKB grid (cost WKB: $2.9 \times 10^{-3}$ vs optimal cost: $3.4 \times 10^{-4}$ for variable $N^2$).

Optimal grids for different Brunt-Väisälä profiles and beta values. Grids optimized for the first eigenvalue.

<figure><img src="image_10.png" /></figure>

**FIGURE 9** Minimisation of the error for the first mode. Optimized grid positions $z(k)$ with $M = 75$ grid points for various values of $\beta$. Each curve corresponds to an independent numerical optimisation of the vertical grid initialized with a different random guess. Constant profile $N^2$ (top row) and variable profile $N^2$ (bottom row). For each case, the minimum and maximum costs for each optimisation results are indicated, as well as, for reference, the cost associated with the use of a uniform grid. In the case $\beta = 2/3$, the optimal grids derived from (53) are also shown (green dashed lines).
Optimal grids for different Brunt-Väisälä profiles and beta values. Grids optimized for the first 4 eigenvalues.

<figure><img src="image_11.png" /></figure>

**FIGURE 10** Minimisation of the error for the fourth modes. Optimized grid positions $z(k)$ with $M = 75$ grid points for various values of $\beta$. Each curve corresponds to an independent numerical optimisation of the vertical grid initialized with a different random guess. Constant profile $N^2$ (top row) and variable profile $N^2$ (bottom row). For each case, the minimum and maximum costs for each optimisation results, as well as, for reference, the cost associated with the use of a uniform grid are indicated.

8 | CONCLUSIONS AND PERSPECTIVES

Making optimal choices for the vertical discretization of primitive ocean equations remains a challenging problem. In this study, we address this issue by considering a simpler yet fundamental problem: the propagation of internal gravity waves. After recalling the continuous theory of vertical mode decomposition, we propose a general discrete framework that unifies discretisation choices (e.g. staggering, position of vertical levels), energy conservation and the associated discrete Sturm-Liouville problem. This formalism can be applied to a wide range of discretisation techniques (including finite-volume and finite-element methods) and to different orders of accuracy. In particular, we give classical examples commonly used in the oceanic and atmospheric communities, namely second-order finite-volume methods on staggered grids (Charney and modified Lorenz grid). For the Charney-Phillips setting, we introduce a $\beta$-indexed family of discretisations which differs from the classical Charney-Phillips scheme by a single parameter $\beta$ involved in the discrete computation of the vertical hydrostatic balance. $\beta = 0$ and $\beta = 1$ correspond to the classical Charney-Phillips and Lorenz schemes, respectively. The leading-order truncation error depends jointly on the vertical grid and on the discretisation parameter $\beta$. This naturally raises the question of whether an optimal discretization through an appropriate choice of the parameter $\beta$ can be identified, and whether an optimal vertical grid can be derived so as to minimize the leading-order eigenvalue error. The main results of the present paper can be summarized as follows:
* Discretization scheme with $\beta = \frac{1}{3}$ can significantly reduces truncation errors (based on error compensation for the discrete eigenvalues systems), yielding fourth-order convergence on uniform grids.

* For a single eigenvalue, we provide a rigorous characterization of the optimal vertical grid in the general case when $\beta = \frac{2}{3}$, which is the only choice which leads a unique solution for the optimization problem. This optimal grid is shown to be a function of the corresponding vertical eigenmode.

* When considering a large number of eigenvalues, the WKB grid (also known as the Charney grid, see Robey and Grooms (2024)) appears to be a near-optimal vertical grid for the same value of $\beta = \frac{2}{3}$, while its performance deteriorates for other choices of the discretization parameter, in particular for the classical Charney-Phillips scheme ($\beta = 0$).

These theoretical results are validated numerically and offer practical guidance for the design of vertical discretisations in ocean models, while motivating further exploration of grid optimization strategies.

Obviously, the present study is a first step towards the design of optimal vertical discretizations for ocean models. The next step would be to extend the present analysis to more complex physical problems, such as the full primitive equations, which would require to account for additional sources of numerical error. In particular those associated with the horizontal pressure gradient as well as the representation of tracer transport and mixing (link with the relative positioning between the vertical coordinate surfaces and isopycnal levels) would need to be considered. Optimization techniques that account for multiple sources of numerical error already exist in the literature Burchard and Beckers (2004); Hofmeister et al. (2010); Kavčič and Thuburn (2018). They rely on numerous user defined parameters that need to be adapted to the particular problem to solve. A particularly interesting avenue would be to investigate, as demonstrated in the present study, how optimal vertical grids depend on the choice of discretisation schemes and how the choice of the discretisation scheme impacts the uniqueness of the resulting optimal grid.

Beyond seeking for a more general optimal vertical coordinate, we stress again that under certain circumstances stated above, a simple WKB based vertical grid could be a valuable alternative. Moreover, the $\beta$-indexed schemes introduced in this paper appear particularly well suited for implementation in realistic ocean models. It does indeed only require a light modification of the discretisation of the hydrostatic balance (the derivation for the classical Lorenz grid remains to be done though).
A | DISCRETE ENERGY CONSERVATION: A CONSTRAINT ON MATRICES

The discrete kinetic energy is defined as

$$ \tilde{E}_C(t) = \int_0^L \frac{1}{2} \mathbf{U}(x,t)^T \mathcal{A}_U \mathbf{U}(x,t) = \int_0^L \frac{1}{2} \langle \mathbf{U}, \mathbf{U} \rangle_{\mathcal{A}_U} $$

where $\mathcal{A}_U$ is a positive-definite matrix. The time derivative of this discrete kinetic energy is then

$$ \partial_t \tilde{E}_C(x, t) = \int_0^L \mathbf{U}(x, t)^T \mathcal{A}_U \partial_t \mathbf{U}(x, t), $$

and, using system (9), we get

$$ \begin{aligned} \partial_t \tilde{E}_C &= \int_0^L \mathbf{U}^T \mathcal{A}_U \partial_t \mathbf{U} \\ &= - \int_0^L \mathbf{U}^T \mathcal{A}_P \partial_x \mathbf{P} \\ &= \int_0^L \partial_x \mathbf{U}^T \mathcal{A}_P \mathbf{P} \\ &= - \int_0^L \mathbf{W}^T \mathcal{D}_W^T \mathcal{A}_U^{-1} \mathcal{A}_P \mathbf{P} \end{aligned} \quad (58) $$

where the third line follows from assuming closed or periodic boundary conditions in the horizontal direction. The potential energy is defined as

$$ \tilde{E}_P(t) = \int_0^L \frac{1}{2} \mathbf{B}(x,t)^T \mathcal{A}_P \mathbf{B}(x,t) = \int_0^L \frac{1}{2} \langle \mathbf{B}, \mathbf{B} \rangle_{\mathcal{A}_B} $$

where $\mathcal{A}_B$ is also a positive-definite matrix. Therefore

$$ \partial_t \tilde{E}_P(t) = \int_0^L \mathbf{B}(x,t)^T \mathcal{A}_B \partial_t \mathbf{B}(x,t), $$

and we get

$$ \begin{aligned} \partial_t \tilde{E}_P(t) &= \int_0^L \mathbf{B}^T \mathcal{A}_B \partial_t \mathbf{B} \\ &= - \int_0^L \mathbf{B}^T \mathcal{A}_W \mathbf{W} \\ &= - \int_0^L \mathbf{W}^T \mathcal{A}_W^T \mathbf{B} \\ &= - \int_0^L \mathbf{W}^T \mathcal{A}_W^T \mathcal{A}_B^{-1} \mathcal{D}_P \mathbf{P} \end{aligned} \quad (59) $$

Using (58) and (59), the time evolution of the discrete energy reads:

$$ \partial_t (\tilde{E}_C + \tilde{E}_P) = - \int_0^L \mathbf{W}^T (\mathcal{D}_W^T \mathcal{A}_U^{-1} \mathcal{A}_P + \mathcal{A}_W^T \mathcal{A}_B^{-1} \mathcal{D}_P) \mathbf{P} $$
In this paper, for all considered discretisations the W points are located at cell interfaces and we have $W_0 = W_{N_W-1} = 0$ (rigid-lid boundary conditions). Therefore the discrete energy is conserved, $\partial_t (\tilde{E}_C + \tilde{E}_P) = 0$, if the matrix $D_W^T A_U^{-1} A_P + A_W^T A_B^{-1} D_P$ has nonzero elements only in its first and last rows. This can be expressed using a projection matrix $P_0$ which removes the first and last rows of a vector of size $N_W$, $P_0 = I_{N_W} - (e_0 e_0^T + e_{N_W-1} e_{N_W-1}^T)$. The discrete energy is conserved if

$$P_0 (D_W^T A_U^{-1} A_P + A_W^T A_B^{-1} D_P) = 0.$$

## B | PERTURBATION ANALYSIS OF THE EIGENVALUE ERROR

The vertical computational grid is defined by a generalised vertical coordinate transformation $z(s), 0 \le s \le 1$. The grid points are defined by

$$z_j = z(s_j) \quad \text{with} \quad s_j = j \Delta s \quad \text{and} \quad \Delta s = \frac{1}{M}, \quad j = 0, \dots, M \qquad (60)$$

where $M+1$ is the number of grid points and $s_j$ are regularly spaced points. The transformation $z(s)$ is such that $z(0) = -H$ and $z(1) = 0$ and that the scale factor $h(s) = z'(s)$ is strictly positive.

Using these notations, $W_j$ is the evaluation of the discrete solution, denoted $W_d$, at $z_j$ points: $W_j = W_d(z_j) = W_d(z(s_j))$. We also note $\lambda_d$ the discrete eigenvalue.

### B.1 | Modified equations associated to the discretizations

For small $\Delta s$, the left hand sides of the discrete equation (36) can be expanded in Taylor series around $s$, with $s = s_j$. One has

$$\frac{W_{j+1} - W_j}{\Delta z_{j+1}} - \frac{W_j - W_{j-1}}{\Delta z_j} = z' \Delta s W_d'' + \frac{1}{12} (\Delta s)^3 (2 z^{(3)} W_d'' + 4 z' z'' W_d^{(3)} + z'^3 W_d^{(4)}) + o((\Delta s)^3)$$

High order derivatives of $W_d$ can be expressed using the continuous equation (7). Indeed, one has

$$\begin{aligned} W_d'' &\approx -\lambda_d N^2 W_d \\ W_d^{(3)} &\approx -\lambda_d [N^2 W_d]' \\ W_d^{(4)} &\approx -\lambda_d [N^2 W_d]'' \end{aligned}$$

Using these expressions, one can write

$$\frac{W_{j+1} - W_j}{\Delta z_{j+1}} - \frac{W_j - W_{j-1}}{\Delta z_j} = z' \Delta s W_d'' - \frac{\lambda_d}{24} (\Delta s)^3 \left( 2 z'^3 [N^2 W_d]'' + 8 z' z'' [N^2 W_d]' + 4 z^{(3)} N^2 W_d \right) + o((\Delta s)^3)$$
Starting from the right hand side of discrete Eq. (36), a Taylor expansion gives

$$z' \Delta s \lambda_d N^2 W_d + \frac{\lambda_d}{24} (\Delta s)^3 \left( 3\beta z'^3 \left[ 2(N^2)' W_d' + W_d(N^2)'' + 2N^2 W_d'' \right] + 12\beta z' z'' (N^2 W_d)' + 4z^{(3)} N^2 W_d \right) + o((\Delta s)^3)$$

The modified equation associated with the discrete equations (36) is then

$$W_d'' + \lambda_d N^2 W_d + \frac{1}{24} (\Delta s)^2 \lambda_d \left( z'^2 \left[ (3\beta - 2)(2(N^2)' W_d' + W_d(N^2)'' + 2(3\beta - 1)N^2 W_d'') \right] \right. \\ \left. + 4(3\beta - 2)z''(N^2 W_d)' \right) = o((\Delta s)^2) \quad (61)$$

## B.2 | Perturbation analysis

One can now write the above modified equations in the general form

$$(A + (\Delta s)^2 \delta A) W_d + \lambda_d (B + (\Delta s)^2 \delta B) W_d + o((\Delta s)^2) = 0$$

where A and B are two differential operators and $(\Delta s)^2 \delta A$ and $(\Delta s)^2 \delta B$ are "small" perturbations of these operators. The eigenpairs of the continuous equation satisfies

$$AW + \lambda BW = 0, \quad \text{with } AW = W'', BW = N^2 W$$

It is well-known Rellich (1969); Duval (2022) that $\mu_n$ and $V_n$, the perturbations of eigenvalues and eigenvectors of index $n$ defined by

$$\lambda_n^d = \lambda_n + (\Delta s)^2 \mu_n + o((\Delta s)^2) \quad \text{and} \quad W_n^d = W_n + (\Delta s)^2 V_n + o((\Delta s)^2)$$

can be written in terms of perturbation of the two initial operators and the unperturbed eigenvectors and eigenvalues as

$$\mu_n = - \frac{\langle W_n, (\delta A + \lambda_n \delta B) W_n \rangle}{\langle W_n, B W_n \rangle} \quad (62)$$

and

$$V_n(z) = -\frac{1}{2} \frac{\langle W_n, \delta B W_n \rangle}{\langle W_n, B W_n \rangle} - \sum_{m \neq n} \frac{\langle W_m, (\delta A + \lambda_m \delta B) W_n \rangle}{(\lambda_n - \lambda_m) \langle W_m, B W_m \rangle} \quad (63)$$

If $W_n$ are chosen to be orthonormal for the scalar product associated with B, i.e. $\langle W_n, B W_m \rangle = \frac{1}{H} \int_{-H}^{0} N^2(z) W_n(z) W_m(z) dz =$ $\delta_{nm}$, then the above expressions simplify to

$$\mu_n = - \langle W_n, (\delta A + \lambda_n \delta B) W_n \rangle \quad (64)$$
and

$$V_n(z) = -\frac{1}{2} \langle W_n, \delta B W_n \rangle W_n(z) - \sum_{m \neq n} \frac{\langle W_m, (\delta A + \lambda_n \delta B) W_n \rangle}{(\lambda_n - \lambda_m)} W_m(z) \quad (65)$$

## Application to the generalised Charney-Phillips discretisations

For the generalised Charney-Phillips discretisations, (61), one has

$$\delta A(W_d) = 0 \\ \delta B(W_d) = \frac{1}{24} \left[ z'^2 \left[ (3\beta - 2)(2(N^2)' W_d' + W_d(N^2)')' + 2(3\beta - 1)N^2 W_d'' \right] + 4(3\beta - 2)z''(N^2 W_d)' \right]$$

In the s coordinate system, the scalar product writes

$$\langle W_n, W_m \rangle = \frac{1}{H} \int_0^1 W_n(z(s)) W_m(z(s)) z'(s) ds$$

thus $\mu_n$ can be written as

$$\mu_n = -\frac{\lambda_n}{H} \int_0^1 W_n(z(s)) \delta B W_n(z(s)) z'(s) ds.$$

Knowing the $\delta BW$ expression, the integrand is composed of terms of the forms $z'(s)^3 f_1(z(s))$ and $z'(s)z''(s)f_2(z(s))$ where $f_1$ and $f_2$ depend on $N^2$ and $W_n$ and their derivatives:

$$\mu_n = -\frac{\lambda_n}{H} \int_0^1 \left( z'(s)^3 f_1(z(s)) + z'(s) z''(s) f_2(z(s)) \right) ds$$

We use the following relation which involves integration by parts:

$$\int_0^1 z'(s) z''(s) f_2(z(s)) ds = \int_0^1 \frac{1}{2} (z'(s)^2)' f_2(z(s)) ds \\ = - \int_0^1 \frac{1}{2} z'(s)^2 \frac{d}{ds} f_2(z) ds + \left[ \frac{1}{2} z'(s)^2 f_2(z(s)) \right] \\ = - \int_0^1 \frac{1}{2} z'(s)^3 f_2'(z(s)) ds$$

since the boundary terms vanish because $f_2(-H) = f_2(0) = 0$.

Using these relations, one can write $\mu_n$ as

$$\mu_n = -\frac{\lambda_n}{H} \int_0^1 z'(s)^3 \left( f_1 - \frac{1}{2} f_2' \right) ds.$$

This leads to the following expressions of $\mu_n$:

$$\mu_n = -\frac{\lambda_n}{24H} \int_0^1 z'(s)^3 \left[ (3\beta - 2)(N^2)' (W_n^2)' + W_n^2(N^2)'' + 2(3\beta - 1)N^2 W_n'' W_n - 2(3\beta - 2)((N^2 W_n)' W_n)' \right] ds$$
To further simplify these integrals, we identify terms that can be written as total derivatives, such as $(N^2W^2)''$. This is useful because, for the special case of a uniform grid, $z'(s)$ constant, the integrals of these total derivatives vanish due to the boundary conditions for $W$. The finally expressions of $\mu_n$ are

$$ \mu_n = \frac{\lambda_n}{24H} \int_0^1 z'(s)^3 [(3\beta - 2)(N^2 W_n^2)'' - 2(3\beta - 1)N^2 W_n W_n'''] ds \quad (66) $$

C | OPTIMAL GRID SPACING VIA THE BELTRAMI IDENTITY

We consider the functional

$$ J[z] = \int_0^1 z'(s)^3 f(z(s)) ds, $$

where $f(z)$ has a constant sign and $z'(s) > 0$. The Beltrami identity states that

$$ \mathcal{L} - z' \frac{\partial \mathcal{L}}{\partial z'} = \text{constant}, $$

where $\mathcal{L}(z, z') = z'(s)^3 f(z(s))$. We compute

$$ \frac{\partial \mathcal{L}}{\partial z'} = 3z'(s)^2 f(z(s)), $$

so that the Beltrami identity becomes

$$ \mathcal{L} - z' \frac{\partial \mathcal{L}}{\partial z'} = z'(s)^3 f(z(s)) - 3z'(s)^3 f(z(s)) = -2z'(s)^3 f(z(s)) = C, $$

where $C$ is a constant. Thus,

$$ z'(s)^3 f(z(s)) = -\frac{C}{2}. $$

Assuming $f(z(s)) > 0$, we set $K = -C/2 > 0$ and obtain

$$ z'(s) = \left( \frac{K}{f(z(s))} \right)^{1/3}. \quad (67) $$

The constant $K$ is determined by the boundary conditions on $z(s)$:

$$ \int_0^1 z'(s) ds = z(1) - z(0) = H. $$

which leads to

$$ K^{1/3} = \int_{-H}^{0} f(z)^{1/3} dz. $$
This gives an ordinary differential equation for $z(s)$, which can be solved numerically for a given $f(z)$. The minimum value of the functional is then

$$J_{\min} = K = \left( \int_{-H}^{0} f(z)^{1/3} dz \right)^{3}.$$

The fact that the optimal solution satisfies $(z')^3 f(z) = \text{const}$ is not accidental. Because the Lagrangian $L(z, z') = (z')^3 f(z)$ is homogeneous of degree three in $z'$ and does not depend explicitly on the independent variable, the Beltrami identity reduces to $(1 - 3)L = \text{const}$. Hence $L$ itself must be constant along any extremal.

## acknowledgements

This research was funded by the French National Research Agency (ANR) under project ANR-23-CE56-0006-01 (MOTIONS). This work was supported by funding from the French government through the National Research Agency as part of France 2030, under reference ANR-22-POCE-0003.

## conflict of interest

The authors declare no conflict of interest.

## references

Arakawa, A. (1988) *Finite-Difference Methods in Climate Modeling*, 79–168. Dordrecht: Springer Netherlands.

Arakawa, A. and Konor, C. S. (1996) Vertical differencing of the primitive equations based on the charney-phillips grid in hybrid and sigma-p vertical coordinates. *Mon. Weather Rev.*, **124**, 511 – 528.

Arbic, B. K., Richman, J. G., Shriver, J. F., Timko, P. G., Metzger, E. J. and Wallcraft, A. J. (2012) *Oceanography*. URL: https://doi.org/10.5670/oceanog.2012.38.

Baer, F. and Ji, M. (1989) Optimal vertical discretization for atmospheric models. *Mon. Weather Rev.*, **117**, 391 – 406.

Barkan, R., Srinivasan, K., Yang, L., McWilliams, J. C., Gula, J. and Vic, C. (2021) Oceanic mesoscale eddy depletion catalyzed by internal waves. *Geophys. Res. Lett.*, **48**, e2021GL094376. URL: https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2021GL094376. E2021GL094376 2021GL094376.

Bell, M. J. and White, A. A. (2017) Analytical approximations to spurious short-wave baroclinic instabilities in ocean models. *Ocean Modell.*, **118**, 31–40.

Bender, C. M. and Orszag, S. A. (1978) *Advanced Mathematical Methods for Scientists and Engineers I: Asymptotic Methods and Perturbation Theory*. New York: McGraw-Hill.

Blayo, E. and Debreu, L. (2005) Revisiting open boundary conditions from the point of view of characteristic variables. *Ocean Modell.*, **9**, 231 – 252.

Buijsman, M. C., Stephenson, G. R., Ansong, J. K., Arbic, B. K., Green, J. M., Richman, J. G., Shriver, J. F., Vic, C., Wallcraft, A. J. and Zhao, Z. (2020) On the interplay between horizontal resolution and wave drag and their effect on tidal baroclinic mode waves in realistic global ocean simulations. *Ocean Modell.*, **152**, 101656. URL: https://www.sciencedirect.com/science/article/pii/S146350032030158X.
Burchard, H. and Beckers, J.-M. (2004) Non-uniform adaptive vertical grids in one-dimensional numerical ocean models. *Ocean Modelling*, **6**, 51–81. URL: https://www.sciencedirect.com/science/article/pii/S1463500302000604.

Charney, J. G. and Phillips, N. A. (1953) Numerical integration of the quasi-geostrophic equations for barotropic and simple baroclinic flows. *J. Atmos. Sci.*, **10**, 71 – 99.

Demange, J. (2014) *Schémas numériques d'advection et de propagation d'ondes de gravité dans les modèles de circulation océanique*. Theses, Université de Grenoble. URL: https://theses.hal.science/tel-01104703.

Duval, E. (2022) *Vers le couplage de modèles de circulation océanique hydrostatique et non hydrostatique*. Ph.D. thesis, Université Grenoble Alpes. URL: http://www.theses.fr/2022GRALM030.

Early, J. J., Lelong, M. P. and Smith, K. S. (2020) Fast and accurate computation of vertical modes. *Journal of Advances in Modeling Earth Systems*, **12**, e2019MS001939. URL: https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2019MS001939. E2019MS001939 10.1029/2019MS001939.

Gill, A. E. (1982) *Atmosphere-Ocean Dynamics*. International geophysics series.

Griffies, S. M., Boning, C., Bryan, F. O., Chassignet, E. P., Gerdes, R., Hasumi, H., Hirst, A., Treguier, A.-M. and Webb, D. (2000) Developments in ocean climate modelling. *Ocean Modell.*, **2**, 123–192.

Hallberg, R. (2013) Using a resolution function to regulate parameterizations of oceanic mesoscale eddy effects. *Ocean Modell.*, **72**, 92–103.

Hofmeister, R., Burchard, H. and Beckers, J.-M. (2010) Non-uniform adaptive vertical grids for 3d numerical ocean models. *Ocean Modelling*, **33**, 70–86. URL: https://www.sciencedirect.com/science/article/pii/S1463500309002248.

Holdaway, D., Thuburn, J. and Wood, N. (2013) Comparison of lorenz and charney-phillips vertical discretisations for dynamics-boundary layer coupling. part i: Steady states. *Quart. J. Roy. Meteorol. Soc.*, **139**, 1073–1086. URL: https://rmets.onlinelibrary.wiley.com/doi/abs/10.1002/qj.2016.

Kavčič, I. and Thuburn, J. (2018) A lagrangian vertical coordinate version of the endgame dynamical core. part i: Formulation, remapping strategies, and robustness. *Quarterly Journal of the Royal Meteorological Society*, **144**, 1649–1666.

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**, 3777–3788.

Leslie, L. and Purser, R. (1992) A comparative study of the performance of various vertical discretization schemes. *Meteorol. Atmos. Phys.*, **50**, 61–73.

Lorenz, E. N. (1960) Energy and numerical weather prediction. *Tellus*, **12**, 364–373.

Nikurashin, M. and Ferrari, R. (2013) Overturning circulation driven by breaking internal waves in the deep ocean. *Geophys. Res. Lett.*, **40**, 3133–3137. URL: http://dx.doi.org/10.1002/grl.50542.

Nycander, J. (2005) Generation of internal waves in the deep ocean by tides. *J. Geophys. Res.: Oceans*, **110**, URL: https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2004JC002487.

Rellich, F. (1969) *Perturbation Theory of Eigenvalue Problems*. New York: Gordon and Breach. Translated from the German original, Springer-Verlag, 1969.

Robey, R. and Grooms, I. (2024) Continuous and discrete baroclinic modes in continuously varying stratification. *SIAM Journal on Applied Mathematics*, **84**, 2502–2521. URL: https://doi.org/10.1137/24M1648181.

Saint-Laurent, L., Alford, M. H. and Paluszkiewicz, T. (2012) An introduction to the special issue on internal waves. *Oceanography*. URL: https://doi.org/10.5670/oceanog.2012.37.
Smith, S. G. L. and Young, W. R. (2002) Conversion of the barotropic tide. *Journal of Physical Oceanography*, **32**, 1554 – 1566.
URL: https://journals.ametsoc.org/view/journals/phoc/32/5/1520-0485_2002_032_1554_cotbt_2.0.co_2.xml.

Stewart, K., Hogg, A., Griffies, S., Heerdegen, A., Ward, M., Spence, P. and England, M. (2017) Vertical resolution of baroclinic modes in global ocean models. *Ocean Modell.*, **113**, 50–65.

Thuburn, J. (2006) Vertical discretizations giving optimal representation of normal modes: Sensitivity to the form of the pressure-gradient term. *Quart. J. Roy. Meteorol. Soc.*, **132**, 2809–2825.

Thuburn, J. and Woollings, T. (2005) Vertical discretizations for compressible euler equation atmospheric models giving optimal representation of normal modes. *J. Comp. Phys.*, **203**, 386–404.

Tokioka, T. (1978) Some considerations on vertical differencing. *J. Meteorol. Soc. Jpn.*, **56**, 98–111.

Ullrich, P. A. (2014) Understanding the treatment of waves in atmospheric models. part 1: The shortest resolved waves of the 1d linearized shallow-water equations. *Quart. J. Roy. Meteorol. Soc.*, **140**, 1426–1440.