Christopher Eldred
email: christopher.eldred@inria.fr

Daniel Y Le Roux

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

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

 that eliminates the spectral gaps. Additionally, the M GD 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.

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 [START_REF] Hiemstra | High order geometric methods with exact conservation properties[END_REF], which are a generalization of compatible finite element methods [START_REF] Arnold | Finite element exterior calculus: from hodge theory to numerical stability[END_REF][START_REF] Arnold | Finite element exterior calculus, homological techniques, and applications[END_REF][START_REF] Cotter | Mixed finite elements for numerical weather prediction[END_REF][START_REF] Natale | Compatible finite element spaces for geophysical fluid dynamics[END_REF] 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 ń Λ 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) [START_REF] Cotter | Mixed finite elements for numerical weather prediction[END_REF][START_REF] Le Roux | Spurious inertial oscillations in shallow-water models[END_REF] 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 posses only three branches in the dispersion relationship, but certain discretizations will have extra branches of either Rossby or inertia-gravity (IGW) modes [START_REF] Cotter | Mixed finite elements for numerical weather prediction[END_REF][START_REF] Danilov | On utility of triangular c-grid type discretization for numerical modeling of large-scale ocean flows[END_REF][START_REF] Nikovi | Geostrophic adjustment on hexagonal grids[END_REF][START_REF] Thuburn | Numerical wave propagation on the hexagonal c-grid[END_REF]. 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 Λ k family on quadrilaterals has spurious Rossby modes. For example, the RT 0 ´P DG 0 and BDM 1 ´P DG 0 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 [START_REF] Ainsworth | Dispersive behaviour of high order finite element schemes for the one-way wave equation[END_REF][START_REF] Ainsworth | Dispersive and dissipative behavior of the spectral element method[END_REF][START_REF] Melvin | Dispersion analysis of the spectral element method[END_REF], 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 ń Λ k and M GD n families are built like this), can ensure a lack of spurious branches [START_REF] Cotter | Mixed finite elements for numerical weather prediction[END_REF]. In this work we will show that the M GD n family additionally avoids spectral gaps in the dispersion relationship for any n, unlike the Q ń Λ k family which has gaps for n ě 2.

Dispersion relationships have been studied since the beginning of geophysical fluid modeling [START_REF] Arakawa | Computational design of the basic dynamical processes of the ucla general circulation model[END_REF], initially using finite differences. Within the class of compatible Galerkin methods, the RT 0 ´P DG 0 and/or BDM 1 ´P DG 0 pairs on triangles were studied in [START_REF] Le Roux | Spurious inertial oscillations in shallow-water models[END_REF][START_REF] Le Roux | Analysis of numerically induced oscillations in twodimensional finite-element shallow-water models part ii: Free planetary waves[END_REF][START_REF] Le Roux | Analysis of numerically induced oscillations in 2d finite-element shallow-water models part i: Inertia-gravity waves[END_REF][START_REF] Rostand | Raviart-thomas and brezzi-douglas-marini finiteelement approximations of the shallow-water equations[END_REF][START_REF] Rostand | Kernel analysis of the discretized finite difference and finite element shallow-water models[END_REF]. These are the lowest-order members of the corresponding P ŕ Λ k (RT 0 ) and P r Λ k (BDM 1 ) families on triangles from finite element exterior calculus [START_REF] Arnold | Finite element exterior calculus, homological techniques, and applications[END_REF]. However, both of these families are known to have spurious branches [START_REF] Cotter | Mixed finite elements for numerical weather prediction[END_REF], and the presence or absence of spectral gaps is unclear, although the results from the P C n ´P DG n´1 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 [START_REF] Thuburn | A primaldual mimetic finite element scheme for the rotating shallow water equations on polygonal spherical meshes[END_REF]) were studied for quadrilaterals and hexagons in [START_REF] Melvin | Wave dispersion properties of compound finite elements[END_REF], and compared to the corresponding C-grid finite difference schemes [START_REF] Nikovi | Geostrophic adjustment on hexagonal grids[END_REF][START_REF] Randall | Geostrophic adjustment and the finite-difference shallow-water equations[END_REF][START_REF] Thuburn | Numerical wave propagation on the hexagonal c-grid[END_REF]. 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 C 2 ´P DG 1 pair [START_REF] Cotter | A mixed discontinuous/continuous finite element pair for shallow-water ocean modelling[END_REF] 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 [START_REF] Staniforth | Analysis of a mixed finite-element pair proposed for an atmospheric dynamical core[END_REF], and was extended to the 2D shallow water equations in [START_REF] Melvin | A two-dimensional mixed finite-element pair on rectangles[END_REF]. However, [START_REF] Eldred | Dispersion analysis of compatible galerkin schemes for the 1d shallow water model[END_REF] showed that this form of lumping does not work for n ě 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 C pairs was studied in [START_REF] Melvin | Choice of function spaces for thermodynamic variables in mixed finite-element methods[END_REF], 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 C n ´P DG n´1 pair in 1D were studied in [START_REF] Melvin | Dispersion analysis of the pnpn1dg mixed finite element pair for atmospheric modelling[END_REF], along with 2D gravity waves using the Q ń Λ k family. Finally, the P C n ´P DG n´1 and GD n ´DGD n´1 pairs were studied for inertia-gravity waves in 1D in [START_REF] Eldred | Dispersion analysis of compatible galerkin schemes for the 1d shallow water model[END_REF]. Both [START_REF] Melvin | Dispersion analysis of the pnpn1dg mixed finite element pair for atmospheric modelling[END_REF] and [START_REF] Eldred | Dispersion analysis of compatible galerkin schemes for the 1d shallow water model[END_REF] showed that the P C n ´P DG n´1 pair suffers from spectral gaps when n ě 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 [START_REF] Eldred | Dispersion analysis of compatible galerkin schemes for the 1d shallow water model[END_REF] 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 [START_REF] Eldred | Dispersion analysis of compatible galerkin schemes for the 1d shallow water model[END_REF] from the 1D shallow water equations to the 2D shallow water equations on quadrilaterals, and investigate the dispersion properties of the Q ń Λ k and M GD n families, which are the 2D extensions of the P C n ´P DG n´1 and GD n DGD n´1 pairs. Additionally, we clarify a partially erroneous claim found at the end of [START_REF] Melvin | Dispersion analysis of the pnpn1dg mixed finite element pair for atmospheric modelling[END_REF], 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 ě 2), such as the Q ń Λ k family when n ě 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 ń Λ k and M GD n families are studied in Section 5 and Section 6. This includes asympototics and a dispersion plot for the Q 2 Λ k family using the mass lumping from [START_REF] Melvin | A two-dimensional mixed finite-element pair on rectangles[END_REF][START_REF] Staniforth | Analysis of a mixed finite-element pair proposed for an atmospheric dynamical core[END_REF] that eliminates the spectral gaps. It is shown that the Q ń Λ k family suffers from spectral gaps when n ě 2, in the form of line discontinuities, while the M GD 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 threedimensional domains with boundaries and variable Coriolis force and the incorporation of time discretization remains to be done.

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 [START_REF] Leblond | Waves in the Ocean[END_REF] as

Bu Bt ´f v `g Bη Bx " 0 , (1)
Bv Bt `f u `g Bη By " 0 , (2)
Bη Bt `H ˆBu Bx `Bv By ˙" 0 , (3)
with relative velocity vpx, y, tq " pu, vq and surface height ηpx, y, tq. 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 ω and spatial wavenumbers k and l, a single Fourier mode is written as η " ηe ´iωt e ikx e ily u " ûe ´iωt e ikx e ily v " ve ´iωt e ikx e ily

where η, û and v are the Fourier amplitudes. Substituting (4) into (1) -(3) gives

´iωû ´f v `ikg η " 0, (5)
´iωv `f û `ilg η " 0, ( 6)

´iω η `ikH û `ilH v " 0. ( 7
)
For non-zero (η, û, v), ( 5) -( 7) can be solved and give a dispersion polynomial of degree 3 in ω. This polynomial has 3 solutions: the geostrophic mode ω " 0 and two inertia-gravity modes:

ω 2 AN " f 2 `gHpk 2 `l2 q " f 2 `gHK 2 (8)
where K 2 " k 2 `l2 . The temporal frequency ω AN can be normalized as σ AN " ω AN f , which gives

σ 2 AN " 1 `λ2 K 2 (9)
where λ 2 " gH f 2 is the squared Rossby radius of deformation λ. Finally, introducing a length scale h to facilitate comparison with numerical results, we have

σ 2 AN " 1 `µp Ą Khq 2 , (10)
with µ " λ 2 h 2 and p Ą Khq 2 " p Ă khq 2 `p r lhq 2 is the non-dimensional wavenumber ( Ă kh P r0, πs and r lh P r0, πs). Plots of σ AN pkh, lhq for λ{ h " 2.0 and λ{ h " 0.1 are found in Figures 1 and2.

Since ω 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 ń Λ k family in [START_REF] Melvin | Dispersion analysis of the pnpn1dg mixed finite element pair for atmospheric modelling[END_REF], and both pure inertia and pure gravity waves for the Q 2 Λ k family are discussed in [START_REF] Melvin | A two-dimensional mixed finite-element pair on rectangles[END_REF]. Figure 1: Plots of the analytic dispersion relationship σ AN versus non-dimensional wavenumbers Ă kh and r lh for λ{ h " 2.0 (left) and λ{ h " 0.1 (right). Note that σ AN is monotonic and non-zero for all Ă kh, r lh ą 0.

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. kh and r lh for λ{ h " 2.0 (left) and λ{ h " 0.1 (right). Note that σ AN is monotonic and non-zero for all Ă kh, r lh ą 0.

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

2D Spaces

Consider a pair of 1D spaces A Ă H 1 and B Ă L 2 that are a partition of unity and form the 1D discrete deRham complex A B{Bx ÝÝÑ B. These are precisely the class of spaces studied in [START_REF] Eldred | Dispersion analysis of compatible galerkin schemes for the 1d shallow water model[END_REF]. A 2D discrete deRham complex is formed by taking tensor products of the 1D spaces:

W 0 " A b A Ă H 1 (11) W 1 " pA b Bq î `pB b Aq ĵ Ă Hpdivq (12) W 2 " B b B Ă L 2 (13)
where î and ĵ 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 [START_REF] Eldred | A Quasi-Hamiltonian Discretization of the Thermal Shallow Water Equations[END_REF][START_REF] Hiemstra | High order geometric methods with exact conservation properties[END_REF][START_REF] Mcrae | Automated generation and symbolic manipulation of tensor product finite elements[END_REF][START_REF] Natale | Compatible finite element spaces for geophysical fluid dynamics[END_REF]. In this paper we consider the Q ń Λ k family from finite element exterior calculus with A " P C n and B " P DG n´1 [START_REF] Arnold | Finite element exterior calculus: from hodge theory to numerical stability[END_REF][START_REF] Arnold | Finite element exterior calculus, homological techniques, and applications[END_REF] and the M GD n family (mimetic Galerkin differences) with A " GD n and B " DGD n´1 [START_REF] Eldred | A Quasi-Hamiltonian Discretization of the Thermal Shallow Water Equations[END_REF][START_REF] Eldred | Dispersion analysis of compatible galerkin schemes for the 1d shallow water model[END_REF]. These are the extension of the element pairs in [START_REF] Eldred | Dispersion analysis of compatible galerkin schemes for the 1d shallow water model[END_REF] to quadrilaterals. Recall that P C n is the order n Lagrange space, P DG n´1 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 [START_REF] Hiemstra | High order geometric methods with exact conservation properties[END_REF] and isogeometric discrete differential forms [START_REF] Buffa | Isogeometric discrete differential forms in three dimensions[END_REF][START_REF] Buffa | Isogeometric methods for computational electromagnetics: B-spline and t-spline discretizations[END_REF][START_REF] Hiemstra | High order geometric methods with exact conservation properties[END_REF]. 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.

Discrete Model

Choose η P W 2 and v " pu, vq P W 1 . This means

u P A b B v P B b A η P B b B ( 14
)
which corresponds with the placement in 1D of η, v P B and u P A [START_REF] Eldred | Dispersion analysis of compatible galerkin schemes for the 1d shallow water model[END_REF]. Given test functions η ˚P W 2 and v ˚" pu ˚, v ˚q P W 1 , (1) -( 3) are discretized as (after appropriate integration by parts)

u ˚, Bu Bt ´f u ˚, v ´g B Bx u ˚, η " 0, (15)
v ˚, Bv Bt `f v ˚, u ´g B By v ˚, η " 0, (16)
η ˚, Bη Bt `H η ˚, B Bx u `H η ˚, B By v " 0, ( 17
)
where a, b is the L 2 inner product. Written in matrix form, (15) - [START_REF] Konor | Impacts of the horizontal and vertical grids on the numerical solutions of the dynamical equations -part 1: Nonhydrostatic inertia-gravity modes[END_REF] are

M u Bu Bt ´f C u v `g G u η " 0, (18)
M v Bv Bt `f C v u `g G v η " 0, (19)
M η Bη Bt `H D u u `H D v v " 0, (20)
where u, v and η denote the basis coefficients in the Galerkin expansion and the matrices are given by

M u Bu Bt " u ˚, Bu Bt , M v Bv Bt " v ˚, Bv Bt , M η Bη Bt " η ˚, Bη Bt , G u η " ´ B Bx u ˚, η , G v η " ´ B By v ˚, η , D u u " η ˚, B Bx u , D v v " η ˚, B By v , C u v " u ˚, v , C v u " v ˚, u .
Use of compatible Galerkin spaces ensures

M u " M u ˚, M v " M v ˚, M η " M η ˚, (21) G u ˚" ´Du , G v ˚" ´Dv , (22)
C u ˚" C v . ( 23
)
with A ˚the Hermitian transpose. These properties are inherited by the matrices r A and p A that arise after performing Fourier analysis. The properties ( 21) -( 23) will be useful in proving some features of the dispersion relationship.

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 ń Λ k family, and d " 1 for the M GD 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 pN dq 2 ˆpN dq 2 . Assume a Fourier dependence in space and time for the basis coefficients: u j ptq " ũj e ´iωt e ikx j e ily j , v j ptq " ṽj e ´iωt e ikx j e ily j , η j ptq " ηj e ´iωt e ikx j e ily j ,

where (x j , y j ) is the nodal location of the jth basis function for the relevant variable and (ũ j , ṽj , ηj ) are the Fourier amplitudes, with j " 1, . . . , pN dq 2 ,. Inserting ( 24) into ( 18) -( 20) yields

´iω Ą M u ũ ´f Ă C u ṽ `g Ă G u η " 0, ( 25
)
´iω Ą M v ṽ `f Ă C v ũ `g Ă G v η " 0, ( 26
)
´iω Ą M η η `H Ă D u ũ `H Ă D v ṽ " 0, (27)
where ũ " " ũ1 , . . . ,

ũpNdq 2 ‰ , ṽ " " ṽ1 , . . . , ṽpNdq 2 ‰ , η " " η1 , . . . , ηpNdq 2 ‰ . (28)
The relationship between A and r A is given by r A mn " A mn e ikpxm´x 0 q e ilpym´y 0 q ,

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 ( N 2 , N 2 ). Using the translational invariance of the basis functions, there are actually only d 2 unique degrees of freedom for ũ, ṽ and η: û " rû 1 , . . . , ûd 2 s T , v " rv 1 , . . . , vd 2 s T , η " rη 1 , . . . , ηd 2 s T .

(

) 30
Thus we can write (28) as ũ " rû 1 , . . . , ûd 2 , û1 , . . . , ûd 2 , . . . s ,

and similarly for ṽ and η. 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ω y M u û ´f x C u v `g x G u η " 0, ( 32
)
´iω y M v v `f x C v û `g x G v η " 0, ( 33
)
´iω y M η η `H x D u û `H x D v v " 0, (34)
The matrices r A and p A are related by

p A mn " ÿ p r A mp , (35)
with m " 1, . . . , d 2 , n " 1, . . . , d 2 , and p is the set of integers belonging to the interval r1, . . . , pN dq 2 s such that xp " xn . The final matrices p A will have the same properties ( 21) -( 23) as the original matrices 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ω x Mx " p Sx, (36)
where x " rû 1 , . . . , ûd 2 , v1 , . . . , vd 2 , η1 , . . . , ηd 2 s T and

x M " ¨y M u 0 0 0 y M v 0 0 0 y M η ‹ ', p S " ¨0 f x C u ´g x G u ´f x C v 0 ´g x G v ´H x D u ´H x D v 0 ‹ '. (37)
Both x M and p S are 3 ˆ3 block matrices of total size 3d 2 ˆ3d 2 , with individual blocks of size d 2 ˆd2 . 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 [START_REF] Thuburn | Numerical wave propagation on the hexagonal c-grid[END_REF] reduces to the eigenvalue problem (since

x M is invertible) ´iω I x " x M ´1 p Sx " p Fx (38)
The eigenvalues of p F " x M ´1 p 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 [START_REF] Thuburn | Numerical wave propagation on the hexagonal c-grid[END_REF]. 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 ω is real.

Proof. Start with [START_REF] Thuburn | Numerical wave propagation on the hexagonal c-grid[END_REF], and note that it is equivalent to

ω x Mx " i p Sx, (39)
Now replace η with b H g η in order to symmetrize the terms in p S involving g and

H p S " ¨0 f x C u ´?gH x G u ´f x C v 0 ´?gH x G v ´?gH x D u ´?gH x D v 0 ‹ ' (40)
where we have used the same symbols η and p S to ease notation. Using ( 21) -( 23), we obtain

x M " x M ˚, (41)
and

p S " ´p S ˚. ( 42
)
and thus pi p Sq ˚" i p S. Since both x M and i p S are Hermitian matrices, the revised generalized eigenvalue problem [START_REF] Ullrich | Impact and importance of hyperdiffusion on the spectral element method: A linear dispersion analysis[END_REF] has real eigenvalues, and thus ω is real.

Conjecture 1. The 3d 2 eigenvalues are divided into d 2 geostrophic modes ω " 0 and d 2 pairs of inertia-gravity waves: d 2 positive modes ω pos,i and d 2 negative modes ω neg,i with ω neg,i " ´ωpos,i (where i " 1, . . . , 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 [START_REF] Eldred | Dispersion analysis of compatible galerkin schemes for the 1d shallow water model[END_REF]. However, we have been unable to find a proof of this conjecture using the matrix properties of x M and p S, although presumably one must exist. A closely related statement has been proven in [START_REF] Cotter | Mixed finite elements for numerical weather prediction[END_REF] 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.

Non-Dimensional Spatial Wavenumbers

It is useful to define an adjusted grid spacing

h " h d , (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 h can be viewed as an average distance. Using h, non-dimensional wavenumbers are naturally defined as

Ă kh " k h, r lh " l h, (44)
with the associated wavelengths

ζ x " 2π k " 2π h Ă kh , ζ y " 2π l " 2π h r lh , (45)
and non-dimensional wavelengths

r ζ x " ζ x h " 2π Ă kh , r ζ y " ζ y h " 2π r lh . ( 46
)
For a grid with N by N elements, the maximal wavelengths are

ζ x " ζ y " N d h " N h, with k " l " 2π N h , Ă kh " r lh " 2π N d and r ζ x " r ζ y " N d;
and the minimum wavelengths are

ζ x " ζ y " 2 h, with k " l " π h , Ă kh " r lh " π and r ζ x " r ζ y " 2.
These are useful in the definition of effective resolution, and also in determining the relevant wavenumber ranges to plot over.

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 ( Ă kh, r lh), when d ě 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 [START_REF] Eldred | Dispersion analysis of compatible galerkin schemes for the 1d shallow water model[END_REF]. 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.

Solution Procedure

Following [START_REF] Eldred | Dispersion analysis of compatible galerkin schemes for the 1d shallow water model[END_REF], two complementary approaches are used to solve [START_REF] Thuburn | Numerical wave propagation on the hexagonal c-grid[END_REF]. 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 σpkh, lh, µq. 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 ˆ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 ń Λ k and M GD 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 " 60 d , the eigenvalues and eigenvectors [START_REF] Thuburn | Numerical wave propagation on the hexagonal c-grid[END_REF] are determined numerically for each pkh, lhq pair and µ of interest. The set of pkh, lhq sampled is a uniformly spaced set belonging to r 2π N d , πs 2 . When d ě 2 (which occurs for the Q ń Λ k family with n ě 2), there are d 2 eigenvalues. Then, using the procedure in Section 3.7, the solution is determined from these d 2 eigenvalues.

Branch determination for d ě 2

The following procedure is used to select the appropriate branch (it is the 2D generalization of the approach in [START_REF] Eldred | Dispersion analysis of compatible galerkin schemes for the 1d shallow water model[END_REF]):

1. Let pkh, lhq P p0, πs 2 and numerically compute the

d 2 eigenvalues ω m (m " 1, . . . , d 2 )
and eigenvectors xm .

2. For each eigenvalue and eigenvector pair, determine the branch (p x , p y ) that it belongs to, with p x " 1, . . . , d and p y " 1, . . . , d. This is done by constructing the discrete solution corresponding to xm and comparing it to the continuous solution e ikx e ily .

3. Determine the effective wavenumber ( Ă kh, r lh) associated with ω px,py using the explicit formula

Ă kh " p´1q αx pkhq `2πβ x d , r lh " p´1q αy plhq `2πβ y d , (47)
where α x " 1 `px , α y " 1 `py , β x " t px 2 u and β y " t py 2 u. This gives an effective wavenumber range of p0, πs 2 .

Exploiting Tensor Product Structure

In [START_REF] Melvin | Dispersion analysis of the pnpn1dg mixed finite element pair for atmospheric modelling[END_REF] it is claimed without proof (see equations ( 80) -(85) in that paper) that the 2D IW, GW and IGW dispersion relationships for the Q ń Λ k family (and other finite element type compatible Galerkin schemes with multiple degrees of freedom when n ě 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 ě 2 we provide counterexamples for the Q ń Λ 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 x M and p S in terms of the 1D matrices.

Kronecker Products

Before continuing, it it useful to recall some basic properties of the Kronecker (or tensor) product [START_REF] Prasolov | Problems and Theorems in Linear Algebra[END_REF]:

• Hermitian Transpose: pA b Bq ˚" A ˚b B ˚(48)
• Determinant: If A is an n x n matrix and B is an m x m matrix then

| A b B | " | A | m | B | n (
)
where I k is the k x k identity matrix.

Eigenvalues of Kronecker Products and Sums

Of particular interest are the eigenvalues of Kronecker products and sums. Let α i denote the eigenvalues of A and β j the eigenvalues of B. Then the eigenvalues of the Kronecker product

C " A b B (53) are c ij " α i β j (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 ' B " A b I `I b B (55) are δ ij " α i `βj (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.

1D IGW Problem

Start by recalling the 1D IGW problem from [START_REF] Eldred | Dispersion analysis of compatible galerkin schemes for the 1d shallow water model[END_REF].

x M " ¨y M x C 0 0 0 y M x D 0 0 0 y M x D ‹ ', p S " ¨0 f x C x u ´g x G x ´f x C x v 0 0 ´H x D x 0 0 ‹ ' (57)
where we have added the superscript x to indicate that the spatial wavenumber is k. The submatrices y [START_REF] Eldred | Dispersion analysis of compatible galerkin schemes for the 1d shallow water model[END_REF], and they satisfy

M x C , y M x D , x G x , x D x , x C x u and x C x v are defined in
y M x C ˚" y M x C , y M x D ˚" y M x D , (58)
x

G x ˚" ´x D x , x C x u ˚" x C x v . (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.

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

y M u " y M x C b y M y D , y M v " y M x D b y M y C , y M η " y M x D b y M y D , (60)
x D u " x D x b y M y D , x D v " y M x D b x D y , (61)
x G u " x G x b y M y D , x G v " y M x D b x G y , (62)
x C u " x C x u b x C y v , x C v " x C x v b x C y u , (63)
and they satisfy

y M u ˚" y M u , y M v ˚" y M v , y M η ˚" y M η , (64)
x

G u ˚" ´x D u , x G v ˚" ´x D v , x C u ˚" x C v . (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:

x M " ˜y M x C 0 0 y M x D ¸, p S " ˜0 f x C x u ´f x C x v 0 ¸. ( 66
)
plus a identical system using y. Inertia waves do not have any dependence on η, and therefore x " rû 1 , . . . , ûd , v1 , . . . , vd s T .

Proposition 2. The eigenvalues of (66) are the solution of the reduced eigenvalue problem

ω 2 1D,x û " f 2 p y M x C q ´1 x C x u p y M x D q ´1 x C x u û " f 2 z F IW,x u û (67) and ω 2 1D,y v " f 2 p y M y D q ´1 x C y v p y M y C q ´1 x C y u v " f 2 z F IW,y v v (68)
Proof. This is obtained by noting that (66) is solving for the determinant of

detp p F `iω Iq " detpp x Mq ´1 p S `iω Iq (69) Now use p F `iω I " « iω I f p y M x C q ´1 x C x u ´f p y M x D q ´1 x C x v iω I ff ( 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 ω ‰ 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 û in x and v in y for the 1D problem. This corresponds with the choice to solve for û in 2D. The same exact expression for the 2D eigenvalues in terms of the 1D eigenvalues is obtained if we solve instead for v in x and û in y for the 1D problem and v in 2D.

2D IW

Letting g " H " 0 in the 2D IGW problem [START_REF] Thuburn | Numerical wave propagation on the hexagonal c-grid[END_REF] gives

x M " ˜y M u 0 0 y M v ¸, p S " ˜0 f x C u ´f x C v 0 ¸. (71)
This has the same form as (66).

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

ω 2 2D û " f 2 p y M u q ´1 x C u p y M v q ´1 x C v û " f 2 y F IW u û (72)
Proof. This again follows by expressing F `iω I as a block matrix

F `iω I " " iω I f pM u q ´1 C u ´f pM v q ´1 C v iω I  ( 73
)
and using the lower right block for the 2 x 2 block matrix determinant.

Proposition 4. The matrix y F IW u can be written as the Kronecker product of the 1D matrices z F IW,x u and z

F IW,y v y F IW u " z F IW,x u b z F IW,y v (74)
Therefore, the eigenvalues of the 2D inertia wave problem (72) can be written as

w 2 2D " w 2 1D,x w 2 1D,y f 2 (75)
Proof. Start by expanding the right hand side y F IW u of (72) in terms of the 1D matrices using the tensor product structure (60) -(63) to get

y F IW u " p y M u q ´1 x C u p y M v q ´1 x C v " p y M x C b y M y D q ´1p x C x u b x C y v qp y M x D b y M y C q ´1p x C x v b x C y u q (76) Applying the inverse property (50) this becomes y F IW u " pp y M x C q ´1 b p y M y D q ´1qp x C x u b x C y v qpp y M x D q ´1 b p y M y C q ´1qp x C x v b x C y u q (77)
Now using the mixed property (51) this is

y F IW u " pp y M x C q ´1 x C x u p y M x D q ´1 x C x v q b pp y M y D q ´1 x C y v p y M y C q ´1 x C y u q " z F IW,x u b z F IW,y v (78)
Therefore, using (54) we can express the eigenvalues of y F IW u in terms of those for z

F IW,x u and z F IW,y v w 2 2D f 2 " w 2 1D,x w 2 1D,y f 4 (79)
which can be simplified to give (75).

Note that (75) is the same as equation (80) (with g " H " 0) from [START_REF] Melvin | Dispersion analysis of the pnpn1dg mixed finite element pair for atmospheric modelling[END_REF].

4.5. 2D Gravity Waves 4.5.1. 1D GW In 1D, gravity waves have no dependence on v and therefore x " rû 1 , . . . , ûd , η1 , . . . , ηd s T . Letting f " 0 in the 1D IGW problem (57) gives the simplified form for 1D GW:

x M " ˜y M x C 0 0 y M x D ¸, p S " ˜0 ´g x G x ´H x D x 0 ¸. ( 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

ω 2 1D,x η " ´gHp y M x D q ´1 x D x p y M x C q ´1 x G x v " ´gH { F GW,x η η (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.

2D GW

Letting f " 0 in the 2D IGW problem [START_REF] Thuburn | Numerical wave propagation on the hexagonal c-grid[END_REF] gives the 2D GW problem

x M " ¨y M u 0 0 0 y M v 0 0 0 y M η ‹ ', S " ¨0 0 ´g x G u 0 0 ´g x G v ´H x D u ´H x D v 0 ‹ '. (82)
Proposition 6. The solutions to 2D eigenvalue problem (82) can be written as

ω 2 2D η " ´gH " p y M η q ´1 x D u p y M u q ´1 x G u `py M η q ´1 x D v p y M v q ´1 x G v ı η " ´gH z F GW η η ( 83
)
plus n geostrophic solutions ω " 0.

Proof. Start by rewriting (82) as

˜iω x M ´g p G ´H p D iω y M η ¸(84) where x M " ˜y M u 0 0 y M v ¸p G " ˜x G u x G v ¸p D " ´x D u x D v ¯. ( 85
)
Now using the 2 x 2 block matrix determinant formula the determinant of (84) is

detpiω x Mqdetpiω y M η ´gH iω p D x M ´1 p Gq. (86)
The first term in (86) gives n geostrophic modes ω " 0 along with a factor of piωq n that can be absorbed into the second term, yielding

detpω 2 y M η `gH p D x M ´1 p Gq, (87)
which is equivalent to detpω 2 `gH y

M η ´1 p D x M ´1 p Gq. (88)
Now using the definitions for p D, p G and x M we have

y M η ´1 p D x M ´1 p G " y M η ´1 ´x D u x D v ¯˜y M u 0 0 y M v ¸´1 ˜x G u x G v ¸ (89)
Expansion using block matrix multiplication then gives (83).

Proposition 7. The 2D matrix z F GW η can be written as the Kronecker sum of the 1D matrices

{ F GW,x η and { F GW,y η z F GW η " { F GW,x η b I `I b { F GW,y η " { F GW,x η ' { F GW,y η (90)
Therefore, the eigenvalues of the 2D gravity wave problem (82) can be written as

w 2 2D " w 2 1D,x `w2 1D,y (91)
Proof. Start by expanding the first term of (83) in terms of the 1D matrices using as (60) -( 63)

p y M η q ´1 x D u p y M u q ´1 x G u " p y M x D b y M y D q ´1p x D x b y M y D qp y M x C b y M y D q ´1p x G x b y M y D q (92)
Use the inverse property (50) to get

pp y M x D q ´1 b p y M y D q ´1qp x D x b y M y D qpp y M x C q ´1 b p y M y D q ´1qp x G x b y M y D q (93)
Now the mixed product property (51) gives

pp y M x D q ´1 x D x p y M x C q ´1 x G x q b pp y M y D q ´1 y M y D p y M y D q ´1 y M y D q " { F GW,x η b I (94)
Identical calculations on the p y

M η q ´1 x D v p y M v q ´1 x G v term yield p y M η q ´1 x D v p y M v q ´1 x G v " I b { F GW,y η (95)
Therefore (90) holds. Since the 2D GW matrix z F GW η is the Kronecker sum of the 1D GW matrices { F GW,x η and { F 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.

Note that (91) is the same as equation (80) (with f " 0) from [START_REF] Melvin | Dispersion analysis of the pnpn1dg mixed finite element pair for atmospheric modelling[END_REF].

2D Inertia-Gravity Waves

In [START_REF] Melvin | Dispersion analysis of the pnpn1dg mixed finite element pair for atmospheric modelling[END_REF] 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 ě 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 ě 2 the 2D IGW eigenvalues cannot be expressed in terms of the 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 [START_REF] Melvin | Dispersion analysis of the pnpn1dg mixed finite element pair for atmospheric modelling[END_REF]. The second is a numerical comparison of the dispersion relationship computed using the approach in this paper and using the approach in [START_REF] Melvin | Dispersion analysis of the pnpn1dg mixed finite element pair for atmospheric modelling[END_REF] , 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 ě 2. Proposition 8. When d " 1, the 2D IGW eigenvalues can be be expressed in terms of the 2D GW and IW eigenvalues (and therefore the 1D GW and IW eigenvalues) as

w 2 IGW " w 2 GW `w2 IW " w 2 1D,IW,x w 2 1D,IW,y f 2 `w2 GW,1D,x `w2 GW,1D,y (96)
Proof. When d " 1, both x M and p S are 3 x 3 matrices, and there exists a direct formula for the eigenvalues of the matrix. The eigenvalue problem [START_REF] Thuburn | Numerical wave propagation on the hexagonal c-grid[END_REF] can be written as

¨iω x M u f x C u ´g x G u ´f x C v iω x M v ´g x G v ´H x D u ´H x D v iω x M η ‹ ', (97)
which has determinant (after factoring out the ω " 0 root)

ω 2 IGW " f 2 x C u x C v x M u x M v `gH ˜x G v x D v x M v x M η `x G u x D u x M u x M η ¸. ( 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

ω 2 GW " gH ˜x G v x D v x M v x M η `x G u x D u x M u x M η ¸, (99)
ω 2 IW " f 2 x C u x C v x M u x M v . ( 100
)
This proves the first part of (96), that w `w2 GW,1D,y . Therefore the second part of (96) holds.

Note that (96) is just equation (80) from [START_REF] Melvin | Dispersion analysis of the pnpn1dg mixed finite element pair for atmospheric modelling[END_REF] with n " 1. Similar arguments can be used to show that the 1D dispersion relationship (when d " 1) can be written as

w 2 1D,IGW " w 2 1D,GW `w2 1D,IW , (101)
in both x and y directions. However, comparing (96) and (101), it is not possible to express w 2 IGW purely in terms of w 2 1D,IGW , even with the ability to express it in terms of w 2 1D,IW and w 2 1D,GW . Since this is not possible when d " 1, it seems very likely that it is also not possible for d ě 2.

Q ń Λ k family

The Q ń Λ k family (also as finite element exterior calculus or compatible finite elements) is constructed by choosing A " P C n and B " P DG n´1 , and is defined for all n ě 1. A diagram of the Fourier degrees of freedom for n " 2 is found in Figure 4. Recall that W 0 Ă H 1 , W 1 Ă Hpdivq and W 2 Ă 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.

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 Ñ 0 we obtain the asymptotic results for n ď 2 :

ω " ω AN ˘1 2 2n`1 n ź j"1 p4j 2 ´1q
p´f 2 pk 2n `l2n q `g H pk 2n`2 `l2n`2 q a f 2 `g H pk 2 `l2 q phq 2n `Oph 2n`2 q.

It is conjectured that this result holds @n ě 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 C n ´P DG n´1 pair from [START_REF] Eldred | Dispersion analysis of compatible galerkin schemes for the 1d shallow water model[END_REF] 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 Oph 2n q order of convergence is superconvergence specific to a regular mesh, and on arbitrary smooth meshes the expected order is Oph n`1 q.

Dispersion Relationships

Plots of the dispersion relationship σ " ω f versus ( Ă kh, r lh) for n " 1, . . . , 4 are found in Figures 5 and6 with λ{ h " 2.0 (a well-resolved Rossby radius) and in Figures 7 and8 with λ{ h " 0.1 (a poorly resolved Rossby radius). The use of λ{ h ensures the same effective resolution of the Rossby radius for all n. The long wavelength part of the spectrum is wellrepresented in all cases, but when λ{ 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 ě 2, spectral gaps where the dispersion relationship is double-valued (in the form of line discontinuities) are found at Ă kh " jπ n (j " 1, . . . , n ´1) and r lh " mπ n (m " 1, . . . , n ´1). As n increases, the jump in frequency at the discontinuities for large Ă kh and r lh 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. [START_REF] Melvin | Dispersion analysis of the pnpn1dg mixed finite element pair for atmospheric modelling[END_REF] In [START_REF] Melvin | Dispersion analysis of the pnpn1dg mixed finite element pair for atmospheric modelling[END_REF], it is claimed (see equation ( 80) in that paper) that

Comparison with the results obtained in

w 2 IGW " w 2 GW `w2 IW " w 2 1D,IW,x w 2 1D,IW,y f 2 `w2 GW,1D,x `w2 GW,1D,y (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 ě 2. The first piece of evidence is found by computing the difference pw tp ´wdirect q{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, . . . , 4. This is shown in Figure 9, which gives the results for λ{ h " 2.0 and n " 1, . . . , 4 using exact quadrature. Similar results (not shown) are found for λ{ h " 0.1. As expected, when n " 1 there is no difference. However, when n ě 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

ω tp " ω AN ˘1 1440 gHpk 6 `l6 q a f 2 `gHpk 2 `l2 q h 4 `Oph 6 q (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, . . . , 4 and the asymptotics for n " 2, we conjecture that (102) holds for 2D IGW only when d " 1.

Mass Lumping for n " 2

As shown in Section 5.2, when n ě 2 the 2D IGW dispersion relationship has line discontinuities. This was also shown for the 2D GW dispersion relationship in [START_REF] Melvin | Dispersion analysis of the pnpn1dg mixed finite element pair for atmospheric modelling[END_REF][START_REF] Melvin | A two-dimensional mixed finite-element pair on rectangles[END_REF]. However, it is possible to eliminate these discontinuities for n " 2 through a careful mass lumping of the M x C and M y C matrices. This was demonstrated for the 1D IGW dispersion relationship in [START_REF] Melvin | A two-dimensional mixed finite-element pair on rectangles[END_REF][START_REF] Staniforth | Analysis of a mixed finite-element pair proposed for an atmospheric dynamical core[END_REF] and for the 2D GW and IW dispersion relationships in [START_REF] Melvin | A two-dimensional mixed finite-element pair on rectangles[END_REF]. 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 ML a to the M x C and M y C matrices that are used in constructing M u and M v , where

ML a " ¨α β ´α β γ β ´α β α ', (104)
with β " γ " 0 and α " 1 60 . The precise value of α will depend on how the basis functions are normalized, and this value corresponds to using a reference element of size r0, 1s 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 x C (M y C ) remains positive-definite, and that the spectral gap at Ă kh " π 2 ( r lh " π 2 ) for 1D IGW is eliminated. As shown in Figure 10, simply using this lumping for M x C and M y C when they are used to construct M u and M v is sufficient to ensure that the 2D IGW gaps are also eliminated. Using this mass lumping and exact quadrature, the asymptotics are ω " ω AN ¯α 4 pk 2 `l2 q a gHpk 2 `l2 q `f 2 h 2 `Oph 4 q.

As in 1D, the convergence rate is reduced from Oph 4 q to Oph 2 q, and when l " 0 the result for the lumped P C

2

´P DG 1 pair from [START_REF] Eldred | Dispersion analysis of compatible galerkin schemes for the 1d shallow water model[END_REF] is recovered.

M GD n family

The M GD n family is constructed by choosing A " GD n [START_REF] Banks | On galerkin difference methods[END_REF] 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 [START_REF] Hiemstra | High order geometric methods with exact conservation properties[END_REF]. More details about this construction can be found in [START_REF] Eldred | A Quasi-Hamiltonian Discretization of the Thermal Shallow Water Equations[END_REF]. 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 M GD n family to avoid spectral gaps.

Exact Quadrature Asymptotics

Just as for the Q ń Λ 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 Ñ 0 we obtain the asymptotic results for n ď 17 :

ω " ω AN ˘αr β r
p´f 2 pk 2n `l2n q `g H pk 2n`2 `l2n`2 q a f 2 `g H pk 2 `l2 q phq 2n `Oph 2n`2 q.

where α e and β e are given in Table 1. It is conjectured that this result holds @n ě 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 [START_REF] Eldred | Dispersion analysis of compatible galerkin schemes for the 1d shallow water model[END_REF] are recovered, albeit with a significant improvement in n from n " 9 to n " 17. Again, the Oph 2n q order of convergence is superconvergence specific to a regular mesh, and on arbitrary smooth meshes the expected order is Oph n`1 q.

2 25 ˆ313 ˆ55 ˆ72 ˆ11 2 ˆ13 2 ˆ17 ˆ19 ˆ23 n " 15 11304410571523 ˆ88334167 2 28 ˆ315 ˆ57 ˆ72 ˆ11 3 ˆ13 2 ˆ19 ˆ23 ˆ29 ˆ31 n " 17 6703742237 ˆ155719 ˆ3541 ˆ2297 2 34 ˆ315 ˆ58 ˆ75 ˆ11 3 ˆ13 2 ˆ17 2 ˆ23 ˆ29 ˆ31

Reduced Quadrature Asymptotics

It is useful to also consider reduced quadrature to ameliorate the computational cost of the M GD n elements. This becomes particularly important as the dimension increases. As proven in [START_REF] Eldred | Dispersion analysis of compatible galerkin schemes for the 1d shallow water model[END_REF], 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 Ñ 0 we obtain the asymptotic results for n ď 19:

ω " ω AN ˘αr β r g H pk 2n `l2n q a f 2 `g H pk 2 `l2 q phq n`1 `Oph n`3 q, ( 105
)
where α r and β r are given in Table 2. It is conjectured that (105) holds for n ě 1.

As in the case of exact quadrature, the results for the GD n ´DGD n´1 pair from [START_REF] Eldred | Dispersion analysis of compatible galerkin schemes for the 1d shallow water model[END_REF] 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. 15 and16 with λ{ h " 0.1. The use of λ{ h ensures the same effective resolution of the Rossby radius for all n. When λ{ h " 2.0, the spectrum is well represented, with the exception of the CD mode that occurs at Ă kh " r lh " π (this is a common feature of all known compatible Galerkin discretizations). When λ{ 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 ń Λ 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 ń Λ k (Figures 7 and8) and M GD n families (Figures 15 and16) for the case λ{ h " 0.1, the long wavelength part of the dispersion relationship for the Q ń rΛ k family is somewhat better, especially near the end of the spectrum before the underestimated branches. However, the Q ń Λ k family has a significant spectral gap in the shortest wavelength part of the spectrum that the M GD 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 [START_REF] Konor | Impacts of the horizontal and vertical grids on the numerical solutions of the dynamical equations -part 1: Nonhydrostatic inertia-gravity modes[END_REF][START_REF] Konor | Impacts of the horizontal and vertical grids on the numerical solutions of the dynamical equations -part 2: Quasi-geostrophic rossby modes[END_REF][START_REF] Nikovi | Geostrophic adjustment on hexagonal grids[END_REF][START_REF] Randall | Geostrophic adjustment and the finite-difference shallow-water equations[END_REF][START_REF] Thuburn | Numerical wave propagation on the hexagonal c-grid[END_REF].

Simulations

To investigate the issues that arise due to the presence of line discontinuities in the dispersion relationship, and to compare the Q ń Λ k and M GD n families, two different simulations were performed. Before proceeding to these tests, let us take a slight digression into a discussion of effective resolution.

Effective Resolution

A measure of effective resolution is obtained by computing the fractional error ω N ´ωAN ω AN in the dispersion relationship as a function of non-dimensional wavelengths ζx and ζy . The calculation was performed for λ{ h " 2 using exact quadrature with n " 1, . . . , 4, in the case of the Q ń Λ k family (Figure 17) and n " 1, 3, 5, 7, in the case of the M GD 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 h in both x and y, with marginal gains after that. However, for the M GD 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 λ{ h " 0.1 (not shown) and when reduced quadrature is used for the M GD n family (not shown). In fact, reduced quadrature for M GD n gives a slight increase in effective resolution. These results are similar to others found in the literature [START_REF] Eldred | Dispersion analysis of compatible galerkin schemes for the 1d shallow water model[END_REF][START_REF] Kent | Determining the effective resolution of advection schemes. part ii: Numerical testing[END_REF][START_REF] Kent | Determining the effective resolution of advection schemes. part i: Dispersion analysis[END_REF][START_REF] Ullrich | Understanding the treatment of waves in atmospheric models. part 1: The shortest resolved waves of the 1d linearized shallow-water equations[END_REF][START_REF] Ullrich | Impact and importance of hyperdiffusion on the spectral element method: A linear dispersion analysis[END_REF] for a variety of finite-difference, finite-element and finite-volume schemes.

Test Cases

Two separate test cases were used to investigate the influence of the line discontinuities and compare the Q ń Λ k and M GD n families. For both test cases, the domain Ω " r0, Ls 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, η 0 " 1 m, α " 1{60, x 0 " L{2, g " 1 ms ´1, H " 1 m and f " 25 s ´1. This gives λ " 1{25, and λ{ h " 2.4, corresponding to a well-resolved Rossby radius. These parameters fit with the tests run in [START_REF] Eldred | Dispersion analysis of compatible galerkin schemes for the 1d shallow water model[END_REF][START_REF] Melvin | Dispersion analysis of the spectral element method[END_REF][START_REF] Staniforth | Analysis of a mixed finite-element pair proposed for an atmospheric dynamical core[END_REF]. The time integrator is implicit midpoint, with ∆t " 1{120 s and taking 100 steps for a total time of T " 100{120 s.

Unsupported Test Case

The first test case sets η " η 0 exp " ´px ´x0 q 2 `py ´y0 q 2 α 2  (106) and u " v " 0. This test was run for n " 1, . . . , 6, with exact quadrature for the Q ń Λ k family and for n " 1, 3, 5, 7, with exact and reduced quadrature for the M GD n family. Additionally, a high resolution run using exact quadrature and M GD 7 with N " 720 (giving a total of 720 2 degrees of freedom) and ∆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 ń Λ k family are found in Figure 20, and those for the M GD 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 inertiagravity mode, which both the Q ń Λ k and M GD n families were able to capture. However, for n ě 2 the Q ń Λ 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 M GD 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 M GD 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 ń Λ 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.

Supported Test Case

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

η " η 0 exp " ´px ´x0 q 2 `py ´y0 q 2 α 2  cospkpx ´x0 qq cosplpy ´y0 qq (107)
where k " l " jπ{pn hq. The test was run for n " 5 and j " 1, 2, 3, 4 with exact quadrature for the Q r ´Λk family and for n " 5 with exact and reduced quadrature for the M GD n family. Results for the Q ń Λ k family are found in Figure 22, and those for the M GD n family in Figure 23. Those for reduced quadrature with the M GD n family are very similar to Figure 23 (not shown). Comparing the two families, it is clear that the Q ń Λ k family has significantly more noise and anomalous dispersion than the M GD 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 ń Λ 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.

Conclusions

This work has shown that the 2D extension of the P C n ´P DG n´1 pair, the Q ń Λ k family, shares its problems: spectral gaps and increasing maximal with n, all of which grow worse with increasing n. When n " 2, the spectral gaps can be removed via mass lumping [START_REF] Melvin | A two-dimensional mixed finite-element pair on rectangles[END_REF][START_REF] Staniforth | Analysis of a mixed finite-element pair proposed for an atmospheric dynamical core[END_REF], 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 M GD n family, also avoids the issues. In particular, for any n the M GD 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 ń Λ 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 ń Λ k family at n " 2. Simulation results support these conclusions, with the M GD n family producing significantly better results than the Q ń Λ k family when n ě 2. Therefore, it appears that the M GD 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 M GD 3 family can be found in [START_REF] Eldred | A Quasi-Hamiltonian Discretization of the Thermal Shallow Water Equations[END_REF] 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. The difference pw tp ´wdirect q{w direct ˚100 between w direct found following Section 3.6, and w tp obtained using (102) for λ{ h " 2.0 and n " 1, . . . , 4 using exact quadrature. As expected, there is no difference when n " 1, but there is an error of around 1% for n ě 2. This is evidence that (102) is not giving the correct dispersion relationship for n ě 2. The initial condition (left) ηp0q and the high resolution reference solution (right) at the end of the simulation ηpT q for the unsupported test case. The reference solution was computed using M GD 7 with N " 720 (720 2 degrees of freedom total) and ∆T " 1{720 s (600 steps total). j " 1 j " 2 j " 3 j " 4

Figure 22: Height ηpT q for the supported test case using the Q ń Λ k family using n " 5 and exact quadrature. This test forces specifically at the gap frequency k " l " jπ{pn hq, 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. j " 1 j " 2 j " 3 j " 4

Figure 2 :

 2 Figure 2: Plots of the analytic dispersion relationship σ AN versus non-dimensional wavenumbers Ăkh and r lh for λ{ h " 2.0 (left) and λ{ h " 0.1 (right). Note that σ AN is monotonic and non-zero for all Ă kh, r lh ą 0.

Figure 3 :

 3 Figure 3: The 2D discrete deRham complex formed by taking tensor products of A and B. The operators ∇ T and ∇¨are strong, while r ∇ and r ∇ T ¨are weak (defined by integration by parts). The degrees of freedom are illustrated for d " 1.

  49) • Inverse: If A and B are invertible then pA b Bq ´1 " A ´1 b B ´1 (50) • Mixed Product: If A, B, C, D are sized such that the products A C and B D can be formed then pA b BqpC b Dq " A C b B D (51) Also recall the Kronecker sum ' for an n x n matrix A and m x m matrix B A ' B " A b I m `In b B (52

Figure 4 :

 4 Figure 4: The Fourier degrees of freedom for the Q ń Λ 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 (ψ, 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 (η) degrees of freedom on the right.

Figure 5 :

 5 Figure 5: Contour plots of the numerical dispersion relationship σ " ω f versus ( Ă kh, r lh) for n " 1, . . . , 4using the Q ń Λ k family and exact quadrature when the Rossby radius is well resolved (λ{ h " 2.0). Note the presence of spectral gaps (in the form of line discontinuities) at Ă kh " jπ n (j " 1, . . . , n ´1) and r lh " mπ n

Figure 9 :

 9 Figure9: The difference pw tp ´wdirect q{w direct ˚100 between w direct found following Section 3.6, and w tp obtained using (102) for λ{ h " 2.0 and n " 1, . . . , 4 using exact quadrature. As expected, there is no difference when n " 1, but there is an error of around 1% for n ě 2. This is evidence that (102) is not giving the correct dispersion relationship for n ě 2.

Figure 12 :Figure 13 :Figure 19 :

 121319 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 M GD 3 family. Unlike the Q ń Λ 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.

Figure 23 :

 23 Figure 23: As Figure 22, but for the M GD n family with n " 5. Unlike the Q ń Λ k family, there is no little anomalous dispersion and almost no noise for any j.

Table 1 :

 1 The coefficients α e and β e for the M GD n family using exact quadrature.

	n	α e

Table 2 :

 2 The coefficients α r and β r for the M GD n family using reduced quadrature.

	n	α r	β r
	n " 3	´1	2 3 ˆ33 ˆ5
	n " 5	´5	2 5 ˆ34 ˆ7
	n " 7	´17	2 8 ˆ33 ˆ5	ˆ7
	n " 9	´133741	2 10 ˆ39 ˆ5 ˆ7 ˆ11
	n " 11	´23 ˆ1069 ˆ4111	2 14 ˆ310 ˆ52 ˆ7 ˆ11 ˆ13
	n " 13	´44346011	2 16 ˆ310 ˆ53 ˆ7 ˆ13
	n " 15	´824695933877	2 19 ˆ314 ˆ53 ˆ72 ˆ13 ˆ17
	n " 17	´45381741803663	2 21 ˆ315 ˆ53 ˆ72 ˆ13 ˆ17 ˆ19
	n " 19 ´3837131 ˆ37661489	2 26 ˆ315 ˆ53 ˆ73 ˆ13 ˆ19
	6.3. Dispersion Relationships
	Plots of the dispersion relationship σ " ω f versus p Ă kh, r lhq for n " 1, 3, 5, 7, are found in
	Figures 13 and 14 with λ{ h " 2.0 and in Figures

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

 5, but for a poorly resolved Rossby radius (λ{ 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.

n " 1 k h 0 13, but with a poorly resolved Rossby radius (λ{ 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.

n " 1 k h 0 13, but with a poorly resolved Rossby radius (λ{ 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. 17, but using the M GD n family with n " 1, . . . , 7. Again, the largest gain in effective resolution comes from increasing n from 1 to 3, giving around 4 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. The final height ηpT q for the unsupported test case using the Q ń Λ k family with exact quadrature for n " 1, . . . , 6. There is significant noise and anomalous dispersion for all n ě 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.

n " 1 n " 3 n " 5 n " 7

Figure 21: As for Figure 20, but using the M GD n family with n " 1, 3, 5, 7. Unlike the Q ń Λ k family there is essentially no noise for any n and reduced anomalous dispersion, especially for higher n.