The DOE E3SM Coupled Model Version 1: Description and Results at High Resolution

Abstract

This study provides an overview of the coupled high‐resolution Version 1 of the Energy Exascale Earth System Model (E3SMv1) and documents the characteristics of a 50‐year‐long high‐resolution control simulation with time‐invariant 1950 forcings following the HighResMIP protocol. In terms of global root‐mean‐squared error metrics, this high‐resolution simulation is generally superior to results from the low‐resolution configuration of E3SMv1 (due to resolution, tuning changes, and possibly initialization procedure) and compares favorably to models in the CMIP5 ensemble. Ocean and sea ice simulation is particularly improved, due to better resolution of bathymetry, the ability to capture more variability and extremes in winds and currents, and the ability to resolve mesoscale ocean eddies. The largest improvement in this regard is an ice‐free Labrador Sea, which is a major problem at low resolution. Interestingly, several features found to improve with resolution in previous studies are insensitive to resolution or even degrade in E3SMv1. Most notable in this regard are warm bias and associated stratocumulus deficiency in eastern subtropical oceans and lack of improvement in El Niño. Another major finding of this study is that resolution increase had negligible impact on climate sensitivity (measured by net feedback determined through uniform +4K prescribed sea surface temperature increase) and aerosol sensitivity. Cloud response to resolution increase consisted of very minor decrease at all levels. Large‐scale patterns of precipitation bias were also relatively unaffected by grid spacing.

1 Introduction

Earth System Models (ESMs) are one of our most important tools available for predicting, reconstructing, and understanding climate change. They can be used both to provide projections of future climate and to test proposed physical explanations for climate system behavior. ESMs are created by coupling global models of the physical processes governing the temporal evolution of the atmosphere, ocean, land, and cryosphere. These components are typically modeled by breaking the planet into a large number of interdependent grid cells from the bottom of the ocean to nearly the top of the atmosphere.

Global coverage, combined with the need to run for centuries at a time, puts severe restrictions on ESM grid resolution. Coarse resolution is a major problem for ESMs because a lot of important climate processes occur on scales smaller than a typical ESM grid cell. Subgrid scale variations must be parameterized in terms of grid‐cell averages. Clouds, for example, are often smaller than the O(100 km) grid spacing typical for GCMs in 2019 (Wood & Field, ). Ocean mesoscale eddies, a critical mechanism for heat advection, occur at scales of less than 10 km in high latitudes (Maslowski et al., ). This is far smaller than the ∼50‐km resolution often used at high latitudes to generate ensembles of 21st century Arctic sea ice loss (Jahn et al., ; Snape & Forster, ). These are just a few examples of resolution‐dependent mechanisms that remain problematic in most fully coupled simulations of Earth's climate.

Recognition of the uncertainties caused by coarse resolution has prompted many studies evaluating the impact of increased spatial resolution. For the atmosphere, resolution sensitivity of many aspects of climate seems to be model dependent. Some commonalities do emerge, however. Precipitation, snowpack, and streamflow in hilly or mountainous terrain improves with increasing horizontal resolution (Delworth et al., ; Duffy et al., ; Iorio et al., ; Pope & Stratton, ) due to better resolution of the terrain itself. Extreme precipitation also improves with increasing horizontal resolution (Iorio et al., ; Kiehl & Williamson, ; Terai et al., ; Wehner et al., , ) because smaller grid cells support more concentrated vertical motions and therefore larger condensation and fallout. The ratio of stratiform‐to‐convective precipitation also increases with horizontal resolution as more scales are explicitly resolved (Bacmeister et al., ; Demory et al., ; Hagemann et al., ; Hertwig et al., ; Pope & Stratton, ; Terai et al., ). In general, precipitation seems to improve more in winter than in summer because summertime rainfall tends to occur in the form of convective events that are smaller than even the high‐resolution grid cells commonly used for GCMs (Duffy et al., ). The impact of model resolution on clouds is more complicated. High resolution has been shown to reduce bias in tropical high cloud but increase bias in low cloud (Wang et al., ). Blocking of flow with resulting air stagnation problems is also better resolved as mountains and valleys come into focus (Berckmans et al., ; Jung et al., ; Matsueda et al., ; Schiemann et al., ). Sea breeze effects, which are important near the coasts in warm areas, also begin to be resolved as horizontal resolution improves (Boyle & Klein, ; Love et al., ). Because average eyewall radius for strong Atlantic hurricanes (Cat 2 or higher) is about 23 km (Knaff et al., ), resolution at least this fine is expected to be a necessary condition for accurately resolving tropical cyclones. Shaevitz et al. () confirmed this expectation in an analysis of 11 simulations from eight distinct climate models ranging in resolution from 130 to 28 km. They found that tropical cyclone intensity was only well captured (relative to observations) by the model with 28‐km resolution. Several other studies have affirmed the ability of ∼25‐km atmosphere models to reproduce the frequency and intensity of tropical storms (Atlas et al., ; Bacmeister et al., ; Wehner et al., ).

Ocean and sea ice model fidelity are particularly affected by the ability of a model to resolve eddies and deformation, respectively. The impact of permitting oceanic mesoscale eddies in fully coupled models is dramatic, affecting, among many things, the spread of near‐inertial energy (Zhai et al., ), water mass modification (Gnanadesikan & Hallberg, ), and modulation of sea ice extent and thickness (McPhee, ; McGeehan & Maslowski, ). Resolving the Gulf Stream improves the simulation of meridional heat transport in the Atlantic (Joyce & Zhang, ; Volkov et al., ), may increase the strength of the Atlantic Meridional Overturning Circulation (AMOC), and may improve simulated AMOC variability (e.g., Cheng et al., ; Hirschi et al., ; Marti et al., ; Stepanov & Haines, ). High resolution also helps to resolve fine‐scale sea ice deformation, which is important for simulating coupled fluxes against a cold atmosphere. Until recently, ocean models were most often coupled daily to atmospheric, land and sea ice components, filtering out the semidiurnal ice‐ocean response to storms, especially in the Southern Ocean (Roberts et al., ), and adversely affecting mixed layer depth in the high north (Holdsworth & Myers, ; Jin et al., ). Even with subinertial (∼hourly) coupling switched on, regions of intense ocean‐atmosphere flux exchange are adversely affected in coarse‐resolution models by horizontal diffusion in the wind stress curl (DuVivier & Cassano, ) and spreading of freshwater fluxes from land models (GarciaQuintana et al., ). Increased horizontal resolution improves upwelling of ocean water off the west coast of continents (Delworth et al., ; McClean et al., ; Small, Bacmeister, et al., ) due to improved surface wind stress from better‐resolved coastal topography (Gent et al., ). This increased upwelling is particularly important because it is expected to increase stratocumulus, the lack of which is a common feature in climate models. Improved atmospheric resolution enhances Antarctic katabatic winds (Lenaerts et al., ) and their forcing of a properly resolved slope current surrounding that continent, which require a minimum oceanic resolution of 0.5° (Mathiot et al., ). The benefits of high‐resolution bathymetry are most keenly felt where narrow passages are critical to water mass exchanges, such as through Nares Strait between Greenland and Ellesmere Island, which can have a significant sea ice flux (Kwok et al., ) and is often too wide or completely closed off in coarse‐resolution ESMs. El Niño Southern Oscillation often also improves with resolution (Delworth et al., ; Sakamoto et al., ). However, we note that a key mechanism for model improvement of ENSO, tropical instability waves (Roberts et al., ), are already well simulated in the standard resolution configuration of E3SMv1 (Golaz et al., ).

While many aspects of climate do improve with increasing horizontal resolution, some seem to be robustly insensitive to resolution. Resolution does not seem to be a panacea for double ITCZ problems (Bacmeister et al., ; McClean et al., ), except perhaps in the Eastern Pacific where water advected by the trade winds has more reasonable SST due to improved upwelling (Delworth et al., ). Resolutions of 25 km in the atmosphere also do not seem to improve simulation of the Madden‐Julian Oscillation (Bacmeister et al., ; Jung et al., ).

It is worth noting that the above discussion focuses on global atmosphere models running in the hydrostatic regime and ocean models able to resolve only the largest eddies. A new class of global atmospheric models that explicitly resolve cloud systems is gaining popularity (Satoh et al., ). Because these models capture such different scales, it is likely that their resolution sensitivity is quite different from that described above. Although E3SM has plans to develop such a model, running these cloud‐resolving models is unfeasible currently for decadal to centennial climate scales and is therefore outside of the scope of the current analysis.

While past studies of the impact of changing horizontal resolution are plentiful, relatively little analysis has been made of the impact of vertical resolution. Parametric assumptions usually are entangled with the vertical grid, so changing vertical resolution convolves discretization error with parametric errors, oftentimes resulting in little payback in model skill from additional computational expense. This is why even though Lindzen and Fox‐Rabinovitz () provide a sound theoretical argument for high vertical resolution in GCMs, we know of no GCMs that actually use the theoretically appropriate resolution. Richter et al. () study the effect of increased vertical resolution in CAM5, finding that an increased model top provides better quasi‐biennial oscillations and improved seasonal cycle of temperature in midlatitude. Partly based on this work, E3SM has switched from 30 to 72 vertical layers and pushed the model top up to 60 km for all model resolutions. The impact of increased vertical resolution on E3SM is discussed in Rasch et al. () and Xie et al. ().

The E3SM project was launched in large part to tackle the challenge of high‐resolution (Bader et al., ). The high‐resolution coupled model configuration described here is a central part of that strategy and is a standard configuration of the v1 release. This paper is to serves as a reference for the coupled high‐resolution configuration to complement the standard resolution configuration of E3SMv1 described in Golaz et al. (). In addition, this paper assesses the impact of horizontal resolution on E3SMv1. While there have already been many papers on this latter topic, the conclusions for E3SMv1 are somewhat novel and therefore of broad interest. Section briefly describes the E3SMv1 model and all differences between its low‐ and high‐resolution configurations, as well as the simulations used in this study. This is followed by an in‐depth evaluation of the model's 1950 control simulation in section . Modes of temporal variability are discussed in section , followed by an analysis of the climatology of extreme weather events in section . This is followed by a short discussion of the high‐resolution model's sensitivity to anthropogenic aerosol and greenhouse gas changes in section . Finally, conclusions are presented in section .

2 Experimental Design

2.1 Model Description

The E3SMv1 Earth System Model is described in Golaz et al. (). The atmospheric component is based on the spectral‐element atmospheric dynamical core (Dennis et al., ) with 72 vertical levels. It features the Zhang‐McFarlane (ZM; Neale et al., ; Richter & Rasch, ; Zhang & McFarlane, ) deep convection scheme and the Cloud‐Layers Unified by Binormals (CLUBB; Golaz et al., ; Larson, ; Larson & Golaz, ) for macrophysics, turbulence, and shallow convection. Morrison‐Gettelman Version 2 (MG2; Gettelman & Morrison, ; Gettelman et al., ) microphysics, four‐mode Modal Aerosol Model (MAM4; Liu et al., ), and Rapid Radiative Transfer Model for general circulation models (RRTMG; Iacono et al., ; Mlawer et al., ) are also used. Additional details are provided in Xie et al. (), Rasch et al. (), Golaz et al. ().

Ocean and sea ice components in E3SMv1 are based on the Model for Prediction Across Scales (MPAS; Petersen et al., ; Ringler et al., ). For the high‐resolution simulations presented here, the Gent‐McWilliams (GM) mesoscale eddy parameterization is not needed and is disabled. For the low‐resolution simulation with high‐resolution tunings, the GM bolus diffusivity is set to a constant 1,800 m²/s (consistent with Golaz et al., ). Neither the high‐resolution or low‐resolution simulations include a parameterization of submesoscale eddy transport. The land model is a slightly revised version of that found in CESM1 (Hurrell et al., ), as described in Golaz et al. ().

The river transport component is the newly developed Model for Scale Adaptive River Transport (MOSART; Li et al., , ). The core of MOSART is the simplified form of the one‐dimensional Saint‐Venant equation, a kinematic routing method that provides explicit simulation of riverine hydrologic variables such as channel velocity, channel water depth, and water surface area. Simulating these variables makes incorporating water management, heat, and biogeochemical processes straightforward. MOSART was developed to work with both grid‐ and vector‐based representations (Li et al., ; Tesfa et al., ). In E3SMv1, the grid‐based representation of MOSART is adapted.

As noted in Golaz et al. (), the ocean model in E3SMv1 accounts for changes in heat content of water with temperature, but the atmosphere does not. Keeping track of the heat content of water in the atmosphere is challenging, so E3SMv1 instead applies an ad hoc correction where the global‐average difference between the temperature at which water evaporates from the ocean and returns back as rain or stream outflow is applied as a globally uniform sensible heat flux correction. This correction was found to have negligible effect on model climate and is described in more detail in Golaz et al. () Appendix A.

2.2 Model Resolution

In this section, we describe the grids used for the high‐resolution (HR) and low‐resolution (LR) configurations of E3SMv1. Details of the HR and LR grid configurations are described below and summarized in Table . To isolate the effect of resolution, we will focus particularly on a LR run using the same parameterization constants (tuning parameters) as HR; we will call this run LRtunedHR to differentiate it from the E3SMv1 standard resolution configuration (hereafter LRv1) described by Golaz et al. ().

Information About Grids Used in This Paper

Grid	∼Δx	# of columns	vertical levels
HR atm	25 km	777,602	72
HR land	25 km	777,602	15
HR ocean	8–16 km	3,693,225	80
HR sea ice	8–16 km	3,693,225	5
HR river	0.125°	4,147,200	1
LR atm	110 km	48,602	72
LR land	110 km	48,602	15
LR ocean	30–60 km	235,160	60
LR sea ice	30–60 km	235,160	5
LR river	0.5°	259,200	1

Note. Grid spacing is given in km for near equal‐distance grids and in degrees for equal‐angle grids. For the sea ice model, the number of vertical levels refers to the number of thickness distribution categories.

At all resolutions, E3SMv1 uses one grid for atmosphere and land components, another grid for ocean and sea ice, and a third grid for streamflow. The atmosphere/land grids are based on a cubed‐sphere topology. Atmospheric parameterizations and the land model are called on each of 16 Gauss‐Lobatto nodes for each spectral element. After accounting for cells on element edges being shared between elements, this results in atmospheric and land processes being called on approximately nine times more columns than used for dynamics. The high‐resolution configuration of E3SMv1 contains 120 quadrilateral spectral elements in both x and y directions of each face of the cube sphere for a total of ∼800,000 columns corresponding to an approximate grid spacing of 25 km. The LR atmosphere, in comparison, has 30² total elements on each side of its cube for a total of ∼50,000 columns (110‐km nominal grid spacing). At all resolutions, E3SMv1 contains 72 vertical layers on a sigma‐pressure coordinate system. As described in Xie et al. (), layers are unevenly spaced with finer resolution near the Earth's surface. Atmosphere and land resolution of  25 km was chosen as a compromise between computational feasibility and improved representation of topography; it is not expected to explicitly capture turbulence or convection.

The MPAS ocean and sea ice models share the same horizontal grid, which at high resolution contains 3.7 million horizontal cells. In the HR configuration, ocean resolution is coarsest in the tropics and becomes finer near the poles. The grid is defined to well resolve the first Rossby Radius of deformation (Chelton et al., ) across most of the globe. This grid resolution is broadly similar with other high‐resolution ocean configurations (e.g., Delworth et al., ; McClean et al., ; Small, Bacmeister, et al., ). The LR MPAS mesh was designed such that the horizontal resolution has a similar meridional distribution (on average) as the standard CESM1 ocean grid, varying from 60 km in the midlatitudes to 30 km near the equator and poles. The high‐resolution version of the E3SMv1 ocean model uses 80 levels in the vertical, with spacings following Stewart et al. (). This vertical grid is designed to resolve the vertical structure of mesoscale eddies. The LR ocean model vertical discretization is very similar to CESM1, with 60 layers varying from 10 m spacing at the surface to 250 m at depth. The sea ice model includes five thickness categories in each grid cell divided at 0.65, 1.39, 2.47, and 4.56 m.

For coupled simulations, MOSART employs a uniform latitude/longitude grid with grid spacing listed in Table . Only grid cells over land are active; there are 1.5 million active cells at high resolution and 0.1 million at low resolution. MOSART is not discretized in the vertical.

Finer grids typically require shorter time steps for stability. This is certainly true for E3SMv1, as evident from the time step information in Table . Simultaneous modification of grid spacing and time step can result in model sensitivity to time step being erroneously attributed to spatial resolution sensitivity (Jung et al., ). Due to time pressures, we do not decouple time step and resolution sensitivity in this study and instead use LR time steps for the LRtunedHR simulation. As a result, what we describe as resolution sensitivity should be interpreted as including all changes required to run at that resolution.

Time Steps (in Min) Used in High‐Resolution and Low‐Resolution Configurations

Time step	HR	LR
atm radiation	60	60
atm phys/dyn coupling	15	30
atm CLUBB+microphysics	5	5
atm dyn remap	2.5	15
atm advection and dyn	1.25	5
atm hyperviscosity	0.417	1.67
ocn	6	10
ocn barotropic	0.2	0.67
ice thermodynamics	15	15
ice dynamics	7.5	15
river	60	60
river coupling	180	180
atm/ice/lnd coupling	15	30
ocn coupling	30	30

2.3 Tuning

The exact values for many parameters in climate models are not precisely known. Values for these parameters are generally chosen to minimize model bias in reproducing current climate or to achieve global top‐of‐atmosphere (TOA) radiation balance for preindustrial conditions. Because empirically chosen tuning parameters can compensate for imperfect process representations, they play a critical role in model development (Hourdin et al., ; Schmidt et al., ). Tuning involves testing a variety of parameter settings and choosing the best configuration. Thorough tuning often requires hundreds of tests, which quickly becomes prohibitively expensive at high resolution. While parameters chosen for LR simulations often behave reasonably well at higher resolution, optimal tuning choices are generally resolution dependent. This is because parameters are used to compensate for unresolved physics, and more physical processes become resolved at higher resolution.

Because of computational expense, tuning of high‐resolution configurations is often accomplished by changing only one or two parameters to remove the most egregious model biases. As a result, it is sometimes hard to disentangle whether lack of improvement at high resolution is due to unimportance of resolution or poor tuning (Bacmeister et al., ). One feature of the current study is that substantial effort was put into tuning the high‐resolution atmosphere component, with the result (shown later) that the high‐resolution version of E3SMv1 has generally lower bias than its low‐resolution counterpart. In order to tune the high‐resolution model in a computationally efficient manner, the Cloud‐Associated Parameterization Testbed (CAPT; Ma, Chuang, et al., ) was extensively employed. CAPT takes advantage of the fact that many atmospheric biases show up in 48 hr forecasts, which are cheaper to run and can be compared against observed weather. Tuning started from the LR parameter settings and targeted mainly the global‐mean TOA radiative budget and surface precipitation. Interestingly, the biases that showed up as a result of increasing horizontal resolution were similar to the biases resulting from increasing the number of vertical layers in the LR configuration. This allowed us to use similar tuning strategies for both problems. These strategies are described in more detail in Xie et al. ().

Another major difference between HR and LRv1 versions of E3SMv1 is in their treatment of heterogeneous ice nucleation. The LRv1 version of E3SMv1 employs a new parameterization of this process based on Classical Nucleation Theory (CNT; Wang et al., ). CNT depends explicitly on aerosol properties (e.g., nucleation efficiency and surface area) and computes ice nucleation separately for immersion, deposition, condensation, and contact modes. CNT replaces the Meyers et al. () parameterization in CESM1. The Meyers scheme only depends on ice supersaturation for liquid‐saturated conditions (and therefore is uniquely specified by temperature). The Meyers scheme was derived from Northern Hemisphere midlatitude measurements. Because aerosol concentrations are lower in polar regions than in Northern Hemisphere midlatitudes, the Meyers scheme is expected to overestimate number concentrations of ice crystals at high latitudes. In contrast, CNT is likely to underestimate ice crystal number at high latitudes because E3SMv1 underestimates aerosol transport to high latitudes (shown later). At high resolution, ice number concentration using the CNT scheme became extremely low, resulting in an unrealistic blanket of supercooled liquid over polar regions during winter. This blanket warmed the wintertime polar regions by more than 10 K, wrecking the model's sea ice distribution. To avoid this problem, we reverted the heterogeneous ice nucleation scheme back to the Meyers scheme as part of the HR tuning. We later realized, however, that CNT was behaving poorly because the rate at which liquid was being converted to ice by the Bergeron process was being scaled by a factor of 0.1 to compensate for the Meyers scheme's propensity to create too many ice crystals and therefore too much ice. Thus, CNT behaved poorly not because of a design flaw, but because it disturbed a balance between compensating errors. In retrospect, a better approach would have been to keep CNT and scale Bergeron to a reasonable value. We will pursue this option in E3SM Version 2.

A downside to our strategy of intensive tuning at high resolution is that our LR and HR configurations differ due to a mixture of resolution and tuning choices. Differences in tuning parameters are listed in Table . It is also worth noting that a few minor bugs were found and final coupled tuning was done between the time the atmosphere‐component papers from Zhang et al. (), Tang et al. (), Xie et al. (), and Rasch et al. () were published and the coupled simulations documented in Golaz et al. () were completed. As a result, LRv1 values of ke, c14, and so4_sz_thresh_icenuc published here differ from those in Xie et al. (). For the most part, these changes did not impact model biases. One exception is an unrealistic lack of high clouds in the Tropical West Pacific in earlier atmosphere‐component papers, which is greatly improved here and in the fixed‐SST simulations shown by Golaz et al. ().

List of Tuning Parameters That Differ Between Low‐ and High‐Resolution Versions of E3SMv1

Parameter	HR	LRv1	Parameterization	Meaning (units)
alfa	0.2	0.1	ZM	ratio of upward versus downward mass flux at downdraft initiation level (‐)
c0_lnd	3.5e–3	7e–3	ZM	Precip efficiency over land (m−1)
c0_ocn	4.3e–3	7e–3	ZM	Precip efficiency over ocean (m−1)
dmpdz	−0.2e–3	−0.7e–3	ZM	parcel fractional mass entrainment rate (m−1)
ke	6e–6	5e–6	ZM	Precip evaporation efficiency ((kg m−2 s−1)−1/2 s−1)
dp1	3.9e–2	4.5e–2	ZM	Deep convective cloud fraction (−)
c1	1.5	1.335	CLUBB	$\overline{w^{\prime 2}}$ dissipation strength (−)
c8	4.73	4.3	CLUBB	$\overline{w^{\prime 3}}$ damping strength (−)
c14	1.75	1.06	CLUBB	Strength of $\overline{u^{\prime 2}}$ and $\overline{v^{\prime 2}}$ dissipation (−)
dust_emis_fact	2.5	2.05	Dust emission	Inverse of dust emission strength (−)

Several more bugs were found and corrected after LRv1 was completed and before these HR simulations started. In particular, a bug was found and corrected in precomputed values of the rrtmg_sw_reftra table lookup (https://github.com/E3SM-Project/E3SM/pull/2202). Another bug related to energy errors in land‐surface water‐phase change and lake/snow interactions (https://github.com/ESCOMP/ctsm/pull/307) was also fixed. Both of these bugs had little effect on model climate.

Unlike the low‐resolution configuration, the coupled version of the HR model uses identical tunings to the fixed‐SST HR simulations described in Zhang et al. (), Tang et al. (), Xie et al. (), and Rasch et al. (). This decision was made because the coupled model seemed to be in relatively good radiative balance without further coupled tuning and time pressures required us to start our simulation campaign. In addition, the efficacy of our tuning parameters seemed to saturate at the values currently used for HR such that further tuning provided little benefit. Another tool for tuning E3SM at high resolution is the ability to use a variable‐resolution mesh in the atmosphere (Tang et al., ). While we did not use this capability much for v1 tuning, we expect this technology to play an important role in the development of future model versions.

2.4 Description of Simulations

Aside from LRv1, all coupled E3SMv1 simulations described in this paper use time‐invariant 1950 forcings and follow the HighResMIP protocol of Haarsma et al. () as closely as possible. Using 1950 conditions instead of the preindustrial conditions typically used for low‐resolution control runs is necessary in order to extend a control simulation into a transient 20th century/future climate simulation—simulating the period between 1850 and 1950 would exhaust our computing allocation, preventing the eventual simulation of the more interesting 1950–2050 period. Although the planet was not strictly in radiative equilibrium in 1950 (as assumed in a control simulation), most anthropogenic greenhouse gas (GHG) emissions and warming have occurred since 1950 (e.g. IPCC, , Fig. SPM.1), so this approach is a reasonable compromise. Using anthropocene‐era control simulations for high‐resolution control runs is a common practice (Delworth et al., ; Sakamoto et al., ).

Forcings for the 1950 coupled simulations are identical to the forcings used for year 1950 of the LRv1 historical simulation (Golaz et al., ) with a few exceptions. The mean of solar insolation from 1939 through 1961 was used in order to average over two periods of the 11‐year solar cycle. Tropospheric aerosol emissions were handled similarly in HR to LR but use higher resolution 0.5° data from input4MIPS. Stratospheric volcanic aerosols are averaged over 1850–2014 to capture the conditions that would be observed if simulating a long period of time. Tropospheric oxidants are specified using 1955 values because data for 1950 are not available and 1955 values should be negligibly different. Dimethyl sulfide is computed by interpolating 1849 and 2000 values to the period of interest; this treatment is identical to Golaz et al. () but was not mentioned in the previous paper. Anthropogenic and fire emissions for aerosols and precursor gases are taken from the CMIP6 emissions data set. One exception to this is secondary organic aerosol (SOA) gas emissions, which are handled specially as described by Rasch et al. () in order to improve the vertical distribution of SOA.

While the focus of this paper is on the high‐resolution configuration of E3SM, it is occasionally useful to compare against the E3SMv1 LR release (LRv1) configuration. In particular, we will occasionally examine the preindustrial (PI) control and the 20th century historical simulations described Golaz et al. ().

To provide crude estimates of model sensitivity to anthropogenic GHG and aerosol changes, we also perform short simulations with prescribed sea surface temperature (SST) and sea ice extent. To estimate the net climate feedback parameter, we perform a pair of 5‐year current‐climate simulations with annually repeating forcings. SST is uniformly raised by 4 K in one simulation and feedback strength is inferred from the difference in resulting TOA radiative imbalance between runs. Aerosol sensitivity is computed from a pair of 15‐month simulations nudged to observations with a 6‐hr relaxation timescale. In one simulation, anthropogenic aerosols are omitted. In the other, they are included. Detailed explanation of these sensitivity studies and their results are provided in section .

For all prescribed‐SST runs, forcings representative of 2010 are used unless explicitly mentioned. As recommended for HighResMIP, these runs use daily‐resolution prescribed SST and sea ice extent at 25‐km resolution from HadISST2 (Haarsma et al., ). GHG and volcanic forcings are taken from year 2010 of the CMIP6 Historical experiment for these runs. Solar forcing is the average of CMIP6 Historical data for the 22 years centered on 2010.

2.5 Initialization

The initial ocean temperature and salinity distributions were interpolated from the polar hydrographic climatology (Steele et al., ). To generate currents and propagate initial gravity waves, MPAS‐Ocean was run in standalone mode for three months, restoring temperature and salinity to the initial climatology and imposing an annually averaged wind stress. The results were used to drive a simulation with interactive ocean and ice components and atmospheric forcing prescribed from the COREv2 interannually varying (Large & Yeager, ) dataset. Ice was originally set to a uniform thickness of 1 m at all locations poleward of 70° for this calculation. Only years 1947–1950 were simulated in this mode. This 3‐year spin‐up period is a compromise between minimizing the amount of model drift inherited by the fully coupled simulation and preventing large shocks upon coupling. The resulting ocean and ice state was combined with an atmosphere initial condition from an earlier experimental HR simulation with fixed SST and a land initial condition interpolated from year 1950 of a LR precursor to the E3SMv1 CMIP LRv1. This configuration was run for 5 years before realizing that polar wintertime surface temperatures were unacceptable and deciding to revert from the CNT mixed‐phase nucleation scheme to the old Meyers et al. () scheme to fix the problem (as described in section ). Instead of discarding those 5 years of simulation with the CNT scheme, we decided to begin our production off the end of the 5‐year CNT run in order to start the model from a slightly more equilibrated state. For consistency, our LR simulation with HR tuning follows an identical spin‐up protocol, including 5 years of coupled simulation with the CNT mixed‐phase nucleation scheme. For both simulations, we denote year 1 as the first year with the Meyers rather than CNT nucleation scheme.

2.6 Performance

The v1 HR simulation was performed at DOE's Argonne Leadership Computing Facility (ALCF) Theta supercomputer, which has 4,400 nodes of Intel KNL (“Knights Landing”) with 64 1.3 GHz cores and 192 GB memory per node. The coupled E3SM model consists of the five components described in section , which exchange fluxes through a top‐level coupling driver (CPL). Overall model throughput is optimized by laying out components across processing cores such that none of the components is a throughput limiter, while at the same time none of the processing cores are idle. Figure shows the layout used for most of our v1‐HR simulation. E3SM's component synchronization approach (RASM_OPTION1) prevents the atmosphere model from running in parallel with sea ice, land, or river runoff models, but allows the ocean to run in parallel on its own computational cores. As evident in this figure, the atmosphere and sea ice models are the main determinants of total run time, with the number of ocean‐model cores chosen to fit within the time taken by these other processes. Inability of the sea ice model to scale to large process counts, combined with its inability to run in parallel with the atmosphere model, results in unavoidable idle time (shown as white space at the center‐bottom of Figure ). CPL sequences, interpolates, and exchanges data among all components and consumes about 20% of overall run time. CPL scaling limitations result in the idle time indicated by white space in the top‐right corner of Figure . Overall throughput for this layout was ∼0.7 simulated years per wall‐clock day. Various performance optimizations (reduced synchronization and improved threading) are being added to ATM and ICE to raise throughput in v2. Using more processor cores would have improved simulation speed but would have been more expensive (because of increased communication lag time) and would have resulted in more time waiting in queue. The 800‐node layout shown was chosen because it is the smallest layout that qualified for the largest production queue and its associated 24‐hr wall‐clock limit.

Typical processor layout and throughput for a 6‐month simulation on Theta. This layout uses 32 MPI processors and 4 OpenMP threads per process for a total of 128 hardware threads per node. In this figure, the atmosphere models is called ATM, the ocean model is OCN, the sea ice model is ICE, the land model is LND, and the river model is the small red bar at the center‐bottom of the plot.

LRtunedHR was run using 54 nodes of Argonne's Anvil cluster, which has 36 64‐GB Intel Xeon Broadwell cores per nodes. It achieved approximately five simulated years per day (SYPD) and consumed ∼9,000 core hours per simulated year. LR runs with prescribed SST were run on 169 nodes of NERSC Cori Haswell nodes, which have two 2.3‐GHz 16‐core Haswell processors per node. It ran at ∼13 SYPD and cost ∼10,000 core hours per simulated year. Prescribed‐SST HR simulations completed ∼0.80 SYPD using ∼640,000 core hours per simulated year on 323 NERSC Cori KNL nodes, which have one 1.4‐GHz 68‐core Intel Xeon Phi Processor 7250 per node.

The HR configuration contains ∼16× more grid points than LR (Table ) and requires a time step which is overall ∼2× shorter (Table ). As a result, we expect (and find) HR runs to be ∼32 times more expensive than LR on a given machine. Resolution increases have less effect on time to solution because of greater opportunities for parallelism.

Source code for the simulations has been taken from a long‐term maintenance branch maint‐1.0(https://github.com/E3SM-Project/E3SM/tree/maint-1.0) that will be maintained to ensure that runs can be exactly (bit for bit) reproduced if run on the same model with the same software environment.

3 Model Climate

3.1 Global‐Average Evaluation

Figure displays the time series over the simulation for several key large‐scale climate parameters. While the LRv1 and LRtunedHR runs are close to TOA radiation balance, the HR run loses energy at the rate of ∼0.5 W/². This slow energy loss rate results in a slight cooling trend to bring the globally averaged surface air temperature closer to the 20th century observed value, whereas the LRtunedHR run is over 0.5° too cold. The small imbalance and slight global trend in our HR run is similar in magnitude to the imbalances found in the first several decades of many fully coupled high‐resolution climate simulations. HR simulations tend to have larger radiative imbalances and corresponding temperature trends than their low‐resolution counterparts because they are too expensive to optimize in a fine‐tuned way. It is recognized in the HighResMip protocol that the 50‐year spin‐up period is a compromise between computer expense and a drift‐free model. The near‐zero energy imbalance in the LRtunedHR simulation was a pleasant surprise, as it took many iterations to achieve an acceptable, drift‐free climate for the standard resolution LRv1 run (Golaz et al., ). Trends in sea ice volume are shown in the lower panels of . In both hemispheres, ice volume is relatively steady across the simulation, indicating that evaluating climatological averages is appropriate.

Time series of key global‐average quantities. Zero TOA net energy balance is depicted as a dashed line in the top panel. In the second panel, horizontal dashed line at 13.9°C is an estimate observed global average temperature over the 20th century (NOAA Global Climate Report, 2018; https://www.ncdc.noaa.gov/sotc/global/201813).

Equilibration of the ocean is shown in Figure . Consistent with the overall loss of energy of the planet, the global ocean is cooling in HR over most of the water column, but especially in the upper 1,000 m. The deeper ocean remains fairly isolated and cools with time. HR is accumulating a small fresh bias in the upper ocean, but this bias is insignificant compared to the fresh bias in the upper 800 m of the water column of LRtunedHR. A separate study is now underway to investigate this accumulating fresh bias at the ocean surface in v1 low‐resolution simulations.

Trends of globally averaged ocean temperature (left panels) and salinity (right panels) versus depth for (a and b) HR and (c and d) LRtunedHR. Trends are computed with respect to annual mean of the first simulated year.

Because of the better simulation of scale interactions, increased resolution is widely assumed to result in better accuracy at scales resolved by lower‐resolution models. We test that assumption by comparing large‐scale measures of model bias among the LRv1, LRtunedHR, and HR model runs.

To get a crude sense of HR bias relative to LRtunedHR, LRv1, and models developed by other groups, Figure provides uncentered global root‐mean‐squared error (RMSE) values for a variety of climatically important variables for all 3 E3SMv1 models and all available CMIP5 models. Prior to analysis, output from all models was spatially filtered through area‐average regridding to 2.5°. We chose uncentered RMSE to account for global‐mean bias in the simulations. As was seen in similar comparisons of the ECHAM atmospheric model by Hertwig et al. () and CESM by Small, Tomas, et al. (), HR has lower error than LRv1 for most variables and seasons. In particular, HR is never in the lowest 25% of CMIP5 models and is among the best CMIP5 models for most variables. Because HR differs from LRv1 due to tuning as well as resolution differences while LRtunedHR has not been tuned to reduce bias at all, it is inappropriate to conclude that HR's superior climatology is due solely to resolution changes. In addition, CMIP5 runs had completed many hundreds of years of simulation by the time we sampled their errors, while HR's spin‐up period was only a few decades. While short spin‐up is typical (and perhaps unavoidable) for high‐resolution simulations, this discrepancy may give high‐resolution simulations an unfair advantage in evaluations of this type. Assessing the impact of initialization approaches on HR evaluation is the subject of important future work.

Comparison of uncentered spatial RMSE (computed following Gleckler et al., ) from years 1981–2005 of 45 CMIP5 historical simulations (box and whiskers with boxes showing 25th–75th percentile values and whiskers showing the span of the data) versus values from years 21–50 of 1950 control runs for HR (cyan) and LRtunedHR (green). Dark blue dots show RMSE values for 1981–2005 of the first member of the E3SMv1 LRv1 historical ensemble. Error for each variable is calculated relative to the relevant observations: CERES‐EBAF for radiation (Loeb et al., ), GPCP for precipitation (Adler et al., ), and ERA‐Interim for the remaining variables (Dee et al., ).

3.2 Clouds and Radiation

Spatial gradients in the amount of radiation absorbed or emitted by the Earth system drive the circulations in both the ocean and the atmosphere, control local temperature and humidity structure, and power the hydrological cycle. As noted above, HR is losing heat at a globally average rate of ∼0.5 W/m² with significant interannual variability. In nature and models, the globally averaged net energy imbalance is the small residual of the highly variable regional TOA net energy fluxes typically 1–2 orders of magnitude larger. As shown in Figure , regional TOA net energy biases are fairly small except for positive bias (too much absorption) in the eastern subtropical ocean basins, over Eastern China, in the Atlantic sector near Antarctica, and (to a lesser extent) over Africa. These are compensated by milder but more widespread negative biases along the west coasts of continents, over the East/Central tropical Pacific (probably associated with intertropical convergence zone biases described later), over the Indian Ocean, and at high northern latitudes. The regional‐scale high‐ and low‐resolution radiative imbalances in both resolutions of E3SMv1 look extremely similar (compare Fig.  with Figure 4 from Golaz et al., ). LRtunedHR also shows very similar pattern (not shown), indicating that this net radiative bias structure is a robust feature of E3SMv1 physics.

TOA radiative imbalance (in W/m2), positive warms planet). Panel a shows the annual‐average climatological distribution from CERES‐EBAFv4.0, while panel b shows bias in HR.

Clouds are the main source of bias in net radiation. This is evident from the fact that bias in shortwave (SW) cloud forcing (TOA SW radiation minus clear‐sky SW radiation, hereafter SWCF) for v1 HR in Figure a has a very similar bias pattern to net radiation in Figure b. From these figures, we conclude that excessive absorption in the eastern subtropical oceans is due to a lack of SW cloud reflection. Longwave (LW) cloud forcing (LWCF) in v1 HR is too weak almost everywhere except at high latitudes. Eliminating the LW bias through tuning proved difficult, forcing us to compensate with a smaller globally averaged SWCF to maintain global‐average net radiative balance. This is yet another example of the limitations of tuning. The problem here is an inability to capture transitions between cloud types such that increasing stratocumulus in the eastern subtropical oceans also increases low clouds in convective regions, exacerbating negative SWCF biases in those regions.

Panels a–d show annual‐average SWCF and LWCF bias in SWCF (left) and LWCF (right) relative to CERES‐EBAFv4.0. Panels e–f show differences between HR and LRtunedHR simulations. Positive SWCF values mean less reflection by clouds, while positive LWCF values mean clouds are trapping more LW radiation. Units for all panels are W/m2.

While net radiation looked very similar between HR and LRtunedHR, more marked differences appear in SWCF and LWCF. In particular, Figures e and f show that finer resolution generally results in SWCF becoming more positive and LWCF becoming less positive. This is fairly uniform except at high latitudes where cloudiness seems to increase. Both of these tendencies are consistent with a reduction in cloudiness with increasing resolution, which is apparent in plots of total cloud fraction (Figure ). This reduction of cloudiness with increasing resolution was also noted in runs with fixed SST and sea ice (Xie et al., ; Zhang et al., ) as well as in other models (Bacmeister et al., ). Reductions in both SWCF and LWCF indicate a reduction of high clouds with resolution, while the ∼1 W m⁻² larger reduction in SWCF versus LWCF indicates reduction from low clouds as well. It is worth noting that although clouds are sensitive to resolution in E3SMv1, resolution sensitivity is greatly decreased relative to many earlier generation models (Bacmeister et al., ).

Comparison of annual mean total cloud fraction between E3SM and Calipso observations. (a) Observed (OBS) total cloud fraction as retrieved from the Calipso satellite with the GOCCP v3.1.2 algorithms (Chepfer et al., ) for the years 2006–2017. (b) Difference between HR and Calipso observations (model minus observed) and (c) the difference between LRtunedHR and observations. Model fields are those after interpretation with the use of a lidar simulator (Chepfer et al., ) and are taken from simulated years 41–50. Number above a indicates the global mean total cloud fraction, whereas numbers above b and c indicate the mean bias (model minus observed) and the root‐mean‐square error (in parenthesis).

Figure reveals that E3SM underestimates cloudiness at low latitudes regardless of resolution. This finding is particularly true in the stratocumulus areas of the eastern subtropical oceans, but also along the equator and over low‐latitude land masses. The marine stratocumulus errors are collocated with large SST errors (Figure ), particularly to the west of Baja California and Namibia/Angola. E3SM atmosphere‐only simulations with observed SSTs (not shown) have very similar but somewhat smaller errors, demonstrating that the stratocumulus error, while originating in the atmosphere model, is amplified by the SST error in the coupled model, similar to the E3SM low‐resolution LRv1 simulations (Golaz et al., ). Other cloud characteristics such as the cloud altitude and optical depth distributions of this E3SM high‐resolution model are quite similar to those found in atmosphere‐only E3SM simulations including those analyzed in Zhang et al. (). An exception is a moderate reduction in the amount of super‐cooled liquid in this high‐resolution model (not shown) due to the use of CNT rather than Meyers ice nucleation at high resolution (see section ). While we have emphasized model errors, it is worth bearing in mind that E3SM's simulated clouds compare favorably to those of other models. Scalar metrics of model performance from Klein et al. () indicate that E3SM would be in the top rank of pre‐CMIP6 models (Zhang et al., ).

To appreciate the value of higher resolution in the climatological distribution of clouds requires closer examination of the fine‐scale structure of simulated fields. Figure shows this structure over South America with model fields plotted close to their native resolution after interpolation to a latitude‐longitude grid. Despite the large‐scale biases present (such as the general underestimate of cloud over South America and in the marine stratocumulus region west of the continent), much of the regional structure—particularly near mountain ranges and coastal boundaries—is improved. In particular, v1 HR is able to capture the sharp gradient between cloudy and cloud‐free conditions off the west coast. The Atacama desert in particular shows up as an area with extremely low cloud fraction, which jumps discontinuously to much larger cloudiness just offshore. Another improvement is the more narrow cross‐width to the northwest‐to‐southeast feature immediately northeast of the Cordillera Oriental portion of the Andes in Peru and Bolivia.

Comparison of annual mean total cloud fraction between E3SM and observations over South America. Total cloud fraction is shown from (a) ISCCP HGM data (Rossow et al., ) for the years 1983–2012 at 1° resolution, (b) the HR at 0.25° resolution, and (c) the LRtunedHR at 1° resolution. Model output is for simulated years 41–50, and modeled total cloud fraction is computed from the ISCCP satellite simulator (Klein and Jakob, ) using optical depths >0.3.

3.3 Aerosols

Because aerosol observations are lacking for the 1950s and aerosol concentrations have changed greatly since then, this section focuses mainly on sensitivity of simulated aerosol life cycle to horizontal resolution. This is one of the first studies that investigates both the individual and total aerosol responses to resolution change, as previous studies have focused on either dust (Ridley et al., ) or carbonaceous aerosols (Ma et al., ; Liu et al., ). We also examine the impact of model tuning (LRtunedHR vs. LRv1) on aerosol simulations.

Table summarizes the global estimates of aerosol budgets. The annual‐mean aerosol optical depth (AOD) is estimated to be 0.123 for the HR run, nearly two thirds of which is contributed by naturally emitted dust and sea‐spray aerosols. In the LRtunedHR simulations, coarser model resolution leads to smaller dust and sea‐spray aerosol emissions, a consequence of the dependence of these emissions on grid‐mean surface wind speeds. The HR model produces more frequent occurrences of high winds due to resolution of fine‐scale topography and dynamics (e.g., Ridley et al., ). Zhang et al. () have developed a parameterization that accounts for the spatial variability in the subgrid winds for dust and sea‐spray aerosol emissions. The impact of applying such a scheme will be explored in future E3SM model development.

Global Aerosol Budgets

Characteristic	HR	LRtunedHR	LRv1	Expected PD value
Global AOD	0.123	0.116	0.14
Dust AOD	0.032	0.025	0.026
Emission (Tg/year)	4,577.	3,881.	4,702.	1,840.±49%
Dry depo (Tg/year)	3,309. (0.33)	3,053. (0.42)	3,653. (0.43)
Wet depo (Tg/year)	1,268. (0.13)	983. (0.12)	1,048. (0.12)
Burden (Tg)	27.2	19.8	23.5	19.2 ±40%
Lifetime (days)	2.17	1.86	1.82	4.14 ±43%
Sea spray AOD	0.046	0.042	0.049
Emission (Tg/year)	3,728.	3,434.	3,601.	16,600. ±199%
Dry depo (Tg/year)	2,320. (1.15)	2,451. (1.46)	2,551. (1.39)
Wet depo (Tg/year)	1,404. (0.70)	983. (0.59)	1,050. (0.57)
Burden (Tg)	5.52	4.60	5.04	7.52 ±54%
Lifetime (days)	0.54	0.49	0.51	0.48 ±58%
Sulfate AOD	0.013	0.014	0.023
Dry depo (Tg/year)	15.9 (0.035)	17.1 (0.036)	27.0 (0.035)
Wet depo (Tg/year)	78.0 (0.173)	78.9 (0.165)	117.2 (0.151)
Burden (Tg)	1.24	1.31	2.06	2.02 ±25%
Lifetime (days)	4.8	5.0	5.39	4.12 ±18%
BC AOD	0.002	0.003	0.005
Dry depo (Tg/year)	1.47 (0.059)	1.69 (0.064)	2.95 (0.056)
Wet depo (Tg/year)	2.84 (0.114)	2.60 (0.099)	4.78 (0.090)
Burden (Tg)	0.068	0.072	0.145	0.24 ±42%
Lifetime (days)	5.8	6.13	6.85	7.12 ±33%
POM AOD	0.005	0.005	0.01
Dry depo (Tg/year)	9.41 (0.047)	10.5 (0.049)	16.3 (0.038)
Wet depo (Tg/year)	21.4 (0.107)	20.2 (0.095)	33.7 (0.079)
Burden (Tg)	0.55	0.58	1.17	1.7 ±27%
Lifetime (days)	6.51	6.9	8.5	6.54 ±27%
SOA AOD	0.025	0.027	0.025
Dry depo (Tg/year)	7.32 (0.009)	7.72 (0.009)	9.83 (0.011)
Wet depo (Tg/year)	61.0 (0.073)	60.5 (0.068)	65.6 (0.076)
Burden (Tg)	2.3	2.46	2.36	0.57 ±117%
Lifetime (days)	12.3	13.2	11.4	6.70 ±115%

Note. The numbers in parentheses are calculated dry (or wet) removal (or deposition rates) in units of day⁻¹, defined as dry (or wet) deposition/(burden*365.) in Textor et al. (). Also shown are the model outputs from LRv1 for simulations with year 2000 aerosol forcing. The means and normalized standard deviations of the Expected PD values are taken from other modeling studies as collected by Liu et al. ().

Increasing resolution also modifies aerosol removal (i.e., deposition) and long‐range transport. To isolate the resolution effect on individual aerosol dry or wet removal pathways, aerosol removal rate, defined as the inverse of the aerosol lifetime (in unit of day⁻¹) for a particular removal pathway (Textor et al., ), is calculated and included in Table (numbers in parentheses). By comparing removal rates rather than the actual values of deposition fluxes, complications due to the different aerosol emissions and burdens in the two simulations are avoided. As a result of better‐resolved subgrid turbulent transport in HR, the dry deposition rates for dust and sea‐spray aerosols are smaller in HR than in LRtunedHR. This result is consistent with previous findings for carbonaceous aerosol (Liu et al., ). If aerosols are removed predominantly by dry deposition, this reduction would result in an increase of aerosol burden and residence time at higher resolution.

In contrast, aerosol wet removal is enhanced in HR due to improved representation of spatial inhomogeneity in clouds and precipitation (Ma, Rasch, et al., ). This enhancement is larger for sea‐spray aerosols, which are more susceptible to wet scavenging than dust. Increased wet removal counteracts but does not overcome increased emissions and decreased dry deposition, with the net effect that increasing resolution leads to increased sea spray and dust burdens, AODs, and lifetime. Dust is more sensitive to increased resolution than sea‐spray aerosols, with the resulting larger enhancement in dust AOD (+28%) versus sea spray (+10%).

In comparison with natural aerosols, the resolution effects on anthropogenic aerosols (sulfate, black carbon (BC), particulate organic matter (POM), and SOA) are less pronounced because the emissions of the anthropogenic aerosol species or their precursors are inventory based and independent of resolution. Anthropogenic aerosols are, however, subject to the same resolution sensitivity of removal mechanisms as described previously for natural aerosols. This is visible in Table in the form of decreased dry deposition and increased wet deposition rates for anthropogenic aerosols. Note, however, that anthropogenic aerosols are removed preferentially by wet deposition (in terms of total mass) rather than dry deposition due to their finer particle sizes. As a result, the net effect of increasing resolution leads to smaller aerosol burden and AOD, and shorter lifetime for sulfate, BC, POM, and SOA.

In our simulations, enhancement in the natural aerosol AOD dominates the reduction in anthropogenic AOD as model resolution increases, resulting in a net increase in total AOD. The estimated global and annual mean total AOD is about 6% higher in HR than in LRtunedHR.

As expected from industrialization over the last 70 years, global AOD in LRv1 is higher than in both HR and LRtunedHR. As mentioned in the beginning of this section, it is not meaningful to compare total and anthropogenic AOD or burdens in LRv1 with HR or LRtunedHR. Aerosol removal rates and lifetime, however, are less dependent on the emission scenarios and therefore worthy of comparison. Differences between LRtunedHR and LRv1 reflect the impact of different tunings since both runs use the same grid spacing. It is interesting to note that the simulated dust and sea‐spray aerosol removal rates and lifetime are very similar between LRtunedHR and LRv1, suggesting that our tuning differences have little impact on aerosol removal and lifetime. Dust and sea‐spray emissions, on the other hand, were purposefully decreased during the HR tuning in order for their AODs to match expected real‐world conditions. This explains why the global emissions of dust and sea‐spray aerosols are similar between HR and the LRv1 despite the resolution differences. Because of the resolution effect on aerosol removal, however, lifetime of natural aerosols is longer in HR than LRv1. For anthropogenic aerosols, the HR simulation results in shorter aerosol lifetime than found in LRv1 or in LRtunedHR. This is the expected result of enhanced aerosol wet removal at higher resolution. Resolution‐induced changes in individual aerosol lifetime may lead to further changes in corresponding aerosol distributions and radiative forcings. These changes are beyond the scope of this overview paper but are explored for dust in Feng et al. (). Comparison with the expected present‐day values in the last column of Table indicates that the simulated aerosol lifetime of natural and anthropogenic aerosols in HR are within the uncertainty range of the AeroCom models (Textor et al., ) as summarized in Liu et al. ().

Resolution effects on aerosol simulations are much more remarkable on regional scales. In particular, the sign of the regional resolution effect changes depending on the predominant aerosol types. The geographical distributions of the annual mean AOD simulated by v1 HR and LRtunedHR are shown in Figures a and b. The aerosol climatology in both figures look as expected, with the lowest AOD in polar regions and highest AOD over the major dust source regions such as the Taklamakan and Gobi deserts in Asia and the Sahara desert in Africa. Compared to the LRtunedHR, HR predicts larger AOD over deserts and over midlatitude oceans, consistent with the stronger dust and sea‐salt emissions at higher resolution noted in Table (Figure c). The maximum increase in regional AOD is much larger than the 6% increase in global AOD, that is, more than a factor of 2 over the Arabian peninsula. Regional reductions in AOD occur mainly over the tropical oceans, associated with the enhanced aerosol wet removal by precipitation and clouds. Changes in sulfate and SOA (not shown) also contribute to the decrease of AODs over some polluted areas, including the Eastern United States and the adjacent North Atlantic, Southern Europe, and Southern Africa. A benefit of HR is the ability to capture highly resolved regional‐scale features, especially over the complex terrains of central Asia and south America. These resolution effects on AOD are particularly clear in Figure c, which shows AOD differences between HR and LRtunedHR simulations.

Annual mean aerosol optical depth (AOD) at 550 nm from the v1: (a) HR and (b) LRtunedHR. Panel (c) shows ΔAOD (HR‐LRtunedHR).

In addition, increased horizontal resolution also results in higher aerosol concentrations near the surface and lower concentrations aloft (not shown). This is consistent with increased sea spray and dust emissions and decreased turbulence removal of aerosols in the boundary layer. Weakened vertical transport by convection is also probably playing a role here; previous studies (Wang et al., ; Xie et al., ) have found that convective transport decreases at higher resolution because more grid‐scale lifting is resolved at a higher resolution. The resultant changes in aerosol vertical distributions are likely to have important implications on aerosol radiative forcing and cloud droplet activation—especially on the regional scales—and will be investigated in further studies.

3.4 Hydroclimate

Global‐average precipitation is slightly lower in HR than in LRv1 (Table ). As discussed in Golaz et al. (), deriving observational uncertainty for global‐average precipitation is challenging. Thus, while E3SMv1 global‐average precipitation is larger than the Rodell et al. () estimate used in Table and the GPCP estimate used in Figure , it is within uncertainty bounds for Stephens et al. () and Wild et al. (). LRtunedHR has lower precipitation, as expected because it is generally colder so has less moisture in the air. The fraction of precipitation or evaporation coming from ocean versus land closely matches observations, and hence moisture flux convergence over land is reasonable, but slightly weaker than observed (Table ).

Land and Ocean Water Budget

	HR	LRtunedHR	DECK	R15 Obs
Global‐ave Pr	2.96	2.82	3.07	2.79±0.15
Ocn Pr	78.6	79.4	77.4	77.7
Ocn Evap	85.6	85.2	84.8	86.5
Land Pr	21.4	20.6	22.6	22.3
Land Evap	14.4	13.9	15.2	13.6
Ocn→Lnd Transport	7.0	5.7	7.4	8.9

Note. Global‐ave Pr values are in mm/day, and all other terms are expressed as percents of global‐average precipitation. Observational estimates (“R15 Obs”) are taken from Rodell et al. (), which provides estimates based mainly on satellite measurements with help from data‐integrating models and combined in a variational framework that enforces several water and energy budget constraints.

Annual surface precipitation (in mm/day) from (a) Global Precipitation Project (GPCP; Adler et al., ) v2.2 observations at 2.5° × 2.5° resolution and model biases for (b) HR and (c) LRtunedHR.

Because precipitation depends strongly on topography and small‐scale baroclinic instabilities, improved precipitation is one of the main motivations for finer resolution. Thus, it is with keen interest that we examine Figure , which provides global maps of precipitation bias. One clear benefit of higher resolution is removal of dry bias to the south of Greenland. This dry bias is caused by erroneous freezing of the Labrador Sea at lower resolutions, which shuts off moisture flux from the ocean surface and therefore starves South Greenland clouds of moisture. See section for further discussion of improved Labrador Sea ice distribution at high resolution. Further corroboration of the importance of sea ice extent comes from Figure of Golaz et al.() that shows that Greenland precipitation bias in LRv1 disappears when SST and sea ice extent are prescribed. Wet biases over mountainous terrain (as described in Golaz et al., ) are slightly improved at large scales, but by no means alleviated. Unfortunately, increased resolution does not appear to correct E3SMv1's other major precipitation biases. The double ITCZ, dry Amazon, and wet Maritime Continent problems are virtually unchanged with changing resolution.

Even though increased resolution does not improve large‐scale features, it has clear benefit at regional scales. This is demonstrated in Figure , which zooms in on the west coast of the United States. Because snowpack is the main source of this region's drinking water, capturing the spatial distribution of precipitation over the U.S. West Coast is critical. While LRtunedHR is unable to capture the maxima of orographic precipitation over the region's mountain ranges, HR is better able to capture these features. Nevertheless, HR overpredicts the maximum value in the Sierras, a problem seen in other high‐resolution simulations of Sierra precipitation (e.g., Caldwell, ; Caldwell et al., ). The benefit of high resolution on stream flow is examined later in this section.

Climatalogical pattern of precipitation (in mm/day) over the west coast of North America as simulated by The North American Regional Reanalysis (NARR; Mesinger et al., ) and E3SM v1 model versions.

In addition to mean precipitation rate, the intensity of precipitation is important for climate impacts. To show the impact of grid resolution on precipitation intensity, the annual mean daily precipitation intensity distributions in the tropics (25°S to 25°N) are examined in Figure with a focus on heavy precipitation. For such histograms, it is more useful to focus on a particular region than on global averages in order to avoid the conflation of disparate changes in different regions. The tropics is a useful region because it dominates global precipitation, it involves significant contributions from both large‐scale and convective precipitation at all seasons, and it displays the largest sensitivity to resolution change. Histograms over the continental United States, however, showed similar but muted responses to resolution change (not shown).

Annual mean precipitation frequency of daily precipitation as a function of precipitation intensity between the models and observations over 25°S to 25° N; black: GPCP1DD; gray: CMORPH; green: TRMM; blue: LRtunedHR; and red: HR. Daily mean precipitation rates (unit: mm/day) gridded at 1° × 1° are used to derive the distribution. Any precipitation rates larger than 100 mm/day are assigned to the last bin that results in an uptick at the end of the plot. (a) Total precipitation, (b) total precipitation over land, (c) total precipitation over ocean, (d) convective precipitation, and (e) stratiform precipitation. Note that a different scale for Y axis is used for convective precipitation.

The observations are from the GPCP 1° daily (GPCP1DD) data (Huffman et al., ), and TRMM (Huffman & Bolvin, ) 3B42 and CPC MORPHing technique (CMORPH; Joyce et al., ) at 0.25° × 0.25°. The differences between the observational data sets provide a crude estimate of potential observational uncertainty. We use years from 1998 to 2013. Before aggregating the distribution, CMORPH/TRMM and modeled precipitation rates are remapped to the same 1° × 1° grids as the GPCP1DD data. For each day, each tropical grid cell is assigned to a precipitation bin using 1 mm/day bin widths. Similar to Wehner et al. (), we placed all precipitation rates larger than 100 mm/day in a single bin. This results in a change in slope at the right side of the plot.

Both LRtunedHR and HR models generally share the same “too frequent, too weak” problem as shown in other climate models (e.g., Stephens et al., ; Trenberth et al., ). They largely overpredict the frequency of precipitation occurrence for rainfall rates less than 15 mm/day and underestimate rainfall between 20 and 70 mm/day (Figure a). The high‐resolution model clearly shows an increase in the frequency of heavy precipitation rates (>50 mm/day) compared to its low‐resolution counterpart. For extreme precipitation rates in the tropics (>70 mm/day), the v1 HR simulation is within the range of the observations while the LRtunedHR simulation is considerably lower than the observed frequency.

The increase of heavy precipitation as model horizontal resolution increases is shown over both land and ocean (Figures b and c), though the increase over land is more significant, bringing the simulated distribution of intense precipitation closer to the observations. This increase is due to resolved‐scale (stratiform) precipitation (Figure e), which increases with increasing resolution for all bins greater than <0.2 mm/day. Convective precipitation, on the other hand, actually decreases with increasing resolution for all bins greater than 10 mm/day.

The diurnal cycle of precipitation is another long‐standing model bias, which one might hope is ameliorated by higher resolution. There is a slight tendency for precipitation to move later in the day over larger islands of the Maritime Continent, but in general, increased resolution does little to improve the model's tendency for precipitation to peak too early in the day (not shown).

In Figure , the land‐model‐simulated mean annual runoff averaged across the year 21–50 simulation period for both HR and LRtunedHR setups is compared with the half degree University of New Hampshire (UNH)/Global Runoff Data Centre (GRDC) Composite Monthly Runoff data (hereafter “GRDC runoff”; Fekete & Vorosmarty, ). The 10‐year‐long (1986‐1995) runoff data are essentially a data assimilation product that preserves the spatial specificity simulated by a large‐scale hydrologic model (Water Balance Model, WBM) while constrained by the in situ streamflow discharge measurements from the Global Runoff Data Centre (GRDC). Therefore, it is useful for spatial pattern comparisons at river basin scale level but may not be suitable as a benchmark product for bias identifications at the gridcell level. Our comparisons (Figures a, b, and c) suggest that both HR and LRtunedHR capture the general spatial pattern of the GRDC runoff. In arid regions such as the western United States, southern Africa, central Asia, and Australia, the model has a wet bias, while in relatively humid areas such as the Amazon basin and central Africa, dry bias is observed. These findings are generally consistent with the precipitation results (Figure ). The difference between HR and LRtunedHR (Figure d) indicates that the HR tends to better represent runoff heterogeneity in mountainous regions (e.g., Andes Mountains and Tibetan Plateau). In addition, the regions noted as having dry bias in LR are noticeably wetter in HR (e.g., central Africa) and regions that are too wet tend to be drier in HR (e.g., western United States).

Mean annual runoff (mm/year) from (a) GRDC, (b) HR, (c) LRtunedHR, and (d) the difference between HR and LRtunedHR.

To evaluate the performance of the E3SM river model (MOSART), we compared the MOSART‐simulated streamflow (years 21 to 50) with the Global Streamflow Indices and Metadata Archive (GSIM; Do et al., ) at 4,515 locations (Figure ). These gauges were selected to accommodate the need for global spatial coverage and for having river gauge drainage areas accurately represented in both HR and LRtunedHR models. Unsurprisingly, Figure shows that biases in runoff carry through into streamflow. Nonetheless, the overall performance of both simulations is very good, with R² over 0.8. In sites with large discharge (e.g., >10⁴ m³ s^‐1), HR's performance is generally better than LRtunedHR.

Mean annual streamflow at 4,515 gauge locations from (a) GSIM, (b) HR, (c) LRtunedHR, and (d) the scatter plot between GSIM observations and simulated streamflow.

In addition to the mean annual streamflow, we also examined the seasonality of the streamflow (Figure ) with color indicating the peak month of the hydrograph and the size of each dot indicating the seasonality index (SI) of the streamflow, which ranges from 1 to 12 and quantifies the level of seasonal variations of the hydrograph. SI equal to one indicates uniformly distributed hydrograph across the year, that is, no seasonal variation; SI equal to 12 indicates that peak streamflow occurs in a single month while the rest months are equal, that is, strongest seasonal variability. The equations for SI calculations can be found in Golaz et al. (). The results suggest that simulated streamflow is generally a good match with observations in terms of both seasonality and peak timing, though simulated SI is slightly overestimated for most rivers (Figure d). This is not surprising because the model did not represent dam operations, which usually lower the seasonality of the hydrographs. We do notice that SI is underrepresented in northern Eurasia in HR for unclear reasons (Figure b). There otherwise does not appear to be any systematic change in SI as a function of resolution.

Seasonality index (SI, size of the dots) (a) for GSIM gauge observations, (b) HR, and (c) LRtunedHR. Panel d provides a scatter plot of observed SI against simulated SI with size of the dots indicating magnitude of the discharge.

3.5 Ocean

A representation of the global upper ocean is shown in Figures as annual averages of Sea Surface Temperature (SST), Sea Surface Salinity (SSS), and Mixed Layer Depth (MLD) computed using the 0.03 kg/m³ density threshold criterion. Moving from low to high horizontal resolution has generally beneficial effects. In particular, the cold bias in the subtropical gyres is greatly reduced at HR, and the extremely cold and fresh biases in the Labrador Sea disappears. This latter bias is most likely related to sea ice, as we shall see in section . Also noticeable is the general improvement in SSS between LRtunedHR and HR, with LRtunedHR showing a fairly strong fresh bias in most of the global ocean. The origins of this fresh bias are being investigated in a separate study but that are believed to be at least partially due to ocean mixing. Biases in LRtunedHR SST and SSS look very similar to those shown for LRv1 in Golaz et al. (), implying that the aforementioned issues are intrinsic to the low‐resolution version of the model rather than the result of poor atmosphere tuning.

Annual‐mean SST from (a) the Hadley/OI observational product (Hurrell et al., ); (b) HR bias with respect to observations (model‐obs); and (c) LRtunedHR bias.

As for Figure but for the annual‐average SSS. Observations are from the Aquarius satellite product (Wentz et al., , years 2012‐2013).

As for Figure but for the annual‐average MLD. Observations are from the Argo‐based climatology of Holte et al. ().

One problem that increased resolution does not fix is the warm SST bias in the Eastern Boundary Current regions (e.g., west coast of North America, off the coast of Peru, and in the Buenguela region of South East Africa). Warm bias is caused by an unrealistic lack of stratocumulus clouds in this region (as shown in Figure ) combined with overly weak alongshore winds. Lack of improvement in eastern boundary current SST with increased resolution is surprising; all other high‐resolution modeling studies we know of found robust improvements in this feature with increasing resolution. Improvement is expected due to improved fidelity of alongshore winds, which is likely the result of a better representation of orography (Gent et al., ). Upwelling does strengthen in HR, stratocumulus do increase, and SST biases do improve with increasing resolution in E3SMv1 (not shown), but only along a very narrow strip just off the coast. Understanding why resolution sensitivity of eastern boundary upwelling has such small areal extent is ongoing research. We also note that while the SST bias patterns suggest a global average warm bias, the years 21–50 average has cooled relative to the initial condition (consistent with Figure a).

Mixed layer depth biases with respect to an Argo‐based climatology (Holte et al., ; Figure ) reflect a variety of upper ocean processes that are not explicitly resolved in models and whose representation depends on mixing parameterizations and shifts in buoyancy fronts (e.g., a shift in the Antarctic Circumpolar Current [ACC] relative to the Argo product). Similar to other Earth System Models, E3SM simulations tend to have larger MLD biases at higher latitudes where strong convection tends to occur. The large positive bias in HR in the Weddell Sea is due to an extensive Weddell Polynya that lasts a little over a decade (see also section ) and is associated with strong oceanic convection, and a warm and salty surface anomaly. The MLD in the Labrador Sea is too shallow in LRtunedHR, but too deep in HR. This is also visible in the regional MLD seasonal cycle in Figure . While the lack of convection in the Labrador Sea in LRtunedHR is most likely due to the unrealistically high sea ice coverage in this region, associated with the accumulation of freshwater and thereby strong stratification probably resulting from too weak Irminger and West Greenland Currents, the overly convective Labrador Sea in HR cannot be easily explained by surface temperature and salinity biases and needs further investigation (a separate study focusing on the Labrador Sea is underway).

March mean MLD (0.03 kg/m3 density threshold) in the Labrador Sea in HR (upper left panel) and LRtunedHR (upper right panel). White and black contours show the 15% and 85% sea ice fraction, respectively. Lower panel shows the seasonal cycle of MLD averaged over the gray box as shown in the upper panels for both simulations. Gray thick curve shows the same averaged MLD from Argo‐based climatology of Holte et al. ().

3.6 Sea Ice

The most important improvement to E3SM's surface polar climate in HR relative to LRtunedHR manifests as a correction to the winter sea ice edge in both the Northern and Southern Hemispheres. Figure summarizes this improvement, where we have used the 1979–1999 Meier et al. () National Oceanic and Atmospheric Administration Climate Data Record (NOAA CDR) as an ice edge benchmark for which there are regular measurements, contrary to the 1950s. At the time of writing the Walsh et al. () industrial‐era ice extent data set was being revised and is not included in this study. Comparing the ice edge of HR with LRtunedHR and observations from the last 20 years of the twentieth century, sea ice extent is more realistic at high resolution owing to reduced bias in the Labrador Sea, the Sea of Okhotsk, and around most of the Southern Ocean cryospheric perimeter (Figure ). Even though there is temporal mismatch between the perpetual 1950 simulations and 1979–1999 measurements, sea ice extent in the 1950 LRtunedHR run is similar to predictions from the late 20th century in the LRv1 simulation, suggesting that the temporal mismatch in observed versus modeled eras is not a cause of LRtunedHR bias (Golaz et al., ). There is also a strong similarity between LRv1 historical ensemble members' ice edge in the Southern Ocean and LRtunedHR.

Seasonal evolution of sea ice for high‐ and low‐resolution simulations for the Northern (a–c) and Southern (d–f) Hemispheres. Sea ice extent may be compared with the Meier et al. () NOAA Climate Data Record (CDR) for the period of 1979–1999 as an approximate guide to real‐world 20th century values. Points indicate monthly median for the designated 30‐year analysis period, and lines indicate the maximum and minimum monthly means. The numeric values on the x axes correspond of consecutive months of the year, where 1 represents January and 12 represents December.

Mean March (a–d) and September (e–h) sea ice thickness for HR and LRtunedHR simulations over the 30‐year climate averaging period. Data are rendered on native model grid cells with opacity (color density) determined by mean sea ice concentration, where fully saturated colors in the legend indicate 100% ice cover. The observed 1979–1999 Meier et al. () NOAA sea ice extent is traced in bold orange, equivalent to the 15% concentration contour. Extent contours appear in thin orange either side of the main measured ice edge for ±1 standard deviation in the observed ice fraction, where a large deviation from the bold orange contour indicates high interannual variability. Native model coastlines indicate the closure of Nares Strait (b and f) and where the top of the Antarctic Peninsula is pinched off (d and h) in the low‐resolution model.

In the case of the Labrador Sea low‐resolution sea ice bias, by far the largest marginal sea bias, the saline ice and its melt water separate early from the west Greenland current and spread toward Newfoundland along the streamlines in Figures e and f. As a consequence, the Labrador Sea becomes stably stratified, and oceanic mixing is shut down (Figure ). By contrast, this does not occur at high resolution, where sea ice in eddies spinning off the West Greenland Current soon melt away, animated in Roberts et al. (). The other significant Northern Hemisphere ice‐edge bias occurs in the Sea of Okhotsk, where low‐resolution bathymetry contribute to poor fidelity and where there is also high observed interannual ice‐edge variability, as indicated by the thin orange ±1 standard deviation 15% concentration measurements in Figure .

Seasonal mean Northern Hemisphere sea ice streamlines (black) and drift speed (shading) for years 21–50 of HR and LRtunedHR, truncated at 15% sea ice concentration. Median sea ice speeds from monthly model means within the respective sea ice zones appear in blue.

High‐resolution topobathymetry, greater extremes in wind and current forcing, and increased oceanic eddy kinetic energy combine to elevate the median and maximum sea ice drift speed in HR relative to LRtunedHR. Figures and demonstrate that sea ice is accelerated in HR relative to its low‐resolution analog. In particular, sea ice reaches much greater drift speeds in coastal or slope currents surrounding Greenland, next to Baffin Island, along the East Antarctica coast and in the Weddell Sea, verified independently of Figures and . These changes manifest from the sharper horizontal gradient ∇_h(f/H) in HR relative to LRtunedHR, for the Coriolis parameter f, and bathymetric depth H. There is also improved steering of the winter Southern Ocean ice edge, including from eddies downstream of Drake Passage, relative to LRtunedHR.

As for Figure but for the Southern Hemisphere.

As seen in the E3SMv0 high‐resolution coupled simulation (Kurtakoti et al., ), the HR simulation evolves intermittent Weddell Polynyas, imprinted in the model climate near Maud Rise (Figures j and k) but not present at low resolution (Figure l). This high‐resolution feature results from Taylor columns and mesoscale eddies around the Maud Rise seamount, where polynyas are initiated, a further ramification of ∇_h(f/H) fidelity; the associated anticlockwise surface circulation around Maud Rise is evident in HR but not LRtunedHR sea ice drift (Figures d and h). Evolution of one such Weddell Polynya may be seen in an animation of the last 10 years of HR sea ice thickness in Roberts et al. (). Improved coastal definition also means that East Antarctic polynyas are better resolved in HR relative to LRtunedHR; an example of the high‐resolution evolution of one such case, the Mertz Polynya, may be seen in Roberts et al. (). Increased ice production in polynyas at high resolution along the Antarctic coast is a function not just of coastal topology, but also of better resolved katabatics. In the Arctic, coastal definitions are important for refining ice and ocean fluxes through the Canadian Archipelago, including in Nares Strait, estimated to have a similar oceanic southward through‐flow as Lancaster Sound (McGeehan and Maslowski, ). It may also have a significant southward freshwater flux in the form of sea ice (Kwok et al., ). However, Nares Strait is completely closed‐off in low‐resolution v1 E3SM simulations, as seen in the native model coastlines in Figure (c) and (i).

Even though key polar surface circulation features including the Beaufort Gyre, Transpolar Drift, Ross Sea Gyre and Weddell Sea circulation are better resolved at high resolution (Figures and ), HR exhibits an ice cover that is too thin. We are working to quantify this bias precisely in ongoing work using ICESat, CryoSat‐2, and ICESat‐2 emulators that compare satellite‐derived and model sea ice freeboard. However, here, we surmise that thin ice is likely the cause of a negative ice extent bias everywhere in summer at high resolution (Figures and ). Mean ice thickness against the Canadian Archipelago never exceeds 3.5 m in HR except in slender coastal margins, even though submarine observations from 1958 to 1987 indicate vast year‐round ice buildup exceeding 5‐m draft, with mean draft of 7–8 m in some Canadian coastal margins (Bourke & McLaren, ). In the Weddell Sea, where Southern Ocean perennial sea ice is most prevalent, upward looking sonar deployed by Alfred Wegener Institute in the early 1990s revealed winter buildup exceeding 3‐m draft against the Antarctic Peninsula (Harms et al., ; Strass and Fahrbach, ; Roberts, ) as compared with HR mean ∼1,950 thickness of about 2 m or less. Ice advecting into the Weddell Gyre in 1992 along the West Antarctic coast near 15° W ranged in draft between about 0.9–1.7 m in March and 1.1–2.3 m in September (Roberts, ), whereas it remains thinner in E3SM in all mid‐20th century high‐resolution cases. By these measures, a much reduced thickness bias exists in LRtunedHR in both hemispheres, pointing to less bias in the polar surface heat budget, especially during summer, at low resolution.

4 Climate Variability

In addition to the mean state of the model, temporal variability is important. In this section we provide a short discussion of the main modes of climate variability.

4.1 SST/SSH Variability

As expected, the HR v1 ocean exhibits a robust mesoscale eddy field, with maximum SSH variability located primarily in western boundary current regions and the ACC (Figure ). The globally averaged SSH variability of the model (4.0 cm) compares well with the satellite altimetry value of 4.2 cm, and the correlation between the model and data is 0.73. Both the Gulf Stream and Kuroshio maxima are shifted about 50 km equatorward in the model when compared to AVISO, whereas the variability pattern in East Australia Current (EAC) matches very well. However, the magnitude of the variability in these three Western Boundary Currents is low by about 25% for the Northern Hemisphere currents, and almost 50% for the EAC. As is typical of this class of coupled models (Delworth et al., ; McClean et al., ; Small, Bacmeister, et al., ), the Gulf Stream Extension reaches too far across the North Atlantic and has a weak representation of eddies hooking around the Northwest Corner. The model Kuroshio exhibits two of the observed modes of multiyear variability: the Large Meander to the south of Japan around 138E (Qiu & Miao, ) and periods of alternating high/low eddy energy in the Kuroshio Extension (Douglass et al., ; Qiu & Chen, ).

Sea surface height variability (m) for (a) E3SMv1 HR (years 21–50 calculated from 5‐day snapshots) and (b) AVISO data (https://www.aviso.altimetry.fr) from satellite altimetry (years 1993–2013).

The low SSH variability in the Gulf of Mexico is due to lack of eddy shedding by the Loop Current for most of the simulation. Soon after the beginning of the run, the stratification in the core of the current drops, resulting in a baroclinically stable flow. It is not until year 38 that an intrusion of cooler water increases the stratification enough for the Loop Current to become unstable, with the first eddy that fully detaches and propagates westward occurring in year 42. The current then continues in a more realistic unstable state for the remainder of the simulation, detaching three more eddies during years 44–50.

Similarly to the Gulf Stream and Kuroshio, the SSH variability in the Agulhas Retroflection is shifted 50–100 km equatorward compared to AVISO. Somewhat remarkably, the meanders in the model match AVISO extremely well in both zonal position and wavelength. Agulhas eddies crossing the South Atlantic do not follow a single dominant path as appears to be the case in, for example, Delworth et al. (), but instead are more realistically distributed over a range of latitudes.

Figure shows the analogous SST variability (with the mean seasonal cycle removed) compared with AVHRR satellite observations. Unlike the SSH, which has no direct influence on the atmosphere, mesoscale SST features do have discernible effects on fluxes of heat and moisture to the atmosphere. SST variability has been shown to modify cloudiness, precipitation, and surface winds (Bryan et al., ; Byrne et al., ), while the location and strength of oceanic fronts can influence the path of storm tracks (Small, Tomas, et al., ). The bias in globally averaged variability is about 40%, or 0.1 °C, though locally the model variations can be as much as 2 °C too large. This suggests the possibility of unrealistically high mesoscale feedbacks in the model, though it is beyond the scope of this paper to quantify this. Despite the large difference in normalized globally averaged variance, the correlation between the model and data for SST (0.72) is essentially the same as for SSH.

Sea surface temperature variability (°C) with the mean seasonal cycle removed for (a) E3SMv1 HR (years 21–50 calculated from 5‐day snapshots) and (b) AVHRR (http://pathfinder.nodc.noaa.gov) satellite data (years 1982–2011 calculated from 5‐day snapshots).

4.2 El Niño‐Southern Oscillation

El Niño Southern Oscillation (ENSO) is a dominant mode of SST variability in the tropical Pacific (e.g., Figure a). Given the short simulation period, ENSO power spectra are not a reliable indicator of model performance (e.g., Wittenberg, ). Although the HR Nino 3.4 index looks reasonable in event magnitude and spacing, a single index only presents a small data point to the fidelity of ENSO in E3SM at HR. To better examine the influence of resolution on ENSO structure, we therefore examine the dominant patterns of SST variability in the HadISST data (Rayner et al., ), v1 HR, and v1 LRtunedHR. As expected, the ENSO SST pattern emerges as the first EOF of each SST data set (shown in Figure ). Broadly, each model captures the ENSO pattern, with subtle but important differences. At low resolution, the warm SST anomalies are too strong and extend too far west relative to observations, which is a common model bias (e.g., Menary et al., ). Surprisingly, increased resolution has little impact on ENSO pattern in E3SMv1. This is in contrast to studies with other models, which found an improvement in westward extent (e.g., Kirtman et al., ; Menary et al., ; Sakamoto et al., ). Across resolutions, the off‐equatorial ENSO response is muted, with the weakest response in HR. While the total tropical SST variance in HR is too large (Figure ), the fraction of SST variance associated with ENSO in HR is too small (as noted in Figure ) with the net result that the magnitude of ENSO‐related SST variations in HR is about right. In other words, tropical SST variance in the HR model run is too large due to sources other than ENSO. This is confirmed in spectra of Nino3, 3.4, and 4 SST anomalies (not shown).

First Empirical Orthogonal Function (EOF) of sea surface temperature for (a) HadISST data (1870–2010), (b) LRtunedHR (years 1–50), and (c) HR (years 1–50). The variance explained by each EOF is shown in the respective panels.

To explore possible mechanisms for the lack of improvement in ENSO in the HR simulation, we examine DJF thermocline depth and SST anomaly variability along the equator (Figures a and b), where we define the thermocline as the depth of the 20 °C isotherm (e.g., Kessler, ; Yang & Wang, ). The LRtunedHR simulation most faithfully represents the ARGO data, while the higher‐resolution simulation has a mean thermocline depth gradient, which is too weak. Previous studies (e.g., Kirtman et al., ) show that SST anomalies are better represented when model resolution is increased. This is also evident in our simulations, although unlike previous work, our LR simulation has more variance along the equator than observed and increasing resolution reduces bias by weakening variability. We also note that the SST anomaly variance observed is consistent with the simulated thermocline depths in Figure a. The HR model has the deepest thermocline depths and a more muted SST anomaly variance than the lower resolution simulation (but higher than observed). Despite the deeper mean thermocline depths, the HR simulation has a larger range of thermocline depth variability than LRtunedHR or observations.

DJF averages of (a) thermocline depth (defined as the depth of the 20 °C isotherm; Z20, where the shading represents the maximum and minimum depths observed during the period. (b) standard deviation of monthly SST anomalies and (c) average maximum vertical viscosity above the thermocline (solid lines) and standard deviation (shading). For all plots, fields have been averaged between 1° S and 1° N. The periods of analysis are years 1–50 for HR and for LRtunedHR. The ARGO data (Gaillard et al., ) in (a) run from 2013 to 2017, and the HadISST (c) is 1870 to 2015.

Deeper mean East Pacific thermocline depths and larger variability in the HR simulation may be responsible for the response seen in Figure (e.g., Jin, ). Large fluctuations in the thermocline depth imply large variations in oceanic heat content (OHC). Excess variation in OHC could reduce the persistence and magnitude of individual ENSO events, which could explain the weakened SST‐precipitation extremes correlation seen in Figure . Physically, we would expect the deeper depths in the East Pacific to result from either vertical advection or diffusion. The vertical velocity at the equator is roughly consistent between high and low resolution (not shown), given that the equatorial resolution in the LR configuration captures essential processes (e.g., Tropical Instability Waves). The vertical diffusivity above the thermocline in the tropical Pacific is shown in Figure c. In the Eastern Pacific, the maximum vertical diffusivity at and above the thermocline is higher throughout the tropical Pacific and increases toward the eastern boundary. It is possible that the increased diffusivity in the cold tongue results in the weaker east/west thermocline gradient, weakening the ENSO response at high resolution. However, it is still unclear what causes this increased diffusivity or if a reduction in diffusivity would improve the model representation of high resolution. This will be investigated in a future publication.

4.3 Atlantic Meridional Overturning Circulation

The Atlantic Meridional Overturning Circulation (AMOC) plays a dominant role in the Earth's climate (e.g., Cheng et al., ), where changes in AMOC can influence, among other things, North Atlantic storm tracks, Northern Hemisphere climate, and CO₂ sequestration. In LRv1 (Golaz et al., ), weak AMOC strength and variability was one of the most prominent model biases. Here we examine the influence of resolution on the modeled AMOC at the location of the RAPID array (Rayner et al., ). A few studies (Hirschi et al., ; Stepanov and Haines, ) have shown a mixed influence of resolution on the simulation of AMOC. These studies have found that there is a clear increase in AMOC variability across resolution, but the influence on the mean AMOC is less clear. Figure shows the time series of the maximum AMOC at 26.5°N, the location of the RAPID array. When resolution is increased the mean AMOC strength increases, with v1 HR being slightly too high.

Maximum AMOC strength (Sv) at 26.5°N, The orange line is HR and teal is LRtunedHR. The gray shaded region represents the mean plus and minus the range of variability observed at the RAPID station (variability is based on monthly observations).

Previous works (e.g., Cheng et al., ; Stepanov and Haines, ) suggest the AMOC strength could change in response to many variables, including surface forcing, ocean stratification, and northward heat transport. Given that the mean LRtunedHR AMOC is similar to LRv1, but with a very different atmosphere, it is unlikely that surface forcing is responsible for the change in magnitude. Figure shows that the near surface stratification is changed dramatically with higher resolution, most notably, the Labrador Sea convection is vigorous in HR and absent at LR. However, it seems unlikely that increased convection in the Labrador Sea is responsible for 10 Sv of increased transport from LRtunedHR to v1 HR. We have also examined mass transport through the Florida Straits in each simulation (not shown), and the HR simulation has much stronger transport through this region, suggesting a possible relationship between increased northward transport (likely due to better simulation of the Gulf Stream as suggested by Joyce and Zhang, in the HR simulation) and AMOC strength. The increased northward heat transport could lead to the improvement in the Labrador Sea ice thickness and extent (consistent with Rugenstein et al. ()). This possible relationship is further explored by examining the Atlantic Meridional Heat Transport (MHT) for each simulation (Figure ). All model results fall within the range of observational variability. Further, consistent with previous results (Volkov et al., ), MHT increases in magnitude due to the resolution of mesoscale eddies. It is possible that the increased heat transport encourages more evaporation in the North Atlantic, increasing deep convection, and hence the AMOC (e.g., Cheng et al., ; Stepanov and Haines, ).

Atlantic Meridional Heat Transport computed from NCEP and ECMWF reanalysis (Trenberth and Caron, ), HR (orange), and v1 LRtunedHR (teal). An estimate of observational error is shown in the shaded envelope. For HR and LRtunedHR, years 21–50 are used.

4.4 Madden Julian Oscillation

The Madden Julian Oscillation (MJO) is the dominant signal of intraseasonal variability in the the tropical atmosphere (e.g., Zhang, ). The MJO has broad reaching impacts through various teleconnections, including North American temperature and precipitation (e.g., Donald et al., ; Lin et al., ). Accurate simulation of the MJO could also lead to improved predictability (e.g., Jones et al., ). Therefore, it is important for E3SM to exhibit skill in simulating the MJO in order to have improved projections of the North American hydrological cycle. In Figure we see a robust simulation of the MJO at both resolutions. This is similar to the LRv1 simulations shown in Golaz et al. (). Increasing model resolution leads to no clear improvement in the modeled MJO, which is consistent with previous work (e.g., Kim et al., ).

Filtered wave number frequency power spectra for total precipitation (mm2/day2), shown as a ratio of the smoothed pseudo‐red background spectra of the same field following Wheeler and Kiladis (1999; 15° S to 15° N), for (a) TRMM, (b) LRtunedHR,and (c) HR with data interpolated to 1°.

Figure also shows a clear underestimate of the tropical wave variability in the HR simulations. In particular the higher‐requency Kelvin wave activity is near absent compared to observations. The tunings in the HR simulations do not themselves explain the lack of Kelvin Wave activity, since when these tunings are applied to a low‐esolution configuration (LRtunedHR; Figure b), the wave activity remains very similar to the default 1° configurations seen in Golaz et al. (), which do retain some Kelvin wave activity, albeit at spuriously high frequencies.

Recent theoretical work on the MJO (Chen & Wang, ; Wang et al., ) have proposed a dynamic moisture mode framework to fully encompass previous MJO theories. This work suggests that models that exhibit robust MJO power preferentially also have strong feedbacks between cloud radiative forcing, precipitation, and moisture convergence (Chen & Wang, ). Chen and Wang () also find that for models to correctly simulate the Kelvin wave component of the MJO, the feedbacks between boundary layer convergence and low level moisture and free tropospheric waves and convective heating must be correctly captured. It is possible that E3SMv1 has biases in one or both of these feedbacks, leading to the anomalously weak Kelvin wave signal seen in Figure . An exploration of these feedbacks is beyond the scope of this overview and is the subject of a future manuscript.

5 Extremes

Changes in the frequency and character of low‐robability events have great impact on people and the environments in which they reside. Scaling and theoretical considerations would favor more accurate representations of these events by higher‐esolution climate models. In this section, we test this assumption for a variety of extreme weather events.

5.1 Extreme Precipitation Over California

In this section, we evaluate the impact of horizontal resolution on extreme precipitation over California. California was chosen because of its large population, its importance to global food production, and its mountainous terrain, which one might hope is better captured at higher resolution.

Wintertime precipitation over California comes from moisture transported by the large‐scale circulation. This leads us to evaluate the water vapor transport and sea surface temperature anomalies, which give rise to wintertime extreme precipitation in this region. To avoid the need for high‐frequency output, this is done by first creating an extreme precipitation index by taking the sum of all daily precipitation values above the 95th percentile (computed separately for each model grid cell) within each month and calculating the DJF mean for each year, then computing the average index over California (defined here as 35–42° N, 125–118° W). The time series of DJF average of SST and vertically integrated water vapor transport (IVT) for every year are then regressed on this extreme precipitation index (separately for each model grid cell) and the resulting slope is plotted. This approach is presented in detail in Dong et al. (). Figure shows the resulting regressions of SST and the vertically integrated water vapor transport between 1,000 and 500 hPa on the winter extreme precipitation averaged over California. Panels a and b show values expected from observations. Both precipitation and water vapor transport in this figure comes from ECMWF Interim Reanalysis Data (ERA‐Interim) at a 1.5° × 1.5° resolution (Dee et al., ) and monthly SST is from the Hadley Centre Global Sea Ice and Sea Surface Temperature (HadISST, 1° × 1°, Rayner et al., ). Observations and reanalysis products from the 1981–2010 period are used. ERA‐Interim data show that extreme precipitation is associated with enhanced cyclonic circulation over the North Pacific, which transports moisture‐laden air from the tropics up the west coast of the United States. Both HR and LRtunedHR simulations are able to capture this basic pattern, though simulated moisture transport is a bit too strong in both models. In addition, moisture destined for California follows a slightly more zonal trajectory in HR due to the cyclonic trajectory over the North Pacific being shifted toward the coast and northward. This flow pattern is associated with an El Niño‐like SST pattern with warm anomalies in the tropical central‐to‐eastern Pacific. This tropical SST pattern is reasonably reproduced at low resolution but is entirely absent at high resolution.

Regression onto the time series of winter extreme precipitation of (left) monthly vertically integrated water vapor transport (kg·m−1 ·s−1) and (right) SST (° C) averaged over DJF from (a) ERA‐Interim, (b) HadISST, (c and d) LRtunedHR, and HR (e and f) runs. See text for details. Hatched areas indicate that the regression coefficients are statistically significant at the 95% level of confidence.

To further explore why the SST anomalies related to extreme precipitation in California are not reproduced in high‐resolution simulations, we examine the atmospheric response to ENSO during boreal winter. Figure shows the regression of precipitation and 500‐hPa geopotential height onto the Nino3.4 index. Observations (Figures a and b) show that El Niño conditions cause a weakened Walker circulation with reduced precipitation in the tropical western Pacific and increased precipitation in the central‐to‐eastern tropical Pacific. The tropical ENSO signal can be propagated to the North Pacific via the Pacific‐North American (PNA) pattern in geopotential height to influence precipitation in California (Figure b), consistent with Dong et al. (). The low‐resolution (LRtunedHR) simulation reproduces the PNA pattern as well as the mean precipitation response to ENSO (Figures c and d). The HR simulation, on the other hand, has a weaker precipitation response over the tropical eastern Pacific and a weaker PNA teleconnection over the North Pacific (Figures e–h), suggesting that the low resolution better reproduces the atmospheric response to ENSO events. These results are consistent with the weak relationship between winter extreme precipitation in California and SST in the high‐resolution simulations (Figure ) as well as the muted ENSO response seen in Figure .

Regression onto time series of the Nino3.4 index during winter (DJF) of (left) monthly mean precipitation (mm/day) and (right) 500‐hPa GPH (m; right) from (a) GPCP, (b) ERA‐Interim, (c and d) low‐resolution E3SM simulation, and (e and f) HR runs. The dotted areas indicate that the regression coefficients are statistically significant at the 95% level of confidence.

5.2 Tropical Cyclones

Tropical Cyclones (TCs) are some of the most destructive weather systems in the global tropics and subtropics with the potential to impact climate (Emanuel, ). It is thus very important to improve our understanding of the response of these features to changes in air/sea temperatures and shifting meteorological patterns. A rigorous analysis of TCs based on observations is often constrained by the length of the satellite record and by the lack of direct measurements of certain ocean‐atmosphere parameters (Henderson‐Sellers et al., ), motivating the use of high‐resolution numerical models to fill this knowledge gap. Consequently, a key driver of high‐resolution climate modeling efforts has been the realistic simulation of TCs and their climatology (Haarsma et al., ; Knutson et al., ; Walsh et al., ). With this consideration, we now examine TC simulation in E3SM.

In Figure , 50 years of TC tracks are shown from observations (Figure a) and E3SMv1 (Figures b and c). Observed TC tracks from the period 1966–2015 are obtained from the International Best Track Archive for Climate Stewardship (IBTrACS, Knapp et al., ). For the low‐resolution TC tracks, we use a 50‐year historical run of E3SM. In the model simulations, TCs are detected using TempestExtremes (Ullrich & Zarzycki, ; Zarzycki & Ullrich, ), a scale‐aware feature tracking software. TC‐like vortices are detected based on local minima in sea‐level pressure, positive temperature anomaly at the 200‐ and 500‐hPa levels to ensure a warm core, and the wind speed at the lowest model level. For further details regarding the TC‐detection algorithm, see Ullrich and Zarzycki (). Not surprisingly, in the low‐resolution version of the model, simulation of TCs is considerably underestimated (Figure b). On average, nearly 15 TC‐like features are generated globally in the low‐resolution version, which is a substantial underestimation when compared to observations (Figure a) where the annual global mean is about 92. However, this is in reasonable agreement with those produced by many other climate models at a similar resolution of about 100 km (Camargo, ). The annual mean global TC frequency is considerably improved in the high‐resolution simulations of E3SM (Figure c). On average, about 64 TC tracks are produced in each model year. When it comes to the distribution of TC tracks in different basins, we find spurious TC activity being produced in the subtropical South Atlantic and Southeast Pacific. The model also produces TC tracks in the Mediterranean region, which are likely related to Medicanes or TC‐like storms in that region (Cavicchia et al., ).

Global distribution of TC tracks from 50 years of (a) observations, (b) low‐resolution E3SM simulations, and (c) high‐resolution E3SM simulations. Tracks are color‐coded based on the along‐track intensity. The legend corresponding to the strength of the TCs is also shown.

Next, when it comes to intensity, the strongest TC produced in the low‐resolution E3SM is a Category 2 on the Saffir‐Simpson scale with a maximum wind speed of about 46 m/s. On the other hand, the most intense TC from the high‐resolution simulations nearly reaches Category 5 with a maximum lifetime wind speed of nearly 70 m/s (Figure a). The track of this TC is reminiscent of Supertyphoon Haiyan, which devastated the Philippines in 2013 (Lin et al., ). The cold SST wake produced by this TC is also shown in Figure a. Consistent with theory, the strongest cooling occurs on the right‐hand side of the track (Price, ). The magnitude of the cold wake increases with the strength of the storm and a maximum cooling of about 3 ° C is produced just before landfall. Similarly, Figure b shows an example of an intense TC in the Atlantic. The track of this TC is almost identical to that of Hurricane Andrew (1992), one of the most destructive hurricanes in Atlantic history (Pielke & Landsea, ). After forming in the Atlantic, the TC goes around the Caribbean islands and makes two landfalls, first as a Category 3 TC over South Florida and subsequently as a high Category 1 TC over the Northern Gulf of Mexico. Again, a cold SST wake can be seen, stronger on the right‐hand side of the track, with a maximum cooling of about 5 ° C in the Gulf of Mexico. While we have only presented an overview of TC simulation in this manuscript, we will follow this with another study focused on TCs with a more detailed description and analysis.

(a) Track of the most intense TC produced in the high‐resolution E3SM simulations, color‐coded by the along‐track intensity. (b) An example of an intense TC produced in the Atlantic basin. In both (a) and (b), the color in the background represents the difference between SST on the day of the TC's formation and SST a week later. The legend shown in (b) corresponds to the along‐track intensity of the TCs.

5.3 Atmospheric Blocking

Atmospheric blocking is a synoptic‐scale weather phenomenon associated with relatively stationary regions of anomalously high pressure. Given that this feature is associated with other forms of extreme weather, including heat waves, cold air outbreaks, and extreme precipitation, it is desirable to capture the frequency and character of these features with high fidelity. Global climate models, while able to approximate broad spatial blocking patterns, have traditionally struggled with correctly replicating observed blocking frequencies (Dunn‐Sigouin & Son, ). Because blocking typically emerges as the result of topographic drivers, improved model resolution (and subsequent improvements in the representation of topography) is expected to give rise to an improvement in model performance. However, Schiemann et al. () studied the sensitivity of blocking quality to resolution and only documented significant improvements in Euro‐Atlantic blocking in GCMs in the summer season, with practically no improvement in the Pacific basin.

The blocking detection algorithm used here is a modification of Dole and Gordon () that is detailed in Pinheiro et al. (). The algorithm searches pointwise for latitude‐adjusted 500‐mb geopotential height (Z500) anomalies that are greater than 1.5 times the standard deviation of the Z500 anomaly at that point, calculated for each day of the year. That is, we define [Image Omitted. See PDF] where Z^* is the Z500 anomaly, Z is the current Z500 value, $\bar{Z}$ is the long term daily mean (LTDM) of the Z500 field at each day of the year, $Z_{adj}^{*}$ is the latitude‐adjusted Z500 anomaly, and ϕ is the latitude. An additional minimum anomaly threshold of 100 m is imposed at coordinates where the calculated threshold value is low (due to anomaly time series with relatively low anomaly magnitudes and small standard deviation value).

Blocks from each data set were identified separately for the Northern and Southern Hemispheres, using the StitchBlobs program in TempestExtremes (Pinheiro et al., ; Ullrich & Zarzycki, ) with the normalized anomalies as the input. The detection algorithm was limited to a latitude range of 25–75° in each hemisphere, with blob persistence of 5 days, minimum size of 10⁶ km², and a minimum 50% overlap of blobs between time steps. The results of this analysis are depicted in Figure . Performance in the Southern Hemisphere was largely unchanged regardless of model resolution, so we focus on analyzing the Northern Hemisphere. Analogous to the results of Schiemann et al. (), the NH summertime blocking density exhibits much smaller biases in HR compared to LRtunedHR across all regions. However, there is little improvement in the NH wintertime blocking density in HR compared with LRtunedHR. Further, a pattern of underestimated blocking density at lower latitudes and overestimated density at higher latitudes emerges in the high‐resolution run. This pattern suggests that blocks are forming too far north, likely because of a bias toward weaker climatological winds that are more susceptible to curving at higher latitudes, and stronger zonally directed winds at lower latitudes. This pattern is particularly prominent over the Northern Pacific. Nonetheless, there is some improvement in the overestimation of wintertime blocking in the Euro‐Atlantic region.

Pointwise frequency of 500‐hPa blocking for (left) December‐January‐February (DJF) and (right) June‐July‐August (JJA). The rows indicate (a) ERA‐Interim reference, and differences between (b) E3SM HR and reference, and (c) E3SM LR and reference.

6 Model Sensitivity

In Golaz et al. (), LRv1 was shown to reproduce the observed total warming over the 20th century. This warming was, however, accomplished through a balance between extremely high equilibrium climate sensitivity (5.3 K) and very strong aerosol response as encapsulated in LRv1's effective radiative forcing (ERF) of −1.65 W/m². These sensitivities were argued to be too large because their combined effect results in underpredicted warming during years when aerosol emissions were particularly strong and overly strong warming in years when aerosol load was reduced. In this study, we will look at the impact of changing resolution and tuning on these sensitivities. Aerosol sensitivity is expected to decrease with resolution (Ma, Rasch, et al., ). Climate sensitivity is also likely to change because more scales are explicitly resolved at higher resolution.

Since we cannot afford to run the fully coupled HR model for the hundreds of years necessary to compute transient or equilibrium climate sensitivity, we instead diagnose the net climate feedback parameter (λ) from the difference between a 5‐year atmosphere‐land simulation with SST and sea ice extent prescribed at current‐climate conditions and a simulation that is identical but has SST increased uniformly by 4 K. This approach, popularized by Cess et al. (), allows one to compute the effective net feedback from the change in TOA radiative imbalance ΔF_TOA caused by forcing global‐average surface temperature T_glob ave to change: [Image Omitted. See PDF]

Ringer et al. () show that λ from Cess experiments provides a reasonable, though slightly high, approximation to values obtained by more expensive methods. Terms on the right‐hand side of equation  are taken from 5‐year simulations with annually repeating forcings representative of the decade centered on 2010. This choice differs from typical Cess experiments (which use preindustrial forcings) and was made because preindustrial fixed‐SST model forcings were not yet available. Simulation year is unlikely to have a leading‐order effect on the sensitivity of λ to resolution. SST and sea ice are actually prescribed from year 2013 because it was more representative of decadal‐average behavior than 2010. For more details on configuration, see section .

The first row of Table provides λ values for HR, LRtunedHR, and LRv1 simulations. The direct effect of increasing resolution causes a slight strengthening of λ, indicating stronger resistance to temperature change and therefore weaker climate sensitivity. Nevertheless, changing parameterization constants has a larger impact than resolution, as can be seen by the differences between the two low‐resolution simulations.

Net Climate Feedback Parameter (λ, unit: W/m²/K) and Aerosol‐Related Effective Radiative Forcing (ERF, unit: W/m²) of Different E3SM Versions

	HR	LRtunedHR	LRv1
λ	−1.29	−1.20	−1.36
ERF	− 1.37	−1.37	−1.64
ERFaci	−1.49	−1.51	−1.76
ERFari	0.04	0.05	0.03
Surf. albedo ERF	0.09	0.08	0.10

Aerosol sensitivity was also computed from 15‐month simulations with horizontal winds nudged toward ERA‐Interim using a nudging timescale of 3 hr. Kooperman et al. () demonstrate that nudged simulations of this type allow for accurate aerosol sensitivity predictions from 1‐year runs by removing weather noise differences between the simulations. Direct aerosol effects are diagnosed by calling radiation twice in each run: once with modeled aerosols and again with clean conditions. Indirect aerosol effects are computed by repeating our present‐day nudged simulation with anthropogenic aerosols and precursor gases replaced with preindustrial values. As explained and defended in Ghan (), this approach allows the total aerosol ERF to be decomposed into ERF due to aerosol‐cloud interaction (ERF_aci), ERF due to aerosol‐radiation interaction (ERF_ari), and ERF due to surface albedo change. Note that changes in land use/land cover may also induce resolution dependence in aerosol changes because higher‐resolution results in stronger wind maxima that lofts more dust (section ; this effect is not included in the present analysis.

ERF, ERF_ari, and ERF_aci are provided for all three model configurations in Table . ERF_ari is positive (warms the planet) due to the strong absorption of anthropogenic black carbon (BC) in the cloudy sky, especially when BC resides above clouds and bright snow/ice. The reduction in surface albedo straightforwardly leads to a positive ERF change. ERF_aci is strongly negative in all simulations due to cloud brightening with increased cloud condensation nuclei. Surprisingly, increased resolution has no effect on ERF. ERF_aci does weaken as expected, but its resolution sensitivity is very small and is counterbalanced by a corresponding weakening of ERF_ari. Tuning has a much larger effect on ERF than resolution change. ERF_aci in particular is weakened by ∼15% due to tuning changes between LRtunedHR and LRv1. This emphasizes the need to understand aerosol sensitivity to parametric choices. More detailed analyses of E3SMv1 aerosol distribution and sensitivity will be presented in a future paper.

7 Conclusions

In this study, we describe the high‐resolution (HR) configuration of E3SMv1 and document its simulation characteristics in a 50‐year‐long climatological control run with time‐invariant 1950 forcings following the HighResMIP protocol. The atmosphere and land components have a horizontal resolution of approximately 0.25° longitude and latitude, which is adequate for realistically simulating most weather phenomena. The variable resolution ocean and sea ice grid employs an 18‐km horizontal grid spacing near the equator, which smoothly reduces to approximately 6 km near the poles. This is sufficient to capture the most energetic motions in the ocean, which are poorly represented in standard resolution coupled climate models.

The HR version of E3SMv1 is superior to its low‐resolution counterpart (and most CMIP5 models) in global RMSE for most important variables. Improvement results from a combination of resolution increases, related tuning changes, and perhaps from using a short initialization period. Improved resolution was particularly beneficial to ocean and ice components. Most strikingly, the excessive sea ice extent in the low‐resolution (LR) simulation is eliminated in our HR run, with a resulting reduction of cold and fresh biases in that region. Increased resolution also strengthened (improved) the Atlantic Meridional Overturning Circulation and associated meridional heat transport, due to the improved Gulf Stream. Reduced ice coverage in the Labrador Sea also has important beneficial impacts on precipitation bias south of Greenland. In this region and elsewhere, improvements come from better‐resolved bathymetry, greater extremes in winds and currents, and strengthened ocean eddy kinetic energy. Unlike its low‐resolution counterpart, the HR model is also able to capture large polynas typical of those observed in the Weddell Sea. Unfortunately, improvements in ice extent do not translate to improvements in ice thickness, which is generally too thin in our HR simulations. Another major improvement in HR relative to LR is in the spatial distribution of SST and salinity.

As expected, regional cloudiness and precipitation patterns are improved by introducing finer‐scale topographic detail. Likewise, hurricanes are better captured, runoff and streamflow are generally improved, extreme precipitation events are better simulated, and the ratio of convective to stratiform precipitation decreases at higher resolution. For aerosol, increased extremes in winds at higher‐resolution result in higher dust and sea‐salt emissions, while enhanced spatial inhomogeneity in clouds and precipitation at finer resolution increases aerosol wet removal. Shift from parameterized convective transport to resolved‐scale vertical advection reduces dry deposition at the surface as more aerosols are lofted into the upper troposphere. The net effect is an increase in global‐average natural aerosol and a decrease in anthropogenic aerosol. The regional picture is more complicated. Clouds are relatively unaffected by resolution, though there is a slight tendency for cloudiness to decrease at all levels as resolution increases.

Interestingly, several features found to improve with resolution in other studies did not do so in our simulation. Most notably, stratocumulus and associated warm SST biases in the eastern subtropical oceans were not substantially improved at higher resolution. This lack of improvement occurs because unlike other studies, increased wind stress‐induced oceanic upwelling with finer resolution is confined in E3SMv1 to an area that is negligibly small compared to the size of the region of deficient low clouds. Another surprise was degradation in the westward extent of ENSO with increasing resolution. Overly deeper East Pacific thermocline depths in the HR simulation offer one possible explanation for the lack of ENSO improvement, but more study is needed. Oddly, MJO‐related Kelvin wave energy also disappears in E3SMv1 as resolution increases. Another surprise was that the diurnal cycle of precipitation was almost completely unaffected by resolution change, even in areas where sea breezes are expected. Additionally, while better topography did improve the frequency of Northern Hemisphere summertime blocking events, wintertime frequency actually got worse.

Another interesting finding is that large‐scale patterns of climatological‐average cloudiness and precipitation were relatively insensitive to resolution change. Climate and anthropogenic‐aerosol sensitivity are also almost completely unaffected by changes in resolution when tuning parameters were held fixed. Lack of aerosol sensitivity change occurred even though ERF_aci is expected to weaken with resolution (Ma, Chuang, et al., ) because that weakening was slight and was completely compensated by increased ERF_ari. The tuning needed to maintain energy balance and model skill with changing resolution, in contrast, had a relatively large effect on climate and aerosol sensitivity. Improved prediction at global scales due to resolution of important fine‐scale processes is a major motivator for the push to higher resolution. In this context, the finding that model sensitivity in E3SMv1 is relatively insensitive to resolution is important and perhaps reassuring. There are, however, several caveats to this result worth keeping in mind. First, this finding is from a single model. Other models may respond differently. In particular, models whose stratocumulus increases with resolution are likely to have stronger changes in climate and aerosol sensitivity with resolution because stratocumulus changes have such a powerful impact on TOA radiation. Second, this finding is specific to the transition from 100‐ to 25‐km atmospheric grid spacing. It could be that even higher resolution is needed in order to capture the processes governing model sensitivity.

In summary, the high‐resolution coupled version of E3SMv1 is a world‐class climate model. Many aspects of model climate benefit from finer resolution. Ocean and sea ice behavior and fine‐scale atmospheric features are particularly improved by increased grid spacing. Surprisingly, several features that improved with resolution in previous high‐resolution studies did not improve in this study. In addition, anthropogenic‐aerosol and climate sensitivity were almost completely unaffected by resolution change. These results suggest that while beneficial, increased resolution does not solve all model problems.

Acknowledgments

This research was supported as part of the Energy Exascale Earth System Model (E3SM) project, funded by the U.S. Department of Energy, Office of Science, Office of Biological and Environmental Research (BER). L. D. and S. A. K. acknowledge support from the Regional and Global Modeling and Analysis Program of BER. The v1 high‐resolution simulation was run using resources of the Argonne Leadership Computing Facility at Argonne National Laboratory, which is supported by the Office of Science of the U.S. Department of Energy under Contract DE‐AC02‐06CH11357. The LRtunedHR simulation was performed on Blues, a high‐performance computing cluster operated by the Laboratory Computing Resource Center at Argonne National Laboratory. This research also used computing resources at the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. DOE under Contract DE‐AC02‐05CH11231. The AVISO (https://www.aviso.altimetry.fr) Ssalto/Duacs altimeter products were produced and distributed by the Copernicus Marine and Environment Monitoring Service (CMEMS). The AVHRR Pathfinder data was obtained from the U.S. National Oceanographic Data Center and GHRSST (http://pathfinder.nodc.noaa.gov). Model output presented in this paper is available online (from https://e3sm.org/data/get-e3sm-data/released-e3sm-data/v1-025-data-hiresmip/). Sandia National Laboratories is a multimission laboratory managed and operated by National Technology and Engineering Solutions of Sandia, LLC, a wholly owned subsidiary of Honeywell International Inc., for the U.S. Department of Energy's National Nuclear Security Administration under Contract DE‐NA0003525. This paper was prepared by LLNL under Contract DE‐AC52‐07NA27344. This paper describes objective technical results and analysis. Any subjective views or opinions that might be expressed in the paper do not necessarily represent the views of the U.S. Department of Energy or the U.S. Government.