Articles | Volume 15, issue 16
Geosci. Model Dev., 15, 6285–6310, 2022
Geosci. Model Dev., 15, 6285–6310, 2022
Development and technical paper
16 Aug 2022
Development and technical paper | 16 Aug 2022

Islet: interpolation semi-Lagrangian element-based transport

Islet: interpolation semi-Lagrangian element-based transport
Andrew M. Bradley, Peter A. Bosler, and Oksana Guba Andrew M. Bradley et al.
  • Sandia National Laboratories, Albuquerque, New Mexico, USA

Correspondence: Andrew M. Bradley (


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.

1 Introduction

An atmosphere model has two parts: the dynamical core, which computes resolved fluid flow and resolved thermodynamics (Lauritzen2011); and the subgrid parameterizations, which compute unresolved chemistry and physics processes (Stensrud2009). 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 Project2018) 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

(1) ρ t + ( ρ u ) = 0 ,

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

(2) ( ρ q ) t + ( ρ q u ) = ρ f ,

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

(3) D q D t q t + u q = f .

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

(4) D q D t = 0 .

At position x and time t1, the exact solution of Eq. (4) in terms of the solution at another time t0 is

(5) q ( x , t 1 ) = q ( X ( t 0 ; x , t 1 ) , t 0 ) .

X(t0;x,t1) is the solution of the ordinary differential equation

(6) d d t X ( t ; x , t 1 ) = u ( X ( t ; x , t 1 ) , t ) ,

with the initial condition X(t1;x,t1)=x. In a semi-Lagrangian method, X(t1;x,t1)=x is called the arrival point and X(t0;x,t1) is called the departure point. In this article, our focus is on the two-dimensional equations on the sphere. Thus, X and u each have two horizontal components and no vertical component. In various settings, the density variable ρ 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 t0, 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 Thuburn2012; Lauritzen2011; 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.

Figure 1Diagram of subelement grids. One spectral element (blue solid line outlining the full square) with dynamics (npv=4, black large circles), tracer (npt=8, small red circles), and physics (nf=6, green dashed lines) subelement grids.


In the horizontal direction, each quadrilateral element has an np×np 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. np is the number of subgrid nodes in a 1D basis function; the degree of the polynomial is np−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 np=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 np. For example, nodes for np=6, 8, and 12 have minimum distances respectively 2.4, 4.3, and 10.0 times smaller than nodes for np=4. Thus, developers of a dynamics solver choose a value for np 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 npv=4, where we use v to refer specifically to the dynamics solver's value of np. 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 nf×nf subgrid of finite-volume quadrilateral subcells. In Fig. 1, the dashed green lines outline the 6×6 physics subgrid, with nf=6. EAM version 2 uses nf=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.

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 npt, where we use t to refer specifically to the transport solver's value of np, a parameter of the tracer transport module and independent of npv, to trade between computational speed (maximum at npt=npv) and higher accuracy (increasing with npt). In Fig. 1, the small red circles are GLL nodes for npt=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 npt relative to npv 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, npv; the physics grid, nf; and the tracer transport grid, npt. For efficiency, EAM version 2 already uses two grids, one with npt=npv=4 and the other with nf=2. This article considers the case npt>npv. 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.

2 Overview of algorithms

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 t1 is advected backward or forward in time to its generally off-grid departure point at time t0. The resulting Lagrangian departure grid then samples the field at time t0 and, by some means, uses these data to construct the advected field at time t1.

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 Nair2014). 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(xi,t1)=q(x¯i,t0), where xi is an Eulerian arrival grid point and x¯i is its departure, generally off-grid, point at t0<t1. q has exact values only at grid points xi; thus, evaluating q(x¯i,t0) 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 L2 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 2Diagram of the linear advection operation in a finite-element interpolation semi-Lagrangian time step. (a) Six spectral elements (blue solid lines) with an npt=6 tracer grid (red small circles). The black large circles in the upper-left element show the npv=4 dynamics grid. The upper-left Eulerian arrival element is advected backward in time to the right and down, resulting in the Lagrangian departure element (blue dashed outline, open circle dynamics grid points, green tracer grid points). The domain of dependence of each green starred tracer grid point is the bottom right Eulerian element, with interpolation grid points (red small circles) associated with this domain of dependence outlined in black. (b) The six 1D npt=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.


Figure 2a illustrates the linear advection step in the Islet method. First, the upper-left element is advected backward in time from time t1 to time t0, using velocity data on the npv=4 dynamics grid to advect each dynamics grid point. Second, the npv=4 isoparametric map interpolates the locations of the Lagrangian npt=6 tracer grid points from the Lagrangian npv=4 dynamics grid points. Third, the tracer data at time t0 in the bottom-right Eulerian element are used to interpolate the tracer mixing ratio at t0 to each green starred Lagrangian point. Interpolation uses the 2D tensor-product npt=6 nodal basis. Fourth, these values are then the values in the upper-left Eulerian element at time t1.

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{0,,np-1} is the mirror image of basis function np-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. 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 xi. 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 qi 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 q(ℰ). When needed, a superscript letter indicates on which grid the field is approximated: t, tracer; v, dynamics; and f, physics.

Arithmetic between vectors applies entry-wise. For two vectors, x and y, we write xy for element-wise multiplication and xy for element-wise division.

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=me(r) and its inverse be r=me-1(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,,(npg)2-1}, where g{v,t}. Global index i is in element elm(i;g)i/(npg)2. The corresponding local index is lcl(i;g)imod(npg)2. Given the local index j and element index e, the global index is glb(j;g)(npg)2e+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 qeq(ℰ(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 rlcl(i;g).

Occasionally, we write the tensor product in the 2D basis explicitly. Let a 1D element basis function set be ψj, j0,,np-1. Let r(r1,r2). Then, ϕi(r)ψi/np(r1)ψimodnp(r2). Figure 2b shows the npt=6 basis functions ψj.

The function normalize(x)x/x2.

Finally, we write the particular map x=me(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 cornere(j;v) be the global index of corner j of the element, here referenced to the dynamics grid for concreteness. Let ϕjbilin(r) be the bilinear basis functions, i.e., the GLL basis with np=2. Then,

(7) x ( r ) = normalize j = 0 3 x corner ( j ; v ) v ϕ j bilin ( r ) .

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 xiv is advected from time t1 to time t0 to give

(8) x ¯ i v X ( t 0 ; x i v , t 1 ) .

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 npv-basis isoparametric map followed by normalization:

(9) x ¯ i t normalize j E v ( elm ( i ; t ) ) x ¯ j v ϕ lcl ( j ; v ) v ( r lcl ( i ; t ) t ) .

Third, an interpolated value from time t0 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.


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

(12) θ ¯ i - 1 1 - 1 1 ϕ lcl ( i ) ( r ) m elm ( i ) ( r ) r d 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 np-basis GLL quadrature. Let


Because ϕlcl(i)(rlcl(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

(13) J i m elm ( i ) ( r lcl ( i ) ) r

and wi as the product of the two 1D basis function integrals. In practice, HOMME modifies each node's Jacobian determinant Ji by multiplying it by a constant very close to 1 such that the global sum iJiwi=4π; that is, the sum of all the products of wi and adjusted Ji 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-t1 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 t1 and grid point i, we are given air density ρi, preliminary tracer mixing ratio 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, qimin and qimax, on the target qi. These are the inputs to CAAS. The output is a set of modified values q(ℰ) that solve the 1-norm minimization problem

miniE|qi-q¯i|subject toiEθiρiqi=b(mass conservation)qiminqiqimaxforiE(shape preservation).

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

(14) i E θ i ρ i q i min b i E θ i ρ i q i max .

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 m^=b, then set qi=q^i for i∈ℰ and terminate. Otherwise, suppose m^<b. Second, compute the capacity:


Third, compute the final values for i∈ℰ:


As a check, note that


If b-m^c, then qiqimax for i∈ℰ; these conditions hold if and only if Eq. (14) does. The case 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 m^<b. ReconstructSafely first tries to compute a one-norm-minimal relaxation to qmax, such that the maximum value maxiEqimax 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: AeiE(e)wiJi. Air density ρe is the element's average density. Let the element's mass be meiE(e)wiJiρi. Then, ρe=me/Ae. The minimum mixing ratio is the minimum tracer mass over the element divided by the total mass: qeminiE(e)wiJiρiqimin/me, and similarly for qemax. Finally, the initial value q¯e is q¯eiE(e)wiJiρiqi/me. 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 qimax-qimin 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 Rastigejev2020); 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 npv-basis to one on the tracer grid represented by the Islet npt-basis, using the natural npv-basis functions to interpolate. Let t→v do the opposite using the Islet npt-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{0,,(npt)2-1},


where we have used the notation fef(ℰ(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 L2 projection. Consider fet=Ivtfev. Then,

(15) i E t ( e ) w i t f i t = i E v ( e ) w i v f i v .

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 Jt using Eq. (13). On the tracer grid, rather than use Eq. (13), we interpolate the Jacobian determinant values from the dynamics grid,

(16) J e t = I v t J e v ,

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 Patera1989; 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 wiJi 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 δijgrid 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 y¯ is written:

(17) y i = DSS ( θ , σ , y ¯ , i ) e j E ( e ) δ i j grid θ j σ j y ¯ j e j E ( e ) δ i j grid θ j σ j ,

where e is the sum over all elements. In our use of the generalized DSS, σ=ρ and y¯=q¯. In the standard DSS, σ is all one and y¯=ρq¯. These two use cases give the same value for qi 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.

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 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, ΔqfΔ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 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 Jρn+1 from the dynamics grid to the tracer grid. In each element e, compute


This grid transfer conserves Jeρen+1 by Eq. (15). Jt(ρt)n+1 is not continuous across element boundaries because Jt 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 Jt according to Eq. (16).

Step::tendency-f2t. Map the tracer tendency Δqn 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 (qef)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 (qet)n. This final extrema update assures that if Δ(qef)n=0, then (q¯et)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 L2 projection is the Islet basis. From the previous step's application of Step::density-d2t, we have ρet. The operator uses this quantity and ρefΔ(qef)n to compute ρetΔ(qet)n.

Step::tendency-f2t::CAAS. In each element e, apply CAAS, described in Sect. 2.2.1, to (q¯et)n(qet)n+Δ(qet)n, using the weight vector θetwetJet and air density (ρet)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, (q¯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 (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 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 npv=4, the npt-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 (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.

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

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 θwtJt; b=iwitJit(ρit)n(q^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 (qt)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 (qt)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 npv-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 (qv)n+1. Third, the standard DSS is applied to (ρv)n+1(qv)n+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 npt=npv, 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.

4 Numerical results

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 npt=4, 6, 8, 9, 12. The npt=npv=4 case provides a reference because it does not use the TTPR algorithms. The npt=6 case has the smallest value of npt providing order of accuracy (OOA) greater than 2, in this case, 4. The npt=8 case provides OOA 5 and has four times as many nodes as the npt=4 basis in two dimensions. The npt=9 case provides OOA 6. Finally, the npt=12 case provides OOA 8 and has four times as many nodes as the npt=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/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 C1 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 Prince1980; Shampine and Reichelt1997) with a tight tolerance (10−8 with an exception noted later) to integrate trajectories very accurately. All tests for npt>4 use TTPR with npv=4 unless we state otherwise. Because npv=4 in these verification problems, the npt=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 npt=4, except where noted, while using CAAS–point for all other npt values, thus maximizing the accuracy obtained with npt=4 to provide the best baseline (npt=npv=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 Δm/atotal to each grid point, where Δm is the global mass discrepancy after advection, and atotal is the total area of the grid. Negative density does not occur in these verification problems. The algorithm for ρv has OOA 2 for npv=4.

We use a quasiuniform, equiangular cubed-sphere element grid. Let a cube face of the cubed-sphere element grid have ne×ne elements. Each element has an npv×npv tensor grid of GLL nodes and, thus, npv-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 6ne steps per T; short, 30ne. 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 39 and 15 show convergence plots, and we explain the format of these figures here. A curve's marker corresponds to npt, as listed in the legend. To maximize font size and minimize notational clutter in figures, in figures, we use np rather than npt , and omit “=”; e.g., “npt=8” is written “np 8”. Additionally, the legend omits “np” 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, ne=5 and npv=4 correspond to the resolution

3604cube  faces1cube  face5elements1element(4-1)intervals=6interval.

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 C1 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 npv=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.

Figure 3Comparison of relative errors calculated at the test simulation's midpoint time of 6 d (1/2 cycle, dashed lines) and endpoint time of 12 d (1 cycle, solid lines). Each number at the right side of the plot is the empirical OOA computed using the final two points of the one-cycle result.


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 npt=12 is omitted because the curves nearly overlap those for npt=9. For fine enough resolution, all curves should have empirical OOA limited to 4, but for npt=8, the full convergence regime is not reached in this plot, thus giving OOA 5 for npt=8. The npt=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 npt=9 and 12. The npt=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 npt, 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.

Figure 4Accuracy diagnostic. Compare with Figs. 1 and 2 in TR14.


Figure 5Accuracy diagnostic. Compare with Figs. 1 and 2 in TR14.


Figure 6Accuracy diagnostic. Compare with Fig. 3 in TR14.


Figure 7Accuracy diagnostic. Compare with Fig. 3 in TR14.


Figure 8Accuracy diagnostic. Compare with Fig. 16 in TR14.


Figure 9Accuracy diagnostic. Compare with Fig. 16 in TR14.


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 npt=12 Islet scheme with the long time step, Fig. 4, is approximately three times more accurate than HEL-ND-CN1.0 in the l2 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 l2 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 npt8. 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 npt6. At 3 resolution, no method other than HEL-ND-CN1.0 provides l2 norm below 10−2; the Islet method does for npt8 with the long time step, and npt10 (only npt=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 C1 IC. The horizontal dash-dotted line provides the relative l2-error-norm value of 0.033 by which the “minimal resolution” diagnostic value is determined; the coordinate of the intersection between the l2-norm curve and this reference line is the value. A larger value is better. For example, with a long time step, for npt=8, this value is a little coarser than 3; for npt=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.

Figure 10Filament diagnostic following Sect. 3.3 of TS12. Compare with Fig. 5 in TR14. The top row shows the diagnostic measured on the npv=4 dynamics grid; the bottom row on the tracer grid. The legend describes the dynamics-grid resolution and the time step length. The prescribed verification problem is the non-divergent 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.


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, τ[0.1,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 τ[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 limits are tighter with increasing npt.

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 npt8, 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.

Figure 11Mixing diagnostic following Sect. 3.5 of TS12. Compare with Figs. 11–14 in TR14. This figure shows results for dynamics-grid resolution of 1.5. lo is exactly 0 in all cases because shape preservation is on, and so is not shown. See text for further details.


Figure 12Same as Fig. 11 but with dynamics-grid resolution 0.75.


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. npt and the time step length are printed in each plot. The diagnostic values lr and lu, 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 lr measures this type of error; smaller is better. In Figs. 11–14 in TR14, the smallest value of lr at 1.5 among the property-preserving methods is 2.15×10-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 npt6; for the short time step, npt8.

Values outside the triangle are overshoots, possible only if the method is not strictly shape preserving. Since our method is, this diagnostic value, lo, 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 lu and, again, smaller is better. This diagnostic is more difficult to compare than lr because very dissipative methods tend to have a large value of lr and, consequently, a very small value for lu. In contrast, a very accurate method, for which lr is small, can have a larger lu value than a very dissipative method. One means of comparison is to consider the best lu values among methods that obtain, say, lr5×10-4. In Figs. 11–14 in TR14, the smallest value of lu 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 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 npt8. However, even with the constraint on lr, comparison is not straightforward as the UCISOM methods are not strictly shape preserving and so have lo>0.

4.2 Slotted cylinders

We observe that in both the filament and mixing diagnostics of TS12, npt6 gives excellent results; npt=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. ϕmin0 and ϕmax0 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 l2, linf, ϕmin, and ϕmax are at least as good as those in Fig. 7 of TS12 for npt6.

Figure 13Images of the slotted cylinders IC advected by the non-divergent flow at the simulation's midpoint. Each column corresponds to a spatial resolution and time step length configuration, as stated at the top of each column. Each row corresponds to a particular value of npt, as stated in the text at the top right of each image. We omit npt=12 results for the 0.75 resolution because they are essentially identical at the resolution of the figure to the npt=8 images.


Figure 14Same as Fig. 13 but for the simulation final point. Error measures are printed at the bottom left of each image; see text for details.


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 q1=s and a manufactured tracer q2=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)=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 nf=npt (dash-dotted lines) and nf=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 nf=npt, the errors are nearly the same as those for s(x,T/2); for npt=nf=4, the curves overlap at the resolution of the plot. As one expects, when nf=2, the error in m(x,t) is much larger than in s(x,T/2), but the OOA remains 2.

Figure 15Verification of the remap of tendencies from physics to tracer grids and state from tracer to dynamics grids. See Sect. 4.3 for a description of the problem.


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, qi=Xi, i=1,2, that interact according to chemical kinetic equations that are nonlinear in one of them: DXi(y,t)/Dt=fi(y,X1,X2), where y is the spatial coordinate on the sphere. (We use y here to avoid overlap with the 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 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 XT without a bar be the corresponding measured quantity. The toy chemistry diagnostics are then c2(t), the l2 norm of XT-X¯T at time t normalized by 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.

Figure 16Toy chemistry diagnostic values as a function of time for 10 cycles of the non-divergent flow. Time is on the x axis and measured in cycles. Diagnostic values c2 (solid lines) and c (dashed lines) are on the y axis. Markers as listed in the bottom legend are placed at the start of each cycle to differentiate the curves.


Figure 17Images of the monatomic tracer at the end of the first cycle. Text at the lower left of each image states the configuration. Text at the upper right reports global extremal values.


Figure 18Same as Fig. 17 but, now, the images are of (XT-X¯T)/X¯T.


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 npt=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 c2(t); dashed, c(t). Markers are placed on the curves at the start of each cycle to differentiate the curves. In each case, nf=npt. We choose this value of nf 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 npt=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, c210-2 and c>10-1 at the end of one cycle, compared with approximately 10−7 and 10−5, respectively, for npt=4 at the end of 10 cycles. For npt>4, the growth in error is very small, with c2<10-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 npt=nf=4 with CAAS–CAAS (left) and npt=nf=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 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 (XT-X¯T)/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 npt=8 at the end of the first cycle in Fig. 16.

5 Conclusions

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 npt and nf 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 Jacob2017). 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 (Kaas2008). 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.

Appendix A: Islet bases

The GLL basis functions yield an unstable ISL method for np≥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 [-1,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 xGnp(i), i{0,,np-1}, be the reference coordinates of the GLL nodes. Let region r{0,,np-2} be the segment [xGnp(r),xGnp(r+1)]. A 1D basis function is a continuous piecewise polynomial over [-1,1]; in each region, it is a polynomial. Each region r has an associated ordered list of support nodes, denoted Irnp. Each list Irnp is a subset of {0,,np-1}. Let |Irnp| be the number of elements in Irnp, let npsubminminr|Irnp|, and let Irnp(j) be the jth element.

In region r and thus for x[xGnp(r),xGnp(r+1)], basis function k is 0 if kIrnp and otherwise has the value given by the Lagrange polynomial


We construct a basis for each np{4,,13}. Each basis is described by a set of support node lists, one Irnp per region. For np=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 Irnp is given by an offset, which is the first index in the list and its size, |Irnp|. For example, np=7, offset 2, and size 4 correspond to support nodes Ir7={2,3,4,5}. In addition, the basis is symmetric, meaning basis function ϕinp(x)=ϕnp-1-inp(-x). Thus, first, support nodes are specified for regions 0 through np/2-1, and the support nodes for the remaining regions are determined by symmetry. Second, if np is even, then the middle region, r=np/2-1, has support nodes Irnp 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 npsubmin-1 if there is no property preservation step and the initial condition is Cnpsubmin-2.

If npsubmin=np, 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 np{5,,13}. In addition, a special basis is constructed for np=4.

Table A1Islet 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.

Download Print Version | Download XLSX

A1np=4 Islet basis

For np=4, we boost the accuracy of the basis by combining an offset nodal basis and the GLL basis. The middle region has support nodes I1={0,1,2,3}, i.e., all four nodes. Region 0 convexly combines basis function k given by I04={0,1,2} with GLL basis function k over the region x[xG4(0),xG4(1)], 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 α(-1)=1, only the GLL basis function k is used; at the right side, where α(xG4(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

The element-local conservation expressed in Eq. (15) holds for the following reasons. First, for any Islet basis for which (npt)subminnpv, the continuum field f is the same on each grid, since the tracer grid can exactly represent polynomials of degree npv-1. Second, for npv even, such as the standard npv=4, if (npt)submin=npv-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(x)i=0daixi, -11f(x)dx = -11(f(-x)+f(x))/2dx, and g(x)(f(-x)+f(x))/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 npv=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 12×100=1200 d); dash-dotted, natural basis for 1 cycle. For each npt, the l2 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 npt=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.

Figure A1Stability of the Islet method with the Islet GLL bases compared with the instability of the method with the natural GLL bases. 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 10l2 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 npt, 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 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

(B1) d d t Ω ( t ) ρ d x = Ω ( t ) ρ t + ( ρ u ) d x .

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

(B2) ρ t + ( ρ u ) = 0 .

In addition, consider a function ϕ(x,t) that satisfies the advection equation,

(B3) D ϕ D t = 0 .

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


which, in turn, implies

(B4) Ω ( t 1 ) ρ ϕ d x = Ω ( t 0 ) ρ ϕ d x

for all times t0, t1.

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 Ω(t1) the arrival Eulerian element. Then,

(B5) Ω ( t ) f ( x , t ) d x = - 1 1 - 1 1 f ( x ( r , t ) , t ) J ( r , t ) d r .

Here, x(r,t)X(t;x(r,t1),t1) 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 np the number of 1D nodes. We apply tensor-product GLL quadrature with this value of np to the right hand side of Eq. (B5). Then,


where i iterates over nodes in the element, ri is the ith GLL node, and wi is the ith GLL weight. Let xi(t)x(ri,t).

Third, we combine the quadrature approximation with a specific form of f, f=ρϕ, where ϕ satisfies Eq. (B3). Let ϕ(x,t1) be an element, GLL nodal, tensor-product np-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 ik. ϕ need not be a natural GLL basis function; for example, it can be an Islet GLL basis function. Because Dϕ/Dt=0, ϕ(xi(t),t)=ϕ(xi(t1),t1) 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



(B6) ρ ( x k ( t 1 ) , t 1 ) = J ( r k , t 0 ) J ( r k , t 1 ) ρ ( x k ( t 0 ) , t 0 ) .

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 np=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 np-basis. The discretization of the continuity equation for ρ then has order of accuracy 2.

Code and data availability

Code and scripts for the algorithms and figures presented in this paper are available at (Bradley2021a) and, alternatively, (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 (Bradley2021b).

Author contributions

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.

Competing interests

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.

Financial support

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.

Review statement

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],, 2021a. a

Bradley, A. M.: Methods data for Islet paper, Zenodo [data set],, 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],, 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,, 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],, 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,, 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,, 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,, 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],, 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,, 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,, 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,, 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,, 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,, 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,, 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,, 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,, 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,, 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

Short summary
Tracer transport in atmosphere models can be computationally expensive. We describe a flexible and efficient interpolation semi-Lagrangian method, the Islet method. It permits using up to three grids that share an element grid: a dynamics grid for computing quantities such as the wind velocity; a physics parameterizations grid; and a tracer grid. The Islet method performs well on a number of verification problems and achieves high performance in the E3SM Atmosphere Model version 2.