## Abstract

Atmospheric modeling with element-based high-order Galerkin methods presents a unique challenge to the conventional physics–dynamics coupling paradigm, due to the highly irregular distribution of nodes within an element and the distinct numerical characteristics of the Galerkin method. The conventional coupling procedure is to evaluate the physical parameterizations (*physics*) on the dynamical core grid. Evaluating the physics at the nodal points exacerbates numerical noise from the Galerkin method, enabling and amplifying local extrema at element boundaries. Grid imprinting may be substantially reduced through the introduction of an entirely separate, approximately isotropic finite-volume grid for evaluating the physics forcing. Integration of the spectral basis over the control volumes provides an area-average state to the physics, which is more representative of the state in the vicinity of the nodal points rather than the nodal point itself and is more consistent with the notion of a “large-scale state” required by conventional physics packages. This study documents the implementation of a quasi-equal-area physics grid into NCAR’s Community Atmosphere Model Spectral Element and is shown to be effective at mitigating grid imprinting in the solution. The physics grid is also appropriate for coupling to other components within the Community Earth System Model, since the coupler requires component fluxes to be defined on a finite-volume grid, and one can be certain that the fluxes on the physics grid are, indeed, volume averaged.

## 1. Introduction

An increasing number of numerical methods publications in the atmospheric science literature concern transport, shallow water, and three-dimensional models employing element-based high-order Galerkin discretizations such as finite-element and discontinuous Galerkin methods [for an introduction to these methods, see, e.g., Durran (2010); Nair et al. (2011); Ullrich (2014)]. Some global models based on Galerkin methods have reached a level of maturity for which they are being considered for next-generation climate and weather models due to their inherent conservation properties, high-order accuracy (for smooth problems), high parallel efficiency, high processor efficiency, and geometric flexibility facilitating mesh-refinement applications. NCAR’s Community Atmosphere Model (CAM; Neale et al. 2012) offers a dynamical core based on continuous Galerkin finite elements (Taylor and Fournier 2010), CAM Spectral Element (CAM-SE; Taylor et al. 2008; Dennis et al. 2012; Lauritzen et al. 2018). CAM-SE is, in particular, being used for high-resolution climate modeling (e.g., Small et al. 2014; Reed et al. 2015; Bacmeister et al. 2018) and static mesh-refinement applications (e.g., Fournier et al. 2004; Zarzycki et al. 2014a,b; Guba et al. 2014b; Rhoades et al. 2016). Other examples of models based on high-order Galerkin methods that are being considered for “operational” weather–climate applications are Giraldo and Restelli (2008), Nair et al. (2009), Brdar et al. (2013), and the Energy Exascale Earth System Model (https://e3sm.org/).

Assumptions inherent to the physical parameterizations (also referred to as *physics*) require the state passed by the dynamical core represent a “large-scale state,” for example, in quasi-equilibrium-type convection schemes (Arakawa and Schubert 1974; Plant and Craig 2008). In finite-volume methods (e.g., Lin 2004), one may think of the dynamical core state as the average state of the atmosphere over a control volume and, for resolutions typical of climate simulations, is entirely consistent with the notion of a large-scale state. For finite-difference methods (e.g., Suarez et al. 1983), the point value is thought of as representative of the atmospheric state in the vicinity of the point value, and one can usually associate a volume with the grid point. Hence, the physics grid (the grid on which the state of the atmosphere is evaluated and passed to physics) and the dynamics grid (the grid the dynamical core uses) coincide. Having the physics and dynamics grids coincide is obviously convenient since no interpolation is needed (which could disrupt conservation properties), and the number of degrees of freedom on both grids is exactly the same.

For the regular latitude–longitude, cubed-sphere, and icosahedral grids, the distance between the grid points is gradually varying for finite-volume/finite-difference discretizations. Examples of models that use these grids are CAM finite volume (CAM-FV; latitude–longitude grid; Lin 2004), FV cubed (FV3; cubed-sphere grid; Putman and Lin 2007), and the Icosahedral Nonhydrostatic model (ICON; icosahedral grid; Wan et al. 2013). For high-order element-based Galerkin methods, the dynamical core grid is defined by the quadrature points. In CAM-SE, these are the Gauss–Lobatto–Legendre (GLL) quadrature nodes. A unique aspect of the high-order quadrature rules is that the nodes within an element are located at the roots of the basis set, which may be irregularly spaced. For example, Fig. 1 shows GLL points on an individual element of a cubed-sphere grid for degree-3 ( quadrature points) and degree-7 ( quadrature points) Lagrange polynomial basis used in CAM-SE. The higher the order of the quadrature rule, the greater variance in distance between GLL quadrature points within an element. GLL quadrature points cluster near the edges and, in particular, the corners of the elements.

The resolved scales of motion are not determined by the distance between quadrature nodes, but rather by the degree of the polynomial basis in each element. The nodes may be viewed as irregularly spaced samples of an underlying spectrally truncated state. From this perspective, one might expect the nodal solutions to be independent of location within an element. While the interior quadrature nodes are in CAM-SE (i.e., the basis representation is infinitely smooth, and all derivatives are continuous), the smoothness of boundary nodes is constrained by the need to patch neighboring solutions together to form the global basis set, an operation known as the direct stiffness summation (DSS; Maday and Patera 1987; Canuto et al. 2007). The DSS operation is attractive because it allows for high-order accuracy with minimal communication between elements, but degrades the solution to at element boundaries (i.e., all derivatives are discontinuous). Through evaluating the physics at the nodal points, strong gridscale forcing or oscillatory behavior near an element boundary may exacerbate the discontinuity, and our initial expectation that the nodal solutions are independent of within-element location is unlikely for nonsmooth problems (e.g., the presence of rough topography or moist physics gridscale forcing).

The purpose of this paper is to document the implementation of an entirely separate, quasi-equal-area finite-volume physics grid into CAM-SE. The use of a separate physics grid is not entirely unheard of; prior studies have utilized the infrastructure developed for global spectral transform methods to experiment with different physics grids (Williamson 1999; Wedi 2014). In our framework, the dynamical core state is integrated over control volumes to provide a volume-averaged state to the physics, thereby minimizing the influence of any one particular nodal value on the physics forcing. Section 2 provides a thorough explanation of how grid imprinting manifests in high-order Galerkin methods for nonsmooth problems. The implementation of the physics grid configuration into CAM-SE is presented in section 3. Results from a hierarchy of idealized model configurations are presented in section 4, illustrating the physics grid is effective at mitigating undesirable grid imprinting in the solution. Section 5 contains a discussion of results and concluding remarks.

## 2. The quadrature node problem

Figure 2 is a schematic illustrating in one dimension how grid imprinting is enabled by the physics, when the dynamical core is built using high-order Galerkin methods. The schematic depicts a time step, starting from smooth initial conditions (Fig. 2a) and subsequently advancing the dynamics one Runge–Kutta time step (Fig. 2b). Since the boundary nodes of adjacent elements overlap one another, there are now two solutions for each boundary node. The DSS operator, effectively a numerical flux applied to the element boundaries such that overlapping nodal values agree, is applied (Fig. 2c), rendering the solutions at element boundaries less smooth than neighboring interior nodes. An element boundary discontinuity may be exacerbated if, for example, the physics updates the state at an element boundary (Figs. 2d,e), resulting in characteristically tighter gradients on the boundary nodes, compared to if the physics forcing were applied to an interior node (Figs. 2g,h).

To test the degree to which nodal solutions depend on within-element position, an aquaplanet simulation (Neale and Hoskins 2000; Medeiros et al. 2016), which consists of an ocean-covered planet in perpetual equinox, with fixed, zonally symmetric sea surface temperatures idealized after the present-day climatology, is carried out using CAM-SE, using CAM, version 4 physics (CAM4; Neale et al. 2010) and run for 1 year. The nominal low-resolution grid is used, pertaining to an average equatorial grid spacing of 111.2 km. The probability density distribution of the upward vertical pressure velocity *ω*, conditionally sampled based on three categories—“interior,” “edge,” and “corner” nodes—is provided in Fig. 3a. The motivation for assessing noise in the *ω* field comes from its connection with the atmosphere’s divergent modes, as follows from the continuity equation in pressure coordinates. These modes are in turn sensitive to the within-element inhomogeneity of the pressure gradient that emerges from high-order Galerkin methods. There is an apparent dependence on nodal location, with interior nodes being characteristically sluggish and corner and edge nodes having systematically larger-magnitude vertical motion. This behavior is consistent with the smoothness properties of the different nodal locations, with discontinuous pressure gradients resulting in greater vertical motion at edge and corner nodes. The main division of solutions shown in Fig. 3a is primarily from whether a node is or is not situated on an element boundary and is a nuanced signature of high-order element-based Galerkin methods for nonsmooth problems.

If the conventional physics–dynamics coupling paradigm is applied to CAM-SE, then the physics are to be evaluated at the GLL nodes, and a volume associated with the quadrature point should be defined. One approach to construct this grid is to decompose each spectral element into subcells and then take the dual grid of this subcell grid. For cubed-sphere meshes, this dual grid will have a control volume associated with each quadrature point. These control volumes will be triangles for the cube corner quadrature points and quadrilaterals for all remaining quadrature points. Newton iteration can then be used to adjust the corners of these control volumes so that their spherical areas exactly match the Gaussian weight multiplied by the metric term (these weights are used for integrating the basis functions over the elements and can therefore, in this context, be interpreted as areas). For cubed-sphere meshes, the Newton iteration can be replaced by a direct method if some of the quadrilaterals are replaced by pentagons, giving additional flexibility in matching the spherical area to the quadrature weights. Such a dual grid is shown in Fig. 4. This grid is used in the NCAR Community Earth System Model (CESM) coupler for passing states among ocean, atmosphere, and land components since the current remapping method is finite-volume based and therefore requires control volumes [it is noted that methods exist that do not require control volumes for conservative interpolation; e.g., Ullrich and Taylor (2015)]. Hence, the components “see” an irregular atmospheric grid. Similarly, the parameterizations in the atmosphere see a state that is anisotropically sampled in space [see Figs. 1 and 5 in Kim et al. (2008)].

The quadrature grid in element-based Galerkin methods is defined to perform mathematical operations on the basis functions (e.g., computing gradients and integrals), rather than evaluating the state variables for physics–dynamics coupling. One may argue that it would be more consistent to integrate the basis functions over quasi-equal-area control volumes within each element and pass those control volume-average values to physics rather than irregularly spaced quadrature point values. In this case, when integrating basis functions over control volumes, a gridcell average value is more representative of the values near the extrema at the element boundary than the quadrature point value. The relationship among the nodal values, the basis functions, and the proposed control volumes is illustrated schematically in one dimension in Figs. 2f and 2i.

## 3. Methods

Here, we focus on CAM-SE; however, in principle, the methods apply to any element-based high-order Galerkin model. The physics grid in CAM-SE is defined by subdividing each element using equiangular gnomonic coordinate lines to define the sides of the physics grid control volumes (see the appendix for details). Note that the element boundaries are defined by equiangular gnomonic grid lines. The notation refers to the configuration where the elements are divided into equiangular physics grid cells (see Fig. 5), resulting in a quasi-equal spherical area grid resembling the cubed sphere. Defining the physics grid by subdividing elements makes it possible to use the same element infrastructure as already used in CAM-SE, thereby facilitating its implementation. Here, we make use of the and grids that use the GLL quadrature point physics grid (physics and dynamics grid coincide) and the same () resolution quasi-equal-area physics grids, respectively. In all configurations, we use degree-3 Lagrange basis () and elements on each cubed-sphere panel.

A consequence of separating physics and dynamics grids is that the atmospheric state must be mapped to the physics grid, and the physics tendencies must be mapped back to the dynamics grid. This is discussed in separate sections below. When separating physics and dynamics grids, it is advantageous to use a vertical coordinate that is static during physics–dynamics coupling. This was one motivation to switch to a dry-mass vertical coordinate in CAM-SE (Lauritzen et al. 2018); since dry mass remains constant throughout physics, the dry-mass vertical coordinate remains fixed during physics–dynamics coupling. The dry-mass coordinate subsequently evolves as floating Lagrangian layers by the dynamics (Lin 2004), periodically mapped back to a reference hybrid-sigma-pressure coordinate after Simmons and Burridge (1981). All variables mapped between grids are collocated, layer-mean values (Lauritzen et al. 2018).

### a. Mapping state from dynamics grid (GLL) to physics grid (pg)

The dynamics state is defined on the GLL grid in terms of temperature , zonal wind component , meridional wind component , and dry pressure-level thickness . In the mapping of the atmospheric state to the physics grid, it is important that the following properties are met:

Conservation of scalar quantities such as mass and dry thermal energy.

For tracers, shape preservation (monotonicity): the mapping method must not introduce new extrema in the interpolated field, in particular, negatives.

Consistency: the mapping preserves a constant.

Linear correlation preservation: if field

*A*is a linear function of*B*, this relationship is still preserved [see, e.g., Eq. (5) in Lauritzen and Thuburn (2012)]

Other properties that may be important, but not pursued here, include total energy conservation and axial angular momentum conservation. Total energy is a quadratic quantity that is inherently difficult to conserve unless one maps total energy, requiring one to diagnose either temperature or momentum components. For example, enforcing total energy conservation locally—using, for example, Lin’s (2004) method where total energy and velocity components are remapped and temperature is a derived variable—has proven problematic (C. Chen 2015, personal communication). Similarly, conservation of axial angular momentum is problematic. Conservation of angular momentum requires one to interpolate the zonal and meridional components of momentum, which creates large errors near the poles. To avoid the pole problem, we interpolate contravariant components of the momentum vector, which violates axial angular momentum conservation.

We argue that the most consistent method for mapping scalar-state variables from the GLL grid to the physics grid is to integrate the Lagrange basis function representation (used by the SE dynamical core) over the physics grid control volumes [i.e., integrate the basis function representation of and over the physics grid control volume; see, e.g., Lauritzen et al. (2017); Ullrich and Taylor (2015)]:

where is the physics grid area. The integrals are numerically computed using the GLL quadrature rule on each physics grid element, which exactly (to machine precision) integrates the basis functions over the control volumes (Lauritzen et al. 2017). Thermal energy and dry air mass is conserved, and the mapping is consistent. For the wind, which is a vector, the zonal and meridional wind components are mapped by transforming to contravariant wind components, evaluating the basis function representation thereof at the equiangular center of the physics grid control volumes and then transformed back to latitude–longitude coordinate system winds. All of the operations are local to the element and do not require communication between elements.

The mapping of tracers is more problematic since the SE basis function representation is oscillatory, although the shape-preserving filter guarantees shape preservation at the GLL nodes (Guba et al. 2014a). To avoid this issue, we use the Conservative Semi-Lagrangian Multitracer transport scheme version of CAM-SE (CAM-SE-CSLAM; Lauritzen et al. 2017), where tracers are advected on the physics grid using the inherently mass- and linear-correlation-preserving CSLAM algorithm. Note that in CAM-SE-CSLAM, the dry mass internally predicted by CSLAM is, by design, equal to integrated over the CSLAM/physics grid control volume (Lauritzen et al. 2017). Since the tracer grid and physics grids are collocated, and , the mass conservation, correlation preservation, consistency, and shape-preservation constraints are inherently fulfilled.

### b. Mapping tendencies from physics grid (pg) to dynamics grid (GLL)

The physics tendencies are computed on the finite-volume physics grid and are denoted , , , and . Note that dry air mass is not modified by physics, and hence, there is no tendency for dry mass . Also, it is important to map tendencies and not state from the physics grid to GLL grid; otherwise, one will get spurious tendencies from mapping errors when the actual physics tendency is zero (unless a reversible map is used).

It is important that this process abides by the following:

For tracers, mass tendency is conserved.

For tracers, in each tracer grid cell, the mass tendency from physics must not exceed tracer mass available in tracer grid cell (it is assumed that the physics tendency will not drive tracer mixing ratio negative on the physics grid).

Linear correlation preservation.

Consistency (i.e., the mapping preserves a constant tendency).

Other properties that may be important, but not pursued here, include total energy conservation (including components of total energy) and axial angular momentum conservation. Scalar variables are mapped from the physics grid to GLL grid using a tensor-product Lagrange interpolation in two dimensions (i.e., we assume that the pressure variations in the vertical are small). The local coordinates on a cubed sphere are discontinuous at the element edges, so the interpolation requires special attention at the cube corners and edges. The details are provided in the appendix. Lagrange interpolation preserves a constant (including zero) and linear correlations. Tracer and physics grids are collocated, so tracer mass, tracer shape, and tracer correlations are trivially preserved on the tracer grid, and the inconsistency in point 2 above will not appear.

Mapping from to GLL grids while conserving mass was found to be difficult without excessive grid imprinting at element edges. Mass conservation (using conventional finite-volume methods) requires a control volume to be defined around the GLL points [see Fig. 4 in this paper or Fig. 8b in Ullrich et al. (2016)]. These volumes are artificial and not consistent with the SE method. Integrating the CSLAM reconstruction of water tracers of such artificial control volumes led to GLL node grid imprinting in the mapping and will not preserve a constant mixing ratio since the mapping of to GLL will not yield the GLL node value for dry pressure-level thickness (i.e., the maps are not reversible). A reversible map requires that the number of degrees of freedom on the source mesh ( has 9 degrees of freedom) equals the number of degrees of freedom on the target mesh ( grid has 16 degrees of freedom). This condition is violated by construction for individual elements.

It was also found important to use an interpolator that is smooth across element boundaries. Using an algorithm that only uses information from an element of control volumes will (at best) be at the element boundaries and therefore lead to boundary node grid imprinting. A stencil that extends beyond one element is therefore necessary. After much experimentation, the best results in terms of grid imprinting were obtained with tensor-cubic interpolation (see the appendix for details) and by using the CAM-SE-CSLAM configuration (which requires the same boundary exchange/communication as used in CSLAM).

### c. Time splitting and physics–dynamics coupling

The physics and dynamics are integrated in time using a sequential-update approach (e.g., Williamson 2002). The dynamical core is subcycled over the (usually) longer physics time step (e.g., the vertical remapping time step is cycled times, totaling to ). In CAM-SE, a fraction of the physics forcing (e.g., ) is applied at the beginning of each vertical remap subcycle, such that the full forcing () is realized over the course of a physics time step. This approach of dribbling the tendencies over subintervals has the advantage of reducing gravity wave noise (Thatcher and Jablonowski 2016), but may disrupt tracer mass conservation (Zhang et al. 2018). In CAM-SE-CSLAM, all but the tracer mass quantities are dribbled, with tracer mass receiving the full physics update (e.g., ), applied only at the beginning of the first remap subcycle, and thereby conserving tracer mass. This is the configuration, described in detail in section 3.6.3 of Lauritzen et al. (2018).

In the SE integration of the equations of motion on the GLL grid, the water species are needed in the computation of the pressure gradient force and generalized expressions for heat capacity at constant pressure and so forth. Hence, the mixing ratios for water vapor and dynamically/thermodynamically active condensates (e.g., cloud liquid and cloud ice) are needed on the GLL grid. We have chosen to advect the water species on the GLL grid using the SE method, as well as on the physics grid using CSLAM. Every time physics updates the water species on the CSLAM grid, a forcing term (equal to the difference between updated CSLAM water variables and the SE values) is applied to the GLL water variables using dribbling so that the CSLAM solution and SE solution for water species are tightly coupled.

## 4. Results

A hierarchy of idealized model configurations is presented in order to elucidate the differences between CAM-SE and CAM-SE-CSLAM (available from the CESM2.1 release; UCAR/NCAR 2018). Here, the configurations are presented in order of increasing complexity, each with a pair of approximately 1° simulations, pertaining to the (CAM-SE) and (CAM-SE-CSLAM) grids, and a = 1800 s.

### a. Moist baroclinic wave

The moist baroclinic wave test case was developed as part of the CESM simpler models project (Polvani et al. 2017) and included in the release of CESM2. It is effectively the dry test case of Ullrich et al. (2014), but initialized with moisture and coupled to the Kessler moist physics routine (Kessler 1969). For more details on this test case (which was part of the 2016 Dynamical Core Model Intercomparison Project; Ullrich et al. 2017), see section 4.1 of Lauritzen et al. (2018). A measure of the uncertainty in the reference solution, the difference norm between two high-resolution solutions using different dynamical cores, was also presented by Lauritzen et al. (2018) and provided again here in Fig. 6. The norm between CAM-SE and CAM-SE-CSLAM lies below the uncertainty of the reference solution, indicating their differences are insignificant.

The flow field of the baroclinic wave test is used to drive the terminator “toy” chemistry test of Lauritzen et al. (2015b, 2017). The terminator test is used to assess linear-correlation preservation using two reactive species advected across the terminator line. The model is initialized with species for which their weighted sum is a constant (constant surface pressure and constant mixing ratio; ), such that if tracer correlations are preserved, then the column-integrated weighted sum of the species should not vary in time. Figure 7 provides a snapshot of the vertically integrated weighted sum of species at day 15. In CAM-SE, the tracer correlations are not preserved at day 15, and the field is populated by overshoots and undershoots. In contrast, by day 15, CAM-SE-CSLAM still conserves tracer correlations to within machine precision, consistent with the previous results of this test case initialized with a dry baroclinic wave (Lauritzen et al. 2017).

### b. Aquaplanets

Two-year-long aquaplanet simulations are performed using CAM-SE and CAM-SE-CSLAM, using the CAM4 physics package (Neale et al. 2010), as discussed in section 2. Away from the grid scale, the mean states in the two models are very similar. Figure 8 shows the zonal-mean climatological precipitation rates in CAM-SE and CAM-SE-CSLAM. Considering how sensitive this aquaplanet configuration is to design choices in CAM-SE (Lauritzen et al. 2018), it is somewhat unexpected that the zonal means look so similar to one another.

A plot similar to Fig. 3a is constructed for the CAM-SE-CSLAM simulation, a probability density distribution of upward *ω* conditionally sampled based on location within the element. Like Fig. 3a, Fig. 3b divides up the control volumes by corner, edge, and interior cells. Through the use of the quasi-equal-area physics grid, the dynamical core state appears more-or-less independent of location within the element, a marked improvement over CAM-SE. Since the state is approximately independent of in-element location, it follows that the physics forcing, which is evaluated from the dynamical core state, may be expected to also show an improvement in grid imprinting.

The low level, mean, and variance of the temperature tendencies from the physics on the GLL grid in the two simulations are shown in Fig. 9. The mean states in the two models resemble one another, consistent with the zonal-mean precipitation rates (Fig. 8). The mean physics tendencies contain modest grid imprinting in CAM-SE (barely visible near the storm track regions), while in the variance field, grid imprinting is both ubiquitous and unmistakable. The variance is larger on boundary nodes, manifesting as a “stitching” pattern resembling the cubed-sphere grid. In CAM-SE-CSLAM, the grid imprinting is all but eliminated, based on the mean and variance of the physics tendencies (Fig. 9), consistent with our expectation.

The global mean and variance of the low-level physics tendencies are marginally lower in CAM-SE-CSLAM, compared with CAM-SE on the GLL grid (by about 1% and 6% for the mean and variance, respectively; Fig. 9). While these differences may be small and potentially insignificant, they are consistent with the state on the GLL grid in the two simulations. Through recreating Fig. 3a, but using the *ω* field on the GLL grid in the CAM-SE-CSLAM run, the frequency of large magnitude *ω* values (less than −1.0 Pa s^{−1}) associated with interior, corner, and edge nodes is slightly lower (not shown). This suggests that the lower-magnitude physics forcing in CAM-SE-CSLAM impacts the state on the GLL grid, albeit modestly. Therefore, the lower frequency of large magnitude *ω* in CAM-SE-CSLAM (Fig. 3) may not be solely due to the smoothing effect of integrating the basis functions over control volumes, but also to the lower-magnitude physics tendencies feeding back onto the dynamical state.

As stated in section 3, the mapping of the state to the physics grid and the reverse interpolation of physics tendencies to the GLL grid is not total energy conserving. CAM has a global energy fixer (Williamson et al. 2015), which can be used to estimate the errors associated with the mapping algorithms. To do so, it is presumed that there are no compensating mapping errors in going to and from the physics and dynamics grids, and CAM-SE-CSLAM and CAM-SE have the same energy dissipation rates. Under these assumptions, the spurious globally integrated total energy errors due to the mapping algorithm are estimated to be approximately 0.0025 W m^{−2} in the aquaplanet simulations. In comparison, the dynamical core total energy dissipation is on the order of 0.1 W m^{−2} (Lauritzen et al. 2018).

### c. Held–Suarez with topography

Grid imprinting associated with the flow around obstacles is more problematic than that encountered on the aquaplanets. To diagnose grid imprinting due to topographic flow, an idealized Held–Suarez configuration (Held and Suarez 1994) is outfitted with real-world topography after Fox-Rabinovitz et al. (2000) and Baer et al. (2006) and run for 2 years. Figure 10 shows the mean *ω* at two different vertical levels in the midtroposphere. The data are presented as a raster plot on their respective unstructured grids in order to delineate whether a particular value is associated with an interior, edge, or element boundary node.

At higher latitudes (e.g., the southern Andes), the flow is smooth, conforming reasonably to the underlying topography. At lower latitudes, over the Andes (between the equator and 20°S) or the Himalayas (from 20° to 30°N), there is a clear preference for extrema to occur at the element boundaries (Fig. 10). The vertical structure of *ω* in regions of strong grid imprinting indicates full-troposphere upward/downward motion (not shown). Grid imprinting is therefore more common in regions of weak stratification, such as occurs in the deep tropics, with forced upslope flow facilitating the release of gravitational instability. Resolved updrafts/downdrafts often align with the element boundaries due to its systematically tighter pressure gradients.

Through the use of the quasi-equal-area physics grid, grid imprinting due to topographic flow is reduced (Fig. 10). The native topography lives on the physics grid, and the topography is mapped to the nodal points at run time in CAM-SE-CSLAM. Mapping topography to the quadrature nodes ensures that no new extrema will be introduced to the boundary nodes, where the solution is least smooth. This effect cannot be very large, since grid noise over topography is similar in CAM-SE and CAM-SE-CSLAM on the GLL grid (not shown). From the perspective of the physics grid, CAM-SE-CSLAM clearly mitigates the influence of grid-induced extrema on the state. This can be seen by comparing Figs. 10a and 10b and their differences (Fig. 10c), which shows that the largest differences coincide with the element boundaries. The reduction in grid imprinting in this modified Held–Suarez configuration appears to be almost entirely a result of the smoothing effect of integrating the basis functions over the control volumes of the physics grid.

### d. AMIP-type simulations

A pair of 20-yr-long AMIP-type simulations is performed using CAM, version 6 physics package (CAM6) and using perpetual year 2000 SST boundary conditions ( compset in CESM2.0; UCAR/NCAR 2018). Figure 11 shows the climatological precipitation fields in CAM-SE (left) and CAM-SE-CSLAM (middle), and over the same mountainous regions as in Fig. 10. The plots have some similar features to the *ω* field in the Held–Suarez runs: the greater variances at lower latitudes and on the windward side of the mountains are broadly similar. CAM-SE-CSLAM has a lower spatial variance (e.g., the lack of extrema over the Andes at about 15°S, compared to CAM-SE, and the gridscale precipitation peak over the Himalayas at about 30°N). The difference plot (Fig. 11; right panel) is more broadly populated by blue, purple, and white contours, indicating that CAM-SE has, in general, larger-magnitude precipitation rates over high topography. The difference plots also highlight a couple of zonally aligned strips of anomalous precipitation, in particular, near the foot of the Himalayas in CAM-SE. These bands are in the same location as the bands of precipitation identified in CAM-SE by Lauritzen et al. (2015a, their Fig. 7), but using CAM, version 5 physics, which they argue are spurious in nature.

To assist in identifying whether a particular precipitation pattern is spurious, an simulation is carried out using the finite-volume dynamical core that uses a regular latitude–longitude 0.9° × 1.25° grid (CAM-FV; grid; Neale et al. 2012). CAM-FV is the default low-resolution model in CESM2.0, and with its smoothly varying grid, it does not suffer from the quadrature node problem (section 2). Figure 12 shows the global precipitation fields in CAM-SE, CAM-SE-CSLAM, and CAM-FV, compared to an observational dataset, the Global Precipitation Climatology Project (GPCP; 1979–2003) gridded dataset (Huffman et al. 2001). The magnitudes of the precipitation rates in all three models are higher than the GPCP dataset, primarily over land in the tropics (note the lack of red contours in the GPCP dataset), which should be interpreted cautiously due to widely accepted issues in constructing a reliable, gridded, global precipitation dataset. At lower latitudes, CAM-FV has lower spatial variance and overall lower magnitudes, compared with CAM-SE. The GPCP dataset indicates that perhaps the precipitation rates in low-latitude mountainous regions in CAM-FV and CAM-SE are larger than in reality. Following suit, the reduction in magnitude and spatial variance in precipitation in these regions in CAM-SE-CSLAM may be interpreted as an improvement over CAM-SE.

## 5. Conclusions

Element-based high-order Galerkin methods possess many of the attractive qualities recommended for next-generation global atmospheric models. Among these, high-order accuracy is achieved with minimal communication between elements, allowing for near-perfect scaling on massively parallel systems. Element communication amounts to a numerical flux applied to the element boundaries, reconciling overlapping solutions of adjacent elements but degrading the smoothness of the boundary nodes in the process (to ). For nonsmooth problems, gradients are systematically tighter at the element boundaries, and local extrema often characterize the boundary nodes. This behavior is illustrated using NCAR’s Community Atmosphere Model Spectral Element (CAM-SE) (CAM-SE) dynamics in an aquaplanet configuration, in a Held–Suarez configuration with real-world topography, and in an AMIP-type configuration.

The authors argue that the conventional physics–dynamics coupling paradigm, in which the physical parameterizations are evaluated on the dynamical core grid, exacerbates grid imprinting. A separate physics grid is proposed and implemented in CAM-SE, referred to as CAM-SE-CSLAM, through dividing the elements into quasi-equal areas with equivalent degrees of freedom. The state is mapped to the physics grid with high-order accuracy through integrating CAM-SE’s Lagrange basis functions over the control volumes. Control volumes near element boundaries now represent a state in the vicinity of the extrema produced through the boundary exchange, as opposed to the nodal value itself. These control volumes are also compatible with a large-scale state, as required by the physical parameterizations. The physical parameterizations are evaluated on the finite-volume grid, and the forcing terms are mapped back to the dynamical core grid using a cubic tensor-product Lagrange interpolation. In aquaplanet simulations, evaluating the parameterizations on the physics grid removes any obvious dependence of proximity to the element boundary, resulting in a more realistic state with negligible grid imprinting. The mapping algorithm does not conserve total energy, but it is estimated that these errors are one to two orders of magnitude less than the total energy dissipation from the dynamical core.

In CAM-SE-CSLAM, the physics grid replaces the default CAM-SE quadrature point-based coupler grid (Fig. 4) to compute fluxes between model components in the Community Earth System Model (CESM). The appeal here is twofold. Through integrating the Lagrange basis functions over control volumes, one can be certain that the fluxes computed from this grid are a volume-averaged flux. The same cannot be said for CAM-SE, where artificial control volumes (with sizes proportional quadrature weights) are constructed around nodal values and assumed to represent the volume-averaged state. The second advantage of the new coupler grid is that extrema occurring on boundary nodes may no longer influence other model components in simulations without rough topography. While grid imprinting is effectively eliminated in the aquaplanets, experiments with real-world topography (Held–Suarez and AMIP-type configurations) reduces, but does not entirely eliminate, imprinting from the mean state. The quasi-equal-area physics grid is nonetheless effective at mitigating numerical nuances associated with high-order element-based Galerkin methods for nonsmooth problems.

Future work will focus on the impact of using a coarser, physics grid configuration. The coarser physics grid may be more effective at reducing spurious noise over regions of rough topography, while potentially reducing the computational overhead. Any advantages of using a coarser-resolution physics grid will be weighed against any potential reduction in a model’s effective resolution.

## Acknowledgments

NCAR is sponsored by the National Science Foundation (NSF). We thank three anonymous reviewers for their helpful comments that significantly improved the clarity of the manuscript. Herrington thanks NCAR’s Computational and Information Systems Laboratory (CISL) and NCAR’s Climate and Global Dynamics division (CGD) for computational resources and technical support. Herrington, Reed, and Lauritzen are grateful to the NCAR Advanced Study Program graduate visitor program for funding Herrington’s 12-month visit. Goldhaber was partially supported by the Accelerated Climate Modeling for Energy (ACME) project and work package 12-015334: “Multiscale Methods for Accurate, Efficient, and the Accelerated Scale-Aware Models,” both funded through the U.S. Department of Energy Office of Biological and Environmental Research. Goldhaber is also partially supported through NSF Award AGS-1500187: “Development and Testing of a Global Quasi 3D multi Modeling Framework.” Model source code is officially released with CESM2.1 (UCAR/NCAR 2018).

### APPENDIX

#### High-Order Mapping to the GLL Grid

The mapping of the physics tendencies from the physics grid to the GLL grid is done with tensor-cubic Lagrange interpolation. The elements of the cubed sphere in SE are created from an equiangular gnomonic projection. Consider one element , where are central angle coordinates, and and are the minimum and maximum central angles in the *α*-coordinate direction, respectively, and similarly for *β*. Let and . The physics gridcell central angle centers are located at

where . The interpolation is performed in central-angle coordinates using tensor-product cubic interpolation. For elements located on a cubed-sphere edge or corner, the coordinate system for neighboring elements may be on a different panel. To take into account this coordinate change, the central-angle locations of physics gridcell centers located in other panels are transformed to the coordinate system of the panel in which the element in question is located [the transformations are given in, e.g., Nair et al. (2005)]. An illustration is given in Fig. A1 for an element located in the lower-left corner of a panel. The element in question is , where, for simplicity, we have transformed the element coordinates into normalized coordinates

[also used internally in the SE dynamical core; see, e.g., section 3.3 in Lauritzen et al. (2018)]. The GLL points are located at −1,, , and 1 in each coordinate direction. Near the edges/corners of an element cubic extrapolation are used if the centered stencil expands beyond the panel.

## REFERENCES

*Spectral Methods: Evolution to Complex Geometries and Applications to Fluid Dynamics*. 1st ed. Springer, 596 pp.

*Numerical Methods for Fluid Dynamics: With Applications to Geophysics*. 2nd ed. Texts in Applied Mathematics, Vol. 32, Springer, 465 pp.

*On the Distribution and Continuity of Water Substance in Atmospheric Circulations*.

*Meteor. Monogr.*, No. 10, Amer. Meteor. Soc., 76 pp.

*State of the Art Surveys on Computational Mechanics*, A. K. Noor, and J. T. Oden, Eds., ASME, 71–143.

*Numerical Techniques for Global Atmospheric Models*, P. H. Lauritzen et al., Eds., Lecture Notes in Computational Science and Engineering Series, Vol. 80, Springer, 251–311.

## Footnotes

© 2018 American Meteorological Society. For information regarding reuse of this content and general copyright information, consult the AMS Copyright Policy (www.ametsoc.org/PUBSReuseLicenses).