This paper is organized as follows. In section 2 , we present an overview of the model components, with a specific focus on new model components (i.e., those not previously used in coupled Earth system modeling), as well as new developments in the atmosphere and land models. In section 3 , we describe the E3SMv1 initialization and spin‐up procedure, including model tuning objectives and simulation campaign. Section 4 provides an overview of the preindustrial control simulation and section 5 an analysis of the E3SMv1 climate in the historical simulations including short‐ and long‐term variability. An extended discussion of E3SMv1's effective radiative forcing (ERF) and climate sensitivity is provided in section 6 . Finally, we offer a summary of E3SMv1 fidelity and discuss future directions for the fully coupled model in section 7 .

Given that E3SMv1 features three specific models of components (ocean, sea ice, and river) that have never been used in a coupled Earth System Model and that there were significant developments in the atmosphere and land components, an examination of the model behavior relative to observations and other CMIP‐class models is needed. We analyze CMIP6 Diagnosis, Evaluation, and Characterization of Klima (DECK) and historical simulations (Eyring et al., 2016 ) performed with E3SMv1. This allows for a rigorous comparison of E3SMv1 behavior against observations and many other models. This work will also provide a baseline for all future E3SM developments and experiments.

The focus here is on the physical model at standard resolution useful for simulations like those specified in Coupled Model Intercomparison Project Phase 6 (CMIP6). This includes a 1° atmosphere and land (equivalent to 110 km at the equator), 0.5° river model (55 km), and an ocean and sea ice with mesh spacing varying between 60 km in the midlatitudes and 30 km at the equator and poles. A higher‐resolution configuration with a 0.25° atmosphere and land, 0.125° river model, and ocean and sea ice with mesh spacing between 18 km at the equator and 6 km at the poles (roughly equivalent to a 0.1° resolution) will be documented in a subsequent paper.

E3SM version 1 (E3SMv1) was branched from the Community Earth System Model (CESM1; Hurrell et al., 2013 ) but has evolved significantly since. E3SMv1 consists of three coupled modeling systems with varying degrees of sophistication. The present work describes the physical Earth system model that represents water and energy cycles in atmosphere, ocean, sea ice, land, and river components. This configuration is aimed at addressing the DOE water cycle science questions relating to interactions between the water cycle and the rest of the human‐Earth system on local to global scales, water availability, and water cycle extremes. This physical model also serves as a foundation for two additional configurations: (i) a biogeochemistry configuration with interactive nitrogen and phosphorous for interactions between biogeochemical cycles and other Earth system components and (ii) a cryosphere configuration with added interactive ice shelf cavities for assessing the impacts of ocean‐ice shelf interactions on Antarctic Ice Sheet dynamics and the implications for sea level rise.

Scientific developments in E3SM are dictated by three science drivers that broadly cover the foundational science for advancing Earth system prediction. Notably, water cycle, biogeochemistry, and cryosphere systems govern variability and changes in water availability and storms, air and river stream temperature, and coastal flooding and sea level rise that are all critical to the energy sector.

In 2013, the U.S. Department of Energy (DOE) developed a report summarizing observed long‐term trends that, if continued for several decades, would have major impacts on the energy sector (U.S. Department of Energy, 2013 ). Among these were regional trends in air and water temperatures, water availability, storms and heavy precipitation, coastal flooding, and sea level rise. The ability to simulate and predict significant, long‐term changes in these environmental variables important to energy‐sector decisions required capabilities beyond the existing state‐of‐the‐science Earth system models. The Energy Exascale Earth System Model (E3SM) project was conceived to meet this mission need (Bader et al., 2014 ).

Typical processor layout, performance, and sequencing for all components in E3SMv1. Components include coupler (CPL), atmosphere (ATM), ocean (OCN), sea ice (ICE), land (LND), and river runoff (ROF). The hatched area represents unoccupied time. OCN runs concurrently on its own nodes and is load‐balanced with ATM + ICE, which run sequentially. Driver overhead is 4%. Other layouts are possible, but we have found this to be the most robust and efficient.

Simulations in this work were performed on DOE's National Energy Research Scientific Computing Center (NERSC) Edison supercomputer ( http://www.nersc.gov/users/computational-systems/edison ). Using 285 nodes (each consisting of two 12‐core Intel “Ivy Bridge" processors at 2.4 GHz) with 24 MPI (Message Passing Interface) tasks per node, the coupled model performance averages 10 simulated years per day. Figure 1 illustrates the sequencing, processor layout and relative cost for all components in E3SMv1.

cpl7 allows several time sequencing options for the components and flexibility for how often the components communicate with the coupler. E3SMv1 uses “RASM_OPTION1" where the sea ice, land, and river runoff models execute simultaneously and in sequence with the atmosphere. The ocean model runs simultaneously with all four of those components. Model time steps are as follows:

cpl7 implements an online‐off‐line method for interpolating values between different grids in E3SMv1. Grid intersections and interpolation weights are calculated with an off‐line tool and then read in at runtime and applied during the coupled integration. cpl7 also allows flux and state variables to be interpolated with different weights. The number of different weight files needed is relatively small because the atmosphere and land models share a grid, as do the ocean and sea ice models. The interpolation weights for nearly all variables are calculated with the ESMF_RegridWeightGen program from ESMF (Collins et al., 2005 ) using the first‐order conservative option. However, interpolation weights for state variables from the atmosphere to the ocean/sea ice grids are calculated using TempestRemap (Ullrich & Taylor, 2015 ). Unlike ESMF, the TempestRemap algorithms have native support for vertex centered finite element grids, such as used by the E3SM atmosphere dycore. Future versions of E3SM will use TempestRemap for all mapping weight calculations.

E3SMv1 uses component coupling software from the Common Infrastructure for Modeling the Earth. The top‐level driver is cpl7 (Craig et al., 2012 ), which provides a main program for forming the single executable of E3SM, directs the time integration of the coupled model, and performs any necessary interpolation or time averaging needed between the components. Communication between the parallel components of E3SM and parallel interpolation is provided by the Model Coupling Toolkit (Jacob et al., 2005 ; Larson et al., 2005 ), which is included in Common Infrastructure for Modeling the Earth.

E3SMv1 replaced the River Transport Model (Branstetter & Erickson, 2003 ), a linear reservoir routing model used in CESM1 and E3SMv0, with the Model for Scale Adaptive River Transport (MOSART), which uses a physically based kinematic wave approach for river routing across local, regional, and global scales (Li et al., 2013 , 2015 ). In the standard E3SMv1 resolution, MOSART uses a regular latitude‐longitude grid with spacing of 0.5°. Surface and subsurface runoff simulated by ELM are mapped from the ELM grid to the 0.5° latitude‐longitude grid as input to MOSART, which routes the runoff and provides freshwater input to the ocean model. MOSART does not exchange water with the atmosphere or return water to the land model. MOSART divides each grid cell into three categories of hydrologic units: hillslopes that contribute both surface and subsurface runoff into tributaries, tributaries that discharge into a single main channel, and the main channel that connects the local grid cell with the upstream/downstream grid cells through the river network. Two simplified forms of the one‐dimensional Saint‐Venant equations are used to represent water flow over hillslopes, in the tributary, or in the main channels. MOSART only routes positive runoff, although spurious negative runoff can be generated occasionally by the land model. Negative runoff is mapped directly from the grid cell where it is generated at any time step to the basin outlet of the corresponding MOSART grid cell. More detailed descriptions of MOSART and its input hydrography data and channel geometry parameters can be found in Li et al. ( 2013 ). When driven by runoff simulated by CLM4 with observed meteorological forcing data, streamflow simulated by MOSART is shown to reproduce the observed annual, seasonal, and daily flow statistics at over 1,600 stations of the world's major rivers reasonably well in terms of the overall bias and the seasonal variation (Li et al., 2015 ).

Several new developments have been made within ELMv0 since its branch point from CLM4.5. First, an extended representation of the influence of aerosols and black carbon (BC) deposition on snow was introduced, based on data summarized by Liu et al. ( 2012 ), causing aerosol optical properties within the snowpack to vary as a function of snow grain size and aerosol/ice mixing state (Flanner et al., 2012 ). This development also fixed a bug related to calculation of snow grain size following snow layer division. Second, a minor modification was made to reduce the rate of evaporation from area characterized as pervious road under dry conditions. Third, the numerical scheme for calculation of leaf stomatal conductance was updated to prevent nonphysical simulation of negative internal leaf CO 2 concentrations. Fourth, albedo calculation was modified to return a land albedo of 1.0 in land cells and time steps when the sun is below the horizon.

The E3SM Land Model, version 0 (ELMv0), is the land component in E3SMv1. ELMv0 adopts many of its capabilities from its parent model, the Community Land Model version 4.5 (CLM4.5). Oleson et al. ( 2013 ) provide a detailed technical description of that parent model, including descriptions of all fundamental equations representing the model's biogeophysical and biogeochemical dynamics. The model describes interactions of the land surface with the near‐surface atmosphere and includes interactions among land subsystems such as vegetation, soil, snow, groundwater, runoff, urban areas, and managed ecosystems. These interactions are represented through conservation equations of state for energy, water, carbon, and nitrogen, including terms that couple these states and fluxes. While ELMv0 includes prognostic carbon and nitrogen biogeochemistry and dynamic ecosystem structure, the E3SMv1 water cycle experiments do not make use of these capabilities. Instead, these simulations use a static representation of land ecosystem structure with prescribed seasonal changes in vegetation canopies based on a climatology of satellite remote sensing data. This simulation mode is referred to as satellite phenology (SP) and is based on methods described by Lawrence and Chase ( 2007 ). In this simulation mode, regardless of long‐term climate changes or short‐term anomalies, the global distribution and seasonal variation of vegetation structure are constant over the course of the simulations. This includes model representation of vegetation height, timing and amount of displayed leaf area, and vegetation contributions to land surface albedo.

Coupling of the sea ice component to the ocean takes advantage of z star ocean coordinates as described by Campin et al. ( 2008 ) and is a departure from the coupling of CICE and POP (Parallel Ocean Program) in CESM1. The weight of sea ice contributes to the ocean's barotropic mode, notably affecting the free surface over continental shelves. In shallow water depths at or less than the floating ice draft, the weight passed to the ocean model is limited to prevent evacuation of the underlying liquid column. When frazil ice forms in the ocean model, the volume of newly formed crystals is passed to the sea ice model with a fixed salinity of 4 PSU, rather than exchanging a freezing potential as in other models. Future versions of E3SM will permit progressive brine drainage to the ocean from the mushy layer physics used in MPAS‐Seaice (Turner & Hunke, 2015 ). For E3SMv1, brine drainage occurs internally in MPAS‐Seaice for thermodynamic calculations, but for the sake of freshwater coupling, the ocean model only receives mass fluxes back from melted sea ice at the fixed salinity that it originally passed to its cryospheric counterpart (4 PSU). The ocean temperature immediately under the ice is the same as the liquid phase in the lowest layer of the sea ice model and is not fixed at −1.8 °C as is typical of previous generation coupled models (Naughten et al., 2017 ). For the current version, we have addressed these long‐standing ocean‐ice coupling issues identified by the modeling community: explicit sea ice mass and salt exchange, a pressure force of the ice on the ocean, a basal sea ice temperature consistent with the ocean model's equation of state, and resolved inertial oscillations (Hibler et al., 2006 ; Lique et al., 2016 ; Schmidt et al., 2004 ).

MPAS‐Seaice is the sea ice component of E3SMv1. MPAS‐Seaice and MPAS‐Ocean share identical meshes, but MPAS‐Seaice uses a B‐grid (Arakawa & Lamb, 1977 ) with sea ice concentration, volume, and tracers defined at cell centers and velocity defined at cell vertices. Velocity components at cell vertices are not aligned with the mesh, as in sea ice models with structured meshes and quadrilateral cells. Instead, the velocity components are aligned with a spherical coordinate system that is locally Cartesian and has its poles on the geographical equator. Velocities are determined by solving the sea ice momentum equation (Hibler, 1979 ; Hunke & Dukowicz, 1997 ) on cell vertices in an identical manner to the CICE model (Hunke et al., 2015 ) except for the internal stress term. MPAS‐Seaice uses the variational formulation of the internal stress term used by CICE (Hunke & Dukowicz, 2002 ) but modified to use the nonregular polygons of MPAS meshes. Instead of the bilinear basis functions used by CICE, MPAS‐Seaice uses Wachspress basis functions (Dasgupta, 2003 ), which are integrated with the quadrature rules of Dunavant ( 1985 ). Horizontal transport of ice concentration, volume, and tracers is achieved with an incremental remapping scheme similar to that described in Dukowicz and Baumgardner ( 2000 ), Lipscomb and Hunke ( 2004 ), and Lipscomb and Ringler ( 2005 ) but adapted to MPAS meshes. MPAS‐Seaice shares the same column physics code as CICE through the Icepack library (Hunke et al., 2018 ). For simulations shown here, MPAS‐Seaice uses the “mushy layer” vertical thermodynamics scheme of Turner et al. ( 2013 ) and Turner and Hunke ( 2015 ), the level ice melt pond scheme of Hunke et al. ( 2013 ), a delta‐Eddington shortwave (SW) radiation scheme (Briegleb & Light, 2007 ; Holland et al., 2012 ), a scheme for transport in thickness space (Lipscomb, 2001 ), and a representation of mechanical redistribution (Lipscomb et al., 2007 ).

E3SMv1 standard resolution simulations employ the classic Gent and McWilliams ( 1990 ) eddy transport parameterization. The Gent‐McWilliams bolus coefficient was tuned, in part, to reduce the transport of heat to depth in the Southern Ocean, to a value of 1,800 m 2 /s for the simulations presented here. The Redi coefficient, which adds diffusion along isopycnal layers, was set to 0 for this set of simulations. In the horizontal, biharmonic viscosity is used for momentum (1.2 × 10 9 m 4 /s 2 ). No explicit horizontal tracer diffusivity is included.

The simulations presented here use a z star vertical coordinate within an arbitrary Lagrangian‐Eulerian scheme, where the layer thicknesses of the full column expand and contract with the sea surface height (Petersen et al., 2015 ; Reckinger et al., 2015 ). The prognostic volume‐based equation of motion includes surface mass fluxes from the coupler; thus, virtual salt fluxes are not needed. Vertical mixing is computed implicitly at the end of each time step using the Community Vertical Mixing project implementation of the K‐profile parameterization as described by Van Roekel et al. ( 2018 ) where our configuration of K‐profile parameterization is based on the results of comparison against large eddy simulations.

MPAS‐Ocean, based on the MPAS framework (Ringler et al., 2010 ), is the ocean component of E3SMv1 (Petersen et al. 2019 ). MPAS‐Ocean uses a mimetic finite volume discretization of the primitive equations and invokes the hydrostatic, incompressible, and Boussinesq approximations on a staggered C‐grid (Arakawa & Lamb, 1977 ; Ringler et al., 2013 ; Thuburn et al., 2009 ). MPAS‐Ocean grid cells for E3SMv1 simulations are near‐hexagons (five or more sides), but the MPAS framework supports cells with any number of sides; the algorithms and code are identical for all cell shapes. The tracer advection scheme is the quasi third‐order flux corrected transport (FCT) scheme (Skamarock & Gassmann, 2011 ) with separate limiting in the horizontal and vertical. The MPAS‐Ocean time stepping method is split‐explicit, where the barotropic component is subcycled within each baroclinic time step.

The computational cost of EAMv1 increased by approximately a factor of 4 relative to CAM5 due to higher vertical resolution and parameterization complexity, along with a larger number of predicted and transported variables (aerosol species and prognostic snow and rain).

The EAM is the atmosphere component of E3SMv1. It is a descendant of the Community Atmosphere Model version 5.3 (CAM5.3). EAM uses a spectral element dynamical core at 110‐km resolution on a cubed sphere geometry. Vertical resolution was increased from 30 layers, with a top at approximately 40 km, in CAM5 to 72 layers with a top at approximately 60 km in EAM. EAM contains many innovations compared to CAM5. While changes in the EAM physics are broadly similar to changes from CAM5.3 to CAM6, many details of tuning, and cloud and aerosols formulations differ in important ways as described in Rasch et al. ( 2019 ). A detailed analysis of its cloud and convective characteristics and the rationale for model tuning are provided in Xie et al. ( 2018 ) and Zhang et al. ( 2019 ). Sensitivity of EAM to a number of its adjustable parameters is explored in Qian et al. ( 2018 ). EAM simulated diurnal cycle of precipitation is analyzed in Xie et al. ( 2019 ).

The project also built a comprehensive infrastructure for code management, development, testing, and analysis to enable development of E3SMv1 and future versions at DOE leadership computing centers. Leveraging DOE investments, a flexible framework provides workflow orchestration, provenance capture and management, simulation analysis and visualization, and automated testing and evaluation capabilities (see Appendix Appendix C for a description of some of the analysis tools).

The subsections below provide a more detailed description of the model components. Coupling of the new MPAS‐Ocean, MPAS‐Seaice, and MOSART models with EAM and ELM provides E3SMv1 with a unique capability for multiresolution modeling using unstructured grids in most of its component models. This capability is critical for future simulation campaigns that have a strong regional focus to meet DOE's needs for Earth system modeling in support of energy‐sector decisions.

E3SM started with a version of the CESM1 (Hurrell et al., http://www.cesm.ucar.edu/models/cesm1.0 ) from which we developed the fully coupled E3SMv1 system. Notable changes between E3SMv1 and CESM1 include the following:

The entire simulation campaign was performed on NERSC Edison. E3SMv1 did not experience any internal failures during the entire 2,930 simulated years. The only failures were system related (generally node failures or file system issues), and the model could always be restarted from the last available set of annual restart files. The source code git hash for the piControl simulation was 2e145acf( https://github.com/E3SM-Project/E3SM/commit/2e145acf ) and for the remaining simulations was 7de18fc7 ( https://github.com/E3SM-Project/E3SM/commit/7de18fc7 ). The difference in source code version arises from code modifications necessary to support transient simulations that were not in place when the piControl simulation was started. The newer 7de18fc7 code bit‐for‐bit reproduces 2e145acf for piControl . A maintenance branch (maint‐1.0; https://github.com/E3SM-Project/E3SM/tree/maint-1.0 ) has also been specifically created to reproduce these simulations. Bit‐for‐bit results on NERSC Edison will be maintained on that branch for as long as the computing environment supports it.

Table 2 summarizes the E3SMv1 simulation campaign. All simulations were configured to adhere to the CMIP6 DECK specifications (Eyring et al., 2016 ) as closely as possible (see Appendix Appendix B for details about input data). piControl spans a total of 500 years. Initial conditions for the idealized 1%/year CO 2 increase ( 1pctCO2 ) and abrupt CO 2 quadrupling ( abrupt‐4xCO2 ) simulations, as well as the first member of the historical simulations ( historical_H1 ) were taken from 1 January of year 101 from piControl . Subsequent members were branched every 50 years for a total of five ensemble members. Atmospheric Model Intercomparison Project (AMIP) simulations (prescribed SST) were also performed to cover the entire period for which CMIP6 provides surface boundary conditions (1870–2014). Atmosphere and land initial conditions for amip_A1 were taken from year 1870 of historical_H1 , and similarly for the other members.

At that point, a decision was made to incorporate a minor code change in the atmospheric dynamical core. This change was a workaround for a compiler bug ( https://github.com/E3SM-Project/E3SM/pull/1922 ). In the interest of time, two new simulations were started simultaneously from the end of beta3rc13, one with identical tuning (“v1a”) and one with a further retuning of the clubb_c14 parameter from 1.08 to 1.06 (“v1b”). After 110 years, both simulations were examined and “v1b” was selected as the official E3SMv1 tuning because it exhibited the smaller drift of the two.

The first of these simulations consisted of 50 years with a prototype E3SMv1 configuration (“FCT2”). Due to the change in forcing from present‐day CORE‐II to preindustrial coupled, the simulation experienced a rapid cooling during the first 20 years followed by a relative stabilization. Ocean and sea ice states after 50 years served as initial conditions for a new coupled simulation (“beta3rc10”) that included bug fixes and associated retuning described in Table 1 , full implementation of the CMIP6 piControl forcings, and a retuning of clubb_c14 from 1.3 to 1.2 based on a few previous tuning attempts with the “FCT2” model configuration. This beta3rc10 configuration was run for a total of 260 years; it exhibited excessive net influx of energy at TOA and a long‐term warm drift with averaged temperature approaching present‐day values near the end, even with preindustrial forcings. As a result, a new configuration (“beta3rc11”) was branched off after 150 years reusing all components states as initial conditions, but with a slight retuning of clubb_c14 (1.2 to 1.1). This new configuration was run for a total of 100 years. Initially, it cooled slightly from the impact of the retuning, then drifted warm, albeit at a smaller pace than its predecessor. The process was repeated once again with “beta3rc13,” further reducing clubb_c14 from 1.1 to 1.08. This last simulation was run for a total of 80 years. The drift was smaller than the preceding attempt, but 80 years was too short to conclude that the drift had been removed entirely.

Since the coupled ocean and sea ice system take centuries to spin‐up, we chose to accelerate the process by simultaneously performing model spin‐up and final coupled tuning. MPAS‐Ocean was initially spun up for one year in stand‐alone mode, from the Polar science center Hydrographic Climatology data set (Steele et al., 2001 ) from rest. During this spin‐up, sea surface temperature (SST) and salinity were relaxed to an annual mean climatology (from Polar science center Hydrographic Climatology) and the currents were forced by an annual averaged CORE‐II normal year wind stress (Large & Yeager, 2009 ). This was followed by a 10‐year Coordinated Ocean‐ice Reference Experiments (CORE‐II), interannually varying, forced (Large & Yeager, 2009 ) ocean and sea ice simulation. The ocean and sea ice state at the end of that simulation served as initial conditions for a series of sequential fully coupled simulations with preindustrial forcing (Figure 2 ).

While most of the coupled tuning was performed under preindustrial conditions, we also performed a few additional simulations to monitor climate sensitivity and total ERF during the developmental phase of E3SMv1. We evaluated sensitivity using idealized +4‐K simulations (Cess et al., 1989 ) and one abrupt quadrupling of CO 2 . Total ERF was estimated using pairs of atmosphere simulations with identical SSTs but differing forcing (climatological 2000 vs. 1850). One historical test simulation was also performed with a near final version of E3SMv1 (“beta3rc10” below). The outcomes of these intermediate simulations were consistent with the final model. A pragmatic, but deliberate, decision was made not to attempt to reduce the high climate sensitivity or reduce the aerosol forcing.

This final tuning required adjustments of less than 1 W/m 2 because the E3SM atmosphere component was already well tuned from simulations with prescribed SSTs and sea ice concentrations. To simplify the process, we chose to focus on a single parameter in the CLUBB atmosphere turbulence parameterization, clubb_c14, which directly impacts dissipation of the horizontal components of the turbulent kinetic energy. Because of the tight coupling between turbulent kinetic energy and boundary layer clouds, clubb_c14 modulates low‐clouds, thus affecting net TOA mostly through SW cloud radiative effects.

Tuning is an integral part of climate model development (e.g., Hourdin et al.,; Schmidt et al.,). While most of the tuning for E3SMv1 was performed at the model component level, some additional retuning was subsequently required for two reasons: (1) Code errors were discovered in the atmosphere model physics, while the coupled system was being developed and (2) the coupled tuning objectives required an additional level of adjustment to accommodate biases present in other model components and interactions between components. Two minor code errors were discovered and addressed in EAMv1:

Additional problems were also identified in the fully coupled system, including nonconservation caused by inconsistent remapping in the exchange of water between model components, excessive storage capacity in the river routing model, missing perched drainage and ponding in the land surface model, and grid definition inconsistencies between different model components. Some of these errors caused nonconservation errors in sea level rise as large as tens of meters per century but all were resolved and corrected. Unfortunately, a smaller nonconservation error was not uncovered until after E3SMv1 was frozen and the DECK simulations were complete. It manifests itself by a steady loss of water in the ocean of 5 cm per century in the piControl simulation. We note that the water loss does not involve any phase change and therefore does not impact TOA energy balance; however, it does impact ocean heat content. The root cause has since been traced to a nonphysical update of water stored in the unconfined aquifer of urban subgrid land units within the land model. This error, which has been corrected for future versions of E3SM ( https://github.com/E3SM-Project/E3SM/pull/2603 ), was partly masked by incomplete and inconsistent water budget checks within the land model and the coupler. These budgets have also been corrected.

Water conservation is also important for a global Earth system model. During the early development of the E3SMv1 model, we identified a number of computational problems, leading to water conservation errors in the atmosphere component (Zhang et al., 2018 ). These problems were all considered and suitable remedies applied (e.g., borrowing from adjacent cells to avoid more drastic nonconservative fixers) so that in the 500‐year piControl simulation, the relative water conservation error (as defined in Zhang et al., 2018 ) in the atmosphere component is about 0.00226%, equivalent to a computational sea level rise of about 2 mm per century.

Among 42 CMIP5 models analyzed by Hawkins and Sutton ( 2016 ), a majority of them have preindustrial global mean surface temperatures between 13 and 14 °C, with some as low at 12 °C and as high as 15 °C. With a long‐term average surface temperature of 13.8 °C, E3SMv1 falls within the range of the bulk of CMIP5 models (tuning objective 3). This value is also consistent with estimated warming and the present‐day global temperature of 14.0 ± 0.5 °C by Jones et al. ( 1999 ) for the period 1961–1990 and with leading reanalyses data sets (14.3 to 14.6 °C) for the period 1979–2008 (Hawkins & Sutton, 2016 ).

The global mean surface air temperature is very stable over the course of the control simulation (tuning objective 2; Figure 3 b). There is a slight positive temperature trend of +0.013 °C/century over the entire 500 years and a much smaller trend of −0.003 °C per century over the last 400 years. We note that all additional simulations were branched from piControl in year 101 or later (Table 2 ), thus ensuring that forcing perturbations are applied atop a nondrifting control simulation. This facilitates subsequent analyses, since it is not necessary to correct for a long‐term drift in the underlying control. Another measure of the stability of the climate in the piControl simulation is provided by the annual maximum and minimum sea ice area in Arctic and Antarctic (Figure 3 c). Both hemispheres have little long‐term drift in their seasonal cycle of sea ice extent.

Time evolution of annual (a) global mean net top‐of‐atmosphere (TOA) radiation (positive down), (b) global mean surface air temperature, and (c) maximum and minimum of total sea ice area for the Arctic and Antarctic in the piControl simulation. Dashed lines in (a) and (b) represent linear trends. The solid straight line in (a) is the mean TOA energy imbalance of −0.011 W/m 2 .

The preindustrial control simulation ( piControl ) is in energy balance at TOA with an average loss of 0.011 W/m 2 over the course of 500 years (tuning objective 1; Figure 3 a) and almost no long‐term linear trend. Among all the model components, the ocean constitutes the largest reservoir of heat. It takes up heat at an average rate of 0.016 W/m 2 , leaving a small net imbalance of 0.027 W/m 2 , either from changes in other components or energy nonconservation. Since this imbalance is sufficiently small compared to anthropogenic forcings of interest, E3SMv1 can be regarded as essentially conserving energy. We note, however, that developmental versions of E3SMv1 suffered from much larger imbalances (on the order of 0.5 W/m 2 ). That imbalance was caused by inconsistent definitions of energy in the ocean and the atmosphere, with the ocean properly accounting for changes in water heat content with temperature while the atmosphere did not. The 0.5‐W/m 2 imbalance was deemed too large to ignore, but rewriting the atmosphere physics more consistently was impractical due to time constraints. As a result, we decided to incorporate an ad hoc correction term. An additional surface (sensible heat) source is introduced to the atmosphere that accounts for the missing energy carried by water molecules leaving the ocean at one temperature and returning (as condensed water) at another temperature. The apparent flux is small (0.4 W/m 2 ) and applied uniformly at every point globally. The correction does not impact the simulated climatology, but its neglect would had led to a long‐term drift in the global mean temperature of a few tens of degree by century. With this correction, both atmosphere and ocean communicate consistent energy loss and gain with each other (see Appendix Appendix A for more detail).

5 Historical Simulations

5.1 Atmosphere Climatology We first evaluate the atmosphere in E3SMv1 with climatologies from the final 30 years (1985–2014) of simulations. Five ensemble members of the historical coupled experiment (H1 to H5) and three ensemble member AMIP simulations (A1 to A3) are evaluated. The net TOA radiative flux (positive down) is shown in Figure 4 compared to CERES‐EBAF Ed4.0 (Loeb et al., 2009). Bias patterns are qualitatively similar between the atmosphere and coupled simulations, indicating a strong imprinting of atmosphere biases on the coupled system, but the biases are frequently larger in the coupled simulations, where component interactions play a role (root‐mean‐square error [RMSE] of 9.13 vs. 7.80 W/m2). Significant positive biases are seen in the subtropical stratocumulus regions off the west coasts of North America, South America, and Africa due to an underestimate of cloudiness. The underpredicted coastal stratocumulus is due to both the use of CLUBB and model retuning (Xie et al., 2018). Southern Oceans are also marked by positive biases due to clouds. The remaining oceanic regions have either neutral or negative biases. Tropical land masses have generally positive biases, while northern high latitudes suffer from negative ones. The global mean net flux is +0.32 W/m2 for the coupled simulations and +0.15 W/m2 for the AMIP simulations, both slightly smaller than the observational estimates of 0.86 W/m2 from CERES‐EBAF Ed 4.0. (uncertainty range of 0.2 to 1.0; Wild et al., 2012). We note that the present‐day energy imbalance was not a target of model tuning. Instead, we tuned the preindustrial control simulation to have a net zero global mean TOA flux. The imbalance for 1985–2014 emerges from the model evolution in response to historical forcings in the case of the coupled simulations or as a result of imposing present‐day SST and sea ice concentrations (as boundary conditions) for the AMIP simulations. Figure 4 Open in figure viewer PowerPoint Annual net top‐of‐atmosphere (TOA) radiative flux: (a) CERES‐EBAF Ed4.0 observational estimate, (b) model bias from the three Atmospheric Model Intercomparison Project (AMIP) ensemble simulations, and (c) model bias from the five ensemble historical coupled simulations (1985–2014). RMSE = root‐mean‐square error. Biases in cloud radiative effects (CRE) at TOA compared to CERES‐EBAF Ed4.0 (Loeb et al., 2009) are shown in Figure 5 for SW and LW. Many regional biases apparent in net TOA (Figure 4) can be traced to SW cloud biases. LW CRE reveal additional biases, such as a lack of LW cloud trapping in the tropics from high clouds and excessive LW trapping from northern high latitudes clouds. As indicated in Xie et al. (2018), high clouds are significantly reduced in the tropical deep convection regions due to the increase of model vertical resolution from 30 levels to 72 levels in EAM, which results in a much weaker LW CRE over these regions. Global mean CRE for both SW and LW are approximately 3 W/m2 below the observational estimates. Figure 5 Open in figure viewer PowerPoint Annual top‐of‐atmosphere cloud radiative effect model biases (W/m2) of the five ensemble historical simulations (1985–2014) against CERES‐EBAF Ed4.0: (a) shortwave (SW) and (b) longwave (LW). RMSE = root‐mean‐square error. Annual precipitation is depicted in Figure 6 compared to Global Precipitation Climatology Project (GPCP) v2.2 (Adler et al., 2003; Huffman et al., 2009). Simulated global mean precipitation is slightly above 3 mm/day for both AMIP and coupled simulations, approximately 15% larger than the GPCP estimate. A comprehensive review of 30 currently available global precipitation data sets found large differences in the magnitude of global land annual precipitation estimates (Sun et al., 2018). In their estimate of the global energy cycle, Wild et al. (2012) note that the GPCP estimate may be too low due to systematic underestimations in the satellite retrievals (Stephens et al., 2012; Trenberth et al., 2009). Wild et al. (2012) estimated the global mean precipitation at 2.94±0.17 mm/day, while the estimate from Stephens et al. (2012) is 3.04 ± 0.35 mm/day. The global mean precipitation from E3SMv1 falls within both estimates. Figure 6 Open in figure viewer PowerPoint Annual precipitation rate (mm/day): (a) Global Precipitation Climatology Project v2.2 observational estimate, (b) model bias from the three Atmospheric Model Intercomparison Project (AMIP) ensemble simulations, and (c) model bias from the five ensemble historical coupled simulations (1985–2014). RMSE = root‐mean‐square error. Major regional precipitation biases in AMIP simulations (Figure 6b) include a dry Amazon, excessive precipitation over elevated terrain (e.g., Andes and Tibetan Plateau), wet biases over tropical Africa and the Indian Ocean, and a dry bias over the central United States (traced to the summer time). Coupled simulations (Figure 6c) tend to amplify regional biases present in AMIP simulations (except the central United States), as well as develop biases typical of coupled simulations: double Intertropical Convergence Zone and excessive precipitation over the maritime continent. Unsurprisingly, the RMSE increases from 0.93 to 1.13 mm/day and the correlation decreases (0.91 vs. 0.86) between AMIP and coupled simulations. Figure 7 shows the annual, zonally averaged temperature for the coupled simulations compared to ERA‐Interim reanalysis (Dee et al., 2011). Overall, the model captures the thermal structure of the atmosphere. There is a tropical cold bias in the upper troposphere, which correlates with insufficient tropical high clouds (Figure 5b). In the northern midlatitudes, the cold bias extends to the surface because of the colder than observed northern latitudes SSTs (see section 5.2 below). There is also a significant warm bias in the lower portion of the troposphere in the southern high latitudes. The corresponding zonal wind is shown in Figure 8. The locations of the maximum jet locations are displaced equatorward, and the wind magnitude is also too large throughout the midlatitude troposphere. Tropical easterlies are too weak in the upper troposphere and too strong near the surface. Figure 7 Open in figure viewer PowerPoint Annual zonally averaged temperature (K): (a) ensemble mean of historical coupled simulations (1985–2014), (b) ERA‐Interim reanalysis, and (c) model bias. RMSE = root‐mean‐square error. Vertical coordinate is pressure level (hPa). Figure 8 Open in figure viewer PowerPoint Annual zonally averaged zonal wind (m/s): (a) ensemble mean of historical coupled simulations (1985–2014), (b) ERA‐Interim reanalysis, and (c) model bias. RMSE = root‐mean‐square error. Vertical coordinate is pressure level (hPa). The evaluation above provides only a limited view of the performance of E3SMv1 from an atmospheric perspective. For a more exhaustive evaluation, we turn to a comparison with an ensemble of 45 CMIP5 models using metrics computed with the PCMDI Metrics Package (PMP; Gleckler et al., 2008, 2016). The comparison covers the period 1981–2005 of the historical simulations. The historical and AMIP ensemble members of E3SMv1 were also processed with PMP. Figure 9 shows global RMSE for the CMIP5 ensemble using box and whiskers plot, as well as the individual E3SMv1 historical simulations with blue dots and AMIP with red dots. Spatial RMSE against observations collected by PMP are shown for nine fields and each one of them for annual and seasonal averages. Lower values are better. For TOA radiation fields (Figures 9a–9c), E3SMv1 coupled (blue) generally falls within the lowest (best) quartile and is even competitive with some of the best CMIP5 models for certain fields and seasons. We note that we are comparing a newer model against older ones, so we do not expect this to necessarily hold for CMIP6. For surface variables, precipitation (Figure 9d), surface air temperature over land (Figure 9e), and zonal wind stress over ocean (Figure 9f), E3SMv1's coupled performance is better than the ensemble median and often falls within the lowest quartile, with the exception of surface air temperature over land during December–February and March–May (MAM). Precipitation during MAM is also notably worse than other seasons relative to the CMIP5 ensemble. For dynamical variables, zonal wind at 200 and 850 hPa (Figures 9g and 9h) and 500‐hPa geopotential height (Figure 9i), E3SMv1 is generally better than the CMIP5 median, except again for MAM. Figure 9 Open in figure viewer PowerPoint 2016 Comparison of RMSE (1981–2005) of an ensemble of 45 Coupled Model Intercomparison Project Phase 5 models (box and whiskers showing 25th and 75 percentiles, minimum and maximum) with the five E3SMv1 historical members (blue dots) and the first member of the Atmospheric Model Intercomparison Project simulations (ret dots). Spatial RMSE against observations are computed for annual and sesonal averages with the PCMDI Metrics Package (Gleckler et al.,). Fields shown include TOA net radiation (a), TOA SW and LW cloud radiative effects (b, c), precipitation (d), surface air temperature over land (e), zonal wind stress over ocean (f), 200‐ and 850‐hPa zonal wind (g, h), and 500‐hPa geopotential height (i). TOA = top‐of‐atmosphere; SW = shortwave; CRE = cloud radiative effects; LW = longwave; DJF = December–February; MAM = March–April; JJA = June–August; SON = September–November; RMSE = root‐mean‐square error. Unsurprisingly, AMIP simulations (red) perform better than their coupled counterparts (blue). However, the relative degradation between AMIP and coupled helps attribute sources of errors. For example, the difference is relatively small for the zonal mean wind, and thus, improving overall performance would likely require atmospheric improvements. On the other hand, surface variables, in particular precipitation and temperature, are much more strongly affected by the coupled model errors that emerge. The MAM seasonal deficiency also appears to be rooted in coupling errors. In summary, E3SMv1 performs better than the median of the of CMIP5 ensemble for most atmospheric fields and seasons, which helps establish the credibility of the E3SMv1 simulated climate.

5.2 Ocean Climatology An annual climatology (1985–2014) of the E3SMv1 ensemble mean of SST is shown in Figure 10, compared to the Hadley‐National Oceanic and Atmospheric Administration (NOAA)/OI merged data product (Hurrell et al., 2008) averaged over the same period. Overall, E3SMv1 captures the observed SST well, with an ensemble mean bias of 0.094 °C and an ensemble mean RMSE of 0.939 °C. A few biases do emerge. First, there is a cold bias in the North Atlantic associated with excessive sea ice in the Labrador Sea. Consistent with this sea ice bias, these cold SST biases are stronger in the first half of the year than in the second half (not shown). While the exact cause of the Labrador Sea ice bias is unknown, it is likely that missing critical heat transports (e.g., from the east and west Greenland currents, Irminger current, and the Northwest corner), perhaps from low resolution or excessive horizontal viscosity (Jochum et al., 2008) play an important role. Second, there are warm biases on the eastern sides of ocean basins, coincident with SW CRE biases discussed in section 5.1. Finally, the Southern Ocean is also ≈2 °C warmer than observed. It is possible this last bias is associated with the relatively large value for the Gent‐McWilliams bolus coefficient used, which was chosen to prevent excess heat transport to the deep ocean (see section 2.2). Figure 10 Open in figure viewer PowerPoint 2008 mean, max are listed for the historical ensemble. Annually averaged sea surface temperature climatology (1985–2014) for (a) E3SMv1 historical simulation (ensemble mean), (b) Hadley‐National Oceanic and Atmospheric Administration/OI merged sea surface temperature data set (1985–2014; Hurrell et al.,), and (c) model bias. In the upper right corner of panel (c), the mean bias is shown in square brackets and the root‐mean‐square error is in angular brackets. For each error, the min,, max are listed for the historical ensemble. The E3SMv1 sea surface salinity bias is shown in Figure 11. The model output is compared to the 2011–2014 National Aeronautics and Space Administration (NASA) Aquarius data (Lagerloef et al., 2015). Overall, the surface ocean is too fresh, with the largest bias in the Labrador Sea. The two minor exceptions to this pattern are the positive salinity bias region over the west Pacific warm pool, which is coincident with the atmospheric precipitation maximum being shifted westward relative to observations (Figure 6) and a positive salinity bias in the Arctic. Figure 11 Open in figure viewer PowerPoint As in Figure 10 but for sea surface salinity (PSU). The data set in (b) is from National Aeronautics and Space Administration Aquarius data (averaged 2011–2014). The ensemble average mixed layer depth (MLD) annual climatology, based on a critical density threshold (σ c = 0.03 kg/m3), is presented in Figure 12. The model output is compared to data described by Holte et al. (2017). The globally averaged model MLD is too shallow relative to observations, with the largest bias coincident with the large fresh bias in the North Atlantic. Unlike many other CMIP5 models (e.g., Sallée et al., 2013), the E3SMv1 MLD in the Southern Ocean is slightly deeper than observed. This is due to the region of deeper mixed layers in E3SMv1 being broader than observed (Figure 12a), possibly due to a positive bias in the Southern Ocean wind stress. The maximum MLD simulated by E3SMv1 is still much shallower than observations, which is consistent with other CMIP5 models. Figure 12 Open in figure viewer PowerPoint 2017 As in Figure 10 but for the annual average mixed layer depth. The data in (b) are from Holte et al. (). Data have been averaged from 2001 to 2017. Figure 13 shows the E3SMv1 maximum Atlantic Meridional Overturning Circulation (AMOC) at 26°N, the site of the RAPID array (Smeed et al., 2017). The mean AMOC simulated in the historical ensemble (approximately 11 Sv) is weaker than the observed mean (16.9 ± 3.35 Sv) and also on the weak end of CMIP model AMOC strength (e.g., Cheng et al., 2013, their Figure 1). There are a number of possible causes for the weak AMOC in E3SMv1. Currently, MPAS‐Ocean utilizes a z star coordinate (Adcroft & Campin, 2004), which is broadly consistent with a traditional z level coordinate in the deeper ocean. Z coordinate models are known to experience spurious diapycnal mixing (e.g., Griffies et al., 2000), which could reduce AMOC strength. Second, recent work has shown that Nordic overflows (Wang et al., 2015) and Arctic freshwater transports (Wang et al., 2018) have a strong impact on AMOC. At low resolution, these processes are poorly represented, where the latter is likely too strong as critical passageways being too wide (e.g., Davis Strait). Finally, since E3SMv1 exhibits excess sea ice in the Labrador sea, the simulated MLD (Figure 12) is reduced, which reduces deep convection and hence AMOC strength. Figure 13 Open in figure viewer PowerPoint Annually averaged maximum Atlantic Meridional Overturning Circulation (MOC) at 26.5°N below 500‐m depth. The solid black line represents the ensemble mean and the spread of the shaded gray region represents the maximum and minimum of the ensemble. The observed AMOC and standard deviation at the RAPID array are shown by the vertical bar (16.9 ± 3.35 Sv.).

5.3 Sea Ice Sea ice in E3SMv1 is too extensive at the end of winter and too confined in late summer, reflecting a heightened model seasonality relative to both the Arctic and Southern Ocean passive microwave record, regardless of the algorithm used to derive ice area from these retrievals (e.g., Cavalieri et al., 1996; Meier et al., 2017). This result is summarized in Figures 14-16. Embedded within this heightened seasonality, the melt season minimum is delayed in the Northern Hemisphere relative to observations. In the high north, the heightened winter seasonal extent relative to observations occurs mainly in the Labrador Sea, as well as the Iceland Sea and Pacific margin of the Sea of Okhotsk, evident in the historical five‐member ensemble means (Figures 14a and 14b). Comparison with satellite‐derived albedo, shown for June to August in Figure 17, suggests that surface radiative processes leading to an albedo bias in the central arctic are unrelated to the physical exchanges responsible for the Labrador sea bias. As a result, the Labrador Sea bias is unlikely to be attributable to MPAS‐Seaice alone and is more likely a result of coupled interactions and unresolved heat advection in the region. Figure 14 Open in figure viewer PowerPoint 2017 Northern Hemisphere ensemble mean March and September sea ice thickness for two decades leading up to year 2000 (a and c), and the first 15 years of the 21st century (b and d). Model ice thickness is truncated at 15% concentration, and magenta represents the Meier et al. () National Oceanic and Atmospheric Administration Climate Data Record ice extent for the same averaging periods. Grid density on the polar stereographic projection is indicated by cell translucence. Numbers in the lower right corner of each panel indicate the mean hemispheric ice thickness thickness for the thinnest ensemble member (blue), the multiensemble mean rendered in the map (black), and the thickest ensemble member (red) for the each period. The model possesses a negative climatological trend in Arctic sea ice thickness that is qualitatively consistent with observed basin‐wide trends (e.g., Kwok & Rothrock, 2009), as quantified in Figure 14. The 1980–1999 hemispheric mean spread among the five ensemble members for March (September) decreases from 1.5–1.71 to 1.16–1.37 m (1.55–2.12 to 0.79–1.35 m) for the period 2000–2014, which includes the thinnest periods known in the Arctic sea ice thickness record. Still, the spatial ice thickness pattern in the central Arctic does not reflect the predominant buildup of thick ice against the Canadian Archipelago seen in observational estimates, including those of Bourke and McLaren (1992), Kwok and Cunningham (2008), and Tilling et al. (2015). This represents a second sea ice bias being addressed in ongoing improvements to the polar components of E3SM. There is no significant climatological trend in mean Southern Ocean sea ice thickness from the historical simulations, in the ensemble mean or spread (Figure 15). Extent bias stemming from the heightened polar seasonality of E3SMv1 is consistent across our chosen 1980–1999 and 2000–2014 analysis periods, seen in Figure 15. While the model does not replicate the observed climatological trend in austral sea ice extent (Figure 16), it does simulate the critical decrease in Northern Hemisphere minimum sea ice extent, which has had a significant impact on planetary albedo in the current century and is therefore an important contribution to the energy balance of E3SMv1 as a whole. Figure 15 Open in figure viewer PowerPoint As for Figure 14 but for the Southern Hemisphere. Figure 16 Open in figure viewer PowerPoint 1996 Maximum (a) and minimum (b) sea ice area observed during the year for the Northern (green) and Southern Hemispheres (blue). The shaded bounds represent the E3SMv1 historical ensemble spread, and the thick colored lines are National Aeronautics and Space Administration TEAM observations (Cavalieri et al.,). Model and observational data are monthly averages. Figure 17 Open in figure viewer PowerPoint Surface albedo bias (E3SMv1‐CERES‐EBAF) for Northern Hemisphere spring, averaged over 1985–2014.

5.4 Land and River Figure 18 shows a comparison of the mean annual total (surface and subsurface) runoff simulated by ELM with the composite runoff map from the Global Runoff Data Center (GRDC; Fekete & Vörösmarty, 2011). Since the GRDC data only provide monthly runoff for 1986–1995, the comparison is shown for the annual total runoff averaged over the same period of the historical_H1 simulation. ELM captures the general spatial distribution of the GRDC runoff, but in relatively arid regions such as Australia and the western United States, ELM has wet biases, while in the Amazon tropical forest, dry biases are notable. These biases are consistent with the precipitation deficiencies shown in Figure 6. Figure 18 Open in figure viewer PowerPoint Mean annual total runoff (mm/day) simulated by (upper panel) E3SMv1 land model (ELM) and from (lower panel) Global Runoff Data Center (GRDC) averaged between 1986 and 1995. (1) (2) Q ij is the monthly streamflow for month i and year j and n is the number of years in the simulation or observation. With this definition, SI is equal to 1 if the monthly streamflow is uniformly distributed throughout the year and SI is equal to 12 if streamflow only occurs in 1 month. The seasonality of the simulated streamflow is generally comparable with that observed in terms of both magnitude and timing. For example, in North America, streamflow seasonality is stronger in the northwest with a peak timing between November and January but in the central and southeastern United States, streamflow seasonality is weaker with a peak timing generally in spring. In Asia, streamflow generally peaks in the late summer. Larger biases in seasonality are found in Australia (Murray Darling River) and central Asia (Yenisey) where biases in the runoff are also more significant (Figure The seasonal cycle of streamflow is an important metric of water availability, as water deficits resulting from a mismatch in the timing of water supply and water demand have important implications for energy and water. Figure 19 compares the seasonal cycle of simulated and observed streamflow at stream gauge stations in major river basins around the world. A seasonality index (SI) is defined as follows:whereis the monthly streamflow for monthand yearandis the number of years in the simulation or observation. With this definition,is equal to 1 if the monthly streamflow is uniformly distributed throughout the year andis equal to 12 if streamflow only occurs in 1 month. The seasonality of the simulated streamflow is generally comparable with that observed in terms of both magnitude and timing. For example, in North America, streamflow seasonality is stronger in the northwest with a peak timing between November and January but in the central and southeastern United States, streamflow seasonality is weaker with a peak timing generally in spring. In Asia, streamflow generally peaks in the late summer. Larger biases in seasonality are found in Australia (Murray Darling River) and central Asia (Yenisey) where biases in the runoff are also more significant (Figure 18 ). Figure 19 Open in figure viewer PowerPoint Seasonality of observed (upper) and simulated (lower) streamflow at stream gauge stations in major river basins around the world. The magnitude of seasonality indicated by the seasonality index (SI) is shown in color, and the timing of the peak streamflow is shown by the position of the arrows. SI is calculated based on monthly streamflow between 1986 and 1995.

5.5 Variability Simulated by E3SMv1 In this section we present a variety of metrics to assess E3SMv1 variability. The long‐term variability of the El Niño–Southern Oscillation (ENSO), as simulated in the piControl and historical ensemble simulations, is examined via wavelet analysis (Torrence & Compo, 1998) of Niño 3.4 SST and is shown in Figure 20. In this analysis, the piControl simulation is subdivided into five 100‐year‐long sections (as in Stevenson, 2012). For reference, we also include the spectrum from HadISST (Rayner et al., 2003) and ERSSTv4 (Huang et al., 2015). In both the piControl and historical ensemble, E3SMv1 ENSO variability is strong with statistically significant peaks near a 3‐year period. We also see a signature of longer term modulation of ENSO variability (6‐ to 9‐year period) in the piControl. This is consistent with other CMIP models (e.g., Stevenson, 2012; Stevenson et al., 2012; Wittenberg, 2009) and proxy data (e.g., McGregor et al., 2013). However, the modulation is a bit weaker than other models. Relative to the CESMv1 Large Ensemble (Kay et al., 2015), Figure 20b, where the CESMv1 PI control is subdivided as in E3SMv1, E3SMv1 variability is slightly closer to observations, but is shifted strongly to a three year period, whereas CESM‐LE is dominant at approximately 4.5 years. In both models, the dominant ENSO peak remains consistent between the PI control and historical ensembles. We also note that the spread of ENSO variability in the E3SMv1 historical simulation is much smaller than in the CESM‐LE. This is likely due to the small E3SMv1 ensemble size relative to the 39 member CESM‐LE (Newman et al., 2018). Figure 20 Open in figure viewer PowerPoint 1998 2012 El Niño–Southern Oscillation (ENSO; Nino3.4) variability of the pre‐industrial (PI) control simulation and historical ensemble. The Morlet wavelet of degree 6 is used (e.g., Torrence & Compo,) and the maximum and minimum wavelet power is used for each period. The PI control is subdivided into 100‐year intervals to form a five‐member ensemble (e.g., Stevenson,). The solid line represents the mean wavelet power for the ensemble. The shading bounds the maximum and minimum power that is above the 90% significance threshold. The black and green lines are two observational data products. (a) E3SMv1 and (b) CESM‐LE. SST = sea surface temperature; CESM = Community Earth System Model. A number of CMIP5 models do not well represent the spatial pattern of ENSO variability, in particular its westward extent (e.g., Menary et al., 2018; Van Oldenborgh et al., 2005). Figure 21 shows the difference of a composite of the strongest El Niño events and strongest La Niña events (following Menary et al., 2018). Broadly, E3SMv1 reproduces the spatial pattern seen in observations well with a few notable differences. E3SMv1 does not capture the signal along the North American coast, suggesting a bias in the coastally trapped Kelvin waves. The westward extent is also larger than observed but better than seen in CESM‐LE. Overall, the comparison with observations is good, with low RMSE and a high correlation coefficient. Figure 21 Open in figure viewer PowerPoint 2018 (a–c) Difference of composite El Ninõ events and composite La Ninã events (1950–2015) for the HadleyISST data set, the E3SMv1 historical ensemble, and the CESM‐LE, respectively. El Ninõ events are defined as periods when the Ninõ 3.4 SST anomaly exceeds 0.8 °C for more than six consecutive months. The La Ninã criterion is Ninõ 3.4 SST anomaly less then −0.8°C for more than 6 months (these definitions are consistent with Menary et al.,). When an El Niño–Southern Oscillation event is identified, the SST is averaged from November to March. For model output, every ensemble member contributes to the mean composite. The Pearson correlation coefficient and RMSE are shown for E3SMv1 and CESM in (b) and (c). SST = sea surface temperature; CESM = Community Earth System Model; RMSE = root‐mean‐square error. On subseasonal time scales the dominant mode of variability and predictability in the tropics is the Madden‐Julian Oscillation (MJO; Waliser et al., 2003). The MJO is generally thought to play a role in ENSO initiation (McPhaden et al., 2006), monsoon active break cycles (Annamalai & Slingo, 2001), tropical cyclogenesis (Sobel & Maloney, 2000), and remote teleconnection effects (Vitart, 2017); therefore, its accurate simulation is key. The simulation of the MJO in E3SMv1 represents a significant improvement in strength, propagation characteristics, and the explained intraseasonal variance compared to CESM1 (Figure 22). However, significant biases in Pacific propagation remain. In CESM1, the intraseasonal propagation is in the wrong direction, westward, in the Indian Ocean. Weak correlations do make it over the Maritime Continent and into the West Pacific, but they are mostly decoupled from the SST signal. In contrast, E3SMv1 has consistent low‐level wind propagation coupled in quadrature with the SSTs from the Indian Ocean to the central Pacific. Figure 22 Open in figure viewer PowerPoint Tropical lag correlation (averaged 10°N to 10°S) of precipitation (colors) and 850‐mb zonal wind (lines) with precipitation in the Indian Ocean region (60°N to 90°E; shown by the vertical dashed green lines) for (a) observed (Tropical Rainfall Measuring Mission for precipitation Modern‐Era Retrospective Analysis for Research and Applications [MERRA] for winds), (b) E3SMv1 preindustrial control, and (c) Community Earth System Model (CESM) preindustrial control. Data used are daily anomalies and band‐pass filtered between 20 and 100 days.