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

HAL
open science

# Dispersion analysis of compatible Galerkin schemes on quadrilaterals for shallow water models

Christopher Eldred, Daniel Le Roux

► To cite this version:

Christopher Eldred, Daniel Le Roux. Dispersion analysis of compatible Galerkin schemes on quadrilaterals for shallow water models. Journal of Computational Physics, 2019, 387, pp.539-568. 10.1016/j.jcp.2019.02.009. hal-01916382v2

HAL Id: hal-01916382

https://hal.science/hal-01916382v2

Submitted on 14 Feb 2019

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.
Dispersion analysis of compatible Galerkin schemes on
quadrilaterals for shallow water models

Christopher Eldred<sup>a,*</sup>, Daniel Y. Le Roux<sup>b</sup>

<sup>a</sup>Univ. Grenoble Alpes, Inria, CNRS, Grenoble INP, LJK, 38000 Grenoble, France

<sup>b</sup>Université de Lyon, CNRS, Université Lyon 1, Institut Camille Jordan, 43, blvd du 11 novembre 1918, 69622 Villeurbanne Cedex, France

**Abstract**

Compatible Galerkin methods (the Galerkin analogue of an Arakawa C-Grid) are growing in popularity for simulating geophysical fluid flows, due to their desirable characteristics, including but not limited to: energy conservation, higher-order accuracy, steady geostrophic modes and the absence of spurious stationary modes, such as pressure modes. However, these characteristics still do not guarantee good wave dispersion properties. In this work, we study the dispersion properties of two compatible Galerkin families for the 2D linear rotating shallow water equations on quadrilaterals: the $Q_n^-\Lambda^k$ family from finite element exterior calculus and the newly developed $MGD_n$ family. These families are the extensions to quadrilaterals of the $P_n^C - P_{n-1}^{DG}$ and $GD_n - DGD_{n-1}$ pairs, respectively, studied for the 1D linear shallow water equations in [13]. A major finding from that paper was that the $P_n^C - P_{n-1}^{DG}$ pair has spectral gaps for $n \ge 2$ and the $GD_n - DGD_{n-1}$ does not. These spectral gaps are non-dimensional wavenumbers where the dispersion relationship is double-valued, and lead to anomalous dispersion and noise in numerical simulations. On quadrilaterals, previous work [24, 26] on the $Q_n^-\Lambda^k$ family for inertia waves with $n = 2$ and for gravity waves for arbitrary $n$ has indicated the presence of spectral gaps, in the form of line discontinuities. The investigation of these gaps for the $Q_n^-\Lambda^k$ family is extended in this paper to inertia-gravity waves for arbitrary $n$, including plots of the dispersion relationship for $n = 2$ when using the lumping developed in [26, 35] that eliminates the spectral gaps. Additionally, the $MGD_n$ family is studied (including the use of reduced quadrature), which is found to be free of spectral gaps. For both families asymptotic convergence rates are established, effective resolutions determined and plots of the dispersion relationships for a range of $n$ and Rossby radii are shown. Finally, a pair of numerical simulations are run to investigate the consequences of the spectral gaps and highlight the main differences between the two families.

**Keywords:**

dispersion relationship, shallow water equations, geophysical fluid dynamics, dynamical core, mixed finite elements, finite element exterior calculus, mimetic Galerkin differences

<footnote>*Corresponding Author: Christopher Eldred (christopher.eldred@inria.fr)</footnote>
# 1. Introduction

Compact Galerkin methods (such as finite element, discontinuous Galerkin and spectral element methods) are finding increasing use in numerical models of geophysical fluid flows. Of particular interest amongst these approaches are compatible Galerkin methods [14], which are a generalization of compatible finite element methods [4, 5, 9, 29] to a more general definition of an element. A compatible Galerkin family is a realization of a discrete deRham complex, and are the compact Galerkin analogue of Arakawa C-grid finite difference schemes (also known as staggered grid methods or Marker and Cell methods). A commonly used compatible Galerkin family on quadrilaterals is the $Q_n^-\Lambda^k$ compatible finite element family from finite element exterior calculus. Like C-grid schemes, they have many desirable properties, including but not limited to energy conservation, steady geostrophic modes and absence of spurious pressure modes or inertial modes; however, all known examples do have an f-mode (also known as a CD or Coriolis mode) [9, 19] due to the discrete Coriolis matrix being rank-deficient. Unfortunately, these desirable properties are not sufficient to ensure that compatible Galerkin methods have good wave dispersion properties. In particular, two types of problems in the dispersion relationship are known to arise: spurious branches and spectral gaps. Considering the linear shallow water equations in 1D or 2D with constant Coriolis parameter $f$, there are three branches of the dispersion relationship: one set of geostrophic modes (these become Rossby modes when $f$ is variable), and two pairs of inertia-gravity waves with opposite sign. Ideally a numerical model would also possess only three branches in the dispersion relationship, but certain discretizations will have extra branches of either Rossby or inertia-gravity (IGW) modes [9, 11, 30, 36]. Typically, these occur on non-quadrilateral grids due to an imbalance between the number of degrees of freedom in the wind field and the mass field, although the $S_n\Lambda^k$ family on quadrilaterals has spurious Rossby modes. For example, the $RT_0 - P_0^{DG}$ and $BDM_1 - P_0^{DG}$ pairs on triangles have spurious branches of IGW and Rossby waves, respectively. The other problem that can arise are spectral gaps. Spectral gaps are non-dimensional wavenumbers where the dispersion relationship is multivalued and the group velocity goes to zero. These gaps seem to be a general feature of high order finite element discretizations [1, 2, 27], not just compatible finite elements. Such gaps, which typically take the form of points in 1D and lines in 2D, result in a piecewise continuous dispersion relationship. A comprehensive understanding of spurious modes and spectral gaps is still lacking, especially for triangular elements. However, it is known that certain choices of spaces, such as quadrilateral elements formed by a tensor product construction (both the $Q_n^-\Lambda^k$ and $MGD_n$ families are built like this), can ensure a lack of spurious branches [9]. In this work we will show that the $MGD_n$ family additionally avoids spectral gaps in the dispersion relationship for any $n$, unlike the $Q_n^-\Lambda^k$ family which has gaps for $n \ge 2$.

Dispersion relationships have been studied since the beginning of geophysical fluid modeling [3], initially using finite differences. Within the class of compatible Galerkin methods, the $RT_0 - P_0^{DG}$ and/or $BDM_1 - P_0^{DG}$ pairs on triangles were studied in [19, 20, 21, 33, 34]. These are the lowest-order members of the corresponding $P_r^-\Lambda^k$ ($RT_0$) and $P_r\Lambda^k$ ($BDM_1$) families on triangles from finite element exterior calculus [5]. However, both of these families
are known to have spurious branches [9], and the presence or absence of spectral gaps is unclear, although the results from the $P_n^C - P_{n-1}^{DG}$ pair in 1D strongly suggest the presence of gaps. A low-order compatible finite element pair on arbitrary polygonal elements (termed compound finite elements, introduced in [37]) were studied for quadrilaterals and hexagons in [28], and compared to the corresponding C-grid finite difference schemes [30, 32, 36]. Both pairs were found to be similar to the finite-difference schemes, and the hexagonal pair possesses a similar spurious branch of Rossby waves. For higher order elements, a spectral gap was found for the $P_2^C - P_1^{DG}$ pair [10] in 1D. A solution to this gap for the 1D gravity wave equations, obtained through partial lumping of the velocity mass matrix, is given in [35], and was extended to the 2D shallow water equations in [26]. However, [13] showed that this form of lumping does not work for $n \ge 3$, and for $n = 2$ it leads to a reduction in the order of convergence. The dispersion relationship for a simplified version of the fully compressible Euler equations, in a doubly periodic slice domain using 2D versions of the $P_2^C - P_1^{DG}$ and $P_1^C - P_0^{DG}$ pairs was studied in [25], where it was found that a vertically continuous and horizontal discontinuous space (the finite element analogue of a C-grid staggering in the horizontal and a Charney-Phillips staggering in the vertical) gave the best results. The pure inertia and pure gravity limits of the $P_n^C - P_{n-1}^{DG}$ pair in 1D were studied in [24], along with 2D gravity waves using the $Q_n^- \Lambda^k$ family. Finally, the $P_n^C - P_{n-1}^{DG}$ and $GD_n - DGD_{n-1}$ pairs were studied for inertia-gravity waves in 1D in [13]. Both [24] and [13] showed that the $P_n^C - P_{n-1}^{DG}$ pair suffers from spectral gaps when $n \ge 2$, and also has a maximal frequency that increases with $n$. Due to the spectral gaps, anomalous wave dispersion occurs when there are wave packets with significant energy close to the frequency of the gaps, and the maximal frequency increase is expected to cause a decrease in maximum stable time step for explicit timestepping schemes due to the CFL limit. It would be possible to partially address these issues through the addition of numerical dissipation, but this is undesirable since such dissipation is not physically motivated. In contrast, the $GD_n - DGD_{n-1}$ pair was found [13] to be free of spectral gaps for any $n$ and has a maximal frequency that does not increase with $n$.

In this work, we extend the study started in [13] from the 1D shallow water equations to the 2D shallow water equations on quadrilaterals, and investigate the dispersion properties of the $Q_n^- \Lambda^k$ and $MGD_n$ families, which are the 2D extensions of the $P_n^C - P_{n-1}^{DG}$ and $GD_n - DGD_{n-1}$ pairs. Additionally, we clarify a partially erroneous claim found at the end of [24], about the expression of the 2D inertia-gravity wave dispersion relationship in terms of the 1D inertia wave and gravity wave dispersion relationships.

The remainder of this paper is structured as follows. In Section 2 we introduce the model problem under consideration: the 2D linear shallow water equations. A general discretization of this problem by compatible Galerkin methods is developed in Section 3, along with the corresponding Fourier analysis to compute the discrete dispersion relationship. The tensor product structure of the spaces is exploited in Section 4 to express the 2D matrices in terms of Kronecker products of the 1D matrices. Using this tensor product structure, it is also shown that the 2D inertia wave (IW) and gravity wave (GW) dispersion relationships can be expressed in terms of the 1D IW and GW dispersion relationships for any $n$. However,
this is not true for 2D inertia-gravity waves (IGW) for compatible Galerkin discretizations that have multiple degrees of freedom per geometric entity (i.e. when $d \ge 2$), such as the $Q_n^- \Lambda^k$ family when $n \ge 2$). Nevertheless, using the general approach from Section 3, and exploiting the tensor product structure discussed in Section 4, the asymptotics and dispersion relationships for the $Q_n^- \Lambda^k$ and $MGD_n$ families are studied in Section 5 and Section 6. This includes asymptotics and a dispersion plot for the $Q_2^- \Lambda^k$ family using the mass lumping from [26, 35] that eliminates the spectral gaps. It is shown that the $Q_n^- \Lambda^k$ family suffers from spectral gaps when $n \ge 2$, in the form of line discontinuities, while the $MGD_n$ family is free of gaps for any $n$. Two test cases are run in Section 7 to investigate the impact of these gaps. This section also includes a study of the effective resolution for the two families. Finally, some conclusions and possible future work are discussed in Section 8. The work in this paper represents a starting point for the analysis of the dispersion properties of compatible Galerkin schemes applied to geophysical fluid models: in particular, an extension to the fully compressible Euler equations (and other relevant geophysical fluid models) in three-dimensional domains with boundaries and variable Coriolis force and the incorporation of time discretization remains to be done.

## 2. Model Problem

We start with the 2D, inviscid, rotating (with a constant Coriolis parameter $f$) shallow water equations in a doubly periodic domain with a flat bottom and uniform gravitational constant $g$. These are linearized about a state of rest ($u = v = 0$) with a constant fluid depth $H$. Using a Cartesian coordinate system the linear shallow water equations are written [22] as

$$ \frac{\partial u}{\partial t} - f v + g \frac{\partial \eta}{\partial x} = 0, \qquad (1) $$

$$ \frac{\partial v}{\partial t} + f u + g \frac{\partial \eta}{\partial y} = 0, \qquad (2) $$

$$ \frac{\partial \eta}{\partial t} + H \left( \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} \right) = 0, \qquad (3) $$

with relative velocity $v(x, y, t) = (u, v)$ and surface height $\eta(x, y, t)$. Since (1) - (3) are a linear system with constant coefficients in a doubly periodic domain, Fourier analysis can be used to determine the linear modes. Introducing temporal frequency $\omega$ and spatial wavenumbers $k$ and $l$, a single Fourier mode is written as

$$ \eta = \hat{\eta} e^{-i\omega t} e^{ikx} e^{ily} \qquad u = \hat{u} e^{-i\omega t} e^{ikx} e^{ily} \qquad v = \hat{v} e^{-i\omega t} e^{ikx} e^{ily} \qquad (4) $$

where $\hat{\eta}$, $\hat{u}$ and $\hat{v}$ are the Fourier amplitudes. Substituting (4) into (1) - (3) gives

$$ -i\omega\hat{u} - f\hat{v} + ikg\hat{\eta} = 0, \qquad (5) $$

$$ -i\omega\hat{v} + f\hat{u} + ilg\hat{\eta} = 0, \qquad (6) $$

$$ -i\omega\hat{\eta} + ikH\hat{u} + ilH\hat{v} = 0. \qquad (7) $$
For non-zero ($\hat{\eta}$, $\hat{u}$, $\hat{v}$), (5) - (7) can be solved and give a dispersion polynomial of degree 3 in $\omega$. This polynomial has 3 solutions: the geostrophic mode $\omega = 0$ and two inertia-gravity modes:

$$ \omega_{AN}^2 = f^2 + gH(k^2 + l^2) = f^2 + gHK^2 \quad (8) $$

where $K^2 = k^2 + l^2$. The temporal frequency $\omega_{AN}$ can be normalized as $\sigma_{AN} = \frac{\omega_{AN}}{f}$, which gives

$$ \sigma_{AN}^2 = 1 + \lambda^2 K^2 \quad (9) $$

where $\lambda^2 = \frac{gH}{f^2}$ is the squared Rossby radius of deformation $\lambda$. Finally, introducing a length scale $\tilde{h}$ to facilitate comparison with numerical results, we have

$$ \sigma_{AN}^2 = 1 + \mu(\tilde{K}\tilde{h})^2, \quad (10) $$

with $\mu = \frac{\lambda^2}{\tilde{h}^2}$ and $(\tilde{K}\tilde{h})^2 = (\tilde{k}\tilde{h})^2 + (\tilde{l}\tilde{h})^2$ is the non-dimensional wavenumber $(\tilde{k}\tilde{h} \in [0, \pi]$ and $\tilde{l}\tilde{h} \in [0, \pi])$. Plots of $\sigma_{AN}(kh, lh)$ for $\lambda/\tilde{h} = 2.0$ and $\lambda/\tilde{h} = 0.1$ are found in Figures 1 and 2.

Since $\omega$ is purely real, all modes are neutrally stable and neither amplify nor decay. Additionally, the dispersion relationship is monotonic and non-zero for all $kh, lh > 0$, provided $g, f$ and $H$ are all non-zero. Other than a discussion of the use of tensor product structure to obtain the 2D dispersion relationship in terms of the 1D dispersion relationship, we do not consider the limiting cases of $f = 0$ (pure gravity waves) or $g = H = 0$ (pure inertia waves) in this paper. The case of pure gravity waves was considered for the $Q_n^-$ family in [24], and both pure inertia and pure gravity waves for the $Q_2^-$ family are discussed in [26].

<figure><img src="image_2.png" /><figcaption>Figure 1: Plots of the analytic dispersion relationship $\sigma_{AN}$ versus non-dimensional wavenumbers $\tilde{kh}$ and $\tilde{l}\tilde{h}$ for $\lambda/\tilde{h} = 2.0$ (left) and $\lambda/\tilde{h} = 0.1$ (right). Note that $\sigma_{AN}$ is monotonic and non-zero for all $\tilde{kh}, \tilde{l}\tilde{h} > 0$.</figcaption></figure>

### 3. General Compatible Galerkin Discretization

We will now discretize (1) - (3) using compatible Galerkin methods, and discuss how to obtain the discrete linear modes corresponding to a particular choice of discretization.
<figure><img src="image_3.png" /><figcaption>Figure 2: Plots of the analytic dispersion relationship $\sigma_{AN}$ versus non-dimensional wavenumbers $\tilde{kh}$ and $\tilde{lh}$ for $\lambda/\tilde{h} = 2.0$ (left) and $\lambda/\tilde{h} = 0.1$ (right). Note that $\sigma_{AN}$ is monotonic and non-zero for all $\tilde{kh}, \tilde{lh} > 0$.</figcaption></figure>

Here we consider only the case when time is continuous. An extension to incorporate time discretization is the subject of future work.

### 3.1. 2D Spaces

Consider a pair of 1D spaces $A \subset H^1$ and $B \subset L_2$ that are a partition of unity and form the 1D discrete deRham complex $A \xrightarrow{\partial/\partial x} B$. These are precisely the class of spaces studied in [13]. A 2D discrete deRham complex is formed by taking tensor products of the 1D spaces:

$$W_0 = A \otimes A \subset H^1 \qquad (11)$$

$$W_1 = (A \otimes B)\hat{\mathbf{i}} + (B \otimes A)\hat{\mathbf{j}} \subset H(div) \qquad (12)$$

$$W_2 = B \otimes B \subset L_2 \qquad (13)$$

where $\hat{\mathbf{i}}$ and $\hat{\mathbf{j}}$ are the unit vectors in the $x$ and $y$ directions on the reference quadrilateral. This complex is illustrated for $d = 1$ in Figure 3. More details about the tensor product construction can be found in [12, 14, 23, 29]. In this paper we consider the $Q_n^- \Lambda^k$ family from finite element exterior calculus with $A = P_n^C$ and $B = P_{n-1}^{DG}$ [4, 5] and the $MGD_n$ family (mimetic Galerkin differences) with $A = GD_n$ and $B = DGD_{n-1}$ [12, 13]. These are the extension of the element pairs in [13] to quadrilaterals. Recall that $P_n^C$ is the order $n$ Lagrange space, $P_{n-1}^{DG}$ the order $n-1$ discontinuous Lagrange space, $GD_n$ the order $n$ Galerkin difference space and $DGD_{n-1}$ the order $n-1$ discontinuous Galerkin difference space. Alternative choices for $A$ and $B$ lead to the mimetic spectral element method [14] and isogeometric discrete differential forms [7, 8, 14]. As shown in Section 4, the tensor product construction can be exploited to simplify the dispersion analysis, especially as compared to triangles or other polygons.
<figure><img src="image_4.png" /><figcaption>Figure 3: The 2D discrete deRham complex formed by taking tensor products of A and B. The operators $\nabla^T$ and $\nabla$ are strong, while $\tilde{\nabla}$ and $\tilde{\nabla}^T$ are weak (defined by integration by parts). The degrees of freedom are illustrated for $d = 1$.</figcaption></figure>

### 3.2. Discrete Model

Choose $\eta \in W_2$ and $\mathbf{v} = (u, v) \in W_1$. This means

$$ u \in A \otimes B \qquad v \in B \otimes A \qquad \eta \in B \otimes B \quad (14) $$

which corresponds with the placement in 1D of $\eta, v \in B$ and $u \in A$ [13]. Given test functions $\eta^* \in W_2$ and $\mathbf{v}^* = (u^*, v^*) \in W_1$, (1) - (3) are discretized as (after appropriate integration by parts)

$$ \left\langle u^*, \frac{\partial u}{\partial t} \right\rangle - f \left\langle u^*, v \right\rangle - g \left\langle \frac{\partial}{\partial x} u^*, \eta \right\rangle = 0, \quad (15) $$

$$ \left\langle v^*, \frac{\partial v}{\partial t} \right\rangle + f \left\langle v^*, u \right\rangle - g \left\langle \frac{\partial}{\partial y} v^*, \eta \right\rangle = 0, \quad (16) $$

$$ \left\langle \eta^*, \frac{\partial \eta}{\partial t} \right\rangle + H \left\langle \eta^*, \frac{\partial}{\partial x} u \right\rangle + H \left\langle \eta^*, \frac{\partial}{\partial y} v \right\rangle = 0, \quad (17) $$

where $\langle a, b \rangle$ is the $L_2$ inner product. Written in matrix form, (15) - (17) are

$$ \mathbf{M}_u \frac{\partial \mathbf{u}}{\partial t} - f \mathbf{C}_u \mathbf{v} + g \mathbf{G}_u \boldsymbol{\eta} = 0, \quad (18) $$

$$ \mathbf{M}_v \frac{\partial \mathbf{v}}{\partial t} + f \mathbf{C}_v \mathbf{u} + g \mathbf{G}_v \boldsymbol{\eta} = 0, \quad (19) $$

$$ \mathbf{M}_{\eta} \frac{\partial \boldsymbol{\eta}}{\partial t} + H \mathbf{D}_u \mathbf{u} + H \mathbf{D}_v \mathbf{v} = 0, \quad (20) $$

where $\mathbf{u}, \mathbf{v}$ and $\boldsymbol{\eta}$ denote the basis coefficients in the Galerkin expansion and the matrices are given by

$$ \begin{gathered} \mathbf{M}_u \frac{\partial \mathbf{u}}{\partial t} = \left\langle u^*, \frac{\partial u}{\partial t} \right\rangle, \quad \mathbf{M}_v \frac{\partial \mathbf{v}}{\partial t} = \left\langle v^*, \frac{\partial v}{\partial t} \right\rangle, \quad \mathbf{M}_{\eta} \frac{\partial \boldsymbol{\eta}}{\partial t} = \left\langle \eta^*, \frac{\partial \eta}{\partial t} \right\rangle, \\ \mathbf{G}_u \boldsymbol{\eta} = - \left\langle \frac{\partial}{\partial x} u^*, \eta \right\rangle, \quad \mathbf{G}_v \boldsymbol{\eta} = - \left\langle \frac{\partial}{\partial y} v^*, \eta \right\rangle, \quad \mathbf{D}_u \mathbf{u} = \left\langle \eta^*, \frac{\partial}{\partial x} u \right\rangle, \quad \mathbf{D}_v \mathbf{v} = \left\langle \eta^*, \frac{\partial}{\partial y} v \right\rangle, \\ \mathbf{C}_u \mathbf{v} = \langle u^*, v \rangle, \quad \mathbf{C}_v \mathbf{u} = \langle v^*, u \rangle. \end{gathered} $$
Use of compatible Galerkin spaces ensures

$$ M_u = M_u^*, \quad M_v = M_v^*, \quad M_\eta = M_\eta^*, \qquad (21) $$

$$ G_u^* = -D_u, \quad G_v^* = -D_v, \qquad (22) $$

$$ C_u^* = C_v. \qquad (23) $$

with $A^*$ the Hermitian transpose. These properties are inherited by the matrices $\tilde{A}$ and $\tilde{A}$ that arise after performing Fourier analysis. The properties (21) - (23) will be useful in proving some features of the dispersion relationship.

### 3.3. Fourier Analysis

Now consider a doubly periodic computational grid made of $N$ by $N$ uniform elements of width $h$. For what follows, $d$ is an integer that depends on the choice of spaces: $d = n$ for the $Q_n^- \Lambda^k$ family, and $d = 1$ for the $MGD_n$ family. In fact, $d$ is the number of distinct Fourier coefficients in 1D for each variable, and therefore the number of distinct Fourier coefficients in 2D are $d^2$. Thus, the matrices in (18) - (20) have size $(Nd)^2 \times (Nd)^2$. Assume a Fourier dependence in space and time for the basis coefficients:

$$ u_j(t) = \tilde{u}_j e^{-i\omega t} e^{ikx_j} e^{ily_j}, \quad v_j(t) = \tilde{v}_j e^{-i\omega t} e^{ikx_j} e^{ily_j}, \quad \eta_j(t) = \tilde{\eta}_j e^{-i\omega t} e^{ikx_j} e^{ily_j}, \qquad (24) $$

where $(x_j, y_j)$ is the nodal location of the $j$th basis function for the relevant variable and $(\tilde{u}_j, \tilde{v}_j, \tilde{\eta}_j)$ are the Fourier amplitudes, with $j = 1, \dots, (Nd)^2$. Inserting (24) into (18) - (20) yields

$$ -i\omega \tilde{M}_u \tilde{u} - f \tilde{C}_u \tilde{v} + g \tilde{G}_u \tilde{\eta} = 0, \qquad (25) $$

$$ -i\omega \tilde{M}_v \tilde{v} + f \tilde{C}_v \tilde{u} + g \tilde{G}_v \tilde{\eta} = 0, \qquad (26) $$

$$ -i\omega \tilde{M}_\eta \tilde{\eta} + H \tilde{D}_u \tilde{u} + H \tilde{D}_v \tilde{v} = 0, \qquad (27) $$

where

$$ \tilde{u} = [\tilde{u}_1, \dots, \tilde{u}_{(Nd)^2}], \quad \tilde{v} = [\tilde{v}_1, \dots, \tilde{v}_{(Nd)^2}], \quad \tilde{\eta} = [\tilde{\eta}_1, \dots, \tilde{\eta}_{(Nd)^2}]. \qquad (28) $$

The relationship between $A$ and $\tilde{A}$ is given by

$$ \tilde{A}_{mn} = A_{mn} e^{ik(x_m - x_0)} e^{il(y_m - y_0)}, \qquad (29) $$

where we have normalized each row by $e^{ikx_0} e^{ily_0}$, with $(x_0, y_0)$ the lower left corner of a specific element: for example, the middle element located at $(\frac{N}{2}, \frac{N}{2})$. Using the translational invariance of the basis functions, there are actually only $d^2$ unique degrees of freedom for $\tilde{u}$, $\tilde{v}$ and $\tilde{\eta}$:

$$ \hat{u} = [\hat{u}_1, \dots, \hat{u}_{d^2}]^T, \quad \hat{v} = [\hat{v}_1, \dots, \hat{v}_{d^2}]^T, \quad \hat{\eta} = [\hat{\eta}_1, \dots, \hat{\eta}_{d^2}]^T. \qquad (30) $$

Thus we can write (28) as

$$ \tilde{u} = [\hat{u}_1, \dots, \hat{u}_{d^2}, \hat{u}_1, \dots, \hat{u}_{d^2}, \dots], \qquad (31) $$
and similarly for $\tilde{\mathbf{v}}$ and $\tilde{\eta}$. Considering now only the $d^2$ rows for the degrees of freedom associated with a specific element (the one used for normalization), the final discrete system is

$$-i\omega \hat{\mathbf{M}}_u \hat{\mathbf{u}} - f \hat{\mathbf{C}}_u \hat{\mathbf{v}} + g \hat{\mathbf{G}}_u \hat{\eta} = 0, \quad (32)$$

$$-i\omega \hat{\mathbf{M}}_v \hat{\mathbf{v}} + f \hat{\mathbf{C}}_v \hat{\mathbf{u}} + g \hat{\mathbf{G}}_v \hat{\eta} = 0, \quad (33)$$

$$-i\omega \hat{\mathbf{M}}_{\eta} \hat{\eta} + H \hat{\mathbf{D}}_u \hat{\mathbf{u}} + H \hat{\mathbf{D}}_v \hat{\mathbf{v}} = 0, \quad (34)$$

The matrices $\tilde{\mathbf{A}}$ and $\hat{\mathbf{A}}$ are related by

$$\hat{\mathbf{A}}_{mn} = \sum_{p} \tilde{\mathbf{A}}_{mp}, \quad (35)$$

with $m = 1, \dots, d^2$, $n = 1, \dots, d^2$, and $p$ is the set of integers belonging to the interval $[1, \dots, (Nd)^2]$ such that $\tilde{x}_p = \hat{x}_n$. The final matrices $\hat{\mathbf{A}}$ will have the same properties (21) - (23) as the original matrices $\mathbf{A}$, which is essential to proving various features of the dispersion relationship. Equations (32) - (34) can be written in the form of a generalized eigenvalue problem (GEP)

$$-i\omega \hat{\mathbf{M}} \hat{\mathbf{x}} = \hat{\mathbf{S}} \hat{\mathbf{x}}, \quad (36)$$

where $\hat{\mathbf{x}} = [\hat{u}_1, \dots, \hat{u}_{d^2}, \hat{v}_1, \dots, \hat{v}_{d^2}, \hat{\eta}_1, \dots, \hat{\eta}_{d^2}]^T$ and

$$\hat{\mathbf{M}} = \begin{pmatrix} \hat{\mathbf{M}}_u & 0 & 0 \\ 0 & \hat{\mathbf{M}}_v & 0 \\ 0 & 0 & \hat{\mathbf{M}}_\eta \end{pmatrix}, \quad \hat{\mathbf{S}} = \begin{pmatrix} 0 & f \hat{\mathbf{C}}_u & -g \hat{\mathbf{G}}_u \\ -f \hat{\mathbf{C}}_v & 0 & -g \hat{\mathbf{G}}_v \\ -H \hat{\mathbf{D}}_u & -H \hat{\mathbf{D}}_v & 0 \end{pmatrix}. \quad (37)$$

Both $\hat{\mathbf{M}}$ and $\hat{\mathbf{S}}$ are $3 \times 3$ block matrices of total size $3d^2 \times 3d^2$, with individual blocks of size $d^2 \times d^2$. As a note, we recover the case of pure gravity waves when $f = 0$ and the case of pure inertia waves when $g = H = 0$. Equation (36) reduces to the eigenvalue problem (since $\hat{\mathbf{M}}$ is invertible)

$$-i\omega \mathbf{I} \hat{\mathbf{x}} = \hat{\mathbf{M}}^{-1} \hat{\mathbf{S}} \hat{\mathbf{x}} = \hat{\mathbf{F}} \hat{\mathbf{x}} \quad (38)$$

The eigenvalues of $\hat{\mathbf{F}} = \hat{\mathbf{M}}^{-1} \hat{\mathbf{S}}$ correspond to the linear modes of the discretized system, and form the dispersion relationship. Unfortunately, unlike the 1D case, it is not possible to obtain a reduced formulation for (36). This is because the reduced formulation relies on the existence of zero blocks in the generalized eigenvalue problem, and there are none for 2D inertia-gravity waves. A reduction is possible in 2D only in the special cases of pure gravity waves (due to the existence of zero blocks) or pure inertia waves (due to the matrix problem being only 2 by 2). See Section 4 for more details.

**Proposition 1.** The generalized eigenvalue problem (36) has $3d^2$ purely imaginary eigenvalues and thus $\omega$ is real.
Proof. Start with (36), and note that it is equivalent to

$$ \omega \hat{\mathbf{M}} \hat{\mathbf{x}} = i \hat{\mathbf{S}} \hat{\mathbf{x}}, \qquad (39) $$

Now replace $\hat{\eta}$ with $\sqrt{\frac{H}{g}}\hat{\eta}$ in order to symmetrize the terms in $\hat{\mathbf{S}}$ involving $g$ and $H$

$$ \hat{\mathbf{S}} = \begin{pmatrix} 0 & f\hat{\mathbf{C}}_u & -\sqrt{gH}\hat{\mathbf{G}}_u \\ -f\hat{\mathbf{C}}_v & 0 & -\sqrt{gH}\hat{\mathbf{G}}_v \\ -\sqrt{gH}\hat{\mathbf{D}}_u & -\sqrt{gH}\hat{\mathbf{D}}_v & 0 \end{pmatrix} \qquad (40) $$

where we have used the same symbols $\hat{\eta}$ and $\hat{\mathbf{S}}$ to ease notation. Using (21) - (23), we obtain

$$ \hat{\mathbf{M}} = \hat{\mathbf{M}}^*, \qquad (41) $$

and

$$ \hat{\mathbf{S}} = -\hat{\mathbf{S}}^*. \qquad (42) $$

and thus $(i\hat{\mathbf{S}})^* = i\hat{\mathbf{S}}$. Since both $\hat{\mathbf{M}}$ and $i\hat{\mathbf{S}}$ are Hermitian matrices, the revised generalized eigenvalue problem (39) has real eigenvalues, and thus $\omega$ is real. $\square$

**Conjecture 1.** The $3d^2$ eigenvalues are divided into $d^2$ geostrophic modes $\omega = 0$ and $d^2$ pairs of inertia-gravity waves: $d^2$ positive modes $\omega_{pos,i}$ and $d^2$ negative modes $\omega_{neg,i}$ with $\omega_{neg,i} = -\omega_{pos,i}$ (where $i = 1, \dots, d^2$).

This has been observed analytically for $d = 1, 2$, and numerically up to $d = 4$ for the case of 2D inertia-gravity waves, and analytically in the case of 1D IGW up to $d = 15$ [13]. However, we have been unable to find a proof of this conjecture using the matrix properties of $\hat{\mathbf{M}}$ and $\hat{\mathbf{S}}$, although presumably one must exist. A closely related statement has been proven in [9] using the Helmholtz decomposition for $W_1$ (see sections 2.5 and 2.8 in that paper) in the context of linear modes for a compatible Galerkin discretization on arbitrary meshes. In this case we deal with tensor-product elements on a doubly periodic domain, which ensures that there are no spurious inertia-gravity or geostrophic modes.

### 3.4. Non-Dimensional Spatial Wavenumbers

It is useful to define an adjusted grid spacing

$$ \tilde{h} = \frac{h}{d}, \qquad (43) $$

For equally-spaced degrees of freedom (nodes) this is the distance between nodes. However, in practice non-uniform nodes are often used (to improve the conditioning of the resulting matrices), and so $\tilde{h}$ can be viewed as an average distance. Using $\tilde{h}$, non-dimensional wavenumbers are naturally defined as

$$ \tilde{kh} = k\tilde{h}, \quad l\tilde{h} = l\tilde{h}, \qquad (44) $$
with the associated wavelengths

$$ \zeta_x = \frac{2\pi}{k} = \frac{2\pi\tilde{h}}{\tilde{k}\tilde{h}}, \quad \zeta_y = \frac{2\pi}{l} = \frac{2\pi\tilde{h}}{\tilde{l}\tilde{h}}, \qquad (45) $$

and non-dimensional wavelengths

$$ \tilde{\zeta}_x = \frac{\zeta_x}{\tilde{h}} = \frac{2\pi}{\tilde{k}\tilde{h}}, \quad \tilde{\zeta}_y = \frac{\zeta_y}{\tilde{h}} = \frac{2\pi}{\tilde{l}\tilde{h}}. \qquad (46) $$

For a grid with N by N elements, the maximal wavelengths are $\zeta_x = \zeta_y = Nd\tilde{h} = Nh$, with $k = l = \frac{2\pi}{Nh}$, $\tilde{k}\tilde{h} = \tilde{l}\tilde{h} = \frac{2\pi}{Nd}$ and $\tilde{\zeta}_x = \tilde{\zeta}_y = Nd$; and the minimum wavelengths are $\zeta_x = \zeta_y = 2\tilde{h}$, with $k = l = \frac{\pi}{\tilde{h}}$, $\tilde{k}\tilde{h} = \tilde{l}\tilde{h} = \pi$ and $\tilde{\zeta}_x = \tilde{\zeta}_y = 2$. These are useful in the definition of effective resolution, and also in determining the relevant wavenumber ranges to plot over.

### 3.5. Branches and Spectral Gaps

The eigenvalue problem in (36) has $d^2$ pairs of solutions corresponding to inertia-gravity modes. For a given non-dimensional wavenumber ($\tilde{k}\tilde{h}$, $\tilde{l}\tilde{h}$), when $d \ge 2$, only one of these solutions should be retained in order to correspond to the analytic dispersion relationship. The procedure for this is given in Section 3.7. Therefore, for each solution of (36) there exists a non-dimensional wavenumber range over which it is retained, and the particular solution of (36) over this range is termed a branch. The remaining parts of each solution are mathematical artifacts that should be disregarded. The union of all branches constitutes the discrete dispersion relationship. This corresponds exactly to the situation in 1D, as discussed in [13]. At the boundaries between two branches, the dispersion relationship is usually discontinuous. This defines a spectral gap as a spatial wavenumber for which the dispersion relationship becomes double valued, and these gaps result in a piecewise continuous dispersion relationship. In 1D, these gaps are point discontinuities. As shown in Section 5, on quadrilaterals in 2D these gaps become line discontinuities.

### 3.6. Solution Procedure

Following [13], two complementary approaches are used to solve (36). In the first, we use a computer algebra system (Maple) to construct the matrices (see Section 4 for how the tensor product structure is exploited to facilitate this) and solve the determinant to obtain $\sigma(kh, lh, \mu)$. The asymptotics are then obtained by performing a series expansion around $h = 0$. This procedure is only tractable at relatively moderate values of $d$ (in fact, only $d = 1, 2$ are solvable), due to the need to compute the determinant of a $3d^2 \times 3d^2$ matrix and solve the resulting polynomial of degree $3d^2$. Some examples of the Maple sheets (run using Maple 2018) employed for the $Q_n^-\Lambda^k$ and $MGD_n$ families can be found in the Supplementary Materials available with the online version of this paper.

Therefore, a second approach is used to make the plots of dispersion relationships and effective resolution. Starting from a computational grid of size N by N with $N = \frac{60}{d}$, the
eigenvalues and eigenvectors (36) are determined numerically for each $(kh, lh)$ pair and $\mu$ of interest. The set of $(kh, lh)$ sampled is a uniformly spaced set belonging to $[\frac{2\pi}{Nd}, \pi]^2$. When $d \ge 2$ (which occurs for the $Q_n^- \Lambda^k$ family with $n \ge 2$), there are $d^2$ eigenvalues. Then, using the procedure in Section 3.7, the solution is determined from these $d^2$ eigenvalues.

### 3.7. Branch determination for $d \ge 2$

The following procedure is used to select the appropriate branch (it is the 2D generalization of the approach in [13]):

1. Let $(kh, lh) \in (0, \pi]^2$ and numerically compute the $d^2$ eigenvalues $\omega_m$ ($m = 1, \dots, d^2$) and eigenvectors $\hat{x}_m$.

2. For each eigenvalue and eigenvector pair, determine the branch $(p_x, p_y)$ that it belongs to, with $p_x = 1, \dots, d$ and $p_y = 1, \dots, d$. This is done by constructing the discrete solution corresponding to $\hat{x}_m$ and comparing it to the continuous solution $e^{ikx}e^{ily}$.

3. Determine the effective wavenumber $(\tilde{k}h, \tilde{l}h)$ associated with $\omega_{p_x, p_y}$ using the explicit formula

$$ \tilde{k}h = \frac{(-1)^{\alpha_x}(kh) + 2\pi\beta_x}{d}, \quad \tilde{l}h = \frac{(-1)^{\alpha_y}(lh) + 2\pi\beta_y}{d}, \qquad (47) $$

where $\alpha_x = 1 + p_x$, $\alpha_y = 1 + p_y$, $\beta_x = [\frac{p_x}{2}]$ and $\beta_y = [\frac{p_y}{2}]$. This gives an effective wavenumber range of $(0, \pi]^2$.

## 4. Exploiting Tensor Product Structure

In [24] it is claimed without proof (see equations (80) - (85) in that paper) that the 2D IW, GW and IGW dispersion relationships for the $Q_n^- \Lambda^k$ family (and other finite element type compatible Galerkin schemes with multiple degrees of freedom when $n \ge 2$, ie $d = n$) can be constructed from the 1D GW and IW relationships by exploiting the tensor product structure. Here we will prove that this claim is true for 2D IW and 2D GW for any $d$ and for 2D IGW when $d = 1$. The case $d = 1$ is special since there is only only degree of freedom per geometric entity, and the resulting determinant matrix is only 3 x 3. However, for 2D IGW with $d \ge 2$ we provide counterexamples for the $Q_n^- \Lambda^k$ family in the form of numerical dispersion relationship calculations and asymptotic convergence results. Nevertheless, the tensor product structure can still be used in the 2D IGW case to determine the matrices $\tilde{M}$ and $\tilde{S}$ in terms of the 1D matrices.

### 4.1. Kronecker Products

Before continuing, it is useful to recall some basic properties of the Kronecker (or tensor) product [31]:

* Hermitian Transpose:

$$ (A \otimes B)^* = A^* \otimes B^* \qquad (48) $$
* Determinant: If A is an n x n matrix and B is an m x m matrix then

$$|A \otimes B| = |A|^m |B|^n \qquad (49)$$

* Inverse: If A and B are invertible then

$$(A \otimes B)^{-1} = A^{-1} \otimes B^{-1} \qquad (50)$$

* Mixed Product: If A, B, C, D are sized such that the products AC and BD can be formed then

$$(A \otimes B)(C \otimes D) = A C \otimes B D \qquad (51)$$

Also recall the Kronecker sum ⊕ for an n x n matrix A and m x m matrix B

$$A \oplus B = A \otimes I_m + I_n \otimes B \qquad (52)$$

where $I_k$ is the k x k identity matrix.

### 4.1.1. Eigenvalues of Kronecker Products and Sums

Of particular interest are the eigenvalues of Kronecker products and sums. Let $\alpha_i$ denote the eigenvalues of A and $\beta_j$ the eigenvalues of B. Then the eigenvalues of the Kronecker product

$$C = A \otimes B \qquad (53)$$

are

$$c_{ij} = \alpha_i \beta_j \qquad (54)$$

or in other words, simply the pairwise product of the eigenvalues of A and B. Similarly, the eigenvalues of the Kronecker sum

$$D = A \oplus B = A \otimes I + I \otimes B \qquad (55)$$

are

$$\delta_{ij} = \alpha_i + \beta_j \qquad (56)$$

or in other words, simply the pairwise sum of the eigenvalues of A and B. These two key facts will be used to show that the 2D IW and GW dispersion relationships can be expressed in terms of the 1D IW and GW dispersion relationships. This, however, is not true for 2D IGW.

### 4.2. 1D IGW Problem

Start by recalling the 1D IGW problem from [13].

$$\hat{M} = \begin{pmatrix} \hat{M}_C^x & 0 & 0 \\ 0 & \hat{M}_D^x & 0 \\ 0 & 0 & \hat{M}_D^x \end{pmatrix}, \quad \hat{S} = \begin{pmatrix} 0 & f \hat{C}_u^x & -g \hat{G}^x \\ -f \hat{C}_v^x & 0 & 0 \\ -H \hat{D}^x & 0 & 0 \end{pmatrix} \qquad (57)$$
where we have added the superscript x to indicate that the spatial wavenumber is k. The submatrices $\hat{M}_C^x$, $\hat{M}_D^x$, $\hat{G}^x$, $\hat{D}^x$, $\hat{C}_u^x$ and $\hat{C}_v^x$ are defined in [13], and they satisfy

$$ \hat{M}_C^* = \hat{M}_C^x, \quad \hat{M}_D^* = \hat{M}_D^x, \qquad (58) $$

$$ \hat{G}^* = -\hat{D}^x, \quad \hat{C}_u^* = \hat{C}_v^x. \qquad (59) $$

When $g = H = 0$, the 1D IW problem is recovered; and when $f = 0$ the 1D GW problem is recovered. There are also similar matrices denoted with a superscript y, where the only change is that spatial wavenumber k has been replaced with l.

### 4.3. 2D Matrices

The tensor product structure of the basis functions means that the submatrices in 2D are Kronecker products of the submatrices in 1D. In fact, each submatrix in 2D is the product of matrix with superscript x (containing only terms with k) and matrix with superscript y (containing only terms with l):

$$ \hat{M}_u = \hat{M}_C^x \otimes \hat{M}_D^y, \quad \hat{M}_v = \hat{M}_D^x \otimes \hat{M}_C^y, \quad \hat{M}_\eta = \hat{M}_D^x \otimes \hat{M}_D^y, \qquad (60) $$

$$ \hat{D}_u = \hat{D}^x \otimes \hat{M}_D^y, \quad \hat{D}_v = \hat{M}_D^x \otimes \hat{D}^y, \qquad (61) $$

$$ \hat{G}_u = \hat{G}^x \otimes \hat{M}_D^y, \quad \hat{G}_v = \hat{M}_D^x \otimes \hat{G}^y, \qquad (62) $$

$$ \hat{C}_u = \hat{C}_u^x \otimes \hat{C}_v^y, \quad \hat{C}_v = \hat{C}_v^x \otimes \hat{C}_u^y, \qquad (63) $$

and they satisfy

$$ \hat{M}_u^* = \hat{M}_u, \quad \hat{M}_v^* = \hat{M}_v, \quad \hat{M}_\eta^* = \hat{M}_\eta, \qquad (64) $$

$$ \hat{G}_u^* = -\hat{D}_u, \quad \hat{G}_v^* = -\hat{D}_v, \quad \hat{C}_u^* = \hat{C}_v. \qquad (65) $$

This tensor product structure can be exploited in the solution of the 2D IW and GW problems, as detailed below. It is also useful in the construction of the matrices used to make dispersion plots and compute asymptotics for 2D IGW.

### 4.4. 2D Inertia Waves

#### 4.4.1. 1D IW

Letting $g = H = 0$ in the 1D IGW (57) gives the simplified form for 1D IW:

$$ \hat{M} = \begin{pmatrix} \hat{M}_C^x & 0 \\ 0 & \hat{M}_D^x \end{pmatrix}, \quad \hat{S} = \begin{pmatrix} 0 & f \hat{C}_u^x \\ -f \hat{C}_v^x & 0 \end{pmatrix}. \qquad (66) $$

plus a identical system using y. Inertia waves do not have any dependence on $\eta$, and therefore $\hat{x} = [\hat{u}_1, \dots, \hat{u}_d, \hat{v}_1, \dots, \hat{v}_d]^T$.
**Proposition 2.** The eigenvalues of (66) are the solution of the reduced eigenvalue problem

$$ \omega_{1D,x}^2 \hat{u} = f^2 (\hat{\mathbf{M}}_C^x)^{-1} \hat{\mathbf{C}}_u^x (\hat{\mathbf{M}}_D^x)^{-1} \hat{\mathbf{C}}_u^x \hat{u} = f^2 \hat{\mathbf{F}}_u^{\text{IW},x} \hat{u} \quad (67) $$

and

$$ \omega_{1D,y}^2 \hat{v} = f^2 (\hat{\mathbf{M}}_D^y)^{-1} \hat{\mathbf{C}}_v^y (\hat{\mathbf{M}}_C^y)^{-1} \hat{\mathbf{C}}_u^y \hat{v} = f^2 \hat{\mathbf{F}}_v^{\text{IW},y} \hat{v} \quad (68) $$

*Proof.* This is obtained by noting that (66) is solving for the determinant of

$$ \det(\hat{\mathbf{F}} + i\omega \mathbf{I}) = \det((\hat{\mathbf{M}})^{-1}\hat{\mathbf{S}} + i\omega \mathbf{I}) \quad (69) $$

Now use

$$ \hat{\mathbf{F}} + i\omega \mathbf{I} = \begin{bmatrix} i\omega \mathbf{I} & f(\hat{\mathbf{M}}_C^x)^{-1} \hat{\mathbf{C}}_u^x \\ -f(\hat{\mathbf{M}}_D^x)^{-1} \hat{\mathbf{C}}_v^x & i\omega \mathbf{I} \end{bmatrix} \quad (70) $$

and apply the 2 x 2 block matrix determinant formula that comes from the Schur complement, using either the upper left or lower right blocks since they are both trivially invertible provided $\omega \neq 0$, which we know is true due to the properties of the matrices. Using the lower right block gives (67), while using the upper left block gives (68). □

Here we have chosen to solve for $\hat{u}$ in $x$ and $\hat{v}$ in $y$ for the 1D problem. This corresponds with the choice to solve for $\hat{u}$ in 2D. The same exact expression for the 2D eigenvalues in terms of the 1D eigenvalues is obtained if we solve instead for $\hat{v}$ in $x$ and $\hat{u}$ in $y$ for the 1D problem and $\hat{v}$ in 2D.

### 4.4.2. 2D IW

Letting $g = H = 0$ in the 2D IGW problem (36) gives

$$ \hat{\mathbf{M}} = \begin{pmatrix} \hat{\mathbf{M}}_u & 0 \\ 0 & \hat{\mathbf{M}}_v \end{pmatrix}, \quad \hat{\mathbf{S}} = \begin{pmatrix} 0 & f \hat{\mathbf{C}}_u \\ -f \hat{\mathbf{C}}_v & 0 \end{pmatrix}. \quad (71) $$

This has the same form as (66).

**Proposition 3.** The solutions to 2D eigenvalue problem (71) can be written as

$$ \omega_{2D}^2 \hat{u} = f^2 (\hat{\mathbf{M}}_u)^{-1} \hat{\mathbf{C}}_u (\hat{\mathbf{M}}_v)^{-1} \hat{\mathbf{C}}_v \hat{u} = f^2 \hat{\mathbf{F}}_u^{\text{IW}} \hat{u} \quad (72) $$

*Proof.* This again follows by expressing $\mathbf{F} + i\omega \mathbf{I}$ as a block matrix

$$ \mathbf{F} + i\omega \mathbf{I} = \begin{bmatrix} i\omega \mathbf{I} & f(\mathbf{M}_u)^{-1} \mathbf{C}_u \\ -f(\mathbf{M}_v)^{-1} \mathbf{C}_v & i\omega \mathbf{I} \end{bmatrix} \quad (73) $$

and using the lower right block for the 2 x 2 block matrix determinant. □
**Proposition 4.** The matrix $\widehat{\mathbf{F}}_{u}^{\text{IW}}$ can be written as the Kronecker product of the 1D matrices $\widehat{\mathbf{F}}_{u}^{\text{IW},x}$ and $\widehat{\mathbf{F}}_{v}^{\text{IW},y}$

$$ \widehat{\mathbf{F}}_{u}^{\text{IW}} = \widehat{\mathbf{F}}_{u}^{\text{IW},x} \otimes \widehat{\mathbf{F}}_{v}^{\text{IW},y} \qquad (74) $$

Therefore, the eigenvalues of the 2D inertia wave problem (72) can be written as

$$ w_{2D}^2 = \frac{w_{1D,x}^2 w_{1D,y}^2}{f^2} \qquad (75) $$

*Proof.* Start by expanding the right hand side $\widehat{\mathbf{F}}_{u}^{\text{IW}}$ of (72) in terms of the 1D matrices using the tensor product structure (60) - (63) to get

$$ \widehat{\mathbf{F}}_{u}^{\text{IW}} = (\widehat{\mathbf{M}}_{u})^{-1} \widehat{\mathbf{C}}_{u} (\widehat{\mathbf{M}}_{v})^{-1} \widehat{\mathbf{C}}_{v} = (\widehat{\mathbf{M}}_{C}^{x} \otimes \widehat{\mathbf{M}}_{D}^{y})^{-1} (\widehat{\mathbf{C}}_{u}^{x} \otimes \widehat{\mathbf{C}}_{v}^{y}) (\widehat{\mathbf{M}}_{D}^{x} \otimes \widehat{\mathbf{M}}_{C}^{y})^{-1} (\widehat{\mathbf{C}}_{v}^{x} \otimes \widehat{\mathbf{C}}_{u}^{y}) \quad (76) $$

Applying the inverse property (50) this becomes

$$ \widehat{\mathbf{F}}_{u}^{\text{IW}} = ((\widehat{\mathbf{M}}_{C}^{x})^{-1} \otimes (\widehat{\mathbf{M}}_{D}^{y})^{-1})(\widehat{\mathbf{C}}_{u}^{x} \otimes \widehat{\mathbf{C}}_{v}^{y})((\widehat{\mathbf{M}}_{D}^{x})^{-1} \otimes (\widehat{\mathbf{M}}_{C}^{y})^{-1})(\widehat{\mathbf{C}}_{v}^{x} \otimes \widehat{\mathbf{C}}_{u}^{y}) \quad (77) $$

Now using the mixed property (51) this is

$$ \widehat{\mathbf{F}}_{u}^{\text{IW}} = ((\widehat{\mathbf{M}}_{C}^{x})^{-1} \widehat{\mathbf{C}}_{u}^{x} (\widehat{\mathbf{M}}_{D}^{x})^{-1} \widehat{\mathbf{C}}_{v}^{x}) \otimes ((\widehat{\mathbf{M}}_{D}^{y})^{-1} \widehat{\mathbf{C}}_{v}^{y} (\widehat{\mathbf{M}}_{C}^{y})^{-1} \widehat{\mathbf{C}}_{u}^{y}) = \widehat{\mathbf{F}}_{u}^{\text{IW},x} \otimes \widehat{\mathbf{F}}_{v}^{\text{IW},y} \quad (78) $$

Therefore, using (54) we can express the eigenvalues of $\widehat{\mathbf{F}}_{u}^{\text{IW}}$ in terms of those for $\widehat{\mathbf{F}}_{u}^{\text{IW},x}$ and $\widehat{\mathbf{F}}_{v}^{\text{IW},y}$

$$ \frac{w_{2D}^2}{f^2} = \frac{w_{1D,x}^2 w_{1D,y}^2}{f^4} \qquad (79) $$

which can be simplified to give (75). □

Note that (75) is the same as equation (80) (with $g = H = 0$) from [24].

## 4.5. 2D Gravity Waves

### 4.5.1. 1D GW

In 1D, gravity waves have no dependence on $v$ and therefore $\hat{\mathbf{x}} = [\hat{u}_1, \dots, \hat{u}_d, \hat{\eta}_1, \dots, \hat{\eta}_d]^T$. Letting $f = 0$ in the 1D IGW problem (57) gives the simplified form for 1D GW:

$$ \widehat{\mathbf{M}} = \begin{pmatrix} \widehat{\mathbf{M}}_{C}^{x} & 0 \\ 0 & \widehat{\mathbf{M}}_{D}^{x} \end{pmatrix}, \quad \widehat{\mathbf{S}} = \begin{pmatrix} 0 & -g\widehat{\mathbf{G}}^{x} \\ -H\widehat{\mathbf{D}}^{x} & 0 \end{pmatrix}. \qquad (80) $$

plus an identical system in $y$ (with $k$ replaced by $l$).

**Proposition 5.** The solutions to 1D eigenvalue problem (80) can be written as

$$ \omega_{1D,x}^2 \hat{\eta} = -gH(\widehat{\mathbf{M}}_{D}^{x})^{-1}\widehat{\mathbf{D}}^{x}(\widehat{\mathbf{M}}_{C}^{x})^{-1}\widehat{\mathbf{G}}^{x}\hat{v} = -gH\widehat{\mathbf{F}}_{\eta}^{\text{GW},x}\hat{\eta} \qquad (81) $$

or the equivalent in $y$.

*Proof.* This is computed similarly to Proposition 2, using the 2x 2 block matrix determinant formula, and the proof is left to the interested reader. □
### 4.5.2. 2D GW

Letting $f = 0$ in the 2D IGW problem (36) gives the 2D GW problem

$$ \widehat{\mathbf{M}} = \begin{pmatrix} \widehat{\mathbf{M}}_u & 0 & 0 \\ 0 & \widehat{\mathbf{M}}_v & 0 \\ 0 & 0 & \widehat{\mathbf{M}}_\eta \end{pmatrix}, \quad \mathbf{S} = \begin{pmatrix} 0 & 0 & -g\widehat{\mathbf{G}}_u \\ 0 & 0 & -g\widehat{\mathbf{G}}_v \\ -H\widehat{\mathbf{D}}_u & -H\widehat{\mathbf{D}}_v & 0 \end{pmatrix}. \qquad (82) $$

**Proposition 6.** The solutions to 2D eigenvalue problem (82) can be written as

$$ \omega_{2D}^2 \hat{\eta} = -gH \left[ (\widehat{\mathbf{M}}_\eta)^{-1} \widehat{\mathbf{D}}_u (\widehat{\mathbf{M}}_u)^{-1} \widehat{\mathbf{G}}_u + (\widehat{\mathbf{M}}_\eta)^{-1} \widehat{\mathbf{D}}_v (\widehat{\mathbf{M}}_v)^{-1} \widehat{\mathbf{G}}_v \right] \hat{\eta} = -gH \widehat{\mathbf{F}}_\eta^{\text{GW}} \hat{\eta} \qquad (83) $$

plus n geostrophic solutions $\omega = 0$.

*Proof.* Start by rewriting (82) as

$$ \begin{pmatrix} i\omega\widehat{\mathbf{M}} & -g\widehat{\mathbf{G}} \\ -H\widehat{\mathbf{D}} & i\omega\widehat{\mathbf{M}}_\eta \end{pmatrix} \qquad (84) $$

where

$$ \widehat{\mathbf{M}} = \begin{pmatrix} \widehat{\mathbf{M}}_u & 0 \\ 0 & \widehat{\mathbf{M}}_v \end{pmatrix} \quad \widehat{\mathbf{G}} = \begin{pmatrix} \widehat{\mathbf{G}}_u \\ \widehat{\mathbf{G}}_v \end{pmatrix} \quad \widehat{\mathbf{D}} = \begin{pmatrix} \widehat{\mathbf{D}}_u & \widehat{\mathbf{D}}_v \end{pmatrix}. \qquad (85) $$

Now using the 2 x 2 block matrix determinant formula the determinant of (84) is

$$ \det(i\omega\widehat{\mathbf{M}})\det(i\omega\widehat{\mathbf{M}}_\eta - \frac{gH}{i\omega}\widehat{\mathbf{D}}\widehat{\mathbf{M}}^{-1}\widehat{\mathbf{G}}). \qquad (86) $$

The first term in (86) gives n geostrophic modes $\omega = 0$ along with a factor of $(i\omega)^n$ that can be absorbed into the second term, yielding

$$ \det(\omega^2 \widehat{\mathbf{M}}_\eta + gH \widehat{\mathbf{D}} \widehat{\mathbf{M}}^{-1} \widehat{\mathbf{G}}), \qquad (87) $$

which is equivalent to

$$ \det(\omega^2 + gH\widehat{\mathbf{M}}_\eta^{-1}\widehat{\mathbf{D}}\widehat{\mathbf{M}}^{-1}\widehat{\mathbf{G}}). \qquad (88) $$

Now using the definitions for $\widehat{\mathbf{D}}$, $\widehat{\mathbf{G}}$ and $\widehat{\mathbf{M}}$ we have

$$ \widehat{\mathbf{M}}_\eta^{-1} \widehat{\mathbf{D}} \widehat{\mathbf{M}}^{-1} \widehat{\mathbf{G}} = \widehat{\mathbf{M}}_\eta^{-1} \begin{pmatrix} \widehat{\mathbf{D}}_u & \widehat{\mathbf{D}}_v \end{pmatrix} \begin{pmatrix} \widehat{\mathbf{M}}_u & 0 \\ 0 & \widehat{\mathbf{M}}_v \end{pmatrix}^{-1} \begin{pmatrix} \widehat{\mathbf{G}}_u \\ \widehat{\mathbf{G}}_v \end{pmatrix} \qquad (89) $$

Expansion using block matrix multiplication then gives (83). $\square$

**Proposition 7.** The 2D matrix $\widehat{\mathbf{F}}_\eta^{\text{GW}}$ can be written as the Kronecker sum of the 1D matrices $\widehat{\mathbf{F}}_\eta^{\text{GW},x}$ and $\widehat{\mathbf{F}}_\eta^{\text{GW},y}$

$$ \widehat{\mathbf{F}}_\eta^{\text{GW}} = \widehat{\mathbf{F}}_\eta^{\text{GW},x} \otimes \mathbf{I} + \mathbf{I} \otimes \widehat{\mathbf{F}}_\eta^{\text{GW},y} = \widehat{\mathbf{F}}_\eta^{\text{GW},x} \oplus \widehat{\mathbf{F}}_\eta^{\text{GW},y} \qquad (90) $$

Therefore, the eigenvalues of the 2D gravity wave problem (82) can be written as

$$ w_{2D}^2 = w_{1D,x}^2 + w_{1D,y}^2 \qquad (91) $$
Proof. Start by expanding the first term of (83) in terms of the 1D matrices using as (60) - (63)

$$ (\widehat{\mathbf{M}}_{\eta})^{-1} \widehat{\mathbf{D}}_{u} (\widehat{\mathbf{M}}_{u})^{-1} \widehat{\mathbf{G}}_{u} = (\widehat{\mathbf{M}}_{\mathrm{D}}^{\mathrm{x}} \otimes \widehat{\mathbf{M}}_{\mathrm{D}}^{\mathrm{y}})^{-1} (\widehat{\mathbf{D}}^{\mathrm{x}} \otimes \widehat{\mathbf{M}}_{\mathrm{D}}^{\mathrm{y}}) (\widehat{\mathbf{M}}_{\mathrm{C}}^{\mathrm{x}} \otimes \widehat{\mathbf{M}}_{\mathrm{D}}^{\mathrm{y}})^{-1} (\widehat{\mathbf{G}}^{\mathrm{x}} \otimes \widehat{\mathbf{M}}_{\mathrm{D}}^{\mathrm{y}}) \quad (92) $$

Use the inverse property (50) to get

$$ ((\widehat{\mathbf{M}}_{\mathrm{D}}^{\mathrm{x}})^{-1} \otimes (\widehat{\mathbf{M}}_{\mathrm{D}}^{\mathrm{y}})^{-1})(\widehat{\mathbf{D}}^{\mathrm{x}} \otimes \widehat{\mathbf{M}}_{\mathrm{D}}^{\mathrm{y}})((\widehat{\mathbf{M}}_{\mathrm{C}}^{\mathrm{x}})^{-1} \otimes (\widehat{\mathbf{M}}_{\mathrm{D}}^{\mathrm{y}})^{-1})(\widehat{\mathbf{G}}^{\mathrm{x}} \otimes \widehat{\mathbf{M}}_{\mathrm{D}}^{\mathrm{y}}) \quad (93) $$

Now the mixed product property (51) gives

$$ ((\widehat{\mathbf{M}}_{\mathrm{D}}^{\mathrm{x}})^{-1} \widehat{\mathbf{D}}^{\mathrm{x}} (\widehat{\mathbf{M}}_{\mathrm{C}}^{\mathrm{x}})^{-1} \widehat{\mathbf{G}}^{\mathrm{x}}) \otimes ((\widehat{\mathbf{M}}_{\mathrm{D}}^{\mathrm{y}})^{-1} \widehat{\mathbf{M}}_{\mathrm{D}}^{\mathrm{y}} (\widehat{\mathbf{M}}_{\mathrm{D}}^{\mathrm{y}})^{-1} \widehat{\mathbf{M}}_{\mathrm{D}}^{\mathrm{y}}) = \widehat{\mathbf{F}}_{\eta}^{\widehat{\mathrm{GW}},x} \otimes \mathbf{I} \quad (94) $$

Identical calculations on the $(\widehat{\mathbf{M}}_{\eta})^{-1} \widehat{\mathbf{D}}_{v} (\widehat{\mathbf{M}}_{v})^{-1} \widehat{\mathbf{G}}_{v}$ term yield

$$ (\widehat{\mathbf{M}}_{\eta})^{-1} \widehat{\mathbf{D}}_{v} (\widehat{\mathbf{M}}_{v})^{-1} \widehat{\mathbf{G}}_{v} = \mathbf{I} \otimes \widehat{\mathbf{F}}_{\eta}^{\widehat{\mathrm{GW}},y} \quad (95) $$

Therefore (90) holds. Since the 2D GW matrix $\widehat{\mathbf{F}}_{\eta}^{\widehat{\mathrm{GW}},x}$ is the Kronecker sum of the 1D GW matrices $\widehat{\mathbf{F}}_{\eta}^{\widehat{\mathrm{GW}},y}$ and $\widehat{\mathbf{F}}_{\eta}^{\widehat{\mathrm{GW}},y}$, using (56) the eigenvalues of (82) are simply the sum of the eigenvalues of the 1D problems (81) in $x$ and $y$, and (91) holds. $\square$

Note that (91) is the same as equation (80) (with $f = 0$) from [24].

## 4.6. 2D Inertia-Gravity Waves

In [24] it is claimed without proof that the 2D IGW dispersion relationship can be constructed from the 1D GW and IW relationships. We have demonstrated this to be true for the 2D IW (Section 4.4) and 2D GW (Section 4.5). For 2D IGW, a proof for $d = 1$ is found below in Section 4.6. However, we have been unable to prove this assertion for $d \ge 2$. The expression of the 2D GW in terms of the 1D GW, and of the 2D IW in terms of the 1D IW rested on the existence of a simplified form and corresponding reduced eigenvalue problem. This simplified form was possible to derive only because there were either zero blocks (GW), or the matrix was only 2 x 2 (IW). In 1D, the reduced forms for IW, GW and IGW also rely on the presence of zero blocks. However, 2D IGW has a 3 x 3 block matrix without any zeros blocks, and therefore it does not appear to be possible to express the determinant in a reduced form.

Two lines of evidence supporting the claim that for $d \ge 2$ the 2D IGW eigenvalues cannot be expressed in terms of the 1D IW and GW eigenvalues are presented in Section 5.3. The first is a demonstration that the asymptotics of the 2D IGW dispersion relationship do not agree with those calculated using equation (80) from [24]. The second is a numerical comparison of the dispersion relationship computed using the approach in this paper and using the approach in [24], where it is found that they differ by around 1%. Another possibility might be that the 2D IGW relationship can be expressed in terms of the 1D IGW relationship. This is demonstrated to be impossible for $d = 1$ below, and therefore it seems very unlikely to hold for $d \ge 2$.
**Proposition 8.** When $d = 1$, the 2D IGW eigenvalues can be expressed in terms of the 2D GW and IW eigenvalues (and therefore the 1D GW and IW eigenvalues) as

$$w_{IGW}^2 = w_{GW}^2 + w_{IW}^2 = \frac{w_{1D,IW,x}^2 w_{1D,IW,y}^2}{f^2} + w_{GW,1D,x}^2 + w_{GW,1D,y}^2 \quad (96)$$

*Proof.* When $d = 1$, both $\hat{M}$ and $\hat{S}$ are 3 x 3 matrices, and there exists a direct formula for the eigenvalues of the matrix. The eigenvalue problem (36) can be written as

$$\left( \begin{array}{ccc} i\omega \hat{M}_u & f \hat{C}_u & -g \hat{G}_u \\ -f \hat{C}_v & i\omega \hat{M}_v & -g \hat{G}_v \\ -H \hat{D}_u & -H \hat{D}_v & i\omega \hat{M}_\eta \end{array} \right), \quad (97)$$

which has determinant (after factoring out the $\omega = 0$ root)

$$\omega_{IGW}^2 = f^2 \frac{\hat{C}_u \hat{C}_v}{\hat{M}_u \hat{M}_v} + gH \left( \frac{\hat{G}_v \hat{D}_v}{\hat{M}_v \hat{M}_\eta} + \frac{\hat{G}_u \hat{D}_u}{\hat{M}_u \hat{M}_\eta} \right). \quad (98)$$

Use of (21) - (23) was made in calculating this. By repeating this procedure for $f = 0$ and $g = H = 0$, it can also be shown that

$$\omega_{GW}^2 = gH \left( \frac{\hat{G}_v \hat{D}_v}{\hat{M}_v \hat{M}_\eta} + \frac{\hat{G}_u \hat{D}_u}{\hat{M}_u \hat{M}_\eta} \right), \quad (99)$$

$$\omega_{IW}^2 = f^2 \frac{\hat{C}_u \hat{C}_v}{\hat{M}_u \hat{M}_v}. \quad (100)$$

This proves the first part of (96), that $w_{IGW}^2 = w_{GW}^2 + w_{IW}^2$. Now use the results of Sections 4.4 and 4.5, which give $w_{IW}^2 = \frac{w_{1D,IW,x}^2 w_{1D,IW,y}^2}{f^2}$ and $w_{GW}^2 = w_{GW,1D,x}^2 + w_{GW,1D,y}^2$. Therefore the second part of (96) holds. $\square$

Note that (96) is just equation (80) from [24] with $n = 1$.

Similar arguments can be used to show that the 1D dispersion relationship (when $d = 1$) can be written as

$$w_{1D,IGW}^2 = w_{1D,GW}^2 + w_{1D,IW}^2, \quad (101)$$

in both $x$ and $y$ directions. However, comparing (96) and (101), it is not possible to express $w_{IGW}^2$ purely in terms of $w_{1D,IGW}^2$, even with the ability to express it in terms of $w_{1D,IW}^2$ and $w_{1D,GW}^2$. Since this is not possible when $d = 1$, it seems very likely that it is also not possible for $d \ge 2$.
### 5. $Q_n^- \Lambda^k$ family

The $Q_n^- \Lambda^k$ family (also as finite element exterior calculus or compatible finite elements) is constructed by choosing $\mathbb{A} = P_n^C$ and $\mathbb{B} = P_{n-1}^{DG}$, and is defined for all $n \ge 1$. A diagram of the Fourier degrees of freedom for $n = 2$ is found in Figure 4. Recall that $W_0 \subset H^1$, $W_1 \subset H(\mathrm{div})$ and $W_2 \subset L_2$. A key point here is that $d = n$, that is, there are $3n^2$ degrees of freedom in the Fourier analysis. This gives rise to multiple branches in the dispersion relationship, and therefore associated line discontinuities.

<figure><img src="image_5.png" /><figcaption>Figure 4: The Fourier degrees of freedom for the $Q_n^- \Lambda^k$ family with $n = 2$. Since $d = 2$ there are $2^2 = 4$ degrees of freedom for each quantity. The $W_0$ degrees of freedom ($\psi$, which are not used) are shown on the left, the two components of $W_1$ ($u$ and $v$) in the middle and the $W_2$ ($\eta$) degrees of freedom on the right.</figcaption></figure>

#### 5.1. Asymptotics

On a uniform grid, the maximal degree of the polynomials that occur in the integrals is $2n$ in $x$ and $y$. Therefore, any quadrature rule of at least that degree will exactly integrate the equations. One example is tensor-product Gauss-Lobatto quadrature with $n + 1$ points in each direction, which exactly integrates polynomials of degree $2n + 1$ in $x$ and $y$. Exact quadrature gives the following results for the asymptotics:

**Theorem 1.** In the limit as mesh spacing $h \to 0$ we obtain the asymptotic results for $n \le 2$:

$$ \omega = \omega_{AN} \pm \frac{1}{2^{2n+1} \prod_{j=1}^{n} (4j^2 - 1)} \frac{(-f^2(k^{2n} + l^{2n}) + g H (k^{2n+2} + l^{2n+2}))}{\sqrt{f^2 + g H (k^2 + l^2)}} (h)^{2n} + O(h^{2n+2}). $$

It is conjectured that this result holds $\forall n \ge 1$.

These results were obtained using Maple and the analytic solution procedure described above, and are at the limits of the software. Solving $n = 3$ would require computing the symbolic determinant of a 27 x 27 matrix, and then finding the roots of the resulting degree 27 polynomial. This proved to be intractable. However, there are two lines of evidence supporting the conjecture. The first is that when $l = 0$, the result for the $P_n^C - P_{n-1}^{DG}$ pair from [13] is recovered, which was shown to hold until $n = 15$. The second is that the 2D GW dispersion relationship (obtained by setting $f = 0$) satisfies this until at least $n = 12$ (an attempt to prove it beyond $n = 12$ was made, but again it was beyond the limits of Maple). The $O(h^{2n})$ order of convergence is superconvergence specific to a regular mesh, and on arbitrary smooth meshes the expected order is $O(h^{n+1})$.
## 5.2. Dispersion Relationships

Plots of the dispersion relationship $\sigma = \frac{\omega}{f}$ versus $(kh, l\tilde{h})$ for $n = 1, \dots, 4$ are found in Figures 5 and 6 with $\lambda/\tilde{h} = 2.0$ (a well-resolved Rossby radius) and in Figures 7 and 8 with $\lambda/\tilde{h} = 0.1$ (a poorly resolved Rossby radius). The use of $\lambda/\tilde{h}$ ensures the same effective resolution of the Rossby radius for all $n$. The long wavelength part of the spectrum is well-represented in all cases, but when $\lambda/\tilde{h} = 0.1$ the short wavelength part of the spectrum is significant underestimated and has a group velocity with the wrong sign. This portion is progressively confined to higher wavenumbers as $n$ increases. When $n \ge 2$, spectral gaps where the dispersion relationship is double-valued (in the form of line discontinuities) are found at $kh = \frac{j\pi}{n}$ ($j = 1, \dots, n-1$) and $l\tilde{h} = \frac{m\pi}{n}$ ($m = 1, \dots, n-1$). As $n$ increases, the jump in frequency at the discontinuities for large $kh$ and $l\tilde{h}$ increase, but the size of the other gaps decreases. It seems likely that numerical dissipation will be required to control these gaps, especially in the short wavelength portion of the spectrum. The use of such dissipation can be undesirable, especially since it is introduced for numerical stability rather than for physically-motivated reasons. In addition to the gaps, there is also significant overestimation of the frequency in the short wavelength part of the spectrum, that gets worse as $n$ increases. This is a plausible explanation of the observation that the CFL condition gets increasingly strict as $n$ increases.

## 5.3. Comparison with the results obtained in [24]

In [24], it is claimed (see equation (80) in that paper) that

$$w_{IGW}^2 = w_{GW}^2 + w_{IW}^2 = \frac{w_{1D,IW,x}^2 w_{1D,IW,y}^2}{f^2} + w_{GW,1D,x}^2 + w_{GW,1D,y}^2 \quad (102)$$

holds for any $d$. In Section 4.6 this was proven to be true for 2D IGW with $d = 1$, and for any $d$ for 2D IW ($g = H = 0$, Section 4.4) and 2D GW ($f = 0$, Section 4.5). However, as discussed previously, here we will present counterexamples for the case of 2D IGW when $d \ge 2$. The first piece of evidence is found by computing the difference $(w_{tp} - w_{direct})/w_{direct} * 100$ between the numerical dispersion relationship $w_{direct}$ found by the procedure in Section 3.6, and $w_{tp}$ obtained by using (102) for $n = 1, \dots, 4$. This is shown in Figure 9, which gives the results for $\lambda/\tilde{h} = 2.0$ and $n = 1, \dots, 4$ using exact quadrature. Similar results (not shown) are found for $\lambda/\tilde{h} = 0.1$. As expected, when $n = 1$ there is no difference. However, when $n \ge 2$, there is a difference of around 1%. The second piece of evidence comes from computing the asymptotics for $w_{tp}$ in the case $n = 2$

$$\omega_{tp} = \omega_{AN} \pm \frac{1}{1440} \frac{gH(k^6 + l^6)}{\sqrt{f^2 + gH(k^2 + l^2)}} h^4 + O(h^6) \quad (103)$$

Such a result is different than found in Section 5.1. Therefore, based on the counterexamples of the numerical dispersion relationships for $n = 1, \dots, 4$ and the asymptotics for $n = 2$, we conjecture that (102) holds for 2D IGW only when $d = 1$.
### 5.4. Mass Lumping for $n = 2$

As shown in Section 5.2, when $n \ge 2$ the 2D IGW dispersion relationship has line discontinuities. This was also shown for the 2D GW dispersion relationship in [24, 26]. However, it is possible to eliminate these discontinuities for $n = 2$ through a careful mass lumping of the $M_C^x$ and $M_C^y$ matrices. This was demonstrated for the 1D IGW dispersion relationship in [26, 35] and for the 2D GW and IW dispersion relationships in [26]. The latter also included numerical evidence in the form of a simulation that the gap was closed for 2D IGW, but no plots of the dispersion relationship were provided.

This lumping is accomplished by adding a term $\mathbf{ML}_a$ to the $M_C^x$ and $M_C^y$ matrices that are used in constructing $\mathbf{M}_u$ and $\mathbf{M}_v$, where

$$ \mathbf{ML}_a = \begin{pmatrix} \alpha & \beta & -\alpha \\ \beta & \gamma & \beta \\ -\alpha & \beta & \alpha \end{pmatrix}, \qquad (104) $$

with $\beta = \gamma = 0$ and $\alpha = \frac{1}{60}$. The precise value of $\alpha$ will depend on how the basis functions are normalized, and this value corresponds to using a reference element of size $[0, 1]^2$ and normalizing the basis functions such that the basis coefficients correspond to integrals over the relevant geometric entities: sampling at points for $W_0$, integrating over lines for $W_1$ and integrating over cells for $W_2$. These lumping parameters are derived by enforcing that $M_C^x$ ($M_C^y$) remains positive-definite, and that the spectral gap at $\widetilde{kh} = \frac{\pi}{2}$ ($\widetilde{lh} = \frac{\pi}{2}$) for 1D IGW is eliminated. As shown in Figure 10, simply using this lumping for $M_C^x$ and $M_C^y$ when they are used to construct $\mathbf{M}_u$ and $\mathbf{M}_v$ is sufficient to ensure that the 2D IGW gaps are also eliminated. Using this mass lumping and exact quadrature, the asymptotics are

$$ \omega = \omega_{AN} \mp \frac{\alpha}{4}(k^2 + l^2)\sqrt{gH(k^2 + l^2) + f^2h^2 + O(h^4)}. $$

As in 1D, the convergence rate is reduced from $O(h^4)$ to $O(h^2)$, and when $l = 0$ the result for the lumped $P_2^C - P_1^{DG}$ pair from [13] is recovered.

## 6. MGD$_n$ family

The $MGD_n$ family is constructed by choosing $A = GD_n$ [6] and $B = DGD_{n-1}$, and is defined for $n$ odd ($n = 1, 3, 5, 7...$). Recall that $GD_n$ is the space of order n Galerkin differences in 1D, and that $DGD_{n-1}$ is the corresponding compatible $L_2$ space in 1D constructed following [14]. More details about this construction can be found in [12]. A diagram of the degrees of freedom for $n = 3$ is found in Figure 11, and the basis functions for $W_0$, $W_1$ and $W_2$ are found in Figure 12. A key point is that $d = 1$ independent of $n$, and therefore there are only 3 degrees of freedom in the Fourier analysis. This is fundamental to the ability of the $MGD_n$ family to avoid spectral gaps.
## 6.1. Exact Quadrature Asymptotics

Just as for the $Q_n^- \Lambda^k$ family, using a quadrature rule of at least degree $2n$ in each direction will ensure that all integrals are exactly integrated. Using exact quadrature yields the following results for the asymptotics:

**Theorem 2.** With exact quadrature, in the limit as mesh spacing $h \to 0$ we obtain the asymptotic results for $n \le 17$:

$$ \omega = \omega_{AN} \pm \frac{\alpha_r (-f^2(k^{2n} + l^{2n}) + g H (k^{2n+2} + l^{2n+2}))}{\beta_r \sqrt{f^2 + g H (k^2 + l^2)}} (h)^{2n} + O(h^{2n+2}). $$

where $\alpha_e$ and $\beta_e$ are given in Table 1. It is conjectured that this result holds $\forall n \ge 1$.

As before, we obtained these results using Maple and the analytic solution procedure described in Section 3.6. When $l = 0$, the results for the $GD_n - DGD_{n-1}$ pair from [13] are recovered, albeit with a significant improvement in $n$ from $n = 9$ to $n = 17$. Again, the $O(h^{2n})$ order of convergence is superconvergence specific to a regular mesh, and on arbitrary smooth meshes the expected order is $O(h^{n+1})$.

<table><caption>Table 1: The coefficients $\alpha_e$ and $\beta_e$ for the $MGD_n$ family using exact quadrature.</caption><thead><tr><th>n</th><th>α<sub>e</sub></th><th>β<sub>e</sub></th></tr></thead><tbody><tr><td>n = 1</td><td>1</td><td>2<sup>3</sup> × 3</td></tr><tr><td>n = 3</td><td>17</td><td>2<sup>6</sup> × 3<sup>3</sup> × 7</td></tr><tr><td>n = 5</td><td>13 × 317</td><td>2<sup>10</sup> × 3<sup>5</sup> × 5<sup>2</sup> × 11</td></tr><tr><td>n = 7</td><td>47 × 318749</td><td>2<sup>13</sup> × 3<sup>6</sup> × 5<sup>3</sup> × 7<sup>2</sup> × 11 × 13</td></tr><tr><td>n = 9</td><td>44734915633</td><td>2<sup>18</sup> × 3<sup>8</sup> × 5<sup>3</sup> × 7<sup>3</sup> × 13 × 17 × 19</td></tr><tr><td>n = 11</td><td>107 × 5693 × 716143 × 6211</td><td>2<sup>21</sup> × 3<sup>11</sup> × 5<sup>5</sup> × 7<sup>3</sup> × 11<sup>2</sup> × 17 × 19 × 23</td></tr><tr><td>n = 13</td><td>631 × 5519 × 145523547391</td><td>2<sup>25</sup> × 3<sup>13</sup> × 5<sup>5</sup> × 7<sup>2</sup> × 11<sup>2</sup> × 13<sup>2</sup> × 17 × 19 × 23</td></tr><tr><td>n = 15</td><td>11304410571523 × 88334167</td><td>2<sup>28</sup> × 3<sup>15</sup> × 5<sup>7</sup> × 7<sup>2</sup> × 11<sup>3</sup> × 13<sup>2</sup> × 19 × 23 × 29 × 31</td></tr><tr><td>n = 17</td><td>6703742237 × 155719 × 3541 × 2297</td><td>2<sup>34</sup> × 3<sup>15</sup> × 5<sup>8</sup> × 7<sup>5</sup> × 11<sup>3</sup> × 13<sup>2</sup> × 17<sup>2</sup> × 23 × 29 × 31</td></tr></tbody></table>

## 6.2. Reduced Quadrature Asymptotics

It is useful to also consider reduced quadrature to ameliorate the computational cost of the $MGD_n$ elements. This becomes particularly important as the dimension increases. As proven in [13], reduced quadrature does not effect any of the desirable properties of the scheme and does not change the matrix sparsity pattern. Using a tensor product, 2 point Gaussian quadrature for all integrals yields the following result for asymptotics:
**Theorem 3.** Using reduced quadrature, in the limit as mesh spacing $h \to 0$ we obtain the asymptotic results for $n \le 19$:

$$ \omega = \omega_{AN} \pm \frac{\alpha_r}{\beta_r} \frac{g H (k^{2n} + l^{2n})}{\sqrt{f^2 + g H (k^2 + l^2)}} (h)^{n+1} + O(h^{n+3}), \quad (105) $$

where $\alpha_r$ and $\beta_r$ are given in Table 2. It is conjectured that (105) holds for $n \ge 1$.

As in the case of exact quadrature, the results for the $GD_n - DGD_{n-1}$ pair from [13] are recovered when $l = 0$, with again a significant improvement in $n$ from $n = 9$ to $n = 19$. It is interesting to note that this reduction in order is simply removing the superconvergence on uniform meshes, and still recovers the expected order on arbitrary smooth meshes. Therefore it is not clear if there will be significant differences between exact and reduced quadrature on more realistic (deformed) meshes.

<table><caption>Table 2: The coefficients $\alpha_r$ and $\beta_r$ for the $MGD_n$ family using reduced quadrature.</caption><thead><tr><th>n</th><th>α<sub>r</sub></th><th>β<sub>r</sub></th></tr></thead><tbody><tr><td>n = 3</td><td>-1</td><td>2<sup>3</sup> × 3<sup>3</sup> × 5</td></tr><tr><td>n = 5</td><td>-5</td><td>2<sup>5</sup> × 3<sup>4</sup> × 7</td></tr><tr><td>n = 7</td><td>-17</td><td>2<sup>8</sup> × 3<sup>3</sup> × 5 × 7</td></tr><tr><td>n = 9</td><td>-133741</td><td>2<sup>10</sup> × 3<sup>9</sup> × 5 × 7 × 11</td></tr><tr><td>n = 11</td><td>-23 × 1069 × 4111</td><td>2<sup>14</sup> × 3<sup>10</sup> × 5<sup>2</sup> × 7 × 11 × 13</td></tr><tr><td>n = 13</td><td>-44346011</td><td>2<sup>16</sup> × 3<sup>10</sup> × 5<sup>3</sup> × 7 × 13</td></tr><tr><td>n = 15</td><td>-824695933877</td><td>2<sup>19</sup> × 3<sup>14</sup> × 5<sup>3</sup> × 7<sup>2</sup> × 13 × 17</td></tr><tr><td>n = 17</td><td>-45381741803663</td><td>2<sup>21</sup> × 3<sup>15</sup> × 5<sup>3</sup> × 7<sup>2</sup> × 13 × 17 × 19</td></tr><tr><td>n = 19</td><td>-3837131 × 37661489</td><td>2<sup>26</sup> × 3<sup>15</sup> × 5<sup>3</sup> × 7<sup>3</sup> × 13 × 19</td></tr></tbody></table>

### 6.3. Dispersion Relationships

Plots of the dispersion relationship $\sigma = \frac{\omega}{f}$ versus $(kh, l\tilde{h})$ for $n = 1, 3, 5, 7$, are found in Figures 13 and 14 with $\lambda/\tilde{h} = 2.0$ and in Figures 15 and 16 with $\lambda/\tilde{h} = 0.1$. The use of $\lambda/\tilde{h}$ ensures the same effective resolution of the Rossby radius for all $n$. When $\lambda/\tilde{h} = 2.0$, the spectrum is well represented, with the exception of the CD mode that occurs at $kh = l\tilde{h} = \pi$ (this is a common feature of all known compatible Galerkin discretizations). When $\lambda/\tilde{h} = 0.1$, the short wavelength part of the spectrum is again underestimated, and has a group velocity with the wrong sign; in a region that grows smaller as $n$ increases. However, unlike the $Q_n^-\Lambda^k$ family, there are no spectral gaps for any $n$ and the maximal frequency does not increase with $n$. Therefore, the CFL condition is not expected to grow worse as $n$ increases. Using reduced quadrature (not shown) does not materially alter these conclusions. Importantly,
there are still no spectral gaps and the maximum frequency is still insensitive to the value of $n$.

Comparing the $Q_n^-\Lambda^k$ (Figures 7 and 8) and $MGD_n$ families (Figures 15 and 16) for the case $\lambda/\tilde{h} = 0.1$, the long wavelength part of the dispersion relationship for the $Q_n^-r\Lambda^k$ family is somewhat better, especially near the end of the spectrum before the underestimated branches. However, the $Q_n^-\Lambda^k$ family has a significant spectral gap in the shortest wavelength part of the spectrum that the $MGD_n$ family avoids. It is unclear which of the two families will perform better in practice, and it seems likely that both will do poorly, just like C-grid finite difference models with a poorly resolved Rossby radius [17, 18, 30, 32, 36].

# 7. Simulations

To investigate the issues that arise due to the presence of line discontinuities in the dispersion relationship, and to compare the $Q_n^-\Lambda^k$ and $MGD_n$ families, two different simulations were performed. Before proceeding to these tests, let us take a slight digression into a discussion of effective resolution.

## 7.1. Effective Resolution

A measure of effective resolution is obtained by computing the fractional error $\frac{\omega_N - \omega_{AN}}{\omega_{AN}}$ in the dispersion relationship as a function of non-dimensional wavelengths $\tilde{\zeta}_x$ and $\tilde{\zeta}_y$. The calculation was performed for $\lambda/\tilde{h} = 2$ using exact quadrature with $n = 1, \dots, 4$, in the case of the $Q_n^-\Lambda^k$ family (Figure 17) and $n = 1, 3, 5, 7$, in the case of the $MGD_n$ family (Figure 18). An error level of 1% was used, which is indicated with the thick black line.

It is clear that for both families, the largest gains in effective resolution come from increasing $n$ from 1 to 3, giving an effective maximal resolution of around $4\tilde{h}$ in both $x$ and $y$, with marginal gains after that. However, for the $MGD_n$ family, despite the marginal gains in effective resolution as $n$ increases, the maximal error decreases significantly and by $n = 7$ it is around 5%. Similar results are obtained for $\lambda/\tilde{h} = 0.1$ (not shown) and when reduced quadrature is used for the $MGD_n$ family (not shown). In fact, reduced quadrature for $MGD_n$ gives a slight increase in effective resolution. These results are similar to others found in the literature [13, 15, 16, 38, 39] for a variety of finite-difference, finite-element and finite-volume schemes.

## 7.2. Test Cases

Two separate test cases were used to investigate the influence of the line discontinuities and compare the $Q_n^-\Lambda^k$ and $MGD_n$ families. For both test cases, the domain $\Omega = [0, L]^2$ with $N = 120/d$ elements in each direction and $h = L/N$, which ensures $120^2$ degrees of freedom independent of $d$; where the parameters are: $L = 2$ m, $\eta_0 = 1$ m, $\alpha = 1/60$, $x_0 = L/2$, $g = 1 \text{ ms}^{-1}$, $H = 1$ m and $f = 25 \text{ s}^{-1}$. This gives $\lambda = 1/25$, and $\lambda/\tilde{h} = 2.4$, corresponding to a well-resolved Rossby radius. These parameters fit with the tests run in [13, 27, 35]. The time integrator is implicit midpoint, with $\Delta t = 1/120$ s and taking 100 steps for a total time of $T = 100/120$ s.
### 7.2.1. Unsupported Test Case

The first test case sets

$$ \eta = \eta_0 \exp \left[ - \frac{(x - x_0)^2 + (y - y_0)^2}{\alpha^2} \right] \quad (106) $$

and $u = v = 0$. This test was run for $n = 1, \dots, 6$, with exact quadrature for the $Q_n^- \Lambda^k$ family and for $n = 1, 3, 5, 7$, with exact and reduced quadrature for the $MGD_n$ family. Additionally, a high resolution run using exact quadrature and $MGD_7$ with $N = 720$ (giving a total of $720^2$ degrees of freedom) and $\Delta T = 1/720$ s for 600 steps was performed. The initial condition along with the high resolution results can be found in Figure 19.

Results for the $Q_n^- \Lambda^k$ family are found in Figure 20, and those for the $MGD_n$ family with exact quadrature are found in Figure 21. The main features of the solution are a steady geostrophic mode in the center of the domain and a propagating, annular inertia-gravity mode, which both the $Q_n^- \Lambda^k$ and $MGD_n$ families were able to capture. However, for $n \ge 2$ the $Q_n^- \Lambda^k$ family shows significant noise in the interior of the annular inertia-gravity mode, along with noise ahead of the wave mode along the $x$ and $y$ axes. There is also significant dispersion, and an underestimation of the magnitude of the propagating mode. These problems lesson somewhat as $n$ increases, but they do not go away. In contrast, the $MGD_n$ family does not have noise at any $n$, captures the amplitude of the propagating mode better and shows significantly lower dispersion, especially at higher $n$. Similar results were obtained using reduced quadrature for the $MGD_n$ family (not shown). A possible hypothesis for this difference in behavior is as follows: the initial condition projects onto a range of spatial wavenumbers, concentrated in the lower-frequency part of the spectrum. For the $Q_n^- \Lambda^k$ family this will include the line discontinuities. These discontinuities lead to noise and enhanced dispersive errors, especially in the high-frequency part of the spectrum. As $n$ is increased, the problematic modes are increasingly localized to the high-frequency part of the spectrum, which has less energy for this initial condition and therefore the problems lessen but are not eliminated. This motivates the second test case, which is designed to study the influence of the line discontinuities directly.

### 7.2.2. Supported Test Case

This test case starts with the initial condition for $\eta$ from the unsupported test case and modulates it by a sinusoidal function, designed to force at the frequency of the line discontinuities:

$$ \eta = \eta_0 \exp \left[ - \frac{(x - x_0)^2 + (y - y_0)^2}{\alpha^2} \right] \cos(k(x - x_0)) \cos(l(y - y_0)) \quad (107) $$

where $k = l = j\pi/(n\tilde{h})$. The test was run for $n = 5$ and $j = 1, 2, 3, 4$ with exact quadrature for the $Q_n^- \Lambda^k$ family and for $n = 5$ with exact and reduced quadrature for the $MGD_n$ family.

Results for the $Q_n^- \Lambda^k$ family are found in Figure 22, and those for the $MGD_n$ family in Figure 23. Those for reduced quadrature with the $MGD_n$ family are very similar to Figure 23 (not shown). Comparing the two families, it is clear that the $Q_n^- \Lambda^k$ family has
significantly more noise and anomalous dispersion than the $MGD_n$ family, even for $j = 1$ where the line discontinuity is very small. It appears than even small discontinuities, such as those that occur in the longer wavelength part of the spectrum for large $n$, can cause issues for the $Q_n^-\Lambda^k$ family. This is anticipated to be particularly relevant for more realistic geophysical fluid models, which have forcings at range of scales due to nonlinear processes, topography and physics parameterizations. It seems likely that some of these issues could be ameliorated through the additional of numerical dissipation, but this is undesirable since such dissipation is not physically motivated.

## 8. Conclusions

This work has shown that the 2D extension of the $P_n^C - P_{n-1}^{DG}$ pair, the $Q_n^-\Lambda^k$ family, shares its problems: spectral gaps and increasing maximal frequency with $n$, all of which grow worse with increasing $n$. When $n = 2$, the spectral gaps can be removed via mass lumping [26, 35], but this degrades the order of convergence and still does not fix the high maximal frequency. In 1D, the alternative pair $GD_n - DGD_{n-1}$ avoids these issues. We have shown that the 2D extension of this pair, the $MGD_n$ family, also avoids the issues. In particular, for any $n$ the $MGD_n$ family pair is free of spectral gaps and the maximal frequency does not increase with $n$. It does have increased computational cost relative to the $Q_n^-\Lambda^k$ family, but this can be partially ameliorated through reduced quadrature without changing any of the features of the dispersion relationship or other desirable properties; except for a loss in order of convergence from $2n$ to $n + 1$. This is still less than the loss seen in the mass lumped $Q_n^-\Lambda^k$ family at $n = 2$. Simulation results support these conclusions, with the $MGD_n$ family producing significantly better results than the $Q_n^-\Lambda^k$ family when $n \ge 2$. Therefore, it appears that the $MGD_n$ family is a strong candidate for the development of geophysical fluid flow models. Study of the effective resolution suggests that $n = 3$ is a sweet spot in terms of maximizing resolvable wavenumber and order of convergence per unit computational cost. Following this idea, a discretization of the thermal shallow water equations using the $MGD_3$ family can be found in [12] and work is ongoing to extend this to the fully compressible Euler equations using Eulerian vertical coordinates. An extension of this study of 2D shallow water dispersion relationships to incorporate time discretization is also underway.

## 9. Acknowledgements

The authors appreciate the thorough and detailed reviews of two anonymous reviewers. Christopher Eldred was supported by the French National Research Agency through contract ANR-14-CE23-0010 (HEAT).

[1] M. Ainsworth. Dispersive behaviour of high order finite element schemes for the one-way wave equation. *Journal of Computational Physics*, 259(Supplement C):1–10, 2014.

[2] M. Ainsworth and H. A. Wajid. Dispersive and dissipative behavior of the spectral element method. *SIAM Journal on Numerical Analysis*, 47(5):3910–3937, 2009.
[3] A. Arakawa and V. R. Lamb. Computational design of the basic dynamical processes of the ucla general circulation model. *General Circulation Models of the Atmosphere*, 17(Supplement C):173–265, 1977.

[4] D. Arnold, R. Falk, and R. Winther. Finite element exterior calculus: from hodge theory to numerical stability. *Bulletin of the American mathematical society*, 47(2):281–354, 2010.

[5] D. N. Arnold, R. S. Falk, and R. Winther. Finite element exterior calculus, homological techniques, and applications. *Acta Numerica*, 15:1155, 2006.

[6] J. Banks and T. Hagstrom. On galerkin difference methods. *Journal of Computational Physics*, 313(Supplement C):310–327, 2016.

[7] A. Buffa, J. Rivas, G. Sangalli, and R. Vazquez. Isogeometric discrete differential forms in three dimensions. *SIAM Journal on Numerical Analysis*, 49(2):818–844, 2011.

[8] A. Buffa, G. Sangalli, and R. Vazquez. Isogeometric methods for computational electromagnetics: B-spline and t-spline discretizations. *Journal of Computational Physics*, 257(Part B):1291–1320, 2014. Physics-compatible numerical methods.

[9] C. Cotter and J. Shipton. Mixed finite elements for numerical weather prediction. *Journal of Computational Physics*, 231(21):7076–7091, 2012.

[10] C. J. Cotter, D. A. Ham, and C. C. Pain. A mixed discontinuous/continuous finite element pair for shallow-water ocean modelling. *Ocean Modelling*, 26(1):86–90, 2009.

[11] S. Danilov. On utility of triangular c-grid type discretization for numerical modeling of large-scale ocean flows. *Ocean Dynamics*, 60(6):1361–1369, Dec 2010.

[12] C. Eldred, T. Dubos, and E. Kritsikis. A Quasi-Hamiltonian Discretization of the Thermal Shallow Water Equations. working paper or preprint, Oct. 2018.

[13] C. Eldred and D. Y. Le Roux. Dispersion analysis of compatible galerkin schemes for the 1d shallow water model. *Journal of Computational Physics*, 371:779–800, 2018.

[14] R. Hiemstra, D. Toshniwal, R. Huijsmans, and M. Gerritsma. High order geometric methods with exact conservation properties. *Journal of Computational Physics*, 257(Part B):1444–1471, 2014. Physics-compatible numerical methods.

[15] J. Kent, C. Jablonowski, J. P. Whitehead, and R. B. Rood. Determining the effective resolution of advection schemes. part ii: Numerical testing. *Journal of Computational Physics*, 278(Supplement C):497–508, 2014.

[16] J. Kent, J. P. Whitehead, C. Jablonowski, and R. B. Rood. Determining the effective resolution of advection schemes. part i: Dispersion analysis. *Journal of Computational Physics*, 278(Supplement C):485–496, 2014.
[17] C. S. Konor and D. A. Randall. Impacts of the horizontal and vertical grids on the numerical solutions of the dynamical equations – part 1: Nonhydrostatic inertia–gravity modes. *Geoscientific Model Development*, 11(5):1753–1784, 2018.

[18] C. S. Konor and D. A. Randall. Impacts of the horizontal and vertical grids on the numerical solutions of the dynamical equations – part 2: Quasi-geostrophic rossby modes. *Geoscientific Model Development*, 11(5):1785–1797, 2018.

[19] D. Y. Le Roux. Spurious inertial oscillations in shallow-water models. *Journal of Computational Physics*, 231(24):7959–7987, 2012.

[20] D. Y. Le Roux and B. Pouliot. Analysis of numerically induced oscillations in two-dimensional finite-element shallow-water models part ii: Free planetary waves. *SIAM Journal on Scientific Computing*, 30(4):1971–1991, 2008.

[21] D. Y. Le Roux, V. Rostand, and B. Pouliot. Analysis of numerically induced oscillations in 2d finite-element shallow-water models part i: Inertia-gravity waves. *SIAM Journal on Scientific Computing*, 29(1):331–360, 2007.

[22] P. H. LeBlond and L. A. Mysak. *Waves in the Ocean*, volume 20 of *Elsevier Oceanography Series*. Elsevier, 1978.

[23] A. McRae, G. Bercea, L. Mitchell, D. Ham, and C. Cotter. Automated generation and symbolic manipulation of tensor product finite elements. *SIAM Journal on Scientific Computing*, 38(5):S25–S47, 2016.

[24] T. Melvin. Dispersion analysis of the pnpn1dg mixed finite element pair for atmospheric modelling. *Journal of Computational Physics*, 355(Supplement C):342–365, 2018.

[25] T. Melvin, T. Benacchio, J. Thuburn, and C. Cotter. Choice of function spaces for thermodynamic variables in mixed finite-element methods. *Quarterly Journal of the Royal Meteorological Society*, 144(712):900–916, 2018.

[26] T. Melvin, A. Staniforth, and C. Cotter. A two-dimensional mixed finite-element pair on rectangles. *Quarterly Journal of the Royal Meteorological Society*, 140(680):930–942, 2014.

[27] T. Melvin, A. Staniforth, and J. Thuburn. Dispersion analysis of the spectral element method. *Quarterly Journal of the Royal Meteorological Society*, 138(668):1934–1947, 2012.

[28] T. Melvin and J. Thuburn. Wave dispersion properties of compound finite elements. *Journal of Computational Physics*, 338:68–90, 2017.

[29] A. Natale, J. Shipton, and C. J. Cotter. Compatible finite element spaces for geophysical fluid dynamics. *Dynamics and Statistics of the Climate System*, 1(1):dzw005, 2016.
[30] S. Nikovi, M. B. Gavrilov, and I. A. Toi. Geostrophic adjustment on hexagonal grids. *Monthly Weather Review*, 130(3):668-683, 2002.

[31] V. Prasolov. *Problems and Theorems in Linear Algebra*. Translations of mathematical monographs. American Mathematical Society, 1994.

[32] D. A. Randall. Geostrophic adjustment and the finite-difference shallow-water equations. *Monthly Weather Review*, 122(6):1371-1377, 1994.

[33] V. Rostand and D. Y. Le Roux. Raviart-thomas and brezzi-douglas-marini finite-element approximations of the shallow-water equations. *International Journal for Numerical Methods in Fluids*, 57(8):951-976, 2008.

[34] V. Rostand, D. Y. Le Roux, and G. Carey. Kernel analysis of the discretized finite difference and finite element shallow-water models. *SIAM Journal on Scientific Computing*, 31(1):531-556, 2008.

[35] A. Staniforth, T. Melvin, and C. Cotter. Analysis of a mixed finite-element pair proposed for an atmospheric dynamical core. *Quarterly Journal of the Royal Meteorological Society*, 139(674):1239-1254, 2013.

[36] J. Thuburn. Numerical wave propagation on the hexagonal c-grid. *Journal of Computational Physics*, 227(11):5836-5858, 2008.

[37] J. Thuburn and C. J. Cotter. A primaldual mimetic finite element scheme for the rotating shallow water equations on polygonal spherical meshes. *Journal of Computational Physics*, 290(Supplement C):274-297, 2015.

[38] P. A. Ullrich. Understanding the treatment of waves in atmospheric models. part 1: The shortest resolved waves of the 1d linearized shallow-water equations. *Quarterly Journal of the Royal Meteorological Society*, 140(682):1426-1440, 2014.

[39] P. A. Ullrich, D. R. Reynolds, J. E. Guerra, and M. A. Taylor. Impact and importance of hyperdiffusion on the spectral element method: A linear dispersion analysis. *Journal of Computational Physics*, 375:427-446, 2018.
<figure><img src="image_6.png" /><figcaption>Figure 5: Contour plots of the numerical dispersion relationship $\sigma = \frac{\omega}{f}$ versus $(\tilde{kh}, \tilde{l}h)$ for $n = 1, \dots, 4$ using the $Q_n^{-} \Lambda^k$ family and exact quadrature when the Rossby radius is well resolved ($\lambda/\tilde{h} = 2.0$). Note the presence of spectral gaps (in the form of line discontinuities) at $\tilde{kh} = \frac{j\pi}{n}$ ($j = 1, \dots, n-1$) and $\tilde{l}h = \frac{m\pi}{n}$ ($m = 1, \dots, n-1$), and the increasing maximal frequency occurring at large $\tilde{kh}$ and $\tilde{l}h$ as $n$ increases.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 6: Surface plots of the numerical dispersion relationship $\sigma = \frac{\omega}{f}$ versus $(kh, \tilde{lh})$ for $n = 1, \dots, 4$ using the $Q_n^{-} \Lambda^k$ family and exact quadrature when the Rossby radius is well resolved ($\lambda/\tilde{h} = 2.0$). Note the presence of spectral gaps (in the form of line discontinuities) at $\tilde{kh} = \frac{j\pi}{n}$ ($j = 1, \dots, n-1$) and $\tilde{lh} = \frac{m\pi}{n}$ ($m = 1, \dots, n-1$), and the increasing maximal frequency occurring at large $kh$ and $\tilde{lh}$ as $n$ increases.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 7: As Figure 5, but for a poorly resolved Rossby radius ($\lambda/\tilde{h} = 0.1$). There are still spectral gaps and an increasing maximal frequency as $n$ increases, but now the short-wavelength part of the spectrum is under-estimated and has a group velocity with the wrong sign.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 8: As Figure 6, but for a poorly resolved Rossby radius ($\lambda/\tilde{h} = 0.1$). There are still spectral gaps and an increasing maximal frequency as $n$ increases, but now the short-wavelength part of the spectrum is under-estimated and has a group velocity with the wrong sign.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 9: The difference $(w_{tp} - w_{direct})/w_{direct} * 100$ between $w_{direct}$ found following Section 3.6, and $w_{tp}$ obtained using (102) for $\lambda/\hbar = 2.0$ and $n = 1, \dots, 4$ using exact quadrature. As expected, there is no difference when $n = 1$, but there is an error of around 1% for $n \ge 2$. This is evidence that (102) is not giving the correct dispersion relationship for $n \ge 2$.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 10: Surface and contour plots of the numerical dispersion relationship $\sigma = \frac{\omega}{f}$ versus $(\tilde{k}h, \tilde{l}h)$ for the lumped version of the $Q_2^-$-$\Lambda^k$ family with exact quadrature when the Rossby radius is well resolved ($\lambda/\tilde{h} = 2.0$) and poorly resolved ($\lambda/\tilde{h} = 0.1$). Note there are no spectral gaps, and the poorly resolved case has a group velocity with the wrong sign.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 11: The Fourier degrees of freedom for the $MGD_n$ family, with $W_0(\psi, \not{v})$ on the left, the two components of $W_1(\not{u}, \not{v})$ in the middle and $W_2(\eta)$ on the right. Note that since $d=1$, there is only one degree of freedom for any $n$.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Figure 12: The basis functions for the $W_0$ (top left), $W_1$ (bottom, $u$ on left and $v$ on right) and $W_2$ (top right) spaces from the $MGD_3$ family. Unlike the $Q_n^- \Lambda^k$ family, the basis functions have extended support beyond neighboring elements. Here a 4 x 4 element mesh with element width equal to 1 is used, and we show the basis function associated with the unique degrees of freedom for the element spanning (2,3) in both $x$ and $y$.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Figure 13: Contour plots of the numerical dispersion relationship $\sigma = \frac{\omega}{f}$ versus $(kh, l\tilde{h})$ for $n = 1, 3, 5, 7$ using the MGD$_n$ family and exact quadrature when the Rossby radius is well resolved ($\lambda/\tilde{h} = 2.0$). Unlike the $Q_n^- \Lambda^k$ family, there are no spectral gaps and the maximal frequency does not increase with $n$.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Figure 14: Surface plots of the numerical dispersion relationship $\sigma = \frac{\omega}{f}$ versus $(kh, l\tilde{h})$ for $n = 1, 3, 5, 7$ using the MGD$_n$ family and exact quadrature when the Rossby radius is well resolved ($\lambda/\tilde{h} = 2.0$). Unlike the $Q_n^-\Lambda^k$ family, there are no spectral gaps and the maximal frequency does not increase with $n$.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Figure 15: As Figure 13, but with a poorly resolved Rossby radius ($\lambda/\tilde{h} = 0.1$). There are still no spectral gaps, but there is a significant underestimation of the frequency in the short wavelength part of the spectrum and a group velocity with the wrong sign.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Figure 16: As Figure 13, but with a poorly resolved Rossby radius ($\lambda/\tilde{h} = 0.1$). There are still no spectral gaps, but there is a significant underestimation of the frequency in the short wavelength part of the spectrum and a group velocity with the wrong sign.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>Figure 17: Plots of the fractional error $\frac{\omega_N - \omega_{AN}}{\omega_{AN}}$ as a function of non-dimensional wavelengths $\tilde{\zeta}_x$ and $\tilde{\zeta}_y$ for the $Q_n^- \Lambda^k$ family using exact quadrature with $n = 1, \dots, 4$ and $\lambda/\tilde{h} = 2$. The thick black line indicates an error level of 0.01 or 1%, and can be used to diagnose an effective resolution. The largest gains in effective resolution occur as $n$ increases from 1 to 3, giving around $4\tilde{h}$ in both $x$ and $y$.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>Figure 18: As for Figure 17, but using the MGD$_n$ family with $n = 1, \dots, 7$. Again, the largest gain in effective resolution comes from increasing $n$ from 1 to 3, giving around $4\tilde{h}$ in both $x$ and $y$. Note also that the maximum error decreases quickly as $n$ increases, and it is around 5% for $n = 7$.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>Figure 19: The initial condition (left) $\eta(0)$ and the high resolution reference solution (right) at the end of the simulation $\eta(T)$ for the unsupported test case. The reference solution was computed using $MGD_7$ with $N = 720$ ($720^2$ degrees of freedom total) and $\Delta T = 1/720$ s (600 steps total).</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>Figure 20: The final height $\eta(T)$ for the unsupported test case using the $Q_n^- \Lambda^k$ family with exact quadrature for $n = 1, \dots, 6$. There is significant noise and anomalous dispersion for all $n \ge 1$, with the effects somewhat improving at $n$ increases but not going away. In particular, the anomalous dispersion along the $x$ and $y$ axes remains strong, even at $n = 6$.</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>Figure 21: As for Figure 20, but using the MGD$_n$ family with $n = 1, 3, 5, 7$. Unlike the $Q_n^- \Lambda^k$ family there is essentially no noise for any $n$ and reduced anomalous dispersion, especially for higher $n$.</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>Figure 22: Height $\eta(T)$ for the supported test case using the $Q_n^- \Lambda^k$ family using $n = 5$ and exact quadrature. This test forces specifically at the gap frequency $k = l = j\pi/(n\tilde{h})$, and was run for $j = 1, 2, 3, 4$. For all $j$ there is significant noise and anomalous dispersion, even for the case $j = 1$ when the gap is very small.</figcaption></figure>
<figure><img src="image_24.png" /><figcaption>Figure 23: As Figure 22, but for the MGD$_n$ family with $n = 5$. Unlike the $Q_n^- \Lambda^k$ family, there is no little anomalous dispersion and almost no noise for any $j$.</figcaption></figure>