The fully coupled regionally refined model of E3SM version 2: overview of the atmosphere, land, and river results

This paper provides an overview of the United States (US) Department of Energy's (DOE's) Energy Exascale Earth System Model version 2 (E3SMv2) fully coupled regionally refined model (RRM) and documents the overall atmosphere, land, and river results from the Coupled Model Intercomparison Project 6 (CMIP6) DECK (Diagnosis, Evaluation, and Characterization of Klima) and historical simulations – a first-of-its-kind set of climate production simulations using RRM. The North American (NA) RRM (NARRM) is developed as the high-resolution configuration of E3SMv2 with the primary goal of more explicitly addressing DOE's mission needs regarding impacts to the US energy sector facing Earth system changes. The NARRM features finer horizontal resolution grids centered over NA, consisting of 25→100 km atmosphere and land, a 0.125^∘ river-routing model, and 14→60 km ocean and sea ice. By design, the computational cost of NARRM is $\sim \mathrm{3}\times$ of the uniform low-resolution (LR) model at 100 km but only ∼ 10 %–20 % of a globally uniform high-resolution model at 25 km.

A novel hybrid time step strategy for the atmosphere is key for NARRM to achieve improved climate simulation fidelity within the high-resolution patch without sacrificing the overall global performance. The global climate, including climatology, time series, sensitivity, and feedback, is confirmed to be largely identical between NARRM and LR as quantified with typical climate metrics. Over the refined NA area, NARRM is generally superior to LR, including for precipitation and clouds over the contiguous US (CONUS), summertime marine stratocumulus clouds off the coast of California, liquid and ice phase clouds near the North Pole region, extratropical cyclones, and spatial variability in land hydrological processes. The improvements over land are related to the better-resolved topography in NARRM, whereas those over ocean are attributable to the improved air–sea interactions with finer grids for both atmosphere and ocean and sea ice. Some features appear insensitive to the resolution change analyzed here, for instance the diurnal propagation of organized mesoscale convective systems over CONUS and the warm-season land–atmosphere coupling at the southern Great Plains. In summary, our study presents a realistically efficient approach to leverage the fully coupled RRM framework for a standard Earth system model release and high-resolution climate production simulations.

1 Introduction

Global Earth system models (ESMs) are fundamental tools for understanding the past evolution of the climate system and projecting future climate changes under various anthropogenic scenarios. High horizontal resolution simulations on climate scales have been recognized as one of the increasingly important directions of ESM development in recent years (Demory et al., 2014; Haarsma et al., 2016). Compared to low-resolution models, high-resolution models show superior fidelity in representing both the large-scale circulation (e.g., meridional ocean heat transport) (Griffies et al., 2015) and small-scale processes (e.g., clouds and streamflow) (Haarsma et al., 2016, and references therein). More importantly, simulations with enhanced horizontal resolution exhibit improved skills in capturing regional climate change signals and facilitating process-level studies, which provide a crucial basis for assessing the impacts of climate extremes with augmented societal implications. However, fine-resolution and multi-century simulations (with ensembles) are competing requirements for climate experiments due to limited computational and human resources. This conflict will likely continue to challenge the climate modeling community, as evidenced by the fact that more than 3 times (72 vs. 23) as many model sources (including different versions of the same model) have published simulations at 100 km than at 25 km nominal resolutions in the current Coupled Model Intercomparison Project 6 (CMIP6) archive (https://esgf-node.llnl.gov/search/cmip6/, last access: 18 August 2022). This suggests that despite the commonly recognized benefits, not many modeling centers can afford to pursue routine high-resolution climate simulations.

The Energy Exascale Earth System Model (E3SM) project (Leung et al., 2020) is supported by the US Department of Energy (DOE) with a primary goal of improving actionable predictions of Earth system variability and change by leveraging advanced DOE computational resources. Scientifically, E3SM development is motivated by modeling requirements in three overarching fields (i.e., water cycle, biogeochemistry, and cryosphere) to address the most critical DOE mission-related questions, such as water availability, wildfires, heat waves, and sea level rise, which all pose challenges to the energy sector with climate change. High-resolution simulations are clearly more desirable to achieve these E3SM objectives since these processes have high spatiotemporal variability. However, uniformly increasing the grid size for climate production simulations is not an easy task even with DOE's world class high-performance computing power. For example, the 25 km simulation is at least 32 times (16× more grid cells, 2× smaller physics time step, and 4× smaller dynamical core time step) more expensive than the 100 km version with the E3SM version 1 (E3SMv1) model (Caldwell et al., 2019), making high-resolution models much more computationally expensive not only to run but also to tune for skillful simulations. With these demands and limitations, a multiscale approach is an attractive avenue for global ESMs to deliver high-resolution production simulations over target areas at a more economical cost.

The multiresolution method (Ringler et al., 2008; Leung et al., 2013), also known as regionally refined model (RRM) or variable-resolution (VR) model, was proposed to alleviate the computational burden of global ESMs by refining a fraction of the globe with higher resolution while keeping (without coarsening) the remaining area at lower resolution. The RRM method is a general tool for all major ESM components, such as atmosphere, land, ocean, and sea ice. With a careful design of the RRM mesh, the high-resolution grids can better represent fine-scale processes over an area of interest at a typical cost of only ∼ 10 %–20 % of a comparable globally uniform high-resolution configuration. Compared to regional or nested climate models, global RRMs by design minimize the impacts from the lack of a two-way dynamical feedback between the refined area and the outside domain.

Recently, an increasing number of studies have successfully applied the RRM technique in global ESMs to tackle a wide range of climate research themes from climatological statistics of idealized aquaplanet (Zarzycki et al., 2014) and mean climate state of more realistic simulations (Sakaguchi et al., 2015, 2016; Gettelman et al., 2018; Tang et al., 2019) to complex terrain climate (Wu et al., 2017; Rhoades et al., 2018c; Rahimi et al., 2019; Bambach et al., 2022) and climate extremes (Huang and Ullrich, 2017; Rhoades et al., 2020a, b; Zarzycki et al., 2021; Reed et al., 2022; Xu et al., 2022). Others leveraged RRM to study specific aspects of climate, such as tropical cyclones (Zarzycki and Jablonowski, 2014, 2015; Hazelton et al., 2018), marine stratocumulus (Bogenschutz et al., 2023), snowpack (Rhoades et al., 2016, 2017), surface energy flux (Burakowski et al., 2019), Greenland surface mass balance (van Kampenhout et al., 2019), irrigation impacts on regional climate (Huang and Ullrich, 2016), and land use and land cover change influence on land–atmosphere coupling and precipitation (Devanand et al., 2020). Lately, the RRM resolution has been pushed to a new limit for watershed-scale hydrology analysis (Xu and Di Vittorio, 2021) and cloud-resolving scale climate simulation (Liu et al., 2022).

The RRM high-resolution results are robust for most places except the Intertropical Convergence Zone (Rauscher et al., 2013; Zarzycki et al., 2014), covering almost all typical climate regimes such as the contiguous US (CONUS) (Gettelman et al., 2018; Tang et al., 2019), the western (Rhoades et al., 2016; Huang et al., 2016; Huang and Ullrich, 2017; Rhoades et al., 2018c) and eastern US (Liu et al., 2022), South America (Sakaguchi et al., 2015, 2016; Bambach et al., 2022), Asia (Sakaguchi et al., 2016), East Asia (Liang et al., 2021), eastern China (Xu et al., 2021), the Tibetan Plateau (Rahimi et al., 2019), the Maritime Continent (Harris and Lin, 2014), Atlantic basin (Zarzycki et al., 2015), the southeastern Pacific (Bogenschutz et al., 2023), Greenland (van Kampenhout et al., 2019), and the Arctic (Veneziani et al., 2022). Furthermore, the RRM capability in representing the general high-resolution climate seems generally acceptable for different models, including the Variable-Resolution Community Earth System Model (VR-CESM) (e.g., Gettelman et al., 2018), the E3SMv1 atmospheric model (EAMv1) (Tang et al., 2019), the Model for Prediction Across Scales-Atmosphere (MPAS-A) (Hagos et al., 2013; Sakaguchi et al., 2015, 2016; Liang et al., 2021), the Geophysical Fluid Dynamics Laboratory finite-volume dynamical core on the cubed-sphere grid (Harris and Lin, 2013, 2014), and the ICOsahedral Non-hydrostatic Earth System Model (ICON-ESM) (Jungclaus et al., 2022).

All of the aforementioned studies utilize RRMs for Atmospheric Model Intercomparison Project (AMIP-type) (Gates et al., 1999) simulations. Although these studies provide valuable experience and important knowledge about RRMs, modeling centers still face the question of how to transform such AMIP-type RRM achievements from individual scientific studies emphasizing specific climate aspects to a standard global ESM release version aiming at a much broader and general scope. At a minimum, the criteria of reasonable global climate should be satisfied for the RRM to be widely adopted for global ESM releases. Most previous AMIP-type RRM studies focus on the regional results within the refined grids without paying much attention to the outside domain. While this might be acceptable for targeted studies, one cannot release a global model without reasonable global results since such a model is expected to address the challenge of a long (multi-century) spinup and demonstrate top-of-atmosphere (TOA) radiative balance in pre-industrial fully coupled simulations. In addition, some physics parameterizations (e.g., deep convection) suffer from poor scale awareness and hence require retuning as the model resolution increases (e.g., Xie et al., 2018). This implies significant model calibration efforts that modeling centers have to seriously consider when planning on releasing the RRM besides the low-resolution model. Furthermore, based on our EAMv1 RRM experience, retuning does not guarantee improved global climate performance. In the present study, building upon the EAMv1 RRM (atmosphere and land area of 25→100 km horizontal resolution with the 25 km mesh over the CONUS) (Tang et al., 2019) plus the E3SMv2 lower resolution configuration (Golaz et al., 2022), we extend the RRM configuration to ocean and sea ice (see grids in Fig. A1) as a fully coupled RRM with fine meshes centered over North America (NA). We propose an innovative RRM strategy (see details in Sect. 2.1) to meet the criteria above with a minimal retuning effort and for the first time to deliver production climate simulations using a fully coupled RRM.

This paper focuses on the atmosphere, land, and river components of the E3SMv2 North American RRM (NARRM), while a companion paper (Luke P. Van Roekel, personal communication, 2023) overviews the NARRM ocean and sea ice. This paper is organized as follows. Section 2 describes the NARRM model, our hybrid time step strategy for the atmospheric component, and key tools and tests used to create its atmospheric configuration. Section 3 summarizes the simulations performed in the present study and reports on the computational cost of the NARRM historical simulation relative to its lower-resolution (LR) counterpart. Analyses of model results start at the global scale in Sect. 4 and then shift to the high-resolution NA region in Sect. 5 for atmosphere, land and river, and land–atmosphere interactions. Conclusions and discussions are presented in Sect. 6.

2 Model description

Except for the mesh and mesh-related settings, E3SMv2 LR and NARRM essentially have the same atmosphere, land, and river components. They are upgraded from E3SMv1 and briefly described here. In the E3SMv2 atmosphere model (EAMv2), the dynamical core uses the High-Order Method Modeling Environment (HOMME) package (Dennis et al., 2005, 2011; Evans et al., 2013) on the spectral element grid (Taylor and Fournier, 2010). HOMME has been updated to use a potential temperature formulation of the equations with a more accurate pressure gradient (Taylor et al., 2020; Herrington et al., 2022) and a new interpolation semi-Lagrangian scheme (Islet) for passive tracer transport (Bradley et al., 2022). The physics operates on a separate finite-volume grid (Hannah et al., 2021), which has four-ninths as many columns as the corresponding spectral element grid (see Table 1) and hence runs about 2× faster than it would on the spectral element grid. The physics parameterization updates include the Cloud Layers Unified By Binormals scheme (Golaz et al., 2002; Larson, 2017) for subgrid turbulent transport and cloud macrophysics, the Zhang–McFarlane (ZM) deep convection scheme (Zhang and McFarlane, 1995) with a new trigger method (Xie et al., 2019), gravity wave parameterizations following Richter et al. (2010) with additional modifications (Beres et al., 2004; Richter et al., 2019), the O3v2 package (Tang et al., 2021) for the prognostic stratospheric ozone, and the four-mode version of Modal Aerosol Module (MAM4) (Liu et al., 2016; Wang et al., 2020) with an updated treatment of dust aerosol (Feng et al., 2022). The same set of EAM physics parameters is used in the LR and NARRM simulations analyzed here. The LR grid is a quasi-uniform 1^∘ cubed sphere grid with an average grid spacing of ∼100 km. The NARRM grid has an average grid spacing of ∼ 25 km over North America, transitioning to match the ∼ 100 km cubed-sphere grid over the rest of the globe. All simulations, except the idealized baroclinic wave simulations described later, utilize E3SM's standard 72 vertical levels (L72).

The E3SMv2 land model (ELMv2) runs on the same grid as the atmospheric physics. ELMv2 upgrades the prescribed vegetation distribution for better consistency between land use and changes in plant functional types across platforms and adopts the new shortwave radiation model SNICAR-AD (Dang et al., 2019) for snow and ice. The land use harmonization version 2f data (LUH2; https://luh.umd.edu/data.shtml, last access: 3 July 2023) (Hurtt et al., 2020) are converted into E3SMv2 plant functional types with an updated version of the land use translator (Di Vittorio et al., 2014). The trajectory of land cover change has also been improved through better tracking of previous land use change. The E3SMv2 river-routing model (Model for Scale Adaptive River Transport, MOSARTv2) utilizes the regular lat–long grid (0.5^∘ for LR and 0.125^∘ for NARRM). MOSARTv2 uses the kinematic wave method to route the runoff from ELM into the ocean model via an eight-direction-based river network (Li et al., 2013). More details about the E3SMv2 model are documented by Golaz et al. (2022).

2.1 EAM hybrid time step strategy for RRM production simulations

In previous RRM studies, including the EAMv1 CONUS RRM (Tang et al., 2019), the atmospheric physics time step is often chosen to be shorter than that of the globally uniform low-resolution model to match the highest-resolution grids in the RRM. However, such treatment faces the challenge of satisfying the criteria above for the purpose of global climate production simulations. Mainly because the ZM deep convection scheme and other cloud parameterizations used by EAM are by design not scale-aware (Xie et al., 2018), if the EAM in NARRM used a shorter physics time step than LR while keeping other physics parameters unchanged, the NARRM results on the unrefined portion of the mesh (covering a larger area than the refined portion) would not match the quality of the LR results and thus undermine the NARRM global performance. Furthermore, even if NARRM used the retuned high-resolution physics parameters along with the shorter physics time step, we would still have degraded global simulation quality over the LR model (see Fig. A2 for the EAMv1 results). With all these considerations, in the present study when employing RRM for climate production campaigns, we opt for a hybrid time step strategy in EAM, which is a combination of an LR physics time step and the high-resolution dynamics time steps (see Table 1). In this way, NARRM retains much of the LR global climate characteristics with possible improvements at the refined area benefiting from the high-resolution dynamics. Moreover, this approach simplifies the RRM development as it naturally avoids further tuning the RRM beyond what was done for LR. This choice also ensures that the physics behaves as similarly as possible between the LR and RRM simulations to facilitate direct comparisons of their climates.

It is worthwhile noting that the hybrid time step strategy is a practical choice before the scale-aware cloud parameterization becomes available. With the coarsened physics time step, NARRM cannot take full advantage of resolved processes (e.g., updrafts) at 25 km because the dynamics at 25 km explicitly resolve greater vertical velocities relative to those at 100 km and hence have faster dynamical timescales, which require the correspondingly shortened physics time step to match the faster-evolving instability. The time-truncation errors of the hybrid time step method are large at 25 km as quantified by a moist bubble test (Herrington et al., 2019).

2.2 EAM running on unstructured meshes

In EAM, the underlying grid is always treated as fully unstructured. EAM can run on any grid that represents a tiling of the sphere with quadrilateral elements. For quasi-uniform grids, EAM relies on cubed-sphere grids since these grids are simple to construct. RRM grids are constructed by external tools as described below. Internally, the code treats all these grids identically, the only difference being the various resolution-dependent parameters. For the dynamical core, these parameters consist of the many time steps in the model (given in Table 1) and the hyperviscosity coefficient. The dynamical core time steps are chosen to ensure stability of the model. For RRM grids, these time steps are chosen to match those that would be used in a global model with the same resolution as the highest resolution contained within the RRM. For the NARRM grid used here, which includes refinement down to 25 km, we use the same time steps as would be used by a global 25 km configuration of EAM.

Table 1 Column numbers and time steps of the atmosphere component used in LR and NARRM simulations.

For hyperviscosity, EAM relies on a resolution-aware tensor hyperviscosity formulation (Guba et al., 2014) applied on each model surface. The tensor coefficients vary spatially based on the two length scales of each spectral element (derived from the eigenvalues of the reference element map). This operator has a built in scaling of Δx³ with strength controlled by a coefficient ν with units of per second. The tensor is designed to have the proper directional resolution dependence for highly distorted elements, while matching the traditional constant-coefficient hyperviscosity on square elements. In EAMv2, we use the tensor hyperviscosity operator with $\mathit{\nu }=\mathrm{3.4}\times {\mathrm{10}}^{-\mathrm{8}}$ s⁻¹ for all grids (cubed-sphere and RRM) and at all resolutions. The only exception is the LR 1^∘ cubed-sphere grid, where for continuity with older simulations we continue to use the constant-coefficient hyperviscosity operator with $\mathit{\mu }=\mathrm{1}\times {\mathrm{10}}^{\mathrm{15}}$ m⁴ s⁻¹. For a uniform degree p spectral element grid with square elements, the tensor operator with coefficient ν is identical to a constant coefficient hyperviscosity operator with coefficient $\mathit{\mu }={\left(\frac{p}{\mathrm{2}}\right)}^{\mathrm{3}}R\mathrm{\Delta }{x}^{\mathrm{3}}\mathit{\nu }$, where Δx is the element edge length divided by p and R is the radius of the sphere. In EAM, we always use p=3.

2.3 Key tools for the RRM configuration

A number of tools have been developed to streamline the workflow for EAM and ELM simulations on RRM grids. These are described as follows, in the approximate order they are employed.

• The Spherical Quadrilateral Grid Generator (SQuadGen). Generation of the atmosphere–land mesh is performed using SQuadGen (Ullrich, 2022; Guba et al., 2014). This tool translates a monochrome PNG image, which denotes the desired level of grid refinement on an equirectangular projection, to a mesh of refined quadrilaterals based on a cubed sphere. The use of quadrilaterals is by necessity for compatibility with the spectral element dynamical core. Transition regions are managed using “paving”, that is, using predefined patterns of quadrilaterals which enable transition between coarse-resolution and fine-resolution regions. Smoothing of the grid is performed via spring dynamics. The spectral elements of the NARRM grid produced with this procedure are shown in Fig. 1. In the spectral element method, each field is represented by polynomials up to degree 3 within each element. The resolution represented by each element (its average length divided by 3) is shown in Fig. 2.

• TempestRemap. The TempestRemap package (Ullrich and Taylor, 2015; Ullrich et al., 2016) is used to generate conservative, consistent, and monotone linear maps between fields stored as volume averages (i.e., updated using the finite-volume methods) and fields stored as spectral elements (i.e., as coefficients of a set of basis functions). The generated maps require the construction of an “overlap mesh”, which is the union of the source and target face; the generation of an approximate map; and subsequent projection of the approximate map onto the linear space of conservative, consistent, and (optionally) monotone maps.

• Topography generation. To generate topography and associated surface roughness fields on the NARRM grid, we rely on the tool chain described in Lauritzen et al. (2015) combined with a topography smoothing tool included with HOMME. The use of HOMME's topography smoothing tool ensures that the smoothing is done with the same discrete Laplace operator used internally in the dynamical core.

• NetCDF Operators (NCOs). NCOs consist of a number of command-line tools that enable manipulation of netCDF files (Zender, 2008). The tools include variable extraction, remapping, and spatial and temporal averaging. Provenance information is preserved within the netCDF files to enable scientific reproducibility.

Figure 1 North American RRM (NARRM) grids for the atmosphere dynamical core shown in (a) a cylindrical equidistant projection and (b) an orthographic projection.

2.4 Idealized test

Before running long coupled NARRM simulations, we first evaluate the dynamical core settings for the NARRM grid using a baroclinic instability test case. This test case establishes that the dynamical core behaves as expected in an idealized setting: the time steps are stable, the model can capture high-resolution features in the high-resolution region, and the presence of the high-resolution and mesh transition regions does not negatively impact the large-scale behavior. For this evaluation, we use an extension of the dry baroclinic wave test case by Ullrich et al. (2014) with two idealized, analytically prescribed mountains (Hughes and Jablonowski, 2023). The latter now serve as the trigger for baroclinic instability. The addition of the two mountains generates a flow with more energy at smaller scales as compared to the original test case, especially downstream of the mountains, making this an attractive test case for studying the impacts of resolution.

For this test case, we run simulations with three different horizontal grids, LR, NARRM, and high-resolution (HR), and 30 hybrid vertical levels (L30), which are specified in Appendix B of Reed and Jablonowski (2012). The LR and NARRM grids are as described above, and we add an HR grid. The HR grid is a global 0.25^∘ grid which matches the high-resolution region of the NARRM grid. All idealized runs use the same settings as in the full model (except L30 instead of L72), with HR and NARRM using identical time steps since they both contain regions of 0.25^∘ resolution. All simulations utilize the EAMv2 tensor hyperviscosity tuning with $\mathit{\nu }=\mathrm{3.4}\times {\mathrm{10}}^{-\mathrm{8}}$ and Δx³ resolution scaling.

The test case is fully described in Hughes and Jablonowski (2023). We use the dry configuration and make one modification to the locations of the mountains. In particular, the center locations of the mountains are shifted longitudinally by 144^∘ to the east in order to place the two mountains within the NARRM's high-resolution region. The new center locations are therefore 144 and 76^∘ W. The peak height of the mountain ranges is 2000 m. Figure 2 illustrates the size and location of the mountains and the NARRM mesh resolution, while Fig. 3 shows the surface pressure at day 6 computed on the same mesh. The latter highlights the topographically generated baroclinic instability in the Northern Hemisphere.

Figure 2 Contour lines of the topographic height with a peak amplitude of 2000 m overlaid on a map of the NARRM grid resolution (square root of element area). The resolution is ∼ 25 km over North America (shown in yellow), transitioning to ∼ 100 km over the rest of the globe (dark blue). The two mountains are mostly contained within the high-resolution region. In the low-resolution region, the faint outline of an inscribed cube shows the slight non-uniformness of the 1^∘ cubed-sphere grid used in that region.

Figure 3 Contours of the surface pressure at day 6 showing the topographically triggered baroclinic instability in the Northern Hemisphere as computed on the NARRM grid. The instability has yet to be triggered in the Southern Hemisphere. The mountain height contours are overlaid. The colors saturate over the mountain ranges with minimum surface pressure values around 750–780 hPa (not shown).

Figure 4 shows contour lines of the 750 hPa temperature field after 6 d on all three grids. The plots are zoomed in over the region with the most activity shown in Fig. 3. We first compare the field in the NARRM's high-resolution region with the HR result and note the remarkable agreement between the two solutions (black contour lines) in the high-resolution region (yellow color). The presence of high resolution in the NARRM simulation allows the model to capture features in that region with finer scales than can be captured by the LR simulation (as expected). Further downstream from the mountains at the right edge of the Fig. 4b, the NARRM resolution has transitioned to match the LR resolution (blue color), and the scales captured by the NARRM solution are no longer as fine as they are in the HR solution. They are somewhat dissipated and fall between the LR and HR results. Thus, the presence of the high-resolution region in the NARRM grid improves some aspects of the solution in the low-resolution region. Finally, we note that there are no visible artifacts from the distorted elements in the mesh transition region. Examination of other fields, such as vorticity (not shown), demonstrate similar results.

Figure 4 Contour lines of the 750 hPa temperature field on day 6 with contour intervals of 5 ^∘C. The temperature contours are overlaid on a map colorized by grid resolution. The data is plotted over a subset of the globe containing the mountains and most of the downstream region affected by the baroclinic instability. Results are shown from the LR grid (a), NARRM (b), and HR grid (c). The NARRM grid shows the transition from high resolution (yellow, ∼ 25 km) to low resolution (blue, ∼ 100 km).

3 Simulations and computational cost

We perform a set of NARRM production simulations parallel to the LR version documented by Golaz et al. (2022) and following the same CMIP6 specifications. The LR and NARRM production simulations analyzed in the present study are summarized in Table 2. These simulations consist of the CMIP6 Diagnosis, Evaluation, and Characterization of Klima (DECK) and historical simulations (Eyring et al., 2016), i.e., one pre-industrial control (piControl, 500 years), two idealized CO₂ runs (1pctCO2 and abrupt-4xCO2, each 150 years), a five-member historical ensemble (historical_N, 1850–2014), and a three-member Atmospheric Model Intercomparison Project (amip) type ensemble (amip_N, 1870–2014). Initial conditions are taken from 1 January of different years of piControl, as indicated in Table 2 for 1pctCO2, abrupt-4xCO2, and historical_N simulations. The amip_N simulations are initialized from the 1870 condition of corresponding historical_N simulations.

In order to estimate the effective radiative forcing of anthropogenic aerosols in LR and NARRM configurations, we perform pairs of nudged simulations with prescribed emissions of aerosols and their precursors for the present-day (PD, year 2010) and pre-industrial (PI, year 1850) values, which are taken from the CMIP6 emission data. Table 3 lists the nudged simulations used to assess the effective radiative forcing of anthropogenic aerosols. Horizontal winds in LR and NARRM are nudged towards wind fields from their respective baseline simulations, with a relaxation timescale of 6 h. These nudged simulations are 15 months long, with the first 3 months discarded as spinup. Previous studies have shown that nudging the horizontal winds can help constrain the large-scale circulation in the model (Zhang et al., 2014; Sun et al., 2019; Tang et al., 2019), meaning that the anthropogenic aerosol effects can be determined with relatively short simulations (K. Zhang et al., 2022; S. Zhang et al., 2022).

Table 2 Summary of E3SMv2 LR (Golaz et al., 2022) and NARRM production simulations analyzed in this study. Numbers in parentheses indicate the simulation year numbers.

Table 3 Nudged LR and NARRM atmospheric model simulations used in this study. All simulations are performed with prescribed sea surface temperature (SST) and sea ice concentration for year 2010. Nudging data are 6-hourly model output saved from the LR and NARRM free-running simulations (middle column). Due to the model instability problem with nudging application in RRM (with a relatively long time step), we use an alternative physics–dynamics coupling approach (see option “se_ftype = 1” in Sect. 3.1 of Zhang et al., 2018) for the NARRM nudged simulations. We find the impact of using different physics–dynamics coupling approaches on the global mean effective aerosol forcing estimate in LR to be small (difference <0.05 W m⁻²).

3.1 Computational performance

A sequence of performance benchmark simulations were run on the Argonne National Laboratory Chrysalis cluster. Chrysalis has 512 compute nodes. Each node has two AMD Epyc 7532 “Rome” 2.4 GHz processors. Each processor has 32 cores, for a total of 64 cores per node. Each node has 256 GB 16-channel DDR4 3200 MHz memory. The interconnect hardware is Mellanox HDR200 InfiniBand and uses the fat tree topology. The model code was compiled using Intel release 20200925 with GCC version 8.5.0 compatibility and run using OpenMPI 4.1.3 provided in the Mellanox HPC-X Software Toolkit.

The simulations are run with one MPI process per core and no OpenMP threading. Throughput values are computed using the maximum wall-clock time (minimum throughput) over all message passing interface (MPI) processes; model initialization time is excluded. A throughput data point corresponds to one simulation run for 90 d. The input/output (I/O) configuration is identical to production simulations. At the end of 90 d, a restart file is written.

Figure 5 Performance of the LR and NARRM historical simulations. (a) Throughput vs. number of computer nodes. Each data point is annotated with its throughput in simulated years per day (SYPD) and computer resource configuration name. The dashed gray line shows the perfect-scaling slope. (b) Computational resource plots for the L process layouts. Each component has one rectangle. A rectangle has the area given by the product of normalized wall-clock time and number of cores, with the NARRM total time normalized to 1.0.

Figure 6 Performance of the atmosphere (Atm.) and ocean components of the NARRM historical simulation. Solid lines show measured performance. Dashed lines show the performance predicted by a simple model that uses the LR simulation with the XS process layout for input data; see the text for a description of the performance model.

Figure 5 summarizes the performance of the LR and NARRM historical_N simulations for several node counts and corresponding process layouts with names T (NARRM only), XS, S, M, and L. Note that while the layout names are shared among models, the specific layout associated with a name differs among models. Each simulation's data point is annotated with its throughput in simulated years per day (SYPD) and process layout name. The highest throughput of the LR simulations is 39.81 SYPD. In Golaz et al. (2022, Fig. 2), the highest throughput is 41.89 for the same node count; historical_N simulations have additional forcings to compute relative to the piControl simulation used in Golaz et al. (2022, Fig. 2). The LR throughput falls off from the perfect scaling slope faster than the NARRM throughput because the LR simulation has less work per node. For the L process layouts, accounting for 105 vs. 100 nodes, the throughput factor difference is 3.14.

Figure 5b shows the wall-clock-time–resource product for each component for the L layouts. A rectangle's width is proportional to the number of cores the component uses; its height is proportional to the wall-clock time to simulate a fixed simulation period, with the time normalized so that the NARRM simulation has a total time of 1.0. The atmosphere (ATM), sea ice (ICE), coupler (CPL), land (LND), and river runoff (ROF; LND and ROF are too small to label) components run on one set of nodes, while the ocean (OCN) component runs on another set. An unfilled rectangle having “LR” or “NARRM” at the top-right corner shows the total product. Because there is no global communication barrier between components run in sequence, the time value of each component is approximate, and thus the filled rectangles do not sum to the total time.

We can understand the NARRM component-level performance as a function of spatial and temporal discretization parameters and one LR simulation to calibrate throughput. The LR calibration simulation should reflect that RRM simulations have a large amount of work per node; thus, we use the LR simulation run with the XS process layout, the left-most LR point in Fig. 5a. We focus on the two most expensive components, the atmosphere and ocean. We start with the ocean, whose performance is simpler to model. For simplicity, we write the formulas in terms of wall-clock time (w.c.t.) for a fixed simulation length, e.g., 90 d. The input measured datum is the top-level ocean component (ocn) wall-clock time in the LR simulation run with the XS process layout, $\mathrm{w}.\mathrm{c}.\mathrm{t}{.}_{\mathrm{LR}}^{\mathrm{ocn}}$. The input parameters are the number of computer cores (n_core) used in the LR (XS) and RRM (variable) simulations, the number of cells (n_cell) in each grid, and the time steps (Δt) in each simulation. For a fixed simulation length, the predicted ocean component RRM performance is then

$$\begin{array}{}\text{(1)}& \begin{array}{rl}\mathrm{w}.\mathrm{c}.\mathrm{t}{.}_{\mathrm{RRM}}^{\mathrm{ocn}}& ={\displaystyle \frac{(n_{\mathrm{core}}{)}_{\mathrm{LR}}^{\mathrm{ocn}}}{(n_{\mathrm{core}}{)}_{\mathrm{RRM}}^{\mathrm{ocn}}}}\cdot {\displaystyle \frac{(n_{\mathrm{cell}}{)}_{\mathrm{RRM}}^{\mathrm{ocn}}}{(n_{\mathrm{cell}}{)}_{\mathrm{LR}}^{\mathrm{ocn}}}}\cdot {\displaystyle \frac{(\mathrm{\Delta }t{)}_{\mathrm{LR}}^{\mathrm{ocn}}}{(\mathrm{\Delta }t{)}_{\mathrm{RRM}}^{\mathrm{ocn}}}}\\ & \cdot \mathrm{w}.\mathrm{c}.\mathrm{t}{.}_{\mathrm{LR}}^{\mathrm{ocn}}.\end{array}\end{array}$$

The performance model for the atmosphere is more complicated because it has two important time steps, one each for the dynamical core (dynamics) and the column parameterizations (physics). Thus, the factor accounting for model time steps is broken into two terms, one each for the physics and dynamics. The predicted atmosphere component RRM performance is then

$$\begin{array}{}\text{(2)}& \begin{array}{rl}\mathrm{w}.\mathrm{c}.\mathrm{t}{.}_{\mathrm{RRM}}^{\mathrm{atm}}& ={\displaystyle \frac{(n_{\mathrm{core}}{)}_{\mathrm{LR}}^{\mathrm{atm}}}{(n_{\mathrm{core}}{)}_{\mathrm{RRM}}^{\mathrm{atm}}}}\cdot {\displaystyle \frac{(n_{\mathrm{cell}}{)}_{\mathrm{RRM}}^{\mathrm{atm}}}{(n_{\mathrm{cell}}{)}_{\mathrm{LR}}^{\mathrm{atm}}}}\\ & \cdot \left({\displaystyle \frac{(\mathrm{\Delta }{t}_{\mathrm{physics}}{)}_{\mathrm{LR}}^{\mathrm{atm}}}{(\mathrm{\Delta }{t}_{\mathrm{physics}}{)}_{\mathrm{RRM}}^{\mathrm{atm}}}}\cdot \mathrm{w}.\mathrm{c}.\mathrm{t}{.}_{\mathrm{LR}}^{\mathrm{atmphysics}}\right.\\ & \left.+{\displaystyle \frac{(\mathrm{\Delta }{t}_{\mathrm{dynamics}}{)}_{\mathrm{LR}}^{\mathrm{atm}}}{(\mathrm{\Delta }{t}_{\mathrm{dynamics}}{)}_{\mathrm{RRM}}^{\mathrm{atm}}}}\cdot \mathrm{w}.\mathrm{c}.\mathrm{t}{.}_{\mathrm{LR}}^{\mathrm{atmdynamics}}\right).\end{array}\end{array}$$

Figure 6 shows the results of these models, where wall-clock time and simulation length have been converted to throughput (SYPD). The solid lines show the measured throughput of each component as a function of number of computer cores. The dashed lines show the corresponding throughput values predicted by Eqs. (1) and (2). The single LR XS layout ocean throughput value is used as the reference for the ocean, and the single LR XS layout atmosphere throughput value is similarly used as the reference for the atmosphere; these are the only measured data inputs to the performance models. The primary error in the performance model is not accounting for a fall-off in scaling at large core counts. Because this fall-off is small for the atmosphere and ocean components, these simple performance models are accurate and can be used to predict the cost of other model configurations. For example, a uniform high-resolution atmosphere model would use the ne120pg2 grid, which has 6⋅120² elements. Using the same time steps and number of vertical levels as in the NARRM configuration, which has 14 454 elements, for fixed computational resources, the high-resolution atmosphere configuration's throughput would be $\mathrm{6}\cdot {\mathrm{120}}^{\mathrm{2}}/\mathrm{14}\mathrm{454}=\mathrm{5.98}$ times smaller than the NARRM configuration's throughput, where this factor is the quotient of the numbers of elements in each of the two grids.

Figure 7 Comparison of the global spatial RMSE of model climatology (annual and seasonal averages of years 1985–2014) vs. observations with the E3SM Diags package (C. Zhang et al., 2022). The model results are from the first historical member of E3SMv2 (0101), LR (blue triangles), and NARRM (red triangles) and 52 CMIP6 models (r1i1p1f1). The boxes and whiskers show the 25th percentile, 75th percentile, and minimum and maximum RMSE of the CMIP6 ensemble. Quantities include (a) TOA net radiation flux, (b, c) TOA SW and LW cloud radiative effects, (d) precipitation, (e) surface air temperature over land, (f) sea level pressure, (g, h) 200 and 850 hPa zonal wind, and (i) 500 hPa geopotential height. TOA is the top of the atmosphere, SW is shortwave, CRE is cloud radiative effects, LW is longwave, ANN is annual, DJF is December–February, MAM is March–April, JJA is June–August, SON is September–November, and RMSE is root-mean-square error. The climatology of the observations and reanalysis data are calculated from CERES-EBAF Ed4.1 (Loeb et al., 2018) (2001–2014) for (a), (b), and (c); GPCP2.3 (Adler et al., 2018) (1985–2014) for (d); and ERA5 (Hersbach et al., 2020) (1985–2014) for (e), (f), (g), and (h).

Figure 8 Surface geopotential height of (a) LR and (b) NARRM over North America.

Figure 9 Time series of global annual mean surface air temperature anomalies from the ensemble mean of LR (blue) and NARRM (red) historical runs and observational datasets (gray) (National Oceanic and Atmospheric Administration (NOAA) National Climatic Data Center (NCDC), National Aeronautics and Space Administration (NASA) GISTEMP, and HadCRUT4). The model ensemble minimum–maximum ranges are shaded, while the observational minimum and maximum numbers are labeled in the parentheses of legend.

Figure 10 Comparison of climate sensitivities between LR (a, c) and NARRM (b, c) derived from idealized CO₂ forcing simulations. (a, b) time series of global annual mean surface air temperature anomaly from the following simulations, abrupt-4xCO2 (red), 1pctCO2 (blue), and the control (piControl; green). The transient climate response (TCR) is computed as a 20-year average around the time of CO₂ doubling (year 70). (c, d) Gregory regression plots. The estimated effective climate sensitivity (ECS) and effective 2× CO₂ radiative forcing (F) are as labeled.

4 Global climate

As described above, the RRM model is expected to simulate a global climate similar to the LR model for production simulation campaigns since most areas are still covered by the same LR grids. In this section, we will examine whether this is the case for the global mean climate, climate sensitivity, and climate feedback.

For the global climatology, we focus on the last 3 decades (years 1985–2014) of historical simulations when more observational datasets are available. Figure 7 provides an overall comparison of the global mean climate among LR (blue triangles), NARRM (red triangles), and CMIP6 (boxes and whiskers) models as quantified by the uncentered spatial root-mean-square error (RMSE) relative to the observations or reanalysis data. The RMSE numbers are calculated with the E3SM Diags package (C. Zhang et al., 2022) for the first historical member (0101 for LR and NARRM, r1i1p1f1 for CMIP6 models). Figure 7 clearly shows that NARRM and LR simulate very similar annual and seasonal averages. NARRM outperforms LR in the June–July–August (JJA) shortwave (SW) cloud radiative effect (CRE) partly because it better represents low clouds in NA (see Fig. 13 for the example in California). NARRM also simulates slightly better December–January–February (DJF) precipitation compared to LR, partly due to its improved topography (Fig. 8) and orographic precipitation in NA (see Fig. 12b, d, f). For other times (e.g., annual mean, Fig. A3) NARRM and LR precipitation results are very similar. On the other hand, NARRM does not perform as well as LR for some other fields, such as the 200 hPa zonal wind in JJA and September–October–November (SON), which are associated with the increased positive biases in the tropical western Pacific and Amazon (not shown).

Figure 9 compares the long time series (years 1850–2014) of global annual average anomalies in the surface air temperature from the ensemble means of LR and NARRM historical simulations and observational datasets (National Oceanic and Atmospheric Administration National Climatic Data Center (Smith et al., 2008; Zhang et al., 2015), National Aeronautics and Space Administration GISTEMP (GISTEMP Team, 2018; Hansen et al., 2010), and HadCRUT4 (Morice et al., 2012)). Over the whole period, NARRM tracks LR closely, including good agreement with observations until the 1930s and low biases afterwards, which are mainly attributed to too strong aerosol-related forcing and feedback (see Golaz et al., 2022, for details). This is further confirmed by the fact that the global mean effective radiative forcing of anthropogenic aerosols in NARRM and LR are very similar (−1.415 W m⁻² vs. −1.421 W m⁻²), as quantified by a pair of nudged simulations (see Sect. 5.1.2).

Following the CMIP6 DECK protocol (Eyring et al., 2016), we quantify the climate sensitivity and feedback with the abrupt quadrupling of CO₂ (abrupt-4xCO2) and the transient climate response (TCR) with a simulation forced by a 1 % yr⁻¹ CO₂ increase (1pctCO2) relative to the pre-industrial control simulation (piControl). The equilibrium climate sensitivity (ECS) is estimated with the linear regression of TOA radiation change against surface temperature change in a 150-year abrupt-4xCO2 simulation (Gregory et al., 2004). The 2xCO₂ effective radiative forcing (ERF) is computed as the y intercept of the Gregory plot divided by two, which measures the energy imbalance caused by doubling the atmospheric CO₂ concentration while keeping the surface temperature unchanged. TCR, which measures the response on shorter timescales, is derived based on its definition – the average surface temperature change in the 20-year period when the CO₂ concentration doubles from a 1pctCO2 experiment.

Figure 10 depicts the annual mean surface temperature change as a function of time and the Gregory plots from the idealized CO₂ experiments with LR and NARRM. The differences in climate sensitivity between LR and NARRM are very subtle as quantified by both ECS (4.00 K vs. 3.94 K) and TCR (2.41 K vs. 2.44 K).

The regression slope in the Gregory plot (Fig. 10c, d) denotes the total radiative feedback caused by the quadrupled CO₂ concentration. We further apply the radiative kernel method (Soden et al., 2008; Held and Shell, 2012) to decompose the total radiative feedback into non-cloud and cloud feedbacks. The cloud feedback is estimated by adjusting the cloud radiative effect anomalies for non-cloud influences. Overall, NARRM shows a slightly larger ERF (3.22 W m⁻² vs. 2.98 W m⁻²), which accompanied with the similar ECS produces a stronger negative total climate feedback in NARRM. The total climate feedback is −0.74 W m⁻² K⁻¹ for LR and −0.82 W m⁻² K⁻¹ for NARRM, which mainly relates to the slightly weaker positive SW cloud feedback in NARRM than in LR (see Fig. 22).

In summary, the results in this section confirm that NARRM with the hybrid time step methodology simulates largely identical global climate as its corresponding LR configuration and hence satisfies the necessary requirement (i.e., good global climate) of global RRM production simulations we proposed in the introduction.

5 North American results

In this section, we will zoom in over the refined region over North America (NA) and emphasize climate aspects most relevant to the E3SM water cycle scientific goals (Leung et al., 2020) as well as some weaknesses in LR revealed by Golaz et al. (2022). The results will be described for the atmosphere, land, and river models, respectively. Moreover, we will analyze interactions between different components (i.e., land–atmosphere coupling) because these interactions are also expected to change with the resolution increase.

5.1 Atmosphere

5.1.1 Hydrology over the CONUS

First, we look at the overall atmospheric results over the CONUS (20–50^∘ N, 65–125^∘ W) by comparing the spatial RMSEs of the historical ensemble means between LR (blue triangles) and NARRM (red triangles) in Fig. 11. The same metric is used in Fig. 7 for the global results, but we adjust the variables to be more relevant to the CONUS. NARRM generally produces better (as quantified by smaller RMSE numbers) results than LR for these annual and seasonal climatologies, such as SW CRE (Fig. 11a), precipitation (Fig. 11c), and 200 hPa zonal wind (Fig. 11f). Because we have not retuned the physics of NARRM, some deteriorations are expected, for example longwave (LW) CRE in DJF and March–April–May (MAM) (Fig. 11b).

Precipitation and clouds, which are obviously important quantities for the water cycle, are largely improved in NARRM compared to LR. The precipitation patterns are better captured by NARRM as observed at the Sierra Madre Occidental in JJA (Fig. 12 left column) and in the western US in DJF (Fig. 12 right column) due to the better-resolved topography in NARRM (see Fig. 8). The poor representation of marine stratocumulus clouds is a long-standing problem that plagues many ESMs (e.g., Bogenschutz et al., 2023). The underestimation of summertime low clouds (manifested as the excessive TOA shortwave CRE) in the California stratocumulus region is substantially improved with NARRM (Fig. 13). The improvement in this bias is likely due to a reduced bias in the simulated sea surface temperature (SST) in the coupled RRM. To constrain the impact of the RRM on the SST bias, we conduct two additional experiments: one where the LR ocean is coupled to the RRM atmosphere, and another where the RRM ocean is coupled to the LR atmosphere. The simulated SST bias averaged over years 51–100 of the piControl is shown in Fig. 14. Comparison of Fig. 14a and b shows a clear reduction in bias for the NARRM simulation. Figure 14c and d illustrate that the regional refinement in the atmosphere (c) is primarily responsible for the reduction in SST bias; however, comparing to the NARRM result, we see that the bias is further reduced when regional refinement is included in both components, highlighting the advantage of coupled RRM over a single-component RRM. Further details will be described in a future paper by Van Roekel et al. (Luke P. Van Roekel, personal communication, 2023).

Another well-known issue of global ESMs is the poorly captured diurnal propagation of organized mesoscale convective systems (MCSs). Over CONUS, such MCSs originate from the front range of the Rockies in the afternoon and propagate eastward, manifesting as a nocturnal precipitation peak in the central US and contributing as much as half of the summertime rainfall in that region (Riley et al., 1987; Jiang et al., 2006). Both LR and NARRM simulate the summertime nocturnal rain peak in the central US (Fig. 15) because of the new convective trigger method for deep convection (Xie et al., 2019). However, the magnitude is weaker and the area of nocturnal peak extends much larger (almost the whole eastern half of the US) than the observation, which could be caused by remaining propagation or convective trigger deficiencies. Nevertheless, NARRM reduces the underestimation of maximum diurnal cycle magnitude with an 80 % greater value (5.61 mm d⁻¹ vs. 3.10 mm d⁻¹) than LR on the same 100 km grids. The NARRM maximum can be as high as 6.71 mm d⁻¹ on the 25 km grids but still biases low compared to the observation (10.89 mm d⁻¹). This result suggests that a resolution of ∼ 25 km is not adequate to capture the physics driving propagating MCSs, which probably require convection-permitting atmospheric simulations to achieve a good agreement with observations (Caldwell et al., 2021).

Figure 11 The same as Fig. 7 but contrasting LR and NARRM historical ensemble means at the refined CONUS area (20–50^∘ N, 65–125^∘ W). Note that variables shown are adjusted to be more appropriate for the CONUS. (a, b) TOA SW and LW cloud radiative effects, (c) precipitation, (d) total precipitable water, (e) surface air temperature, (f, g) 200 and 850 hPa zonal wind, and (h, i) 500 hPa geopotential height and vertical velocity (pressure). The same reference climatology data are used for the variables also shown in Fig. 7, whereas ERA5 (Hersbach et al., 2020) (1985–2014) is used for (d), (e), and (i).

Figure 12 Comparison of CONUS JJA (a, c, e) and DJF (b, d, f) precipitation geographic patterns from ERA5 reanalysis (a, b), LR (c, d), and NARRM (e, f) historical ensemble means.

Figure 13 Mean TOA shortwave cloud radiative effects at California in JJA of (a) observations (CERES-EBAF Ed4.1), (b) LR (H1-5) minus observation, and (c) NARRM (H1-5) minus observation.

Figure 14 Sea surface temperature (SST) bias (model–observations) simulated by four configurations of E3SMv2 (a) NARRM, (b) LR, and (c) RRM atmosphere coupled to LR ocean and (d) LR atmosphere coupled to RRM ocean. The data are averaged over years 51–100 of the respective piControl simulations.

Figure 15 Mean diurnal phase (local time, colors) and magnitude (color density) of the maximum precipitation in JJA calculated from the first harmonic of 3-hourly total precipitation (mm d⁻¹) for (a) Tropical Rainfall Measuring Mission (TRMM) observations (Huffman et al., 2007), (b) LR (H1-5), (c) NARRM (H1-5) regridded to the same 100 km grids as (b), and (d) NARRM (H1-5) on 25 km grids.

5.1.2 Aerosols

The E3SMv2 LR model (Golaz et al., 2022) simulates too strong aerosol-related forcing, which has been identified as the primary cause of the underestimated warming in the later portion of historical period in Fig. 9. We will examine here if the NARRM configuration helps bring down the biases in aerosols and anthropogenic forcing by better resolving the meteorological and climate fields. In addition, since NARRM employs the hybrid time step approach that eliminates retuning the scale-dependent aerosol parameters used in LR, e.g., the global scaling factor used to constrain the total emission fluxes of natural aerosols (dust and sea salt), which depend non-linearly on the model-resolved small-scale surface winds, we will also discuss the impact of increasing model horizontal resolution on the natural aerosols and total aerosol optical depth (AOD) in NARRM.

1. Impact on anthropogenic aerosols.

   Aerosols in the NARRM configuration are represented in the same manner as in the LR with the enhanced MAM4 (Wang et al., 2020) and improved dust aerosol properties (Feng et al., 2022). Anthropogenic and wildfire emissions used in the LR and NARRM experiments are also from the same input datasets. However, cloud microphysical processes, horizontal advection, and convection can affect aerosol loading within the NARRM high-resolution domain if wet deposition and/or transport are substantially different from the LR configuration (Caldwell et al., 2019). We first compare modeled surface mass concentrations of SO₂, sulfate, black carbon, and organic carbon between the ensemble means of LR and NARRM historical simulations, and evaluate them against ground-based observations from the Clean Air Status and Trends Network (CASTNET) and the Interagency Monitoring of Protected Visual Environments (IMPROVE). The results are shown in Fig. 16. In general, both sets of simulations show a strong correlation with measurements, and biases are very similar between the two for the anthropogenic aerosol species. NARRM simulations slightly increase (less than 7 %) the concentrations of the four species shown here, compared to LR simulations, which may result from less wet removal or vertical transport in the refined mesh.

2. Impact on natural aerosols.

   In addition to aerosol removal, emissions of natural aerosols such as dust and sea salt are highly sensitive to resolution changes due to their strong dependence on the resolved surface wind speeds in the model. Increasing model horizontal resolution normally requires retuning the dust and sea salt aerosol emission factors. For E3SMv1, Feng et al. (2022) showed that without retuning, an increase in the horizontal resolution by a factor of 4 (i.e., from ∼ 100 km in LR to ∼ 25 km in HR) results in about 29 % increase in global dust emissions and an even larger increase in dust AOD of 42 % due to the combined effects from the weakened removal. In contrast, as shown in Table 4, NARRM historical runs simulate nearly the same global mean AODs as LR for all the aerosol species including dust and sea salt, without changing their emission factors. Over the regionally refined CONUS, the mean dust and sea salt AODs are slightly increased (< 5 %). This suggests that NARRM largely retains the performance of LR for the aerosol simulations on the global and regional mean basis without requiring additional retuning of the scale-dependent emission factors.

3. Aerosol spatial variability and extremes.

   On the other hand, NARRM shows improvement over LR in representing aerosol spatial variability and extreme values over the refined mesh region. Figure 17 compares the simulated AOD (550 nm) distributions between LR (0101) and NARRM (0101) historical simulations for the present-day time period of 2000–2014. While both depict a similar general geographical pattern, e.g., higher AODs over the more polluted eastern US than the western part of the country, NARRM captures greater and finer detail in spatial variability than LR, e.g., over the mountainous areas along the Rockies, Sierra Nevada, and Appalachians. The better-resolved AOD variability in NARRM results from the spatial refinement of the resolution-dependent aerosol emission fluxes (natural species), transport, and removal, as discussed above. Compared to the ground-based AOD measurements at the 37 AERONET (Holben et al., 1998) sites (2006–2015), NARRM shows stronger spatial correlation with the observations than LR (Fig. 17c). Both configurations overestimate the mean AOD averaged over the AERONET sites, possibly linked to the weak wet removal in E3SMv2 (Golaz et al., 2022).

   In addition to the improved spatial variability, higher resolution in NARRM also leads to more frequent occurrences of large AOD predictions over CONUS than the LR model, especially over the regions dominated by wind-driven dust or sea salt aerosols. Figure 18 shows an example of the calculated probability density function (PDF) for dust AOD over the major dust source region in the US (32–42^∘ N, 118–108^∘ W; indicated by the $\mathrm{10}{}^{\circ }\times \mathrm{10}{}^{\circ }$ box in Fig. 17b), from both the LR (0101) and NARRM (0101) simulations in 2000–2014, which are remapped to the same 0.25^∘ grid resolution. It is worth noting that the remapping of the LR results to the finer resolution leads to little improvement in the resolved spatial variability in dust AOD. Clearly, NARRM predicts more occurrences of high dust AOD over this region than LR, e.g., 22 % of the dust AODs predicted by NARRM exceed 0.015, which is the top 98th percentile of the LR model predictions remapped to the same resolution. This suggests that LR may significantly underestimate the occurrences of large dust outbreaks in the southwestern US region relative to NARRM due to the unresolved surface winds for dust mobilization in the model. Similarly, NARRM would be more suitable for urban climate or air quality studies for capturing the extremely polluted cases occurring at finer spatial or temporal scales.

4. Effective radiative forcing of anthropogenic aerosols.

   Figure 19 shows the effective radiative forcing of anthropogenic aerosols (ΔF) over CONUS and adjacent ocean areas estimated using nudged LR and NARRM simulations. ΔF is overall negative in both LR and NARRM and dominated by the shortwave component (ΔF_SW in Figure 19b, e). The regional mean ΔF_SW and ΔF_LW are both slightly stronger (more negative for ΔF_SW and more positive for ΔF_LW) in NARRM compared to LR. Over the Pacific Ocean near 20^∘ N and 120^∘ W, ΔF_SW and ΔF in NARRM are much stronger than in LR. This is mainly caused by larger low cloud fraction simulated in NARRM (see Appendix Fig. A4), which causes a larger contrast in droplet number concentration and liquid water path between the PD and PI simulations compared to LR (Fig. A5).

Figure 16 Scatter plots of modeled annual mean surface concentrations of (a) SO₂, (b) sulfate, (c) black carbon, and (d) organic carbon (POM+SOA) from LR (H1-5) and NARRM (H1-5) compared to observations at CASTNET and IMPROVE network surface sites during 2005–2014. The numbers are mean concentration and correlation coefficient (R) for data at the individual sites.

Table 4 Comparison of simulated annual mean AOD (550 nm) between LR (H1-5) and NARRM (H1-5) for the time period of 1985–2014. POM stands for particulate organic matter, BC stands for black carbon, and SOA stands for secondary organic aerosol.

Figure 17 Aerosol optical depth (AOD) at 550 nm from (a) LR (0101) and (b) NARRM (0101) historical simulations averaged over 2000–2014. Panel (c) shows the AOD comparison of the two model simulations with the AERONET observations during 2006–2015. The site locations of AERONET are denoted by the gray dots in panel (a). The gray box in panel (b) denotes the dust region referenced in Fig. 18.

Figure 18 Calculated probability density function (PDF) of the dust AOD predictions from LR (0101) and NARRM (0101) between 2000–2014 remapped to the same 0.25^∘  grid resolution, over the major dust source region in the US (32–42^∘ N, 118–108^∘ W; indicated by the $\mathrm{10}{}^{\circ }\times \mathrm{10}{}^{\circ }$ box in Fig. 17b).

Figure 19 Anthropogenic aerosol effects simulated by the nudged LR and NARRM simulations.

5.1.3 Cloud and cloud feedback

Here we examine the impact of increased horizontal resolution over NA on the simulated clouds and their radiative effects with the LR and NARRM historical simulations and on cloud feedback changes with the quadrupling 4xCO₂ simulations.

Figure 20 compares the cloud cover between the E3SM CALIPSO (Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observation) simulator output and the GCM-Oriented CALIPSO Cloud Product (CALIPSO-GOCCP) (Zhang et al., 2023). Cloud cover and cloud thermodynamic phase are diagnosed with the same algorithm in the CALIPSO simulator and CALIPSO-GOCCP data, facilitating consistent model–observation comparisons (Chepfer et al., 2008, 2010; Cesana and Chepfer, 2013). Observed total cloud cover is larger over the NA polar region than the CONUS in the CALIPSO-GOCCP data. Total cloud cover is larger than 60 % over the eastern Pacific Ocean, northern Atlantic Ocean, and Arctic Ocean. Strong land–ocean contrast is observed – liquid phase clouds dominate over the ocean, while ice phase clouds prevail over the land in areas such as Greenland and CONUS. Compared to CALIPSO-GOCCP, LR overestimates total cloud cover at NA high latitudes and the western CONUS and underestimates it in Greenland, particularly over Baffin Bay and near the Greenland coast. The underestimated cloud cover over the western coast of the CONUS is also notable, which is consistent with the previous discussion on Fig. 13. The excessive modeled total cloud cover is primarily attributed to the positive biases in the liquid cloud over the polar region. The positive biases of ice cloud cover contributes to the biases over the mountainous regions in the western NA. On the other hand, ice clouds are underestimated over Greenland and northern Canada.

Over land, NARRM displays improvements relative to LR in western NA and the Arctic for both cloud phases. For instance, the ice cloud biases are significantly reduced from Alaska to the western CONUS (Fig. 20f, i), and the liquid cloud deficiencies over Alaska and Greenland are generally decreased (Fig. 20e, h). The better represented topography in NARRM (Fig. 8) is probably the key factor of these NARRM improvements. The impact of increased horizontal resolution on simulated cloud phase is also noted in the E3SMv1 model with the CALIPSO simulator (Zhang et al., 2019), where increased horizontal resolution also slightly decreases simulated liquid and ice clouds at temperatures warmer than −40 ^∘C in the Arctic region.

Over ocean, NARRM substantially improves the stratocumulus clouds to the west of coastal regions. NARRM also moderately outperforms LR in representing liquid clouds over the North Atlantic to the west of Greenland. This is related to the warmer (∼ 1.5 ^∘C) NARRM surface air temperature over the Labrador Sea. This warmer NARRM surface air temperature is consistent with the decreased sea ice concentration in that region (not shown), which is somewhat expected as an advantage of refining grids for both atmosphere and ocean and sea ice (see Fig. A1). Further process-level analysis is necessary to fully understand this LR-NARRM model behavior change and will be reported in separate papers.

Figure 21 compares the simulated SW and LW CRE with the CERES-EBAF Ed4.1 observations. Large negative SW CRE biases and positive LW CRE biases are shown over the Arctic land area and western coast of NA (i.e., Alaska to Oregon) in LR, which mainly result from the overestimated cloud cover. These biases are substantially reduced in NARRM, primarily owing to the improved cloud cover (Fig. 20). Given the reduced negative bias of marine stratocumulus clouds in NARRM near the western coasts of the CONUS, simulated SW CRE is also largely improved. As discussed by Golaz et al. (2022), sea ice concentration is largely overestimated over the North Atlantic Ocean in LR. The too large sea ice extent leads to weaker SW and LW CRE than observed in the Labrador Sea. This is primarily because of the brighter and colder sea ice surface in LR that reflects more SW radiative fluxes and emits less LW radiative fluxes than the observations under clear-sky conditions (not shown). Compared to LR, the maximum positive bias in NARRM sea ice extent in Labrador Sea is greatly alleviated. The better simulated sea ice extent reduces the biases of overly reflective clear-sky SW radiation and the insufficient clear-sky outgoing LW radiation. With generally comparable all-sky SW and LW radiative fluxes between LR and NARRM, those reduced clear-sky biases thus lead to a better CRE in NARRM over Labrador Sea.

Given the improved historical cloud cover and cloud radiative effects over NA, we further examine the regional climate feedbacks over this region in Fig. 22. Relative to the global mean value, the total climate feedback over NA is more negative from NARRM than from LR. This mainly results from the more negative Planck feedback and less positive SW cloud feedbacks. The more negative Planck feedback is related to the stronger surface warming over the northeastern Pacific (not shown).

Figure 23 shows the spatial distribution of cloud feedback of LR and NARRM in the NA region. Due to the different cloud types over land and ocean, we report their regional averages separately. Notably, the total land cloud feedback is 0.41 W m⁻² K⁻¹ smaller in NARRM (0.26 W m⁻² K⁻¹) than in LR (0.67 W m⁻² K⁻¹), which is dominated by the reduced SW cloud feedback over the northeastern US (Fig. 23e). Further examination indicates this reduction is mainly related to the weaker reduction in low cloud cover under warming in NARRM. Figure 20g shows that the overestimated cloud cover is slightly alleviated over the Arctic land region in NARRM, implying that the lower mean state cloud cover might contribute to a weaker cloud reduction under warming there. Over ocean, NARRM presents a stronger SW cloud feedback and a weaker LW cloud feedback over the marine low cloud regime, leading to a small change in total cloud feedback. Across the CSS/WGNE Pacific Cross-Section Intercomparison (GPCI) transect (Teixeira et al., 2011), NARRM tends to show a weaker positive cloud feedback near the coast and more positive cloud feedback off the coast. These factors suggest that the regional refinement can significantly affect the regional cloud responses under warming and the predictability of regional climate.

Figure 20 Spatial distribution of annual mean total cloud cover (a) and cloud cover in liquid phase (b) and ice phase (c) from the CALIPSO-GOCCP data. The cloud cover biases in LR (H1-5) and NARRM (H1-5) historical simulations (1985–2014) are shown in (d)–(f) and (g)–(i), respectively. Simulated cloud cover and cloud thermodynamic phase are derived by the CALIPSO simulator. Climatology data of CALIPSO-GOCCP version 3.1.2 (Chepfer et al., 2010) from 2006–2018 are used in the model evaluation.

Figure 21 Spatial distribution of the observed shortwave cloud radiative effect (a), the longwave cloud radiative effect (b), and the simulated cloud radiative effect biases in LR (H1-5) (c, d) and NARRM (H1-5) (e, f) historical simulations (1985–2014). The observed cloud radiative effect is from the CERES-EBAF Ed4.1.

Figure 22 Mean global and NA climate feedbacks of LR and NARRM decomposed using radiative kernels (e.g., Soden et al., 2008; Held and Shell, 2012).

Figure 23 Spatial distribution of total (a, d, g), SW (b, e, h), and LW (c, f, i) North American cloud feedbacks for LR (a–c), NARRM (d–f), and the difference between NARRM and LR (g–i). The GPCI transect is denoted by the dashed black line. The average values over land and ocean are labeled in the brackets.

5.1.4 Extratropical cyclone

One of the primary motivations for pushing climate simulation resolution is to potentially better capture extremes. Extratropical cyclones (ETCs) are a major weather extreme phenomenon at middle and high latitudes, bringing with them strong winds and precipitation that can exert substantial societal impacts along their pathways over days and over hundreds to thousands of kilometers. Climate changes are likely to induce changes to the dynamical and physical characteristics of ETCs as well as their geospatial distribution (e.g., Bengtsson et al., 2006; Ulbrich et al., 2009). Projections of such future changes rely heavily on numerical climate and ESM models. Their skills in simulating major weather systems like ETC have been carefully scrutinized by modeling centers and by the climate science community in conjunction with the major intercomparison campaigns (Greeves et al., 2007; Chang, 2013). While the conventional climate models with grid resolution around 100 km show reasonable skill in producing ETC frequency and spatial track density, it has also been found that higher-resolution models are better capable of capturing more intense ETCs (e.g., Jung et al., 2006), which is critical for using ESM to project future climates, as growing evidence shows that global warming tends to shift the weather spectrum to the more extreme end (Melillo et al., 2014). Here, we will demonstrate the benefits of higher resolution in simulating ETCs in a regionally refined setting by comparing the results from NARRM with those from LR simulations against the ETC activities derived from the ERA5 reanalysis.

The ETC tracks and statistics can be obtained using automated identification and tracking algorithms (e.g., Blender and Schubert, 2000; Geng and Sugi, 2001; Bengtsson et al., 2006; Jung et al., 2006; Ullrich and Zarzycki, 2017; Ullrich et al., 2021). The objective identification and tracking also make it suitable to compare ETC activities and statistics derived from different data sources in particular for model evaluations. The algorithms usually identify and track the spatial features of a meteorological variable, such as mean sea level pressure (MSLP) or 850 hPa vorticity, that can characterize the structure of cyclones and their movements. In this work, we use a community feature detection and tracking framework, TempestExtremes (Ullrich and Zarzycki, 2017; Ullrich et al., 2021), to derive ETC activities from 6-hourly MSLP data during the period of 1985–2014 from the E3SM simulations and the ERA5 reanalysis. Considering higher-resolution data can more accurately identify the storms and their tracks (Blender and Schubert, 2000; Geng and Sugi, 2001), all the model and reanalysis data are placed on $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ grids to feed the tracking software. The algorithm takes two steps. First, a candidate cyclone is detected when a minimum MSLP feature is enclosed by a contour of 200 Pa interval within 6^∘ of the center. Candidates within 6^∘ of one another are merged, with the lower center pressure taking precedence. The candidates are then stitched together to define the tracks if the features persist for at least 60 h with a maximum gap of at most 18 h. From the start to the end, a candidate cyclone must travel at least 12^∘ great circle distance to qualify as an ETC.

Over the NARRM high-resolution domain, ETCs are most active during winter. Figure 24 shows the mean DJF track density for the models and the analysis derived by casting the computed ETC track data onto $\mathrm{5}{}^{\circ }\times \mathrm{5}{}^{\circ }$ grids. The tracks are mostly concentrated over the northeastern Pacific and northwestern Atlantic that form the well-known storm tracks. Both LR and NARRM simulations capture these main features to a large extent. There are, however, notable differences between LR and NARRM over these oceanic storm tracks. The track densities are clearly underestimated in the LR simulations inside the refined region, except for the Atlantic storm track in the coupled mode. NARRM clearly produces higher ETC track density than LR does for both sections of the oceanic storm tracks and mostly agrees better with the ERA5 data, although in the coupled mode the track density tends to be overestimated. The shape and orientation of both the Pacific and Atlantic storm tracks are much better produced by NARRM in the coupled mode. Given the significant differences from both coupled LR (Fig. 24d) and uncoupled NARRM (Fig. 24c), it is reasonable to believe that the better captured storm track shapes in the coupled NARRM are due to interactions with the refined ocean (Fig. A1). Several secondary centers of active ETCs in the ERA5 over land are also reproduced by the models, including the active regions over the Great Lakes and Hudson Bay, although the densities are overestimated in the models (more so in the NARRM). It is worthwhile mentioning that the NARRM simulations are able to produce the chain of secondary centers to the east of the Rocky Mountains that are also present in the ERA5 reanalysis but are largely missing in the LR simulations. This is presumably due to the NARRM's better-resolved mountainous terrain features, a benefit by design.

The benefit of grid refinement can be further seen in Fig. 25, which shows the histograms of the ETC as a function of the minimum center pressure and the maximum deepening rate during its lifetime. All events within the refined region bounded by the dashed black lines as shown in Fig. 24 are used to compute these statistics. The maximum deepening rate is defined as the maximum 6-hourly center pressure drop, normalized by $\sin \left({\mathit{\phi }}_{\mathrm{ref}}\right)/\sin \left(\mathit{\phi }\right)$, with φ being the latitude and φ_ref the reference latitude at 45^∘ (see also in Jung et al., 2006). Clearly the NARRM very closely reproduces the number of intense cyclones (minimum center MSLP <960 hPa), while unsurprisingly the LR model underproduces. This is true in coupled and uncoupled modes. Both LR and NARRM simulations overestimate the number of weaker ETCs. On a similar note, the observed number of rapid-growth cyclones is closely reproduced by the NARRM simulations but is clearly underestimated by a large margin in the LR simulations.

Figure 24 Mean extratropical cyclone track density in the DJF season between 1985 and 2014 from (a) ERA5, (b, c) AMIP, and (d, e) historical simulations of LR and NARRM. Dashed black lines denote the western and eastern boundary of the refined region.

Figure 25 Histogram of extratropical cyclone minimum center sea level pressure and maximum 6 h deepening rate in the DJF season between 1985 and 2014 (a, b) AMIP and (c, d) historical simulations. LR is shown in dotted, NARRM in dashed, and ERA5 in solid black lines.

5.2 Land and river

5.2.1 Snowpack

Natural storage provided by mountain snowpack is central to water supply reliability in the western US (Siirila-Woodburn et al., 2021). To evaluate model skill in representing this critical hydroclimate benchmark variable, intra-annual snowpack dynamics are evaluated using the methodology known as the snow water equivalent (SWE) triangle (Rhoades et al., 2018a, b). The seven metrics that make up the SWE triangle attempt to distill management-relevant aspects of the accumulation and ablation of snowpack (e.g., peak water volume and snowmelt rate) for any arbitrary gridded SWE dataset. Five HUC2 basins of the mountainous western US are used to derive five-member ensemble and basin average evaluations of LR and NARRM fully coupled historical simulations and are compared with ERA5. All datasets are bi-linearly regridded using the Earth System Modeling Framework (ESMF) to 0.25^∘ resolution prior to masking and computing the basin-average SWE triangle metrics.

NARRM provides enhanced winter (DJF) climatological representation of the spatial variability of SWE across the CONUS relative to LR (Fig. 26a, b). This is seen through higher SWE magnitudes and more granular spatial structures in NARRM compared with LR, particularly in coastal mountain ranges such as the Cascades and Sierra Nevada, and corroborates a long history of ESM studies that highlight the critical importance of horizontal resolution (≤ 0.25^∘) in properly representing the mountainous hydrologic cycle (Demory et al., 2014; Rhoades et al., 2017; Kapnick et al., 2018; Palazzi et al., 2019; Bambach et al., 2022; Rhoades et al., 2022). As shown through the more granular intra-seasonal perspective of the SWE triangle metrics, certain aspects in the snowpack dynamics are improved with NARRM (e.g., peak water volume), namely in the Pacific Northwest and California (Fig. 26c). With that said, some E3SM SWE biases are not ameliorated with horizontal resolution and may arise due to the combination of higher winter season precipitation (Fig. 12) and a general cool bias (Fig. 9) in both the LR and NARRM fully coupled historical simulations.

Figure 26 Climatological DJF average snow water equivalent (SWE) as simulated by E3SMv2 with (a) LR and (b) NARRM over the 1985–2014 period. (c) LR (blue) and NARRM (red) SWE triangle metrics for five HUC2 basins within the mountainous western US compared with ERA5 (gray). Black bars at the end of each histogram represent the mean 95 % confidence intervals.

5.2.2 Runoff and evapotranspiration

NARRM better captures spatial variability in land hydrologic processes, as indicated by two most important land hydrologic variables, total runoff (Fig. 27) and evapotranspiration (ET) (Fig. 28). For instance, in the coastal Pacific regions, the Rocky Mountains block atmospheric moisture from ocean to inland areas and lead to two distinct hydrologic regimes: a wet regime in the western mountains and a dry regime in the eastern mountains. This abrupt spatial shift from a wet to dry hydrologic regime can be clearly seen in the composite runoff map from the Global Runoff Data Center (GRDC) (Fekete et al., 2011), Fig. 27a, and the observed evapotranspiration map from the Moderate Resolution Imaging Spectroradiometer (MODIS) satellite observations (Running et al., 2017), Fig. 28a. Note that the GRDC runoff map is not completely based on the observational data since runoff measurements are not available at the regional or global scales due to technical and economic limitations. It is nevertheless a more realistic estimate than any model simulations because it was first generated with a monthly hydrologic model (hence producing spatiotemporal variability) and then corrected for bias against discharge measurements at thousands of river gauges (Fekete et al., 2011). This abrupt shift of hydrologic regime around the Rocky Mountains, along with the other spatial variations, is much better resolved in the NARRM simulation than LR, which is the case for both simulated runoff and evapotranspiration, as shown in Figs. 27b, c and 28b, c. The NARRM simulated spatial patterns are thus more realistic than the LR ones over NA. Over the remaining regions of the globe, the NARRM and LR simulated spatial patterns are quite similar to each other in terms of both runoff and evapotranspiration (not shown).

The simulation biases in runoff and evapotranspiration are further examined in terms of absolute biases, i.e., the absolute difference between the simulated and “benchmark” values. Here the GRDC runoff and MODIS ET data are used as the benchmark data. Figures 27d and 28d show the maps of absolute bias difference, i.e., the difference between the absolute biases in the LR simulation and those in the NARRM simulation (former subtracting latter), for annual mean runoff and evapotranspiration, respectively. For a specific grid cell in these two maps, a positive difference means the absolute bias in the LR simulation is larger than that in the NARRM and vice versa. It appears that there are more absolute biases in LR than NARRM over both the western and eastern US. Using the longitude 100^∘ W as the divide, the average absolute bias differences (positive indicates NARRM has less overall absolute bias than LR) are 22.8 and 0.9 mm yr⁻¹ over the western and eastern US, respectively, for annual mean runoff, and are 21.6 and 18.5 mm yr⁻¹ over the western and eastern US, respectively, for annual mean evapotranspiration. When compared to LR, NARRM can thus help reduce simulation biases in hydrologic variables.

Figure 27 Annual mean runoff from (a) GRDC, (b) LR, and (c) NARRM simulations and (d) the differences in absolute biases between LR and NARRM (LR−NARRM); a positive value suggests the absolute bias in the LR simulation is larger than that in the NARRM, while a negative value indicates the opposite.

Figure 28 Annual mean evapotranspiration from (a) MODIS, (b) LR, and (c) NARRM simulations and (d) the differences in absolute biases between LR and NARRM (LR−NARRM); a positive value suggests the absolute bias in the LR simulation is larger than that in the NARRM, while a negative value indicates the opposite.

5.2.3 Streamflow

Streamflow simulations are typically affected by multiple sources of uncertainties, such as the biases in the simulated runoff, the uncertainties in the river model parameters (e.g., river network topology, channel geometry, Manning's roughness coefficients), and water demand data (Li et al., 2013, 2015a, b; Zhou et al., 2020). For river network topology, the 0.5^∘ and 0.125^∘ resolution river network data are used for the LR and NARRM simulations, respectively, as shown in Fig. 29a, b. It is expected that a higher-resolution river network data can represent rivers more smoothly and hence more realistically. Another benefit of higher resolution river network data is to enable more extensive streamflow validation. Terrestrial water fluxes, particularly surface runoff and streamflow, are dominated by gravity and controlled by topography and hence mostly follow irregular watershed boundaries. In most land surface and ESMs, including E3SM, regular lat–long grids are used to resolve spatial heterogeneity for both runoff and river processes to be compatible with the other land and atmospheric components (Lawrence et al., 2019; Golaz et al., 2022). Both the magnitude and timing of streamflow at each river gauge are dominated by the corresponding upstream drainage area. Streamflow simulations are thus largely affected by the discrepancies between the watershed boundaries and regular lat–long grids. These discrepancies can be significantly reduced with higher-resolution river network data. For example, in this study a 10 % discrepancy threshold is used to select the river gauges for validating streamflow simulations; i.e., the relative difference between the real upstream drainage area of a river gauge and that estimated from a lat–long grid-based river network should not exceed 10 %. Over the NA domain, 615 river gauges satisfy the requirement for the 0.5^∘ resolution river network, whilst 2924 river gauges satisfy the requirement for the 0.125^∘ resolution river network. There are 563 river gauges that simultaneously satisfy the requirement for both resolutions.

Figure 29 Simulated annual mean streamflow at 563 river gauges in NA compared against The Global Streamflow Indices and Metadata (GSIM) database. (a) River network used in LR, demonstrated by the mean annual discharge. (b) River network used in NARRM, demonstrated by the mean annual discharge. (c) Simulated annual mean streamflow against observed for LR and NARRM.

Figure 30 Simulated JJA (a, c, e) and DJF (b, d, f) mean streamflow at CONUS river gauges compared against The Global Streamflow Indices and Metadata (GSIM) database. (a, b) Relative bias of LR. (c, d) Relative bias of NARRM. (e) The difference in absolute relative bias between LR (absolute value of a) and NARRM (absolute value of c) for the JJA season; a positive value (purple) indicates LR has greater absolute bias than NARRM. (f) The same as (e) but for the DJF season

Figure 29c shows the comparison between the annual mean observed and simulated streamflow over these 563 river gauges. Overall, both simulations produce the long-term average streamflow reasonably well across these gauges. NARRM performs noticeably better (closer to the red 1:1 line) for the top four gauges with the largest discharges. An additional analysis (figure not shown) indicates that LR produces greater absolute bias than NARRM in 330 out of 563 gauges (about 60 %) in the streamflow simulation. Interestingly, it appears that the overestimation and underestimation of JJA and DJF streamflow are concentrated in the western and eastern US, respectively, for both LR and NARRM, as shown in Fig. 30a–d. Figure 30e, f displays the difference in absolute biases between LR and NARRM (the former subtracted from the latter) at individual river gauges for the JJA and DJF seasons. Positive differences (indicating greater bias in LR than in NARRM, purple color) dominate over most gauges in the eastern US during JJA and over the CONUS during DJF. Taken together, Figs. 29 and 30 suggest an overall better performance of NARRM despite all the uncertainties.

5.3 Land–atmosphere coupling

Accurate representation of the interactive processes between the land surface, planetary boundary layer (PBL), and clouds and precipitation is an ongoing challenge for current state-of-art climate models. Here we assess the land–atmosphere (L-A) coupling in LR (H1-5), LR (A1-3), NARRM (H1-5), and NARRM (A1-3) using the 9-year warm-season observations at the Atmospheric Radiation Measurement (ARM) Southern Great Plains (SGP) site following Tao et al. (2021). Before the detailed analysis of L-A coupling, we first examine the seasonal variations of daytime mean surface heat fluxes from May to August during 2004–2012. As shown in Fig. 31, the evaporative fraction (EF) is in general underestimated in model simulations except for LR (H1-5). The much lower simulated EF compared with the ARM observations is mainly attributed to a large negative bias in surface latent heat fluxes (LH). Different from the other simulations, the surface sensible heat flux (SH) is significantly underestimated in LR (H1-5) from May to early July. As a result, the simulated daytime mean EF is higher than that from the observations. Overall, the surface state and fluxes are better reproduced in the historical runs than in the AMIP runs, where both LR (A1-3) and NARRM (A1-3) show a significant negative bias in LH and EF persisting since July. This is surprising and requires further analyses, which is beyond the scope of this paper. In the following, we focus on two local convective regimes and diagnose model behaviors using the local coupling metrics (Santanello et al., 2018).

During the selected 9-year period, 165 and 154 clear-sky days are classified from LR (A1-3) and NARRM (A1-3), respectively (Table A1). This is double the 66 clear-sky days identified from ARM observations. However, the occurrence frequency of shallow cumulus (ShCu) days is much lower in these AMIP runs compared with that observed. For ShCu, only 6 and 5 d are identified in LR (A1-5) and NARRM (A1-5), respectively (not shown). For the historical runs, the number of selected clear-sky days from both LR (H1-5) and NARRM (H1-5) are comparable to that observed, but the occurrence frequency of ShCu days is still low. As we are targeting a statistical and climatological comparison between the long-term ARM data and climate model simulations, we extend the analysis period to 1980–2012 for model simulations on ShCu days due to the limited sample size between 2004 and 2012.

Figure 32 shows the composite clear-sky day mixing diagrams (Santanello et al., 2009), which relates the conservative variables, potential temperature (θ), and total water-specific humidity (q) to the water and energy budgets and the growth of planetary boundary layer (PBL). The coevolution of Lvq and Cpθ (07:30 to 17:30 LST) is decomposed by vector components that represent the integrated fluxes of heat and moisture from the land surface (V_sfc), the advection (V_adv), and the entrainment at the PBL top (V_ent as a residual). Six metrics are derived from these diagrams and summarized in Table 5. In general, although the differences among various model simulations are minor, several common model biases are noted when compared with the ARM observations. For example, the model-simulated clear-sky days are featured with too warm and too dry conditions in the early morning (07:30 LST). The θ_sfc in the two AMIP runs and historical runs, for both LR and NARRM, is about 3 and 2 times that which was observed, respectively. The high θ_sfc indicates that more energy at the surface goes to heating rather than moistening. Moreover, the E_LH is significantly overestimated in the model simulations, which is about 5 (3) times that observed in the two AMIP runs (historical runs). The much higher simulated E_LH suggests that the entrainment heating and drying dominates the surface fluxes on the simulated clear-sky days, which supports rapid and deep PBL growth in models. Different from the observations, the simulated advection tends to cool and dry the mixed layer, but the overall impact is much smaller compared to those from the surface and entrainment.

Figure 33 shows the daytime evolution composites of PBL, lifting condensation level (LCL), and LCL deficit (PBL top height minus LCL) on clear-sky and ShCu days. Both PBL and LCL on model-simulated clear-sky days are much higher than on the ARM-observed clear-sky days. The model behaviors are in general consistent among different simulations, except that the bias of LCL is significantly lower in LR (H1-5) compared with the others. In models, the PBL grows rapidly after sunrise on clear-sky days, corresponding to the large warm and dry air entrainments that dominate the PBL budget (Fig. 32). But the too warm and too dry early morning surface conditions lead to an even higher LCL on model-simulated clear-sky days. The PBL never reaches the LCL, with a negative LCL deficit throughout the day, which supports clear skies. The diurnal evolution of LCL on the ARM-observed ShCu days is similar to that on clear-sky days, but the development of the PBL is much more vigorous. As a result, the PBL is deep enough to touch the LCL for cloud formation around noon. Different from the observed results, the daytime evolution of PBL is much weaker on ShCu days than on clear-sky days in all model simulations. However, the decrease in LCL from clear-sky days to ShCu days is even greater, where the growth of PBL is high enough to touch the LCL for cloud formation. Note that the models simulate a positive LCL deficit at around 09:00 LST, a few hours earlier than that in the observations. To summarize, ShCu forms as a result of strong surface SH fluxes that drives the rapid development of PBL in observations, while in models ShCu results from a relatively more humid lower troposphere that leads to a lowered LCL. Differences among various model simulations are pretty minor.

Figure 31 The seasonal variation of 2004–2012 daytime mean (06:00–18:00 LST): (a) surface sensible heat flux (SH), (b) surface latent heat flux (LH), and (c) surface evaporative fraction (EF, defined as [LH/(LH+SH)]) from ARM observations (black), LR (A1-3) (green), LR (H1-5) (blue), NARRM (A1-3) (orange), and NARRM (H1-5). A moving average of 30 d is applied to smooth out short-term fluctuations and highlight longer-term trends.

Figure 32 Clear-sky-day mixing diagram of the PBL conservative variables, Lvq vs. Cpθ, during the daytime evolution from ARM observations (black), LR (A1-3) (green), LR (H1-5) (blue), NARRM (A1-3) (orange), and NARRM (H1-5) (red). Dots denote the composite hourly means from 07:30 to 17:30 LT. The text annotations depict the vector component contributions from surface (V_sfc), advection (V_adv), and entrainment fluxes (V_ent) to the evolution.

Table 5 The surface (β_sfc) and entrainment (β_ent) Bowen ratios, the entrainment ratio of heat (E_SH) and moisture (E_LH), and the advective flux ratio of heat (A_SH) and moisture (A_LH) from the ARM observations, LR (A1-3), LR (H1-5), NARRM (A1-3), and NARRM (H1-5) on clear-sky days. The flux values (W m⁻²) are derived using the mixing diagram theory and surface, advection, and entrainment flux vectors depicted in Fig. 32.

Figure 33 Composite daytime evolution of (a) PBL, (b) LCL, and (c) LCL deficit (PBL minus LCL) from ARM observations (black), LR (A1-3) (green), LR (H1-5) (blue), NARRM (A1-3) (orange), and NARRM (H1-5) (red) on clear-sky days. Panels (d)–(f) are the same as panels (a)–(c) but on shallow cumulus days.

6 Conclusions and discussion

A primary Earth system model (ESM) advancement is to represent the spatially continuous world more realistically on discretized grids, which often requires constantly increasing the finest scale of explicitly resolved processes within the computational limit. Before uniformly high-resolution global models solve their severe computational challenge for climate simulation campaigns, the multiresolution ESM (e.g., regionally refined model (RRM)) is a natural alternative for these campaigns. Nevertheless, it has been over a decade since such a multiresolution method (e.g., Ringler et al., 2008) was proposed.

To our knowledge, this is the first study with a global ESM that has accomplished the CMIP6 climate simulation campaign with a fully coupled RRM configuration – a potentially significant step in the long journey of improving the explicitly resolved resolution of climate simulations. The key to this success is the application of the hybrid time step strategy (i.e., merging the high-resolution dynamics time step with the low-resolution physics time step) in the atmosphere model, which mitigates the negative impacts caused by the persistent poor scale-aware problem of atmospheric physics in a multi-scale framework (e.g., RRM). The powerful aspect of RRM is that it typically only costs ∼ 10 %–20 % of the globally uniform high-resolution model, substantially reducing the computational burden of production simulations. This is particularly important for high-resolution ensemble simulations, which are necessary to account for the internal variability of the climate system, but whose cost would otherwise be prohibitive. On the global scale, we show that NARRM reproduces the LR climate well. Within the high-resolution domain (i.e., North America), NARRM displays more improvements than deteriorations relative to LR. Furthermore, some of the NARRM improvements (e.g., marine shallow cumulus clouds over California and mixed-phase clouds near the Arctic) are attributable to the better-captured coupling processes, highlighting the strength of refining multiple components over a single component. The main detailed findings are as follows.

• The new dry baroclinic idealized test (Hughes and Jablonowski, 2023) allows us to test the NARRM grid with the stand-alone atmospheric dynamical core and confirms that the NARRM mesh is numerically stable and that the results are reasonable compared to the LR and HR grids (Fig. 4).

• By employing the EAM hybrid time step method, NARRM successfully matches the global climate (including climatology, time series, and climate sensitivity and feedback) simulated by LR (Figs. 7, 9, 10, 22) without retuning physics parameters.

• Within the high-resolution region over the CONUS, precipitation and clouds are largely improved in NARRM compared to LR (Figs. 11a, c, 12, 13) due to the better topography in NARRM (Fig. 8) and/or reduced sea surface temperature biases.

• Refining the atmospheric grid spacing from 100 to 25 km is not adequate to improve the diurnal propagation of organized MCSs over the CONUS (Fig. 15).

• NARRM retains the LR performance of aerosol simulations on a global scale and regional mean basis without retuning of the scale-dependent aerosol emissions. Over the refined mesh, NARRM improves the simulated aerosol spatial variability and predictions of extreme polluted cases (e.g., the upper tail of AOD distribution) (Fig. 18). On the other hand, the refined grid resolution does not eliminate the high biases in aerosol loadings and effective radiative forcing inherited from the LR model.

• NARRM generally simulates better cloud cover than LR for both liquid and ice phase clouds. Over land (e.g., western NA and Greenland), this improvement is likely related to topography, whereas over ocean it is attributed to air–sea interactions.

• NARRM produces a comparable global mean cloud feedback to LR but a less positive cloud feedback over the NA (Fig. 22). The reduction in cloud feedback there mainly relates to the shortwave component. The total cloud feedback over coastal California does not change much due to the compensation between the shortwave and longwave components (Fig. 23).

• While both the LR and NARRM simulations are to a large extent able to capture the spatial and statistical distributions of the observed extratropical cyclone (ETC) activities, NARRM shows a particularly improved skill when simulating the ETC activities along the oceanic storm tracks and over the mountain range to the east of the Rocky Mountains (Fig. 24). NARRM in coupled mode outperforms all other configurations (LR and uncoupled NARRM) when simulating the shape and orientation of the oceanic storm tracks within the NARRM high-resolution domain due to the coupling with the refined ocean surrounding North America. NARRM in general produces more ETCs than LR and overestimates the total number of cyclones compared to the ERA5 reanalysis. More importantly, for intense or rapidly developing cyclones, the NARRM simulations are in close agreement with the observations, whereas the LR simulations are mismatched by a significant margin (Fig. 25).

• NARRM appears to better represent the spatial variability in land hydrologic processes by resolving the land features more realistically over the western US (Figs. 27, 28). With higher grid resolution, NARRM can better capture surface topography that dominates surface water flows across hillslopes and through rivers, and hence not only improves the river model performance but also provides more precise river gauge geo-referencing information for streamflow validations (Figs. 29, 30).

• NARRM provides enhanced winter (DJF) climatological representation of the spatial variability of snow water equivalent (SWE) across the CONUS relative to LR (Fig. 26ab) as a result of higher SWE magnitudes and more granular spatial structures in NARRM. Certain biases (e.g., peak water volume) in snowpack are reduced in NARRM compared with LR (Fig. 26c).

• Over the ARM SGP site during warm seasons, the surface conditions are warm and dry on the model-simulated clear-sky days, with overestimation in both the PBL height and the LCL (Figs. 32, 33), while the ShCu days in models result from a much moister environment compared with that in the observations (Fig. 33). In general, the surface properties and fluxes are better reproduced in the historical runs than in the AMIP runs (Fig. 31), showing only limited impact due to resolution.

Besides the NARRM configuration illustrated in the present study, E3SMv2 has been successfully run with RRM meshes with finer grids located in other regions (Antarctic, Arctic, and southeastern Pacific). We expect that the hybrid time step strategy is a general approach that can be applied to these RRMs to simulate high-resolution climate in different areas. With that in mind, we streamlined the process of creating new RRM configurations to facilitate broader RRM applications in the next phase of the E3SM project. Depending on the goal of RRM simulations, further improvements over the refined domain can be achieved via additional parameter tuning. In that case, nudging the outside coarser domain may be necessary to avoid severe degradations of the climate there. Such nudging capability is available in E3SM (e.g., Tang et al., 2019), and one has the option to nudge towards the data from the reanalysis product or low-resolution E3SM simulation. Lastly, we highlight that this paper serves as an overview of the NARRM atmosphere, land, and river models. More in-depth analysis is planned to be reported in follow-up papers.

Appendix A

Figure A1 North American RRM (NARRM) grids for ocean and sea ice.

Figure A2 The same as Fig. 11 but for the global spatial RMSE of model climatology from the EAMv1 LR (blue triangles) and the EAMv1 RRM (red triangles) with high-resolution grids over the CONUS and both physics parameters and time steps tuned for high resolution. The details of these two EAMv1 simulations are documented in Tang et al. (2019).

Figure A3 Comparison of global annual mean precipitation geographic patterns from (a) GPCP2.3, (b) LR, and (c) NARRM historical ensemble means.

Figure A4 Cloud fractions over North America simulated by the nudged LR (a, b, c) and NARRM (c, d, e) for low clouds (CLDLOW; a, d), middle clouds (CLDMED; b, e), and high clouds (CLDHGH; c, f).

Figure A5 Present-day (PD) minus pre-industrial (PI) changes in cloud droplet number concentration (a, d), liquid water path (LWP), and ice water path (IWP) in North America calculated from the nudged LR (a, b, c) and NARRM (d, e, f) simulations.

Table A1 Definition criteria and sample size of the clear-sky regime (Clear) and shallow cumulus regime (ShCu) based on the ARM observations and four different model simulations.