Impact of Grids and Dynamical Cores in CESM2.2 on the Surface Mass Balance of the Greenland Ice Sheet

Abstract

Six different configurations, a mixture of grids and atmospheric dynamical cores available in the Community Earth System Model, version 2.2 (CESM2.2), are evaluated for their skill in representing the climate of the Arctic and the surface mass balance of the Greenland Ice Sheet (GrIS). The finite-volume dynamical core uses structured, latitude-longitude grids, whereas the spectral-element dynamical core is built on unstructured meshes, permitting grid flexibility such as quasi-uniform grid spacing globally. The 1°–2° latitude-longitude and quasi-uniform unstructured grids systematically overestimate both accumulation and ablation over the GrIS. Of these 1°–2° grids, the latitude-longitude grids outperform the quasi-uniform unstructured grids because they have more degrees of freedom to represent the GrIS. Two Arctic-refined meshes, with 1/4° and 1/8° refinement over Greenland, were developed for the spectral-element dynamical core and are documented here as newly supported configurations in CESM2.2. The Arctic meshes substantially improve the simulated clouds and precipitation rates in the Arctic. Over Greenland, these meshes skillfully represent accumulation and ablation processes, leading to a more realistic GrIS surface mass balance. As CESM is in the process of transitioning away from conventional latitude-longitude grids, these new Arctic-refined meshes improve the representation of polar processes in CESM by recovering resolution lost in the transition to quasi-uniform grids, albeit at increased computational cost.

Introduction

General Circulation Models (GCMs) are powerful tools for understanding the meteorology and climate of the high latitudes, which are among the most sensitive regions on Earth to global and environmental change. GCMs differ vastly in their numerical treatment of polar regions because of the so-called pole problem (Williamson, 2007). The pole problem refers to numerical instability arising from the convergence of meridian lines into polar singularities on latitude-longitude grids (e.g., Figure 1a, hereafter referred to as lat-lon grids). Depending on the numerics, methods exist to suppress this instability, and lat-lon grids may be advantageous for polar processes by representing structures with finer resolution than elsewhere in the computational domain. With the recent trend toward quasi-uniform unstructured grids, any potential benefits of lat-lon grids in polar regions may be lost (hereafter, quasi-uniform refers to approximately isotropic grids, in contrast to lat-lon grids, which are highly anisotropic due to the polar singularity). In this study, we evaluate a number of grids and dynamical cores (hereafter referred to as dycores) available in the Community Earth System Model, version 2.2 (CESM2.2; Danabasoglu et al., 2020), including new variable-resolution grids (i.e., grids with enhanced resolution over a particular region), to understand their impacts on the simulated Arctic climate. We focus specifically on the climate and surface mass balance of the Greenland Ice Sheet.

Figure 1. Computational grids for the 1°–2° lat-lon and quasi-uniform unstructured grids in this study, using a northern hemisphere polar projection. Grids names after Table 1, (a) 2° lat-lon grid (f19) $(\mathtt{f}19)$, (b) 1° lat-lon grid (f09) $(\mathtt{f}09)$, (c) 1° quasi-uniform unstructured grid with reduced physics resolution (ne30pg2) $(\mathtt{n}\mathrm{e}30\mathrm{pg}2)$ and (d) 1° quasi-uniform unstructured grid (ne30pg3) $(\mathtt{n}\mathrm{e}30\mathrm{pg}3)$.

In the 1970s, the pole problem was largely defeated through the adoption of efficient spectral transform methods in GCMs (see Williamson, 2007, and references therein). These methods transform grid point fields into a global, isotropic representation in wave space, where linear operators (e.g., horizontal derivatives) in the (truncated) equation set can be solved exactly. While spectral transform methods are still used today, local numerical methods have become desirable for their ability to run efficiently on massively parallel systems. The pole problem has thus re-emerged in contemporary climate models that use lat-lon grids, and some combination of reduced grids (modified lat-lon grids in which cells around the polar singularity are elongated in the zonal direction) and polar filters are necessary to ameliorate this numerical instability (Jablonowski & Williamson, 2011). Polar filters subdue the growth of unstable computational modes by applying additional damping to the numerical solution over polar regions. This damping reduces the effective resolution in polar regions such that the resolved scales are approximately the same everywhere on the grid. We emphasize approximately, since it is conceivable that marginal increases in effective resolution occur over polar regions in lat-lon grids, despite polar filtering, since resolved waves can be represented with more grid points than at lower latitudes.

Dycores built on lat-lon grids have some advantages over dycores built on unstructured grids. Lat-lon coordinate lines are orthogonal, and aligned with zonally symmetric circulations that characterize many large-scale features of Earth's atmosphere. Lauritzen et al. (2010) has experimented with rotating lat-lon models such that their coordinate lines no longer align with an idealized, zonally balanced circulation. For the finite-volume lat-lon dycore considered in this paper (hereafter FV), numerical errors were shown to be largest when the polar singularity is rotated into the baroclinic zone (45°N latitude), generating spurious wave growth much earlier in the simulation than for other rotation angles. This illustrates the advantages of coordinate surfaces aligned with latitude bands, albeit an extreme example where the polar singularity and the polar filter are also contributing to the spurious wave growth. The unstructured grids all generate spurious baroclinic waves earlier in the simulations than the (unrotated) lat-lon models, although the unstructured model considered in this paper, the spectral-element dycore (hereafter SE), holds a balanced zonal flow without spurious wave growth appreciably longer than the rotated FV experiments (Lauritzen et al., 2010). And unlike FV, the SE dycore has the same error characteristics regardless of how the grid is rotated.

The polar filter in the FV model impedes efficiency at large processor (CPU) counts because it requires a spectral transform, which has a large communication overhead (Dennis et al., 2012; Suarez & Takacs, 1995). Unstructured grids support quasi-uniform grid spacing globally, and there is no pole problem (e.g., Figure 1c). This is in part why unstructured grids are becoming more common; their improved performance on massively parallel systems and lack of constraints on grid structure (Putman & Lin, 2007; Taylor et al., 1997; Wan et al., 2013). This increased grid flexibility allows for the adoption of variable-resolution grids (e.g., Figure 2; hereafter abbreviated as VR), sometimes referred to as regional grid refinement. In principle, grid refinement over polar regions can make up for any loss of resolution in transitioning away from lat-lon grids (e.g., Figure 2). However, local grid refinement comes at the cost of a smaller CFL-limited time step in the refined region; the CFL-condition—short for Courant–Friedrichs–Lewy condition—is a necessary condition for numerical stability when using discrete data in time and space.

Figure 2. Variable-resolution grids available in CESM2.2, (a) Arctic $\mathtt{A}\mathrm{rctic}$, (b) Arctic−GrIS $\mathtt{A}\mathrm{rctic}-\mathrm{GrIS}$ and (c) CONUS $\mathtt{C}\mathrm{ONUS}$. Note what is shown is the element grid; the computational grid has 3 × 3 independent grid points per element.

We emphasize that the pole problem is a distinctive feature of the dycore in atmospheric models. Polar filters do not directly interfere with the physical parameterizations, nor do they have any bearing on the surface models; for example, the land model can take full advantage of the greater number of grid cells in polar regions on lat-lon grids. This is particularly relevant for the surface mass balance (SMB; the integrated sum of precipitation and runoff) of the Greenland Ice Sheet, which relies on hydrological processes represented in the land model.

The SMB of the Greenland Ice Sheet (hereafter GrIS) is determined by processes occurring over a range of scales (Fyke et al., 2018) that are difficult to represent in GCMs (Pollard, 2010). GrIS precipitation is concentrated at the ice-sheet margins, where steep topographic slopes drive orographic precipitation. The truncated topography used by low resolution GCMs enables moisture to penetrate well into the GrIS interior, manifesting as a positive precipitation bias in the interior (Pollard & Groups, 2000; Van Kampenhout et al., 2019). GrIS ablation areas (marginal regions where the annual SMB is negative) are typically less than 100 km wide and are confined to low-lying areas or regions with low precipitation. These narrow ablation zones are not fully resolved in low-resolution GCMs, and may further degrade the simulated SMB. More recently, GCMs such as CESM and the NASA Goddard Institute for Space Studies GCM (Alexander et al., 2019), have implemented an elevation class downscaling scheme for computing the SMB. The downscaling helps to resolve these narrow ablation zones in GCMs (Sellevold et al., 2019), but large SMB biases remain. For example, CESM, version 2.0 (CESM2) underestimates ablation in the northern GrIS, leading to unrealistic ice advance when run with an interactive ice sheet component (Lofverstrom et al., 2020).

Regional climate models (RCMs) are commonly relied upon to provide more accurate SMB estimates. The limited area domain used by RCMs permits the use of high-resolution grids capable of resolving SMB processes, and can skillfully simulate the GrIS SMB (Box et al., 2004; Fettweis et al., 2013; Mottram et al., 2017; Noël et al., 2018; Rae et al., 2012; Van Angelen et al., 2012). However, unlike GCMs, RCMs are not a freely evolving system, and the atmospheric state must be prescribed at the lateral boundaries of the model domain. The inability of the RCM solution to influence larger-scale dynamics outside the RCM domain (due to the prescribed boundary conditions) severely limits this approach from properly representing the role of the GrIS in the climate system. In addition, the boundary conditions are derived from a separate host model, which introduces inconsistencies due to differences in model design between the host model and the RCM.

In order to retain the benefits of RCMs in a GCM, Van Kampenhout et al. (2019) used the VR capabilities of the SE dycore in CESM, generating grids where Greenland is represented with up to 1/4° resolution, and elsewhere with the more conventional 1° resolution. The simulated SMB compared favorably to the SMB from RCMs and observations. The VR approach is therefore emerging as a powerful tool for simulating and understanding the GrIS and its response to different forcing scenarios.

The SE dycore has been included in the model since CESM version 1, but has been under active development ever since. This includes the switch to a dry-mass vertical coordinate (Lauritzen et al., 2018) and incorporation of an accelerated multi-tracer transport scheme (Lauritzen et al., 2017), made available in CESM2. This paper documents several additional enhancements to the SE dycore as part of the release of CESM2.2. These include three new VR configurations (Figure 2), two Arctic meshes and a Contiguous United-States mesh ($$\mathtt{C}\mathrm{ONUS}$$; featured in Pfister et al. [2020]). While there are dozens of published studies using VR in CESM (e.g., Bambach et al., 2021; Burakowski et al., 2019; Gettelman et al., 2017; Rhoades et al., 2016; Zarzycki et al., 2014), these studies either used development code or collaborated closely with model developers. CESM2.2 is the first code release that contains out-of-the-box VR functionality.

This study compares the representation of Arctic regions using the SE and FV dycores in CESM2.2 (see description below), as these two dycores treat high latitudes (i.e., the pole problem) in different ways. Section 2 documents the grids, dycores, and physical parameterizations used in this study, and also describes the experiments, data sets, and evaluation methods. Section 3 analyzes the results of the experiments, and Section 4 provides a general discussion and conclusions.

Methods

CESM2.2 is a CMIP6-class (Coupled Model Intercomparison Project Phase 6; Eyring et al., 2016) Earth System Model maintained by the National Center for Atmosphere Research. CESM2.2 contains sub-component models for the atmosphere, land, ocean, sea-ice, and land-ice, and can be configured to run with varying degrees of complexity. All simulations described in this study use an identical transient 1979–1998 Atmospheric Model Intercomparison Project (AMIP) configuration, with prescribed monthly sea-surface temperature and sea ice following Hurrell et al. (2008). In CESM terminology, AMIP simulations use the $$\mathtt{F}\mathrm{HIST}$$ computational set and run out of the box in CESM2.2. The land-ice component is not dynamically active in $$\mathtt{F}\mathrm{HIST}$$. However, the surface mass balance is computed by the land model before being passed to the land-ice component; $$\mathtt{F}\mathrm{HIST}$$ includes the necessary functionality to simulate the SMB of the Greenland Ice Sheet.

Dynamical Cores

The atmospheric component of CESM2.2 (Danabasoglu et al., 2020), the Community Atmosphere Model, version 6.3 (CAM6; Craig et al., 2021; Gettelman et al., 2019), supports several different atmospheric dynamical cores. These include dycores on lat-lon grids, such as finite-volume (FV; Lin, 2004) and Eulerian spectral transform (EUL; Collins et al., 2006) models, and dycores built on unstructured grids, including spectral-element (SE; Lauritzen et al., 2018) and finite-volume 3 (FV3; Putman & Lin, 2007) models. This study compares the performance of the SE and FV dycores, omitting the EUL and FV3 dycores. CESM2 runs submitted to CMIP6 used the FV dycore, whereas the SE dycore is often used for global high-resolution simulations (e.g., Bacmeister et al., 2016; Chang et al., 2020; Small et al., 2014) due to its higher throughput on massively parallel systems (Dennis et al., 2012).

Finite-Volume (FV) Dynamical Core

The FV dycore integrates the hydrostatic equations of motion using a finite-volume discretization on a spherical lat-lon grid (Lin & Rood, 1997). The 2D dynamics evolve in floating Lagrangian layers that are periodically mapped to an Eulerian reference grid in the vertical (Lin, 2004). Hyperviscous damping is applied to the divergent modes, and is increased in the top few layers (referred to as a sponge layer) to prevent undesirable interactions with the model top, such as wave reflection (Lauritzen et al., 2011). A polar filter damps computational instability due to the convergence of meridians, permitting a longer time step. It takes the form of a Fourier filter in the zonal direction, with the damping coefficients increasing monotonically in the meridional direction (Suarez & Takacs, 1995). The form of the filter is designed to slow down the propagation of large zonal wave-numbers to satisfy the CFL condition of the shortest resolved wave at some reference latitude.

Spectral-Element (SE) Dynamical Core

The SE dycore integrates the hydrostatic equations of motion using a high-order continuous Galerkin method (Taylor & Fournier, 2010; Taylor et al., 1997). The computational domain is a cubed-sphere grid (obtained by projecting each face of a cube onto a sphere) tiled with quadrilateral elements (see Figure 2). Each element contains a fourth-order basis set in each horizontal direction, with the solution defined at the roots of the basis functions, the Gauss-Lobatto-Legendre quadrature points. This results in 16 nodal points per element, with 12 of the points lying on the (shared) element boundary. Communication between elements uses the direct stiffness summation (Canuto et al., 2007), which applies a numerical flux to the element boundaries to reconcile overlapping nodal values and produce a continuous global basis set.

As with the FV dycore, the dynamics evolve in floating Lagrangian layers that are subsequently mapped to an Eulerian reference grid. A dry mass vertical coordinate was recently implemented for thermodynamic consistency with condensates (Lauritzen et al., 2018). The 2D dynamics have no implicit dissipation, and so hyperviscosity operators are applied to all prognostic variables to remove spurious numerical errors (Dennis et al., 2012). Laplacian damping is applied in the sponge layer.

SE is a next-generation dycore, and is less mature than the FV dycore due to its shorter history. In CESM2.2, the SE numerics have been enhanced relative to CESM2.0 to mitigate spurious noise over topography. These algorithmic changes are described in Appendix A. Future versions of CESM will likely include further optimizations and enhancements to the (more junior) SE dycore.

The SE dycore supports regional grid refinement via its VR configuration, requiring two enhancements over quasi-uniform resolution grids. First, as the numerical viscosity increases with resolution, explicit hyperviscosity relaxes according to the local element size, reducing in strength by an order of magnitude per halving of the grid spacing. A tensor-hyperviscosity formulation is used (Guba et al., 2014), which adjusts the coefficients in two orthogonal directions to more accurately target highly distorted quadrilateral elements. Second, the topography boundary conditions are smoothed in a way that does not excite grid scale modes, and so the NCAR topography software (Lauritzen et al., 2015) has been modified to scale the smoothing radius by the local element size, resulting in rougher topography in the refinement zone.

When using the SE dycore with quasi-uniform grid spacing, the SE tracer transport scheme is replaced with the Conservative Semi-Lagrangian Multi-tracer transport scheme (CSLAM) (Lauritzen et al., 2017). Atmospheric tracers have large, nearly discontinuous horizontal gradients that are difficult to represent with spectral methods, which are prone to oscillatory “Gibbs-ringing” errors (Rasch & Williamson, 1990). CSLAM has improved tracer property preservation and accelerated multi-tracer transport. It uses a separate grid from the spectral-element dynamics, dividing each element into 3 × 3 control volumes with quasi-equal area. The physical parameterizations are computed from the state on the CSLAM grid, which has clear advantages over the original SE dycore in which the physics are evaluated at Gauss-Lobatto-Legendre points (Herrington et al., 2018). CSLAM advection is not available in the VR configuration, which instead uses the standard SE tracer transport scheme with the physics evaluated at Gauss-Lobatto-Legendre points.

Physical Parameterizations

All simulations in this study use the CAM6 physical parameterization package (hereafter referred to as the physics; Gettelman et al., 2019). The physics in CAM6 differs from its predecessors through the incorporation of high-order turbulence closure, Cloud Layers Unified by Binormals (CLUBB; Bogenschutz et al., 2013; Golaz et al., 2002), which jointly acts as a planetary boundary layer, shallow convection, and cloud macrophysics scheme. CLUBB is coupled with the MG2 microphysics scheme (Gettelman & Morrison, 2015; Gettelman et al., 2015), which computes prognostic precipitation and uses classical nucleation theory to represent cloud ice for improved cloud-aerosol interactions. Deep convection is parameterized using a convective quasi-equilibrium mass flux scheme (Neale et al., 2008; Zhang & McFarlane, 1995) and includes convective momentum transport (Richter et al., 2010). Boundary layer form drag is modeled after Beljaars et al. (2004), and orographic gravity wave drag is represented with an anisotropic method informed by the orientation of topographic ridges at the sub-grid scale (the ridge orientation is derived from a high-resolution, global topography data set [Danielson & Gesch, 2011]).

Initial simulations with the SE dycore produced weaker shortwave cloud forcing relative to the tuned finite-volume dycore in the standard CESM2 configuration. The SE dycore in CESM2.2 therefore has two CLUBB parameter changes to provide more realistic cloud forcing and top-of-atmosphere radiation balance. We reduced the width of the sub-grid distribution of vertical velocity $$(\mathtt{c}\mathrm{lubb}\text{\_}\mathrm{gamma}=0.308\to 0.270)$$ and also reduced the strength of the damping for horizontal component of turbulent energy $$(\mathtt{c}\mathrm{lubb}\text{\_}c14=2.2\to 1.6)$$ to increase cloudiness. For a description of how CLUBB parameters impact the simulated climate, see Guo et al. (2015).

Grids

We evaluate model simulations on six different grids in this study (Table 1). The FV dycore is run with nominal 1° and 2° grid spacing, referred to as $$\mathtt{f}09$$ and $$\mathtt{f}19$$, respectively (Figures 1a and 1b). We also run the 1° equivalent of the SE-CSLAM grid, referred to as $$\mathtt{n}\mathrm{e}30\mathrm{pg}3$$ (Figure 1d), where ne refers to a grid with ne × ne quadrilaterial elements per cubed-sphere face, and pg denotes that there are pg × pg control volumes per element for computing the physics. We run an additional 1° SE-CSLAM simulation with the physical parameterizations computed on a grid with 2 × 2 control volumes per element, $$\mathtt{n}\mathrm{e}30\mathrm{pg}2$$ (Figure 1c; Herrington et al., 2019, note CSLAM is still run on the 3 × 3 control volume grid).

Table 1 Grids and Dycores Used in This Study

grid name	Dycore	Δxeq (km)	Δxfine (km)	Δtphys (s)	Cost (25 nodes)	Cost (50 nodes)
$$\mathtt{f}19$$	FV	278	-	1,800	436.66	-
$$\mathtt{f}09$$	FV	139	-	1,800	1534.57	2024.24
$$\mathtt{n}\mathrm{e}30\mathrm{pg}2$$	SE-CSLAM	167	-	1,800	1497.26	1683.97
$$\mathtt{n}\mathrm{e}30\mathrm{pg}3$$	SE-CSLAM	111	-	1,800	1890.48	2090.43
$${\mathtt{n}\mathtt{e}\mathtt{3}\mathtt{0}\mathtt{p}\mathtt{g}\mathtt{3}}^{\ast }$$	SE-CSLAM	111	-	450	-	-
$$\mathtt{A}\mathrm{rctic}$$	SE	111	28	450	15947.41	16675.45
$$\mathtt{A}\mathrm{rctic}-\mathrm{GrIS}$$	SE	111	14	225	40305.03	41036.67

Note. Δx_eq is the average equatorial grid spacing, Δx_fine is the grid spacing in the refined region (if applicable), and Δt_phys is the physics time step. FV refers to the finite-volume dycore, SE the spectral-element dycore, and SE-CSLAM the spectral-element dycore with CSLAM tracer advection. The FV dycore uses lat-lon grids, whereas the SE and SE-CSLAM dycores run on unstructured grids. We use the $$\mathtt{n}\mathrm{e}30\mathrm{pg}3$$ grid for two runs with different values of Δt_phys. The last columns provide the computational costs in core hours per simulated year (CHPSY). The costs are from single-month runs using 25 nodes and 50 nodes on the Cheyenne supercomputer (Computational and Information Systems Laboratory, 2017).

Three VR meshes were developed for the CESM2.2 release to support grid refinement over the Arctic and the United States (Figure 2). This paper serves as the official documentation of these grids. The VR meshes were developed using the software package SQuadgen (https://github.com/ClimateGlobalChange/squadgen). The $$\mathtt{A}\mathrm{rctic}$$ grid is a 1° grid with 1/4° regional refinement over the broader Arctic region (when using degrees to describe the resolution of unstructured grids, we refer to degrees at the equator). The $$\mathtt{A}\mathrm{rctic}-\mathrm{GrIS}$$ grid is identical to the $$\mathtt{A}\mathrm{rctic}$$ grid, but with an additional patch covering the island of Greenland with 1/8° resolution. The $$\mathtt{C}\mathrm{ONUS}$$ grid contains 1/8° refinement over the United States, and 1° everywhere else. The $$\mathtt{C}\mathrm{ONUS}$$ grid is not discussed any further in this paper; see Pfister et al. (2020) for simulations with the $$\mathtt{C}\mathrm{ONUS}$$ grid.

The accuracy of the simulated surface mass balance is expected to be sensitive to grid resolution. Figure 3a shows the average grid spacing over the Greenland Ice Sheet (GrIS hereafter) in all six grids, as well as two grids pertaining to the Regional Atmospheric Climate Model (RACMO; Noël et al., 2018, 2019), which are used for validation (Table 2). The $$\mathtt{n}\mathrm{e}30\mathrm{pg}2$$ grid has the coarsest representation with an average grid spacing (Δx) of Δx = 160 km, and the $$\mathtt{A}\mathrm{rctic}-\mathrm{GrIS}$$ grid has the highest resolution with an average grid spacing of Δx = 14.6 km, similar to the 11 km grid spacing of the RACMO2.3 grid. The $$\mathtt{n}\mathrm{e}30\mathrm{pg}3$$ grid has an average Δx = 111.2 km, substantially coarser than the $$\mathtt{f}09$$ grid, with an average Δx = 60 km. Although $$\mathtt{n}\mathrm{e}30\mathrm{pg}3$$ and $$\mathtt{f}09$$ have similar average grid spacing over the entire globe, and comparable computational costs (Table 1), the convergence of meridians on the FV grid enhances the resolution over the GrIS. The $$\mathtt{A}\mathrm{rctic}$$ grid has an average grid spacing of Δx = 27.8 km.

Figure 3. The spatial properties of the GrIS as represented by different grids in this study. (Left) approximate average grid spacing over GrIS, (right) GrIS area error, computed as the relative differences from the reference data set described in the text. Also included are the two RACMO grids used for validation (see Table 2).

Table 2 Description of Validation Data Sets Used in This Study

Data product	Years used in this study	Resolution	Citation
ERA5	1979–1998	1/4°	Copernicus (2019)
CERES-EBAF ED4.1	2003–2020	1°	Loeb et al. (2018)
CALIPSO-GOCCP	2006–2017	2°	Chepfer et al. (2010)
MODIS-MAR	2000–2013	25 km	Tedesco and Alexander (2013)
RACMO2.3	1979–1998	11 km	Noël et al. (2015)
RACMO2.3p2	1979–1998	5.5 km	Noël et al. (2019)

The physics time step depends on the grid resolution. Increased horizontal resolution permits faster vertical velocities that reduce characteristic time scales, so the physics time step is reduced to avoid large time truncation errors (Herrington & Reed, 2018). The $$\mathtt{A}\mathrm{rctic}$$ and $$\mathtt{A}\mathrm{rctic}-\mathrm{GrIS}$$ grids are run with a 4× and 8× reduction in physics time step relative to the default 1,800 s time step used in the 1° and 2° grids (Table 1).

All grids and dycores in this study use 32 hybrid pressure-sigma levels in the vertical, with a model top of 2 hPa or about 40 km. However, any grid or dycore can in principle be run with a higher model top or finer vertical resolution.

Computational Costs

The last columns of Table 1 provide cost estimates for the different grids and dycores. The costs, expressed as core hours per simulated year (CHPSY), are taken from single-month runs of $$\mathtt{F}\mathrm{HIST}$$, with no i/o, and using 25 nodes (900 tasks) and 50 nodes (1,800 tasks) on the Cheyenne supercomputer (Computational and Information Systems Laboratory, 2017). For a typical multi-decadal climate simulation at 1° resolution, 25 nodes are on the low side. However, this is the largest number of tasks supported by the $$\mathtt{f}19$$ grid, and we fixed the number of tasks across all grids for the purpose of comparing costs. We also show costs for 50 nodes, excluding $$\mathtt{f}19$$, to provide a more practical estimate for longer climate integrations. Other approaches for comparing costs across different grids and dycores, for example, holding fixed the number of grid columns per task, are beyond the scope of this study.

The cheapest grid is the $$\mathtt{f}19$$ grid at 436.66 CHPSY, as this is the only grid with 2° dynamics. The $$\mathtt{f}09$$ grid costs 1534.57 CHPSY using 25 nodes, which is noticeably cheaper than $$\mathtt{n}\mathrm{e}30\mathrm{pg}3$$ at 1890.48 CHPSY. The $$\mathtt{n}\mathrm{e}30\mathrm{pg}2$$ grid is 20% cheaper than the $$\mathtt{n}\mathrm{e}30\mathrm{pg}3$$ grid, in both the 25 node and 50 node runs, consistent with previous estimates (Herrington & Reed, 2018). The FV model is known to be cheaper than SE at small core counts, whereas SE is more efficient than FV at large core counts due to its improved scalability (Dennis et al., 2005, 2012). In the more conventional 50-node runs, $$\mathtt{f}09$$ costs are more similar to $$\mathtt{n}\mathrm{e}30\mathrm{pg}3$$, due to a 30% cost increase in $$\mathtt{f}09$$ relative to the 25-node run (Table 1). The $$\mathtt{A}\mathrm{rctic}$$ grid is an order of magnitude more expensive than the lat-lon and quasi-uniform grids, at about 16k CHPSY, whereas the $$\mathtt{A}\mathrm{rctic}-\mathrm{GrIS}$$ grid adds another factor of ∼2 (40k CHPSY). All timing numbers are from runs without threading. The $$\mathtt{f}09$$ grid is the only grid that runs out-of-the-box with threading; holding the number of tasks fixed reduces CHPSY by 4%–6.5% reduction compared to runs without threading.

Simulated Surface Mass Balance (SMB)

CESM simulates the GrIS SMB as the sum of ice accumulation and ice ablation. The latter includes sublimation, as well as liquid runoff from ice melt. Liquid precipitation and runoff may also contribute to ice accumulation by penetrating pore spaces in the snowpack/firn layer and freezing, forming ice lenses. The relevant SMB processes are represented by different CESM components, but it is the Community Land Model, version 5 (CLM; Lawrence et al., 2019), that aggregates these processes and computes the SMB.

CLM runs on the same grid as the atmosphere, and uses a downscaling technique to account for sub-grid variability in SMB. In short, the ice sheet patch in a CLM grid cell is subdivided into 10 elevation classes (ECs), each with a distinct surface energy balance and SMB. The area fraction of each EC is computed from the CISM initial conditions, which are based on a high-resolution data set of the observed, modern extent and thickness of the GrIS (Morlighem et al., 2014). For configurations with a dynamically active ice sheet, the area fractions are continuously updated throughout the run to reflect the evolving ice sheet geometry in CISM. The near-surface air temperature, humidity, and air density are calculated for each EC using an assumed lapse rate and the elevation difference from the grid-cell mean. Precipitation from CAM is repartitioned into solid or liquid based on the surface temperature of the EC; precipitation falls as snow for temperatures between T < −2°C, as rain for T > 0°C, and as a linear combination of rain and snow for temperatures between −2° and 0°C.

Changes in ice depth, not snow depth, count toward the SMB. Snow accumulation in each EC is limited to a depth of 10 m liquid water equivalent. Any snow above the 10 m cap contributes to ice accumulation, and refreezing of liquid water within the snowpack is an additional source of ice. Surface ice melting (after melting of the overlying snow) yields a negative SMB. Integrating over all ECs, weighting by the area fractions, provides a more accurate SMB than would be found using the grid-cell mean elevation. For a more detailed description of how the SMB is computed in CESM, we refer the reader to Lipscomb et al. (2013), Sellevold et al. (2019), van Kampenhout et al. (2020), and Muntjewerf et al. (2021).

Since snow in the accumulation zone must reach the cap to simulate a positive SMB, the snow depths on the VR grids were spun up by forcing CLM in standalone mode, cycling over data from a 20-year $$\mathtt{A}\mathrm{rctic}\ \mathtt{F}\mathrm{HIST}$$ simulation (with perpetual 1979 boundary conditions) for about 500 years. The 1°–2° lat-lon and quasi-uniform unstructured grids are initialized with the SMB from an existing $$\mathtt{f}09$$ spun-up initial condition.

Validation Data Sets

We use several validation data sets (Table 2) to assess the performance of the simulations. The ERA5 reanalysis product (Copernicus, 2019) is used to validate the large-scale circulation and extreme precipitation events. Clouds and radiation fields are validated using remote sensing products, the CERES-EBAF ED4.1 (Loeb et al., 2018) and the CALIPSO-GOCCP (Chepfer et al., 2010) data sets, respectively. We also use a MODIS GrIS surface albedo product, mapped to a 25 km regional climate model grid (Tedesco & Alexander, 2013).

SMB data sets are gathered from multiple sources. RACMO, version 2.3 11 km (RACMO23; Noël et al., 2015) and version 2.3p2 5.5 km (RACMO2.3p2; Noël et al., 2018, 2019) are RCM simulations targeting Greenland, forced by ERA renalyses products at the domain's lateral boundaries. The RACMO simulations have been shown to perform skillfully against observations and are often used as modeling targets (e.g., Evans et al., 2019; van Kampenhout et al., 2020). The RACMO data sets are used along with the CERES-EBAF product to validate the radiative fluxes around Greenland.

In-situ SMB (snow pit and ice cores) and remote sensing data sets (e.g., IceBridge radar accumulation data set) are maintained in The Land Ice Verification and Validation toolkit (LIVVkit), version 2.1 (Evans et al., 2019). However, these point-wise measurements are difficult to compare to model output spanning several different grids, especially since the SMB from each elevation class is not available from the model output. We used a nearest-neighbor technique for an initial analysis, which showed that the model biases are similar to those computed using the RACMO data sets. Because of the uncertainty of comparing gridded fields to point-wise measurements, and the lack of information added with regard to model biases, we omitted these data sets from our analysis.

SMB Analysis

We seek to integrate SMB components over a GrIS ice mask and to diagnose their contributions to the GrIS mass budget. However, the ice masks vary across the grids, especially in comparison to the RACMO3.2 ice mask, whose total area is about 3% less than that of the reference data set (i.e., the GrIS initial conditions in CISM; Figure 3b). CLM's data set creation tool generates the model ice mask by mapping the reference data set to the target grid using the Earth System Modeling Framework (ESMF) first-order conservative remapping algorithm (Team et al., 2021). The figure suggests that the mapping errors are less than 1.5% across the CESM2.2 grids. The area errors in Figure 3b may seem small, but even 1–2% area differences can lead to large differences in integrated SMB (Hansen et al., 2022).

We have taken a common-ice-mask approach by mapping all model fields to the lowest-resolution grids, that is, the $$\mathtt{f}19$$ and $$\mathtt{n}\mathrm{e}30\mathrm{pg}2$$ grids, and integrating over these low-resolution ice masks. The use of low-resolution common ice masks is a conservative decision, and is justified because we seek to use first-order remapping algorithms to map fields to the common ice mask, which is not generally reliable when mapping to a higher-resolution grid than the source grid. We use two remapping algorithms: ESMF first-order conservative and the TempestRemap (Ullrich & Taylor, 2015) high-order monotone algorithm. Since mapping errors are sensitive to grid type, we evaluate all quantities on both common ice masks, the $$\mathtt{f}19$$ and $$\mathtt{n}\mathrm{e}30\mathrm{pg}2$$ masks. Thus, we evaluate an integrated quantity on a given grid up to four times to estimate the uncertainty due to differences in grid type and remapping algorithms.

The SMB is expressed in a form that is agnostic of water phase, a total water mass balance, to facilitate comparisons across different grids with different ice masks and to increase consistency with the variables available in the RACMO data sets. The SMB for total water can be expressed as: [Image Omitted. See PDF]where all terms have consistent sign conventions (positive values contribute mass, and negative values reduce mass). The accumulation source term refers to the combined solid and liquid precipitation, runoff refers to the liquid water sink, and evaporation/sublimation is the vapor sink. Since the runoff term aggregates many processes, we isolate the melting contribution by also tracking the combined melt of snow and ice.

The total water SMB (Equation 1) is different from the SMB internally computed by CLM and described in Section 2.5, which only tracks ice mass. We do not use CLM's internally computed SMB in this study. Rather, we use the components of the internally computed SMB to construct the total water SMB.

We consider two approaches for mapping and integrating the SMB components over the common ice masks:

1. Map the grid-cell mean quantities to the common grid, and integrate the mapped fields over the common ice masks.

2. Map the patch-level quantities (i.e., the state over the ice fractional component of the grid cell) to the common grid, and integrate the mapped fields over the common ice masks.

Note that we are mapping to low-resolution grids that have larger GrIS areas than the source grids (Figure 3b). Since the components of Equation 1 are not confined to the ice mask, method 1 reconstructs the SMB over the portion of the common ice mask that is outside the ice mask on the source grid. While this may be a an acceptable way to reconstruct the mass source terms over different ice masks, ice melt is zero outside the source ice mask, and so method 1 will underestimate the mass sink term. This underestimation is systematic in method 2, where all variables are exclusive to the ice mask; mapping to a lower-resolution grid will dilute a field of non-zero values over the ice mask with a field of zeros outside the ice mask. However, patch-level values for processes exclusive to the ice mask (e.g., ice melt) will on average have larger magnitudes than the grid-mean quantities used in method 1.

The different error characteristics of the two methods are used to further diversify the ensemble. Each of the four regridding combinations (with conservative and high-order remapping to the $$\mathrm{f}19$$ and $$\mathtt{n}\mathrm{e}30\mathrm{pg}2$$ grids) are repeated with each method, resulting in (up to) eight values for each integrated quantity. Unfortunately, the patch-level values of evaporation/sublimation are not available from the model output, and we estimate their contribution by zeroing out the field for grid cells that have no ice, prior to mapping to the common ice mask. This will degrade the SMB estimates using method 2, but we are more interested in characterizing the behavior of individual processes across grids and dycores, expressed as the components of the SMB, rather than the SMB itself.

Results Tropospheric Temperatures

Before delving into the simulated Arctic climate conditions, we assess the global mean differences between the various grids and dycores. Figure 4 shows 1979–1998 annual mean, zonal mean height plots of temperature, expressed as differences between 1° and 2° lat-lon and quasi-uniform unstructured grids and dycores. The $$\mathtt{f}09$$ grid is warmer than the $$\mathtt{f}19$$ grid, primarily in the mid-to-high latitudes throughout the depth of the troposphere. This is a common response to increasing horizontal resolution in GCMs (Pope & Stratton, 2002; Roeckner et al., 2006). Herrington and Reed (2020) have shown that this occurs in CAM due to higher resolved vertical velocities which, in turn, generate more condensational heating in the CLUBB macrophyiscs. The right panel in Figure 4a supports this interpretation, showing an increase in the climatological CLUBB heating at all latitudes in the $$\mathtt{f}09$$ grid, but with the largest increase in the mid-latitudes.

Figure 4. 1979–1998 annual mean temperature (left column) and CLUBB temperature tendencies (right column) in zonal mean height space, expressed as differences between the various 1°–2° grids. The thick gray and magenta lines are the tropopause for the control run and the test run, respectively.

As the SE dycore is less diffusive than the FV dycore, the resolved vertical velocities are larger in the SE dycore, and so the $$\mathtt{n}\mathrm{e}30\mathrm{pg}3$$ troposphere is modestly warmer than $$\mathtt{f}09$$ (Figure 4b). The stratosphere responds differently, with $$\mathtt{n}\mathrm{e}30\mathrm{pg}3$$ much cooler than $$\mathtt{f}09$$ in the mid-to-high latitudes. Figure 4c also shows small temperature differences between $$\mathtt{n}\mathrm{e}30\mathrm{pg}3$$ and $$\mathtt{n}\mathrm{e}30\mathrm{pg}2$$, with $$\mathtt{n}\mathrm{e}30\mathrm{pg}3$$ slightly warmer near the tropopause at high latitudes. This is consistent with the similar climates found for these two grids by Herrington et al. (2019).

Comparing the VR grids to the lat-lon and quasi-uniform grids is complicated because we simultaneously increase the resolution and reduce the physics time-step, both of which influence the solution (Williamson, 2008). We therefore run an additional $$\mathtt{n}\mathrm{e}30\mathrm{pg}3$$ simulation with the shorter physics time step used in the $$\mathtt{A}\mathrm{rctic}$$ grid (450 s), referred to as $${\mathtt{n}\mathtt{e}\mathtt{3}\mathtt{0}\mathtt{p}\mathtt{g}\mathtt{3}}^{\ast }$$ (Table 1). Figure 5a shows the difference between $${\mathtt{n}\mathtt{e}\mathtt{3}\mathtt{0}\mathtt{p}\mathtt{g}\mathtt{3}}^{\ast }$$ and $$\mathtt{n}\mathrm{e}30\mathrm{pg}3$$ for climatological summer temperatures in zonal-mean height space. The troposphere is warmer with the reduced time step, and the mechanism is similar in that the shorter time step increases resolved vertical velocities (not shown) and CLUBB heating (right panel in Figure 5a). Figure 5b shows the difference in climatological summer temperature between the $$\mathtt{A}\mathrm{rctic}$$ grid and the $${\mathtt{n}\mathtt{e}\mathtt{3}\mathtt{0}\mathtt{p}\mathtt{g}\mathtt{3}}^{\ast }$$ grid. With the same physics time step, the greater condensational heating and warmer temperatures are confined to the refined Arctic region.

Figure 5. As in Figure 4 but for the short-time-step experiment and the VR grids. The fields plotted are for the climatological northern hemisphere summer. We focus on summer because that is when the resolution response is largest, and the refined regions are located in the northern hemisphere.

There is not perfect alignment of the CLUBB heating anomalies and the temperature anomalies in Figures 4 and 5. This leaves some uncertainty as to whether CLUBB is solely responsible for the temperature anomalies caused by changing grid resolution or dycore, as asserted above. In the deep tropics, gravity waves are expected to propagate temperature anomalies far from their heating source, and so some misalignment is expected for this region. In the mid-latitudes, we speculate that the misalignment is due to averaging the CLUBB heating over all time. Herrington and Reed (2020) have shown that only the CLUBB heating coinciding with upward resolved vertical velocities tends to align spatially with the temperature response, whereas heating anomalies associated with descending grid columns have no alignment with the temperature response. Therefore the average over all time will dilute the heating signal, potentially leading to misalignments.

Figure 5c shows that the $$\mathtt{A}\mathrm{rctic}-\mathrm{GrIS}$$ grid is much warmer than the $$\mathtt{A}\mathrm{rctic}$$ grid in the Arctic summer. This may be due, in part, to the shorter physics time step, but the temperature response is too large to be explained by enhanced condensational heating from CLUBB alone. This summer warming appears to be a result of variations in the stationary wave pattern, with a swath of anomalous southerly winds to the west of Greenland (not shown). This dynamic response is interesting, because other than the physics time step, the only difference between the $$\mathtt{A}\mathrm{rctic}-\mathrm{GrIS}$$ and $$\mathtt{A}\mathrm{rctic}$$ runs is the doubling of resolution over Greenland. This behavior will be explored further in a future study.

Keeping our focus on the Arctic region, and in particular Greenland, it is useful to understand summer temperature biases due to their control on ice and snow melt over the GrIS (Ohmura, 2001). Figure 6 shows the 1979–1998 lower-troposphere summer temperature bias relative to ERA5, computed by equating a layer mean virtual temperature with the 500–1,000 hPa geopotential thickness. The results are consistent with the zonal mean height plots; increasing resolution from $$\mathtt{f}19$$ to $$\mathtt{f}09$$ warms the climate, and the 1° SE grids are warmer than the FV grids. The FV summer temperatures are persistently colder than ERA5, whereas the 1° SE grids are not as cold, and are actually warmer than ERA5 at high-latitudes, north of 75°. All grids show a north-south gradient in bias over Greenland, with the summer temperature bias more positive for the northern part of the ice sheet. This pattern is also evident in the near-surface temperature bias over Greenland (not shown).

Figure 6. 1979–1998 lower troposphere, northern hemisphere summer virtual temperature biases, computed as the difference from ERA5. Lower troposphere layer mean virtual temperature is derived from the 1,000–500 hPa geopotential thickness, using the hypsometric equation. Differences are computed after mapping the ERA5 data to the finite-volume grids since the geopotential field is only available on the output tapes in the spectral-element runs that have been interpolated to the f09 $\mathtt{f}09$ grid, inline.

The $$\mathtt{A}\mathrm{rctic}$$ grid has summer temperatures similar to the 1° SE grids, but is slightly warmer over northern Eurasia and the North Pole (Figure 6). An anomalous cooling patch forms to the west of Greenland, centered over Baffin Island. The $$\mathtt{A}\mathrm{rctic}-\mathrm{GrIS}$$ grid is warmer than the $$\mathtt{A}\mathrm{rctic}$$ grid over most of the Arctic, but with a similar spatial pattern of summer temperature bias.

Some of these temperature differences may be related to differences in summer cloudiness. Figure 7 shows the summer shortwave cloud forcing bias in the six runs, using the CERES product. Shortwave cloud forcing quantifies the impact of clouds on shortwave radiation, taken as the difference between all-sky and clear-sky shortwave radiative fluxes at the top of the atmosphere. A negative bias corresponds to excessive reflection and cooling. The lat-lon and quasi-uniform grids have similar biases, with the clouds reflecting 20–40 W/m² too much shortwave radiation over a wide swath of the Arctic, primarily the land masses. There is also a halo of positive bias (clouds not reflective enough) around the ocean perimeter of Greenland. The $$\mathtt{A}\mathrm{rctic}$$ grid has much smaller cloud forcing biases over the Arctic land masses, but is still too reflective over Alaska, the Canadian Archipelago, and parts of Eurasia. Compared to the $$\mathtt{A}\mathrm{rctic}$$ grid, the $$\mathtt{A}\mathrm{rctic}-\mathrm{GrIS}$$ grid vastly reduces the cloud forcing bias over Eurasia, and also improves the bias over North America. In both VR grids, the halo of positive shortwave cloud forcing bias around the perimeter of Greenland is absent.

Figure 7. 1979–1998 Northern Hemisphere summer shortwave cloud forcing bias, relative to the CERES-EBAF gridded data set. Shortwave cloud forcing is defined as the difference between all-sky and clear-sky net shortwave fluxes at the top of the atmosphere. Differences are computed after mapping all model output to the 1° CERES-EBAF grid.

The summer cloud forcing biases are consistent with the summer temperature biases in Figure 6—regions where clouds are too reflective coincide with regions that are too cold. While we have not quantified the contribution of cloud biases to the cooler Arctic temperatures, shortwave radiation is a crucial component of the Arctic energy budget during summer. The shortwave cloud forcing biases are on the order of 10 W/m2, which is a significant fraction of the total absorbed shortwave radiation during Arctic summer (Serreze et al., 2007) and is therefore likely a factor contributing to the cooler temperatures.

Clouds and Precipitation Over Greenland

In addition to summer temperatures, shortwave radiation is an important determinant of snow and ice melt. Figure 8 shows the summer incident shortwave radiation bias at the surface over Greenland and surrounding seas. The top panel shows the bias relative to the RACMO2.3p2 data set, and the middle panel relative to the CERES data set. The halo of excessive incident shortwave radiation around the coasts of Greenland is apparent for both data sets in relation to the coarser grids, consistent with the shortwave cloud forcing biases in Figure 7.

Figure 8. 1979–1998 northern hemisphere summer surface incident shortwave radiation bias (W/m2), computed as the difference (top) from RACMO2.3p2, and (middle) the CERES-EBAF data set, and the (bottom) total cloud fraction bias relative to the CALIPSO data set. CERES and CALIPSO differences are found by mapping the model output to their 1° and 2° grids, respectively, and differences in the bottom panel are computed after mapping the RACMO2.3p2 data set to the individual model grids. Note that the averaging period for the CALIPSO-GOCCP and CERES-EBAF panels, 2006–2017 and 2003–2020, respectively, are different from the averaging period for the model results.

The ice sheet interior receives too little shortwave radiation in the coarser grids. On the VR grids, both the interior shortwave deficit and the excessive shortwave around the ocean perimeter are improved. This suggests that the coarse-grid clouds are too thick in the interior of Greenland and too thin around the perimeter, which is consistent with the total summer cloud fraction bias, computed from the CALIPSO cloud data set and shown in the bottom panel of Figure 8. Note that total cloud fraction characterizes the cloud field at all vertical levels, and attenuates the changes arising from any single layer due to the occurrence of overlapping clouds at other levels. The VR grids exhibit an overall improvement in total cloud fraction bias, relative to the coarse grids.

The top panel of Figure 9 shows the annual climatological mean precipitation bias over the GrIS, expressed as the fractional difference from the RACMO2.3p2 solution. The coarse 1°–2° grids have large, positive biases centered over the southern dome. The $$\mathtt{A}\mathrm{rctic}$$ grid reduces this bias substantially, and the $$\mathtt{A}\mathrm{rctic}-\mathrm{GrIS}$$ grid reduces it further, with precipitation centers migrating from the interior toward the margins.

Figure 9. 1979–1998 (top) annual precipitation and (bottom) ice/snow melt biases relative to RACMO2.3p2, evaluated on the native model grids. The precipitation biases are expressed as fractional changes, whereas the melt biases are absolute changes (mm/yr). In the bottom panel, the Rignot and Mouginot (2012) basin boundaries are shown in gray for each model grid. Note that Figure 11 uses the basin boundaries for the two common ice masks, shown in the f19 $\mathtt{f}19$ and ne30pg2 $\mathtt{n}\mathrm{e}30\mathrm{pg}2$ panels, in computing the basin-scale integrals. Blue lines in the f19 $\mathtt{f}19$ panel show the location of the two transects plotted in Figure 12.

The more accurate representation of orographic precipitation in the VR grids is consistent with the cloud and radiation biases, cf. Figures 7–9. The agreement of the cloud, radiation and precipitation biases in and around Greenland from multiple independent data sets indicates that the biases are a robust feature of the coarser grids. The reduced biases in the VR grids suggest that the deficiencies of the coarse grids are due to insufficient horizontal resolution, consistent with previous findings that coarse GCMs have large, positive precipitation biases over Greenland (Pollard & Groups, 2000; Van Kampenhout et al., 2019).

Greenland Surface Mass Balance

Table 3 shows the 1979–1998 climatological SMB components for each grid, compared with RACMO. The values in the table are averages over the ensemble of mapping methods to the common ice masks described in Section 2.7, and the RACMO values refer to the average of both RACMO ensembles. Table 3 also contains (in parenthesis) the SMB components derived from evaluating the integrals on each model's native grid and ice mask. For the accumulation term in the quasi-uniform SE grids, the native grid values are actually less than common ice mask values, likely due to the inflationary impacts of method 1 discussed in Section 2.7. For the total melt term, the common ice mask values aren't nearly as dissipative as expected, with only the $$\mathtt{A}\mathrm{rctic}-\mathrm{GrIS}$$ grid showing a large reduction in melt rates relative to the native grid values.

Table 3 1979–1998 Surface Mass Balance of the Greenland Ice Sheet in Gt/yr

Grid name	Accumulation	Total melt	Runoff	Evap + subl	SMB
RACMO	681.7 (733.5)	−318.6 (−436.4)	−189.1 (−258.5)	−34.5 (−38.8)	458.1 (436.2)
$$\mathtt{n}\mathrm{e}30\mathrm{pg}2$$	1013.6 (986.9)	−560.6 (−561.4)	−410.0 (−380.6)	−39.0 (−37.2)	564.6 (569.1)
$$\mathtt{n}\mathrm{e}30\mathrm{pg}3$$	930.6 (909.7)	−583.1 (−543.6)	−403.2 (−337.9)	−39.1 (−37.6)	488.3 (534.2)
$$\mathtt{f}19$$	888.2 (919.7)	−446.3 (−472.1)	−304.5 (−307.6)	−41.9 (−42.9)	541.8 (569.2)
$$\mathtt{f}09$$	878.6 (891.6)	−419.5 (−401.5)	−274.4 (−230.6)	−42.3 (−42.4)	561.9 (618.6)
$$\mathtt{A}\mathrm{rctic}$$	786.3 (823.5)	−361.2 (−350.5)	−231.4 (−199.1)	−47.0 (−48.7)	507.9 (575.7)
$$\mathtt{A}\mathrm{rctic}-\mathrm{GrIS}$$	695.9 (751.6)	−471.6 (−513.7)	−297.7 (−301.8)	−53.2 (−57.0)	345.0 (392.8)

Note. Values shown are using the common ice mask approach described in the methods section, whereas values in parentheses are from integrating over the native grid and ice mask. Evap + subl refers to the combined evaporation and sublimation term.

The coarse grids are characterized by too much precipitation and too much melting and runoff, compared with RACMO. The total SMB on coarse grids therefore has smaller errors than the individual components (Table 3), because large errors in the source and sink terms offset one another when added together. Such compensating errors underscore the importance of understanding whether a model is getting the right SMB for the right reasons.

Figure 10 shows time series of annually integrated precipitation and snow/ice melt over the GrIS for the various different grids and dycores, and RACMO in black. The 1979–1998 climatological mean values from Table 3 are shown as circles on the right side of the panels. The 1°–2° lat-lon and quasi-uniform grids have positive precipitation biases, whereas the VR grids have the smallest biases, with precipitation comparable to RACMO. The $$\mathtt{f}19$$ and $$\mathtt{f}09$$ grids perform similarly, with +200 Gt/yr bias, whereas $$\mathtt{n}\mathrm{e}30\mathrm{pg}3$$ is biased by about +250 Gt/yr and $$\mathtt{n}\mathrm{e}30\mathrm{pg}2$$ by +330 Gt/yr.

Figure 10. Time-series of annual (solid + liquid) precipitation (top) and annual runoff (bottom) integrated over the Greenland Ice Sheet for all six simulations and compared to the RACMO data sets. The time-series were generated using the common ice mask approach, which results in up to four ensembles, with the mean value given by the solid line and shading spanning the extent of the ensemble members.

The combined annual snow/ice melt shown in the bottom panel of Figure 10 indicates that the $$\mathtt{A}\mathrm{rctic}$$ grid simulates the most realistic melt rates, with the other grids having more melt than RACMO. The $$\mathtt{A}\mathrm{rctic}-\mathrm{GrIS}$$ grid over-predicts melting by about 150 Gt/yr. This is likely due to an anomalously warm lower troposphere during the summer, relative to the $$\mathtt{A}\mathrm{rctic}$$ run (Figure 6). The $$\mathtt{f}19$$ and $$\mathtt{f}09$$ melting rates are improved over $$\mathtt{A}\mathrm{rctic}-\mathrm{GrIS}$$, overestimating melt by only 100–120 Gt/yr. The SE grids have the largest positive melt bias, between 240 and 265 Gt/yr.

To illustrate the regional behavior of the SMB components, Figure 11 shows the precipitation and combined snow/ice melt integrated over the basins defined by Rignot and Mouginot (2012). The uncertainty due to differences in basin area is larger than for GrIS-wide integrals, owing to the differences in basin boundaries on the common ice masks, which are shown in the $$\mathtt{f}19$$ and $$\mathtt{n}\mathrm{e}30\mathrm{pg}2$$ panels of Figure 9. Nonetheless, the regional totals in Figure 11 correctly show that the southeast and southwest basins have the most accumulation. In all basins, accumulation decreases monotonically with increasing grid resolution, though with some exceptions. The $$\mathtt{A}\mathrm{rctic}-\mathrm{GrIS}$$ grid simulates less precipitation than RACMO in the central-east and southeast basins, and is closest of all grids to RACMO in the large southwest basin.

Figure 11. 1979–1998 basin integrated components of the SMB; (top) precipitation, (middle) ice/snow melt and (bottom) ice/snow melt estimated from the PDD method. Whiskers span the max/min of the four ensemble members generated from the common-ice-mask approach. Basin definitions are after Rignot and Mouginot (2012), and are found on the common ice masks using a nearest neighbor approach, and shown in Figure 9.

The basin-integrated melt rates in Figure 11 depend on the dycore. The quasi-uniform SE grids have the largest positive biases in all basins. The $$\mathtt{A}\mathrm{rctic}-\mathrm{GrIS}$$ grid is a close second, while the FV grids have systematically smaller melt-rates and melt-rate biases. The “second-place” standing of $$\mathtt{A}\mathrm{rctic}-\mathrm{GrIS}$$ is somewhat unexpected, as this grid has the warmest lower-troposphere summer temperatures (Figure 6) and greatest incident shortwave radiation (Figure 8), yet it has less melting than the quasi-uniform SE grids.

Lower-troposphere temperature is not a strict proxy for melting; for example, it may not capture microclimate effects as a result of a better representation of the low-elevation ablation zones. The Positive Degree-Day temperature-based melt index (PDD; Braithwaite, 1984), which accumulates the near-surface temperature in °C for days with temperature above freezing, is a more accurate proxy for melting. PDD is nonlinear in mean monthly temperature (Reeh, 1991). We compute PDD from monthly mean 2-m temperature using the method of Calov and Greve (2005), assuming a fixed monthly mean standard deviation of 3°C and a degree-day factor of 5 mm d⁻¹ °C⁻¹. This specific degree-day factor lies between typical values reported for snow and ice, in order to apply the PDD method to estimate the combined snow and ice melt.

Figure 11c shows the basin-integrated PDD melt estimate. In the large southeast and southwest basins (and all the other western basins), the $$\mathtt{n}\mathrm{e}30\mathrm{pg}3$$ grid has larger PDD-based melt than the $$\mathtt{A}\mathrm{rctic}-\mathrm{GrIS}$$ grid. The FV grids also have large PDD-based melt in the southwest basin, relative to $$\mathtt{A}\mathrm{rctic}-\mathrm{GrIS}$$. The PDD plots indicate that the relationship between temperature and melt is not well approximated by the summer lower-troposphere temperatures in Figure 6.

The bottom panel of Figure 9 presents the biases in the combined ice/snow melt as map plots. These plots show that the largest melt biases are on the southeast and northwest coasts, where large coarse-grid cells overlap with the ocean. One possibility is that these problematic grid cells are situated at lower elevations than the true ice sheet surface, leading to a warm bias and too much melt. Figure 12 shows the representation of the ice sheet surface along two transects on the different grids, compared to the high-resolution data set used to generate CAM topographic boundary conditions (Danielson & Gesch, 2011; Lauritzen et al., 2015). The two transects are shown in Figure 9: the east-west “K-transect” in southwest Greenland and a transect extending from the central dome down to the Kangerlussuaq glacier on the southeast coast. The 1°–2° grids are noticeably coarse, with only a handful of grid cells populating the transect. The $$\mathtt{f}09$$ grid is a bit of an exception since the grid cells narrow in the zonal direction at high latitudes, and so more grid cells fit into the east-west transects. The VR grids more accurately reproduce the steep margins of the ice sheet, capturing the characteristic parabolic shape of the GrIS margin.

Figure 12. Model surface elevation along the (top) K-transect, and (bottom) a transect spanning the central dome down to the Kangerlussuaq glacier in southeast Greenland, for all model grids. The GMTED reference surface is a 1 km surface elevation data set (Danielson & Gesch, 2011) used for generating the CAM topographic boundary conditions.

The transects in Figure 12 show that the ice sheet surface on the coarse grids is not systematically lower than the true surface in ablation zones. Rather, the smoothing and flattening of the raw topography, necessary to prevent the model from exciting grid-scale numerical modes, causes the lower-elevation ablation zones to extend beyond the true ice sheet margin, causing the modeled ablation zones (which must reside within the ice sheet mask) to be elevated relative to the actual ice surface. The $$\mathtt{f}19$$ grid has both the smoothest topography and the flattest ice sheet since its dynamics are coarsest, whereas the $$\mathtt{f}09$$, $$\mathtt{n}\mathrm{e}30\mathrm{pg}2$$ and $$\mathtt{n}\mathrm{e}30\mathrm{pg}3$$ grids have similar dynamical resolution and use identical smoothing. This suggests that coarser models will tend to elevate the ablation zones relative to where they should be, which may be expected to cause anomalous (adiabatic) cooling and depressed melt rates. This is opposite the melt bias in the coarse grid simulations. Also, the EC downscaling should be able to correct for situations where model ice surface is elevated relative to the true ice surface, although for very large elevation differences or errors, the EC scheme may not be adequate.

Figure 12 also shows the ice margin boundary, illustrating that the ablation zone lies in a narrow horizontal band where the ice sheet rapidly plunges to sea-level. Due to this abrupt transition, coarse grids will commonly represent the ablation zone with grid cells containing mixtures of ice-covered and ice-free regions. We hypothesize that coarser models have larger melt biases because summer melting is confined to these mixed ice/land/ocean grid cells. CLM deals with land heterogeneity in a complex and sophisticated manner, but CAM only sees the homogenized state after area-averaging over the sub-grid mixture. Thus, warm ice-free land patches in a grid cell may unduly influence the climate over the entire grid cell, causing a warm bias over the ice-covered patch and more melting.

Figure 13 shows the mean melt bias, relative to both RACMO data sets, conditionally sampled based on grid cell ice fraction in the GrIS region. Errors are computed after mapping the melt rates to the common ice masks using different methods, described in Section 2.7. The grid cell ice fractions therefore pertain to ice fractions on the low-resolution common ice masks. Also shown are the ±1 standard deviation of the biases for each bin. The figure shows that coarser grids can be characterized by a monotonic increase in melt bias as the grid cell ice fraction decreases. The VR grids have the smallest melt biases for small grid cell ice fractions (less than 50%), the quasi-uniform SE grids and $$\mathtt{f}19$$ have the largest melt biases, and the $$\mathtt{f}09$$ grid melt biases lie between these two cases. Figure 13 generally supports our hypothesis that the prevalence of mixed-grid cells in the ablation zone of coarse grids is responsible for their large melt bias.

Figure 13. Average grid cell mean melt bias over the GrIS, computed relative to the RACMO data sets using the common ice mask approach, and conditionally sampled by grid cell ice fraction provided by the common ice masks. Solid lines are the mean of the distribution with ± one standard deviation expressed by shading.

Changes in surface albedo may also explain the melt differences in the grids and dycores. The top and middle rows of Figure 14 show maps of summer surface albedo biases in the simulations, relative to the RACMO2.3p2 and the processed MODIS albedo data set of Tedesco and Alexander (2013). The model surface albedo is diagnosed from the incident and net shortwave radiative fluxes at the surface. All grids and dycores have been mapped to the $$\mathtt{f}19$$ grid in the figure, and only cells overlapping with the $$\mathtt{f}19$$ GrIS ice mask are shown. Both RACMO2.3p2 and MODIS data sets indicate that the albedo in the quasi-uniform grids is lower in the central-west and northwest region, relative to the other grids, whereas the albedo is anomalously low at the southwest and southeast coastal margins in both lat-lon and quasi-uniform grids, compared with the VR grids. All grids and dycores have a positive albedo anomaly in northeast Greenland, which may explain the depressed melt rates and unrealistic ice sheet expansion in this region when coupled to the dynamic ice sheet (Lofverstrom et al., 2020).

Figure 14. (Middle and bottom) 1979–1998 northern hemisphere summer albedo bias (fraction), computed as the difference (middle) from RACMO2.3p2, and (bottom) the MODIS data set of Tedesco and Alexander (2013). (Bottom) 1979–1998 bias in summer snow cover (fraction) relative to ERA5. All model and validation data have been mapped to the f19 $\mathtt{f}19$ grid to facilitate the comparison across the different grids.

The lower albedo values in the lat-lon and quasi-uniform grids is consistent with their larger melt rates. Since bare ice has a lower albedo than snow, lower snow fractions over the GrIS will expose larger areas of bare ice and reduce the surface albedo. The bottom panel of Figure 14 shows the summer snow fraction bias, computed relative to the ERA5 reanalysis data set. There is broad agreement with the albedo patterns. The quasi-uniform grids have low snow fractions in the central-west and north-west coastal regions, and both lat-lon and quasi-uniform grids have anomalously low snow fractions along the southeast and southwestern coasts. In contrast, the $$\mathtt{A}\mathrm{rctic}$$ grid has larger snow fractions around the south and western coasts, in better agreement with ERA5. The $$\mathtt{A}\mathrm{rctic}-\mathrm{GrIS}$$ grid has low snow fractions along the western coasts, similar to the lat-lon and quasi-uniform grids.

To the extent that the snow fraction changes are responsible for the albedo anomalies in the lat-lon and quasi-uniform grids, we seek to understand the snow cover changes. As the lat-lon and quasi-uniform grids accumulate too much precipitation in the GrIS interior, it is possible that the coastal ice margins—where seasonal melt typically occurs—receive too little snow, which may be more easily melted away to expose the underlying glacial ice. However, Figure 9 does not support a reduction in precipitation along the coastal margins, even when only plotting the solid precipitation (not shown). Instead it seems more likely that the warmer near surface environment, especially in the quasi-uniform grids, leads to anomalously low snow cover and consequent albedo changes.

Precipitation Extremes

Synoptic storms are tracked using TempestExtremes atmospheric feature detection software (Ullrich et al., 2021). As the $$\mathtt{A}\mathrm{rctic}$$ grid contains 1/4° refinement north of about 45° latitude, the storm tracker is applied to this region for the $$\mathtt{A}\mathrm{rctic}$$ and $$\mathtt{n}\mathrm{e}30\mathrm{pg}3$$ runs to identify differences in storm characteristics due to horizontal resolution.

Figure 15 shows monthly PDF (probability density function) of the precipitation rates associated with storms. The PDFs are constructed by sampling all the precipitation rates within 30° of the storm center, for each point on the storm track and for all storms. The PDFs are evaluated on an identical composite grid for all runs, and so storm statistics are not impacted by differences in output resolution. The $$\mathtt{A}\mathrm{rctic}$$ run has larger extreme precipitation rates compared to $$\mathtt{n}\mathrm{e}30\mathrm{pg}3$$ in every month, but the increase is greatest in the summer months, which coincides with the most extreme events of the year. This is primarily due to increased resolution and not the reduced physics time step; the $${\mathtt{n}\mathtt{e}\mathtt{3}\mathtt{0}\mathtt{p}\mathtt{g}\mathtt{3}}^{\ast }$$ run only marginally increases the extreme precipitation rates compared with $$\mathtt{n}\mathrm{e}30\mathrm{pg}3$$ (Figure 15).

Figure 15. PDFs of the total precipitation rate associated with tracked storms, by month, in the ne30pg3 $\mathtt{n}\mathrm{e}30\mathrm{pg}3$, ne30pg3∗ ${\mathtt{n}\mathtt{e}\mathtt{3}\mathtt{0}\mathtt{p}\mathtt{g}\mathtt{3}}^{\ast }$ and Arctic $\mathtt{A}\mathrm{rctic}$ runs, and compared with the ERA5 data set.

The extreme precipitation rates in the $$\mathtt{A}\mathrm{rctic}$$ run are closer than $$\mathtt{n}\mathrm{e}30\mathrm{pg}3$$ to the ERA5 reanalysis (Figure 15). It is difficult to know how much the extreme precipitation rates in ERA5 are constrained by data assimilation, or whether these precipitation rates are due to using a similar 1/4° model as the $$\mathtt{A}\mathrm{rctic}$$ grid. However, it is well documented that 1/4° models are more skillful at simulating extreme events (Bacmeister et al., 2013; Obrien et al., 2016). A more realistic representation of extreme precipitation events is an additional benefit of the VR grids.

Conclusions

Running CESM2.2 in an AMIP style configuration, we have evaluated six grids from two dynamical cores for their performance over the Arctic and in simulating the GrIS SMB. The 1–2° FV grids have enhanced resolution over polar regions due to their convergence of meridian lines, although a polar filter is used to prevent spurious atmospheric features from forming in these regions. SE grids comparable to the resolution of the FV grids have an isotropic grid structure where the grid resolution is similar over the entire model domain. We developed two VR grids and introduced them into CESM2.2 as part of this work. Both use the SE dycore; the $$\mathtt{A}\mathrm{rctic}$$ grid has 1/4° refinement over the broader Arctic, whereas the $$\mathtt{A}\mathrm{rctic}-\mathrm{GrIS}$$ grid is identical except for a 1/8° patch of refinement over Greenland. A third VR grid, $$\mathtt{C}\mathrm{ONUS}$$, with 1/8° refinement over the contiguous US, has also been made available in CESM2.2.

In general, the FV grids simulate colder summer temperatures over the Arctic compared with the SE grids (including the VR grids). The cloud biases in all the lat-lon and quasi-uniform resolution grids (FV and SE) are similar, in general being too cloudy over Arctic land masses. It should be emphasized that our analysis is specific to the Arctic summer because of its relevance to GrIS melt rates. An improved representation of clouds in the Arctic does not imply improved clouds at lower latitudes.

At the regional level, there is a halo of negative cloudiness bias around the ocean perimeter of Greenland on all 1–2° grids, but not the VR grids. This negative cloud bias contrasts with a positive cloud bias over the ice sheet interior. This anomaly pattern is attributed to deficient orographic precipitation in the coarser model grids. With overly smooth topography on the 1–2° grids, synoptic systems moving into Greenland are not sufficiently lifted when encountering the steep ice margins. As a result, excess precipitation falls in the GrIS interior, instead of being concentrated on the steep coastal margins as shown by observations (Pollard & Groups, 2000; Van Kampenhout et al., 2019). This results in a positive precipitation and cloud bias in the ice sheet interior, and a halo of low cloud bias about the perimeter. The agreement of different observational data products on this bias lends confidence in the attribution of causes. The VR grids compare better to the observations and show that orographic precipitation in Greenland is largely resolved when the horizontal resolution is increased sufficiently.

We integrated the primary source and sink terms of the SMB equation over the GrIS for each of the six grids. The 1°–2° lat-lon and quasi-uniform grids have large positive accumulation biases because they fail to resolve orographic precipitation. The quasi-uniform SE grids have larger accumulation biases, suggesting that the FV grids are more skillful for precipitation due to finer resolution over Greenland, and despite a polar filter. The VR grids have the most accurate accumulation rates of all the grids. The primary mass sink term of the GrIS, ice/snow melt, has similar biases; the coarse grids melt too much, with a greater bias for quasi-uniform SE grids. In general, on coarse grids, errors in the individual SMB terms are larger than the errors in the SMB itself, due to compensating errors. This observation serves as a precaution; projecting mass-loss from a glacier or ice sheet cannot be reliable if the processes representing the components of the SMB are incorrect from the start, even if the total SMB has the right magnitude.

The $$\mathtt{A}\mathrm{rctic}-\mathrm{GrIS}$$ grid has the warmest summer lower troposphere of all grids, yet it has less melting than the quasi-uniform resolution SE grids. This suggests that grid resolution (in coarse grids) is contributing to the melt biases in a way that is not obvious from the large-scale dynamics. The mechanism we propose is that coarse grids represent ablation zones using grid cells with mixed surface types, ice-covered and ice-free. The warmer ice-free patches may largely determine the mean state, leading to a warm bias over the ice-covered patches of the grid cell. This mechanism is supported by analysis of melt biases binned by grid-cell ice fraction. We leave further analysis of this hypothesis to future work.

The $$\mathtt{A}\mathrm{rctic}$$ grid substantially improves the simulated Arctic climate, including precipitation extremes and the GrIS SMB, compared to the 1°–2° lat-lon and quasi-uniform grids. The $$\mathtt{A}\mathrm{rctic}-\mathrm{GrIS}$$ grid has the most realistic cloud and precipitation fields, but the summer temperatures are too warm. The 1° FV model simulates a surprisingly realistic SMB, likely due to the relatively fine resolution of Greenland on lat-lon grids (but perhaps also because it is the most heavily tuned model configuration in CESM). In particular, a greater number of grid cells in the ablation zone reduces the influence of mixed ice-covered/ice-free grid cells that represent ablation poorly on the other (coarser) lat-lon and quasi-uniform grids.

As modeling systems move away from lat-lon grids toward quasi-uniform unstructured grids, it is worth taking stock of whether this will degrade the simulated polar climate. We have found that the 1° FV model has clear advantages over the 1° SE model for simulating the GrIS SMB. The SE dycore is still under active development (e.g., Appendix A) compared to the more mature FV dycore, and future algorithmic improvements may reduce the FV-SE GrIS SMB skill gap. However, such developments are unlikely to eliminate this skill gap entirely, because quasi-uniform unstructured grids have fewer grid cells representing high-latitude structures. Thus, the simulated GrIS SMB will likely be adversely impacted in future CESM versions, after the FV dycore is phased out. This finding will not interrupt the ongoing transition toward unstructured grids in CESM however, which is largely driven by gains in computational efficiency and grid refinement capabilities. We therefore provide the Arctic refined-meshes to the community by way of CESM2.2, providing users the option to simulate a more realistic GrIS SMB, although at a substantial computational premium relative to conventional 1°–2° grids.

We are working to develop a configuration of the $$\mathtt{A}\mathrm{rctic}$$ grid that is fully coupled with the CESM ocean and sea ice components and the Community Ice Sheet Model (CISM), to provide multi-century projections of the state of the GrIS and its contribution to sea-level rise. We have also developed a visualization of the $$\mathtt{A}\mathrm{rctic}-\mathrm{GrIS}$$ run, now available on YouTube (see link in Acknowledgments) to increase awareness of the capabilities of CESM2.2. Figure 16 shows a snapshot of this visualization, illustrating mesoscale katabatic winds descending the southeastern slopes of the GrIS. These new grids and configurations will provide new opportunities for CESM polar science, and they aim to contribute to an improved understanding of the polar environment. However, we recognize the potentially prohibitive costs for some users, and so will continue to explore different grids, parameterizations and workflows that can provide some of the same benefits of the VR grids, but at a lower cost.

Figure 16. Snapshot of the lowest model level streamlines from the Arctic−GrIS $\mathtt{A}\mathrm{rctic}-\mathrm{GrIS}$ visualization, with color shading denoting the wind magnitude.

A AppendixDetails on Spectral-Element Dynamical Core Improvements Since the CESM2.0 Release

Since the CESM2.0 release of the spectral-element dynamical core documented in Lauritzen et al. (2018) some important algorithmic improvements have been implemented and released with CESM2.2. These pertain mainly to the flow over orography that, for the spectral-element dynamical core, can lead to noise aligned with the element boundaries (Herrington et al., 2018).

A1 Reference Profiles

Significant improvement in removing noise for flow over orography can be achieved by using reference profiles for temperature and pressure [Image Omitted. See PDF] [Image Omitted. See PDF](Simmons & Jiabin, 1991) where $${T}_{1}={{\Gamma}}_{0}{T}_{\mathit{ref}}{c}_{p}^{(d)}/g\approx 192K$$, with gravity g, and standard lapse rate Γ₀ ≡ 6.5 K/km and T₀ ≡ T_ref − T₁ ≈ 97 K; T_ref = 288 K ($${c}_{p}^{(d)}$$ specific heat of dry air at constant pressure; R^(d) gas constant for dry air), and Φ_s is the surface geopotential. The reference Exner function is [Image Omitted. See PDF]where $$\kappa =\frac{{R}^{(d)}}{{c}_{p}^{(d)}}$$. The reference surface pressure p₀ = 1,000 hPa and at each model level the reference pressure p^(ref) is computed from $${p}_{s}^{(\mathrm{ref})}$$ and the standard hybrid coefficients [Image Omitted. See PDF]where A and B are the standard hybrid coefficients (using a dry-mass generalized vertical mass coordinate η). These reference profiles are subtracted from the prognostic temperature and pressure-level-thickness states before applying hyperviscosity: [Image Omitted. See PDF] [Image Omitted. See PDF]

This reduces spurious transport of temperature and mass up/down-slope due to the hyperviscosity operator.

A2 Rewriting the Pressure Gradient Force (PGF)

In the CESM2.0 the following (standard) form of the pressure gradient term was used: [Image Omitted. See PDF]where Φ is geopotential and $$\rho =\frac{{R}^{(d)}{T}_{v}}{p}$$ is density (for details see Lauritzen et al. [2018]). To alleviate noise for flow over orography, we switched to an Exner pressure formulation following Taylor et al. (2020), which uses that Equation A7 can be written in terms of the Exner pressure [Image Omitted. See PDF]where the Exner pressure is [Image Omitted. See PDF]and virtual temperature is [Image Omitted. See PDF]where m^(ℓ) is dry mixing ratio of component of moist air ℓ; $${\mathcal{L}}_{\mathit{all}}$$ is the set of all components of moist air and, in particular, “wv” is water vapor.

The derivation showing that Equations A7 and A8 are equivalent is given here: [Image Omitted. See PDF]where virtual potential temperature is T_v/Π. Using the reference states from Simmons and Jiabin (1991), [Image Omitted. See PDF] [Image Omitted. See PDF]we can define a geopotential as a function of Exner pressure [Image Omitted. See PDF]

This “balanced” geopotential obeys [Image Omitted. See PDF]for any Exner pressure. Subtracting this “reference” profile from the PGF yields [Image Omitted. See PDF]

In the continuum, the two formulations (left and right-hand side of Equation A15) are identical. But under discretization, the second formulation can have much less truncation error.

A3 Results

One year averages of vertical pressure velocity at 500 hPa (OMEGA500) have been found to be a useful quantity to detect spurious up or down-drafts induced by steep orography (Figure A1). While the true solution is not known, strong vertical velocities aligned with element edges that are not found in the CAM-FV reference solution (Figure A1a) are likely not physical (spurious). The older CESM2.0 version of SE (Figure A1d) using the “traditional” discretization of the PGF, (A15), exhibits significant spurious noise patters around steep orography compared to CAM-FV (e.g., around Himalayas and Andes). This is strongly alleviated by switching to the Exner formulation of the PGF (A8; Figure A1c). By also subtracting reference profiles from pressure-level thickness and temperature, Equation A5 and A6 respectively, reduces strong up-down drafts further (Figure A1d). Switching to the CAM-SE-CSLAM version where physics tendencies are computed on an quasi-equal area physics grid and using the CSLAM transport scheme, marginal improvements are observed in terms of a smoother vertical velocity field (Figures A1e and A1f). The configuration shown in Figure A1d is used for the simulations shown in the main text of this paper.

Figure A1. One year averages of vertical pressure velocity at 500 hPa (OMEGA500) $(\mathtt{O}\mathrm{MEGA}500)$ using (a) CAM-FV (Finite-Volume dynamical core) and (b–f) various versions of the spectral-element (SE) dynamical core at approximately 1° horizontal resolution and using 32 levels. Plot (b) is equivalent to the CESM2.0 version of the SE dynamical core using the “traditional”/“old” discretization of the pressure-gradient force (PGF). Plot (c) is equivalent to configuration (b) but using the Exner form of the PGF. Plot (d) is the same as configuration (c) but also subtracting reference profiles from pressure and temperature before applying hyperviscosity operators (which is equivalent to the CESM2.2 version of SE in terms of the dynamical core). Plots (e and f) are equivalent to (c and d), respectively, by using the SE-CSLAM (ne30pg3) $(\mathtt{n}\mathrm{e}30\mathrm{pg}3)$ version of the SE dynamical core (i.e., separate quasi-uniform physics grid and CSLAM transport scheme).

It is interesting to note that the noise issues and algorithmic remedies found in the real-world simulations discussed above, can be investigated by replacing all of physics with a modified version of the Held-Suarez forcing (Held & Suarez, 1994). The original formulation of the Held-Suarez idealized test case used a flat Earth (Φ_s = 0) and a dry atmosphere. By simply adding the surface topography used in “real-world” simulations and removing the temperature relaxation in the lower part of domain (σ > 0.7; see Held and Suarez [1994] for details), surprisingly realistic vertical velocity fields (in terms of structure) result (see Figure A2). Since this was a very useful development tool it is shared in this manuscript.

Figure A2. Same as Figure A1 but using modified Held-Suarez forcing and the average is over 18 months (excl. spin-up).

Acknowledgments

This material is based upon work supported by the National Center for Atmospheric Research (NCAR), which is a major facility sponsored by the NSF under Cooperative Agreement 1852977. Computing and data storage resources, including the Cheyenne supercomputer (Computational and Information Systems Laboratory, 2017), were provided by the Computational and Information Systems Laboratory (CISL) at NCAR. A. Herrington was funded through the NCAR Advanced Study Program Postdoctoral Fellowship. A. Herrington thanks Matt Rehme (NCAR/CISL) for his role in generating the $$\mathtt{A}\mathrm{rctic}-\mathrm{GrIS}$$ visualization available on youtube (https://www.youtube.com/watch?v=YwHgqDu75s8&t=4s&abchannel=NCARVisLab). To download CESM2.2, go to https://github.com/ESCOMP/CESM and check-out the latest $\tt cesm2_2$ tag.

Data Availability Statement

The data presented in main part of this manuscript is available at https://github.com/adamrher/2020-arcticgrids. The source code and data for the Appendix is available at https://github.com/PeterHjortLauritzen/CAM/tree/topo-mods.