E3SMv0‐HiLAT: A Modified Climate System Model Targeted for the Study of High‐Latitude Processes

We document the configuration, tuning, and evaluation of a modified version of the Community Earth System Model version 1 (Hurrell et al., 2013, https://doi.org/10.1175/BAMS-D-12), introduced here as E3SMv0‐HiLAT, intended for study of high‐latitude processes. E3SMv0‐HiLAT incorporates changes to the atmospheric model affecting aerosol transport to high northern latitudes and to reduce shortwave cloud bias over the Southern Ocean. An updated sea ice model includes biogeochemistry that is coupled to an extended version of the ocean model's biogechemistry. This enables cloud nucleation to depend on the changing marine emissions of aerosol precursors, which may be expected in scenarios with strongly changing sea ice extent, oceanic stratification and associated nutrient availability, and atmospheric state. An evaluation of the basic preindustrial state of E3SMv0‐HiLAT is presented in order to ensure that its climate is adequate to support future experimentation. Additional capability is not achieved without some cost, relative to the extraordinarily well‐tuned model from which it was derived. In particular, a reduction of bias in cloud forcing achieved over the Southern Hemisphere also allows for greater Southern Ocean sea ice extent, a tendency that has been partially but not fully alleviated through experimentation and tuning. The most interesting change in the behavior of the model may be its response to greenhouse gas forcing: While the climate sensitivity is found to be essentially unchanged from that of Community Earth System Model version 1, the adjusted radiative forcing has increased from within one standard deviation above that of Coupled Model Intercomparison Project Phase 5 models to nearly two standard deviations.


Introduction
The high latitudes, where snow and ice coverage persist over much of the year, have seen especially large changes in their climate. In the Arctic, the regional rise in temperature has been approximately twice the global average, based on annual mean estimates (Serreze et al., 2009;Walsh, 2014). A dramatically increased seasonal cycle of sea ice extent means that much more open water is now exposed to the atmosphere during summer and autumn seasons (Comiso, 2012), and loss of multiyear ice leaves much of the year-round Arctic sea ice vulnerable to further retreat (Peng & Meier, 2017).
The Antarctic, more isolated in terms of its meteorology and heavily influenced by ocean currents, has not experienced such region-wide polar amplification of atmospheric warming. Annual maximum sea ice extents had actually increased until recently (Comiso et al., 2017;Simmonds, 2015). But changing conditions have led to loss of major portions of the ice shelves along the more exposed Antarctic Peninsula (Cook et al., 2014), and it has been suggested that the Thwaites and Pine Island Glaciers have passed tipping points in terms of their retreat (Joughin et al., 2014;Rignot et al., 2014), with great consequence for global sea level over the longer term.
These remarkable changes at the high latitudes pose challenges for Earth system models, including and going beyond the widely appreciated issues around the modeling of sea ice of the past 10 or 20 years, such as slow model loss of Arctic ice (Stroeve et al., 2007) and the difficulty models have in capturing the south-

Model Description
E3SMv0-HiLAT is a fully coupled Earth System Model based on the CESM1. The atmosphere component is the Community Atmosphere Model version 5 with the Spectral Element dynamical core (Dennis et al., 2012), the land component is the Community Land Model version 4.5 (Lawrence et al., 2011), the ocean component consists of the Parallel Ocean Program version 2 model (Smith, 2010), and the sea ice component is the CICE5 model . The code base is the E3SMv0.1 model, which in turn is based upon the CESM 1.3beta10 version (as used at higher resolution in Kurtakoti et al., 2018). Relative to the code that ran the CESM1 LENS simulations, E3SMv0.1 mostly differs in its atmospheric model component, which is CAM5.3 using the spectral element core (the CAM5.3 version includes a bug fix related to how a cloud radiation parameter was calculated). E3SMv0-HiLAT is significantly more expensive than CESM1, primarily due to our use of the spectral element core (subsequent performance improvements in the new CESM2 version of the core, as documented in Lauritzen et al., 2018, have largely eliminated the computational penalty of using this core, even at our relatively low resolution). Although this limits the number of years of simulation we were able to conduct, it was considered important to support the development path of subsequent versions of E3SM (Golaz et al., 2019;Petersen et al., 2019). Secondarily, although the errors associated with modern atmospheric dynamical cores tend to be smaller than those associated with the physical parameterizations, it was considered preferable for our purposes to avoid the filtering at polar latitudes required by the finite volume core used in LENS (the impact of this filtering at high latitudes is discussed by Lauritzen et al., 2010, andEvans et al., 2013).
An important research front in high-latitude climate modeling that we do not specifically address here is the suitability for use with an ice sheet model. While the Community Ice Sheet Model Lipscomb et al., 2019) may be coupled into E3SMv0-HiLAT and the initial work of Vizcaíno et al. (2013) is accounted for in our model, the subsequent work of van Kampenhout et al. (2017) to further improve the mass balance of CESM is not and so would require further effort to port to either E3SMv0-HiLAT or to the substantially different system described by Golaz et al. (2019), where much of the effort on ice sheet modeling has been aimed at support for the strong regional refinement allowed by the use of a Voronoi tessellation (Hoffman et al., 2018).
An issue we do take up here, standing as a notable deficiency of the current generation of climate models even as it represents a significant control in the climate system (e.g., Wang, Maltrud, Burrows, et al.,2018), is the impact of high-latitude ocean/sea ice ecosystems on cloud nucleation. Our effort aims to improve on this situation through changes in our model formulation that are detailed in the subsections below, including (i) modifications to the atmospheric model to allow for prognostic marine aerosol precursors and to reduce SWCF biases, (ii) replacing version 4 of the sea ice model with a more recent version that includes biogeochemistry, and (iii) extending the ocean model's biogeochemistry.
While our aim is to establish a greater capability for the study of high-latitude climate processes, we mean to do so without too much loss of fidelity elsewhere, relative to LENS. This will be a concern throughout section 3 but also to some extent in this section.

Atmosphere
The atmospheric component of E3SMv0-HiLAT is based on the standard CAM5.3 using the spectral element dynamical core. Following Zhang et al. (2018), a hybrid numerical coupling between the resolved dynamics and subgrid physics that adopts different methods for fluid dynamics and tracer transport is used in E3SMv0-HiLAT to minimize water conservation error in the default model. The atmospheric component is tuned to run at an effective horizontal resolution of 1 • and 30 vertical layers (up to 3.6 hPa). In addition to the standard treatments of transport and physical/chemical processes (i.e., tracer transport, turbulence, shallow/deep convection, radiation, ozone/oxidants, aerosols, cloud macrophysics, and microphysics), as used in the LENS simulations (Kay et al., 2015) and described by Neale et al. (2010), we have included some improved treatments of aerosols and clouds. More detailed descriptions of new treatments are given below.

Aerosol Emissions, Aging, Transport, and Removal
The standard CAM5.3 treats major aerosol types (i.e., mineral dust, sea salt, sulfate, black carbon [BC], primary organic matter, and secondary organic aerosol) using a three-mode modal aerosol module, MAM3. The size distributions of particles in MAM3 are described by three sets of lognormal modes, including Aitken, accumulation, and coarse mode (Liu et al., 2012). Different types of aerosols within each of the three modes are internally mixed, while particles from the different modes are externally mixed. In E3SMv0-HiLAT, we include an additional mode (i.e., primary carbon mode) to represent freshly emitted BC and primary organic matter particles, following Liu et al. (2016). Particles in the primary carbon mode cannot have microphysical effects on clouds by acting as CCN until they age into the accumulation mode by being coated with sulfate or secondary organic aerosol. Under present-day conditions, this new treatment was shown to improve the representation of the aging process of BC particles and long-range transport of BC into the Arctic (Wang et al., 2013). Inclusion of this additional aerosol mode is therefore a significant enhancement to the model's capability for subsequent study of 20th and 21st Century climate, when BC transport into the Arctic becomes substantial. Aerosol transport and wet removal will affect snow albedo through the deposition of light-absorbing particles such as BC and dust upon the snow surface. Although this effect is likely to be minimal in our preindustrial simulations, the impact on Arctic snow could be discernible under present-day conditions.
Along with the standard treatments of aerosol-cloud interactions, E3SMv0-HiLAT also includes several modifications, following Wang et al. (2013), that improve the simulation of aerosol wet removal and convective transport. Although aerosols have no direct microphysical impact on convective clouds in the model, convective precipitation can scavenge and remove aerosols. In the default model, a constant fraction of aerosols entering convective cloud base is removed, while the remaining is lifted up. In the new treatment, aerosol activation in convective cloud is calculated at cloud base using the scheme of Abdul-Razzak and Ghan (2000), as for activation in stratiform cloud. In addition, the secondary activation of aerosols (e.g., laterally entrained particles into convective plumes) above convective cloud base is also explicitly treated in convective updrafts. As a result, the unified treatment of convective transport and the removal of aerosols by convective precipitation becomes more physically based. These new treatments also improve the representation of aerosol vertical distributions and long-range transport (Wang et al., 2013).
This model also includes improved treatments of the emissions of dimethyl sulfide (DMS), an important aerosol precursor gas that is emitted to the atmosphere almost entirely from marine ecosystems. After emission to the atmosphere, DMS is chemically oxidized to form SO 2 and other products, with a reactive lifetime on the order of one to several days. Over the Southern Ocean and other remote marine regions, DMS-derived SO 2 is the most important source of CCN in the atmosphere (Carslaw et al., 2013;McCoy et al., 2015;Yang et al., 2017;Yang et al., 2018). These remote marine regions are particularly sensitive to changes in CCN number (Carslaw et al., 2013;Karydis et al., 2012;Moore, Karydis, et al., 2013); thus, it is necessary to adequately represent DMS emissions and DMS-derived CCN in atmospheric models that simulate high-latitude clouds.
Historically, global Earth System Models have typically used observationally based climatologies of oceanic DMS concentrations (Kettle & Andreae, 2000;Lana et al., 2011), in combination with calculations of the piston velocity to parameterize the rate of air-sea transfer. These oceanic inventories neglected important differences in DMS emission rates by different phytoplankton species. In particular, Phaeocystis, a family of phytoplankton prevalent at high latitudes, produces DMS at higher rates compared to phytoplankton species dominant in the tropics and midlatitudes. In previous work by a member of this team and coworkers, an explicit representation of DMS emissions by Phaeocystis was developed and shown to improve agreement of modeled DMS with observations (Wang et al., 2015).
The E3SMv0-HiLAT model configuration couples this improved representation of oceanic DMS concentrations (Wang et al., 2015) to DMS fluxes used by the atmosphere model. Fluxes are calculated in the ocean model using atmospheric surface winds and transferred to the atmosphere model at the ocean-atmosphere daily coupling frequency.
This coupling enables high-latitude aerosols and clouds to respond to simulated changes in oceanic DMS production associated with sea ice retreat and to shift in the simulated distribution of Phaeocystis and other phytoplankton species associated with warming and freshening of the surface ocean, increasing ocean stability, and injection of nutrients into the ocean surface from melting sea ice. Shifts in community composition under future climate change have the potential to change the regional distribution of DMS production, with possible regional-scale impacts on atmospheric surface temperature and sea ice .

Modifications to Cloud Physics
After making all of the adjustments to our model, the default parameterization values for the stratiform cloud macrophysics scheme (Park et al., 2014) would result in overestimation of liquid cloud fraction. To counteract this effect, the threshold value for the minimum grid-mean relative humidity for low liquid stratus was adjusted, as discussed in section 2.5. Also, the "freeze-dry" scheme (Vavrus & Waliser, 2008) is enabled to reduce the fractional areal extent of stratiform liquid-containing clouds in cold and dry environments (e.g., high latitudes). Liquid cloud fraction is scaled by a factor of q v q −1 v 0 when the ambient specific humidity (q v ) is smaller than the threshold (q v0 = 3 g/kg). The freeze-dry scheme is expected to also reduce the low bias in Arctic aerosol loading and liquid water content (Wang et al., 2013).
In many cloud schemes, when condensates are detrained from convective cloud in a supercooled environment, they need to be partitioned into ice and liquid before allocating condensate mass to preexisting liquid or ice stratus. In the shallow convective parameterization of the standard model, the ratio of ice to total condensate detrained from convective cloud is 1 when the ambient temperature is below −35 • C (238.15 K) and 0 for temperature above −5 • C (268.15 K), and the ratio is a linear function of temperature in between the two threshold values. Kay et al. (2016) found that the bias toward a low value of SWCF over the Southern Ocean is associated with insufficient supercooled liquid in CAM5. Considering that shallow convection is the primary source of liquid condensate over the midlatitude oceans, they changed the threshold of 268.15 to 253.15 K and found a large reduction in midlatitude SWCF bias, especially over the Southern Ocean. We adopt their strategy, but to a more moderate degree, by adjusting the threshold from 268.15 to 261.15 K in E3SMv0-HiLAT. Our use of a change in threshold of about half that of Kay et al. (2016) results in a consider-  (Hunke & Lipscomb, 2008) and E3SMv0-HiLAT's Version 5    Table 2 for details of their adjustment during the course of integration.
able bias reduction in SWCF, as discussed in section 2.5 below, while avoiding the more substantial retuning of other parameters that would be required to reestablish radiative balance if the full change in threshold were used.

Sea Ice
The E3SMv0-HiLAT model uses an early implementation of CICE version 6 , largely equivalent to CICE version 5 , but with the single-column physics and biogeochemical aspects of CICE, such as ice ridging, thermodynamics, and hydrology, separated into a submodule ("Icepack"; Roberts et al., 2018). We choose here to treat sea ice as an ice-brine "mushy layer" (Turner et al., 2013;, with prognostic temperature and salinity in seven vertical layers within each of the five thickness categories, which facilitates the inclusion of sea ice biogeochemistry. A delta-Eddington radiation scheme determines albedo and radiative fluxes throughout the snow-ice column for snow, ice, and ponded surfaces (Briegleb & Light, 2007;Hunke et al., 2013). Ice dynamics include an elastic-viscous-plastic rheology (Hunke, 2001;Hunke & Dukowicz, 1997), incremental-remapping advection (Lipscomb & Hunke, 2004), and ridging (Lipscomb et al., 2007). In contrast, CICE version 4  was used for the LENS simulations, including the thermodynamics of Bitz and Lipscomb (1999) in four vertical layers and five thickness categories. In the LENS simulations, aerosols were treated using the original (Community Climate System Model; CCSM) scheme for vertical aerosol cycling (Holland et al., 2012), including their influence on the shortwave absorption and transmission through the ice, while E3SMv0-HiLAT uses a full vertical description of biogeochemical cycling within the ice but without feedbacks with the shortwave parameterization (which is also delta-Eddington). While the lack of aerosol feedback on shortwave absorption represents a reduction in capability, we will first assess the impact of the brine-dependent transport scheme, used for all sea ice biochemical tracers including aerosols, before allowing aerosols to determine albedo in modern-day application of this model.
Melt ponds are included through the "level-ice" scheme , for which ponds pool on undeformed ice, while the LENS simulations employ the more empirical CCSM melt pond scheme (Holland et al., 2012).
Sea ice model parameter differences between E3SMv0-HiLAT and LENS are listed in Table 1. The first three parameters in this table were chosen to calibrate CICE, as described in the next section.

Sea Ice Parameter Refinement From Statistical Estimation
The differences in parameterization of sea ice physics outlined above, relative to LENS, require a thorough consideration of parameter settings. We chose to use the results of a statistically based parameter refinement exercise, shown to produce realistic estimates of sea ice in stand-alone simulations using prescribed atmospheric forcing, as a starting point. Urrego-Blanco et al. (2017) produced an ensemble of 397 stand-alone model simulations spanning different combinations of CICE model parameters, resulting in a variety of responses of sea ice volume and extent. They ranked the configurations using a variance-based distance metric, to capture the spatial character of model skill as compared to different observational data sets of sea ice concentration and thickness. Using this distance metric and building on the work of Urrego-Blanco et al. (2016), who identified parameters for which the simulated sea ice response is more sensitive, the baseline parameter values for the sea ice model component in the coupled E3SMv0-HiLAT system were taken to be those that had allowed the simulation to best match observations of sea ice concentration (EUMET-SAT, NASA Team 1 and 2, and Bootstrap data sets) and sea ice thickness (Unified Sea Ice Thickness Data Record), as documented in Urrego-Blanco et al. (2017).
Important parameters include the snow thermal conductivity (ksno), the maximum snow grain size during melting (rsnw_mlt), and a coefficient R_snw that scales the nominal standard deviation of snow grain size (taken to be 250 μm, based on observations from the Surface Heat Budget of the Arctic Ocean (SHEBA) project, as discussed in Briegleb & Light, 2007). The first parameter is critical for the amount of sea ice produced by congelation at the bottom of the sea ice pack in winter, while the latter two determine snow grain radius, which is considered to grow linearly, once melting commences, up to a maximum value, and which affects surface albedo through the delta-Eddington radiation scheme. Rsnw_mlt effectively puts a lower bound on the albedo by limiting the maximum snow grain size, and R_snw reduces the snow grain size by R_snw × if R_snw > 0 (as it is here) or would increase it by the same amount if R_snw were set to be < 0, where is the standard deviation of the snow grain size distribution.
After the first stage of parameter optimization was conducted in the stand-alone model, parameters with large first-order effects were further adjusted in the coupled E3SMv0-HiLAT model to better match climatological observations (Urrego-Blanco et al., 2016). The CICE model parameters involved in the model tuning are presented in Table 1, along with the parameter values used in the coupled model.

Sea Ice Biogeochemistry
In order to support a more complete representation of marine fluxes of DMS to the atmosphere, we treat sea ice here as a biogeochemically active medium. In doing so, one must acknowledge its multiphase nature-ice crystals are interspersed with brine pockets and channels, creating a protective habitat for algae and, also, because of its location at the air-sea interface, allowing accessibility to both ocean nutrients and light. Modeling this environment comes with challenges. Ice algae combat extremes in salinity and low temperatures with disproportionately high production of the cryoprotectant (Tison et al., 2010) and osmotic regulator (Kirst et al., 1991) dimethylsulfoniopropionate (DMSP), the precursor of DMS. In order to include the impact of ice algae on the marine production of this climatically important trace gas, we have extended CICE5 with a vertically resolved, sea ice biogeochemical submodule (zBGC).
The goal of zBGC is to simulate seasonal patterns of ice algal chlorophyll concentrations and carbon and DMS production in a global and climate evolving context. That is, the component tracer equations satisfy a single set of biochemical and physical constraints, which are in theory applicable to both polar regions as they undergo changing conditions through the current century and beyond. The main assumptions are as follows.
1. Ice algal growth rate is primarily determined by local brine concentrations of macronutrients (nitrate, ammonium, and silicate), dissolved iron and photosynthetically active radiation based on the law of limiting factors (Liebig's Law of the Minimum; Sprengel, 1828). In general, ice algae are not observed to be iron limited, but we include dissolved iron for its potentially important role in the ocean biogeochemistry of iron replete polar regions. In addition, the maximal growth rate satisfies an Arrhenius's (1889) equation dependence on temperature and salinity (also available in Arrhenius, 1967). 2. Biogeochemistry occurs in the liquid volume fraction of the sea ice with no a priori assumptions restricting the vertical distribution to, for example, the skeletal layer in the Arctic or the freeboard in the Southern Ocean. 3. Biogeochemical tracers are present in the ice in two possible phases: mobile and attached. Mobile phase tracers move freely with the brine whose motion is governed by external pressure gradients and desalination processes, which are functions of the temperature and bulk salinity profiles of the ice. In the former case, vertical transport of mobile tracers is expressed as advection by Darcy flow, while in the latter, mixing generated by gravity drainage is parameterized using the mixing length diffusivity of Jeffery et al. (2011). Attached tracers, on the other hand, are fixed to the ice crystals and immune to brine motion. Exchange between phases is governed by prescribed adsorption/desorption time scales and the assumption that ice growth enables adsorption to ice crystals while ice melt facilitates desorption. Both phases are assumed to participate in biogeochemical processes equivalently.
A complete description of the zBGC submodule is given in Hunke et al. (2018). We note that a dedicated paper describing the sea ice model's biogeochemisty and mushy layer physics, in coupled application, has yet to be written.

Ocean
The ocean model is identical to that of the LENS simulations of Kay et al. (2015), as detailed by Danabasoglu et al. (2012), except for its biogeochemistry (as in their simulations, even when biogeochemistry (BGC) is included in the model, it is not allowed to feed back on the ocean physics). The model has 60 vertical levels, with resolution varying between 10 m at the surface and 250 m near the bottom. Notable parameterizations include versions of the eddy parameterization scheme of Gent and McWilliams (1990) and the K-profile boundary layer scheme of Large et al. (1994), as well as the dense water overflow scheme of Briegleb et al. (2010) applied within the Denmark Strait, Faroe Bank Channel, and Ross and Weddell Seas.

Extended Ocean Biogeochemistry
Our focus on high-latitude climate and the impact of marine DMS on cloud forcing require extension to the ocean model biogeochemistry, composed here of a modified version of Biogeochemical Elemental Cycling (BEC; Moore et al., 2001;Moore et al., 2004) and tracegas modules (Elliott, 2009;Wang et al., 2015). The standard BEC includes three explicit phytoplankton functional groups (diatoms, diazotrophs, and smaller phytoplankton) and an implicit group of coccolithophores, along with one general zooplankton functional group. Growth and distribution of phytoplankton and zooplankton are determined by modeled growth-limiting nutrients (nitrate, ammonium, phosphate, iron, and dissolved silicon), photosynthetically active radiation, and local temperatures. Our modified version of BEC adds two phaeocystis groups, one for each hemisphere, implemented based on their nutrient and temperature preferences (Wang et al., 2015). Calculations of the production and removal of DMS and DMSP are driven by a complex set of biogeochemical fields and ecological processes simulated by BEC. Parameters in the extended ocean biogeochemistry modules are set according to the available data from laboratory measurements or field studies (Wang et al., 2015).
The major global ocean ecological fields including nutrients, light, DMS, and biomass have previously been shown to be well reproduced (Moore et al., 2001(Moore et al., , 2004Moore & Braucher, 2008; J. K. Moore, Lindsay, et al., 2013;Wang et al., 2015). The oceanic carbon cycle has also been thoroughly validated in previous studies . Simulated macronutrient distributions are in good agreement with the World Ocean Atlas data sets, with correlation coefficient r > 0.8 and comparable variability between simulated and observed values (J. K. Moore, Lindsay, et al., 2013). With the extended version of the BEC model, individual phytoplankton biomass have been validated against the MAREDAT database and showed reasonable agreement (Wang et al., 2015). In general, the simulated phytoplankton biomass were found to have lower variability than observational data, likely due to the coarseness of model resolution (the same as used in this model), which cannot resolve local spikes of nutrient concentrations. The simulated annual mean DMS concentration of the surface ocean was, in these previous studies, found to be 2.3 nM under the contemporary climate scenario (Wang et al., 2015;, in close agreement with the observation-based estimate of 2.34 nM (Lana et al., 2011). Simulated zonal mean DMS was also found to be in good agreement with the observational estimate.

Initial Conditions
The ocean, which presents the longest time scales of adjustment, was initialized from an observational climatology, while the other components were initialized from previous fully coupled CESM1 simulations: The atmosphere and land started with preindustrial (1850) climate conditions of 1 January; having no prior preindustrial simulation of CICE5, the sea ice was initialized from 1 January of a fully coupled modern-day simulation (at the Year 2006). In any case, after several decades of simulation, the state of the sea ice will have lost most memory of its initial condition. For the ocean, January mean temperature and salinity was used from the Polar Science Center Hydrographic Climatology (PHC2), as was also done with LENS. PHC2 combines the global climatology of Levitus et al. (1998) with Arctic climatology of Steele et al. (2001). The cooling seen in the first 15 or so years in Figure 1 is an indication of a slight mismatch between the clima-  Table 2). Note that the first two of the five periods over which the global imbalance has been quantified include additional years beyond those that comprise the final preindustrial control (piControl) simulation. SW = shortwave; DMS = dimethyl sulfide.
tology of the initial conditions and the internal model state. When biogeochemistry was added to the ocean model after 161 years (see Table 2), the biogeochemical tracers were initialized with fields from a previous fully coupled preindustrial simulation, taken after 75 years of integration; the initial fields for that simulation with extended ocean biogeochemistry had been taken, in turn, from a simulation with standard ocean biogeochemistry after 580 years of integration (as documented by Wang et al., 2015).
The initial condition for the sea ice physical state was obtained from a fully coupled simulation without active ice biogeochemistry. We tested several initialization options for the ice biogeochemical tracer profiles Note. Changes to BGC parameters only affect the physics after the atmospheric model begins using marine model-produced DMS, rather than an input climatology, as an aerosol precursor in cloud nucleation. DMS = dimethyl sulfide; LENS = Large Ensemble Community Project; RH = relative humidity.
for grid cells with nonzero initial condition ice volumes: (1) zero initial values throughout the ice, (2) fixed ocean tracer concentrations throughout the ice, and (3) scaling of the tracers based on their surface ocean value, scaled by the sea ice salinity profile. In general, the different procedures yielded very small differences after about 5 years, though the zero initial condition produced the greatest initial adjustment. Given that the concentration of macronutrients without active biochemistry is expected to scale with salinity, we chose the third option as representing the most physical upper bound profile.
Before proceeding, we note that land use change may have been appreciable before 1850, and coal was already being burned to power steam engines. A more precise definition, which also takes into considera-   tion the variability in natural forcings, has been suggested as the period 1720-1800 (Hawkins et al., 2017). The state supporting the near perfect balance that would be reached after a millenium or more of model integration would depend on the precise definition. Here, however, we use an approximate definition of preindustrial as before 1850 and consider a simulation of several centuries to be adequate for our purpose.

Tuning and Changes in Configuration
During the course of the preindustrial control simulation (abbreviated piControl), a number of relatively small changes were made that affected the physical solution, as noted in Table 2 and sketched in Figure 2. This followed a period of experimentation in which larger changes in configuration were considered. Along with the table, any figures that show full time series have these transition points marked.
The first changes to configuration, after 35 years, had the greatest impact on the model's climate, allowing for more supercooled liquid water to be detrained from convective cloud in the atmospheric column. Following Kay et al. (2016), as discussed above, we allowed for supercooled water down to a temperature of 261.15 K (only about half as far from the freezing point as the value used by Kay et al., 2016). Simultaneously, the minimum relative humidity for low clouds (cldfrc_rhminl) was increased from 0.8975 to 0.91 to compensate for the increase in SWCF. As can be seen in Figure 1, these changes do compensate to approximately maintain overall radiative balance while substantially reducing the bias in SWCF in the tropics and the Northern Hemisphere (Figure 3). At this same point, at the beginning of Year 36, an additional modification was made in one sea ice model parameter, reducing the thermal conductivity of snow from the value that had been used in stand-alone simulation to one that may be more appropriate to coupled simulation (parameter ksno reduced from 0.357 to 0.3). After integrating this configuration for 114 years, a 30-year simulation was branched off Year 150 that involves instantaneously quadrupling the CO 2 concentration in the model from its preindustrial value (hereafter referred to as 4xCO2). The results of the 4xCO2 simulations (see below for a mention of the second 4xCO2 simulation) are discussed briefly in section 3.5 and in much greater detail in section 4.
From the overlapping time series of model integrations included in Figure 4, it would seem that sea ice area in the Southern Hemisphere ceases its expansion after Year 35 (as seen in an unused branch of the integration, extended through 50). But this is also when the correction to SWCF was made to the main branch of the piControl simulation, and as one might expect, the reduction in short wave radiation reaching the marine surface over the Southern Ocean allows for continued expansion in ice area in the Southern Hemisphere.
The sea ice area and volume in the Southern Ocean are clearly becoming excessive, even for preindustrial conditions. While the correction in SWCF would be adding to this bias, it is also sufficiently important  Table 2; overlapping periods of integration are included to provide an indication of the impact of sea ice model parameter changes, first after 35 years of simulation (change of one parameter to a value believed to be more appropriate for fully coupled use), then after the 110th year (in order to reduce sea ice area and volume; see Table 2 and discussion thereof). NH = Northern Hemisphere; SH = Southern Hemisphere.

10.1029/2018MS001524
that we were not willing to forsake this improvement. Two further changes, however, were made to sea ice model parameters after 110 years, aimed at reducing sea ice area in the Southern Hemisphere. This involved further reduction of the thermal conductivity of snow (parameter ksno reduced from 0.3 to 0.27) and an increase in the maximum grain radius of melting snow (parameter rsnw_mlt increased from 800 to 1,500; both of these changes were made globally, affecting sea ice in both hemispheres). As discussed in section 2.2.1, the decision to adjust these two parameters, rather than others that could also affect a reduction in sea ice volume, was based on their relatively high degree of sensitivity as determined by Urrego-Blanco et al. (2016).
Marine biogeochemistry was added after 161 years of simulation, with changes made to the ocean model's BGC after another 10 years of simulation. At this point, however, the marine BGC was still passive.
Marine BGC was only allowed to affect the physics when, after 204 years of model simulation, the atmospheric model began using the marine model-produced DMS field, rather than a prescribed climatology, as an aerosol precursor in cloud formation. Finally, after 228 years of simulation, changes were made in the sea ice model's BGC parameters based on more recent experience gained with our still-new sea ice BGC system in stand-alone simulations. This last change does have some small effect on the physical climate through modification of DMS fluxes. After making this final change, the model was run for another 26 years, out to a total of 253 years of preindustrial simulation. A second 4xCO2 simulation was branched at Year 236 and run for 50 years; this simulation is discussed briefly in section 3.5 and in much greater detail in section 4.

Assessment of Resulting State
In assessing the resulting state, comparisons are made here to observations and to the extraordinarily well balanced LENS control simulation. While remembering that our aim is to establish a greater capability for the study of high-latitude climate processes, we do so while attempting to minimize loss of fidelity elsewhere.

Remaining Trends in Energy Balance
Returning to the global radiation budget of Figure 1, there is some evidence of multiple time scales of adjustment. To see this, one may consider the second of the five periods over which the global radiation imbalance has been quantified, from Years 36 through 150 (the preindustrial Control simulation, piControl, includes this period through Year 110, with the years from 111 to 150 unused except within this plot). The imbalance is at a minimum during this period, at a level well below 0.1 W/m 2 , before rising back up to a level of around 0.2 W/m 2 over the third and fourth periods of the simulation, from Years 111 through 227. Some explanation may be derived from the other heat-related time series plot of Figure 1, showing the evolution of heat content in different depth ranges of the ocean. Up until around Year 110, a secular warming of the deep ocean is offset by cooling of the upper 700 m. The changes made at this point to sea ice model parameters, intended to reduce excess Southern Hemisphere sea ice extent, appear to have been partially responsible for the reduction of upper ocean cooling. But there is the suggestion of some reduction in the trend of upper ocean cooling appearing around that time even without the parameter change, as seen in the overlapping (and otherwise unused) years of integration from Years 111 through 150. Together, the change in sea ice parameters after Year 110, along with the difference in adjustment time scales for deep and upper ocean, explains the emergence of the higher ∼0.2 W/m 2 level of imbalance recorded during these middle years, with the more persistent deep ocean trend only allowing the overall radiative balance to drop back down to the ∼0.1 W/m 2 level over the last few decades of the piControl simulation.
There is a suggestion of a resumption in upper ocean cooling over the last 50 years, with the change from climatological to marine model-based DMS fields, for use in the nucleation of stratiform clouds. It is also possible, however, that this resumption of cooling is due to internal variability. In either case, an analysis of the time series of surface air temperature has been done, based on fits to a semi-infinite diffusion model (see Appendix A; MacMynowski et al., 2011); that analysis indicates that the changes made to model physics after 36, 110 and 204 years did not alter global mean surface air temperature in a way that appreciably affected the results.

Atmosphere
Our E3SMv0-HiLAT simulation is generally colder than LENS, as seen in Figure 5. The comparison indicates warmer temperatures over some land masse, and also broadly over the Arctic during the JJA months. During the DJF months, however, there is a reversal of this difference over much of the Eurasian sector of the Arctic, where colder temperatures are evident in our E3SMv0-HiLAT simulation. There is also a seasonal reversal of temperature difference over North America. Another region in which our E3SMv0-HiLAT simulation tends to be warmer, contrary to its cooler global tendency, is over the tropics.
Our E3SMv0-HiLAT simulation also tends to be drier than LENS, if far from uniformly so, as seen in Figure 6. This overall tendency toward less precipitation is to be expected from the cooler marine surface air temperatures. Differences in precipitation show up more strongly around the intertropical convergence zone (ITCZ), where the strongest magnitudes of precipation are found, but the patterns in the two models are broadly similar. The double ITCZ is somewhat more pronounced in E3SMv0-HiLAT.

Ocean
A time series of annual mean, global mean sea surface temperature (SST) drifts over the entire 253-year period, for both E3SMv0-HiLAT and LENS, is included within Figure 7a. While SST cools over the first 35 years in both models, the cooling is greater in E3SMv0-HiLAT by around a factor of 2, as also anticipated by the surface air temperature comparison of Figure 5. After a short warming adjustment to the changes of radiation and cloud parameters made at Year 36, the E3SMv0-HiLAT SST continues to cool with time toward the end of the simulation. In contrast, the LENS SST undergoes a cooling phase for roughly 130 years, after which it displays a gradual rise for the remainder of the simulation. This leads to a −1.0 • C drift in E3SMv0-HiLAT relative to the initial state, while the drift is −0.3 • C in LENS. We further compare the model results with the Hadley Center's SST estimation of the preindustrial period 1870-1900 (Rayner et al., 2003), whose global mean difference from the modern value estimated in PHC2 is −0.25 • C (Figure 8).
In this sense, the simulated cooling is to be expected, since the initial condition is based on PHC2, but the cooling in E3SMv0-HiLAT exceeds that of LENS.

Journal of Advances in Modeling Earth Systems
10.1029/2018MS001524 Even before the initial change in model configuration is made after 35 years, it is notable in Figure 1 that global ocean heat contents drops and then recovers over that initial period, while upper ocean heat content drops steadily. With the initial drop in total ocean heat content being about 10 23 J over 15 years, this corresponds to an average air/sea heat flux anomaly of around 1 2 W/m 2 . The recovery in global heat content over the next 15 or 20 years corresponds to a heat flux anomaly of similar magnitude but opposite sign. From an inspection of the various terms contributing to globally integrated surface heat over the initial 35 years (not shown here), the latent heat flux was noted to display a downward trend amounting to around 2 W/m 2 over 30 years. A latent heat flux that is anomalously high by 1 2 W/m 2 over the first 15 years would potentially explain the initial loss of global ocean heat content. And a latent heat flux that is anomalously low by a similar amount over the subsequent 15 to 20 years would be sufficient to explain the recovery in global ocean heat content. It is plausible then to suggest that our model does not find its radiatively balanced state until the SST has developed the cold bias seen in Figure 7. We also note that much of the early tendency toward greater cooling of SSTs in E3SMv0-HiLAT may be attributed to the North Pacific Subtropical Gyre (not shown). If the suggestion above is correct, then the North Pacific Subtropical Gyre would be a key region in which our model, as configured, requires the development of a cold bias, allowing a drop in latent heat flux, before it will establish a radiative balance.
Looking again at the later more equilibrated solutions, SST bias patterns for both models relative to the Hadley Center's SST preindustrial estimation are shown in Figure 8. E3SMv0-HiLAT has an overall cold bias globally, with an enhanced cold bias over much of the subtropical and subpolar gyres. In comparison, LENS has a similar cold bias pattern over tropical regions but a warm bias over subpolar regions. The excessive cold bias over midlatitudes in E3SMv0-HiLAT is mostly due to reduced surface shortwave heat flux over these regions (not shown here), partially compensated by reduced latent heat loss in response to the surface cooling; some of the tendency toward reduced shortwave heat flux can be ascribed to our use of the Kay et al. (2016) allowance for detrainment of supercooled liquid water, as described above in section 2.1.2. The cooling bias over the Southern Ocean is clearly associated with the expansion of sea ice extents in both the winter and summer seasons, as will be shown below; this cold bias would also be seen in LENS, though to a smaller degree.
As discussed above, the time series of global ocean heat content shows nearly zero change for the first 110 years (Figure 1b). In Figure 7c, one sees cooling over the upper 500 m that slightly overcomes the warming in the middepth and deep ocean. After 110 years, the middle-to-deep ocean continues to warm and becomes dominant in the trend of ocean heat content, indicating net ocean heat uptake, and accordingly, top-of-model  heat gain in the coupled system. The LENS simulation shows greater ocean heat uptake, since the upper ocean cooling is weaker (Figure 7e). Both simulations show a net ocean heat gain.
The cool temperatures over the uppermost 500 m of the ocean appear to be a consequence of the E3SMv0-HiLAT model's tendency to produce low SSTs. An examination of zonally averaged temperature biases within basins (not shown) finds the greatest cool biases of the upper ocean within the latitudinal range of 20 • S to 50 • S. Some of the coldest biases in SST occur over this same range of latitude (Figure 8), with subduction of these cold surface waters producing the strong upper layer temperature drift.
While the LENS simulation undergoes less vertical redistribution of heat and salinity, the bias of cooling above the thermocline and warming below is common in Earth system models (Kuhlbrodt & Gregory, 2012;Kuhlbrodt et al., 2018). The small overall oceanic heat trend is less significant and could readily be reversed-at least through the duration of the few centuries that we have run-by a slight retuning of parameters; one should only expect the overall radiative imbalance to drop below 0.1 W/m 2 in a millenial scale simulation.
The sea surface salinity (SSS) drifts are similar in the two models, showing a freshening trend of similar magnitude over the 253 years in both models (Figure 7b). Given the very different feedback mechanisms impacting SST versus SSS, it is not surprising to see different responses in these two fields. By the end of the simulations, the drifts from the initial state are −0.30 and −0.34 psu in E3SMv0-HiLAT and LENS, respectively. The slightly greater freshening drift in LENS is mostly caused by stronger freshening over the tropical Indian Ocean. The similar evolution of SSSs are indicative of similar hydrological boundary forcing and transport through surface ocean currents in these two models. It further confirms that the difference in evolution of SSTs is not caused by the surface ocean circulation but by the boundary radiative forcing.
The global volume mean salinity falls only very slightly from the PHC2 value in both models, even while salinity is redistributed within the ocean, with the upper 700 m getting fresher and the ocean below 700 m getting saltier (Figures 7d and 7f). These global mean and depth-dependent trends are consistent with previous CCSM3 and CCSM4 results (Danabasoglu et al., 2012).
The change in SWCF does lead to some change in ocean mixed layer depths; whereas the maximum in the zonal average of mixed layer depth over the Southern Ocean was around 90 m in Figure 9, after the improvement was made to SWCF the Southern Ocean maximum of zonal mean mixed layer depth deepens to around 97 m, which would be an improvement from the shallow bias to Southern Ocean mixed layer depths reported for this ocean model (Danabasoglu et al., 2012).
Global and Atlantic overturning circulations in the E3SMv0-HiLAT model are very similar to those of LENS, despite the difference in surface conditions between the two models. Referring back to Figure 5, the warmer Figure 10. Seasonal concentrations of sea ice, averaged over Years 234-253 of the piControl run, and compared with climatological seasonal cycle derived from modern-day satellite observations (Cavalieri et al., 1996). DJF temperatures off southern Greenland, in E3SMv0-HiLAT relative to LENS, may be more a function of warmer winter time conditions over the land mass of Greenland itself than of any systematic difference in the overturning circulations. The circulations of Antarctic Bottom Waters are weak in both versions of the model, as is typical of non-eddy resolving ocean climate models (Danabasoglu et al., 2014). Drake Passage mass transport (not shown) is slightly weaker in our model, at around 142 Sv, as compared with 155 in LENS.

Sea Ice Extent and Volume
Even with the modification to sea ice parameters after 110 years, the annual maximum area of sea ice in the Southern Hemisphere remains around 50% higher than seen in present-day observations, as seen in Figure 4. This is particularly evident in the summer sea ice area. More impact from the adjustment to parameters is apparent on annual minimum areal extent than on the maximum. There is also a significant reduction in Arctic volumes, which was not the intention of the parameter change but which may result in reduced bias. Indeed, sea ice area and volume in the Arctic are similar to LENS's results, in both winter and summer seasons.
Spatial comparisons of sea ice concentration and volume distributions, relative to modern-day observed fields, are seen in Figures 10 and 11. Southern ice concentrations are more extensive in most regions, while Austral Spring thicknesses are particularly great around the prime meridian and off of the Ross and Bellingshausen Seas. Most likely, this increase in sea ice coverage represents an increased bias, one indicative of the correction of one contribution to compensating biases, with the improved SWCF (Figure 3). In any case, Figure 11. Seasonal thickness of sea ice, averaged over Years 234-253 of the piControl run, and compared with climatological seasonal cycle derived from modern-day satellite observations (Kurtz & Markus, 2012;Yi & Zwally, 2014). Ice thickness is only shown where ice concentration exceeds 15%. even while our initial set of experiments will be done from the piControl state, it will be informative to learn how the sea ice responds through transient twentieth Century simulation.
We also note that one of the next planned additions to E3SMv0-HiLAT is an iceberg model. Although not primarily intended for this purpose, one potential benefit to future inclusion of the iceberg model may be a reduction in overall extent and volume of sea ice over the Southern Ocean, as found recently by Storkey et al. (2018).

DMS and Clouds
The performance of DMS-related marine biogeochemistry in the model is summarized in Figures 12 and  13. Near-surface concentrations of DMS from the piControl and 4xCO2 simulations, considered in their annual and zonal means (Figure 12a), are consistent with present-day observations from Lana et al. (2011). DMS concentrations are sparsely sampled in natural sea ice: Modeled concentrations are consistent (within error bars) with available in situ observations in the Southern Ocean (Trevena & Jones, 2006;Tison et al., 2010) but higher in the Arctic (Levasseur, 2013). Overall, Chlorophyll-a (Figure 12b) also compares favorably with observations, but surface concentrations are significantly underestimated in the subpolar Northern Hemisphere. Figure 13 shows that the simulated zonal mean DMS fluxes for both the piControl and 4xCO2 simulations are more consistent with the updated climatology of Lana et al. (2011) than with the more limited DMS flux climatology of Kettle et al. (1999) and Kettle and Andreae (2000), which is the default forcing data set for CESM1 and the standard for the Coupled Model Intercomparison Project Phase 5 (CMIP5)  (Gradinger et al., 2012), Barrow (Deal et al., 2011), and Baffin Bay (Michel et al., 2002), with those sea ice observations weighted by the annual mean ice concentration for the given latitude band from passive microwave data (Cavalieri et al., 1996). DMS = dimethyl sulfide; Chl a = chlorophyll-a.
Atmospheric Chemistry and Climate Model Intercomparison Project. It is encouraging to see that our model produces DMS fluxes that are closer to the new observational estimate than to the old while establishing the new capability for this aspect of cloud forcing to respond to the prognostic marine field.
The inclusion of this capability reveals a significant sensitivity of the fluxes to the climate state. The most prominent differences between the DMS fluxes of 4xCO2 and piControl simulations are an overall increase in emissions south of 30 • S and in the high Arctic and a decrease in fluxes between 30 • S and 40 • N but with large geographic variations ( Figure 14). The increased fluxes in the Southern Ocean, the sea ice zone south of 65 • S, and the Arctic are consistent with increases in surface DMS concentration from oceanic phytoplankton activity, while the decreases elsewhere are consistent with decreased surface DMS concentrations . The changes in DMS concentrations occur despite relatively modest changes in phytoplankton biomass, as DMS surface concentrations reflect the phytoplankton community structure and not strictly the chlorophyll biomass (Wang et al., 2015;Wang, Maltrud, Burrows, et al., 2018). The sea ice algal contribution to DMS production decreases by almost a factor of 2 between the piControl and warmer 4xCO2 climates.
The change in sulfate aerosol burden with use of prognostic marine DMS is qualitatively similar in piControl and 4xCO2 states, as seen in Figure 15 (top middle and top right panels, respectively; note that because the model records sulfate aerosol burden only where the planet is seasonally sunlit, that we have not calculated a statistical significance of the difference in sulfate burdens at polar latitudes, and so there is no stippling there). Changes in SWCF (middle panels) are evident in several regions, with opposite responses seen between piControl and 4xCO2 states over portions of the high latitudes of both hemispheres and also over much of the (split) ITCZ of the Pacific. Stippling indicates much of this difference to be significant, consistent with the finding of  of a statistically significant impact of DMS on cloud forcing under RCP8.5 (Representative Concentration Pathway 8.5, a greenhouse gas concentration trajectory resulting in an approximate 8.5 W/m 2 radiative forcing by the year 2100) conditions, which has a greenhouse gas forcing similar to that in our 4xCO2 scenario. Changes in total cloud fraction (bottom panels) are again opposite in sign between piControl and 4xCO2 states over much of the ITCZ, with stippling there indicating significance. Elsewhere, changes in total cloud fraction appear to be less systematic.
The impact of the new DMS formulation on DMS emissions was discussed in detail by Wang and Moore (2011) and Wang et al. (2015), and the subsequent impact of the new emissions formulation on cloud properties (aerosol indirect effects) was described in greater detail in .

Radiative Forcing and Climate Response
Introduction of a new model begs the question as to the performance of the model in representing response to radiative forcing. There is no benchmark value of climate sensitivity (the equilibrium climate response to a doubling of the atmospheric CO 2 concentration) to which the model can be compared. However, here we evaluate the global mean response of the model to abruptly quadrupling the model's CO 2 concentration from its preindustrial value (4xCO2 simulations). This is a standard simulation under the CMIP5 , allowing us to compare results with those of other models. In particular, a "Gregory plot" (Gregory et al., 2004) is a useful tool for evaluating metrics akin to climate sensitivity. This involves creating a scatter plot of annual mean change (difference from piControl) in surface air temperature (x axis) and net Top-of-Atmosphere (TOA) radiative flux (y axis). Regression through the points yields estimates of adjusted radiative forcing (y intercept) and equilibrium climate sensitivity (half of x intercept, corresponding to a doubling of the CO 2 concentration).
As is indicated in Figure 2, we performed two 4xCO2 simulations, one prior to making parameter value modifications to ksno and rsnw_mlt and one after all modifications were made. These will subsequently be referred to as the Y150 and Y236 versions of the model, respectively, in reference to the preindustrial control year from which the 4xCO2 runs were begun. The Gregory plot for the two 4xCO2 simulations conducted with the HiLAT model is shown in Figure 16. The Y150 and Y236 versions of the model have climate sensitivities of approximately 3.77 and 3.50 • C, respectively, based on moderate length simulations of 30 and 50 years (a more definitive assessment would require 150 years of simulation, as discussed by Gent, 2009, andby Andrews et al., 2012). These values are not statistically significantly different from each other, and the uncertainty estimates for the Y236 version of the model nearly fall entirely within the error bars of the Y150 model version. The mean values are essentially identical to that of a simulation we also performed with CESM 1.0.3 (3.49 • C; not shown), from which E3SMv0-HiLAT is derived. A sample of the CMIP5 models (Andrews et al., 2012) has values  Kettle et al., 1999, and ;Kettle & Andreae, 2000). Both modeled and prescribed DMS fluxes assume the piston velocity parameterization of Nightingale et al. (2000). (right column) Contoured differences of annual average DMS fluxes between the piControl and Prescribed (top right), 4xCO2 simulation and Prescribed (middle right), and between 4xCO2 and piControl simulations (lower right). Stippling in the panel at lower right indicates statistical significance at the 95% level according to a Welch's t test. DMS = dimethyl sulfide. that range between 2.08 and 4.67 • C, with a mean of 3.37 • C, a median of 3.45 • C, and a standard deviation of 0.83 • C, putting E3SMv0-HiLAT only very slightly above this mean value.
Comparing adjusted radiative forcing from the 4xCO2 runs, E3SMv0-HiLAT is at the upper end of the pack at 8.03 and 9.04 W/m 2 , for initial and final evaluations, respectively. In this respect, E3SMv0-HiLAT differs from CESM 1.0.3, which reported a more midrange value of 7.65 W/m 2 . The same CMIP5 model set ranges between 5.17 and 8.62 W/m 2 , with a mean of 6.89 W/m 2 , a median of 6.49 W/m 2 , and a standard deviation of 1.12 W/m 2 . The Y236 version of E3SMv0-HiLAT is nearly two standard deviations above the model mean value and unique in being that far on the high side of the distribution. The high value of adjusted radiative forcing in the Y236 4xCO2 simulation could either be explained by the change in two sea ice parameters (see Table 2) or by the use of marine model-produced DMS in nucleation of stratiform clouds. However, as is illustrated in Appendix A, these changes have minimal effects on global mean temperature change in the piControl run and estimates of equilibrium climate sensitivity. Note that sulfate aerosol burden as logged by the model (panels a, d, and g) is only defined during sunlit months, such that the field at polar latitudes represents a partial annual average. Panels b, e, and h show shortwave cloud forcing, while panels c, f, and i show total cloud fraction. Differences are between cases in which the atmosphere receives model-produced DMS and the corresponding case in which a DMS flux climatology is used; the middle column is the difference under preindustrial conditions, the right-hand column the difference under 4xCO2 conditions. Stippling in these difference plots indicates statistical significance at the 95% level according to a Welchs t test. Because of the partial record of sulfate burden at polar latitudes, we have not calculated a statistical significance of the difference in sulfate burdens there, and so there is no stippling at the high latitudes of (d) and (g). DMS = dimethyl sulfide.
Globally, temperature change in E3SMv0-HiLAT is slightly less sensitive to global radiative forcing than CESM 1.0.3. However, the differences are minor, leading to the conclusion that the modifications made to produce E3SMv0-HiLAT did not substantially alter global climate sensitivity. Nevertheless, the lack of global average change could be due to compensating regional changes. Figure 17 shows net TOA radiative flux and surface air temperature changes in the two 4xCO2 simulations. Responses in both simulations are qualitatively similar, and quantitative differences between the two different responses are small. The largest radiative flux differences are in the Equatorial and Northern Pacific, the Arctic, and the Southern Ocean. These are known regions of persistent cloud cover, so it is difficult to ascertain whether any of these changes were due to alterations to the model physics, differences in model state (as was seen in Figure 7, the Year 150 base period is warmer than the Year 236 period), or noise due to cloud field variability. For surface air temperature, again, the qualitative responses in the two simulations are similar. The largest quantitative differences are increases in temperature in areas with known monsoon precipitation (the Indian subcontinent, East Asia, and the American Southwest), decreases in temperature over nearly all other land regions, increases in temperature in the Southern Ocean (in the same regions of increased radiative flux, suggesting reduced cloud cover), and a decrease in temperature in the North Atlantic (corresponding to a known response of the Atlantic Meridional Overturning Circulation; e.g., Drijfhout et al., 2012). Although the results discussed here are suggestive, because the analysis was performed only over 10-year periods, we cannot rule out the possibility that these results are also influenced by natural variability.

Summary and Outlook
We have presented our modified version of the CESM1 Earth System Model, E3SMv0-HiLAT, as an enhanced model for the study of high-latitude processes and as a context in which to further develop improved representations of those processes. The advances include extended marine biogeochemistry and dependence of stratiform cloud nucleation on marine emissions, along with improved aerosol transport to high latitudes and reduced bias in SWCF over the Southern Ocean.
The inclusion of biogeochemistry in a coupled simulation with CICE, facilitated by an ice-brine mushy layer, is a first. Along with the extended ocean model BGC, the marine components were shown to produce a DMS flux at the air/sea interface that resembles the zonal mean of the observed flux. Use of these prognostic marine model-produced DMS fluxes allows one to consider changes in atmospheric model cloud nucleation, under scenarios such as the 4xCO2 presented here (which is not necessarily beyond the range of radiative forcing scenarios one must consider toward the end of our present century).
A number of biases known from the CESM1 persist in E3SMv0-HiLAT, as was discussed in section 3. While one reduction in bias is especially substantial that of SWCF, this correction has consequences for other fields: SST anomalies, evaluated relative to the Hadley Center's estimate of the preindustrial state (Rayner et al., 2003), show a stronger cold bias than that of LENS. Sea ice extent in the Southern Hemisphere, as shown in Figures 4 and 10, would appear to be excessive, even after partial alleviation through the adjustment of two sea ice model parameters. That adjustment, as well as the initial setting of sea ice model parameters, were guided by parameter estimation (Urrego-Blanco et al., 2016;Urrego-Blanco et al., 2017).
While climate sensitivity is essentially unchanged from that of CESM1, the adjusted radiative forcing of our new model is nearly two standard deviations above that of a sampling of CMIP5 models (Andrews et al., 2012). Being an outlier does not necessarily mean that E3SMv0-HiLAT is worse in this respect than LENS. We intend to perform additional analyses in the near future to understand whether the unusually strong initial radiative response of the final configuration of E3SMv0-HiLAT, as evaluated in the Y236 4xCO2 simulation, is due to the atmosphere's dependence on marine model-produced DMS or to the last change made in sea ice model parameters. This will include a more thorough investigation of how the modifications to the model affected feedback strength, especially cloud feedbacks.
We plan to use E3SMv0-HiLAT to pursue additional analysis and experiments focused on high-latitude science problems, making use of the model's additional capability to resolve changing marine biogeochemical activity and its impact on the atmosphere and clouds. More specifically, we will pursue a set of experiments to evaluate climate response, as mediated by biogeochemistry originating in ocean and sea ice, to changing freshwater and nutrient inputs from Antarctica. While the initial set of experiments will be done from the piControl state established here, we will also perform a historical simulation, with attention to the sea ice distributions over the late 20th and early 21st Centuries. Figure 17. Maps of net TOA radiative flux (left; W/m 2 ) and surface air temperature (right; K) change in the two 4xCO2 simulations from their corresponding piControl simulations. All 4xCO2 values were averaged over Years 21-30, which is the latest 10-year period shared by both simulations. Bottom panel is the difference between the two 4xCO2 simulations to show the effects of model parameter and physics changes on model response.

Appendix A: Semi-Infinite Diffusion Fits
Figure 7 provides a time series of global, annual mean SST evolution over the entire 253-year piControl run, including all of the model changes described in Table 2 and Figure 2. The SST change in Figure 7 shows evidence of response to an abrupt forcing, likely due to adjustment to a mismatch between the internal model state and the forcing files, which were generated with a different model. Figure A1. Global, annual mean surface air temperature change in the piControl run (black) and semi-infinite diffusion fits to that data (colors) using the leave-one-out procedure described in Appendix A. Dashed lines indicate the equilibrium global mean temperature for each of the fits, calculated as the asymptotic value of semi-infinite diffusion (equation (A1)).
Here we consider a time series of global mean surface air temperature and the question of whether changes in model physics have had an appreciable impact on this one measure of the global mean climate. The global mean temperature response to a global-scale step forcing can be described via semi-infinite diffusion (MacMynowski et al., 2011): with parameters and , where erfc is the complementary error function: Figure A2. Parameter values for (top) and (bottom) for the four different semi-infinite diffusion fits (equation (A1)) to the global, annual mean surface air temperature model output for the piControl simulation. Blue dots indicate the best fit values from nonlinear least squares regression, and red dots indicate 95% confidence intervals for those parameters. 10.1029/2018MS001524 Using Matlab's native nonlinear least squares fitting algorithm lsqcurvefit, we performed four semi-infinite diffusion fits to the global, annual mean temperature time series: the entire time series (Years 1-253), omitting Years 37-110, omitting Years 111-204, and omitting Years 205-254. These omitted periods reflect changes in the model's physics, with the aim of understanding when a change in the physics might have led to statistically significantly different semi-infinite diffusion fits. Figure A1 shows the results of the fitting procedure under this "leave-one-out" trial, and Figure A2 shows the semi-infinite diffusion parameter values for and , as well as their 95% confidence intervals from the nonlinear fitting procedure. From these figures, there are several apparent conclusions. As was hypothesized, semi-infinite diffusion fits the model time series rather well. There is no appreciable difference in the different fits. None of the leave-one-out fits robustly indicates that any change in the model physics resulted in a substantial change in the model's global mean climate. This may be a result of the short duration of the different periods involved in the leave-one-out procedure. 89233218CNA000001). This research also used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract DE-AC02-05CH11231. This paper has been greatly improved thanks to the efforts of two anonymous reviewers and the Editors. The ship-based sea ice and snow thickness data were provided by the SCAR Antarctic Sea Ice Processes and Climate (ASPeCt) program (aspect.antarctica.gov.au). Public release of the E3SMv0-HiLAT model source code will be made available, pending review (https:// github.com/lanl/E3SMv0-HiLAT). Climatological data from the final 20 years of the preindustrial control simulation is available online: atmosphere (https://doi.org/10.5281/ zenodo.2554182), ocean (https://doi. org/10.5281/zenodo.2554321), and sea ice (https://doi.org/10.5281/zenodo. 2554240).