Islet: Interpolation semi-Lagrangian element-based transport

. 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 efﬁcient 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 5 class of interpolation semi-Lagrangian (ISL) methods contains potentially extremely efﬁcient SL methods. We describe a ﬁnite-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 conﬁgurable number of ﬁnite-volume subcells per element; and a tracer grid supporting use of 10 Islet bases , with particular basis again conﬁgurable. This method provides extremely accurate tracer transport and excellent diagnostic values in a number of veriﬁcation problems.

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 30 tracer transport and physics parameterizations for a given, fixed dynamical core configuration.

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 35 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 (4) in terms of the solution at another time t 0 is q(x, t 1 ) = q(X(t 0 ; x, t 1 ), t 0 ).
X(t 0 ; x, t 1 ) is the solution of the ordinary differential equation d dt X(t; x, t 1 ) = u(X(t; x, t 1 ), t), Figure 1. Diagram of subelement grids. One spectral element (blue solid line outlining the full square) with dynamics (n v p = 4, black large circles), tracer (n t p = 8, small red circles), and physics (n f = 6, green dashed lines) subelement grids.
The 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 v p = 4, where we use 90 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 95 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 4/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. 100

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 t p , 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 v p , to trade between computational speed (maximum at n t p = n v p ) and higher accuracy (increasing with n t p ). In Fig. 1, the small red 105 circles are GLL nodes for n t p = 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 t p relative to n v p 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 v p ; the physics grid, n f ; and the tracer transport grid, n t p . For efficiency, EAM version 2 already uses two grids, one with n t p = n v p = 4 110 and the other with n f = 2. This article considers the case n t p > n v p . 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 115 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.

Outline
The rest of this article is structured as follows. The Islet method depends on a number of lower-level algorithms. Section 2 and 120 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.

Overview of algorithms
This section describes the details of low-level algorithms on which the Islet method depends, associated concepts, and related 125 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.

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 . 135 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. 140 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 Sects. 2.1.4 and 2.2, respectively.

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

150
An interpolation SL (ISL) method directly discretizes the tracer transport equation, Eq. (5): q(x i , t 1 ) = q(x i , t 0 ), where x i is an Eulerian arrival grid point andx i is its departure, generally off-grid, point at t 0 < t 1 . q has exact values only at grid points 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 155 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 160 method to second-order accuracy; however, Ullrich et al. (2013) describe a higher-order edge reconstruction that yields a thirdorder 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. 165 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 170 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.

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 175 maximize computational efficiency. We develop an ISL method that uses the Islet bases, summarized in Appendix A and derived in Bradley et al. (2021, Sect. 2 and 3) and a forthcoming article based on that material, to satisfy a necessary condition for stability. n v p = 4 isoparametric map interpolates the locations of the Lagrangian n t p = 6 tracer grid points from the Lagrangian n v p = 4 dynamics grid points. Third, the tracer data at time t 0 in the bottom-right Eulerian element is used to interpolate the tracer mixing ratio at t 0 to each green-starred Lagrangian point. Interpolation uses the 2D tensor-product n t p = 6 nodal basis. Fourth, these values are then the values in the upper-left Eulerian element at time t 1 .
Figure 2(b) shows two bases, the GLL basis, called the natural basis when a distinguishing name is needed, and the Islet 185 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 . 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. 2(b), the blue basis function is associated with the third node of six. The basis 190 is symmetric; basis function k ∈ {0, . . . , n p − 1} is the mirror image of basis function n p − k − 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.

195
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. The six 1D n t p = 6 element basis functions, a different color for each, for the natural GLL (dashed curves) and Islet GLL (solid curves) bases. In the bottom-right Eulerian element, a set of 1D basis functions is conceptually associated with each of the red horizontal and green vertical rectangles outlining an element side, and a 2D basis function over the element is a tensor product of one 1D basis function from each set.
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 200 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 205 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.

210
The tracer mixing ratio degrees of freedom are written q i ; 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 E, is written q(E). When needed, a superscript letter indicates on which grid the field is approximated: t, tracer; v, dynamics; f, physics.
Arithmetic between vectors applies entry-wise. For two vectors x and y, we write x y for element-wise multiplication and x y for element-wise division. 215 We need to define a map between a reference element with coordinates r ∈ [−1, 1] 2 and element e on the manifold with coordinates x, e.g., the sphere. Let this map be x = m e (r) and its inverse be r = m −1 e (x). Additionally, let e = E(x) return the element index for the element that contains e. In this work, if x is on the boundary of multiple elements, E(x) can be any of the corresponding element indices without affecting the solution.
Similarly, we need to map between (element-)local and global indices. Local indices are in the set {0, . . . , (n g p ) 2 − 1}, where 220 g ∈ {v, t}. Global index i is in element elm(i; g) ≡ i/(n g p ) 2 . The corresponding local index is lcl(i; g) ≡ i mod (n g p ) 2 . Given the local index j and element index e, the global index is glb(j; g) ≡ (n g p ) 2 e + j. Let E(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(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) .

225
Occasionally we write the tensor product in the 2D basis explicitly. Let a 1D element basis function set be ψ j , j ∈ 0, . . . , n p − 1.
Finally, we write the particular map x = m e (r) we use. The map uses just the corner nodes of an element; it is the isoparametric map for a bilinear basis. In element e, let corner e (j; v) be the global index of corner j of the element, here referenced 230 to the dynamics grid for concreteness. Let φ bilin j (r) 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).

The linear advection operator
Now we write the steps illustrated in Fig. 2. First, each grid point x v i 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 v p -basis isoparametric map followed by normalization: Third, an interpolated value from time t 0 is assigned to each advected tracer grid point:

Property preservation
Shape preservation and mass conservation on their own are not necessarily difficult to achieve. Rather, the hard problem in 245 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)  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 265 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.
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 i J i w i = 4π; 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.

290
Now we describe CAAS, independent of two sets of details: first, the means by which time-t 1 data are obtained; second, the index set E 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 ratioq i , and weight θ i . In addition, we are given a target total mass b over all i ∈ E. Finally, we have lower and upper bounds, q min i and q max i , on the target q i . These are the inputs to CAAS.
The output is a set of modified values q(E) that solves the 1-norm minimization problem This problem has a solution, i.e., the constraint set is non-empty, if and only if i∈E 300 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 ∈ E: Let the preliminary mass be Third, compute the final values for i ∈ E: As a check, note that carry out these steps.
CAAS is a CEDR because the global version, i.e., the problem in which E is over all grid points, can be implemented using a single global reduction. For details, see Bradley et al. (2019)[Sect. 6].

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 nonempty 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 andm < b. RECONSTRUCTSAFELY first tries to compute a one-norm-minimal relaxation to q max such that the maximum value max i∈E q max i 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(E) with uniform values that satisfy the mass conservation constraint.

Local and global problems
We use CAAS to solve problems at two levels: first, within each element separately; 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.

335
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 E.
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 340 are formulated as follows. Weight e is element e's area: A e ≡ i∈E(e) w i J i . Air density ρ e is the element's average density.
Let the element's mass be m e ≡ i∈E(e) w i J i ρ i . Then ρ e = m e /A e . The minimum mixing ratio is the minimum tracer mass over the element divided by the total mass: q min e ≡ i∈E(e) w i J i ρ i q min i /m e , and similarly for q max e . Finally, the initial valuē 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.

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 interelement, 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 max i − q min i 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.

Grid remap
Remapping data among multiple component grids is common in many applications of 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 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 365 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 370 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 I v→t interpolate a field in an element on the dynamics grid, represented by the natural GLL n v p -basis, to one on the tracer grid, represented by the Islet n t p -basis, using the natural n v p -basis functions to interpolate. Let I t→v do the opposite, using the Islet n t p -basis functions. 375 We can write these linear operators explicitly using the notation introduced in Sect.
where we have used the notation f e ≡ f (E(e)). I 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 380 each element e despite not being in all cases an L 2 projection.
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 385 determinant values from the dynamics grid, for two reasons. First, this operation conserves the area of the element,

Direct stiffness summation
In a spectral element method, most operations are performed independently in each element, often leading to discontinuities where e is the sum over all elements. In our use of the generalized DSS, σ = ρ andȳ =q. In the standard DSS, σ is all one andȳ = ρ 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 q continuous but leave ρ discontinuous and, in particular, unmodified from its original value.

405
3 The Islet method 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 410 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 f or ∆q is sometimes called a tendency. Third, the remapped tracer tendencies are added to the tracer-grid tracers. Fourth, the 415 tracers are advected on the tracer grid. Fifth, tracer states are remapped to the physics and, optionally, dynamics grids.

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 J ρ n+1 from the dynamics grid to the tracer grid. In each element e, compute 420 This grid transfer conserves J e ρ n+1 e by Eq. (15). J t (ρ t ) n+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).

425
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 (q f e ) 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 val-430 ues of the tracer-grid mixing ratio state (q t e ) n . This final extrema update assures that if ∆(q f e ) n = 0, then (q t e ) 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 ρ t e . The operator uses this quantity and ρ f e ∆(q f e ) n to compute ρ t e ∆(q t e ) n .
Step::tendency-f2t::CAAS. In each element e, apply CAAS, described in Sect. 2.2.1, to (q t e ) n ≡ (q t e ) n + ∆(q t e ) n , using the weight vector θ t e ≡ w t e J t e and air density (ρ t e ) 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 440 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, (q t ) n is discontinuous across element boundaries. Neither order of accuracy nor 445 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 (q t ) n , using (ρ t ) n , to obtain the continuous field (q 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 450 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 I 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 455 are normalized to the sphere. For n v p = 4, the n t p -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 460 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 (q 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.

465
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 470 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.

475
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 θ ≡ w t J t ; , 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 480 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.

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

490
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, I t→v interpolates the tracer-grid mixing ratio to the n v p -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 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.

500
If n t p = n v p , 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.

Numerical results
This section presents results for a number of verification problems. Except in Sect. 4.3, the equation is the sourceless advection 505 equation, Eq. (4). Two-dimensional, time-dependent flow, u(x, t), is prescribed on the sphere.
In most figures, we show results for n t p = 4, 6, 8, 9, 12. The n t p = n v p = 4 case provides a reference because it does not use the TTPR algorithms. The n t p = 6 case has the smallest value of n t p providing order of accuracy (OOA) greater than 2, in this case 4. The n t p = 8 case provides OOA 5 and has four times as many nodes as the n t p = 4 basis in two dimensions. The n t p = 9 case provides OOA 6. Finally, the n t p = 12 case provides OOA 8 and has four times as many nodes as the n t p = 6 basis.

510
Tests follow the procedures detailed in . 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 515 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 nondivergent 520 one and a divergent one. Each prescribes a flow that lasts for T = 12 days and such that at time T , the exact solution is the same as the initial condition. This 12-day prescribed flow can be run for multiple cycles to lengthen the simulation. The nondivergent flow creates a filament of maximum aspect ratio at time T /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 525 the behavior of a transport method on tracers having various degrees of smoothness is important because atmosphere tracers can be smooth or nonsmooth.
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 530 tests for n t p > 4 use TTPR with n v p = 4 unless we state otherwise. Because n v p = 4 in these verification problems, the n t p = 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 535 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 t p = 4, except where noted, while using CAAS-point for all other n t p values, thus maximizing the accuracy obtained with n t p = 4 to provide the best baseline (n t p = n v p = 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 540 occurs when transport is coupled to a full dynamical core. For the nondivergent 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 ∆m/a 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 545 does not occur in these verification problems. The algorithm for ρ v has OOA 2 for n v p = 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 v p ×n v p tensor grid of GLL nodes and thus n v p − 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 6n e steps per T ; short, 30n e .

550
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 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. Fig. 15 show convergence plots, and we explain the format of these figures here. A curve's marker corresponds to n t p , as listed in the legend. To maximize font size and minimize notational clutter in figures, in figures we use n p rather than n t p and omit "="; e.g., "n t p = 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 v p = 4 correspond to the resolution 560 360 • 4 cube faces · 1 cube face 5 elements · 1 element (4 − 1) intervals = 6 • interval .

Figures 3-9 and
The y-axis is log 10 of the relative error, with the norm or norms indicated in each figure

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 v p = 4 and I 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 570 than would occur in a more realistic problem. Of particular concern is the symmetry of the flow in time around the midpoint time, 6 days. 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 nondivergent flow, the Gaussian hills IC, and the long time step. The midpoint reference solution is computed using one 6-day step and the natural GLL basis. The time integrator's relative error tolerance is set to 10 −14 rather 575 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-day 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 580 almost exactly parallel lines for resolution at least as fine as 0.75 • . In addition, the empirical OOA for n t p = 12 is omitted because the curves nearly overlap those for n t p = 9. For fine enough resolution, all curves should have empirical OOA limited to 4, but for n t p = 8, the full convergence regime is not reached in this plot, thus giving OOA 5 for n t p = 8. The n t p = 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 t p = 9 and 12. The n t p = 9 curve shows at coarse resolution higher OOA, governed by the spatial error, and then a drop to 4 as the temporal 585 error becomes dominant. Importantly, the half-and full-cycle errors are very close for each value of n t p , demonstrating that the endpoint error metrics are valid measurements for the Islet method.

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 590 cases.  Fig. 1, bottom right, of TR14). The n t p = 12 Islet scheme with the long time step, Fig. 4, is approximately 595 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 t p ≥ 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 t p ≥ 6. At 3 • resolution, no 600 method other than HEL-ND-CN1.0 provides l 2 norm below 10 −2 ; the Islet method does for n t p ≥ 8 with the long time step and n t p ≥ 10 (only n t p = 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   Fig. 13 for an illustration of the filamentary structure at the simulation midpoint, 620 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, τ ∈ [0.1, 1]. The diagnostic is then this area divided by the correct area, which for any nondivergent flow is the area at the initial time. The perfect diagnostic value is 100% for all τ ∈ [0.1, 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 625 limits are tighter with increasing n t p . 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 630 n t p ≥ 8, the resolution 1.5 • , and on the tracer grid are among the best of those in TR14 at the resolution 1.5 • . The test uses the nondivergent flow with two ICs. The diagnostics assess preservation of the nonlinear correlation of two 635 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

Mixing diagnostic
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 t p and the time step length are printed For the long time step, the Islet method gives at least as small a value for n t p ≥ 6; for the short time step, n t p ≥ 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 650 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 ≤ 5 × 10 −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 655 usable, unlike HEL-ND, and are designed to preserve tracer correlations exactly. Other than the HEL methods, the next best value is 4.80 × 10 −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 t p ≥ 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. 660 We observe that in both the filament and mixing diagnostics of TS12, n t p ≥ 6 gives excellent results; n t p = 12, nearly perfect. Figures 13 and 14 show solution quality using latitude-longitude images and further underscore these observations. The problem is the nondivergent flow with the slotted cylinders IC, at resolutions 1.5 • and 0.75 • , with the long and short time steps.

Slotted cylinders
The text in the individual images in Fig. 14 provides normwise accuracy at the end of one cycle and deviation from the initial extrema. φ min ≥ 0 and φ max ≤ 0 are consistent with no global extrema violation. Compare Fig. 13 with Figs. 7-10 in TR14 665 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 t p ≥ 6.

Source term
Now we move to verification problems that include a source term: f = 0 in Eq. (3). The legend describes the dynamics-grid resolution and the time step length. The prescribed verification problem is the nondivergent flow with cosine bells IC. Property preservation is on. The x-axis is τ , the mixing ratio threshold. The y-axis is the percent area having mixing ratio at least τ relative to that at the initial time.

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 x is position on the sphere. A tendency ∆m is applied to m on the physics grid: ∆m(x, t) ≡ −[cos(2π(t + ∆t)/T ) − cos(2πt/T )]s(x, t)/2, so that the exact solution is m(x, t) = (1 − cos(2πt/T ))s(x, t)/2 and, in particular, m(x, T /2) = 675 s(x, T /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 /2, on the dynamics grid as usual, as in Fig. 3. Figure   15 shows the results. We run this test with n f = n t p (dash-dotted lines) and n f = 2 (dashed lines). The solid lines show the error in s(x, T /2) as a reference. The dotted line provides the OOA-2 reference. We see that when n f = n t p , the errors are nearly the same as those for s(x, T /2); for n t p = n f = 4, the curves overlap at the resolution of the plot. As one expects, when 680 n f = 2, the error in m(x, t) is much larger than in s(x, T /2), but the OOA remains 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 = 1, 2, that interact according to chemical kinetic equations that are nonlinear in one of them: DX i (y, t)/Dt = f i (y, X 1 , X 2 ), where y is the spatial coordinate on the sphere. (We use y here to avoid overlap with the 685 species symbol, 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 690 constant,X T . The source terms have this property, too, since they model chemical reactions. Thus, in the exact solution,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 −X T at time t normalized byX T , and c ∞ (t), the same but for the l ∞ norm.
As explained in the context of equation 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, 695 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 705 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.
710 Figure 16 shows the diagnostic values for the case of nondivergent flow, 1 • dynamics-grid resolution, and a 30-minute time step, where these configuration details are prescribed in TC15. For the n t p = 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 days) of the flow, but it is useful to view it over multiple cycles. Figure 16 shows ten cycles on the x-axis, for a total of 120 days. 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 f = n t p . 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 t p = 4 case with CAAS-point is greatly improved relative to the Eulerian spectral element results shown in Fig. 7 of TC15, even after ten cycles instead of the one cycle shown in that figure. In Fig. 7 of TC15, c 2 ≈ 10 −2 and c ∞ > 10 −1 at the end of one cycle, compared with approximately 10 −7 and 10 −5 , respectively, for n t p = 4 at 720 the end of ten cycles. For n t p > 4, the growth in error is very small, with c 2 < 10 −10 through ten cycles. To illustrate what these diagnostics measure, Figure 17 shows latitude-longitude images of the monatomic tracer at the end of the first cycle for n t p = n f = 4 with CAAS-CAAS (left) and n t p = n f = 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.

725
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 4×10 −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 −X T )/X T at the end of the first cycle. The correct Toy chemistry: 1 • , ∆t 30min, nondivergent flow 0.0e+00 5.0e-07 1.0e-06 1.5e-06 2.0e-06 2.5e-06 3.0e-06 3.5e-06 4.0e-06 4.5e-06 -7.5e-02 -3.7e-02 0.0e+00 3.7e-02 7.5e-02 n p 8 n f 8 CAAS-point min -3.4e-13 max 2.6e-12 Toy chemistry: 1 • , ∆t 30min, nondivergent flow -2.6e-12 -1.3e-12 0.0e+00 1.3e-12 2.6e-12 Figure 18. Same as Fig. 17, but now the images are of (XT −XT )/XT . 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 t p = 8 at the end of the first cycle in Fig. 16.

Conclusions
We described the Islet method, a property preserving tracer transport method to support a three-grid atmosphere model, one 735 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.

740
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 745 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 t p 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 755 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.
Code and data availability. Code and scripts for the algorithms and figures presented in this paper are available at https://doi.org/10.5281/ zenodo.5595499 and, alternatively, https://github.com/E3SM-Project/COMPOSE/releases/tag/v1.1.2. In this repository, read methods/-We construct a basis for each n p ∈ {4, . . . , 13}. Each basis is described by a set of support node lists, one I np r 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 I np r is given by an offset, which is the first index in the list, and its size, |I np r |. For example, n p = 7, offset 2, and size 4 correspond to support nodes I 7 r = {2, 3, 4, 5}. In addition, the basis is symmetric, meaning basis function φ np i (x) = φ np np−1−i (−x). Thus, first, support nodes are specified for regions 0 through n p /2 − 1, and the support nodes for the remaining regions are determined by symmetry. Second, if n p is even, then the 785 middle region, r = n p /2 − 1, has support nodes I np r 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 submin p − 1 if there is no property preservation step and 790 the initial condition is C n submin p −2 .
If n submin p = n 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 795 nodal bases, most of them offset nodal bases, for n p ∈ {5, . . . , 13}. 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 I 1 = {0, 1, 2, 3}, i.e., all four nodes. Region 0 convexly combines basis function k given by I 4 0 = {0, 1, 2} with GLL basis function k over the region x ∈ [x 4 G (0), x 4 G (1)], with convex combination function α(x). α(x) is a quadratic 800 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 α(−1) = 1, only the GLL basis function k is used; at the right side, where α(x 4 G (1)) = 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 805
The element-local conservation expressed in Eq. (15) holds for the following reasons. First, for any Islet basis for which (n t p ) submin ≥ n v p , the continuum field f is the same on each grid, since the tracer grid can exactly represent polynomials of degree n v p − 1. Second, for n v p even, such as the standard n v p = 4, if (n t p ) submin = n v p − 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 −1 (f (−x) + f (x))/2 dx, and g(x) ≡ (f (−x) + f (x))/2 has degree at least one 810 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, offsets {0, 0, 0, 0, 0, 1} Table A1. Islet GLL nodal subset bases. Each row provides a formula for the row's np value. Columns are np, order of accuracy (OOA), the support sizes |Ir| for each region ordered left to middle, and the supports. For offset nodal subset bases, supports are given by offsets. For general nodal subset bases, supports are given by nodal subsets, again ordered from left region to middle. The case np = 4 is described in Appendix A1. In all cases, the support points are GLL nodes.
is the same on the tracer grid as on the dynamics grid. These two reasons assure that I v→t is conservative for the Islet bases in Table A1 when n v p = 4. Figure A1, which has the format of the figures described in Sect. 4, compares accuracy and stability between the natural 815 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 12 × 100 = 1200 days); dash-dotted, natural basis for 1 cycle. For each n t p , 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 t p = 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. The x-axis is the average dynamics-grid grid point spacing at the equator in degrees for the quasiuniform cubed-sphere grid. The y-axis is log 10 l2 relative error. A curve's line pattern corresponds to the basis type and number of cycles, as listed in the top legend. A curve's marker corresponds to n t p , as listed in the bottom legend.

Appendix B: Discretization of the continuity equation
Our primary interest is the advection equation and not the continuity equation. However, to run tests with divergent flow and 825 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.
Here, x(r, t) ≡ X(t; x(r, t 1 ), t 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, |det(∂x(r, t)/∂r)|.
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 ith GLL node, and w i is the ith GLL weight. Let x i (t) ≡ x(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 855 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 Dφ/Dt = 0, φ(x i (t), t) = φ(x i (t 1 ), t 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 ρ(x k (t 1 ), t 1 ) J(r k , t 1 ) w k = ρ(x k (t 0 ), t 0 ) J(r k , t 0 ) w k , or ρ(x k (t 1 ), t 1 ) = J(r k , t 0 ) J(r k , t 1 ) ρ(x k (t 0 ), t 0 ).
Eq. (B6) differs from Eq. (5) in the appearance of the density factor, the quotient of the two Jacobian determinants. In this 865 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 870 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 nondivergent 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).

875
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.
Author contributions. AMB developed the algorithms, wrote the software, ran the numerical studies, and wrote the manuscript, all with