**Development and technical paper**| 16 Aug 2022

# Islet: interpolation semi-Lagrangian element-based transport

Andrew M. Bradley, Peter A. Bosler, and Oksana Guba

**Andrew M. Bradley et al.**Andrew M. Bradley, Peter A. Bosler, and Oksana Guba

- Sandia National Laboratories, Albuquerque, New Mexico, USA

- Sandia National Laboratories, Albuquerque, New Mexico, USA

**Correspondence**: Andrew M. Bradley (ambradl@sandia.gov)

**Correspondence**: Andrew M. Bradley (ambradl@sandia.gov)

Received: 27 Aug 2021 – Discussion started: 27 Sep 2021 – Revised: 19 Jul 2022 – Accepted: 19 Jul 2022 – Published: 16 Aug 2022

Advection of trace species, or *tracers*, also called tracer transport,
in models of the atmosphere and other physical domains
is an important and potentially computationally expensive part of a model's dynamical core.
Semi-Lagrangian (SL) advection methods are efficient because
they permit a time step much larger than the advective stability limit
for explicit Eulerian methods
without requiring the solution of a globally coupled system of equations as implicit Eulerian methods do.
Thus, to reduce the computational expense of tracer transport,
dynamical cores often use SL methods to advect tracers.
The class of interpolation semi-Lagrangian (ISL) methods contains potentially extremely efficient SL methods.
We describe a finite-element ISL transport method that we call the
interpolation semi-Lagrangian element-based transport (Islet) method,
such as for use with atmosphere models discretized using the spectral element method.
The *Islet method* uses three grids that share an element grid:
a dynamics grid supporting, for example, the Gauss–Legendre–Lobatto basis of degree three;
a physics parameterizations grid with a configurable number of finite-volume subcells per element;
and a tracer grid supporting use of *Islet bases* with particular basis again configurable.
This method provides extremely accurate tracer transport
and excellent diagnostic values in a number of verification problems.

An atmosphere model has two parts:
the dynamical core, which computes resolved fluid flow and resolved thermodynamics (Lauritzen, 2011);
and the subgrid parameterizations, which compute unresolved chemistry and physics processes
(Stensrud, 2009).
In turn, the dynamical core also has two parts:
first, the dynamics solver, which solves the equations of fluid motion;
second, the tracer transport solver, which advects trace atmosphere species, or *tracers*,
using the air density and flow fields from the dynamics solver.

Because of the large number of tracers in climate models, tracer transport can be computationally very expensive. For example, in the Dept. of Energy's Energy Exascale Earth System Model (E3SM) (E3SM Project, 2018) Atmosphere Model (EAM) version 1 (Golaz et al., 2019), configured with the default 40 tracers, tracer transport takes approximately 75 % of the total dynamical core wall clock time and approximately 23 % of the total atmosphere model wall clock time on a typical computer cluster (Golaz et al., 2022, Fig. 3). For EAM version 2, we developed a new tracer transport method that is 6.5 to over 8 times faster than in EAM version 1 in the cases of, respectively, low and high workload per computer node (Golaz et al., 2022, Fig. 3). In addition, we developed remap operators to remap data between separate grids for physics parameterizations and dynamics, permitting the physics parameterization computations to run on a coarser grid and thus 1.6 to 2.2 times faster in version 2 than in version 1 (Hannah et al., 2021; Golaz et al., 2022).

In this article, we describe the interpolation semi-Lagrangian element-based transport method, with the acronym stylized as “Islet” rather than “ISLET” to minimize distracting all-capitalized text. The Islet method extends those we developed for EAM version 2 by providing the discretizations with natural parameters to trade between computational cost and accuracy in tracer transport and physics parameterizations for a given, fixed dynamical core configuration.

## 1.1 Problem statement

The mass continuity equation for the air density *ρ* is

where ** u** is the flow velocity.
The tracer transport equation in continuity form and with a source term for a tracer mixing ratio

*q*and corresponding tracer density

*ρ*

*q*is

where *f* is a source term for *q*.
Combining Eqs. (1) and (2) gives
the advective form of the tracer transport equation,

Because of time splitting in the atmosphere model, our focus is often on the sourceless advection equation,

At position ** x** and time

*t*

_{1}, the exact solution of Eq. (4) in terms of the solution at another time

*t*

_{0}is

$\mathit{X}({t}_{\mathrm{0}};\mathit{x},{t}_{\mathrm{1}})$ is the solution of the ordinary differential equation

with the initial condition $\mathit{X}({t}_{\mathrm{1}};\mathit{x},{t}_{\mathrm{1}})=\mathit{x}$.
In a semi-Lagrangian method,
$\mathit{X}({t}_{\mathrm{1}};\mathit{x},{t}_{\mathrm{1}})=\mathit{x}$ is called the *arrival* point
and $\mathit{X}({t}_{\mathrm{0}};\mathit{x},{t}_{\mathrm{1}})$ is called the *departure* point.
In this article, our focus is on the two-dimensional equations on the sphere.
Thus, ** X** and

**each have two horizontal components and no vertical component. In various settings, the density variable**

*u**ρ*can be, for example, mass per unit area, layer thickness, or layer pressure increment.

An approximate numerical solution for *q* of Eq. (4)
is said to be *property preserving* if
(possibly just a subset of) properties that hold for the exact solution
also hold for the approximate one.
Eq. (5) implies that advection cannot introduce new extrema in the mixing ratio;
advection is said to be *shape preserving*.
Equation (2) with *f*=0 implies the global mass is conserved.
Although the focus of this article is not the continuity equation,
we note that the Lagrangian form of the continuity equation,
Eq. (B4) in Appendix B,
implies that the total mass in a Lagrangian parcel,
which is a parcel of fluid that moves with the flow, is constant.
A final property that is a special case of the shape-preserving property
relates to coupling a solver for Eq. (4)
to a dynamics solver: mass-tracer consistency.
This property means that if *q* is constant in space at time *t*_{0},
then it remains constant in space at every other time.
In other words, the dynamics solver and transport solver use the same air density.
The methods in this article conserve global mass, do not introduce new nodal extrema,
and provide mass-tracer consistency when coupled to a dynamics solver.

There is a large body of work on tracer transport methods. As a result, researchers have developed a large number of accuracy, property preservation, and other diagnostics for use in comparisons of methods (Lauritzen et al., 2012, 2014; Lauritzen and Thuburn, 2012; Lauritzen, 2011; Lauritzen et al., 2015).

Another quality of a transport method is important: its computational efficiency. Computational efficiency is some measure of solution accuracy for a given set of computational resources.

Thus, when developing a new transport method, the objective is to obtain high efficiency, as measured by diagnostic values and computational cost, constrained by the need to couple to specific dynamics and physics grids. Our objective in this article is to extend our highly efficient tracer transport method for EAM version 2 by introducing, first, parameters whose values trade between solution accuracy and computational cost and, second, algorithms associated with these parameters.

Our motivation is to support the computationally efficient simulation of strongly tracer-dependent models, such as those of aerosols, by enabling extremely highly resolved tracer filamentary structure for a fixed dynamics resolution. A recent report from the National Academies of Sciences, Engineering, and Medicine (Field et al., 2021) includes the recommendation to “explore whether global aerosol optical depth (AOD) distribution is significantly affected by plume-scale effects,” and asks: “Are nested grids needed to represent plume processes? What spatial resolution is needed to faithfully represent the radiative forcing and impact outcomes?” In this article, we present algorithms that can help to address these questions.

## 1.2 Overview of the E3SM Atmosphere Model

EAM uses the High-Order Method Modeling Environment (HOMME) dynamical core (Dennis et al., 2005, 2012).
HOMME uses the spectral element method to discretize the equations of atmospheric flow.
HOMME's grid consists of
horizontally unstructured quadrilateral elements
extruded as columns in the vertical dimension.
*Element-based* methods permit extremely flexible discretization of domains:
for example, E3SM's regionally refined model (RRM) configurations (Tang et al., 2019),
in which the atmosphere element grid is refined in a particular region of the earth.

In the horizontal direction, each quadrilateral element has an *n*_{p}×*n*_{p} subgrid.
HOMME uses the Gauss–Legendre–Lobatto (GLL) spectral element polynomial basis and nodes.
The two-dimensional (2D) basis functions are tensor products
of one-dimensional (1D) GLL basis functions.
*n*_{p} is the number of subgrid nodes in a 1D basis function;
the degree of the polynomial is *n*_{p}−1.
Figure 1 shows one spectral element outlined in blue.
The large black dots are the 4×4 tensor-product grid of GLL nodes with *n*_{p}=4.
The Courant–Friedrichs–Lewy (CFL) condition for the advection operator in spectral element fluid dynamics scales inversely
with the minimum distance between GLL nodes.
In turn, this minimum distance scales roughly quadratically in *n*_{p}.
For example, nodes for *n*_{p}=6, 8, and 12 have minimum distances
respectively 2.4, 4.3, and 10.0 times smaller than nodes for *n*_{p}=4.
Thus, developers of a dynamics solver choose a value for *n*_{p} that is large enough
to provide good accuracy and wave representation
but small enough to permit a large time step,
i.e., one determined by accuracy rather than stability considerations.
EAM always sets ${n}_{\mathrm{p}}^{\mathrm{v}}=\mathrm{4}$, where we use v to refer specifically to the dynamics solver's value of *n*_{p}.
We refer to the spectral element grid, including the spectral element nodes,
as the *dynamics grid*,
and the element grid excluding subgrid nodes as the *element grid*.

For efficiency, the EAM version 2 physics and chemistry parameterizations run
on a different grid than the dynamics;
we call this grid the *physics grid*.
To couple efficiently to the dynamics grid,
the physics grid has the same element grid as the dynamics grid
but its own subelement grid (Herrington et al., 2019; Hannah et al., 2021).
The physics grid divides each horizontal quadrilateral element into
a regular *n*_{f}×*n*_{f} subgrid of finite-volume quadrilateral subcells.
In Fig. 1,
the dashed green lines outline the 6×6 physics subgrid, with *n*_{f}=6.
EAM version 2 uses *n*_{f}=2.
Compared with EAM version 1, in which the physics and dynamics grids were the same,
EAM version 2 has $\mathrm{4}/\mathrm{9}$ as many physics grid points,
speeding up the parameterizations part of the model by up to 2.25 times
with a small cost for remapping between dynamics and physics grids.
The shared element grid and the element-local remap operators
minimize interprocess communication costs.

## 1.3 Overview of the Islet method

While our tracer transport method is based on spectral elements
and, like the physics grid, uses the same element grid as the dynamics grid,
it is not restricted by the CFL condition.
This important fact motivates making ${n}_{\mathrm{p}}^{\mathrm{t}}$,
where we use t to refer specifically to the transport solver's value of *n*_{p},
a parameter of the tracer transport module and independent of ${n}_{\mathrm{p}}^{\mathrm{v}}$,
to trade between computational speed (maximum at ${n}_{\mathrm{p}}^{\mathrm{t}}={n}_{\mathrm{p}}^{\mathrm{v}}$)
and higher accuracy (increasing with ${n}_{\mathrm{p}}^{\mathrm{t}}$).
In Fig. 1,
the small red circles are GLL nodes for ${n}_{\mathrm{p}}^{\mathrm{t}}=\mathrm{8}$.
We refer to this approach as *tracer transport p-refinement*.
In the finite element method, *p*-refinement means increasing the basis polynomial degree.
In the Islet method, we increase ${n}_{\mathrm{p}}^{\mathrm{t}}$ relative to ${n}_{\mathrm{p}}^{\mathrm{v}}$ to represent the mixing ratio fields
at higher resolution.

In summary, we use one element grid, three subelement grids, and three corresponding parameters:
the dynamics grid, ${n}_{\mathrm{p}}^{\mathrm{v}}$; the physics grid, *n*_{f}; and the tracer transport grid, ${n}_{\mathrm{p}}^{\mathrm{t}}$.
For efficiency, EAM version 2 already uses two grids, one with ${n}_{\mathrm{p}}^{\mathrm{t}}={n}_{\mathrm{p}}^{\mathrm{v}}=\mathrm{4}$ and the other with *n*_{f}=2.
This article considers the case ${n}_{\mathrm{p}}^{\mathrm{t}}>{n}_{\mathrm{p}}^{\mathrm{v}}$.
Figure 1 shows the three subgrids for one spectral element.

The tracer transport solver couples to the dynamics solver through the flow velocity ** u**,
the air density

*ρ*, and possibly a small subset of tracers such as specific humidity. It couples to the parameterizations through the tracers. Coupling requires remapping a field between grids. It is natural to give responsibility for these remap procedures to the tracer transport solver module as the remap operators share many of the same requirements and computations as the transport solver. As in the case of remapping quantities between the physics grid and the dynamics grid in EAM version 2, the shared element grid and element-local remap procedures minimize communication when remapping among the three grids. For brevity, we refer to this overall spectral element, three-grid scheme as the

*Islet method*.

## 1.4 Outline

The rest of this article is structured as follows. The Islet method depends on a number of lower-level algorithms. Section 2 and Appendix A describe these; mathematical details can be skipped on a first reading. Section 3 describes each step of the Islet method over the course of one time step. Descriptions are compact, relying on the details in Sect. 2. Section 4 describes and presents the results of a number of verification problems. Finally, Sect. 5 concludes.

This section describes the details of low-level algorithms on which the Islet method depends, associated concepts, and related work for context. Our objective is to isolate the notation for and details of these low-level algorithms to this section, then focus on the higher-level details of the Islet method in Sect. 3 without lower-level clutter. Section 2.1 discusses semi-Lagrangian transport, the finite-element interpolation semi-Lagrangian method in particular, and the bases that the Islet method uses. Section 2.2 describes our approach to solving the property preservation problem. Section 2.3 describes algorithms to remap data between grids. Section 2.4 describes the spectral element direct stiffness summation operator.

## 2.1 Semi-Lagrangian transport

The Islet method uses semi-Lagrangian (SL) transport.
In an SL method,
following Eqs. (5) and (6),
each Eulerian arrival grid point at time *t*_{1}
is advected backward or forward in time to its
generally off-grid departure point at time *t*_{0}.
The resulting Lagrangian departure grid then samples the field at time *t*_{0}
and, by some means, uses these data to construct the advected field at time *t*_{1}.

In the problem of tracer transport, for which the flow velocity ** u** is provided,
SL methods can take much longer time steps than Eulerian methods.
Data dependencies are local to a small domain of dependence;
thus, unlike an implicit method,
SL tracer transport does not require the solution of a global linear equation.
This long time step requires
many fewer interprocess communication rounds
and sometimes fewer floating point operations
per unit of simulation time
than Eulerian methods,
making SL methods substantially more efficient than Eulerian methods for tracer transport.

Transport methods often are the composition of a number of operators. Two key operators are the linear advection operator and the nonlinear property preservation operator; we describe ours in Sect. 2.1.4 and 2.2, respectively.

### 2.1.1 Types of SL methods

For a broad survey of atmosphere tracer transport methods, including a large number of Lagrangian and semi-Lagrangian methods, we refer the reader to Lauritzen et al. (2014, Sect. 2). In addition, Giraldo et al. (2003), Natarajan and Jacobs (2020), and Bradley et al. (2019) review types of SL methods. This article focuses on remap-form interpolation methods.

A method in *remap form* directly remaps
the tracer field on the Eulerian grid at the previous time step to
the Lagrangian grid,
in contrast to flux-form methods that integrate over swept regions (Lauritzen et al., 2010; Lee et al., 2016)
or characteristic curves (Erath and Nair, 2014).
The computational and communication costs of a remap-form method
are very nearly independent of time step
whereas a flux-form method's cost grows roughly linearly with the time step.

An interpolation SL (ISL) method directly discretizes the tracer transport equation,
Eq. (5):
$q({\mathit{x}}_{i},{t}_{\mathrm{1}})=q({\overline{\mathit{x}}}_{i},{t}_{\mathrm{0}})$,
where *x*_{i} is an Eulerian arrival grid point and
${\overline{\mathit{x}}}_{i}$ is its departure, generally off-grid, point at *t*_{0}<*t*_{1}.
*q* has exact values only at grid points *x*_{i};
thus, evaluating $q({\overline{\mathit{x}}}_{i},{t}_{\mathrm{0}})$ requires interpolation.

Interpolation is in contrast to *exactly cell-integrated* methods,
which accurately integrate the basis of a target (e.g., Lagrangian) element against those of the source;
see, e.g., Bosler et al. (2019).
(In some cases, an inaccurate cell-integrated method can be interpreted as an interpolation method;
see Appendix B for an example.)
Exactly cell-integrated methods have substantially greater cost than interpolation methods for three reasons.

First, to obtain smoothness in the integrand, integration is over facets computed by geometric intersection of a target element against source elements; intersection calculations are not needed in interpolation methods. Typically, to minimize computational geometry complexity, departure cell edges are approximated by great arcs rather than flow-distorted curves, limiting the method to second-order accuracy; however, Ullrich et al. (2013) describe a higher-order edge reconstruction that yields a third-order accurate advection method. In contrast, achieving arbitrarily high order in an ISL method's linear advection operator does not entail any additional complexity.

Second, accurate integration has a larger computational cost because it requires sphere-to-reference point calculation and interpolant evaluations at many quadrature points.

Third, an exactly cell-integrated method requires a larger communication volume because all data from a source element are used to integrate against each target basis function.

In trade for these additional costs,
exactly cell-integrated methods are locally mass conserving,
and the fact that they are *L*^{2} projections can be used to prove stability.
Local mass conservation means that
one can identify numerical, possibly Lagrangian, fluid parcels on the grid
that have constant tracer mass.
*Local* is in contrast to *global* mass conservation;
the latter means that the mass of the tracer fluid is conserved over the whole domain
but not necessarily in any identifiable parcels smaller than the domain.
Although an exactly cell-integrated method is locally mass conserving,
coupling it to a dynamics solver
still generally requires additional measures to obtain mass-tracer consistency.

### 2.1.2 Spectral element ISL transport

Our objective in this work is to use the design freedom
provided by giving up local, but not global, mass conservation
to maximize computational efficiency.
We develop an ISL method that uses
the *Islet bases*,
summarized in Appendix A
and derived in Bradley et al. (2021, Sects. 2 and 3)
and a forthcoming article based on that material,
to satisfy a necessary condition for stability.

Figure 2a illustrates the
linear advection step in the Islet method.
First, the upper-left element is advected backward in time from time *t*_{1} to time *t*_{0},
using velocity data on the ${n}_{\mathrm{p}}^{\mathrm{v}}=\mathrm{4}$ dynamics grid to advect each dynamics grid point.
Second, the ${n}_{\mathrm{p}}^{\mathrm{v}}=\mathrm{4}$ isoparametric map interpolates
the locations of the Lagrangian ${n}_{\mathrm{p}}^{\mathrm{t}}=\mathrm{6}$ tracer grid points
from the Lagrangian ${n}_{\mathrm{p}}^{\mathrm{v}}=\mathrm{4}$ dynamics grid points.
Third, the tracer data at time *t*_{0} in the bottom-right Eulerian element
are used to interpolate the tracer mixing ratio at *t*_{0} to each green starred Lagrangian point.
Interpolation uses the 2D tensor-product ${n}_{\mathrm{p}}^{\mathrm{t}}=\mathrm{6}$ nodal basis.
Fourth, these values are then the values in the upper-left Eulerian element at time *t*_{1}.

Figure 2b shows two bases,
the GLL basis, called the *natural* basis when a distinguishing name is needed,
and the Islet GLL basis.
The ISL method using the natural GLL basis is unstable
for the test problem of uniform flow on a uniform grid,
but the ISL method using the Islet GLL basis is stable for this problem (Bradley et al., 2021).
In this article, all of the algorithms accommodate either basis type,
and we use the Islet GLL bases to obtain stability;
see Appendix A for details.

Each basis is a *nodal* basis:
a basis function has value 1 at one node and 0 at every other node.
Thus, each basis function is associated with a node.
For example, in Fig. 2b,
the blue basis function is associated with the third node of six.
The basis is symmetric;
basis function $k\in \mathit{\{}\mathrm{0},\mathrm{\dots},{n}_{\mathrm{p}}-\mathrm{1}\mathit{\}}$ is the mirror image of basis function ${n}_{\mathrm{p}}-k-\mathrm{1}$.
Thus, the blue and cyan functions are mirror images around reference coordinate 0.

There are two conventional ways to enumerate
finite element nodal degrees of freedom and basis functions,
where for a scalar field, one degree of freedom is associated with each basis function.
First, in this article, by a basis we always mean an *element* basis
unless stated otherwise.
A basis function in an element basis set has non-zero values only within an element.
As a result,
nodes in the interior of an element have single-valued field values
while nodes at an element boundary may have multi-valued field values.

Second, the viewpoint of *global* continuous degrees of freedom and *global* basis functions
always assigns only one degree of freedom of the field to each grid point.
The basis function associated with a grid point in the interior of an element
is the same as in the element basis.
A grid point on the element boundary has a basis function
that spans all the elements that share that grid point.
This global basis function is the sum of the element basis functions associated with that grid point.

In the case of continuous finite element methods, the two viewpoints are equivalent, and convenience determines which viewpoint is adopted. In the case of a discontinuous finite element method, only the element viewpoint is possible.

The Islet method accommodates both continuous and discontinuous variants;
in this article, we present only the continuous variant.
We always take the element viewpoint,
regardless of the continuity of a field.
Thus, for a scalar field,
there is one degree of freedom associated with each grid point in an element
and, thus, *n* degrees of freedom at a grid point shared by *n* elements.
Index *i* is associated with grid point location *x*_{i}.
Again, to be clear, for a grid point at an element edge,
multiple index values are associated with the same location ** x**.
Some steps in the Islet method introduce discontinuities in the field across elements.
Section 2.4 describes the operator that restores continuity to a discontinuous field.

### 2.1.3 Definitions and notation

The tracer mixing ratio degrees of freedom are written *q*_{i} and
the air density, *ρ*_{i}.
The vector of all scalar values is bold-faced, e.g., ** q**.
A subset of these, with the index set given by ℰ,
is written

**(ℰ). When needed, a superscript letter indicates on which grid the field is approximated: t, tracer; v, dynamics; and f, physics.**

*q*Arithmetic between vectors applies entry-wise.
For two vectors, ** x** and

**, we write**

*y***⊙**

*x***for element-wise multiplication and**

*y***⊘**

*x***for element-wise division.**

*y*We need to define a map between
a reference element with coordinates $\mathit{r}\in [-\mathrm{1},\mathrm{1}{]}^{\mathrm{2}}$
and element *e* on the manifold with coordinates ** x**,
e.g., the sphere.
Let this map be

**=**

*x**m*

_{e}(

**) and its inverse be $\mathit{r}={m}_{e}^{-\mathrm{1}}\left(\mathit{x}\right)$. Additionally, let**

*r**e*=

*E*(

**) return the element index for the element that contains**

*x**e*. In this work, if

**is on the boundary of multiple elements,**

*x**E*(

**) can be any of the corresponding element indices without affecting the solution.**

*x*Similarly, we need to map between (element-)local and global indices.
Local indices are in the set $\mathit{\{}\mathrm{0},\mathrm{\dots},({n}_{\mathrm{p}}^{g}{)}^{\mathrm{2}}-\mathrm{1}\mathit{\}}$,
where $g\in \mathit{\{}\mathrm{v},\mathrm{t}\mathit{\}}$.
Global index *i* is in element $\mathrm{elm}(i;g)\equiv \lfloor i/({n}_{\mathrm{p}}^{g}{)}^{\mathrm{2}}\rfloor $.
The corresponding local index is $\mathrm{lcl}(i;g)\equiv i\mathrm{mod}({n}_{\mathrm{p}}^{g}{)}^{\mathrm{2}}$.
Given the local index *j* and element index *e*,
the global index is $\mathrm{glb}(j;g)\equiv ({n}_{\mathrm{p}}^{g}{)}^{\mathrm{2}}e+j$.
Let ℰ(*e*;*g*) be the set of global indices for the grid points in element *e*.
We omit *g* when it is not necessary to specify a grid.
For convenience, let *q*_{e}≡** q**(ℰ(

*e*)).

Element basis functions and reference-coordinate nodes are indexed locally.
Thus, global grid point *i* corresponds to
the 2D tensor-product element basis function *ϕ*_{lcl(i;g)} and
reference-coordinate node *r*_{lcl(i;g)}.

Occasionally, we write the tensor product in the 2D basis explicitly.
Let a 1D element basis function set be *ψ*_{j}, $j\in \mathrm{0},\mathrm{\dots},{n}_{\mathrm{p}}-\mathrm{1}$.
Let $\mathit{r}\equiv ({r}_{\mathrm{1}},{r}_{\mathrm{2}})$.
Then, ${\mathit{\varphi}}_{i}\left(\mathit{r}\right)\equiv {\mathit{\psi}}_{\lfloor i/{n}_{\mathrm{p}}\rfloor}\left({r}_{\mathrm{1}}\right){\mathit{\psi}}_{i\mathrm{mod}{n}_{\mathrm{p}}}\left({r}_{\mathrm{2}}\right)$.
Figure 2b shows the ${n}_{\mathrm{p}}^{\mathrm{t}}=\mathrm{6}$ basis functions *ψ*_{j}.

The function $\mathrm{normalize}\left(\mathit{x}\right)\equiv \mathit{x}/\Vert \mathit{x}{\Vert}_{\mathrm{2}}$.

Finally, we write the particular map ** x**=

*m*

_{e}(

**) we use. The map uses just the corner nodes of an element; it is the isoparametric map for a bilinear basis. In element**

*r**e*, let corner

_{e}(

*j*;v) be the global index of corner

*j*of the element, here referenced to the dynamics grid for concreteness. Let ${\mathit{\varphi}}_{j}^{\mathrm{bilin}}\left(\mathit{r}\right)$ be the bilinear basis functions, i.e., the GLL basis with

*n*

_{p}=2. Then,

This is the map described in Guba et al. (2014, Appendix A).

### 2.1.4 The linear advection operator

Now we write the steps illustrated in Fig. 2.
First, each grid point ${\mathit{x}}_{i}^{\mathrm{v}}$ is advected from time *t*_{1} to time *t*_{0}
to give

In this article, we assume this step is performed nearly exactly, i.e., with numerical error substantially smaller than any other source of error because this article focuses on the spatial part of the advection discretization. Second, the tracer grid points are computed from the dynamics grid points using the ${n}_{\mathrm{p}}^{\mathrm{v}}$-basis isoparametric map followed by normalization:

Third, an interpolated value from time *t*_{0} is assigned to each advected tracer grid point:

## 2.2 Property preservation

Shape preservation and mass conservation on their own are not necessarily difficult to achieve. Rather, the hard problem in the development of an efficient transport method is the consistent coupling of the transport scheme to the dynamical core, or mass-tracer consistency. Mass-tracer consistency means that the transport solver and dynamics solver agree exactly on the air density field. Achieving both mass-tracer consistency and a long transport time step requires redistributing tracer mass.

It is possible to couple a tracer transport method to a dynamical core consistently and preserve properties using only local operations. Lauritzen et al. (2017) describe the coupling of a finite-volume, cell-integrated SL method to the HOMME spectral element dynamical core using Lagrangian grid adjustments, where each adjustment is a function of only data in the adjusted entity's local neighborhood. However, in their method, the time step is limited to a Courant number of at most 1, limiting efficiency.

Our approach to consistent coupling is to use what is sometimes called a *global fixer*,
although, often, this term applies to only the problem of global mass conservation.
In Bradley et al. (2019), we describe several classes of methods that we call
Communication-Efficient Density Reconstructors (CEDR).
Each CEDR can solve a variety of property preservation problems,
including shape preservation and global mass conservation.
Each CEDR provides a tight and practically meaningful bound on mass redistribution.
A solution is assured if the input data meet necessary and sufficient conditions.
Depending on the algorithm, a CEDR uses either one or two global reductions or reduction-like communication rounds.
See Bradley et al. (2019, Sect. 1.2) for a discussion of how CEDRs relate to previous global fixers,
in particular, with respect to solution guarantees and number of communication rounds.

Because we use a CEDR for consistent coupling, we are then free to exploit it in other parts of our transport method to maximize efficiency. For example, we can use an interpolation SL method, which is not mass conserving, because the mass will be corrected at the end of each time step.

Exact shape preservation uses mixing ratio bounds computed from each arrival grid point's domain of dependence and does not permit any violation of a domain's nodal extrema. As a result, exact shape preservation clips smooth extrema, limiting the overall transport method to an order of accuracy generally not more than a little above two. Nonetheless, higher-order linear advection is still useful because it increases the absolute accuracy of the overall method even if not its order of accuracy, as examples in Sect. 4 will demonstrate.

### 2.2.1 ClipAndAssuredSum

In this article, we use the simplest CEDR: ClipAndAssuredSum, subsequently CAAS, algorithm 3.1 in Bradley et al. (2019).

To explain CAAS,
we need to define the weight associated with a grid point
since weights are needed to define the mass.
Conceptually, a weight is an integral over the element of the product of
the grid point's basis function and the Jacobian determinant of ** x**=

*m*(

**):**

*r*
Details of the discretization modify this value.
Our manifold is the sphere,
and we use HOMME's definition of a weight.
HOMME discretizes the weight integral using the *n*_{p}-basis GLL quadrature.
Let

Because *ϕ*_{lcl(i)}(*r*_{lcl(j)})=1 when *j*=*i* and is 0 for every other value of *j*,
GLL quadrature applied to Eq. (12) gives

In the final line, we define

and *w*_{i} as the product of the two 1D basis function integrals.
In practice, HOMME modifies each node's Jacobian determinant *J*_{i}
by multiplying it by a constant very close to 1 such that
the global sum ${\sum}_{i}{J}_{i}{w}_{i}=\mathrm{4}\mathit{\pi}$;
that is, the sum of all the products of *w*_{i} and adjusted *J*_{i} gives exactly the area of the unit sphere.
This is an optional and minor implementation detail,
and we neglect this scaling factor subsequently.
Note that, in general, ** J** is discontinuous across elements.

Now we describe CAAS
independent of two sets of details:
first, the means by which time-*t*_{1} data are obtained and,
second, the index set ℰ over which CAAS is applied.
Later, we will specify these details for each application of CAAS.
At time *t*_{1} and grid point *i*,
we are given air density *ρ*_{i},
preliminary tracer mixing ratio ${\overline{q}}_{i}$,
and weight *θ*_{i}.
In addition, we are given a target total mass *b* over all *i*∈ℰ.
Finally, we have lower and upper bounds, ${q}_{i}^{min}$ and ${q}_{i}^{max}$, on the target *q*_{i}.
These are the inputs to CAAS.
The output is a set of modified values ** q**(ℰ) that solve the 1-norm minimization problem

This problem has a solution, i.e., the constraint set is non-empty, if and only if

For details, see Bradley et al. (2019, Sect. 2).

When there is a solution, there are usually an infinite number,
and CAAS efficiently finds one as follows.
First, clip the mixing ratios for *i*∈ℰ:

Let the preliminary mass be

If $\widehat{m}=b$, then set ${q}_{i}={\widehat{q}}_{i}$ for *i*∈ℰ and terminate.
Otherwise, suppose $\widehat{m}<b$.
Second, compute the *capacity*:

Third, compute the final values for *i*∈ℰ:

As a check, note that

If $b-\widehat{m}\le c$, then ${q}_{i}\le {q}_{i}^{max}$ for *i*∈ℰ;
these conditions hold if and only if Eq. (14) does.
The case $\widehat{m}>b$ is similar.
Let

carry out these steps.

CAAS is a CEDR because the global version, i.e., the problem in which ℰ is over all grid points, can be implemented using a single global reduction. For details, see Bradley et al. (2019, Sect. 6).

### 2.2.2 Infeasible problems

In some situations, the constraint set is empty.
There is a sequence of up to two relaxations to the constraints to
produce a non-empty constraint set.
Algorithm 3.4 in Bradley et al. (2019), ReconstructSafely,
which wraps a lower-level algorithm like CAAS,
formalizes this sequence.
Suppose the original constraint set is empty and $\widehat{m}<b$.
ReconstructSafely first tries to compute a one-norm-minimal relaxation to *q*^{max},
such that the maximum value ${max}_{i\in \mathcal{E}}{q}_{i}^{max}$ is not exceeded:

This relaxed problem provides mass-tracer consistency and mass conservation
but does not solve the exact shape preservation problem.
If this relaxation has an empty constraint set,
then ReconstructSafely returns ** q**(ℰ)
with uniform values that satisfy the mass conservation constraint.

### 2.2.3 Local and global problems

We use CAAS to solve problems at two levels: first, within each element separately and, second, over the whole grid. In addition, at the global level, there are two natural choices for scalar values in the vectors passed to CAAS: grid point values and integrals over elements. In practice, some local problems may not have assuredly non-empty constraint sets, but the global problem always has a non-empty constraint set.

In the first approach to the global problem,
which we refer to as *CAAS–point*,
the inputs to CAAS are simply the grid point values.
Global and local applications of CAAS are then identical
except for the index set ℰ.

In the alternative approach,
which we refer to as *CAAS–CAAS*,
the global problem is solved over element-level scalar values,
and then a second step applies CAAS separately to each element.
The global step's task is to provide each element with a sufficient mass adjustment
so that each element's local CAAS step has a non-empty constraint set.
The global problem's inputs are formulated as follows.
Weight *e* is element *e*'s area: ${A}_{e}\equiv {\sum}_{i\in \mathcal{E}\left(e\right)}{w}_{i}{J}_{i}$.
Air density *ρ*_{e} is the element's average density.
Let the element's mass be
${m}_{e}\equiv {\sum}_{i\in \mathcal{E}\left(e\right)}{w}_{i}{J}_{i}{\mathit{\rho}}_{i}$.
Then, ${\mathit{\rho}}_{e}={m}_{e}/{A}_{e}$.
The minimum mixing ratio is the minimum tracer mass over the element
divided by the total mass:
${q}_{e}^{min}\equiv {\sum}_{i\in \mathcal{E}\left(e\right)}{w}_{i}{J}_{i}{\mathit{\rho}}_{i}{q}_{i}^{min}/{m}_{e}$,
and similarly for ${q}_{e}^{max}$.
Finally, the initial value ${\overline{q}}_{e}$ is
${\overline{q}}_{e}\equiv {\sum}_{i\in \mathcal{E}\left(e\right)}{w}_{i}{J}_{i}{\mathit{\rho}}_{i}{q}_{i}/{m}_{e}$.
This global CAAS step returns element-level mass adjustments.
The adjusted mass in an element is then *b* in the input to the local CAAS step.

### 2.2.4 Relaxed problems

In the Islet method,
multiple local and one global property preservation problems are solved in each time step.
Not every problem must enforce strict property preservation;
strictness is required only in certain outputs,
e.g., at the end of a time step.
Relaxing the extremal bounds on the mixing ratio grid point values
reduces the amount of mass redistribution.
More subtly, solving relaxed element-local problems before solving the strict global problem
increases intra-element and decreases inter-element mass redistribution,
thus, increasing the locality of mass redistribution.
Finally, problems such as the toy chemistry problem described in Sect. 4.3.2
that are sensitive to finite-precision round-off errors benefit from
relaxed bounds because the perturbation to the mixing ratio is decreased.
Thus, in some element-local CAAS applications in a time step,
for each index *i*, we widen each pair of bounds by 1 % of ${q}_{i}^{max}-{q}_{i}^{min}$ in each direction.
The global CAAS–point application uses strict bounds,
assuring the tracer field is property-preserving at the end of a time step.

## 2.3 Grid remap

Remapping data among multiple component grids is common in many applications of partial differential equation (PDE)-based modeling because it is a direct means to permit each component to run in its most efficient configuration. For example, many whole-earth, fully coupled earth system models use different grids for ocean, land, and atmosphere components. Remapping among subcomponent grids is less common. Recent examples include separate physics parameterizations and dynamics grids in the atmosphere (Herrington et al., 2019; Hannah et al., 2021); adaptive mesh refinement (AMR) of a tracer (Chen et al., 2021; Semakin and Rastigejev, 2020); and local vertical refinement in physics parameterizations relative to the shared background vertical grid, the Framework for Improvement by Vertical Enhancement (FIVE) (Yamaguchi et al., 2017; Lee et al., 2020).

In the Islet method, grid remap operators transfer data among the dynamics, tracer, and physics grids. In all of these grid transfers, there is one key property: the linear remap operators use only element-local data. Thus, remap between grids is extremely efficient.

Remapping a field between a GLL grid, dynamics or tracer, and the physics grid is described in detail in Hannah et al. (2021, Sect. 2), and these details are not reproduced here. Remapping from the dynamics grid to the physics grid is simple: a physics grid subcell's value is assigned the average value of the density over the subcell. Remapping from the physics grid to the dynamics grid requires a high-order reconstruction of the data on the physics grid in addition to maintaining an element-local mass-conservation constraint, and most of the details in Hannah et al. (2021, Sect. 2) focus on this high-order reconstruction.

Remapping a field between the dynamics and tracer grids in either direction is simple;
each grid hosts a field already in a high-order representation
and, thus, no high-order reconstruction is needed.
The first step is interpolation.
Let ℐ^{v→t} interpolate a field in an element on the dynamics grid
represented by the natural GLL ${n}_{\mathrm{p}}^{\mathrm{v}}$-basis
to one on the tracer grid
represented by the Islet ${n}_{\mathrm{p}}^{\mathrm{t}}$-basis,
using the natural ${n}_{\mathrm{p}}^{\mathrm{v}}$-basis functions to interpolate.
Let ℐ^{t→v} do the opposite
using the Islet ${n}_{\mathrm{p}}^{\mathrm{t}}$-basis functions.
We can write these linear operators explicitly using the notation
introduced in Sect. 2.1.2.
Let ** f** be the grid-point values of a scalar field.
Consider element

*e*. Then, for local index $i\in \mathit{\{}\mathrm{0},\mathrm{\dots},({n}_{\mathrm{p}}^{\mathrm{t}}{)}^{\mathrm{2}}-\mathrm{1}\mathit{\}}$,

where we have used the notation *f*_{e}≡** f**(ℰ(

*e*)). ℐ

^{t→v}is written the same way but with t and v switched.

The interpolation of a field *f* from the dynamics grid to the tracer grid
has the useful property that it is conservative within each element *e*
despite not being in all cases an *L*^{2} projection.
Consider ${\mathit{f}}_{e}^{\mathrm{t}}={\mathcal{I}}^{\mathrm{v}\to \mathrm{t}}{\mathit{f}}_{e}^{\mathrm{v}}$.
Then,

See Appendix A2 for a proof of Eq. (15).

Remapping between grids requires, as a second step, applying CAAS for property preservation;
thus, we need grid point weight data on each grid.
We define *J*^{t} using Eq. (13).
On the tracer grid,
rather than use Eq. (13),
we interpolate the Jacobian determinant values from the dynamics grid,

for two reasons. First, this operation conserves the area of the element,

by Eq. (15).
We describe the second reason when discussing *ρ* in Sect. 3.1.
On the physics grid,
we follow a similar procedure, explained in Hannah et al. (2021, Sect. 2.2.1),
that is specialized to the finite-volume physics grid.
Here again, the area of an element on the physics grid is the same as on the dynamics grid.
Thus, all three subgrid definitions of the element area agree on element area values.

## 2.4 Direct stiffness summation

In a spectral element method, most operations are performed independently in each element,
often leading to discontinuities in a field across element boundaries.
The global direct stiffness summation (DSS) operator (Fischer and Patera, 1989; Dennis et al., 2012) restores continuity.
At each grid point on an edge of an element,
the multi-valued solution is restored to a single value
by weighted summation of contributions from each element sharing the grid point.
An element's weight at grid point *i*
is the value *w*_{i}*J*_{i} normalized by the sum of all contributing elements' values.
We will need a generalization of this operator
in which an extra factor *σ*_{i} is absorbed into the weight in the weighted sum.
Let ${\mathit{\delta}}_{ij}^{\mathrm{grid}}$ be
1 if global indices *i* and *j* are associated with the same grid point
and 0 otherwise.
For global index *i*, the generalized DSS applied to a preliminary field $\overline{\mathit{y}}$
is written:

where ∑_{e} is the sum over all elements.
In our use of the generalized DSS,
** σ**=

**and $\overline{\mathit{y}}=\overline{\mathit{q}}$. In the standard DSS,**

*ρ***is all one and $\overline{\mathit{y}}=\mathit{\rho}\odot \overline{\mathit{q}}$. These two use cases give the same value for**

*σ**q*

_{i}if

**is continuous across elements. We use the generalized DSS when we want to make**

*ρ***continuous but leave**

*q***discontinuous and, in particular, unmodified from its original value.**

*ρ*Now that we have described the problem setting and the core algorithms we use,
we describe each algorithm that runs in one time step of the Islet method,
advancing the simulation from time step *n* to *n*+1.

As part of an earth system model,
the advection equation Eq. (3) has a source term.
*f* is a function of possibly all the variables in a simulation.
The physics parameterizations compute *f* on the physics grid.
In practice, there are multiple equations of the form Eq. (3) to solve,
e.g., 40 in the E3SM version 2 water cycle simulations (Golaz et al., 2022);
because they decouple during one transport time step,
we can focus on just one in this section.

In summary,
first, air density ** ρ** and trajectory data are remapped from the dynamics grid to the tracer grid.
Second, Δ

**≡**

*q***Δ**

*f**t*is remapped from the physics grid to the tracer grid, where Δ

*t*is the physics parameterization time step. Either of

**or Δ**

*f***is sometimes called a**

*q**tendency*. Third, the remapped tracer tendencies are added to the tracer-grid tracers. Fourth, the tracers are advected on the tracer grid. Fifth, tracer states are remapped to the physics and, optionally, dynamics grids.

## 3.1 Algorithms

There are a number of algorithmic steps in a time step.
To organize these algorithms,
we name each algorithm using the format **Step::algorithm-name** and sometimes **Step::algorithm-name::sub-algorithm-name**.

**Step::density-d2t**.
Interpolate $\mathit{J}\odot {\mathit{\rho}}^{n+\mathrm{1}}$ from the dynamics grid to the tracer grid.
In each element *e*, compute

This grid transfer conserves ${\mathit{J}}_{e}\odot {\mathit{\rho}}_{e}^{n+\mathrm{1}}$ by Eq. (15).
${\mathit{J}}^{\mathrm{t}}\odot ({\mathit{\rho}}^{\mathrm{t}}{)}^{n+\mathrm{1}}$ is not continuous across element boundaries
because *J*^{t} is not,
but continuity is not needed.
Equation (16) implies that
if *ρ*_{e} is constant on the dynamics grid,
then it is constant on the tracer grid too.
Mapping a constant air density exactly is not necessary,
but since it is possible, we do;
this property is the second reason to define *J*^{t} according to Eq. (16).

**Step::tendency-f2t**.
Map the tracer tendency Δ*q*^{n} from the physics grid to the tracer grid.
This step involves multiple sub-algorithms.

**Step::tendency-f2t::bounds**.
In each element *e*, compute and store the minimum and maximum values of $({\mathit{q}}_{e}^{\mathrm{f}}{)}^{n}$,
subsequently, *extrema*.
Let an element's neighborhood contain itself
and every other element that shares a vertex with it.
Augment an element's extrema with the extrema over its neighborhood.
Finally, augment element *e*'s extrema again with the extremal values of the tracer-grid mixing ratio state $({\mathit{q}}_{e}^{\mathrm{t}}{)}^{n}$.
This final extrema update assures that if $\mathrm{\Delta}({\mathit{q}}_{e}^{\mathrm{f}}{)}^{n}=\mathrm{0}$,
then $({\overline{\mathit{q}}}_{e}^{\mathrm{t}}{)}^{n}$ is unmodified in **Step::tendency-f2t::CAAS**.
These final extrema
are the bounds used in subsequent property preservation corrections.

**Step::tendency-f2t::linear-remap**.
In each element *e*,
apply the linear, element-local, conservative, panel-reconstruction remap operator
described in Hannah et al. (2021, Sect. 2.2.3)
to map the tendency from the physics grid to the tracer grid.
In the Islet method, but unlike in Hannah et al. (2021),
the basis used in the mass matrix of the *L*^{2} projection is the Islet basis.
From the previous step's application of **Step::density-d2t**, we have ${\mathit{\rho}}_{e}^{\mathrm{t}}$.
The operator uses this quantity and ${\mathit{\rho}}_{e}^{\mathrm{f}}\odot \mathrm{\Delta}({\mathit{q}}_{e}^{\mathrm{f}}{)}^{n}$ to compute ${\mathit{\rho}}_{e}^{\mathrm{t}}\odot \mathrm{\Delta}({\mathit{q}}_{e}^{\mathrm{t}}{)}^{n}$.

**Step::tendency-f2t::CAAS**.
In each element *e*, apply CAAS,
described in Sect. 2.2.1,
to $({\overline{\mathit{q}}}_{e}^{\mathrm{t}}{)}^{n}\equiv ({\mathit{q}}_{e}^{\mathrm{t}}{)}^{n}+\mathrm{\Delta}({\mathit{q}}_{e}^{\mathrm{t}}{)}^{n}$,
using the weight vector ${\mathit{\theta}}_{e}^{\mathrm{t}}\equiv {\mathit{w}}_{e}^{\mathrm{t}}\odot {\mathit{J}}_{e}^{\mathrm{t}}$
and air density $({\mathit{\rho}}_{e}^{\mathrm{t}}{)}^{n}$.
Target mass *b* is set to the current total element mass
because **Step::tendency-f2t::linear-remap** is conservative;
thus, no adjustment to the total element mass is needed.
The same upper and lower bounds apply to each GLL node in the element.
In this step, the bounds that **Step::tendency-f2t::bounds** computed
are relaxed in each direction by 1 % of the difference between upper and lower bounds,
as discussed in Sect. 2.2.4.
The exact bounds will be enforced in **Step::CEDR::global**.
The constraint set of mass conservation and bounds nonviolation is non-empty
because the constant mixing ratio is a solution,
as explained in Hannah et al. (2021, Sect. 2.3).

**Step::tendency-f2t::DSS**.
At this point, $({\overline{\mathit{q}}}^{\mathrm{t}}{)}^{n}$ is discontinuous across element boundaries.
Neither order of accuracy nor property preservation requires continuity,
but we find that continuity improves the solution quality.
(*ρ*^{t})^{n} can remain discontinuous.
Thus, we apply the generalized DSS, Eq. (17), to $({\overline{\mathit{q}}}^{\mathrm{t}}{)}^{n}$,
using (*ρ*^{t})^{n}
to obtain the continuous field $({\widehat{\mathit{q}}}^{\mathrm{t}}{)}^{n}$.

**Step::advect-interp**. Compute and apply the linear advection ISL operator
described in Sect. 2.1.4.
This step involves two substeps:
computing the grid point trajectories and
computing the interpolants.

**Step::advect-interp::trajectory**.
Compute Eq. (8).
The dynamics component supplies velocity data at the dynamics GLL grid points.
Any of a number of algorithms can compute departure points at time *n* backward in time from dynamics-grid arrival GLL points at time *n*+1.
The Islet method takes as input these departure points as 3D Cartesian departure points.
Next, compute Eq. (9).
In each element, the natural GLL interpolant ℐ^{v→t}
is applied to each of the Cartesian components separately
to provide departure points on the tracer grid.
This procedure implies that
adjacent elements compute identical departure points at shared boundaries.
Finally, for simplicity in the subsequent sphere-to-reference map computations,
the 3D Cartesian points are normalized to the sphere.
For ${n}_{\mathrm{p}}^{\mathrm{v}}=\mathrm{4}$, the ${n}_{\mathrm{p}}^{\mathrm{t}}$-basis departure points are obtained at order of accuracy 4.

**Step::advect-interp::mixing-ratio**.
Compute Eq. (10).
Each departure tracer-grid GLL node is mapped to the containing element, subsequently the *source* element.
Details of finding the source element depend on host-model implementation details and are omitted here;
possibilities include
octree search,
*O*(1) arithmetic for quasiuniform cubed-sphere element grids having certain reference-to-sphere maps,
and search within a predefined element neighborhood
whose size is proportional to maximum wind speed times advection time step.
Then, the corresponding reference coordinates within the source element are computed using Newton's method.
Next, compute Eq. (11).
The mixing ratio value is computed at the departure point using the Islet basis interpolant and the source element's $({\widehat{\mathit{q}}}^{\mathrm{t}}{)}^{n}$ values.
In addition, the source element's stored extrema are associated to this point as bounds.
Finally, the mixing ratio value and bounds are assigned to the target GLL node on the arrival tracer grid.
Departure points and interpolant weights are calculated once and then are reused for each tracer.

**Step::CEDR**.
Apply Communication-Efficient Density Reconstructors (CEDR)
on the tracer grid,
first to each element and then globally.
In this work, we use CAAS, described in Sect. 2.2.1.
Node weights are *w*^{t}⊙*J*^{t}.

**Step::CEDR::local**.
First, apply the element-local CAAS algorithm.
This step is not required
but it reduces the amount of global mass redistribution in **Step::CEDR::global**,
is computationally inexpensive,
involves no interprocess communication,
and, thus, is worth including.
As in **Step::tendency-f2t::CAAS**,
we relax bounds in each direction by 1 % of the difference between upper and lower bounds.
Unlike in **Step::tendency-f2t::CAAS**, in this application of the element-local CAAS algorithm,
any two target GLL nodes in an element may have different bounds;
the bounds depend on the source element for the target GLL node,
as detailed in **Step::advect-interp::mixing-ratio**.
At this point, the global mass has not yet been corrected;
thus, this local CAAS application's mass constraint is to maintain the element's current tracer mass.
The constraint set is not assuredly non-empty.
Thus, ReconstructSafely,
described in Sect. 2.2.2,
wraps the call to CAAS.

**Step::CEDR::global**.
Apply CAAS–point, described in Sect. 2.2.3, on the global tracer grid.
The exact rather than relaxed bounds are applied to each node.
The inputs to CAAS–point are as follows:
the air density (*ρ*^{t})^{n+1};
the weights $\mathit{\theta}\equiv {\mathit{w}}^{\mathrm{t}}\odot {\mathit{J}}^{\mathrm{t}}$;
$b={\sum}_{i}{w}_{i}^{\mathrm{t}}{J}_{i}^{\mathrm{t}}\left({\mathit{\rho}}_{i}^{\mathrm{t}}{)}^{n}\right({\widehat{\mathit{q}}}^{\mathrm{t}}{)}^{n}$,
the global tracer mass after the tendency update;
the mixing ratio bounds computed in **Step::tendency-f2t::bounds**;
and the current mixing ratio values from **Step::CEDR::local**.
Let the output values be (*q*^{t})^{n+1}.

Continuity across element boundaries was restored to the mixing ratio field in **Step::tendency-f2t::DSS**.
Subsequent steps maintain it in exact arithmetic.
However, in finite precision, continuity in (*q*^{t})^{n+1} holds only to a little above machine precision and not exactly.
This numerical discontinuity does not grow in time because in each time step,
both **Step::tendency-f2t::DSS** and **Step::advect-interp::mixing-ratio**
restore exact continuity in finite precision,
although only one restoration per time step is necessary to prevent growth of the numerical discontinuity.
No step of the overall algorithm is sensitive
to this small and roughly temporally constant numerical discontinuity.

At this point, the time step is complete; the remaining computations remap the tracer grid data between grids.

**Step::state-t2f**.
Remap the mixing ratio state to the physics grid.
This is a purely element-local operation.
First, the linear operator described in Hannah et al. (2021, Sect. 2.2.1) is applied to the tracer density.
Second, the element-local CAAS algorithm is applied on the physics grid,
with the extremal mixing ratio values in the element on the tracer grid as the bounds.
For the same reason as in **Step::tendency-f2t::CAAS**,
the constraint set is assuredly non-empty.

**Step::state-t2v**.
The dynamics solver needs one or more mixing ratios on the dynamics grid, e.g.,
specific humidity.
In addition, in our numerical results in Sect. 4,
we compute all errors, except as indicated, on the dynamics grid,
so we use this step to obtain those errors.
First, in an element, ℐ^{t→v} interpolates the tracer-grid mixing ratio to the ${n}_{\mathrm{p}}^{\mathrm{v}}$-basis.
Second, the element-local CAAS algorithm is applied on the dynamics grid
to preserve shape and conserve mass;
the details are as in **Step::state-t2f** but with GLL nodes instead of finite-volume subcells.
The result is (*q*^{v})^{n+1}.
Third, the standard DSS is applied to $({\mathit{\rho}}^{\mathrm{v}}{)}^{n+\mathrm{1}}\odot ({\mathit{q}}^{\mathrm{v}}{)}^{n+\mathrm{1}}$,
where (*ρ*^{v})^{n+1} is the continuous air density from the dynamical core
to obtain continuous tracer density and mixing ratio fields.

**Step::state-v2t**.
In verification problems in Sect. 4,
we need to remap a mixing ratio initial condition from the dynamics grid to the tracer grid.
The algorithm is the same as **Step::state-t2v**
except that the DSS follows the procedure in **Step::tendency-f2t::DSS**
because the air density on the tracer grid is and remains discontinuous at element boundaries.

If ${n}_{\mathrm{p}}^{\mathrm{t}}={n}_{\mathrm{p}}^{\mathrm{v}}$, then tracer transport *p*-refinement (TTPR) is not enabled.
In this case, identity maps replace a subset of the algorithms described in this subsection:
**Step::density-d2t**, the interpolation part of **Step::advect-interp::trajectory**, **Step::state-t2v**, and **Step::state-v2t**.
When describing numerical experiments, we indicate when TTPR is not enabled.

This section presents results for a number of verification problems.
Except in Sect. 4.3,
the equation is the sourceless advection equation, Eq. (4).
Two-dimensional, time-dependent flow, ** u**(

**,**

*x**t*), is prescribed on the sphere.

In most figures, we show results for ${n}_{\mathrm{p}}^{\mathrm{t}}=\mathrm{4}$, 6, 8, 9, 12. The ${n}_{\mathrm{p}}^{\mathrm{t}}={n}_{\mathrm{p}}^{\mathrm{v}}=\mathrm{4}$ case provides a reference because it does not use the TTPR algorithms. The ${n}_{\mathrm{p}}^{\mathrm{t}}=\mathrm{6}$ case has the smallest value of ${n}_{\mathrm{p}}^{\mathrm{t}}$ providing order of accuracy (OOA) greater than 2, in this case, 4. The ${n}_{\mathrm{p}}^{\mathrm{t}}=\mathrm{8}$ case provides OOA 5 and has four times as many nodes as the ${n}_{\mathrm{p}}^{\mathrm{t}}=\mathrm{4}$ basis in two dimensions. The ${n}_{\mathrm{p}}^{\mathrm{t}}=\mathrm{9}$ case provides OOA 6. Finally, the ${n}_{\mathrm{p}}^{\mathrm{t}}=\mathrm{12}$ case provides OOA 8 and has four times as many nodes as the ${n}_{\mathrm{p}}^{\mathrm{t}}=\mathrm{6}$ basis.

Tests follow the procedures detailed in Lauritzen et al. (2012).
Results can be compared with those from many models described in Lauritzen et al. (2014).
We refer to these articles frequently and thus abbreviate them as TS12 (“test suite”) and TR14 (“test results”), respectively.
Initial conditions are generated on the dynamics grid.
Similarly, error diagnostics are computed on the dynamics grid in most cases;
we state the exceptions when they occur.
In most cases, we omit results for simulations without property preservation
as tracer transport modules in earth system models are expected to be property preserving.
We have not attempted to make this section self-contained
as describing details of the large number of verification problems would take too much space.
We recommend the reader not familiar with these problems read TS12.
In addition, we refer to specific figures in TR14 and sometimes TS12
so the reader can compare our results with those from previously documented methods.
We briefly summarize the key characteristics of the verification problems.
There are two prescribed flows: a non-divergent one and a divergent one.
Each prescribes a flow that lasts for *T*=12 d
and such that at time *T*, the exact solution is the same as the initial condition.
This 12 d prescribed flow can be run for multiple cycles to lengthen the simulation.
The non-divergent flow creates a filament of maximum aspect ratio at time $T/\mathrm{2}$.
The divergent flow tests treatment of divergence.
There are four initial conditions that share the feature of placing
two circular shapes at two points along the equator:
the *C*^{∞} Gaussian hills,
the *C*^{1} cosine bells,
the correlated cosine bells used in the mixing diagnostic,
and the discontinuous slotted cylinders.
Assessing the behavior of a transport method on tracers having various degrees of smoothness
is important because atmosphere tracers can be smooth or non-smooth.

This article does not study time integration methods to generate the dynamics-grid departure points.
Thus, to remove temporal errors due to time integration algorithms,
we use an adaptive Runge–Kutta method (Dormand and Prince, 1980; Shampine and Reichelt, 1997) with a tight tolerance (10^{−8} with an exception noted later) to integrate trajectories very accurately.
All tests for ${n}_{\mathrm{p}}^{\mathrm{t}}>\mathrm{4}$ use TTPR with ${n}_{\mathrm{p}}^{\mathrm{v}}=\mathrm{4}$ unless we state otherwise.
Because ${n}_{\mathrm{p}}^{\mathrm{v}}=\mathrm{4}$ in these verification problems, the ${n}_{\mathrm{p}}^{\mathrm{t}}=\mathrm{4}$ configuration does not use TTPR.

As we discussed in Sect. 2.2.3, there are two natural approaches when applying CAAS globally, which we refer to as CAAS–CAAS and CAAS–point. CAAS–CAAS tends to give slightly more accuracy than CAAS–point. However, it does not permit the bound relaxations in the element-local CAAS applications that improve the solution quality when simulating source terms sensitive to round-off errors, as discussed in Sect. 2.2.4. However, because it provides slightly more accurate results, we use CAAS–CAAS for ${n}_{\mathrm{p}}^{\mathrm{t}}=\mathrm{4}$, except where noted, while using CAAS–point for all other ${n}_{\mathrm{p}}^{\mathrm{t}}$ values, thus maximizing the accuracy obtained with ${n}_{\mathrm{p}}^{\mathrm{t}}=\mathrm{4}$ to provide the best baseline (${n}_{\mathrm{p}}^{\mathrm{t}}={n}_{\mathrm{p}}^{\mathrm{v}}=\mathrm{4}$) accuracy.

The air density on the dynamics grid, *ρ*^{v},
is not needed for the advection of mixing ratios,
but it and its remapped versions are needed for property preservation.
In these verification problems, we have no independent computation of density *ρ*^{v}
as occurs when transport is coupled to a full dynamical core.
For the non-divergent flow problem, we could set *ρ*^{v} to a constant,
but for the divergent flow problem, a constant is incorrect.
Thus, within our standalone verification code,
we need a means to compute *ρ*^{v} as a surrogate for a dynamics solver.
Appendix B describes the linear advection part of computing *ρ*^{v}.
After the linear advection operator is applied,
because it is not mass-conserving,
the density is corrected by adding $\mathrm{\Delta}m/{a}_{\mathrm{total}}$ to each grid point,
where Δ*m* is the global mass discrepancy after advection,
and *a*_{total} is the total area of the grid.
Negative density does not occur in these verification problems.
The algorithm for *ρ*^{v} has OOA 2 for ${n}_{\mathrm{p}}^{\mathrm{v}}=\mathrm{4}$.

We use a quasiuniform, equiangular cubed-sphere element grid.
Let a cube face of the cubed-sphere element grid have *n*_{e}×*n*_{e} elements.
Each element has an ${n}_{\mathrm{p}}^{\mathrm{v}}\phantom{\rule{-0.125em}{0ex}}\times \phantom{\rule{-0.125em}{0ex}}{n}_{\mathrm{p}}^{\mathrm{v}}$ tensor grid of GLL nodes
and, thus, ${n}_{\mathrm{p}}^{\mathrm{v}}-\mathrm{1}$ intervals between GLL nodes along each direction of an element.
GLL nodes are mapped to the sphere using the isoparametric map for a bilinear element,
as described in Appendix A of Guba et al. (2014) and Eq. (7).
Long tracer time steps correspond to 6*n*_{e} steps per *T*;
short, 30*n*_{e}.
For the test flows, these correspond to, respectively,
approximately 5.5 and 1 times the maximum Courant number.
These are the same time step settings as the two Conservative Semi-Lagrangian Multi-tracer (CSLAM) (Lauritzen et al., 2010) model configurations used in TR14.
Essentially all SL methods, particularly when given exact trajectory data,
exhibit greater error with smaller time step, e.g., CSLAM in TR14.
This is because the only source of error given exact trajectories
is the remap error.
Smaller time steps correspond to more remaps to reach a fixed simulation time.

Figures 3–9 and 15
show convergence plots,
and we explain the format of these figures here.
A curve's marker corresponds to ${n}_{\mathrm{p}}^{\mathrm{t}}$, as listed in the legend.
To maximize font size and minimize notational clutter in figures,
in figures, we use *n*_{p} rather than ${n}_{\mathrm{p}}^{\mathrm{t}}$ , and omit “=”; e.g., “${n}_{\mathrm{p}}^{\mathrm{t}}=\mathrm{8}$” is written “*n*_{p} 8”.
Additionally, the legend omits “*n*_{p}” entirely.
The *x* axis is the average dynamics-grid node spacing at the equator in degrees for a cubed-sphere grid
with north and south cube faces centered at the poles.
Thus, for example, *n*_{e}=5 and ${n}_{\mathrm{p}}^{\mathrm{v}}=\mathrm{4}$ correspond to the resolution

The *y* axis is log _{10} of the relative error, with the norm or norms indicated in each figure.
The title of a plot lists details of the configuration:
the test flow, the initial condition (IC), the time step (short or long),
and, if property preservation is disabled, an extra line stating that.
Many plots have a dotted line showing the slope for OOA 2.

## 4.1 Accuracy for *C*^{∞} and *C*^{1} tracers

### 4.1.1 Accuracy limited only by trajectories

The first experiment tests the OOA limit due to computing trajectories
on the tracer grid using velocities from the dynamics grid.
Because the dynamics grid uses ${n}_{\mathrm{p}}^{\mathrm{v}}=\mathrm{4}$ and ℐ^{v→t} uses the natural GLL basis,
we expect this OOA limit to be 4.
Property preservation is turned off to expose this OOA limit.

Sometimes verification problems have unexpected features that interact with a method to produce higher-accuracy solutions than would occur in a more realistic problem. Of particular concern is the symmetry of the flow in time around the midpoint time, 6 d. To be sure our trajectory interpolation procedure is not interacting with this symmetry to produce artificially higher accuracy, in this experiment, we compute the solution error at the midpoint time as well as the final time.

The test uses the non-divergent flow, the Gaussian hills IC, and the long time step.
The midpoint reference solution is computed using one 6 d step and the natural GLL basis.
The time integrator's relative error tolerance is set to 10^{−14} rather than its usual 10^{−8} for this step.
One remap step with the natural GLL basis provides a midpoint solution much more accurate than the time-stepped case,
thus serving as an appropriate reference.

Figure 3 shows the results.
The dashed curves show the errors measured at the midpoint, half of a cycle of the 12 d problem;
the solid curves, the usual endpoint.
The numbers at the right side of the plot show empirical OOA
computed using the final two points of the one-cycle results.
Numbers for the half-cycle curves are omitted
because each pair of curves has almost exactly parallel lines
for resolution at least as fine as 0.75^{∘}.
In addition, the empirical OOA for ${n}_{\mathrm{p}}^{\mathrm{t}}=\mathrm{12}$ is omitted because the
curves nearly overlap those for ${n}_{\mathrm{p}}^{\mathrm{t}}=\mathrm{9}$.
For fine enough resolution, all curves should have empirical OOA limited to 4,
but for ${n}_{\mathrm{p}}^{\mathrm{t}}=\mathrm{8}$, the full convergence regime is not reached in this plot,
thus giving OOA 5 for ${n}_{\mathrm{p}}^{\mathrm{t}}=\mathrm{8}$.
The ${n}_{\mathrm{p}}^{\mathrm{t}}=\mathrm{4}$ curves have OOA less than 4 due to a spatial OOA limit of 2 in the full convergence regime.
We see OOA at about 4 for ${n}_{\mathrm{p}}^{\mathrm{t}}=\mathrm{9}$ and 12.
The ${n}_{\mathrm{p}}^{\mathrm{t}}=\mathrm{9}$ curve shows at coarse resolution higher OOA,
governed by the spatial error,
and then a drop to 4 as the temporal error becomes dominant.
Importantly, the half- and full-cycle errors are very close for each value of ${n}_{\mathrm{p}}^{\mathrm{t}}$,
demonstrating that the endpoint error metrics are valid measurements
for the Islet method.

### 4.1.2 Accuracy data for other standard configurations

The next figures show accuracy for standard configurations used in TR14. Although we explain how each figure can be compared with corresponding ones in TR14, these figures also stand on their own as simply convergence plots for various test cases.

Figures 4 and 5 can be compared with Figs. 1 and 2 in TR14.
They evaluate error on an infinitely smooth IC.
The Islet method with a high-order basis compares extremely favorably with the methods in TR14.
For example, the most accurate shape-preserving method in TR14
for the non-divergent flow with the Gaussian hills IC is the Hybrid Eulerian Lagrangian–Non-Diffusive (HEL-ND) method run with Courant number 1.0
(HEL-ND-CN1.0) by a substantial margin
(cyan curves in Fig. 1, bottom right, of TR14).
The ${n}_{\mathrm{p}}^{\mathrm{t}}=\mathrm{12}$ Islet scheme with the long time step, Fig. 4,
is approximately three times more accurate than HEL-ND-CN1.0 in the *l*_{2} norm
at resolution 0.375^{∘}, and approximately twice as accurate at resolution 3^{∘}.
Yet, HEL-ND is, quoting TR14, an “unphysical” method.
It is run for comparison with the practically useful HEL scheme.
After HEL-ND, the next most accurate method in the *l*_{2} norm at 0.375^{∘} resolution is CSLAM-CN5.0.
The Islet method with the long time step, the same as that of CSLAM-CN5.0,
is at least as accurate for ${n}_{\mathrm{p}}^{\mathrm{t}}\ge \mathrm{8}$.
With the short time step, the same as that of CSLAM-CN1.0,
the Islet method is at least as accurate as CSLAM-CN1.0 for ${n}_{\mathrm{p}}^{\mathrm{t}}\ge \mathrm{6}$.
At 3^{∘} resolution, no method other than HEL-ND-CN1.0 provides *l*_{2} norm
below 10^{−2}; the Islet method does for ${n}_{\mathrm{p}}^{\mathrm{t}}\ge \mathrm{8}$ with the long time step, and
${n}_{\mathrm{p}}^{\mathrm{t}}\ge \mathrm{10}$ (only ${n}_{\mathrm{p}}^{\mathrm{t}}=\mathrm{12}$ is shown) with the short time step.

Now we must remind the reader that, in this article, we use nearly exact trajectories on the dynamics grid because the description of practical trajectory methods is outside the scope of this article. However, first, many of the schemes in TR14 use nearly exact trajectories, as well, e.g., CSLAM. Second, when using a practical trajectory algorithm, highly accurate, even if not exact, trajectories are possible because a trajectory over a tracer time step can be computed from multiple dynamics time steps. Thus, the use of exact trajectories is only slightly unrealistic. Third, the diagnostic values in Sect. 4.1.3, 4.1.4, and 4.3.2 are roughly independent of temporal errors.

Figures 6 and 7 provide data that can be compared with the top panel of Fig. 3 in TR14.
They evaluate error on a *C*^{1} IC.
The horizontal dash-dotted line provides the relative *l*_{2}-error-norm value of 0.033 by which the
“minimal resolution” diagnostic value is determined;
the coordinate of the intersection between the *l*_{2}-norm curve and this reference line is the value.
A larger value is better.
For example, with a long time step,
for ${n}_{\mathrm{p}}^{\mathrm{t}}=\mathrm{8}$, this value is a little coarser than 3^{∘};
for ${n}_{\mathrm{p}}^{\mathrm{t}}=\mathrm{12}$, approximately 6^{∘}.
For comparison, no model in TR14 reports a value larger than 2.5^{∘}.

Figures 8 and 9 provide data that can be compared with the top two panels in Fig. 16 of TR14 given the minimum resolutions plotted in Fig. 3 of TR14. They are like Figs. 6 and 7 but, here, the divergent flow is used.

### 4.1.3 Filament diagnostic

Figure 10 shows results for the filament diagnostic
described in Sect. 3.3 of TS12
to compare with Fig. 5 in TR14.
The two dynamics grid resolutions are as prescribed in TR14.
We used the code distributed with TS12 to compute the results.
The diagnostic uses the non-divergent flow and the cosine bells IC.
In this test, the midpoint solution is analyzed
to determine the quality of the filamentary structure.
See Fig. 13 for an illustration
of the filamentary structure at the simulation midpoint, although with the slotted cylinders IC.
For each possible value *τ* of the tracer mixing ratio at the initial time,
the area over which the mixing ratio is at least *τ* at the midpoint time is computed.
For the cosine bells IC, $\mathit{\tau}\in [\mathrm{0.1},\mathrm{1}]$.
The diagnostic is then this area divided by the correct area,
which for any non-divergent flow is the area at the initial time.
The perfect diagnostic value is 100 *%* for all $\mathit{\tau}\in [\mathrm{0.1},\mathrm{1}]$ and 0 otherwise.
In each plot in Fig. 10,
the *x* axis is *τ*
and the *y* axis is the diagnostic value.
The diagnostic is computed for the dynamics-grid resolutions and time step lengths listed in the legend.
Note that the *y* axis limits are tighter with increasing ${n}_{\mathrm{p}}^{\mathrm{t}}$.

A subtlety with this diagnostic is that the area calculation must use the quadrature method of the discretization. On the dynamics grid, the resulting curve is noisier than is implied by the underlying solution on the tracer grid. Thus, in Fig. 10, we show the diagnostic as computed on the dynamics grid in the top row of plots; on the tracer grid, in the bottom row.

There is no summary number that can be compared directly with the results in Fig. 5 of TR14,
but, visually, the curves for ${n}_{\mathrm{p}}^{\mathrm{t}}\ge \mathrm{8}$, the resolution 1.5^{∘}, and on the tracer grid
are among the best of those in TR14 at the resolution 1.5^{∘}.

### 4.1.4 Mixing diagnostic

Figures 11 and 12 show results for the mixing diagnostics described in Sect. 3.5 of TS12 to compare with Figs. 11–14 in TR14. Again, the two dynamics grid resolutions are as prescribed in TR14.

The test uses the non-divergent flow with two ICs.
The diagnostics assess preservation of the nonlinear correlation of two tracers;
one mixing ratio is a nonlinear function of the other.
The mixing ratio of each tracer is on an axis,
cosine bells on the *x* axis,
and the correlated cosine bells field on the *y* axis.
Each dot is a grid-point sample from the dynamics grid of the two mixing ratios.
Like the filament diagnostic, the analysis is done at the solution midpoint rather than the endpoint.
Also as for the filament diagnostic, we used the code distributed with TS12 to compute the results.
${n}_{\mathrm{p}}^{\mathrm{t}}$ and the time step length are printed in each plot.
The diagnostic values *l*_{r} and *l*_{u}, explained in a moment, are also printed in each plot;
the parenthesized “(v)” suffix means the value is measured on the dynamics grid,
while the un-suffixed values are measured on the tracer grid.

At the initial time, all points are on the curve.
Perfect nonlinear correlation corresponds to staying on the upper curve
as the simulation proceeds.
Points that drop into the convex hull of the starting points
can be interpreted as physical mixing of air parcels,
since such mixing results in a linear combination of two points on the curve.
The diagnostic value *l*_{r} measures this type of error; smaller is better.
In Figs. 11–14 in TR14,
the smallest value of *l*_{r} at 1.5^{∘} among the property-preserving methods
is $\mathrm{2.15}\times {\mathrm{10}}^{-\mathrm{4}}$ by the University of California, Irvine, Second-Order Moments (UCISOM) method run with Courant number 5.5 (UCISOM-CN5.5),
except for a value of 0 by HEL-ND,
which, again, cannot be used in practice.
For the long time step, the Islet method gives at least as small a value for ${n}_{\mathrm{p}}^{\mathrm{t}}\ge \mathrm{6}$;
for the short time step, ${n}_{\mathrm{p}}^{\mathrm{t}}\ge \mathrm{8}$.

Values outside the triangle are overshoots,
possible only if the method is not strictly shape preserving.
Since our method is, this diagnostic value, *l*_{o}, is always 0,
and is not displayed in the figures.

Values outside the convex hull cannot be described as physical mixing of parcels
and, thus, are purely numerical artifacts of a method.
The corresponding diagnostic is *l*_{u} and, again, smaller is better.
This diagnostic is more difficult to compare than *l*_{r}
because very dissipative methods tend to have a large value of *l*_{r}
and, consequently, a very small value for *l*_{u}.
In contrast, a very accurate method, for which *l*_{r} is small,
can have a larger *l*_{u} value than a very dissipative method.
One means of comparison is to consider the best *l*_{u} values
among methods that obtain, say, ${l}_{r}\le \mathrm{5}\times {\mathrm{10}}^{-\mathrm{4}}$.
In Figs. 11–14 in TR14,
the smallest value of *l*_{u} at 1.5^{∘} under this restriction is 0
obtained by the HEL-CN1.0 and HEL-CN5.5 methods.
These HEL variants are practically usable
unlike HEL-ND, and are designed to preserve tracer correlations exactly.
Other than the HEL methods, the next best value is $\mathrm{4.80}\times {\mathrm{10}}^{-\mathrm{5}}$,
again by the UCISOM-CN5.5 method.
For both the long and short time steps, at 1.5^{∘},
the Islet method gives at least as small a value for ${n}_{\mathrm{p}}^{\mathrm{t}}\ge \mathrm{8}$.
However, even with the constraint on *l*_{r}, comparison is not straightforward
as the UCISOM methods are not strictly shape preserving
and so have *l*_{o}>0.

## 4.2 Slotted cylinders

We observe that in both the filament and mixing diagnostics of TS12,
${n}_{\mathrm{p}}^{\mathrm{t}}\ge \mathrm{6}$ gives excellent results;
${n}_{\mathrm{p}}^{\mathrm{t}}=\mathrm{12}$, nearly perfect.
Figures 13 and 14 show solution quality using latitude–longitude images
and further underscore these observations.
The problem is the non-divergent flow with the slotted cylinders IC,
at resolutions 1.5 and 0.75^{∘}, with the long and short time steps.
The text in the individual images in Fig. 14
provides normwise accuracy at the end of one cycle
and deviation from the initial extrema.
${\mathit{\varphi}}_{min}\ge \mathrm{0}$ and ${\mathit{\varphi}}_{max}\le \mathrm{0}$ are consistent with no global extrema violation.
Compare Fig. 13 with Figs. 7–10 in TR14
and both figures with Fig. 7 in TS12.
Although TR14 does not provide error norm values for this problem,
those in Fig. 7 of TS12 can be compared with the
Islet method's values at 1.5^{∘} resolution and the long time step (left column of Fig. 14).
The Islet method's values of *l*_{2}, *l*_{inf}, *ϕ*_{min}, and *ϕ*_{max}
are at least as good as those in Fig. 7 of TS12 for ${n}_{\mathrm{p}}^{\mathrm{t}}\ge \mathrm{6}$.

## 4.3 Source term

Now we move to verification problems that include a source term:
*f*≠0 in Eq. (3).

### 4.3.1 Accuracy

The first test verifies the property-preserving remaps between the physics and tracer grids.
The test is constructed as follows.
Two tracer mixing ratios, a source *q*_{1}=*s* and a manufactured tracer *q*_{2}=*m*, are paired.
At time *t*=0, *m*(** x**,

*t*) is set to 0, where

**is position on the sphere. A tendency Δ**

*x**m*is applied to

*m*on the physics grid: $\mathrm{\Delta}m(\mathit{x},t)\equiv -\left[\mathrm{cos}\right(\mathrm{2}\mathit{\pi}(t+\mathrm{\Delta}t)/T)-\mathrm{cos}(\mathrm{2}\mathit{\pi}t/T\left)\right]s(\mathit{x},t)/\mathrm{2}$, so that the exact solution is $m(\mathit{x},t)=(\mathrm{1}-\mathrm{cos}(\mathrm{2}\mathit{\pi}t/T\left)\right)s(\mathit{x},t)/\mathrm{2}$ and, in particular, $m(\mathit{x},T/\mathrm{2})=s(\mathit{x},T/\mathrm{2})$. To compute the tendency, the state

*s*(

**,**

*x**t*) must be remapped to the physics grid. Thus, the results for this test depend on both grid-transfer directions. We measure the error at time $T/\mathrm{2}$, on the dynamics grid as usual, as in Fig. 3. Figure 15 shows the results. We run this test with ${n}_{\mathrm{f}}={n}_{\mathrm{p}}^{\mathrm{t}}$ (dash-dotted lines) and

*n*

_{f}=2 (dashed lines). The solid lines show the error in $s(\mathit{x},T/\mathrm{2})$ as a reference. The dotted line provides the OOA-2 reference. We see that when ${n}_{\mathrm{f}}={n}_{\mathrm{p}}^{\mathrm{t}}$, the errors are nearly the same as those for $s(\mathit{x},T/\mathrm{2})$; for ${n}_{\mathrm{p}}^{\mathrm{t}}={n}_{\mathrm{f}}=\mathrm{4}$, the curves overlap at the resolution of the plot. As one expects, when

*n*

_{f}=2, the error in

*m*(

**,**

*x**t*) is much larger than in $s(\mathit{x},T/\mathrm{2})$, but the OOA remains 2.

### 4.3.2 Toy chemistry diagnostic

The toy chemistry verification problem is described in Lauritzen et al. (2015), subsequently TC15.
The problem consists of two tracer mixing ratios, *q*_{i}=*X*_{i}, $i=\mathrm{1},\mathrm{2}$, that interact according to chemical kinetic equations that are nonlinear in one of them:
$\mathrm{D}{X}_{i}(\mathit{y},t)/\mathrm{D}t={f}_{i}(\mathit{y},{X}_{\mathrm{1}},{X}_{\mathrm{2}})$,
where ** y** is the spatial coordinate on the sphere.
(We use

**here to avoid overlap with the species symbol,**

*y**X*, that is used in TC15.) The tracers are composed of a monatomic and a diatomic molecule, respectively, of the same atomic species. The reactions are extremely sensitive to solar insolation. The sun's position is held fixed with respect to the grid. As a result, the largest-scale spatial pattern one sees in the fields is the boundary dividing nonzero (day) and zero (night) solar insolation, the solar

*terminator*; this boundary is particularly visible in the right image of Fig. 17, a figure that we will describe subsequently. The ICs are designed so that the sum over the atomic mixing ratio at each point in space is a constant, ${\overline{X}}_{T}$. The source terms have this property too, since they model chemical reactions. Thus, in the exact solution, ${\overline{X}}_{T}$ is maintained at every point in space and time. Let

*X*

_{T}without a bar be the corresponding measured quantity. The toy chemistry diagnostics are then

*c*

_{2}(

*t*), the

*l*

_{2}norm of ${X}_{T}-{\overline{X}}_{T}$ at time

*t*normalized by ${\overline{X}}_{T}$, and

*c*

_{∞}(

*t*) the same but for the

*l*

_{∞}norm.

As explained in the context of Eq. (14) in TC15,
any advection operator that is *semi-linear* will produce a perfect diagnostic value of 0 when using exact arithmetic.
Linear operators are semi-linear; the CEDR algorithms we use are as well, as explained in Bradley et al. (2019);
and a composition of semi-linear operators is too.
Thus, the Islet method is semi-linear,
and deviation from 0 in the diagnostic values is wholly due to the effects of finite precision arithmetic.

We compute the diagnostics, as usual, on the dynamics grid. Following the Islet tracer transport method described in Sect. 3.1, the source terms are computed on the physics grid using states remapped from the tracer grid, and then the computed tendencies are remapped to the tracer grid.

It is already known that the Eulerian spectral element tracer transport method yields poor values for this diagnostic due to finite-precision effects of the limiter (Lauritzen et al., 2017). The Islet method with property preservation using CAAS–CAAS does as well. Again, in exact arithmetic, each of these methods would produce perfect values. The poor diagnostic values are due to quickly accumulating machine-precision round-off errors that break semi-linearity in finite precision. The interaction of the chemistry source term with exact bounds in element-local limiter applications is responsible for this fast accumulation. In contrast, the Islet method using CAAS–point and relaxed-bound, element-local CAAS applications produces good diagnostic values. The relaxed bounds in the element-local part make unnecessary many of the mixing ratio adjustments that lead to loss of semi-linearity in finite precision. Recall that CAAS–point, applied at the end of the step, imposes the exact bounds, so at the end of a time step, shape preservation still holds to machine precision. Clips to bounds must still occur, but adjustments to other grid points to compensate are smaller because the adjustments are spread over many more grid points.

Figure 16 shows the diagnostic values
for the case of non-divergent flow, 1^{∘} dynamics-grid resolution, and a 30 min time step,
where these configuration details are prescribed in TC15.
For the ${n}_{\mathrm{p}}^{\mathrm{t}}=\mathrm{4}$ case, we use CAAS–point rather than CAAS–CAAS as previously,
since we already know that CAAS–CAAS will produce poor values.
The diagnostic is usually plotted over the course of one cycle (12 d) of the flow,
but it is useful to view it over multiple cycles.
Figure 16 shows 10 cycles on the *x* axis for a total of 120 d.
The *y* axis is the diagnostic value.
Solid lines plot *c*_{2}(*t*); dashed, *c*_{∞}(*t*).
Markers are placed on the curves at the start of each cycle to differentiate the curves.
In each case, ${n}_{\mathrm{f}}={n}_{\mathrm{p}}^{\mathrm{t}}$.
We choose this value of *n*_{f} because
the toy chemistry source term has an extremely large gradient at the terminator and,
thus, it makes sense to compute the physics tendencies at high spatial resolution.
The ${n}_{\mathrm{p}}^{\mathrm{t}}=\mathrm{4}$ case with CAAS–point is greatly improved relative to the Eulerian spectral element results shown in Fig. 7 of TC15,
even after 10 cycles instead of the one cycle shown in that figure.
In Fig. 7 of TC15, ${c}_{\mathrm{2}}\approx {\mathrm{10}}^{-\mathrm{2}}$ and ${c}_{\mathrm{\infty}}>{\mathrm{10}}^{-\mathrm{1}}$ at the end of one cycle,
compared with approximately 10^{−7} and 10^{−5}, respectively, for ${n}_{\mathrm{p}}^{\mathrm{t}}=\mathrm{4}$ at the end of 10 cycles.
For ${n}_{\mathrm{p}}^{\mathrm{t}}>\mathrm{4}$, the growth in error is very small, with ${c}_{\mathrm{2}}<{\mathrm{10}}^{-\mathrm{10}}$ through 10 cycles.

To illustrate what these diagnostics measure,
Fig. 17 shows latitude–longitude images of the monatomic tracer at the end of the first cycle
for ${n}_{\mathrm{p}}^{\mathrm{t}}={n}_{\mathrm{f}}=\mathrm{4}$ with CAAS–CAAS (left)
and ${n}_{\mathrm{p}}^{\mathrm{t}}={n}_{\mathrm{f}}=\mathrm{8}$ with CAAS–point (right).
Note that the images in TS12 and TR14 are plotted with longitude ranging from 0 to 2*π*;
those in TC15, from −*π* to *π*.
In our latitude–longitude figures so far,
we have chosen the convention used in TS12 and TR14, and we continue to use it in these toy-chemistry images.
Thus, these images are circularly shifted horizontally by half the image width relative to those in TC15.
The globally extremal tracer values are printed in the upper right quadrant of each image.
The correct maximum is $\mathrm{4}\times {\mathrm{10}}^{-\mathrm{6}}$ and the correct minimum is at least 0.
The right image is free of noise and satisfies these bounds.
The left image shows substantial noise, as we expect when using exact bounds in the local property preservation problems,
and consistent with previous observations about spectral element transport.
Other than noise and some filaments that grow from the noise, the two images are qualitatively similar.
Figure 18 shows images in the same format,
but the quantity is now $({X}_{T}-{\overline{X}}_{T})/{\overline{X}}_{T}$ at the end of the first cycle.
The correct value is 0 everywhere.
In the right image, the pointwise relative error is a little better than 10^{−11},
consistent with the *l*_{∞}-norm diagnostic value for ${n}_{\mathrm{p}}^{\mathrm{t}}=\mathrm{8}$ at the end of the first cycle in Fig. 16.

We described the Islet method, a property-preserving tracer transport method to support a three-grid atmosphere model, one with a shared element grid but separate subelement grids for physics parameterizations, dynamics, and tracer transport. This configuration permits the modeler to create a dynamics grid with a tolerable CFL-limited time step independent of the other two subcomponents, while physics parameterizations and tracer transport can run at resolutions potentially substantially higher. The shared element grid minimizes communication during remaps between grids, and almost all operations are local to each element, making the Islet method extremely efficient.

Section 4 presented a number of verification problems with diagnostics assessing accuracy, order of accuracy, numerical mixing, filament preservation, and nonlinear correlation preservation. The corresponding figures and diagnostic values can be compared directly with those of several other methods. The Islet method performs well in a number of detailed comparisons.

Possible future work includes the following.
First, because all operations except the DSS are either local to the element or act on element-level scalars,
both the tracer and the physics grids permit various types and subsets of
spatially and temporally adaptive and tracer-dependent refinement and derefinement,
possibly in combination with already existing E3SM regionally refined models (RRM).
Future work could develop procedures to permit spatially and temporally varying ${n}_{\mathrm{p}}^{\mathrm{t}}$ and *n*_{f} values,
and incorporate these into the Islet method.

Second, this article focuses on the horizontal dimensions or, more generally, 2D Lagrangian levels in a three-dimensional discretization. However, it is necessary to increase simultaneously both the horizontal and vertical resolutions of a plume to capture and maintain its structure (Eastham and Jacob, 2017). Future work using methods similar to those in Sect. 3 should address the vertical dimension.

Third, it is possible to recover local mass conservation in an ISL method by modifying the linear advection operator's coefficients so that the sum of the mass over all the target points associated with a source element is consistent with the source element's total mass (Kaas, 2008). This coefficient modification step might be compatible with the Islet method and should be explored.

Finally, we predict that applications related to, in particular, aerosols will benefit from the Islet method. In future work, we intend to integrate the Islet method into the E3SM Atmosphere Model to investigate its impact on science applications.

The GLL basis functions yield an unstable ISL method for *n*_{p}≥4.
We developed a set of bases on GLL nodes
that yield ISL methods that meet a necessary condition for stability
on the test problem of uniform periodic flow.
We refer to these as the *Islet bases*
and use these bases in this article.
Bradley et al. (2021, Sects. 2 and 3) and a forthcoming article based on that material
detail the derivation of the Islet bases,
while this appendix summarizes the results of that derivation for completeness.
This article can be understood equally well by assuming the standard GLL bases are used;
only the numerical results depend on the details of the basis functions.

A 1D reference element has domain $[-\mathrm{1},\mathrm{1}]$.
Higher-dimensional basis functions are tensor products of 1D basis functions.
The bases are nodal, meaning a basis function has the value 1 at one node and 0 at all others.
The nodes are GLL.
Let ${x}_{G}^{{n}_{\mathrm{p}}}\left(i\right)$, $i\in \mathit{\{}\mathrm{0},\mathrm{\dots},{n}_{\mathrm{p}}-\mathrm{1}\mathit{\}}$, be the reference coordinates of the GLL nodes.
Let *region* $r\in \mathit{\{}\mathrm{0},\mathrm{\dots},{n}_{\mathrm{p}}-\mathrm{2}\mathit{\}}$ be the segment $\left[{x}_{G}^{{n}_{\mathrm{p}}}\right(r),{x}_{G}^{{n}_{\mathrm{p}}}(r+\mathrm{1}\left)\right]$.
A 1D basis function is a continuous piecewise polynomial over $[-\mathrm{1},\mathrm{1}]$;
in each region, it is a polynomial.
Each region *r* has an associated ordered list of *support nodes*, denoted ${\mathcal{I}}_{r}^{{n}_{\mathrm{p}}}$.
Each list ${\mathcal{I}}_{r}^{{n}_{\mathrm{p}}}$ is a subset of $\mathit{\{}\mathrm{0},\mathrm{\dots},{n}_{\mathrm{p}}-\mathrm{1}\mathit{\}}$.
Let $\left|{\mathcal{I}}_{r}^{{n}_{\mathrm{p}}}\right|$ be the number of elements in ${\mathcal{I}}_{r}^{{n}_{\mathrm{p}}}$,
let ${n}_{\mathrm{p}}^{\mathrm{submin}}\equiv {min}_{r}\left|{\mathcal{I}}_{r}^{{n}_{\mathrm{p}}}\right|$,
and let ${\mathcal{I}}_{r}^{{n}_{\mathrm{p}}}\left(j\right)$ be the *j*th element.

In region *r* and thus for $x\in \left[{x}_{G}^{{n}_{\mathrm{p}}}\right(r),{x}_{G}^{{n}_{\mathrm{p}}}(r+\mathrm{1}\left)\right]$,
basis function *k* is 0 if $k\notin {\mathcal{I}}_{r}^{{n}_{\mathrm{p}}}$
and otherwise has the value given by the Lagrange polynomial

We construct a basis for each ${n}_{\mathrm{p}}\in \mathit{\{}\mathrm{4},\mathrm{\dots},\mathrm{13}\mathit{\}}$.
Each basis is described by a set of support node lists, one ${\mathcal{I}}_{r}^{{n}_{\mathrm{p}}}$ per region.
For *n*_{p}=4, there are additional details that we describe in Appendix A1.
There are two types of bases,
each corresponding to a method of describing the support node lists.

The first type of basis is the *offset nodal basis*.
Each list ${\mathcal{I}}_{r}^{{n}_{\mathrm{p}}}$ is given by an *offset*, which is the first index in the list
and its size, $\left|{\mathcal{I}}_{r}^{{n}_{\mathrm{p}}}\right|$.
For example, *n*_{p}=7, offset 2, and size 4 correspond to
support nodes ${\mathcal{I}}_{r}^{\mathrm{7}}=\mathit{\{}\mathrm{2},\mathrm{3},\mathrm{4},\mathrm{5}\mathit{\}}$.
In addition, the basis is symmetric, meaning basis function ${\mathit{\varphi}}_{i}^{{n}_{\mathrm{p}}}\left(x\right)={\mathit{\varphi}}_{{n}_{\mathrm{p}}-\mathrm{1}-i}^{{n}_{\mathrm{p}}}(-x)$.
Thus, first, support nodes are specified for regions 0 through $\lfloor {n}_{\mathrm{p}}/\mathrm{2}\rfloor -\mathrm{1}$,
and the support nodes for the remaining regions are determined by symmetry.
Second, if *n*_{p} is even, then the middle region, $r={n}_{\mathrm{p}}/\mathrm{2}-\mathrm{1}$,
has support nodes ${\mathcal{I}}_{r}^{{n}_{\mathrm{p}}}$ that are symmetric
around reference coordinate 0.

The second type of basis is the *nodal basis*.
An offset nodal basis is a nodal basis
whose support nodes all permit the more compact description of offset and size.
A general nodal basis has at least one support nodal list that cannot be described
by just an offset and a size.
Instead, the list is stated explicitly.
Again, all nodal bases are symmetric.

The order of accuracy of the ISL method using an Islet GLL basis is ${n}_{\mathrm{p}}^{\mathrm{submin}}-\mathrm{1}$ if there is no property preservation step and the initial condition is ${C}^{{n}_{\mathrm{p}}^{\mathrm{submin}}-\mathrm{2}}$.

If ${n}_{\mathrm{p}}^{\mathrm{submin}}={n}_{\mathrm{p}}$, then the basis is the natural GLL basis. Thus, any statement about an Islet GLL basis also holds for the natural GLL basis.

A search procedure searches for nodal bases that provide ISL methods
that are stable on the test problem of uniform flow on a uniform grid.
An accuracy heuristic is used to find the most accurate basis satisfying this stability condition.
Table A1 lists nodal bases, most of them offset nodal bases,
for ${n}_{\mathrm{p}}\in \mathit{\{}\mathrm{5},\mathrm{\dots},\mathrm{13}\mathit{\}}$.
In addition, a special basis is constructed for *n*_{p}=4.

## A1 *n*_{p}=4 Islet basis

For *n*_{p}=4, we boost the accuracy of the basis by combining an offset nodal basis and the GLL basis.
The middle region has support nodes ${\mathcal{I}}_{\mathrm{1}}=\mathit{\{}\mathrm{0},\mathrm{1},\mathrm{2},\mathrm{3}\mathit{\}}$, i.e., all four nodes.
Region 0 convexly combines basis function *k* given by ${\mathcal{I}}_{\mathrm{0}}^{\mathrm{4}}=\mathit{\{}\mathrm{0},\mathrm{1},\mathrm{2}\mathit{\}}$
with GLL basis function *k* over the region $x\in \left[{x}_{G}^{\mathrm{4}}\right(\mathrm{0}),{x}_{G}^{\mathrm{4}}(\mathrm{1}\left)\right]$,
with convex combination function *α*(*x*).
*α*(*x*) is a quadratic polynomial having values 1, 0.306, and 0
at, respectively, the left side, middle, and right side of region 0.
The combination is defined so that at the left side, where $\mathit{\alpha}(-\mathrm{1})=\mathrm{1}$,
only the GLL basis function *k* is used;
at the right side, where $\mathit{\alpha}\left({x}_{G}^{\mathrm{4}}\right(\mathrm{1}\left)\right)=\mathrm{0}$,
only the offset nodal basis function *k* is used.
Region 2 is constructed by symmetry.
*α*(*x*) is determined to maximize the accuracy of the resulting ISL method
subject to the constraint of stability on the test problem.

## A2 Conservation when interpolating

The element-local conservation expressed in Eq. (15)
holds for the following reasons.
First, for any Islet basis for which $({n}_{\mathrm{p}}^{\mathrm{t}}{)}^{\mathrm{submin}}\ge {n}_{\mathrm{p}}^{\mathrm{v}}$,
the continuum field *f* is the same on each grid,
since the tracer grid can exactly represent polynomials of degree ${n}_{\mathrm{p}}^{\mathrm{v}}-\mathrm{1}$.
Second, for ${n}_{\mathrm{p}}^{\mathrm{v}}$ even, such as the standard ${n}_{\mathrm{p}}^{\mathrm{v}}=\mathrm{4}$,
if $({n}_{\mathrm{p}}^{\mathrm{t}}{)}^{\mathrm{submin}}={n}_{\mathrm{p}}^{\mathrm{v}}-\mathrm{1}$, the continuum field *f* is, in general, different on each grid,
but the integral of each over an element has the same value.
This is because for polynomial degree *d* odd and $f\left(x\right)\equiv {\sum}_{i=\mathrm{0}}^{d}{a}_{i}{x}^{i}$,
${\int}_{-\mathrm{1}}^{\mathrm{1}}f\left(x\right)\phantom{\rule{0.33em}{0ex}}\mathrm{d}x$ = ${\int}_{-\mathrm{1}}^{\mathrm{1}}\left(f\right(-x)+f(x\left)\right)/\mathrm{2}\phantom{\rule{0.33em}{0ex}}\mathrm{d}x$,
and $g\left(x\right)\equiv \left(f\right(-x)+f(x\left)\right)/\mathrm{2}$ has degree at least one less than *f*(*x*).
*g*(*x*) is exactly represented on the tracer grid
and, thus, the integral of *f*, which is the same as the integral of *g*,
is the same on the tracer grid as on the dynamics grid.
These two reasons assure that ℐ^{v→t} is conservative
for the Islet bases in Table A1 when ${n}_{\mathrm{p}}^{\mathrm{v}}=\mathrm{4}$.

## A3 Unstable and stable integration

Figure A1,
which has the format of the figures described in Sect. 4,
compares accuracy and stability between the natural and Islet bases.
The divergent flow, Gaussian hills IC, property preservation, TTPR, and long time steps are used.
To test the stability of the ISL linear advection operator using the Islet bases,
we run the verification problem for both 1 cycle and 100 cycles.
A curve's line pattern corresponds to basis type and number of cycles:
solid, Islet basis for 1 cycle; dashed, Islet basis for 100 cycles (or $\mathrm{12}\times \mathrm{100}=\mathrm{1200}$ d);
dash-dotted, natural basis for 1 cycle.
For each ${n}_{\mathrm{p}}^{\mathrm{t}}$, the *l*_{2} norm of the solution using the natural basis diverges with increasing resolution within the first cycle,
demonstrating that the basis leads to an unstable ISL linear advection operator.
In the case of ${n}_{\mathrm{p}}^{\mathrm{t}}=\mathrm{4}$, we see the start of the curve's divergence, but further element-grid refinement is needed to see the curve fully diverge.
In contrast, the curves for the Islet method with the Islet bases converge at order of accuracy 2.

Our primary interest is the advection equation and not the continuity equation.
However, to run tests with divergent flow
and property preservation enabled,
we need to solve the continuity equation for the air density
to provide the ** ρ** argument to CAAS.
This appendix describes a discretization of the continuity equation
that differs in only one factor from the discretization of the advection equation.

We start by following, e.g., Giraldo (1997, Sect. 3.3) and Bosler et al. (2019, Sect. 3).
The Reynolds transport theorem for a Lagrangian fluid parcel Ω(*t*)
containing a fluid having density *ρ*(** x**,

*t*) is

To balance clarity and concision,
we include or omit the position and time arguments depending on the context.
Suppose *ρ* satisfies the continuity equation,

In addition, consider a function *ϕ*(** x**,

*t*) that satisfies the advection equation,

Then, by Eqs. (B1), (B2), and (B3),

which, in turn, implies

for all times *t*_{0}, *t*_{1}.

Now we discretize Eq. (B4) in space.
First, we write an integral over a scalar *f* in terms of the element reference coordinate system.
Let Ω(*t*) be a quadrilateral element,
with Ω(*t*_{1}) the arrival Eulerian element.
Then,

Here, $\mathit{x}(\mathit{r},t)\equiv \mathit{X}(t;\mathit{x}(\mathit{r},{t}_{\mathrm{1}}),{t}_{\mathrm{1}})$
is the time-dependent map from the reference element to the manifold,
and satisfies Eq. (6).
*J*(** r**,

*t*) is the absolute value of the Jacobian determinant of this map, $\left|\mathrm{det}\right(\partial \mathit{x}(\mathit{r},t)/\partial \mathit{r}\left)\right|$.

Second, let the element have a tensor-product grid of GLL nodes, with *n*_{p} the number of 1D nodes.
We apply tensor-product GLL quadrature with this value of *n*_{p}
to the right hand side of Eq. (B5).
Then,

where *i* iterates over nodes in the element,
*r*_{i} is the *i*th GLL node,
and *w*_{i} is the *i*th GLL weight.
Let ${\mathit{x}}_{i}\left(t\right)\equiv \mathit{x}({\mathit{r}}_{i},t)$.

Third, we combine the quadrature approximation with a specific form of *f*,
*f*=*ρ**ϕ*, where *ϕ* satisfies Eq. (B3).
Let *ϕ*(** x**,

*t*

_{1}) be an element, GLL nodal, tensor-product

*n*

_{p}-basis function on an Eulerian arrival element. Because

*ϕ*is a nodal basis function with the same nodes as the element,

*ϕ*=1 at one node, say node

*k*, and is 0 at all other nodes

*i*≠

*k*.

*ϕ*need not be a natural GLL basis function; for example, it can be an Islet GLL basis function. Because $\mathrm{D}\mathit{\varphi}/\mathrm{D}t=\mathrm{0}$, $\mathit{\varphi}\left({\mathit{x}}_{i}\right(t),t)=\mathit{\varphi}\left({\mathit{x}}_{i}\right({t}_{\mathrm{1}}),{t}_{\mathrm{1}})$ for any time

*t*. Thus, the quadrature sum can be simplified:

Fourth, we apply quadrature to each side of Eq. (B4) to obtain the discretization

or,

Equation (B6) differs from Eq. (5) in the appearance
of the *density factor*, the quotient of the two Jacobian determinants.
In this equation, both *ρ* values and the denominator of the density factor are grid point values.
Note that this simple form depends on three different node sets being the same:
those for the element, the GLL quadrature, and the basis function *ϕ*.

Fifth, to complete the discretization, the discretization steps in Sect. 2.1.4 are applied to Eq. (B6) as they were to Eq.(5), but now with the additional factor present.

Since Eq. (B6) is not the focus of this article,
we omit detailed analysis of the discretization.
However, we note some basic facts.
Importantly, the Jacobian determinant of the Lagrangian element appears in the numerator
rather than the denominator of the density factor;
thus, the denominator is always well-behaved if the Eulerian grid has reasonable quality.
Second, for certain simple flows such as uniform translational flow on the plane,
*J* is constant in time,
and Eq. (B6) reduces to the discretization of the advection equation,
Eq. (5).
For general non-divergent flow, *J* is time dependent because of the spatial discretization error.
Third, only the density factor can make
the order of accuracy of the discretization Eq. (B6)
differ from that of Eq. (5).
Details concerning the manifold on which flow occurs and
the map from the reference element to the manifold
affect the order of accuracy.
In the verification problems in this article,
the density *ρ* is approximated using the Islet *n*_{p}=4 basis,
the manifold is the sphere,
and the map from the reference element to the sphere
is the isoparametric map of the finite element method for the *n*_{p}-basis.
The discretization of the continuity equation for *ρ* then has order of accuracy 2.

Code and scripts for the algorithms and figures presented in this paper are available at https://doi.org/10.5281/zenodo.5595499 (Bradley, 2021a) and, alternatively, https://github.com/E3SM-Project/COMPOSE/releases/tag/v1.1.2 (last access: 3 September 2021).
In this repository, read `methods/islet/readme.txt`

for further instructions;
in particular, this file points to detailed instructions to recreate the data in this article
in the file `methods/islet/figures/figs.tex`

.
These data are also available at https://doi.org/10.5281/zenodo.5595518 (Bradley, 2021b).

AMB developed the algorithms, wrote the software, ran the numerical studies, and wrote the manuscript, all with contributions from PAB and OG. PAB administers the project that funded most of this work.

The contact author has declared that none of the authors has any competing interests.

Sandia National Laboratories is a multimission laboratory managed and operated by National Technology and Engineering Solutions of Sandia, LLC., a wholly owned subsidiary of Honeywell International, Inc., for the US Department of Energy's National Nuclear Security Administration under contract DE-NA-0003525. This paper describes objective technical results and analysis. Any subjective views or opinions that might be expressed in the paper do not necessarily represent the views of the US Department of Energy or the United States Government. SAND2022-9879 J.

Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

We thank Mark A. Taylor for valuable discussions regarding material in this article and Topical Editor James Kelly and the anonymous reviewers for valuable suggestions that improved this article's exposition.

This work was supported by the US Department of Energy (DOE) Office of Science's Advanced Scientific Computing Research (ASCR) and Biological and Environmental Research (BER) Programs under the Scientific Discovery through Advanced Computing (SciDAC 4) ASCR/BER Partnership Program.

This paper was edited by James Kelly and reviewed by three anonymous referees.

Bosler, P. A., Bradley, A. M., and Taylor, M. A.: Conservative multimoment transport along characteristics for discontinuous Galerkin methods, SIAM J. Sci. Comput., 41, B870–B902, 2019. a, b

Bradley, A. M.: COMPOSE (v1.1.2): Methods for Islet paper, Zenodo [code], https://doi.org/10.5281/zenodo.5595499, 2021a. a

Bradley, A. M.: Methods data for Islet paper, Zenodo [data set], https://doi.org/10.5281/zenodo.5595518, 2021b. a

Bradley, A. M., Bosler, P. A., Guba, O., Taylor, M. A., and Barnett, G. A.: Communication-efficient property preservation in tracer transport, SIAM J. on Sci. Comput., 41, C161–C193, 2019. a, b, c, d, e, f, g, h

Bradley, A. M., Bosler, P. A., and Guba, O.: Islet: Interpolation semi-Lagrangian element-based transport, Geosci. Model Dev. Discuss., [preprint], https://doi.org/10.5194/gmd-2021-296, 2021. a, b, c

Chen, Y., Simon, K., and Behrens, J.: Extending legacy climate models by adaptive mesh refinement for single-component tracer transport: a case study with ECHAM6-HAMMOZ (ECHAM6.3-HAM2.3-MOZ1.0), Geosci. Model Dev., 14, 2289–2316, https://doi.org/10.5194/gmd-14-2289-2021, 2021. a

Dennis, J., Fournier, A., Spotz, W. F., St-Cyr, A., Taylor, M. A., Thomas, S. J., and Tufo, H.: High-resolution mesh convergence properties and parallel efficiency of a spectral element atmospheric dynamical core, Int. J. High Perform C., 19, 225–235, 2005. a

Dennis, J., Edwards, J., Evans, K. J., Guba, O., Lauritzen, P., Mirin, A. A., St-Cyr, A., Taylor, M. A., and Worley, P. H.: CAM-SE: A scalable spectral element dynamical core for the Community Atmosphere Model, Int. J. High Perform C., 26, 74–89, 2012. a, b

Dormand, J. R. and Prince, P. J.: A family of embedded Runge-Kutta formulae, J. Comput. Appl. Math., 6, 19–26, 1980. a

E3SM Project: Energy Exascale Earth System Model (E3SM), E3SM Project [computer software], https://doi.org/10.11578/E3SM/dc.20180418.36, 2018. a

Eastham, S. D. and Jacob, D. J.: Limits on the ability of global Eulerian models to resolve intercontinental transport of chemical plumes, Atmos. Chem. Phys., 17, 2543–2553, https://doi.org/10.5194/acp-17-2543-2017, 2017. a

Erath, C. and Nair, R. D.: A conservative multi-tracer transport scheme for spectral-element spherical grids, J. Comput. Phys., 256, 118–134, 2014. a

Field, C., Cheung, W., Dilling, L., Frumhoff, P., Greely, H., Hordequin, M., Hurrell, J., Light, A., Lin, A., MacMartin, D., McHenry, R., Moreno-Cruz, J., Ricke, K., Russell, L., Sagar, A., and Wennberg, P.: Reflecting Sunlight: Recommendations for Solar Geoengineering Research and Research Governance, The National Academies Press, Washington, DC, https://doi.org/10.17226/25762, 2021. a

Fischer, P. F. and Patera, A. T.: Parallel spectral element methods for the incompressible Navier-Stokes equations, in: Solution of Superlarge Problems in Computational Mechanics, Springer, 49–65, https://doi.org/10.1007/978-1-4613-0535-4_3, 1989. a

Giraldo, F. X.: Lagrange–Galerkin methods on spherical geodesic grids, J. Comput. Phys., 136, 197–213, 1997. a

Giraldo, F. X., Perot, J. B., and Fischer, P. F.: A spectral element semi-Lagrangian (SESL) method for the spherical shallow water equations, J. Comput. Phys., 190, 623–650, 2003. a

Golaz, J.-C., Caldwell, P. M., Van Roekel, L. P., Petersen, M. R., Tang, Q., Wolfe, J. D., Abeshu, G., Anantharaj, V., Asay-Davis, X. S., Bader, D. C., et al.: The DOE E3SM coupled model version 1: Overview and evaluation at standard resolution, J. Adv. Model Earth Sy., 11, 2089–2129, 2019. a

Golaz, J.-C., Roekel, L. P. V., Zheng, X., Roberts, A., Wolfe, J. D., Lin, W., Bradley, A., Tang, Q., Maltrud, M. E., Forsyth, R. M., Zhang, C., Zhou, T., Zhang, K., Zender, C. S., Wu, M., Wang, H., Turner, A. K., Singh, B., Richter, J. H., Qin, Y., Petersen, M. R., Mametjanov, A., Ma, P.-L., Larson, V. E., Krishna, J., Keen, N. D., Jeffery, N., Hunke, E. C., Hannah, W. M., Guba, O., Griffin, B. M., Feng, Y., Engwirda, D., Di Vittorio, A. V., Dang, C., Conlon, L. M., Chen, C.-C.-J., Brunke, M. A., Bisht, G., Benedict, J. J., Asay-Davis, X. S., Zhang, Y., Zhang, M., Zeng, X., Xie, S., Wolfram, P. J., Vo, T., Veneziani, M., Tesfa, T. K., Sreepathi, S., Salinger, A. G., Eyre, J. E. J. R., Prather, M. J., Mahajan, S., Li, Q., Jones, P. W., Jacob, R. L., Huebler, G. W., Huang, X., Hillman, B. R., Harrop, B. E., Foucar, J. G., Fang, Y., Comeau, D. S., Caldwell, P. M., Bartoletti, T., Balaguru, K., Taylor, M. A., McCoy, R. B., Leung, L. R., and Bader, D. C.: The DOE E3SM Model Version 2: Overview of the physical model and initial model evaluation, Earth and Space Science Open Archive [preprint], https://doi.org/10.1002/essoar.10511174.2, 2022. a, b, c, d

Guba, O., Taylor, M. A., Ullrich, P. A., Overfelt, J. R., and Levy, M. N.: The spectral element method (SEM) on variable-resolution grids: evaluating grid sensitivity and resolution-aware numerical viscosity, Geosci. Model Dev., 7, 2803–2816, https://doi.org/10.5194/gmd-7-2803-2014, 2014. a, b

Hannah, W. M., Bradley, A. M., Guba, O., Tang, Q., Golaz, J.-C., and Wolfe, J.: Separating Physics and Dynamics grids for Improved Computational Efficiency in Spectral Element Earth System Models, J. Adv. Model Earth Sy., 13, e2020MS002419, https://doi.org/10.1029/2020MS002419, 2021. a, b, c, d, e, f, g, h, i, j

Herrington, A. R., Lauritzen, P. H., Reed, K. A., Goldhaber, S., and Eaton, B. E.: Exploring a lower resolution physics grid in CAM-SE-CSLAM, J. Adv. Model Earth Sy., 11, 2019MS001684, https://doi.org/10.1029/2019MS001684, 2019. a, b

Kaas, E.: A simple and efficient locally mass conserving semi-Lagrangian transport scheme, Tellus A, 60, 305–320, 2008. a

Lauritzen, P. H.: Atmospheric Transport Schemes: Desirable Properties and a Semi-Lagrangian View on Finite-Volume Discretizations, in: Numerical Techniques for Global Atmospheric Models, edited by: Lauritzen, P. H., Jablonowski, C., Taylor, M. A., and Nair, R. D., vol. 80 of Lecture Notes in Computational Science and Engineering, Springer, Berlin, 185–250, 2011. a, b

Lauritzen, P. H. and Thuburn, J.: Evaluating advection/transport schemes using interrelated tracers, scatter plots and numerical mixing diagnostics, Q. J. Roy. Meteor. Soc., 138, 906–918, 2012. a

Lauritzen, P. H., Nair, R. D., and Ullrich, P. A.: A conservative semi-Lagrangian multi-tracer transport scheme (CSLAM) on the cubed-sphere grid, J. Comput. Phys., 229, 1401–1424, 2010. a, b

Lauritzen, P. H., Ullrich, P. A., Jablonowski, C., Bosler, P. A., Calhoun, D., Conley, A. J., Enomoto, T., Dong, L., Dubey, S., Guba, O., Hansen, A. B., Kaas, E., Kent, J., Lamarque, J.-F., Prather, M. J., Reinert, D., Shashkin, V. V., Skamarock, W. C., Sørensen, B., Taylor, M. A., and Tolstykh, M. A.: 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, 105–145, https://doi.org/10.5194/gmd-7-105-2014, 2014. a, b

Lauritzen, P. H., Ullrich, P. A., Jablonowski, C., Bosler, P. A., Calhoun, D., Conley, A. J., Enomoto, T., Dong, L., Dubey, S., Guba, O., Hansen, A. B., Kaas, E., Kent, J., Lamarque, J.-F., Prather, M. J., Reinert, D., Shashkin, V. V., Skamarock, W. C., Sørensen, B., Taylor, M. A., and Tolstykh, M. A.: 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, 105–145, https://doi.org/10.5194/gmd-7-105-2014, 2014. a, b, c

Lauritzen, P. H., Conley, A. J., Lamarque, J.-F., Vitt, F., and Taylor, M. A.: The terminator “toy” chemistry test: a simple tool to assess errors in transport schemes, Geosci. Model Dev., 8, 1299–1313, https://doi.org/10.5194/gmd-8-1299-2015, 2015. a, b

Lauritzen, P. H., Taylor, M. A., Overfelt, J., Ullrich, P. A., Nair, R. D., Goldhaber, S., and Kelly, R.: CAM-SE–CSLAM: Consistent Coupling of a Conservative Semi-Lagrangian Finite-Volume Method with Spectral Element Dynamics, Mon. Weather Rev., 145, 833–855, 2017. a, b

Lee, D., Lowrie, R., Petersen, M., Ringler, T., and Hecht, M.: A high order characteristic discontinuous Galerkin scheme for advection on unstructured meshes, J. Comput. Phys., 324, 289–302, 2016. a

Lee, H.-H., Bogenschutz, P., and Yamaguchi, T.: The Implementation of Framework for Improvement by Vertical Enhancement (FIVE) into Energy Exascale Earth System Model (E3SM), J. Adv. Model Earth Sy., 13, e2020MS002240, https://doi.org/10.1029/2020MS002240, 2020. a

Natarajan, H. and Jacobs, G. B.: An explicit semi-Lagrangian, spectral method for solution of Lagrangian transport equations in Eulerian-Lagrangian formulations, Comput. Fluids, 207, 104526, https://doi.org/10.1016/j.compfluid.2020.104526, 2020. a

Semakin, A. N. and Rastigejev, Y.: Optimized wavelet-based adaptive mesh refinement algorithm for numerical modeling of three-dimensional global-scale atmospheric chemical transport, Q. J. Roy. Meteor. Soc., 146, 1564–1574, 2020. a

Shampine, L. F. and Reichelt, M. W.: The Matlab ODE suite, SIAM J. Sci. Comput., 18, 1–22, 1997. a

Stensrud, D. J.: Parameterization Schemes: Keys to Understanding Numerical Weather Prediction Models, Cambridge University Press, ISBN 9780521126762, 2009. a

Tang, Q., Klein, S. A., Xie, S., Lin, W., Golaz, J.-C., Roesler, E. L., Taylor, M. A., Rasch, P. J., Bader, D. C., Berg, L. K., Caldwell, P., Giangrande, S. E., Neale, R. B., Qian, Y., Riihimaki, L. D., Zender, C. S., Zhang, Y., and Zheng, X.: Regionally refined test bed in E3SM atmosphere model version 1 (EAMv1) and applications for high-resolution modeling, Geosci. Model Dev., 12, 2679–2706, https://doi.org/10.5194/gmd-12-2679-2019, 2019. a

Ullrich, P. A., Lauritzen, P. H., and Jablonowski, C.: Some considerations for high-order “incremental remap”-based transport schemes: edges, reconstructions, and area integration, Int. J. Numer. Meth. Fl., 71, 1131–1151, 2013. a

Yamaguchi, T., Feingold, G., and Larson, V. E.: Framework for improvement by vertical enhancement: A simple approach to improve representation of low and high-level clouds in large-scale models, J. Adv. Model Earth Sy., 9, 627–646, 2017. a

- Abstract
- Introduction
- Overview of algorithms
- The Islet method
- Numerical results
- Conclusions
- Appendix A: Islet bases
- Appendix B: Discretization of the continuity equation
- Code and data availability
- Author contributions
- Competing interests
- Disclaimer
- Acknowledgements
- Financial support
- Review statement
- References

- Abstract
- Introduction
- Overview of algorithms
- The Islet method
- Numerical results
- Conclusions
- Appendix A: Islet bases
- Appendix B: Discretization of the continuity equation
- Code and data availability
- Author contributions
- Competing interests
- Disclaimer
- Acknowledgements
- Financial support
- Review statement
- References