# Multidimensional method-of-lines transport for atmospheric flows over steep
terrain using arbitrary meshes

James Shaw<EMAIL_ADDRESS>Hilary Weller  John Methven  Terry Davies
Department of Meteorology, University of Reading, Reading, United Kingdom  Met
Office, Exeter, United Kingdom

###### Abstract

Including terrain in atmospheric models gives rise to mesh distortions near
the lower boundary that can degrade accuracy and challenge the stability of
transport schemes. Multidimensional transport schemes avoid splitting errors
on distorted, arbitrary meshes, and method-of-lines schemes have low
computational cost because they perform reconstructions at fixed points.

This paper presents a multidimensional method-of-lines finite volume transport
scheme, “cubicFit”, which is designed to be numerically stable on arbitrary
meshes. Constraints derived from a von Neumann stability analysis are imposed
during model initialisation to obtain stability and improve accuracy in
distorted regions of the mesh, and near steeply-sloping lower boundaries.
Reconstruction calculations depend upon the mesh only, needing just one vector
multiply per face per time-stage irrespective of the velocity field.

The cubicFit scheme is evaluated using three, idealised numerical tests. The
first is a variant of a standard horizontal transport test on severely
distorted terrain-following meshes. The second is a new test case that
assesses accuracy near the ground by transporting a tracer at the lower
boundary over steep terrain on terrain-following meshes, cut-cell meshes, and
new, slanted-cell meshes that do not suffer from severe time-step constraints
associated with cut cells. The third, standard test deforms a tracer in a
vortical flow on hexagonal-icosahedral meshes and cubed-sphere meshes. In all
tests, cubicFit is stable and largely insensitive to mesh distortions, and
cubicFit results are more accurate than those obtained using a
multidimensional linear upwind transport scheme. The cubicFit scheme is
second-order convergent regardless of mesh distortions.

###### keywords:

Finite volume, unstructured mesh, atmospheric modelling, least-squares
approximation, von Neumann stability analysis

##  1 Introduction

Numerical simulations of atmospheric flows solve equations of motion that
result in the transport of momentum, temperature, water species and trace
gases. The numerical representation of Earth’s terrain complicates the
transport problem because the mesh is necessarily distorted next to the lower
boundary. As new atmospheric models use increasingly fine mesh spacing, meshes
are able to resolve steep, small-scale slopes. Numerical schemes in
operational weather forecast models can perform poorly over large mountain
ranges, exhibiting small-scale numerical noise in momentum [1], temperature,
humidity [2] and potential vorticity fields [3], or even violating the
Courant–Friedrich–Lewy stability constraint resulting in so-called ‘grid-point
storms’ [4]. A transport scheme is desired that yields stable and accurate
solutions, particularly near the surface where the weather affects us
directly. We present a new transport scheme which is numerically stable on
arbitrary meshes and which is generally insensitive to mesh distortions
created by steep slopes. It has a low computational cost since most
calculations are not repeated every time-step because they depend upon the
mesh geometry only.

There are two main methods for representing terrain in atmospheric models:
terrain-following layers and cut cells. Both methods modify regular meshes to
produce distorted meshes with cells that are aligned in columns. Most
operational models use terrain-following layers in which horizontal mesh
surfaces are moved upwards to accommodate the terrain. A vertical decay
function is chosen so that mesh surfaces slope less steeply with increasing
height. The most straightforward is the linear decay function used by the
basic terrain-following transform [5] (also called the sigma coordinate), but
many atmospheric models suffer from large numerical errors on such meshes [2,
6, 7]. To reduce such errors, more complex decay functions have been developed
so that mesh surfaces are smoother [8, 2, 9, 6].

An alternative to terrain-following layers is the cut cell method. Cut cell
meshes are constructed by ‘cutting’ a regular mesh with a piecewise-linear
representation of the terrain. New vertices are created where the terrain
intersects mesh edges, and cell volumes that lie beneath the ground are
removed. Cut cell meshes can have arbitrarily small cells that impose severe
time-step constraints on explicit transport schemes. Several techniques have
been developed to alleviate this problem, known as the ‘small-cell problem’:
small cells can be merged with adjacent cells [10], cell volumes can be
artificially increased [11], or an implicit scheme can be used near the ground
with an explicit scheme used aloft [12].

Another method for avoiding the small-cell problem was proposed by Shaw and
Weller [13] in which cell vertices are moved vertically so that they are
positioned at the terrain surface. We refer to this alternative method as the
slanted cell method in order to distinguish it from the traditional cut cell
method. Slanted cell meshes do not suffer from arbitrarily small cells because
the horizontal cell dimensions are not modified by the presence of terrain.

Smoothed terrain-following layers, cut cells and slanted cell methods all
reduce the amount of mesh distortion but any mesh that represents sloping
terrain must necessarily be distorted, at least near the ground. Even when
distortions are minimal, transport across mesh surfaces tends to be more
common near steep slopes, and this misalignment between the flow and mesh
surfaces increases numerical errors [14, 2, 13]. A huge variety of transport
schemes have been developed for atmospheric models, but few are able to
account for distortions associated with steep terrain because they treat
horizontal and vertical transport separately [15], resulting in numerical
errors called ‘splitting errors’. Such errors can be reduced by explicitly
accounting for transverse fluxes when combining fluxes [16], but splitting
errors are still apparent in flows over steep terrain where meshes are highly
distorted and metric terms in a terrain-following coordinate transform are
large [17].

Transport schemes are often classified as dimensionally-split or
multidimensional. Dimensionally-split schemes such as [18, katta2015]
calculate transport in each dimension separately before the flux contributions
are combined. Such schemes are computationally efficient and allow existing
one-dimensional high-order methods to be used. Dimensionally-split schemes
have only been used with quadrilateral meshes where dimensions are inherently
separable. Special treatment is required at the corners of cubed-sphere panels
where local coordinates differ [20, katta2015]. For similar reasons,
dimensionally-split schemes have only been used with terrain-following
coordinate transforms and not cut cells. Perhaps confusingly, dimensionally-
split schemes are sometimes called multidimensional, too, because they use
one-dimensional techniques for multidimensional transport.

Unlike dimensionally-split schemes, multidimensional schemes consider
transport in two or three dimensions together. There are several subclasses of
multidimensional schemes that include semi-Lagrangian finite volume schemes
(also called conservative mesh remapping), swept-area schemes (also called
flux-form semi-Lagrangian, incremental remapping, or forward-in-time), and
method-of-lines schemes (also called Eulerian schemes). Two-dimensional semi-
Lagrangian finite volume schemes such as [21, 22] integrate over departure
cells that are found by tracing backward the trajectories of cell vertices.
These schemes are conservative because departure cells are constructed so that
there are no overlaps or gaps, which requires that cell areas are simply-
connected domains [23]. SLICE-3D is a three-dimensional semi-Lagrangian finite
volume scheme for latitude-longitude meshes that applies separate conservative
remappings in each dimension [24]. Swept area schemes such as [25, 26, 27, 28]
calculate the flux through a cell face by integrating over the upstream area
that is swept out over one time-step. Such schemes differ in their choice of
area approximation, sub-grid reconstruction, and spatial integration method.
Because swept area schemes integrate over the reconstructed field, they
typically require a matrix-vector multiply per face per time-stage [28, 26].
Method-of-lines schemes such as [29, 30] use a spatial discretisation to
reduce the transport PDE to an ODE that is typically solved using a multi-
stage time-stepping method. A method-of-lines scheme using a spectral element
reconstruction was recently developed to achieve accurate solutions near the
surface of cut cell meshes [31]. Unlike semi-Lagrangian finite volume schemes,
swept-area and method-of-lines schemes achieve conservation for small-scale
rotational flows. Such flows can twist the departure domain to such an extent
that the domain intersects itself [27]. In two dimensions, a self-intersecting
departure domain has a bowtie or hourglass shape. There are many more types of
atmospheric transport schemes, but all can be classified according to their
treatment of the three spatial dimensions. A more comprehensive overview is
presented by Lauritzen et al. [32].

For transport schemes that are ordinarily classified as ‘multidimensional’, a
further distinction ought to made between horizontally-multidimensional and
three-dimensional schemes. Most multidimensional schemes are only
horizontally-multidimensional because, while the two horizontal dimensions are
considered together, horizontal and vertical transport are still treated
separately. This separate treatment becomes less justifiable as atmospheric
models are using increasingly fine horizontal mesh spacings that resolve
small-scale steep slopes, resulting in greater mesh distortion and possible
splitting errors [15]. Three-dimensional schemes avoid any splitting errors
over steep slopes, but only a few conservative three-dimensional schemes have
been used in atmospheric models. The multi-moment constrained finite volume
scheme [33] is a three-dimensional scheme that has been used to simulate
nonhydrostatic flows over orography with terrain-following coordinates on a
$x$–$z$ plane [34]. Simulations of subcritical flow around a cylinder have
also been performed on a three-dimensional hexahedral-prismatic hybrid mesh
[35]. The Multidimensional Positive Definite Advection Transport Algorithm
(MPDATA) is another three-dimensional scheme that is suitable for arbitrary
meshes. It has been used on triangular unstructured meshes to simulate two-
dimensional nonhydrostatic flows over orography [36], and in three-dimensional
transport tests [37]. Most recently, MPDATA has been extended to enable semi-
implicit integrations of the compressible Euler equations on arbitrary meshes
[38]. The three-dimensional method-of-lines scheme developed by Weller and
Shahrokhi [39] has been used in two-dimensional flows over orography on
Cartesian $x$–$z$ planes with distorted meshes [13, 17]. This finite volume
scheme uses a moving least-squares reconstruction that makes it suitable for
arbitrary meshes. This least-squares approach has been applied previously to
shallow water flows [40], aeronautic [41] and porous media [42] simulations.

In this paper, we present a new multidimensional method-of-lines scheme,
‘cubicFit’, that improves the stability of the Weller and Shahrokhi scheme
[39] and avoids all splitting errors. To reconstruct values at cell faces, the
scheme fits a multidimensional cubic polynomial over an upwind-biased stencil
using a least-squares approach. The implementation uses constraints derived
from a von Neumann stability analysis to select appropriate polynomial fits
for stencils in highly-distorted mesh regions. Almost all of the least-squares
procedure depends upon the mesh geometry only and reconstruction weights can
be precomputed without knowledge of the velocity field or tracer field. Hence,
the computational cost of the cubicFit scheme is lower than most swept-area
schemes, being more comparable to dimensionally-split schemes, requiring only
$n$ multiplies per cell face per time-stage where $n$ is the size of the
stencil. Based on numerical experiments, the scheme is found to be
conditionally stable up to maximum Courant numbers of about $1.3$ to $3.3$.

The remainder of this paper is organised as follows. Section 2 starts by
discretising the transport equation using a method-of-lines approach before
describing the cubicFit transport scheme and a multidimensional linear upwind
transport scheme. Section 3 evaluates the cubicFit scheme using three
idealised numerical tests. The first test follows Schär et al. [2],
transporting a tracer horizontally above steep mountains on two-dimensional,
highly-distorted terrain-following meshes. The second is a new test case
designed to assess numerical accuracy next to a mountainous lower boundary. In
this test, a tracer placed at the ground is transported over steep slopes by a
terrain-following velocity field on terrain-following, cut cell and slanted
cell meshes. The third is a standard test of deformational flow on a single-
layer spherical Earth, specified by Lauritzen et al. [43], which we use to
assess the cubicFit transport scheme on hexagonal-icosahedral meshes and
cubed-sphere meshes. Concluding remarks are made in section 4.

##  2 Transport schemes for arbitrary meshes

The transport of a dependent variable $\phi$ in a prescribed, non-divergent
velocity field $\mathbf{u}$ is given by the equation

$\displaystyle\frac{\partial\phi}{\partial
t}+\nabla\cdot\left(\mathbf{u}\phi\right)=0\text{ .}$ (1)

The time derivative is discretised using a two-stage, second-order Heun
method,

$\displaystyle\phi^{\star}$ $\displaystyle=\phi^{(n)}+\Delta t\>g(\phi^{(n)})$
(2a) $\displaystyle\phi^{(n+1)}$ $\displaystyle=\phi^{(n)}+\frac{\Delta
t}{2}\left[g(\phi^{(n)})+g(\phi^{\star})\right]$ (2b)

where $g(\phi^{(n)})=-\nabla\cdot(\mathbf{u}\phi^{(n)})$ at time level $n$.
The same time-stepping method is used for both the cubicFit scheme and the
multidimensional linear upwind scheme. Although the Heun method is unstable
for a linear oscillator [44] and for solving the transport equation using
centred, linear differencing, it is stable when it is used for transport
schemes with sufficient upwinding.

Using the finite volume method, the velocity field is prescribed at face
centroids and the dependent variable is stored at cell centroids. The
divergence term in equation (1) is discretised using Gauss’s theorem:

$\displaystyle\nabla\cdot\left(\mathbf{u}\phi\right)\approx\frac{1}{\mathcal{V}_{c}}\sum_{f\in\>c}\mathbf{u}_{f}\cdot\mathbf{S}_{f}\phi_{F}$
(3)

where subscript $f$ denotes a value stored at a face and subscript $F$ denotes
a value approximated at a face from surrounding values. $\mathcal{V}_{c}$ is
the cell volume, $\mathbf{u}_{f}$ is a velocity vector prescribed at a face,
${\mathbf{S}_{f}}$ is the surface area vector with a direction outward normal
to the face and a magnitude equal to the face area, $\phi_{F}$ is an
approximation of the dependent variable at the face, and $\sum_{f\in\>c}$
denotes a summation over all faces $f$ bordering cell $c$. Note that equation
(3) is a second-order approximation of the divergence term which limits the
cubicFit transport scheme to second-order numerical convergence.

This discretisation is applicable to arbitrary meshes. A necessary condition
for stability is given by the multidimensional Courant number,

$\displaystyle\mathrm{Co}_{c}=\frac{\Delta
t}{2\mathcal{V}_{c}}\sum_{f\in\>c}\lvert\mathbf{u}\cdot\mathbf{S}_{f}\rvert$
(4)

such that $\mathrm{Co}_{c}\leq 1$ for all cells $c$ in the domain. Hence,
stability is constrained by the maximum Courant number of any cell in the
domain.

The accurate approximation of the dependent variable at the face, $\phi_{F}$,
is key to the overall accuracy of the transport scheme. The cubicFit scheme
and multidimensional linear upwind scheme differ in their approximations, and
these approximation methods are described next.

###  2.1 Cubic fit transport scheme

Figure 1: Upwind-biased stencils for faces far away from the boundaries of
two-dimensional (a) rectangular and (b) hexagon meshes. The stencil is used to
fit a multidimensional polynomial to cell centre values, $\phi_{c}$, marked by
grey circles, in order to approximate the value $\phi_{F}$ at the face
centroid marked by an open circle. $\phi_{u}$ and $\phi_{d}$ are the values at
the centroids of the upwind and downwind cells neighbouring the target face,
drawn with a heavy line. The velocity vector $\mathbf{u}_{f}$ is prescribed at
face $f$ and determines the choice of stencil at each time-step.

The cubicFit scheme approximates the value of the dependent variable at the
face, $\phi_{F}$, using a least-squares fit over a stencil of surrounding
known values. To introduce the approximation method, we will consider how an
approximate value is calculated for a face that is far away from the
boundaries of a two-dimensional uniform rectangular mesh. For any mesh, every
interior face connects two adjacent cells. The velocity direction at the face
determines which of the two adjacent cells is the upwind cell. Since the
stencil is upwind-biased and asymmetric, two stencils must be constructed for
every interior face, and the appropriate stencil is chosen depending on the
velocity direction at each face for every time-step.

The upwind-biased stencil for a face $f$ is shown in figure 1a. The wind at
the face, $\mathbf{u}_{f}$, is blowing from the upwind cell $c_{u}$ to the
downwind cell $c_{d}$. To obtain an approximate value at $f$, a polynomial
least-squares fit is calculated using the stencil values. The stencil has $4$
points in $x$ and $3$ points in $y$, leading to a natural choice of polynomial
that is cubic in $x$ and quadratic in $y$,

$\displaystyle\phi=a_{1}+a_{2}x+a_{3}y+a_{4}x^{2}+a_{5}xy+a_{6}y^{2}+a_{7}x^{3}+a_{8}x^{2}y+a_{9}xy^{2}\text{
.}$ (5)

A least-squares approach is needed because the system of equations is
overconstrained, with $12$ stencil values but only $9$ polynomial terms. The
stencil geometry is expressed in a local coordinate system with the face
centroid as the origin so that the approximated value $\phi_{F}$ is equal to
the constant coefficient $a_{1}$. The stencil is upwind-biased to improve
numerical stability, and the multidimensional cubic polynomial is chosen to
improve accuracy in the direction of flow [14].

The remainder of this subsection generalises the approximation technique for
arbitrary meshes and describes the methods for constructing stencils,
performing a least-squares fit with a suitable polynomial, and ensuring
numerical stability of the transport scheme.

####  2.1.1 Stencil construction

For every interior face, two stencils are constructed, one for each of the
possible upwind cells. Stencils are not constructed for boundary faces because
values of $\phi$ at boundaries are calculated from prescribed boundary
conditions. For a given interior face $f$ and upwind cell $c_{u}$, we find
those faces that are connected to $c_{u}$ and ‘oppose’ face $f$. These are
called the opposing faces. The opposing faces for face $f$ and upwind cell
$c_{u}$ are determined as follows. Defining $G$ to be the set of faces other
than $f$ that border cell $c_{u}$, we calculate the ‘opposedness’,
$\mathrm{Opp}$, between faces $f$ and $g\in G$, defined as

$\displaystyle\mathrm{Opp}(f,g)\equiv-\frac{\mathbf{S}_{f}\cdot\mathbf{S}_{g}}{|\mathbf{S}_{f}|^{2}}$
(6)

where $\mathbf{S}_{f}$ and $\mathbf{S}_{g}$ are the surface area vectors
pointing outward from cell $c_{u}$ for faces $f$ and $g$ respectively. Using
the fact that
$\mathbf{a}\cdot\mathbf{b}=|\mathbf{a}|\>|\mathbf{b}|\cos(\theta)$ we can
rewrite equation (6) as

$\displaystyle\mathrm{Opp}(f,g)=-\frac{|\mathbf{S}_{g}|}{|\mathbf{S}_{f}|}\cos(\theta)$
(7)

where $\theta$ is the angle between faces $f$ and $g$. In this form, it can be
seen that $\mathrm{Opp}$ is a measure of the relative area of $g$ and how
closely it parallels face $f$.

The set of opposing faces, $\mathrm{OF}$, is a subset of $G$, comprising those
faces with $\mathrm{Opp}\geq 0.5$, and the face with the maximum opposedness.
Expressed in set notation, this is

$\displaystyle\mathrm{OF}(f,c_{u})\equiv\\{g:\mathrm{Opp}(f,g)\geq
0.5\\}\cup\\{g:\max_{g\>\in\>G}(\mathrm{Opp}(f,g))\\}\text{ .}$ (8)

On a rectangular mesh, there is always one opposing face $g$, and it is
exactly parallel to the face $f$ such that $\mathrm{Opp}(f,g)=1$.

Once the opposing faces have been determined, the set of internal and external
cells must be found. The internal cells are those cells that are connected to
the opposing faces. Note that $c_{u}$ is always an internal cell. The external
cells are those cells that share vertices with the internal cells. Note that
$c_{d}$ is always an external cell. Finally, the stencil boundary faces are
boundary faces having Dirichlet boundary conditions111Boundary faces with
Neumann boundary conditions would require extrapolated boundary values to be
calculated. This would create a feedback loop in which boundary values are
extrapolated from interior values, then interior values are transported using
stencils that include boundary values. We have not considered how such an
extrapolation could be made consistent with the multidimensional polynomial
reconstruction. Hence, boundary faces with Neumann boundary conditions are
excluded from the set of stencil boundary faces. that share a vertex with the
internal cells. Having found these three sets, the stencil is constructed to
comprise all internal cells, external cells and stencil boundary faces.

Figure 2: A fourteen-point, upwind-biased stencil for face $f$ connecting the
pentagonal upwind cell, $c_{u}$, and the downwind cell $c_{d}$. The dashed
lines denote the two faces of cell $c_{u}$ that oppose $f$, and black circles
mark the centroids of the internal cells that are connected to these two
opposing faces. The stencil is extended outwards by including cells that share
vertices with the three internal cells, where black squares mark these
vertices. Four stencil boundary faces, marked by black triangles, are also
included. The local coordinate system $(x,y)$ has its origin at the centroid
of face $f$, marked by an open circle, with $x$ normal to $f$ and $y$
perpendicular to $x$.

Figure 2 illustrates a stencil construction for face $f$ connecting upwind
cell $c_{u}$ and downwind cell $c_{d}$. The two opposing faces are denoted by
thick dashed lines and the centres of the three adjoining internal cells are
marked by black circles. The stencil is extended outwards by including the
external cells that share vertices with the internal cells, where the vertices
are marked by black squares. A boundary at the far left has Dirichlet boundary
conditions, and so the four stencil boundary faces are also included in the
stencil, where the boundary face centres are marked by black triangles. The
resultant stencil contains fourteen points.

####  2.1.2 Least-squares fit

To approximate the value of $\phi$ at a face $f$, a least-squares fit is
calculated from a stencil of surrounding known values. First, we will show how
a polynomial least-squares fit is calculated for a face on a rectangular mesh.
Second, we will make modifications to the least-squares fit that are necessary
for numerical stability.

For faces that are far away from the boundaries of a rectangular mesh, we fit
the multidimensional polynomial given by equation (5) that has nine unknown
coefficients, $\mathbf{a}=a_{1}\ldots a_{9}$, using the twelve cell centre
values from the upwind-biased stencil, $\bm{\phi}=\phi_{1}\ldots\phi_{12}$.
This yields a matrix equation

$\displaystyle\begin{bmatrix}1&x_{1}&y_{1}&x_{1}^{2}&x_{1}y_{1}&y_{1}^{2}&x_{1}^{3}&x_{1}^{2}y_{1}&x_{1}y_{1}^{2}\\\
1&x_{2}&y_{2}&x_{2}^{2}&x_{2}y_{2}&y_{2}^{2}&x_{2}^{3}&x_{2}^{2}y_{2}&x_{2}y_{2}^{2}\\\
\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\\
1&x_{12}&y_{12}&x_{12}^{2}&x_{12}y_{12}&y_{12}^{2}&x_{12}^{3}&x_{12}^{2}y_{12}&x_{12}y_{12}^{2}\\\
\end{bmatrix}\begin{bmatrix}a_{1}\\\ a_{2}\\\ \vdots\\\
a_{9}\end{bmatrix}=\begin{bmatrix}\phi_{1}\\\ \phi_{2}\\\ \vdots\\\
\phi_{12}\end{bmatrix}$ (9)

which can be written as

$\displaystyle\mathbf{B}\mathbf{a}=\bm{\phi}\text{ .}$ (10)

The rectangular matrix $\mathbf{B}$ has one row for each cell in the stencil
and one column for each term in the polynomial. $\mathbf{B}$ is called the
stencil matrix, and it is constructed using only the mesh geometry. A local
coordinate system is established in which $x$ is normal to the face $f$ and
$y$ is perpendicular to $x$. The coordinates $(x_{i},y_{i})$ give the position
of the centroid of the $i$th cell in the stencil. A two-dimensional stencil is
also used for the tests on spherical meshes in section 3.3. In these tests,
cell centres are projected perpendicular to a tangent plane at the face
centre. Previous studies found that results were largely insensitive to the
projection method [30, 25].

The unknown coefficients $\mathbf{a}$ are calculated using the pseudo-inverse,
$\mathbf{B}^{+}$, found by singular value decomposition,

$\displaystyle\mathbf{a}=\mathbf{B}^{+}\bm{\phi}\text{ .}$ (11) Recall that
the approximate value $\phi_{F}$ is equal to the constant coefficient $a_{1}$,
which is a weighted mean of $\bm{\phi}$, $\displaystyle
a_{1}=\begin{bmatrix}b_{1,1}^{+}\\\ b_{1,2}^{+}\\\ \vdots\\\
b_{1,12}^{+}\end{bmatrix}\cdot\begin{bmatrix}\phi_{1}\\\ \phi_{2}\\\ \vdots\\\
\phi_{12}\end{bmatrix}$ (12)

where the weights $b_{1,1}^{+}\ldots b_{1,12}^{+}$ are the elements of the
first row of $\mathbf{B}^{+}$. Note that the majority of the least-squares fit
procedure depends on the mesh geometry only. An implementation may precompute
the pseudo-inverse for each stencil during model initialisation, and only the
first row needs to be stored. Since each face has two possible stencils
depending on the orientation of the velocity relative to the face, the
implementation stores two sets of weights for each face. Knowledge of the
values of $\bm{\phi}$ is only required to calculate the weighted mean given by
equation (12), which is evaluated once per face per time-stage.

In the least-squares fit presented above, all stencil values contributed
equally to the polynomial fit. It is necessary for numerical stability that
the polynomial fits the cells connected to face $f$ more closely than other
cells in the stencil, as shown by [25, 26]. To achieve this, we allow each
cell to make an unequal contribution to the least-squares fit. We assign an
integer multiplier to each cell in the stencil, $\mathbf{m}=m_{1}\ldots
m_{12}$, and multiply equation (10) to obtain

$\displaystyle\mathbf{\tilde{B}}\mathbf{a}=\mathbf{m}\cdot\bm{\phi}$ (13)

where $\mathbf{\tilde{B}}=\mathbf{M}\mathbf{B}$ and
$\mathbf{M}=\mathrm{diag}(\mathbf{m})$. The constant coefficient $a_{1}$ is
calculated from the pseudo-inverse, $\mathbf{\tilde{B}}^{+}$,

$\displaystyle a_{1}=\mathbf{\tilde{b}_{1}^{+}}\cdot\mathbf{m}\cdot\bm{\phi}$
(14)

where
$\mathbf{\tilde{b}_{1}^{+}}=\tilde{b}_{1,1}^{+}\ldots\tilde{b}_{1,12}^{+}$ are
the elements of the first row of $\mathbf{\tilde{B}}^{+}$. Again, $a_{1}$ is a
weighted mean of $\bm{\phi}$, where the weights are now
$\mathbf{\tilde{b}_{1}^{+}}\cdot\mathbf{m}$. Values for $\mathbf{m}$ are
chosen so that the cells connected to face $f$ make a greater contribution to
the least-squares fit, as discussed later in section 2.1.4.

For faces of a non-rectangular mesh, or faces that are near a boundary, the
number of stencil points and number of polynomial terms may differ: a stencil
will have one or more cells and, for two-dimensional meshes, its polynomial
will have between one and nine terms. Additionally, the polynomial cannot have
more terms than its stencil has cells because this would lead to an
underconstrained system of equations. The procedure for choosing suitable
polynomials is discussed next.

####  2.1.3 Polynomial generation

Figure 3: Upwind-biased stencils for faces near the lower boundary of a
rectangular $x$–$z$ mesh, with (a) a $3\times 2$ stencil for the face
immediately adjacent to the lower boundary, and (b) a $3\times 3$ stencil for
the face immediately adjacent to the face in (a). Each stencil belongs to the
face marked by a thick line. The local coordinate system is shown, having an
$x$ direction normal to the face a $y$ direction tangent to the face. For both
stencils, attempting a least-squares fit using the nine-term polynomial in
equation (5) would result in an underconstrained problem. There is no normal
flow at the lower boundary.

The majority of faces on a uniform two-dimensional mesh have stencils with
more than nine cells. For example, a rectangular mesh has 12 points (figure
1a), and a hexagonal mesh has 10 points (figure 1b). In both cases,
constructing a system of equations using the nine-term polynomial in equation
(5) leads to an overconstrained problem that can be solved using least-
squares. However, this is not true for faces near boundaries: stencils that
have fewer than nine cells (figure 3a) would result in an underconstrained
problem, and stencils that have exactly nine cells may lack sufficient
information to constrain high-order terms. For example, the stencil in figure
3b lacks sufficient information to fit the $x^{3}$ term. In such cases, it
becomes necessary to perform a least-squares fit using a polynomial with fewer
terms.

For every stencil, we find a set of candidate polynomials that do not result
in an underconstrained problem. In two dimensions, a candidate polynomial has
some combination of between one and nine terms from equation (5). There are
two additional constraints that a candidate polynomial must satisfy.

First, high-order terms may be included in a candidate polynomial only if the
lower-order terms are also included. More precisely, let

$\displaystyle M(x,y)={x^{i}y^{j}:i,j\geq 0\text{ and }i\leq 3\text{ and
}j\leq 2\text{ and }i+j\leq 3}$ (15)

be the set of all monomials of degree at most $3$ in $x,y$. A subset $S$ of
$M(x,y)$ is “dense” if, whenever $x^{a}y^{b}$ is in $S$, then $x^{i}y^{j}$ is
also in $S$ for all $0\leq i\leq a$, $0\leq j\leq b$. For example, the
polynomial $\phi=a_{1}+a_{2}x+a_{3}y+a_{4}xy+a_{5}x^{2}+a_{6}x^{2}y$ is a
dense subset of $M(x,y)$, but $\phi=a_{1}+a_{2}x+a_{3}y+a_{4}x^{2}y$ is not
because $x^{2}y$ can be included only if $xy$ and $x^{2}$ are also included.
In total there are 26 dense subsets of the two-dimensional polynomial in
equation (5).

Second, a candidate polynomial must have a stencil matrix $\mathbf{B}$ that is
full rank. The matrix is considered full rank if its smallest singular value
is greater than $1\text{\times}{10}^{-9}$. Using a polynomial with all nine
terms and the stencil in figure 3b results in a rank-deficient matrix and so
the nine-term polynomial is not a candidate polynomial.

The candidate polynomials are all the dense subsets of $M(x,y)$ that have a
cardinality greater than one with a stencil matrix that is full rank. The
final stage of the cubicFit transport scheme selects a candidate polynomial
and ensures that the least-squares fit is numerically stable.

####  2.1.4 Stabilisation procedure

So far, we have constructed a stencil and found a set of candidate
polynomials. Applying a least-squares fit to any of these candidate
polynomials avoids creating an underconstrained problem. The final stage of
the transport scheme chooses a suitable candidate polynomial and appropriate
multipliers $\mathbf{m}$ so that the fit is numerically stable.

The approximated value $\phi_{F}$ is equal to $a_{1}$ which is calculated from
equation (14). The value of $a_{1}$ is a weighted mean of $\bm{\phi}$ where
$\mathbf{w}=\mathbf{\tilde{b}_{1}^{+}}\cdot\mathbf{m}$ are the weights. If the
cell centre values $\bm{\phi}$ are assumed to approximate a smooth field then
we expect $\phi_{F}$ to be close to the values of $\phi_{u}$ and $\phi_{d}$,
and expect $\phi_{F}$ to be insensitive to small changes in $\bm{\phi}$. When
the weights $\mathbf{w}$ have large magnitude then this is no longer true:
$\phi_{F}$ becomes sensitive to small changes in $\bm{\phi}$ which can result
in large, numerically unstable departures from the smooth field $\bm{\phi}$.

To avoid numerical instabilities, a simplified, one-dimensional von Neumann
analysis was performed, presented in appendix A. The analysis is used to
impose three stability conditions on the weights $\mathbf{w}$,

$\displaystyle 0.5\leq w_{u}\leq 1$ (16a) $\displaystyle 0\leq w_{d}\leq 0.5$
(16b) $\displaystyle w_{u}-w_{d}\geq\max_{p\>\in\>P}(|w_{p}|)$ (16c)

where $w_{u}$ and $w_{d}$ are the weights for the upwind and downwind cells
respectively. The peripheral points $P$ are the cells in the stencil that are
not the upwind or downwind cells, and $w_{p}$ is the weight for a given
peripheral point $p$. The upwind, downwind and peripheral weights sum to one
such that $w_{u}+w_{d}+\sum_{p\in P}w_{p}=1$.

The stabilisation procedure comprises three steps. In the first step, the set
of candidate polynomials is sorted in preference order so that candidates with
more terms are preferred over those with fewer terms. If there are multiple
candidates with the same number of terms, the minimum singular value of
$\mathbf{B}$ is calculated for each candidate, and an ordering is imposed such
that the candidate with the larger minimum singular value is preferred. This
ordering ensures that the preferred candidate is the highest-order polynomial
with the most information content.222Note that singular values are used for
two purposes: first, to test if the matrix $\mathbf{B}$ is full-rank and,
second, to impose an ordering on candidates. We have used the minimum singular
value, $\sigma_{\mathrm{min}}(\mathbf{B})$, for both purposes. Alternatively,
we could use the condition number, $\mathrm{cond}(\mathbf{B})$, which is the
ratio of smallest to largest singular value. Experiments revealed that only
the candidate ordering was sensitive to the choice of $\sigma_{\mathrm{min}}$
or $\mathrm{cond}$. The most suitable choices of singular value calculations
could be explored in future.

In the second step, the most-preferred polynomial is taken from the list of
candidates and the multipliers are assigned so that the upwind cell and
downwind cell have multipliers $m_{u}=2^{10}$ and $m_{d}=2^{10}$ respectively,
and all peripheral points have multipliers $m_{p}=1$. These multipliers are
very similar to those used by [25], leading to a well-conditioned matrix
$\mathbf{\tilde{B}}$ and a least-squares fit in which the polynomial passes
almost exactly through the upwind and downwind cell centre values.

In the third step, we calculate the weights $\mathbf{w}$ and evaluate them
against the stability conditions given in equation (16). If any condition is
violated, the value of $m_{d}$ is halved and the conditions are evaluated with
the new weights. This step is repeated until the weights satisfy the stability
conditions, or $m_{d}$ becomes smaller than one. In practice, the conditions
are satisified when $m_{d}$ is either small (between 1 and 4) or equal to
$2^{10}$. The upwind multiplier $m_{u}$ is fixed at $2^{10}$ and the
peripheral multipliers $m_{p}$ are fixed at $1$. If the conditions are still
not satisfied, then we start again from the second step with the next
polynomial in the candidate list.

Finally, if no stable weights are found for any candidate polynomial, we
revert to an upwind scheme such that $w_{u}=1$ and all other weights are zero.
In our experiments we have not encountered any stencil for which this last
resort is required. Furthermore, our experiments show that the stabilisation
procedure only modifies the least squares fit for stencils near boundaries and
for stencils in distorted mesh regions. For stencils in the interior of a
uniform rectangular mesh, the least squares fit includes all terms in equation
(5) with $m_{u}=m_{d}=2^{10}$.

Figure 4: One-dimensional least-squares fits with a stencil of five points
using (a) a cubic polynomial with multipliers $m_{u}=1024$, $m_{d}=1024$ and
$m_{p}=1$, (b) a quadratic polynomial with the same multipliers, and (c) a
quadratic polynomial with multipliers $m_{u}=1024$, $m_{d}=1$ and $m_{p}=1$.
Notice that the curves in (a) and (b) fit almost exactly through the upwind
and downwind points immediately adjacent to the $y$-axis, but in (c) the curve
fits almost exactly only through the upwind point immediately to the left of
the $y$-axis. The point data are labelled with their respective weights.
Points that have failed one of the stability conditions in equation (16) are
marked in red with italicised labels. The upwind point is located at
$(-1,1.8)$ and the downwind point at $(0.62,1.9)$, and the peripheral points
are at $(-2.8,2.4)$, $(-1.6,2.7)$ and $(-1.2,2.2)$. The stabilisation
procedure (section 2.1.4) calculates weights using only $x$ positions, and
values of $\phi$ are included here for illustration only.

To illustrate the stabilisation procedure, figure 4a presents a one-
dimensional example of a cubic polynomial fitted through five points, with the
weight at each point printed beside it. The stabilisation procedure only uses
the $x$ positions of these points and does not use the values of $\phi$
themselves. The $\phi$ values are included here for illustration only. Hence,
for a given set of $x$ positions, the same set of weights are chosen
irrespective of the $\phi$ values.

For a one-dimensional cubic polynomial fit, the list of candidate polynomials
in preference order is

$\displaystyle\phi$ $\displaystyle=a_{1}+a_{2}x+a_{3}x^{2}+a_{4}x^{3}\text{
,}$ (17) $\displaystyle\phi$ $\displaystyle=a_{1}+a_{2}x+a_{3}x^{2}\text{ ,}$
(18) $\displaystyle\phi$ $\displaystyle=a_{1}+a_{2}x\text{ ,}$ (19)
$\displaystyle\phi$ $\displaystyle=a_{1}\text{ .}$ (20)

We begin with the cubic equation (17). The multipliers are chosen so that the
polynomial passes almost exactly through the upwind and downwind points that
are immediately to the left and right of the $y$-axis respectively. The
stability condition on the upwind point is violated because $w_{u}=1.822>1$
(equation 16a). Reducing the downwind multiplier does not help to satisfy the
stability condition, so we start again with the quadratic equation (18), and
the new fit is presented in figure 4b. Again, the multipliers are chosen to
force the polynomial through the upwind and downwind points, but this violates
the stability condition on the downwind point because $w_{d}=0.502>0.5$
(equation 16b). This time, however, stable weights are found by reducing
$m_{d}$ to one (figure 4c) and these are the weights that will be used to
approximate $\phi_{F}$, where the polynomial intercepts the $y$-axis.

####  2.1.5 Future extension to three dimensions

All the procedures used in the cubicFit scheme generalise to three dimensions.
The stencil construction procedure described in section 2.1.1 creates a
stencil with $12$ cells for a face in the interior of a two-dimensional
rectangular mesh. In three dimensions, the same procedure creates a stencil
with $3\times 12=36$ cells. A three-dimensional stencil has three times as
many cells as its two-dimensional counterpart if the mesh has prismatic cells
arranged in columns. Hence, the computational cost during integration
increases three-fold when moving from two dimensions to three dimensions.

To extend the least squares fit to three dimensions, the two-dimensional
polynomial in equation (5) is replaced with its three-dimensional counterpart,

$\phi=a_{1}+a_{2}x+a_{3}y+a_{4}z+a_{5}x^{2}+a_{6}xy+a_{7}y^{2}+a_{8}xz+a_{9}yz+a_{10}z^{2}+\\\
a_{11}x^{3}+a_{12}x^{2}y+a_{13}xy^{2}+a_{14}x^{2}z+a_{15}xz^{2}+a_{16}yz^{2}+a_{17}y^{2}z+a_{18}xyz\text{
.}$ (21)

The procedure for generating candidate polynomials described in section 2.1.3
results in 26 dense subsets in two dimensions and 842 dense subsets in three
dimensions. Note that the combinatorial explosion of dense subsets in three
dimensions does not increase the computational cost during integration.

The stabilisation procedure described in section 2.1.4 requires further
numerical experiments to verify that it is sufficient for three-dimensional
flows and arbitrary polyhedral meshes. An initial three-dimensional test with
uniform flow and a uniform Cartesian mesh obtained a numerically stable
result. For stencils in the interior of the domain, the least squares fit
includes all polynomial terms in equation (21) with $m_{u}=m_{d}=2^{10}$. The
stabilisation procedure does not modify the least squares fit for these
stencils, but we have not explored the three-dimensional extension of cubicFit
any further.

###  2.2 Multidimensional linear upwind transport scheme

The multidimensional linear upwind scheme, called “linearUpwind” hereafter, is
documented here since it provides a baseline accuracy for the experiments in
section 3. The approximation of $\phi_{F}$ is calculated using a gradient
reconstruction,

$\displaystyle\phi_{F}$
$\displaystyle=\phi_{u}+\nabla_{c}\>\phi\cdot\left(\mathbf{x}_{f}-\mathbf{x}_{c}\right)$
(22)

where $\phi_{u}$ is the upwind value of $\phi$, and $\mathbf{x}_{f}$ and
$\mathbf{x}_{c}$ are the position vectors of the face centroid and cell
centroid respectively. The gradient $\nabla_{c}\>\phi$ is calculated using
Gauss’ theorem:

$\displaystyle\nabla_{c}\>\phi=\frac{1}{\mathcal{V}_{c}}\sum_{f\in\>c}\widetilde{\phi}_{F}\mathbf{S}_{f}$
(23)

where $\widetilde{\phi}_{F}$ is linearly interpolated from the two
neighbouring cells of face $f$. The resulting stencil comprises all cells
sharing a face with the upwind cell, including the upwind cell itself. For a
face in the interior of a two-dimensional rectangular mesh, the stencil for
the linearUpwind scheme is a ‘$+$’ shape with 5 cells. On the same mesh, the
stencil for the cubicFit scheme is more than twice the size with 12 cells. For
cells adjacent to boundaries having zero gradient boundary conditions, the
boundary value is set to be equal to the cell centre value before equation
(23) is evaluated. This implementation of the multidimensional linear upwind
scheme is included in the OpenFOAM software distribution [45].

##  3 Results

Three idealised numerical tests are performed to compare the accuracy of the
cubicFit transport scheme with the multidimensional linear upwind scheme and
with other transport schemes in the literature. The first test transports a
tracer horizontally on two-dimensional, highly-distorted terrain-following
meshes, following Schär et al. [2]. The second is a new test case that
modifies the velocity field and tracer position from the first test in order
to challenge the stability and accuracy of the transport schemes near
mountainous lower boundaries. The third test evaluates the cubicFit scheme on
hexagonal-icosahedral meshes and cubed-sphere meshes with deformational flow
on a spherical Earth, as specified by Lauritzen et al. [43].

We have implemented the cubicFit transport scheme and the numerical test cases
using the OpenFOAM CFD library because it enables a like-for-like comparison
between mesh types and transport schemes. We provide source code archives for
the OpenFOAM implementation of the cubicFit scheme [46], the ASAM cut cell
mesh generator [47] and associated OpenFOAM converter [48], and the hexagonal-
icosahedral mesh generator [49]. For the numerical test cases presented here
we also supply the source code [50] and result data [51].

###  3.1 Horizontal transport over mountains

A two-dimensional transport test was developed by Schär et al. [2] to study
the effect of terrain-following coordinate transformations on numerical
accuracy. In this standard test, a tracer is positioned aloft and transported
horizontally over wave-shaped mountains. The test challenges transport schemes
because the tracer must cross mesh layers, which acts to reduce numerical
accuracy [2, 13]. Here we use a more challenging variant of this test that has
steeper mountains and highly-distorted terrain-following meshes. Convergence
results are compared using the linearUpwind and cubicFit transport schemes.

The domain is defined on a rectangular $x$–$z$ plane that is
$301\text{\,}\mathrm{km}$ wide and $25\text{\,}\mathrm{km}$ high as measured
between parallel boundary edges. Boundary conditions are imposed on the tracer
density $\phi$ such that
$\phi=$0\text{\,}\mathrm{kg}\text{\,}{\mathrm{m}}^{-3}$$ at the inlet
boundary, and a zero normal gradient $\partial\phi/\partial
n=$0\text{\,}\mathrm{kg}\text{\,}{\mathrm{m}}^{-4}$$ is imposed at the outlet
boundary. There is no normal flow at the lower and upper boundaries. A range
of mesh spacings are chosen so that $\Delta x\mathbin{:}\Delta
z=2\mathbin{:}1$ to match the original test specification from Schär et al.
[2].

The terrain is wave-shaped, specified by the surface elevation $h$ such that

$\displaystyle h(x)$ $\displaystyle=h^{\star}\cos^{2}(\alpha x)$ (24a) where
$\displaystyle h^{\star}(x)$ $\displaystyle=\left\\{\begin{array}[]{l
l}h_{0}\cos^{2}(\beta x)&\quad\text{if $|x|<a$}\\\
0&\quad\text{otherwise}\end{array}\right.$ (24d)

where $a=$25\text{\,}\mathrm{km}$$ is the mountain envelope half-width,
$h_{0}=$6\text{\,}\mathrm{km}$$ is the maximum mountain height,
$\lambda=$8\text{\,}\mathrm{km}$$ is the wavelength, $\alpha=\pi/\lambda$ and
$\beta=\pi/(2a)$. Note that, in order to make this test more challenging, the
mountain height $h_{0}$ is double the mountain height used by [2].

A basic terrain-following (BTF) mesh is constructed by using the terrain
profile to modify the uniform mesh. The BTF method uses a linear decay
function so that mesh surfaces become horizontal at the top of the model
domain [5],

$z(x)=\left(H-h(x)\right)\left(z^{\star}/H\right)+h(x)$ (25)

where $z$ is the geometric height, $H$ is the height of the domain, $h(x)$ is
the surface elevation and $z^{\star}$ is the computational height of a mesh
surface. If there were no terrain then $h=0$ and $z=z^{\star}$.

A velocity field is prescribed with uniform horizontal flow aloft and zero
flow near the ground,

$\displaystyle u(z)=u_{0}\left\\{\begin{array}[]{l l}1&\quad\text{if $z\geq
z_{2}$}\\\
\sin^{2}\left(\frac{\pi}{2}\frac{z-z_{1}}{z_{2}-z_{1}}\right)&\quad\text{if
$z_{1}<z<z_{2}$}\\\ 0&\quad\text{otherwise}\end{array}\right.$ (29)

where $u_{0}=$10\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$$,
$z_{1}=$7\text{\,}\mathrm{km}$$ and $z_{2}=$8\text{\,}\mathrm{km}$$.

A tracer with density $\phi$ has the shape

$\displaystyle\phi(x,z)$ $\displaystyle=\phi_{0}\left\\{\begin{array}[]{l
l}\cos^{2}\left(\frac{\pi r}{2}\right)&\quad\text{if $r\leq 1$}\\\
0&\quad\text{otherwise}\end{array}\right.$ (30c) with radius $r$ given by
$\displaystyle r$
$\displaystyle=\sqrt{\left(\frac{x-x_{0}}{A_{x}}\right)^{2}+\left(\frac{z-z_{0}}{A_{z}}\right)^{2}}$
(30d)

where $A_{x}=$25\text{\,}\mathrm{km}$$, $A_{z}=$3\text{\,}\mathrm{km}$$ are
the horizontal and vertical half-widths respectively, and
$\phi_{0}=$1\text{\,}\mathrm{kg}\text{\,}{\mathrm{m}}^{-3}$$ is the maximum
density of the tracer. At $t=$0\text{\,}\mathrm{s}$$, the tracer is centred at
$(x_{0},z_{0})=($-50\text{\,}\mathrm{km}$,$12\text{\,}\mathrm{km}$)$ so that
the tracer is upwind of the mountain, in the region of uniform flow above
$z_{2}$.

Figure 5: Numerical convergence of the two-dimensional horizontal transport
test over mountains. $\ell_{2}$ errors (equation 31) and $\ell_{\infty}$
errors (equation 32) are marked at mesh spacings between
$5000\text{\,}\mathrm{m}$ and $250\text{\,}\mathrm{m}$ using linearUpwind and
cubicFit transport schemes on basic terrain-following meshes.

Tests are integrated for $10\,000\text{\,}\mathrm{s}$ using time-steps chosen
for each mesh so that the maximum Courant number is about $0.4$. This choice
yields a time-step that is well below any stability limit, as recommended by
Lauritzen et al. [43]. By the end of integration the tracer is positioned
downwind of the mountain. The analytic solution at
$t=$10\,000\text{\,}\mathrm{s}$$ is centred at
$(x_{0},z_{0})=($50\text{\,}\mathrm{km}$,$12\text{\,}\mathrm{km}$)$. Error
norms are calculated by subtracting the analytic solution from the numerical
solution,

$\displaystyle\ell_{2}$
$\displaystyle=\sqrt{\frac{\sum_{c}\left(\phi-\phi_{T}\right)^{2}\mathcal{V}_{c}}{\sum_{c}\left(\phi_{T}^{2}\mathcal{V}_{c}\right)}}$
(31) $\displaystyle\ell_{\infty}$
$\displaystyle=\frac{\max_{c}|\phi-\phi_{T}|}{\max_{c}|\phi_{T}|}$ (32)

where $\phi$ is the numerical value, $\phi_{T}$ is the analytic value,
$\sum_{c}$ denotes a summation over all cells $c$ in the domain, and
$\max_{c}$ denotes a maximum value of any cell.

Tests were performed using the linearUpwind and cubicFit schemes at mesh
spacings between $\Delta x=$250\text{\,}\mathrm{m}$$ and $\Delta
x=$5000\text{\,}\mathrm{m}$$. Numerical convergence in the $\ell_{2}$ and
$\ell_{\infty}$ norms is plotted in figure 5. The linearUpwind and cubicFit
schemes are second-order convergent at all but the coarsest mesh spacings
where errors are saturated for both schemes.

The cubicFit scheme achieves a given $\ell_{2}$ error using a mesh spacing
that is almost twice as coarse as that needed by the linearUpwind scheme.
Doubling the mesh spacing results in a coarser mesh with four times fewer
cells because the $\Delta x\mathbin{:}\Delta z$ aspect ratio is fixed. Recall
that the stencil for the cubicFit scheme has about twice as many cells as the
stencil for the linearUpwind scheme. Hence, for a given $\ell_{2}$ error, the
computational cost during integration of the cubicFit scheme is about half the
computational cost of the linearUpwind scheme.

This test demonstrates that cubicFit is second-order convergent in the domain
interior irrespective of mesh distortions. The cubicFit scheme achieves In the
next test, we assess the numerical accuracy of the transport schemes near a
distorted, mountainous lower boundary.

###  3.2 Transport over a mountainous lower boundary

The horizontal transport test in the previous section is useful for assessing
numerical accuracy on terrain-following meshes, but it presents no particular
challenge on cut cell meshes because there is no flow through the distorted
cut cells near the ground [52]. Here we present another variant of the
standard horizontal transport test that challenges transport schemes on all
mesh types. By positioning the tracer next to the ground and modifying the
velocity field, we can assess the accuracy of the cubicFit scheme near the
lower boundary. Results using the cubicFit scheme are compared with the
linearUpwind scheme on basic terrain-following, cut cell and slanted cell
meshes.

Cut cell meshes are constructed using the ASAM grid generator [53, 54].
Slanted cell meshes are constructed following the approach by Shaw and Weller
[13]: vertices that are underground are moved up to the surface and zero-area
faces and zero-volume cells are removed. Unlike [13], vertices are never moved
downwards.

Following Schär et al. [2], the domain is $301\text{\,}\mathrm{km}$ wide and
$25\text{\,}\mathrm{km}$ high as measured between parallel boundary edges,
with a mesh spacing of $\Delta x=$1000\text{\,}\mathrm{m}$$ and $\Delta
z=$500\text{\,}\mathrm{m}$$. The boundary conditions are the same as those
used in section 3.1. Cell edges in the central region of the domain are shown
in figure 6 for each of the three mesh types. Cells in the BTF mesh are highly
distorted over steep slopes (figure 6a) while the cut cell mesh (figure 6b)
and slanted cell mesh (figure 6c) are orthogonal everywhere except for cells
nearest the ground.

Figure 6: Two dimensional $x$-$z$ meshes created with the (a) basic terrain-
following, (b) cut cell, and (c) slanted cell methods, and used for the tracer
transport tests in section 3.2. Cell edges are marked by thin black lines. The
peak mountain height $h_{0}=$5\text{\,}\mathrm{km}$$. The velocity field is
the same for all mesh types with streamlines marked on each panel by thick red
lines. The velocity field (equation 33) follows the lower boundary and becomes
entirely horizontal above $H_{1}=$10\text{\,}\mathrm{km}$$. Only the lowest
$10\text{\,}\mathrm{km}$ for the central region of the domain is shown. The
entire domain is $301\text{\,}\mathrm{km}$ wide and $25\text{\,}\mathrm{km}$
high.

Similar to the approach by [13], a velocity field is chosen that follows the
terrain at the surface and becomes entirely horizontal aloft. A streamfunction
$\Psi$ is used so that the discrete velocity field is non-divergent, such that

$\Psi(x,z)=-u_{0}H_{1}\frac{z-h}{H_{1}-h}$ (33)

where $u_{0}=$10\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$$, which is the
horizontal velocity where $h(x)=0$. There is no normal flow at the lower and
upper boundaries. The velocity field becomes purely horizontal above
$H_{1}=$10\text{\,}\mathrm{km}$$ and, elsewhere, the flow is predominantely
horizontal, with non-zero vertical velocities only above sloping terrain. The
horizontal and vertical components of velocity, $u$ and $w$, are given by

$\displaystyle u$ $\displaystyle=-\frac{\partial\Psi}{\partial
z}=u_{0}\frac{H_{1}}{H_{1}-h},\quad w=\frac{\partial\Psi}{\partial
x}=u_{0}H_{1}\frac{\mathrm{d}h}{\mathrm{d}x}\frac{H_{1}-z}{\left(H_{1}-h\right)^{2}}\text{
,}$ (34) $\displaystyle\frac{\mathrm{d}h}{\mathrm{d}x}$
$\displaystyle=-h_{0}\left[\beta\cos^{2}\left(\alpha x\right)\sin\left(2\beta
x\right)+\alpha\cos^{2}\left(\beta x\right)\sin\left(2\alpha
x\right)\right]\text{ .}$ (35)

Unlike the horizontal transport test in [2], the velocity field presented here
extends from the top of the domain all the way to the ground.

The flow is deliberately misaligned with the BTF, cut cell and slanted cell
meshes away from the ground (figure 6) to ensure that flow always crosses mesh
surfaces in order to challenge the transport scheme. The value of $H_{1}$ is
chosen to be much smaller than the domain height $H$ in equation (25) so that
flow crosses the surfaces of the BTF mesh. This is evident in figure 6a where
the the velocity streamlines are tangential to the mesh only at the ground.

| Peak mountain height $h_{0}$ ($\mathrm{km}$)
---|---
Mesh type | $0$ | $3$ | $4$ | $5$ | $6$
BTF | $40$ | $16$ | $10$ | $8$ | $5$
Cut cell | $40$ | $1.6$ | $1.6$ | $0.5$ | $1.5$
Slanted cell | $40$ | $8$ | $6.25$ | $5$ | $4$
Table 1: Time-steps ($\mathrm{s}$) for the two-dimensional transport test over
a mountainous lower boundary. The time-steps were chosen so that the maximum
Courant number was between $0.36$ and $0.46$.

The tracer is again defined by equation (30) but is now positioned at the
ground with
$(x_{0},z_{0})=($-50\text{\,}\mathrm{km}$,$0\text{\,}\mathrm{km}$)$ with half-
widths $A_{x}=$25\text{\,}\mathrm{km}$$ and $A_{z}=$10\text{\,}\mathrm{km}$$.
Tests are integrated forward for $10\,000\text{\,}\mathrm{s}$. The time-step
was chosen for each mesh so that the maximum Courant number was about $0.4$
(table 1). An analytic solution at $10\,000\text{\,}\mathrm{s}$ is obtained by
calculating the new horizontal position of the tracer. Integrating along the
trajectory yields $t$, the time taken to move from the left side of the
mountain to the right [13]:

$\displaystyle\mathrm{d}t$ $\displaystyle=\mathrm{d}x/u(x)$ (36)
$\displaystyle t$
$\displaystyle=\int_{0}^{x}\frac{H-h(x)}{u_{0}H}\>\mathrm{d}x$ (37)
$\displaystyle t$
$\displaystyle=\frac{x}{u_{0}}-\frac{h_{0}}{16u_{0}H}\left[4x+\frac{\sin
2(\alpha+\beta)x}{\alpha+\beta}\frac{\sin
2(\alpha-\beta)x}{\alpha-\beta}+2\left(\frac{\sin 2\alpha
x}{\alpha}+\frac{\sin 2\beta x}{\beta}\right)\right]$ (38)

By solving this equation we find that
$x(t=$10\,000\text{\,}\mathrm{s}$)=$54\,342.8\text{\,}\mathrm{m}$$ when
$h_{0}=$5\text{\,}\mathrm{km}$$.

The tracer density boundary conditions are the same as those in the previous
test. Since the cubicFit transport scheme uses values at boundaries with
Dirichlet boundary conditions, the cubicFit scheme uses only inlet boundary
values in this test case.

Figure 7: Evolution of the tracer in the two-dimensional transport test over a
mountainous lower boundary. The tracer is transported to the right over the
wave-shaped terrain. Tracer contours are every
$0.1\text{\,}\mathrm{kg}\text{\,}{\mathrm{m}}^{-3}$. The result obtained using
the cubicFit scheme on the basic terrain-following mesh is shown at
$t=$0\text{\,}\mathrm{s}$$, $t=$5000\text{\,}\mathrm{s}$$ and
$t=$10\,000\text{\,}\mathrm{s}$$ with solid black contours. The analytic
solution at $t=$10\,000\text{\,}\mathrm{s}$$ is shown with dotted contours.
The shaded box indicates the region that is plotted in figure 8.

Three series of tests were performed using similar configurations. The first
series uses a peak mountain height of $h_{0}=$5\text{\,}\mathrm{km}$$ to
examine errors on different mesh types using the two transport schemes. The
second series varies the peak mountain height to examine the sensitivity of
the transport schemes to mesh distortions. The third series verifies accuracy
at Courant numbers close to the limit of stability, and examines the longest
stable time-step for different mesh types.

For the first series of tests with $h_{0}=$5\text{\,}\mathrm{km}$$, tracer
contours at the initial time $t=$0\text{\,}\mathrm{s}$$, half-way time
$t=$5000\text{\,}\mathrm{s}$$, and end time $t=$10\,000\text{\,}\mathrm{s}$$
are shown in figure 7 using the cubicFit scheme on the BTF mesh. As apparent
at $t=$5000\text{\,}\mathrm{s}$$, the tracer is distorted by the terrain-
following velocity field as it passes over the mountain, but its original
shape is restored once it has cleared the mountain by
$t=$10\,000\text{\,}\mathrm{s}$$. A small phase lag is apparent when the
numerical solution marked with solid contour lines is compared with the
analytic solution marked with dotted contour lines.

Figure 8: Tracer contours at $t=$10\,000\text{\,}\mathrm{s}$$ for the two-
dimensional tracer transport tests over a mountainous lower boundary. A region
in the lee of the mountain is plotted corresponding to the shaded area in
figure 7. Results are presented on BTF, cut cell and slanted cell meshes
(shown in figure 6) using the linearUpwind and cubicFit transport schemes. The
numerical solutions are marked by solid black lines. The analytic solution is
marked by dotted lines. Contours are every
$0.1\text{\,}\mathrm{kg}\text{\,}{\mathrm{m}}^{-3}$.

Numerical errors are more clearly revealed by subtracting the analytic
solution from the numerical solution. Error fields are compared between BTF,
cut cell and slanted cell meshes using the linearUpwind scheme (figures 8a, 8b
and 8c respectively) and the cubicFit scheme (figures 8d, 8e and 8f
respectively). Results are least accurate using the linearUpwind scheme on the
slanted cell mesh (figure 8c). The final tracer is slightly distorted and does
not extend far enough towards the ground. The $\ell_{\infty}$ error magnitude
is reduced by using the linearUpwind scheme on the cut cell mesh (figure 8b),
but the shape of the error remains the same. The cubicFit scheme is less
sensitive to the choice of mesh with similar error magnitudes on the BTF mesh
(figure 8d), cut cell mesh (figure 8e) and slanted cell mesh (figure 8f).
Errors using the cubicFit scheme on cut cell and slanted cell meshes are much
smaller than the errors using the linearUpwind scheme on the same meshes.

Figure 9: Error measures for the two-dimensional tracer transport tests over a
mountainous lower boundary. Peak mountain heights $h_{0}$ are from
$0\text{\,}\mathrm{km}6\text{\,}\mathrm{km}$. Results are compared on BTF, cut
cell and slanted cell meshes using the linearUpwind and the cubicFit schemes.
At $h_{0}=$0\text{\,}\mathrm{km}$$ the terrain is entirely flat and the BTF,
cut cell and slanted cell meshes are identical. At
$h_{0}=$6\text{\,}\mathrm{km}$$ the linearUpwind scheme is unstable on the
slanted cell mesh.

To further examine the performance of the cubicFit scheme in the presence of
steep terrain, a second series of tests were performed in which the peak
mountain height was varied from $0\text{\,}\mathrm{km}6\text{\,}\mathrm{km}$
keeping all other parameters constant. Results were obtained on BTF, cut cell
and slanted cell meshes using the linearUpwind scheme and cubicFit scheme.
Again, the time-step was chosen for each test so that the maximum Courant
number was about $0.4$ (table 1). The $\ell_{2}$ error was calculated by
subtracting the analytic solution from the numerical solution (figure 9). Note
that the analytic solution is a function of mountain height, with the tracer
travelling farther over higher mountains due to non-divergent flow through a
narrower channel. In all cases, error increases with increasing mountain
height because steeper slopes lead to greater mesh distortions. Errors are
identical for a given transport scheme when $h_{0}=$0\text{\,}\mathrm{km}$$
and the ground is entirely flat because the BTF, cut cell and slanted cell
meshes are identical. Compared with the cubicFit scheme, the linearUpwind
scheme is more sensitive to the mesh type and mountain height. The
linearUpwind scheme is unstable on the slanted cell mesh with a peak mountain
height $h_{0}=$6\text{\,}\mathrm{km}$$ despite using a Courant number of
$0.428$. In contrast, the cubicFit scheme is less sensitive to the mesh type
and errors grow more slowly with increasing mountain height. The cubicFit
scheme yields stable results in all tests.

Figure 10: (a) Longest stable time-steps, $\Delta t_{\mathrm{max}}$, and (b)
largest stable maximum Courant numbers, $\max(\mathrm{Co})$, for the two-
dimensional tracer transport test over a mountainous lower boundary. Results
were obtained on basic terrain-following, cut cell and slanted cell meshes at
mesh spacings between $\Delta x=$5000\text{\,}\mathrm{m}$$ and $\Delta
x=$250\text{\,}\mathrm{m}$$. The largest stable maximum Courant numbers were
calculated from the corresponding longest stable time-steps using equation
(4).

A final series of tests were performed to determine the stability limit of the
cubicFit scheme with the two-stage Heun time-stepping scheme (equation 2). The
tracer was transported on BTF, slanted cell and cut cell meshes with a variety
of mesh spacings between $\Delta x=$5000\text{\,}\mathrm{m}$$ and $\Delta
x=$125\text{\,}\mathrm{m}$$. $\Delta z$ was chosen so that a constant aspect
ratio is preserved such that $\Delta x/\Delta z=2$. For each test, the time-
step was increased until the result became unstable. The largest stable time-
steps, $\Delta t_{\mathrm{max}}$, are presented in figure 10a. BTF meshes
permit the longest time-steps of all three mesh types since cells are almost
uniform in volume. As expected, the longest stable time-step scales linearly
with BTF mesh spacing. There is no such linear scaling on cut cell meshes
because these meshes can have arbitrarily small cells. The time-step
constraints on cut cell meshes are the most severe of the three mesh types.
Slanted cell meshes have a slightly stronger time-step constraint than BTF
meshes but still exhibit similar linear scaling with mesh spacing.
Furthermore, a dynamical model that uses slanted cell meshes instead of BTF
meshes is expected to calculate pressure gradients more accurately [13].

Figure 10b presents the largest stable maximum Courant numbers,
$\max(\mathrm{Co})$, which were calculated by substituting $\Delta t=\Delta
t_{\mathrm{max}}$ into equation (4). On basic terrain following meshes, the
maximum Courant number tends towards about $1.3$ with finer mesh spacings. No
such trend is found on cut cell or slanted cell meshes. Cut cell meshes permit
the largest maximum Courant numbers of around $3$, but the largest stable
time-steps on cut cell meshes are still smaller than corresponding time-steps
on basic terrain following and slanted cell meshes.

This paper focuses on the spatial discretisation of the cubicFit scheme, but
the stability limit depends also upon the choice of time-stepping. As such, we
have not calculated a theoretical Courant number limit, although such an
analysis should be possible using the techniques in [55].

The transport tests presented in this section demonstrate that the cubicFit
scheme is suitable for flows over very steep terrain on two-dimensional
terrain-following, cut cell and slanted cell meshes. The cubicFit scheme is
less sensitive to the mesh type and mountain steepness compared to the
linearUpwind scheme. The linearUpwind scheme becomes unstable over very steep
slopes but the cubicFit scheme is stable for all tests. The accuracy of the
cubicFit scheme was largely insensitive to the choice of time-step. In the
next section, we evaluate the cubicFit scheme using more complex,
deformational flows on icosahedral meshes and cubed-sphere meshes.

###  3.3 Deformational flow on a sphere

The tests so far have used flows that are mostly uniform on meshes that are
based on rectangular cells. To ensure that the cubicFit transport scheme is
suitable for complex flows on a variety of meshes, we use a standard test of
deformational flow on a spherical Earth, as specified by Lauritzen et al.
[43]. Results are compared between linearUpwind and cubicFit schemes using
hexagonal-icosahedral meshes and cubed-sphere meshes. Hexagonal-icosahedral
meshes are constructed by successive refinement of a regular icosahedron
following the approach by [28, 56, 57] without any mesh twisting. Cubed-sphere
meshes are constructed using an equi-distant gnomic projection of a cube
having a uniform Cartesian mesh on each panel [58].

Following appendix A9 in [32], the average equatorial spacing $\Delta\lambda$
is used as a measure of mesh spacing. It is defined as

$\displaystyle\Delta\lambda=$$\frac{\overline{\Delta x}}{2\pi R_{e}}$ (39)

where $\overline{\Delta x}$ is the mean distance between cell centres and
$R_{e}=$6.3712\text{\times}{10}^{6}\text{\,}\mathrm{m}$$ is the radius of the
Earth.

The deformational flow test specified by Lauritzen et al. [43] comprised six
elements:

  1. 1.

a convergence test using a Gaussian-shaped tracer

  2. 2.

a “minimal” resolution test using a cosine bell-shaped tracer

  3. 3.

a test of filament preservation

  4. 4.

a test using a “rough” slotted cylinder tracer

  5. 5.

a test of correlation preservation between two tracers

  6. 6.

a test using a divergent velocity field

We assess the cubicFit scheme using the first two tests only. We do not
consider filament preservation, correlation preservation, or the transport of
a “rough” slotted cylinder because no shape-preserving filter has yet been
developed for the cubicFit scheme. Stable results were obtained when testing
the cubicFit scheme using a divergent velocity field, but no further analysis
is made here.

The first deformational flow test uses an infinitely continuous initial tracer
that is transported in a non-divergent, time-varying, rotational velocity
field. The velocity field deforms two Gaussian ‘hills’ of tracer into thin
vortical filaments. Half-way through the integration the rotation reverses so
that the filaments become circular hills once again. The analytic solution at
the end of integration is identical to the initial condition. A rotational
flow is superimposed on a time-invariant background flow in order to avoid
error cancellation. The non-divergent velocity field is defined by the
streamfunction $\Psi$,

$\displaystyle\Psi(\lambda,\theta,t)=\frac{10R_{e}}{T}\sin^{2}\left(\lambda^{\prime}\right)\cos^{2}\left(\theta\right)\cos\left(\frac{\pi
t}{T}\right)-\frac{2\pi R_{e}}{T}\sin\left(\theta\right)$ (40)

where $\lambda$ is a longitude, $\theta$ is a latitude,
$\lambda^{\prime}=\lambda-2\pi t/T$, and $T=$12$$ days is the duration of
integration. The time-step is chosen such that the maximum Courant number is
about 0.4.

The initial tracer density $\phi$ is defined as the sum of two Gaussian hills,

$\displaystyle\phi=\phi_{1}(\lambda,\theta)+\phi_{2}(\lambda,\theta)\text{ .}$
(41)

An individual hill $\phi_{i}$ is given by

$\displaystyle\phi_{i}(\lambda,\theta)=\phi_{0}\exp\left(-b\left(\frac{|\mathbf{x}-\mathbf{x}_{i}|}{R_{e}}\right)^{2}\right)$
(42)

where $\phi_{0}=$0.95\text{\,}\mathrm{kg}\text{\,}{\mathrm{m}}^{-3}$$ and
$b=5$. The Cartesian position vector $\mathbf{x}=(x,y,z)$ is related to the
spherical coordinates $(\lambda,\theta)$ by

$\displaystyle(x,y,z)=(R_{e}\cos\theta\cos\lambda,R_{e}\cos\theta\sin\lambda,R_{e}\sin\theta)\text{
.}$ (43)

The centre of hill $i$ is positioned at $\mathbf{x}_{i}$. In spherical
coordinates, two hills are centred at

$\displaystyle(\lambda_{1},\theta_{1})$ $\displaystyle=(5\pi/6,0)$ (44)
$\displaystyle(\lambda_{2},\theta_{2})$ $\displaystyle=(7\pi/6,0)$ (45)

Figure 11: Tracer fields for the deformational flow test using initial
Gaussian hills. The tracer is deformed by the velocity field before the
rotation reverses to return the tracer to its original distribution: (a) the
initial tracer distribution at $t=$0\text{\,}\mathrm{s}$$; (b) by $t=T/2$ the
Gaussian hills are stretched into a thin S-shaped filament; (c) at $t=T$ the
tracer resembles the initial Gaussian hills except for some distortion and
diffusion due to numerical errors. Results were obtained with the cubicFit
scheme on a hexagonal-icosahedral mesh with an average equatorial mesh spacing
of $\Delta\lambda=$$$.

The results in figure 11 are obtained using the cubicFit scheme on a
hexagonal-icosahedral mesh with $\Delta\lambda=$$$. The initial Gaussian hills
are shown in figure 11a. At $t=T/2$ the tracer has been deformed into an
S-shaped filament (figure 11b). By $t=T$ the tracer has almost returned to its
original distribution except for some slight distortion and diffusion that are
the result of numerical errors (figure 11c).

Figure 12: Numerical convergence of the deformational flow test on the sphere
using initial Gaussian hills. $\ell_{2}$ errors (equation 31) and
$\ell_{\infty}$ errors (equation 32) are marked at mesh spacings between
$None$ and $None$ using the linearUpwind scheme (dotted lines) and the
cubicFit scheme (solid lines) on hexagonal-icosahedral meshes and cubed-sphere
meshes.

To determine the order of convergence and relative accuracy of the
linearUpwind and cubicFit schemes, the same test was performed at a variety of
mesh spacings betweeen $\Delta\lambda=$$$ and $\Delta\lambda=$$$ on hexagonal-
icosahedral meshes and cubed-sphere meshes. The results are shown in figure
12. The solution is slow to converge at coarse resolutions, and this behaviour
agrees with the results from Lauritzen et al. [43]. Both linearUpwind and
cubicFit schemes achieve second-order accuracy at smaller mesh spacings. For
any given mesh type and mesh spacing, the cubicFit scheme is more accurate
than the linearUpwind scheme. Results are more accurate using hexagonal-
icosahedral meshes compared to cubed-sphere meshes. It is not known whether
the larger errors on cubed-sphere meshes are due to mesh non-uniformities at
panel corners but there is no evidence of grid imprinting in the error fields
(not shown).

A slightly more challenging variant of the same test is performed using a
quasi-smooth tracer field defined as the sum of two cosine bells,

$\displaystyle\phi=\begin{cases}b+c\phi_{1}(\lambda,\theta)&\quad\text{if
$r_{1}<r$,}\\\ b+c\phi_{2}(\lambda,\theta)&\quad\text{if $r_{2}<r$,}\\\
b&\quad\text{otherwise.}\end{cases}$ (46)

The velocity field is the same as before. This test is used to determine the
“minimal” resolution, $\Delta\lambda_{m}$, which is specified by Lauritzen et
al. [43] as the coarsest mesh spacing for which $\ell_{2}\approx 0.033$.

Transport scheme | Mesh type | Minimal resolution ($\mathrm{\SIUnitSymbolDegree}$)
---|---|---
linearUpwind | Cubed-sphere | $0.15$
FARSIGHT, grid-point semi-Lagrangian [59] | Cubed-sphere | $0.1875$
linearUpwind | Hexagonal-icosahedral | $0.2$
SLFV-SL, swept-area scheme [60] | Hexagonal-icosahedral | $0.25$
cubicFit | Cubed-sphere | $0.25$
cubicFit | Hexagonal-icosahedral | $0.3$
ICON-FFSL, swept-area scheme [60] | Triangular-icosahedral | $0.42$
Table 2: Minimal resolutions for the cubicFit and linearUpwind schemes in the
test of deformational flow using cosine bells. Italicised values have been
extrapolated using the second-order convergence obtained at coarser mesh
spacings. For comparison with existing models, some results are also included
for unlimited versions of the transport schemes from the intercomparison by
Lauritzen et al. [32].

The minimal resolution for the cubicFit scheme on a hexagonal-icosahedral mesh
is about $\Delta\lambda_{m}=$$$. Tests were not performed at mesh spacings
finer than $\Delta\lambda=$$$ but approximate minimal resolutions have been
extrapolated from the second-order convergence that is found at fine mesh
spacings. These minimal resolutions are presented in table 2 along with a
selection of transport schemes having similar minimal resolutions from the
model intercomparison by Lauritzen et al. [32].

The series of deformational flow tests presented here demonstrate that the
cubicFit scheme is suitable for transport on spherical meshes based on
quadrilaterals and hexagons. The cubicFit scheme is largely insensitive to the
mesh type, and results are more accurate compared to the linearUpwind scheme
for a given mesh type and mesh spacing. Neither scheme requires special
treatment at the corners of cubed-sphere panels.

##  4 Conclusion

Atmospheric models are using increasingly fine horizontal mesh spacings that
resolve steep slopes in terrain resulting in highly-distorted meshes,
increased numerical errors and numerical instabilities. We have presented a
new multidimensional method-of-lines transport scheme, cubicFit, that applies
constraints derived from a von Neumann stability analysis to make the scheme
stable over steep terrain on highly-distorted, arbitrary meshes. The scheme
has a low computational cost at runtime, requiring only $n$ multiplies per
face per time-stage using a stencil with $n$ cells. Stability constraint
calculations are pre-computed during model initialisation since they depend
upon the mesh geometry only.

The cubicFit scheme was compared to a multidimensional linear upwind scheme
using three idealised numerical tests. The first test transported a tracer
horizontally above steep slopes on highly-distorted, two-dimensional terrain-
following meshes. The cubicFit scheme was second-order convergent regardless
of mesh distortions. The second test transported a tracer over a mountainous
lower boundary using terrain-following, cut cell and slanted cell meshes. The
cubicFit scheme was generally insensitive to the type of mesh and less
sensitive to terrain steepness compared to the multidimensional linear upwind
scheme, and the scheme maintained accuracy up to its stability limit. The
third test evaluated the transport schemes in a standard deformational flow
field on hexagonal-icosahedral meshes and cubed-sphere meshes. In all tests,
compared to the multidimensional linear upwind scheme, the cubicFit transport
scheme was more stable and more accurate.

##  5 Acknowledgements

James Shaw acknowledges support from a PhD studentship funded jointly by NERC
grant NE/K500860/1 and the University of Reading with CASE support from the
Met Office. We are grateful to the Leibniz Institute for Tropospheric Research
for providing their cut cell mesh generator, to the three anonymous reviewers
for their helpful questions, and to Dr Shing Hing Man for his assistance with
candidate polynomial generation. We also thank Dr Tristan Pryer at the
University of Reading for useful discussions about the cubicFit transport
scheme.

## Appendix A: One-dimensional von Neumann stability analysis

Two analyses are performed in order to find stability constraints on the
weights $\mathbf{w}=\mathbf{\tilde{b}_{1}^{+}}\cdot\mathbf{m}$ as appear in
equation (14). The first analysis uses a two-cell approximation to derive
separate constraints on the upwind weight $w_{u}$ and downwind weight $w_{d}$.
The second analysis uses three cells to derive a constraint that considers all
weights in a stencil.

### Two-cell analysis

We start with the conservation equation for a dependent variable $\phi$ that
is discrete-in-space and continuous-in-time

$\displaystyle\frac{\partial\phi_{j}}{\partial t}$
$\displaystyle=-v\frac{\phi_{R}-\phi_{L}}{\Delta x}$ (47) where $v$ is the
velocity, and the left and right fluxes, $\phi_{L}$ and $\phi_{R}$, are
weighted averages of the neighbouring cell centres. Assuming that $v$ is
positive $\displaystyle\phi_{L}$
$\displaystyle=\alpha_{u}\phi_{j-1}+\alpha_{d}\phi_{j}$ (48)
$\displaystyle\phi_{R}$ $\displaystyle=\beta_{u}\phi_{j}+\beta_{d}\phi_{j+1}$
(49)

where $\phi_{j-1},\phi_{j},\phi_{j+1}$ are cell centre values, and $j$ denotes
a cell centre position $x=j\Delta x$ where $\Delta x$ is a uniform mesh
spacing. $\alpha_{u}$ and $\beta_{u}$ are the upwind weights and $\alpha_{d}$
and $\beta_{d}$ are the downwind weights for the left and right fluxes
respectively, and $\alpha_{u}+\alpha_{d}=1$ and $\beta_{u}+\beta_{d}=1$.

At a given time $t=n\Delta t$ at time-level $n$ and with a time-step $\Delta
t$, we assume a wave-like solution with an amplification factor $A$, such that

$\displaystyle\phi_{j}^{(n)}$ $\displaystyle=A^{n}e^{{i\mkern 1.0mu}jk\Delta
x}$ (50)

where $\phi_{j}^{(n)}$ denotes a value of $\phi$ at position $j$ and time-
level $n$. Using this to rewrite the left-hand side of equation (47)

$\displaystyle\frac{\partial\phi_{j}}{\partial t}$
$\displaystyle=\frac{\partial}{\partial t}\left(A^{t/\Delta
t}\right)e^{ijk\Delta x}=\frac{\ln A}{\Delta t}A^{n}e^{ikj\Delta x}$ (51)
hence equation (47) becomes $\displaystyle\frac{\ln A}{\Delta t}$
$\displaystyle=-\frac{v}{\Delta x}\left(\beta_{u}+\beta_{d}e^{ik\Delta
x}-\alpha_{u}e^{-ik\Delta x}-\alpha_{d}\right)$ (52) $\displaystyle\ln A$
$\displaystyle=-c\left(\beta_{u}-\alpha_{d}+\beta_{d}\cos k\Delta x+{i\mkern
1.0mu}\beta_{d}\sin k\Delta x-\alpha_{u}\cos k\Delta x+{i\mkern
1.0mu}\alpha_{u}\sin k\Delta x\right)$ (53) where the Courant number
$c=v\Delta t/\Delta x$. Let $\Re=\beta_{u}-\alpha_{d}+\beta_{d}\cos k\Delta
x-\alpha_{u}\cos k\Delta x$ and $\Im=\beta_{d}\sin k\Delta x+\alpha_{u}\sin
k\Delta x$, then $\displaystyle\ln A$ $\displaystyle=-c\left(\Re+{i\mkern
1.0mu}\Im\right)$ (54) $\displaystyle A$ $\displaystyle=e^{-c\Re}e^{-{i\mkern
1.0mu}c\Im}$ (55) and the complex modulus of $A$ is $\displaystyle|A|$
$\displaystyle=e^{-c\Re}=\exp\left(-c\left(\beta_{u}-\alpha_{d}+\left(\beta_{d}-\alpha_{u}\right)\cos
k\Delta x\right)\right)\text{ .}$ (56)

For stability we need $|A|\leq 1$ and, imposing the additional constraints
that $\alpha_{u}=\beta_{u}$ and $\alpha_{d}=\beta_{d}$, then

$\displaystyle\left(\alpha_{u}-\alpha_{d}\right)\left(1-\cos k\Delta x\right)$
$\displaystyle\geq 0\quad\forall k\Delta x$ (57) and, given $0\leq 1-\cos
k\Delta x\leq 2$, then $\displaystyle\alpha_{u}-\alpha_{d}\geq 0\text{ .}$
(58)

Additionally, we do not want more damping than a first-order upwind scheme
(where $\alpha_{u}=\beta_{u}=1$, $\alpha_{d}=\beta_{d}=0$), having an
amplification factor, $A_{\mathrm{up}}$, so we need $\left\lvert
A\right\rvert\geq\left\lvert A_{\mathrm{up}}\right\rvert$, hence

$\displaystyle\exp\left(-c\left(\alpha_{u}-\alpha_{d}\right)\left(1-\cos
k\Delta x\right)\right)$ $\displaystyle\geq\exp\left(-c\left(1-\cos k\Delta
x\right)\right)\quad\forall k\Delta x$ (59) therefore
$\displaystyle\alpha_{u}-\alpha_{d}$ $\displaystyle\leq 1\text{ .}$ (60)

Now, knowing that $\alpha_{u}+\alpha_{d}=1$ (or $\alpha_{d}=1-\alpha_{u}$)
then, using equations (58) and (60),

$\displaystyle 0.5\leq\alpha_{u}$ $\displaystyle\leq 1\text{ and}$ (61)
$\displaystyle 0\leq\alpha_{d}$ $\displaystyle\leq 0.5\text{ .}$ (62)

### Three-cell analysis

We start again from equation (47) but this time approximate $\phi_{L}$ and
$\phi_{R}$ using three cell centre values,

$\displaystyle\phi_{L}$
$\displaystyle=\alpha_{uu}\phi_{j-2}+\alpha_{u}\phi_{j-1}+\alpha_{d}\phi_{j}$
(63) $\displaystyle\phi_{R}$
$\displaystyle=\alpha_{uu}\phi_{j-1}+\alpha_{u}\phi_{j}+\alpha_{d}\phi_{j+1}$
(64)

having used the same weights $\alpha_{uu}$, $\alpha_{u}$ and $\alpha_{d}$ for
both left and right fluxes. Substituting equation (50) into equation (47) we
find

$\displaystyle A=\exp\left(-c\left[\alpha_{uu}\left(e^{-ik\Delta
x}-e^{-2ik\Delta x}\right)+\alpha_{u}\left(1-e^{-ik\Delta
x}\right)+\alpha_{d}\left(e^{ik\Delta x}-1\right)\right]\right)$ (65) so that,
if the complex modulus $\left\lvert A\right\rvert\leq 1$ then
$\displaystyle\alpha_{u}-\alpha_{d}+\left(\alpha_{uu}-\alpha_{u}+\alpha_{d}\right)\cos
k\Delta x-\alpha_{uu}\cos 2k\Delta x\geq 0\text{ .}$ (66)

If $k\Delta x=\pi$ then $\cos k\Delta x=-1$ and $\cos 2k\Delta x=1$ and
$\alpha_{u}-\alpha_{d}\geq\alpha_{uu}$. If $k\Delta x=\pi/2$ then $\cos
k\Delta x=0$ and $\cos 2k\Delta x=-1$ and
$\alpha_{u}-\alpha_{d}\geq-\alpha_{uu}$. Hence we find that

$\displaystyle\alpha_{u}-\alpha_{d}$
$\displaystyle\geq\left\lvert\alpha_{uu}\right\rvert\text{ .}$ (67) When the
same analysis is performed with four cells, $\alpha_{uuu}$, $\alpha_{uu}$,
$\alpha_{u}$ and $\alpha_{d}$, by varying $k\Delta x$ we find that equation
(67) holds replacing $\left\lvert\alpha_{uu}\right\rvert$ with
$\max(\left\lvert\alpha_{uu}\right\rvert,\left\lvert\alpha_{uuu}\right\rvert)$.
Hence, we generalise equation (67) to find the final stability constraint
$\displaystyle\alpha_{u}-\alpha_{d}$
$\displaystyle\geq\max_{p\>\in\>P}|\alpha_{p}|$ (68)

where the peripheral cells $P$ is the set of all stencil cells except for the
upwind cell and downwind cell, and $\alpha_{p}$ is the weight for a given
peripheral cell $p$. We hypothesise that the three stability constraints
(equations 61, 62 and 67) are necessary but not sufficient for a transport
scheme on arbitrary meshes.

The stability of the one-dimensional transport equation discretised in space
and time could be analysed using existing techniques [55], but we have only
analysed the spatial stability of the cubicFit scheme. Numerical experiments
presented in section 3.2 demonstrate that the cubicFit scheme is generally
insensitive to the time-step, provided that it is below a stability limit.

## Appendix B: Mesh geometry on a spherical Earth

The cubicFit transport scheme is implemented using the OpenFOAM CFD library.
Unlike many atmospheric models that use spherical coordinates, OpenFOAM uses
global, three-dimensional Cartesian coordinates with the $z$-axis pointing up
through the North pole. In order to perform the experiments on a spherical
Earth presented in section 3.3, it is necessary for velocity fields and mesh
geometries to be expressed in these global Cartesian coordinates.

### Velocity field specification

The non-divergent velocity field in section 3.3 is specified as a
streamfunction $\Psi(\lambda,\theta)$. Instead of calculating velocity
vectors, the flux $\mathbf{u}_{f}\cdot\mathbf{S}_{f}$ through a face $f$ is
calculated directly from the streamfunction,

$\displaystyle\mathbf{u}_{f}\cdot\mathbf{S}_{f}=\sum_{e\>\in\>f}\mathbf{e}\cdot\mathbf{x}_{e}\Psi(e)$
(69)

where $e\in f$ denotes the edges $e$ of face $f$, $\mathbf{e}$ is the edge
vector joining the two vertices of the edge, $\mathbf{x}_{e}$ is the position
vector of the edge midpoint, and $\Psi(e)$ is the streamfunction evaluated at
the same position. Edge vectors are directed in a counter-clockwise
orientation.

### Spherical mesh construction

Since OpenFOAM does not support two-dimensional spherical meshes, instead, we
construct meshes that have a single layer of cells that are
$2000\text{\,}\mathrm{m}$ deep, having an inner radius
$r_{1}=R_{e}-$1000\text{\,}\mathrm{m}$$ and an outer radius
$r_{2}=R_{e}+$1000\text{\,}\mathrm{m}$$. By default, OpenFOAM meshes comprise
polyhedral cells with straight edges and flat faces. This is problematic for
spherical meshes because face areas and cell volumes are too small. For tests
on a spherical Earth, we override the default configuration and calculate our
own face areas, cell volumes, face centres and cell centres that account for
the mesh curvature. Note that the new centres are no longer centroids, but
they are consistent with the horizontal transport tests on a sphere presented
in section 3.3.

A face is classified as either a surface face or radial face. A surface face
has any number of vertices, all of equal radius. A radial face has four
vertices with two different radii, $r_{1}$ and $r_{2}$, and two different
horizontal coordinates, $(\lambda_{1},\theta_{1})$ and
$(\lambda_{2},\theta_{2})$. A radial face centre is modified so that it has a
radius $R_{e}$. The latitudinal and longitudinal components of a radial face
centre need no modification. The face area $A_{f}$ for a radial face $f$ is
the area of the annular sector,

$\displaystyle A_{f}=\frac{d}{2}\left\lvert r_{2}^{2}-r_{1}^{2}\right\rvert$
(70)

where $d$ is the great-circle distance between $(\lambda_{1},\theta_{1})$ and
$(\lambda_{2},\theta_{2})$.

To calculate the centre of a surface face $f$, a new vertex is created that is
positioned at the mean of the face vertices. Note that this centre position,
$\mathbf{\tilde{c}}_{f}$, is used in intermediate calculations and it is not
the face centre position. Next, the surface face is subdivided into spherical
triangles that share this new vertex [61]. The face centre direction and
radius are calculated separately. The face centre direction $\mathbf{\hat{r}}$
is the mean of the spherical triangle centres weighted by their solid angle,

$\displaystyle\mathbf{\hat{r}}=\frac{\sum_{t\>\in\>f}{\Omega_{t}\left(\mathbf{x}_{t,1}+\mathbf{x}_{t,2}+\mathbf{\tilde{c}}_{f}\right)}}{\left\lvert\sum_{t\>\in\>f}{\Omega_{t}\left(\mathbf{x}_{t,1}+\mathbf{x}_{t,2}+\mathbf{\tilde{c}}_{f}\right)}\right\rvert}$
(71)

where $t\in f$ denotes the spherical triangles $t$ of face $f$, $\Omega_{t}$
is spherical triangle’s solid angle which is calculated using l’Huilier’s
theorem, $\mathbf{x}_{t,1}$ and $\mathbf{x}_{t,2}$ are the positions of the
vertices shared by the face $f$ and spherical triangle $t$, and
$\mathbf{\tilde{c}}_{f}$ is the position of the centre vertex shared by all
spherical triangles of face $f$. The face centre radius $r$ is the mean radius
of the face vertices, again weighted by the solid angle of each spherical
triangle,

$\displaystyle r=\frac{\sum_{t\in
f}{\Omega_{t}\left(\left\lvert\mathbf{x}_{t,1}\right\rvert+\left\lvert\mathbf{x}_{t,2}\right\rvert\right)/2}}{\Omega_{f}}$
(72)

where the solid angle $\Omega_{f}$ of face $f$ is the sum of the solid angles
of the constituent spherical triangles,

$\displaystyle\Omega_{f}=\sum_{t\in f}{\Omega_{t}}\text{ .}$ (73)

We use equations (71) and (72) to calculate the centre $\mathbf{c}_{f}$ of the
face $f$,

$\displaystyle\mathbf{c}_{f}=r\>\mathbf{\hat{r}}$ (74)

The area vector $\mathbf{S}_{f}$ of the surface face $f$ is the sum of the
spherical triangle areas [61],

$\displaystyle\mathbf{S}_{f}=r^{2}\Omega_{f}\mathbf{\hat{r}}\text{ .}$ (75)

Cell centres and cell volumes are corrected by considering faces that are not
normal to the sphere such that

$\displaystyle\frac{\left(\mathbf{S}_{f}\cdot\mathbf{c}_{f}\right)^{2}}{\left\lvert\mathbf{S}_{f}\right\rvert^{2}\left\lvert\mathbf{c}_{f}\right\rvert^{2}}>0\text{
.}$ (76)

Let $\mathcal{F}$ be the set of faces satisfying equation (76). Then, the cell
volume $\mathcal{V}_{c}$ is

$\displaystyle\mathcal{V}_{c}=\frac{1}{3}\sum_{f\>\in\>\mathcal{F}}\mathbf{S}_{f}\cdot\mathbf{c}_{f}$
(77)

which can be thought of as the area $A$ integrated between $r_{1}$ and $r_{2}$
such that
$\int_{0}^{R}{A(r)\>\mathrm{d}r}=\int_{r_{1}}^{r_{2}}{r^{2}\Omega\>\mathrm{d}r}=\frac{1}{3}\Omega\left(r_{2}^{3}-r_{1}^{3}\right)$.
The cell centre is modified so that it has a radius $R_{e}$, which is
consistent with radial faces.

Edges can be classified in a similar manner to faces where surface edges are
tangent to the sphere and radial faces are normal to the sphere. The edge
midpoints $\mathbf{x}_{e}$ are used to calculate the face flux for non-
divergent velocity fields (equation 69). For transport tests, corrections to
edge midpoints are unnecessary. Due to the choice of $r_{1}$ and $r_{2}$
during mesh construction, the midpoint of a radial edge is at a radial
distance of $R_{e}$ which is necessary for the correct calculation of non-
divergent velocity fields. The position of surface edge midpoints is
unimportant because these edges do not contribute to the face flux since
$\mathbf{e}\cdot\mathbf{x}_{e}=0$. Edge lengths are the straight-line distance
between the two vertices and not the great-circle distance. Again, the edge
lengths are not corrected because it makes no difference to the face flux
calculation.

## References

  * [1] R. L. Walko, R. Avissar, The Ocean-Land-Atmosphere Model (OLAM). Part II: Formulation and tests of the nonhydrostatic dynamic core, Mon. Wea. Rev. 136 (11) (2008) 4045–4062.  doi:10.1175/2008MWR2523.1.
  * [2] C. Schär, D. Leuenberger, O. Fuhrer, D. Lüthi, C. Girard, A new terrain-following vertical coordinate formulation for atmospheric prediction models, Mon. Wea. Rev. 130 (10) (2002) 2459–2480.  doi:10.1175/1520-0493(2002)130<2459:ANTFVC>2.0.CO;2.
  * [3] K. P. Hoinka, G. Zängl, The influence of the vertical coordinate on simulations of a pv streamer crossing the alps, Mon. Wea. Rev. 132 (7) (2004) 1860–1867.  doi:10.1175/1520-0493(2004)132<1860:TIOTVC>2.0.CO;2.
  * [4] S. Webster, A. Brown, D. Cameron, C. Jones, Improvements to the representation of orography in the Met Office Unified Model, Quart. J. Roy. Meteor. Soc. 129 (591) (2003) 1989–2010.  doi:10.1256/qj.02.133.
  * [5] T. Gal-Chen, R. C. Somerville, On the use of a coordinate transformation for the solution of the Navier-Stokes equations, J. Comp. Phys. 17 (2) (1975) 209–228.  doi:10.1016/0021-9991(75)90037-6.
  * [6] J. B. Klemp, A terrain-following coordinate with smoothed coordinate surfaces, Mon. Wea. Rev. 139 (7) (2011) 2163–2169.  doi:10.1175/MWR-D-10-05046.1.
  * [7] S. D. Eckermann, J. P. McCormack, J. Ma, T. F. Hogan, K. A. Zawdie, Stratospheric analysis and forecast errors using hybrid and sigma coordinates, Mon. Wea. Rev. 142 (1) (2014) 476–485.  doi:10.1175/MWR-D-13-00203.1.
  * [8] A. J. Simmons, D. M. Burridge, An energy and angular-momentum conserving vertical finite-difference scheme and hybrid vertical coordinates, Mon. Wea. Rev. 109 (4) (1981) 758–766.  doi:10.1175/1520-0493(1981)109<0758:AEAAMC>2.0.CO;2.
  * [9] D. Leuenberger, M. Koller, O. Fuhrer, C. Schär, A generalization of the SLEVE vertical coordinate, Mon. Wea. Rev. 138 (9) (2010) 3683–3689.  doi:10.1175/2010MWR3307.1.
  * [10] H. Yamazaki, T. Satomura, N. Nikiforakis, Three-dimensional cut-cell modelling for high-resolution atmospheric simulations, Quart. J. Roy. Meteor. Soc. 142 (696) (2016) 1335–1350.  doi:10.1002/qj.2736.
  * [11] J. Steppeler, H.-W. Bitzer, M. Minotte, L. Bonaventura, Nonhydrostatic atmospheric modeling using a $z$-coordinate representation, Mon. Wea. Rev. 130 (8) (2002) 2143–2149.  doi:10.1175/1520-0493(2002)130<2143:NAMUAZ>2.0.CO;2.
  * [12] S. Jebens, O. Knoth, R. Weiner, Partially implicit peer methods for the compressible Euler equations, J. Comp. Phys. 230 (12) (2011) 4955–4974.  doi:10.1016/j.jcp.2011.03.015.
  * [13] J. Shaw, H. Weller, Comparison of terrain following and cut cell grids using a non-hydrostatic model, Mon. Wea. Rev. 144 (6) (2016) 2085–2099.  doi:10.1175/MWR-D-15-0226.1.
  * [14] B. Leonard, M. MacVean, A. Lock, Positivity-preserving numerical schemes for multidimensional advection, Tech. Rep. 106055, NASA (1993).
  * [15] J. Kent, P. A. Ullrich, C. Jablonowski, Dynamical core model intercomparison project: Tracer transport test cases, Quart. J. Roy. Meteor. Soc. 140 (681) (2014) 1279–1293.  doi:10.1002/qj.2208.
  * [16] B. Leonard, A. Lock, M. MacVean, Conservative explicit unrestricted-time-step multidimensional constancy-preserving advection schemes, Mon. Wea. Rev. 124 (11) (1996) 2588–2606.  doi:10.1175/1520-0493(1996)124<2588:CEUTSM>2.0.CO;2.
  * [17] Y. Chen, H. Weller, S. Pring, J. Shaw, Dimension splitting errors and a long time-step multi-dimensional scheme for atmospheric transport, Quart. J. Roy. Meteor. Soc.In revision; preprint, http://arxiv.org/abs/1701.06907.
  * [18] S.-J. Lin, R. B. Rood, Multidimensional flux-form semi-Lagrangian transport schemes, Mon. Wea. Rev. 124 (9) (1996) 2046–2070.  doi:10.1175/1520-0493(1996)124<2046:MFFSLT>2.0.CO;2.
  * [19] W. Guo, R. D. Nair, J.-M. Qiu, A conservative semi-Lagrangian discontinuous Galerkin scheme on the cubed sphere, Mon. Wea. Rev. 142 (1) (2014) 457–475.  doi:10.1175/MWR-D-13-00048.1.
  * [20] W. M. Putman, S.-J. Lin, Finite-volume transport on various cubed-sphere grids, J. Comp. Phys. 227 (1) (2007) 55–78.  doi:10.1016/j.jcp.2007.07.022.
  * [21] A. Iske, M. Käser, Conservative semi-Lagrangian advection on adaptive unstructured meshes, Numer. Methods. Partial Differ. Equ. 20 (3) (2004) 388–411.  doi:10.1002/num.10100.
  * [22] P. H. Lauritzen, R. D. Nair, P. A. Ullrich, A conservative semi-lagrangian multi-tracer transport scheme (CSLAM) on the cubed-sphere grid, J. Comp. Phys. 229 (5) (2010) 1401–1424.  doi:10.1016/j.jcp.2009.10.036.
  * [23] P. H. Lauritzen, C. Jablonowski, M. A. Taylor, R. D. Nair, Numerical techniques for global atmospheric models, Vol. 80, Springer Science & Business Media, 2011\.  doi:10.1007/978-3-642-11640-7.
  * [24] M. Zerroukat, T. Allen, A three-dimensional monotone and conservative semi-Lagrangian scheme (SLICE-3D) for transport problems, Quart. J. Roy. Meteor. Soc. 138 (667) (2012) 1640–1651.  doi:10.1002/qj.1902.
  * [25] R. K. Lashley, Automatic generation of accurate advection schemes on unstructured grids and their application to meteorological problems, Ph.D. thesis, University of Reading (2002).
  * [26] W. C. Skamarock, M. Menchaca, Conservative transport schemes for spherical geodesic grids: High-order reconstructions for forward-in-time schemes, Mon. Wea. Rev. 138 (12) (2010) 4497–4508.  doi:10.1175/2010MWR3390.1.
  * [27] P. H. Lauritzen, C. Erath, R. Mittal, On simplifying ‘incremental remap’–based transport schemes, J. Comp. Phys. 230 (22) (2011) 7957–7963.  doi:10.1016/j.jcp.2011.06.030.
  * [28] J. Thuburn, C. Cotter, T. Dubos, A mimetic, semi-implicit, forward-in-time, finite volume shallow water model: comparison of hexagonal-icosahedral and cubed-sphere grids, Geosci. Model Dev. 7 (3) (2014) 909–929.  doi:10.5194/gmd-7-909-2014.
  * [29] H. Weller, H. G. Weller, A. Fournier, Voronoi, Delaunay, and block-structured mesh refinement for solution of the shallow-water equations on the sphere, Mon. Wea. Rev. 137 (12) (2009) 4208–4224.  doi:10.1175/2009MWR2917.1.
  * [30] W. C. Skamarock, A. Gassmann, Conservative transport schemes for spherical geodesic grids: High-order flux operators for ODE-based time integration, Mon. Wea. Rev. 139 (9) (2011) 2962–2975.  doi:10.1175/MWR-D-10-05056.1.
  * [31] J. Steppeler, J. Klemp, Advection on cut-cell grids for an idealized mountain of constant slope, Mon. Wea. Rev.doi:10.1175/MWR-D-16-0308.1.
  * [32] P. Lauritzen, P. Ullrich, C. Jablonowski, P. Bosler, D. Calhoun, A. Conley, T. Enomoto, L. Dong, S. Dubey, O. Guba, et al., A standard test case suite for two-dimensional linear transport on the sphere: results from a collection of state-of-the-art schemes, Geosci. Model Dev. 7 (1) (2014) 105–145.  doi:10.5194/gmd-7-105-2014.
  * [33] S. Ii, F. Xiao, High order multi-moment constrained finite volume method. Part I: Basic formulation, J. Comp. Phys. 228 (10) (2009) 3669–3707.  doi:10.1016/j.jcp.2009.02.009.
  * [34] X. Li, C. Chen, X. Shen, F. Xiao, A multimoment constrained finite-volume model for nonhydrostatic atmospheric dynamics, Mon. Wea. Rev. 141 (4) (2013) 1216–1240.  doi:10.1175/MWR-D-12-00144.1.
  * [35] B. Xie, F. Xiao, A multi-moment constrained finite volume method on arbitrary unstructured grids for incompressible flows, J. Comp. Phys. 327 (2016) 747–778.
  * [36] P. K. Smolarkiewicz, J. Szmelter, A nonhydrostatic unstructured-mesh soundproof model for simulation of internal gravity waves, Acta Geophys. 59 (6) (2011) 1109–1134.  doi:10.2478/s11600-011-0043-z.
  * [37] P. K. Smolarkiewicz, J. Szmelter, MPDATA: An edge-based unstructured-grid formulation, J. Comp. Phys. 206 (2) (2005) 624–649.  doi:10.1016/j.jcp.2004.12.021.
  * [38] C. Kühnlein, P. K. Smolarkiewicz, An unstructured-mesh finite-volume MPDATA for compressible atmospheric dynamics, J. Comp. Phys. 334 (2017) 16–30.  doi:10.1016/j.jcp.2016.12.054.
  * [39] H. Weller, A. Shahrokhi, Curl free pressure gradients over orography in a solution of the fully compressible Euler equations with implicit treatment of acoustic and gravity waves, Mon. Wea. Rev. 142 (12) (2014) 4439–4457.  doi:10.1175/MWR-D-14-00054.1.
  * [40] L. Cueto-Felgueroso, I. Colominas, J. Fe, F. Navarrina, M. Casteleiro, High-order finite volume schemes on unstructured grids using moving least-squares reconstruction. Application to shallow water dynamics, Int. J. Numer. Meth. Engng 65 (3) (2006) 295–331.  doi:10.1002/nme.1442.
  * [41] L. Cueto-Felgueroso, I. Colominas, X. Nogueira, F. Navarrina, M. Casteleiro, Finite volume solvers and moving least-squares approximations for the compressible Navier–Stokes equations on unstructured grids, Comput. Methods Appl. Mech. Engrg 196 (45) (2007) 4712–4736.  doi:10.1016/j.cma.2007.06.003.
  * [42] R. P. Laurent White, D. Trenev, Flow simulation in heterogeneous porous media with the moving least-squares method, SIAM J. Sci. Comput. 39 (2) (2017) B323–B351.  doi:10.1137/16M1070840.
  * [43] P. H. Lauritzen, W. C. Skamarock, M. Prather, M. Taylor, A standard test case suite for two-dimensional linear transport on the sphere, Geosci. Model Dev. 5 (3) (2012) 887–901.  doi:10.5194/gmd-5-887-2012.
  * [44] D. R. Durran, Numerical methods for wave equations in geophysical fluid dynamics, Vol. 32, Springer Science & Business Media, 2013.  doi:10.1007/978-1-4419-6412-0.
  * [45] CFD Direct, OpenFOAM user guide: Numerical schemes, http://cfd.direct/openfoam/user-guide/fvschemes/ (2016).
  * [46] H. Weller, J. Shaw, AtmosFOAM: OpenFOAM applications and libraries for atmospheric modelling (2017).  doi:10.5281/zenodo.247327.
  * [47] Leibniz Institute for Tropospheric Research, J. Shaw, ASAM cut cell grid generator (2017).  doi:10.5281/zenodo.242374.
  * [48] J. Shaw, GMV-to-OpenFOAM file format converter (2017).  doi:10.5281/zenodo.242387.
  * [49] J. Thuburn, C. J. Cotter, T. Dubos, H. Weller, J. Shaw, Hexagonal/triangular geodesic mesh generator (2017).  doi:10.5281/zenodo.245327.
  * [50] J. Shaw, AtmosTests: Idealised atmospheric test suite (2017).  doi:10.5281/zenodo.247480.
  * [51] J. Shaw, Test result data for idealised numerical atmospheric transport (2017).  doi:10.5281/zenodo.259328.
  * [52] B. Good, A. Gadian, S.-J. Lock, A. Ross, Performance of the cut-cell method of representing orography in idealized simulations, Atmos. Sci. Lett. 15 (1) (2014) 44–49.  doi:10.1002/asl2.465.
  * [53] M. Jähn, O. Knoth, M. König, U. Vogelsberg, ASAM v2.7: a compressible atmospheric model with a Cartesian cut cell approach, Geosci. Model Dev. 8 (2) (2015) 317–340.  doi:10.5194/gmd-8-317-2015.
  * [54] U. Vogelsberg, ASAM wiki: Grid generator, http://asamwiki.tropos.de/index.php?title=Grid_Generator (2010).
  * [55] M. Baldauf, Stability analysis for linear discretisations of the advection equation with Runge–Kutta time integration, J. Comp. Phys. 227 (13) (2008) 6638–6659.  doi:10.1016/j.jcp.2008.03.025.
  * [56] R. Heikes, D. A. Randall, Numerical integration of the shallow-water equations on a twisted icosahedral grid. Part I: Basic design and results of tests, Mon. Wea. Rev. 123 (6) (1995) 1862–1880.  doi:10.1175/1520-0493(1995)123<1862:NIOTSW>2.0.CO;2.
  * [57] R. Heikes, D. A. Randall, Numerical integration of the shallow-water equations on a twisted icosahedral grid. Part II: A detailed description of the grid and an analysis of numerical accuracy, Mon. Wea. Rev. 123 (6) (1995) 1881–1887.  doi:10.1175/1520-0493(1995)123<1881:NIOTSW>2.0.CO;2.
  * [58] A. Staniforth, J. Thuburn, Horizontal grids for global weather and climate prediction models: a review, Quart. J. Roy. Meteor. Soc. 138 (662) (2012) 1–26.  doi:10.1002/qj.958.
  * [59] J. White, J. J. Dongarra, High-performance high-resolution semi-Lagrangian tracer transport on a sphere, J. Comp. Phys. 230 (17) (2011) 6778–6799.  doi:10.1016/j.jcp.2011.05.008.
  * [60] H. Miura, An upwind-biased conservative advection scheme for spherical hexagonal-pentagonal grids, Mon. Wea. Rev. 135 (12) (2007) 4038–4044.  doi:10.1175/2007MWR2101.1.
  * [61] G. Van Brummelen, Heavenly mathematics: The forgotten art of spherical trigonometry, Princeton University Press, 2013.