Effective radiative forcing of anthropogenic aerosols in E3SM version 1: historical changes, causality, decomposition, and parameterization sensitivities

. The effective radiative forcing of anthropogenic aerosols (ERF aer ) is an important measure of the anthropogenic aerosol effects simulated by a global climate model. Here we analyze ERF aer simulated by the E3SM version 1 (E3SMv1) atmospheric model using both century-long free-running atmosphere–land simulations and short nudged simulations. We relate the simulated ERF aer to characteristics of the aerosol composition and optical properties, and we evaluate the relationships between key aerosol and cloud properties. In terms of historical changes from the year 1870 to 2014, our results show that the global mean anthropogenic aerosol burden and optical depth increase during the simulation period as expected, but the regional averages show large differences in the temporal evolution. The largest regional differences are found in the emission-induced evolution of the burden and optical depth of the sulfate aerosol: a strong decreasing trend is seen in the Northern Hemisphere high-latitude region after around 1970, while a continued increase is simulated in the tropics. The relationships between key aerosol and cloud properties (relative changes between pre-industrial and present-day conditions) also show evident changes after 1970, diverging from the linear relationships exhibited for the period of 1870–1969. In addition to the regional differences in the simulated relationships, a reduced sensitivity in cloud droplet number and other cloud properties to aerosol perturbations is seen when the aerosol perturbation is large. Consequently, the global annual mean ERF aer magnitude does not increase after around 1970.The ERF aer in E3SMv1 is relatively large compared to the recently published multi-model estimates; the primary reason is the large indirect aerosol effect (i.e., through aerosol–cloud interactions). Compared to other models, E3SMv1 features large relative changes in the cloud droplet effective radius in response to aerosol perturbations. Large sensitivity is also seen in the liquid cloud optical depth, which is determined by changes in both the effective radius and liquid water path. Aerosol-induced changes in liquid and ice cloud properties in E3SMv1 are found to have a strong correlation, as the evolution of anthropogenic sulfate aerosols affects both the liquid cloud formation and the homogeneous ice nucleation in cirrus clouds (that causes a large effect on longwave ERF aer ). As suggested by a previous study, the large ERF aer appears to be one of the reasons why the model cannot reproduce the observed global mean temperature evolution in the second half of the 20th century. Sensitivity simulations are performed to understand which parameterization and/or parameter changes have a large impact on the simulated ERF aer . The ERF aer estimates in E3SMv1 for the shortwave and longwave components are sensitive to the parameterization changes in both liquid and ice cloud processes. When the parameterization of ice cloud processes is modiﬁed, the top-of-model forcing changes in the shortwave and longwave components largely offset each other, so the net effect is negligible. This suggests that, to reduce the magnitude of the net ERF aer , it would be more effective to reduce the anthropogenic aerosol effect through liquid or mixed-phase clouds.


Introduction
Variations in atmospheric composition have played an important role in causing the observed (historical) climate changes since the pre-industrial era. The concentrations of both greenhouse gases and aerosols have changed, with partially offsetting impacts on the Earth's radiative balance. Compared to greenhouse gases, the lifetime of aerosols is much shorter and the processes involved in the aerosol life cycle are more complex. In addition, aerosols affect the Earth's radiative balance not only by directly scattering and absorbing shortwave and longwave radiation but also by acting as cloud condensation nuclei and ice nuclei, hence changing the cloud formation process and cloud radiative forcing. Consequently, aerosols have been considered to be a group of the most uncertain forcing agents in the climate system (Penner et al., 2001;Forster et al., 2007;Boucher et al., 2013;Szopa et al., 2022;. Due to our incomplete knowledge of the aerosol-related processes and the limited spatial and temporal resolutions used by Earth system models or global climate models, most such models have rather simplified representations of aerosol lifecycles and their interactions with clouds, precipitation, and radiation. The complex and multi-scale aerosol-related processes are simplified and represented differently in the models, leading to large difference in the simulated aerosol properties  and their impacts on clouds . As a result, the estimated anthropogenic aerosol effects in global climate models are quite uncertain for both the aerosol-radiation interactions Myhre et al., 2013) and aerosol-cloud interactions (Quaas et al., 2009;Ghan et al., 2016).
In an earlier study, Golaz et al. (2019) report that the present-day (average of 1995-2014) net (shortwave + longwave) effective radiative forcing of anthropogenic aerosols (ERF aer ) in the atmosphere component of the Energy Exascale Earth System Model version 1 (E3SMv1) is about −1.65 W m −2 . This value is outside the likely range of −0.63 to −1.37 W m −2 estimated from models participating in the Radiative Forcing Model Intercomparison Project (RFMIP, Smith et al., 2020) for the World Climate Research Programme Coupled Model Intercomparison Project Phase 6 (CMIP6, Eyring et al., 2016). Furthermore, E3SMv1's present-day ERF aer is close to the lower bound of the multimodel estimate from Bellouin et al. (2020), who report a 68 % likelihood range of −1.6 to −0.6 W m −2 and a 90 % likelihood range of −2.0 to −0.4 W m −2 .
We note that E3SMv1 considers complex interactions between aerosols, clouds, and radiation (see Sect. 2.1.2), so the estimated ERF aer might be larger than those estimated in models with simplified aerosol treatment (Fiedler et al., 2019;Shi et al., 2019).
The large aerosol forcing in E3SMv1 appears to be one of the reasons why the model cannot reproduce the observed global mean temperature evolution in the second half of the 20th century, as suggested by the analysis of Golaz et al. (2019) based on a two-layer energy balance model. The present study complements that analysis by providing a comprehensive investigation of the aerosol forcing characteristics.
This study evaluates the ERF aer simulated by the E3SMv1 Atmosphere Model (EAMv1, Rasch et al., 2019) and relates the ERF aer to characteristics of the aerosol composition and optical properties. We also evaluate the relationships between key aerosol and cloud properties. Our analysis focuses on the following aspects.
-The historical changes in aerosol properties and the regional differences in ERF aer are analyzed in Sect. 3. We investigate whether the historical changes in aerosol properties are consistent with the changes in the aerosol forcing and explore the temporal variations in different regions.
-Second, following previous studies, we identify the most important aerosol, cloud, and radiation flux quantities and analyze the relationship between them (Sect. 4). The analysis focuses on the relative changes in these quantities based on the causality between them.
-Possible causes of the strong ERF aer in E3SMv1 are explored in Sect. 5. The decomposition method of Ghan (2013) is used to quantify the contributions of individual forcing mechanisms (i.e., aerosol-radiation interaction, aerosol-cloud interaction, and residual effect) to the top-of-model (TOM) and surface ERF aer . The contribution of individual aerosol species to the overall ERF aer is quantified by perturbing one type of anthropogenic aerosol emission at a time.
-The sensitivities of ERF aer to parameterization changes, including both formulation changes and parameter tuning, are discussed in Sect. 6. These parameterization changes were applied during the E3SMv1 development and some of them have been further tuned for an alternate configuration of E3SM .
-Section 7 evaluates the impact of anthropogenic aerosols on the near-surface temperature through fast processes (in contrast to the slow processes involving ocean-sea ice responses to aerosol forcing).
The goal of the study is to better characterize and understand the strong ERF aer in E3SMv1 and hence provide insights to help constrain its value in future versions of the model. Descriptions of the model and simulations are provided in Sect. 2. The conclusions from the study are summarized in Sect. 8.

Model and simulations
The model used in this study is the atmospheric component of the E3SMv1, often referred to as EAMv1. A comprehensive description of EAMv1 can be found in , and the representation of aerosol processes is documented in Wang et al. (2020); hence we only provide brief overviews here (Sect. 2.1). The simulations conducted and analyzed in this study are described in Sect. 2.2.

Overview
EAMv1 uses the High-Order Methods Modeling Environment (HOMME) spectral element dynamical core (Taylor and Fournier, 2010;Dennis et al., 2012) and a cubed-sphere mesh. The standard model resolution is used, which has approximately 100 km horizontal resolution, with 30 spectral elements ("ne30"). The model has 72 vertical layers (see Fig.  1 of Xie et al., 2018) and a top at approximately 0.1 hPa, with about 21 levels below 700 hPa and 26 levels above 100 hPa. The key physical processes considered include deep convection (Zhang and McFarlane, 1995), turbulence, shallow convection, and cloud macrophysics that are handled by CLUBB Larson et al., 2002;Larson, 2017), as well as cloud microphysics (Gettelman and Morrison, 2015;Wang et al., 2014), aerosol processes (see details in Sect. 2.1.2), and radiation (Iacono et al., 2008). More information about aerosol and cloud microphysics parameterizations is provided in the next subsection. Most parameterized physical processes use a time step of 30 min, which is also the time step for physics and dynamics coupling. A static (fixed sub-step length) or dynamic sub-stepping method is used by some parameterizations. For example, droplet activation and ice nucleation are coupled with CLUBB and cloud microphysical processes, with a coupling frequency of 5 min. The radiation time step is 1h. More details about process coupling and time stepping can be found in Zhang et al. (2018) and Wan et al. (2021). The evaluations of simulated largescale circulation and clouds in EAMv1 are documented in Rasch et al. (2019), Xie et al. (2018), and Zhang et al. (2019). As summarized in Zhang et al. (2019), the major cloud biases in EAMv1 include underestimated marine stratocumulus clouds at subtropical latitudes near the west coast of major continents, underestimated thin cirrus clouds that result in a smaller longwave cloud radiative effect, excessive low clouds over the Arctic, and too much supercooled cloud liquid fraction.
All of these biases might affect the effective aerosol forcing estimate. EAMv1 is interactively coupled with a land surface model (Oleson et al., 2013), which also includes the dust emission calculation and the parameterization that considers the impacts of light-absorbing impurities (from black carbon and dust deposition) on surface snow-ice albedo (Wang et al., 2020).

Parameterization of aerosols and their interactions with radiation and clouds
Aerosol processes in EAMv1 (Wang et al., 2020) are represented by the Modal Aerosol Module (MAM) (Liu et al., 2012a(Liu et al., , 2016. Major aerosol compositions are considered, including sulfate (SU), black carbon (BC), primary organic matter (POM), secondary organic aerosol (SOA), dust (DU), sea salt (SS), and marine organic aerosol (MOM). The aerosol size distribution is represented by four lognormal modes, with three soluble modes with mixed aerosol species and one primary carbon mode for BC and POM (see Fig. 2 of Wang et al., 2020). Aerosols are assumed to be internally mixed within each mode but externally mixed between different modes. MAM prognoses both the interstitial aerosols (particles suspended in the air and not in cloud droplets) and cloud-borne aerosols (particles in cloud droplets). The model evaluation of present-day aerosol properties in EAMv1 is documented in Wang et al. (2020). The model representation and evaluation of dust and marine organic aerosols are documented in Burrows et al. (2022) and Feng et al. (2022). The aerosol optical properties are parameterized based on the study by Ghan and Zaveri (2007), which approximates the calculation using analytic functions of the mode radius and refractive index of wet particles. The refractive indices for individual aerosol compositions are from Hess et al. (1998), except that the refractive index for BC is taken from Bond and Bergstrom (2006). The volume mixing rule is assumed when calculating the volume-mean wet refractive index for aerosol modes with mixed compositions. The Köhler-theory-based water uptake scheme is used to calculate the wet particle size in each mode based on the predicted aerosol composition and hygroscopicity, as well as the relative humidity with respect to water (RHw) in the clear-sky portion of the grid box. The hygroscopicity values for individual compositions can be found in Table 1 of Burrows et al. (2022). In MAM, an RHw ceiling of 98 % is applied for the aerosol water uptake calculation.
For liquid cloud formation, soluble aerosols can act as cloud condensation nuclei and get activated to form cloud droplets in existing or newly formed clouds. The Abdul-Razzak and Ghan (2000) activation scheme is used, which parameterizes the aerosol activation (or droplet nucleation) based on the subgrid updraft velocity and the aerosol properties (especially for dry size and hygroscopicity) in each mode. Instead of calculating the aerosol activation with a spectrum of updraft velocities as in Ghan et al. (2001), MAM uses a so-called characteristic subgrid updraft velocity to simplify the calculation, which is parameterized as a function of the turbulent kinetic energy (TKE) prognosed by CLUBB. A lower bound of 0.2 m s −1 (inherited from the CAM5.4) is applied to the characteristic subgrid updraft velocity 1 in order to compensate for the potentially underestimated turbulence strength due to unresolved processes (e.g., cloud-top radiative cooling). No minimum cloud droplet number concentration is applied in E3SMv1. The activation process is coupled with an explicit vertical diffusion scheme with dynamical sub-stepping, which handles the vertical diffusion transport of all aerosol species and the cloud droplet number. 2 For mixed-phase clouds, E3SMv1 uses a parameterization based on classical nucleation theory (CNT). The deposition, 1 Note that the 0.2 m s −1 lower bound might be too high to be applied to all conditions, but removing it in E3SMv1 will degrade the cloud simulation if the model is not further tuned. In newer versions of E3SM Golaz et al., 2022), this lower bound is reduced to 0.1 m s −1 , in combination with other tunings in the cloud and turbulence parameterizations.
2 CLUBB calculates the vertical diffusion transport of other quantities, including water vapor, cloud hydrometeors (excluding cloud droplet number), and other trace gases. immersion, and contact freezing of natural dust and BC particles are considered based on a PDF (probability density function) representation of contact angles for those particles (Wang et al., 2014). The initial ice particle size is assumed to be the same as the mean cloud droplet size. Among the different heterogeneous ice formation mechanisms considered in E3SMv1 for mixed-phase clouds, the immersion freezing of dust is the dominant source for ice formation. The contributions from other ice nucleation pathways (the deposition and contact freezing of dust and the immersion freezing of BC) are much smaller. In addition, all the droplets are assumed to homogeneously freeze when T < −40 • C (note that this is different from the homogeneous nucleation of sulfate aerosols), with an initial size of 25 µm. The previously used Meyers scheme (Meyers et al., 1992) for deposition and condensation freezing as well as the contact nucleation scheme (Young, 1974) are optional and are explored as a sensitivity test in this study (see Sect. 6).
For cirrus cloud formation (T < −37 • C), the homogeneous freezing of sulfate aerosols and heterogeneous immersion freezing on mineral dust are considered using the  ice nucleation scheme. For both nucleation mechanisms, an initial ice particle size of 10 µm is assumed. As for droplet nucleation, the ice nucleation scheme also uses the characteristic updraft velocity (as a function of TKE), with a lower bound of 0.2 m s −1 . Sulfate aerosols (or sulfate solution droplets) in the Aitken mode with a diameter larger than a threshold are considered ice-nucleating aerosols for the homogenous ice nucleation. This size threshold was first introduced in the CAM5 model to better reproduce observations (Neale et al., 2010, page 135), and it was set differently in various modeling studies. For example, a threshold size of 100 nm was used in the CAM5 model, while in some sensitivity studies (e.g., Liu et al., 2012b;Zhang et al., 2013;Shi et al., 2015), all Aitken-mode particles are considered to be potential ice-nucleating aerosols in cirrus clouds. E3SMv1 inherited this treatment and further tuned the threshold mainly based on the evaluation of cloud fraction and cloud radiative forcing. For the E3SMv1 CMIP6 experiments and the reference model used in this study, a threshold size of 50 nm is used . We discuss this parameter sensitivity in Sect. 6.

Simulations
Four groups of simulations are analyzed in this study (see Table 1). First, we further analyze the atmosphere-only simulations from Golaz et al. (2019) with prescribed historical sea surface temperature (SST) and sea ice concentration for the years 1870-2014. In the reference simulation ensemble (AMIP), transient-historical anthropogenic and biomass burning-biofuel emissions as well as other forcings (e.g., concentrations of greenhouse gases) are prescribed using the CMIP6 emission data (Hoesly et al., 2018) with additional corrections by Feng et al. (2020). In the second simulation Table 1. List of simulation configurations. Note that for SOA, the 3D formation rates calculated based on Shrivastava et al. (2015) are applied as emissions in E3SMv1.

Simulations Description
Group 1 AMIP Free-running simulation with prescribed SST, 1870-2014, with CMIP6 historical emissions of aerosols and their precursors AMIP1850aer As AMIP, but with CMIP6 pre-industrial (1850)  Group 4 BERG07 (PD and PI) as CTRL or E1850, but the Bergeron process rate tuning parameter is changed from 0.1 to 0.7 following Ma et al. (2022) KK2000 (PD and PI) as CTRL or E1850, but the autoconversion tuning parameters are reset to Khairoutdinov and Kogan (2000) values HOM100 (PD and PI) as CTRL or E1850, but the cutoff size of Aitken-mode sulfate aerosols for homogeneous ice nucleation is changed from 50 to 100 nm MEY (PD and PI) as CTRL or E1850, but the Meyers et al. (1992) heterogeneous ice nucleation scheme is used instead of the CNT-based scheme (Wang et al., 2014). ensemble (AMIP1850aer), the anthropogenic and biomass burning-biofuel emissions are kept at the 1850 level (also from CMIP6). Dimethyl sulfide (DMS) emissions for three time slices (1850, 2000, and 2100) are available from a coupled model simulation, with detailed representation of DMS formation in seawater , and are linearly interpolated for other periods. Dust, sea salt, and marine organic aerosol emissions are calculated interactively in all the simulations analyzed in this study based on the simulated surface winds (frictional velocity for dust) and other surface quantities from the land and ocean components. We take the difference between AMIP and AMIP1850aer to obtain the temporal evolution of anthropogenic aerosol properties, including mass burden, cloud condensation nuclei (CCN) number, and optical depth, as well as the transient ERF aer . Each set of simulations has three ensemble members, and the ensemble-averaged fields are used for analysis. The AMIP-type simulations from Golaz et al. (2019) only include limited monthly mean output for aerosol-related quantities, so additional output data are needed for more detailed analysis. Due to limited computational resources, we rely on shorter simulations that are nudged towards reanalysis data for the additional three simulation groups. Nudged simulations help identify the aerosol forcing signal (i.e., spa-tial pattern) over shorter time periods (e.g, a specific year). See further discussions later in this section.
The second group includes nudged simulations with emissions of anthropogenic and biomass burning-biofuel aerosols as well as their precursors fixed to a specific year. We picked six time slices to represent pre-industrial (E1850), early industrial (E1900), mid-industrial (E1950), European peak industrial (E1970), near-present (E2000), and present (E2010, also as the nudged CTRL) conditions. By taking the difference between E1850 and all other simulations in this group, we derive the anthropogenic aerosol effects in each (emission) period and investigate the regional shifts of the aerosol forcing pattern. In addition to using the CMIP6 emissions, we also performed simulations using the CMIP5 emissions . As shown in Table 2 and Fig. S5 (in the Supplement), using the CMIP5 emissions results in a slightly smaller ERF aer estimate in E3SMv1.
The third group consists of simulations with 1850 aerosol emissions, but with individual aerosol species emission (or 3D production rates for SOA) set to the present-day (2010) value. Again, by taking the difference between E1850 and simulations in this group, the aerosol forcing by individual aerosol species can be derived. For DMS, the global mean difference between DMS emissions in 1850 and 2010 (in-terpolated) is very small (< 1%), and sensitivity simulations show that it has a negligible impact on the simulated ERF aer . Therefore, the result for DMS is not included here. Dust, sea salt, and marine organic aerosols are not considered to be forcing agents, since their emissions are calculated online and/or interactively in the model and are mainly affected by the natural variability and feedback of the perturbed climate and Earth system. The last group includes simulations similar to the nudged CTRL and E1850, but with the cloud microphysics parameterization switched back to older schemes or tuned with a different parameter. Several pairs of simulations were performed with both present-day (PD) and pre-industrial (PI) aerosol emissions to estimate the anthropogenic aerosol effects.
Since the large-scale circulation is constrained, nudging can help distinguish the signals (e.g., caused by changes in external forcings) from the noises caused by chaotic model responses to forcing or parameterization changes (Kooperman et al., 2012;Zhang et al., 2014). We only nudge horizontal winds in our simulations because it is better not to constrain the thermodynamical response of meteorological fields to aerosol effects. Also, earlier studies show a negative impact of temperature and humidity nudging on the modeled climate (Zhang et al., 2014;Lin et al., 2016;Sun et al., 2019), which will make the forcing estimate inconsistent with that estimated using free-running simulations. For all the nudged runs, we nudge the horizontal winds towards the ERA-Interim reanalysis data, with a relaxation timescale of 6 h. Since the vertical interpolation-extrapolation of reanalysis data to model vertical grids might have large errors near the surface, nudging is not applied in the lowest two layers. In the layer next to the no-nudging layers, the nudging relaxation time is doubled (weaker nudging compared to the layers above with full nudging strength). Each nudged simulation is 15 months long and the results from a whole year (after a 3-month spin-up) are used for analysis.
We note that since the horizontal winds are constrained in the nudged simulations, the anthropogenic aerosol effect on the large-scale circulation is greatly diminished. As we will show later, for global mean or regional mean (over a large area) estimates, the nudged simulations provide similar ERF aer estimates as the free-running simulations.

Historical changes and regional shifts of aerosol compositions, optical depth, and forcings
Before showing the effective aerosol forcing, we first analyze the simulated historical changes in anthropogenic aerosol burdens (vertical integration of aerosol mass) over the simulation period . Here, we define the anthropogenic change in a certain quantity (denoted as , e.g., aerosol burden, optical depth, and radiative fluxes) as the differences between the modeled quantities in simulations with transient historical forcing and the pre-industrial forcing (year 1850) for aerosols and their precursors. Note that the anthropogenic change we define here also includes the contribution from biomass burning emissions, which is consistent with the IPCC (Intergovernmental Panel on Climate Change) definition of "anthropogenic" in a radiative forcing context. Figure 1 shows the historical changes in global and regional mean anthropogenic aerosol burden ( burden) for different compositions. In addition to the results from AMIP-type simulations, we also show the nudged simulations for 5 selected years (shown as dots). Overall, the results from the nudged simulations agree very well with those from the AMIP-type simulations in terms of both historical variations and regional differences of burden. For aerosols with anthropogenic sources (i.e., sulfate, BC, POM, SOA), the global mean anthropogenic burden ( Fig. 1a) shows an overall increase during the simulation period (1870 to 2014). The burden of anthropogenic sulfate shows largest change, followed by that of SOA and POM (note that the BC burden is scaled by a factor of 5 in Fig. 1). There are large regional differences in burden. For example, the anthropogenic sulfate burden in the Northern Hemisphere (NH) polar region (60-90 • N) peaks in the 1970s and decreases by about 50% in the year 2010 compared to the peak value ( Fig. 1b). This is consistent with the emission reduction of air pollutants in Europe and North America where the legislation for clear air was implemented from the 1970s (Sliggers and Kakebeeke, 2004), while in the tropics (20 • S-20 • N), the sulfate burden continues to increase (Fig. 1d).

Changes in anthropogenic aerosol burden and optical depth
The carbonaceous aerosols are not only affected by the anthropogenic aerosol sources (e.g., from industrial and transport sectors), but also by the biomass burning emissions. In addition to the gradual increase in time, there are several intermittent spikes in the anthropogenic aerosol burden change. The change in dust aerosol burden caused by the anthropogenic aerosol effect has large interannual variabilities, and there is a small decreasing trend after about 1950, indicating that anthropogenic aerosols can indirectly affect the dust life cycle, especially in the tropics and NH midlatitude region (Fig. 1d, c). To understand what causes the decreasing trend in the dust aerosol burden, the dust mass budget is eval- uated using the nudged simulations (the impact of interannual variability in the AMIP simulations can be avoided). We find that the small decreasing trend is mainly caused by slightly weakened dust emission in simulations with increased anthropogenic emissions. Compared to the simulation with 1850 emissions, the global mean dust emission rate decreased by 2 %-3 % in the simulations with emissions after 1970. In the nudged simulations, we only weakly constrain the large-scale horizontal winds, so the near-surface winds can still be affected by stability changes in the lower troposphere. The slightly weakened dust emission is very likely due to the changed surface winds caused by anthropogenic aerosols through changes in the atmospheric stability in the lower troposphere (Jacobson and Kaufman, 2006;Baró et al., 2017). The dust wet removal rate also decreases with increased anthropogenic emissions, suggesting that the impact of anthropogenic aerosols on dust wet scavenging (through coating) is less important in our model compared to the dust emission changes.
Changes in sea salt and marine organic aerosol burdens are small (only evident in SH midlatitude and polar regions) and also have large interannual variations. Both sea salt and marine organic aerosol emissions are dependent on surface winds, SST, and the ocean fraction. The large interannual variations (from AMIP simulations) are mainly caused by changes in surface winds in combination with sea ice concentration and SST changes in different years. For the nudged simulations, since the large-scale winds are constrained and single-year SST and sea ice concentrations are used, the changes in sea salt and marine organic aerosol burden are very small. Figure 2 shows the historical changes in global and regional mean anthropogenic aerosol optical depth ( τ a ) for the total and individual aerosol components. Again, we see consistent results between the estimates by AMIP-type simulations and nudged simulations constrained to the same largescale wind fields but with different years of emission input. The total τ a (black) shows an overall increase during 1870-2014, but with large interannual variations. The interannual variability in total τ a is caused by changes in dust aerosols and by changes in carbonaceous aerosols (BC, POM, SOA) affected by fluctuations in biomass burning emissions and by changes in dust aerosols. Compared to the multi-model estimates for τ a (0.015-0.04) by Bellouin et al. (2020) based on multiple global model estimates (vertical bar in Fig. 2a), the estimated present-day τ a in E3SMv1 (about 0.03 around the year 2010) is well within the range.
τ a of individual aerosol compositions shows similar temporal variations and regional differences as seen in burden (Fig. 1). It can be clearly seen that after 1970, both sulfate burden and sulfate τ a over NH polar and NH midlat-itude regions decrease significantly, while for carbonaceous aerosols, there is an overall increase in burden and τ a . Figure 3 provides a more quantitative comparison of the contributions of sulfate, BC, POM, and SOA to the total anthropogenic τ a (excluding dust and sea salt). The contribution of carbonaceous aerosols to total anthropogenic τ a increases from 27 % (BC 8 %, POM 6 %, and SOA 13 %) in 1970 to 48 % (BC 11 %, POM 13 %, and SOA 24 %) in 2000 and 50 % (BC 12 %, POM 10 %, and SOA 28 %) in 2010. Since carbonaceous aerosols are less hygroscopic compared to sulfate in the model, their increased contribution to aerosol mass (burden) will decrease the mean hygroscopicity of the aerosol particles and thus affect the optical property and the ability to be activated to cloud droplets. As will be shown later, the changes in the relative contribution of different types of anthropogenic aerosol components will have an impact on the relationship between τ a , CCN (cloud condensation nuclei) concentrations, and droplet concentrations (N d ). Figure 4 shows the historical changes in global and regional mean ERF aer ( F) for the net (dark grey), shortwave (blue), and longwave (red) components. The global mean shortwave F is mostly negative and the largest negative forcing (about −2.6 W m −2 ) first appears around 1980. In contrast, the longwave F is overall positive, with the largest forcing of about 1.0 W m −2 . After 1980, the global mean F does not show a clear increasing or decreasing trend, but there are fluctuations caused by interannual variations. Similar to the burden and τ a results, we can see large regional differences in the temporal evolution of F. The overall increasing or decreasing trend in F resembles the changes in sulfate burden and τ a , especially for NH midlatitude and polar regions.

Historical changes in ERF aer
Compared to F estimated by CMIP6 RFMIP models , which adopt year 2014 emissions as present-day conditions, the shortwave and longwave F estimated by E3SMv1 are barely within the max-to-min range in the CMIP6 RFMIP model estimates (thin vertical bars) but outside the 1 standard deviation around the mean estimate (thick vertical bars). For the net F, the E3SMv1 estimate is even larger than the max-min range of the multi-model estimates, since the large shortwave and longwave F estimates from CMIP6 are from the same model and the net change is relatively small.
As for burden and τ a , the F values estimated by the nudged simulations (dots in Fig. 4) are very similar to the estimates from the AMIP-type free-running simulations (lines). Therefore, in the following we mainly use results from the nudged simulations to discuss the regional forcing shifts, forcing decomposition (Sect. 5), and sensitivity to parameterization changes (Sect. 6). Note that since the large-scale circulation is constrained, the estimated ERF aer is mainly af-fected by fast processes associated with interactions between aerosol, cloud, precipitation, and radiation. Figure 5 shows the geographical distributions of F for different time slices estimated based on several pairs of nudged simulations (group 3 in Table 1) with emissions set to 1850 or a specific year (1900, 1950, 1970, 2000, or 2010). The corresponding changes in aerosol optical depth (AOD) and column-integrated number concentrations of CCN (at 0.1% supersaturation) and cloud droplets are shown in Figs. S8, S9, and S10.
In 1900, relatively large F appears over western and northern Europe, North America, and the downwind regions. We see an overall increase in the magnitude of F in 1950 (about a factor of 3 larger than the 1900 value), but the regions with large F are still very similar as in 1900. The F over western and northern Europe and North America peaks in 1970. Afterwards, the largest F gradually shifts to East Asia. The estimated F values for the years 2000 and 2010 are very similar, although the forcing over western and northern Europe and North America is further reduced in 2010. Figure 6 clearly shows that F has strong seasonal variations. The global mean F is larger in boreal summer (June-July-August) and spring (March-April-May), mainly due to stronger oxidations of SO 2 and higher sulfate concentrations in these seasons. Also, both the annual and seasonal averaged zonal mean F fields have a regional shift of larger aerosol forcing from NH high-latitude regions in 1970 to the tropics and subtropics in 2000 and 2010, mainly due to regional changes in anthropogenic emissions. The largest changes in the NH polar region appear in boreal summer, from about −8 W m −2 in 1970 to −4 W m −2 in 2010.
In winter, shortwave F in the NH polar regions is close to zero due to very small incoming shortwave radiation. As is shown later, the increased liquid water path by anthropogenic aerosol effects results in a positive (longwave) radiative forcing at the surface, which causes surface temperature to increase in addition to the greenhouse warming.

Relationships between changes in aerosol, cloud, and radiative flux quantities
The impacts of aerosol perturbations on the top-of-model radiative flux changes involve many processes related to aerosol and cloud microphysics. Previous studies (Quaas et al., 2009;Ghan et al., 2016;Bellouin et al., 2020) have summarized the key causal relationships that determine the aerosol effect through interactions with radiation and clouds. These relationships are often expressed in the logarithmic form of d ln Y /d ln X, since most processes involved in aerosol-cloud interactions are more sensitive to a relative change (between pre-industrial and present-day conditions) in quantity than the absolute change (Carslaw et al., 2013;Bellouin et al., 2020). We follow the conventional expressions as used in Ghan et al. (2016) and Bellouin et al. (2020). In the following discussions, ln X is the relative change in X (spatial and temporal mean) between PD and PI simulations, and it is calculated as 3 X/X PD . In earlier studies, some aerosol and cloud quantities are sampled at certain locations or under selected conditions based on their role in the physical processes or the consideration for comparing with satellite retrievals. For example, the cloud-base values of CCN (e.g., 1 km above the surface) and the cloud-top values of cloud droplet number concentrations (N d ) are often used (e.g., Quaas et al., 2009;Ghan et al., 2016;Zhang et al., 2016). Furthermore, to focus on the analysis on certain types of clouds, in some studies (e.g., Ghan et al., 2016;Gryspeerdt et al., 2020) these quantities are conditionally sampled only when some criteria are met  Table 1, and "E2010" represents the CTRL simulation. See Sect. 3.2 for details.
(e.g., when cloud-top temperature is larger than a threshold value). For such analysis, high-frequency output data for the related aerosol and cloud properties are needed unless the model can conditionally sample data online and process them accordingly .
It should be noted that in our AMIP simulations, not all the quantities needed for such analysis are available or available at the required frequency (e.g., 3-hourly). Therefore, for some quantities (such as for CCN and N d ), we use the available column-integrated quantities sampled as monthly mean fields for the analysis using AMIP simulations. For the nudged simulations, both monthly mean output and high-frequency data sampled under required conditions (as in Ghan et al., 2016) are available, so we use both and evaluate whether the results are sensitive to the data sampling methods.
4.1 Relationships between ∆ ln τ a (or ∆τ a ) and other quantities Figure 7 shows the relationships between simulated changes in AOD ( τ a , ln τ a ) and key cloud and radiative flux properties. The global annual mean total AOD (τ a , Fig. 7a) from the AMIP simulation increased (almost linearly) from about 0.110 to 0.145 during the simulation period. Compared to the AMIP simulation (small dots), τ a is consistently smaller in the nudged simulations with anthropogenic emissions for different years (shown as bigger markers). This is mainly caused by the slightly weakened dust emission when reanalysis nudging is applied (Timmreck and Schulz, 2004;Sun et al., 2019). Although we do not apply nudging in the surface layer, the simulated surface winds and dust emission fluxes are still affected by the nudging in upper layers (see Fig. 18 in Sun et al., 2019). The estimated present-day (year 2010) τ a is well within the range of multi-model estimates (0.12-0.16 from Ghan et al., 2016, and0.13-0.17 from Bellouin et al., 2020, shaded area in Fig. 7). The simulated anthropogenic AOD ( τ a ) shows a similar increase during the simulation period. The estimated present-day τ a (∼ 0.03) in E3SM is also well within the range of multi-model estimates.
The relationship between ln τ a ( τ a /τ a ) and lnCCN 0.1% is shown in Fig. 7b. They have a quasilinear correlation before 1970 (grey dots). As indicated by the changed slope, the increase in lnCCN 0.1 % is weakened as ln τ a further increases after 1970 (red dots). Results from the nudged simulations show very similar changes. When sampling the data conditionally (cloud-top temperature > −10 • C) at high frequency, the result (shown as an orange square in Fig. 7b, c, d) is similar to that derived using the monthly mean data available, and it is within the multi-model estimated range in Ghan et al. (2016). ln τ a and lnN d also show a quasi-linear correlation before 1970, but the slope changes sign for samples after 1970. This is likely related to the fact that there is a continuous increase in anthropogenic carbonaceous aerosols in most regions (that increases τ a ) but decreased or stabilized anthropogenic sulfate aerosols in the NH high-latitude and midlatitude regions. We note that the N d changes in response to CCN changes are different between individual regions (see Figs. S9 and S10 and discussions in Sect. 4.2). In addition, the droplet nucleation is also dependent on the maximum supersaturation in a rising cloud parcel, which is parameter- ized based on the ambient meteorological conditions (e.g., vertical mixing and subgrid updraft velocity).
Following Bellouin et al. (2020), Fig. 7d shows the relationship between τ a and the changes in clear-sky shortwave radiation ( F SW,clear ). F SW,clear magnitude increases with increasing τ a almost linearly, but the increase is much weaker after 1970. This is likely also related to the fact that there are relatively more carbonaceous aerosols but fewer and/or stabilized sulfate aerosols, which increase the absorption and offset the cooling effect by scattering. The estimated F SW,clear in E3SMv1 is well within the multi-model estimates presented in Bellouin et al. (2020). Figure 8 shows the relationships between simulated relative changes in column-integrated CCN at 0.1 % supersaturation ( lnCCN 0.1 % ) and in cloud microphysical quantities, including cloud droplet number concentration ( lnN d ), liquid water path of stratiform clouds ( lnLWP strat ), and cloud-top droplet effective radius ( lnRe liq ). The global annual mean CCN 0.1 % increases by more than 80 % from the pre-industrial era to present-day conditions. Compared to the year 1970, the global annual mean CCN 0.1 % increases by more than 60 % in the year 2010. However, it shows opposite trends in the NH polar region and in the tropics (see Figs. S1 and S3).

Susceptibility of N d , LWP, and R e,liq to changes in CCN
Before 1970, there is strong positive correlation between lnN d and lnCCN 0.1 % (grey dots in Fig. 8b) and between lnLWP strat and lnCCN 0.1 % (Fig. 8c) but a negative correlation between ln Re liq and lnCCN 0.1% (Fig. 8d). This is consistent with the causal relationship between these quantities simulated in the model (more CCN leads to higher N d and LWP but smaller Re liq ) and the fact that all types of anthropogenic aerosol emissions increase (so does CCN) from the pre-industrial era to about 1970 (see Fig. 1).
After 1970, the correlation between these global quantities become weaker (red dots in Fig. 8b-d). This is mainly caused by three factors. First, more primary (carbonaceous) aerosols and less secondary (sulfate) aerosols decrease the particle hy-groscopicity and size, which subsequently reduce the droplet nucleation even if CCN 0.1 % does not change much (e.g., Figs. S1b and S2b). The impact of carbonaceous aerosol increase is most evidently seen from the differences between 2000 and 2010 (Fig. 9a, b). lnCCN 0.1 % slightly increases in the NH polar and NH midlatitude regions, while lnN d decreases in both regions.
Second, there are regional differences in the simulated relationships. The large regional differences in lnN d and lnCCN 0.1 % contribute to the temporal variations of global mean values. In the NH polar region, CCN 0.1 % increases before 1970 but decreases afterwards (Figs. 9a and S1a). In midlatitudes, CCN 0.1% is stabilized after 1970 as a result of the increase in Asia and the decrease in North America. In the tropics, CCN 0.1 % increases continuously during the whole simulation period (Figs. 9a and S3a). Consequently, there is an increase in the global mean lnCCN 0.1 % . In contrast, the global mean lnN d is greatly affected by the decreases in both the NH polar and NH midlatitude regions, so it decreases slightly after 1970. Although lnN d still increases continuously in the tropics, the value is much smaller than those in the NH polar and NH midlatitude regions.
Third, the CCN-N d relationship also appears to have a "saturation effect" (Martin et al., 1994;Boucher and Lohmann, 1995;, which means a saturation in N d with increasing aerosol number or CCN. The saturation effect can be caused by (1) reduced wet growth by carbonaceous aerosols (Martin et al., 1994), which has a similar effect as the first factor, and (2) a suppression in the maximum supersaturation (that the air parcel can have) caused by more CCN; this will lead to a reduced fraction of CCN that is activated (Chuang et al., 2000). As shown in Fig. 9c, the slope ( lnN d / lnCCN 0.1 % ) calculated using regional mean fields slightly decreases from 1900 to 2010 for all the inspected regions.
Compared to the multi-model estimates reported in Ghan et al. (2016), the estimated values of lnCCN 0.1 % , ln N d , and lnLWP strat in E3SMv1 are all well within the min-tomax range (almost right in the middle). However, we find that the estimated ln Re liq magnitude in E3SMv1 is relatively large compared to the multi-model estimates.
4.3 Susceptibility of Re liq , LWP, τ liq , and cloud radiative forcing to changes in N d Figure 10 shows the relationships between simulated relative changes in column-integrated N d ( ln N d ) and in other cloud microphysical quantities, including cloud-top droplet effective radius ( ln Re liq ), liquid and ice water path of stratiform clouds ( lnLWP and lnIWP), cloud optical depth ( ln τ liq ), liquid cloud fraction ( ln f liq ), and shortwave cloud radiative forcing ( lnCF SW ). For the sampling using monthly model output, ln N d , lnLWP, and lnIWP (for stratiform clouds) are calculated using grid mean values. ln Re liq and ln τ liq are calculated using in-cloud values, which are sampled using the COSP (Cloud Feedback Model Intercomparison Project Observation Simulator Package) simulator and derived during post-processing. For the sampling using high-frequency data (orange square in Figs. 7-10) following Ghan et al. (2016), they are in-cloud values calculated using the grid box mean values divided by the cloud cover after averaging in time and space. Because they are relative changes in global annual mean quantities, using in-cloud values (derived from grid mean values and cloud cover after averaging in time and space) or grid average values gives very similar results. The global annual mean N d increases by about 50 % from the pre-industrial era to present-day conditions (blue square in Fig. 10a). Compared to the year 1970, the global annual mean N d is smaller in the year 2010, which is different from the results we see in τ a (Fig. 7a) and CCN 0.1 % (Fig. 8a). More aerosols simulated in the model (in terms of global annual mean) do not necessarily lead to an increase in N d , since there are regional decreases in the sulfate burden that reduce the aerosol hygroscopicity and subsequent aerosol activation. On the other hand, similar opposite trends (as for lnCCN 0.1 % ) appear in the NH polar region (decreases after 1970) and in the tropics (increase continuously after 1970).
There is a strong negative correlation between ln N d and ln Re liq (grey and red dots in Fig. 10b). There is also strong positive correlation between ln N d and other quantities including lnLWP, ln τ liq , and lnCF SW (Fig. 10c, d, f). The results are consistent with the causal relationship between these quantities simulated in the model (more N d leads to larger LWP but smaller Re liq ). Unlike lnCCN 0.1 % , after 1970, the correlation between N d and other global quantities is very similar to those before 1970. Since the relations shown in Fig. 10 characterize cloud properties (instead of the relations between aerosols and clouds), the lack of a change in slope (as in Fig. 8) is expected.

Relationships between changes in liquid and ice cloud properties
As described in Sect. 2.1.2, E3SMv1 considers the activation of aerosol particles to cloud droplets and the freezing of cloud droplets to ice crystals in mixed-phase clouds when ice nuclei are available. It also considers the homogeneous ice nucleation of sulfate aerosols and heterogeneous ice nucleation of dust aerosols in cirrus clouds. Because dust is natural and the ice nucleation rate caused by BC is rather small (compared to the contribution by dust), sulfate is the dominant anthropogenic aerosol type that affects both the droplet (and the subsequent freezing) and ice crystal formation processes. Figure 11 shows the relationships between simulated relative changes in liquid and ice cloud properties. There are strong positive correlations between the changes in the liquid and ice cloud quantities for the global annual mean columnintegrated particle number ( ln N d vs. ln N i ), cloud-top effective radius ( ln Re liq vs. ln Re ice ), column-integrated water path in stratiform clouds ( lnLWP vs. lnIWP), and cloud optical depth ( ln τ liq vs. ln τ ice ). The ice properties shown here do not consider the contribution from snow. On the other hand, the relative changes in liquid and ice cloud fractions are not well-correlated. Compared to liquid clouds, the ice cloud properties have a similar but slightly weaker sensitivity to changes in anthropogenic aerosols (Fig. S4).
The relative change in global annual mean columnintegrated droplet number concentration ( ln N d in Fig. 11a, also in Fig. 10) increases from the pre-industrial era but starts to decrease after 1970. This is not the case for the relative changes in the column-integrated ice crystal number, for which we see a continuous increase. This difference is mainly caused by the fact that the major changes in N d and N i appear in different regions. As shown in Fig. 12, the largest changes in N d are in NH midlatitude and polar regions, where anthropogenic aerosols affect the low liquid cloud most. In contrast, the largest changes in N i are mainly in the tropics and subtropics regions, where the cirrus clouds are affected by anthropogenic sulfate particles through homogeneous ice nucleation. In our model, the homogeneous ice nucleation in cirrus clouds is parameterized based on  and Gettelman et al. (2010). We use the aerosol number concentration in the Aitken mode truncated at a specified cutoff diameter (50 nm in CTRL) to estimate the number of sulfate aerosol particles that can initiate homogeneous ice nucleation. Similar to findings in Gettelman et al. (2010), the parameterized homogeneous ice nucleation in Figure 11. As Fig. 7, but for relationships between changes in liquid and ice cloud properties as well as in SW and LW cloud radiative forcing. ln Re ice , lnIWP ice , ln τ ice , lnf ice , and lnCF LW are the relative changes in the effective radius of cloud ice at the cloud top, ice water path in stratiform clouds, ice cloud optical depth, ice cloud fraction, and longwave cloud radiative forcing, respectively. E3SMv1 is sensitive to the sulfate aerosol number, and there is a large difference between PD and PI simulations (this is further discussed in Sect. 6). The anthropogenic aerosol effects on N i below 300 hPa (mainly in mixed-phase clouds) are very small.

Forcing decomposition
In this section, we decompose the effective aerosol forcing based on the method proposed by Ghan (2013) and based on the individual aerosol forcing and emission types. These decompositions help us understand the role of various forcing mechanisms and the impacts of individual aerosol compositions on the overall ERF aer .

Decomposition of ERF aer due to the aerosol-radiation interaction, aerosol-cloud interaction, and surface albedo effect
We use the Ghan (2013) method to decompose the total ERF aer ( F) into individual forcings caused by aerosolradiation interactions (ARI, direct aerosol effect), aerosolcloud interactions (ACI, mainly from indirect aerosol effect), and the residual term (mainly from surface albedo effect). To do this, an additional call to the radiation calculation is added in the model to calculate the radiative fluxes without considering the aerosol scattering and absorption of shortwave radiation in the calculation for all-sky (F clean ) and clear-sky (F clear,clean ) conditions. The decomposed forcings are calculated as follows.
-Aerosol-radiation interaction ("direct radiative forcing" in Ghan, 2013): -Aerosol-cloud interaction ("cloud radiative forcing" in Ghan, 2013): (F clean − F clear,clean ) -Residual forcing ("surface albedo forcing" in Ghan, 2013): F clear,clean We note that F and F clean are calculated in a single simulation. They both include the impact of ACI, which is considered the same for the two terms in the decomposition calculation. In addition to the indirect aerosol effect, the ACI term ("cloud radiative forcing" in Ghan, 2013) also includes the semi-direct aerosol effect (light-absorbing aerosols can heat the atmosphere and affect clouds), although the estimated ef- fect is much smaller than the indirect aerosol effect . Figure 13 shows the geographical distributions of the decomposed top-of-model F for the year 2010 using the nudged simulations with E3SMv1 (CTRL -E1850). The forcing caused by aerosol-cloud interactions has the largest contribution to the overall F. In contrast, the global annual mean forcing caused by aerosol-radiative interactions is very small (0.04 W m −2 ), but it is positive over subtropical land and stratocumulus-prevalent regions, where BC absorption of shortwave radiation is strong (especially above clouds). The large positive forcing due to aerosol-radiative interactions off the Peruvian coast is mainly caused by absorption enhancement caused by the light-absorbing aerosols (mainly by biomass burning aerosols) above highly reflective clouds (Haywood and Shine, 1997). The clear-sky aerosol forcing is overall negative, except over North Africa and the Arabian Peninsula, where there are significant changes in absorbing aerosols. The residual forcing ( F clear,clean ) is negative over the NH polar region, which is mainly caused by anthropogenic sulfate aerosols. Positive forcing appears over subtropical land areas and is caused by both shortwave (surface albedo effect) and longwave (water vapor) changes.
At the surface, the anthropogenic aerosol effect through the interactions with radiation contributes to more than 40 % of the total forcing (Fig. 14g). Large surface longwave ACI appears in the NH high-latitude and polar regions (Fig. 14f), which is caused by more low liquid cloud induced by anthropogenic aerosols. The thicker liquid clouds can trap more outgoing longwave radiation and emit it back to the surface. As discussed in Sect. 7, this causes a strong increase in the near-surface temperature during the boreal winter season.

Decomposition of forcing by individual compositions
To investigate the impact of changes in individual aerosol compositions, we decompose F based on the different anthropogenic aerosol emission types that contribute to the total forcing. To do this, we change one type of aerosol emission in E1850 at a time from the pre-industrial condition to the 2010 level so that the contribution of this single emission type can be estimated. All other aerosol emissions are set to pre-industrial conditions. Note that we take the 1850 aerosol condition as the background. Due to the nonlinear effect (e.g., BC mixing state in 1850 is different from that in 2010), the sum of the estimated forcings by individual aerosol species would not be the same as the total forcing estimated using CTRL and E1850. Figure 15 shows the geographical distributions of the decomposed F for individual aerosol emission types. The Figure 13. Geographical distributions of decomposed net (left column), shortwave (middle column), and longwave (right column) F at the top of the model (TOM). The decomposition calculation follows Ghan (2013). ALL indicates the total F calculated from the difference between CTRL and E1850 (ALL = ARI + ACI + RES). ACI indicates the F caused by aerosol-cloud interactions (second row), ARI the F caused by aerosol-radiation interactions (third row), and RES (bottom row) the residual forcing. The clear-sky direct aerosol effect (fourth row) is also shown. See Sect. 5.1 for details.
anthropogenic sulfur emission has a large impact on the F (Fig. 15d). The global annual mean net forcing is about −1.66 W m −2 , which is very close to the F estimate due to all aerosol emission changes. However, both the shortwave and longwave F magnitudes increase significantly (by 0.25 W m −2 ), and they compensate for each other. The anthropogenic BC emission has an overall positive F (Fig. 15g), with a global annual mean net forcing of 0.27 W m −2 . The positive F over land is mainly caused by the shortwave absorption of BC (Fig. 15h).
The compensating changes over the Tropical Warm Pool region ( F SW and F LW ) are mainly caused by reduced high clouds in that area. The impacts of the anthropogenic POM and SOA emissions on F are overall negative (−0.3 to −0.4 W m −2 , Fig. 15j, m). The impacts of biomass burning emission changes (BB) on F are mostly in the shortwave, and they have a distinct regional pattern (Fig. 15m, n). Over Siberia, BB F is positive (mainly in JJA and MAM), while over the west of Columbia, Ecuador, and Peru, there is a significant negative forcing (in JJA and SON). This suggests that there are large differences of BB emissions near (or on the upwind side) these two regions between the preindustrial (1850) and present-day (2010) emissions. Figure 14. Geographical distributions of decomposed net (left column), shortwave (middle column), and longwave (right column) F at the surface. The decomposition calculation follows Ghan (2013). ALL indicates the total F calculated from the difference between CTRL and E1850 (ALL = ARI + ACI + RES). ACI indicates the F caused by aerosol-cloud interactions (second row), ARI the F caused by aerosol-radiation interactions (third row), and RES (bottom row) the residual forcing. The clear-sky ARI (fourth row) is also shown. See Sect. 5.1 for details.

Sensitivity of effective aerosol forcing estimation to parameterization changes
Based on the sensitivity simulations performed during the model development, we have identified several parameterization changes that have important impacts on the anthropogenic aerosol effect in the model. These include changes in cloud microphysics that affect the autoconversion in warm clouds, the Bergeron process and heterogeneous ice nucleation in mixed-phase clouds, and homogeneous ice nucleation in cirrus clouds (see Table 1). As described in Rasch et al. (2019), E3SMv1 uses a modified autoconversion scheme based on Khairoutdinov and Ko-gan (2000) (hereafter KK2000), wherein the parameters have been tuned so that the autoconversion process converts less liquid to rain compared to the standard KK2000 treatment. As a result, the liquid water content increases less rapidly (especially under more pristine conditions) when the aerosol and cloud droplet number concentrations are small. With the parameters reset to KK2000 standard values, the liquid water path is decreased significantly in most regions (Fig. S6e). In contrast, the changes in liquid water path between the PD and PI simulations become larger (Figs. 16c and S6f). Consequently, the simulated anthropogenic aerosol effect (−2.03W m −2 ) is much larger than the reference model (Table 2 and Fig. 16). The Bergeron process rate determines how much liquid is converted to ice through the Wegener-Bergeron-Findeisen process. In E3SMv1, a scaling factor was set to make the Bergeron process rate 10 % of its nominal value. This was originally intended to compensate for an overestimated ice formation rate when using the Meyers scheme (Meyers et al., 1992). The overly small Bergeron process rate produces too much supercooled liquid cloud in the simulation. By changing the scaling factor from 0.1 to 0.7, the liquid and ice water path is decreased in most places, especially in midlatitude and polar regions (Figs. S6c and S7c). The changes in liquid and ice water paths between the PD and PI simulations also become smaller (Figs. S6d and S7d). As a result, the simulated anthropogenic aerosol effect is significantly weaker than the reference model for both the shortwave and longwave components (Table 2 and Fig. 16). E3SMv1 uses the CNT-based heterogeneous ice nucleation scheme to consider the ice formation in mixed-phase clouds (Wang et al., 2014). Compared to the Meyers scheme (Meyers et al., 1992) used in E3SMv0, the CNT scheme predicts less heterogeneous ice nucleation and freezes less supercooled liquid (Figs. 16c and S6h), especially in the NH polar region. The increased supercooled liquid amount leads to optically thick liquid clouds in the lower troposphere, which subsequently increase the downward longwave radiation and warm the surface in the NH polar region (this is further explained in the next section). A much larger impact has been seen in the high-resolution simulations , which reported severely overestimated Arctic surface temperature when the CNT scheme was used along with the 0.1 Bergeron process rate scaling factor. This might be related to the fact that the high-resolution model can better resolve the heterogeneity in cloud phase (Tan and Storelvmo, 2016). As mentioned in Caldwell et al. (2019) and above, it is more reasonable to use a larger scaling factor (e.g., 0.7 instead of 0.1) when the parameterization uses the CNT scheme. When the Meyers scheme is used, the relative changes in the column-integrated mass and number are reduced for both cloud liquid and ice, leading to a reduced ERF aer in both the shortwave and longwave components. Similar sensitivities were found in the CAM6/CESM2 model simulations .
In E3SMv1, the threshold size of Aitken sulfate aerosols for homogeneous ice nucleation was adjusted to 50 nm in the DECK simulation setup from the original 100 nm value used in CAM5. This allows more Aitken-mode sulfate particles to nucleate ice homogeneously, which greatly increases the ice crystal number concentration and their subsequent growth. When the threshold size is reset to 100 nm, the ice crystal number and ice water content decrease by about a factor of 2 in most regions (Fig. S7i), and the changes in ice crystal number and ice water content between the PD and PI simulations become smaller (Figs. 16d and S7j). As a result, the simulated anthropogenic aerosol effects in both the shortwave and longwave are reduced compared to the reference model (Table 2). Due to the strong compensation between the changes in SW and LW, the change in the net effect is fairly small (−1.62 vs. −1.64 Wm −2 in the reference model).
We note that the sensitivities discussed above are likely dependent on the quantity of simulated background (preindustrial) aerosols and cloud hydrometeor concentrations in the model. For example, the autoconversion rates estimated by the original and modified KK2000 parameterizations depend on the droplet number concentrations, and the difference between them is large under pristine conditions (in which the cloud droplet number is small, see

Impact on near-surface temperature over land and ice-covered regions through fast processes
The near-surface temperature is strongly affected by the radiative heating and cooling caused directly or indirectly by anthropogenic aerosols through fast processes. The "fast processes" defined here are processes that do not involve the ocean feedback to the atmosphere (in response to the anthropogenic aerosol forcing). Our definition of "fast processes" is slightly different from the "rapid adjustment" by Sherwood et al. (2015), which means "changes that occur directly due to the forcing, without mediation by the global-mean temperature". As pointed out by Sherwood et al. (2015), the rapid adjustment may include fast changes in SST patterns caused by SST changes in areas with relatively shallow mixed layers. In this case, the global mean temperature change is negligible, but these SST changes can significantly affect the cloud and circulation patterns. Such rapid adjustments are not considered to be fast processes for our case.
Since the sea surface temperature is prescribed in our AMIP and nudged simulations, the impact of anthropogenic aerosols on sea surface temperature and the associated slow processes, such as the ocean circulation and its feedback to atmosphere, cannot be considered. Nevertheless, it is useful to estimate the impact of anthropogenic aerosols on the surface temperature through fast processes. The impact through slow processes can be inferred if both the AMIP and coupled simulation results are available.
We first look at the results from the nudged simulations. Since the large-scale circulation and SST are constrained, the local changes in diabatic heating and cooling (affected mainly by fast processes) have a dominant influence on the near-surface temperature over land. Fig. 17 (first row) shows the near-surface temperature (at 2 m) changes ( TREFHT) caused by the emission changes between preindustrial and present-day conditions. As expected, near the emission sources and the downwind regions (especially in NH midlatitudes), anthropogenic aerosols in general have strong cooling effects. However, over the NH polar region, we see warmer surface temperature associated with re- Figure 17. Global annual and seasonal mean changes in near-surface temperature at 2 m ( TREFHT, K), decomposed ACI forcing at the surface ( F ACI,SRF , W m −2 ), and downward longwave radiation flux ( FLDS, W m −2 ) due to the anthropogenic aerosol effects. See Sect. 7 for details. sponses over sea ice from increased anthropogenic aerosol emissions.
The increase in near-surface temperature mainly happens in the boreal winter season (Fig. 17b) and in late autumn and early spring, when the shortwave radiation is weak (due to the large solar zenith angle) and the longwave radiation plays a more important role in determining the near-surface temperature. Between simulations using 2010 and 1850 emissions, there is about a 0.5-2 K increase in December-January-February mean near-surface temperature over land and areas covered by sea ice near the Arctic. During the boreal summer season, the shortwave cooling effect is usually stronger than the longwave warming effect. Also, during the summer season the sea ice cover is smaller than in winter, and the near-surface temperature is less affected by the surface radiation changes. Therefore, the overall effect is cooling over land areas in the NH polar region (Fig. 17c), and over the Arctic Ocean the near-surface temperature change is not large.
In order to determine what processes cause the temperature increase in the NH polar region, we compared the seasonal mean near-surface temperature changes with the decomposed surface anthropogenic aerosol forcing (see previous section). The interaction between anthropogenic aerosols and clouds (ACI) has a strong warming effect at the surface (Fig. 17g, h). Furthermore, both the spatial pattern and the seasonal variations agree well with temperature increase. The strong ACI forcing is mainly contributed by the large increase in liquid water path (not shown). The thicker liquid clouds trap more outgoing longwave radiation and reflect it back to the surface ( FLDS in Fig. 17j, k), which causes the strong temperature increase in the boreal winter season. The results are consistent with previous observational studies (Garrett and Zhao, 2006;Lubin and Vogelmann, 2006), which reported a significant anthropogenic aerosol effect in the longwave over the Arctic through aerosol-cloud interactions.
Results from the AMIP simulations (Fig. 18) show a similar relationship between the changes in near-surface temperature and in the downward longwave radiation due to the anthropogenic aerosol effects. Note that we do not have the decomposed radiative forcing diagnostics output for the AMIP simulation, so here only the relationship between TREFHT and FLDS is shown. There is a strong correlation between TREFHT and FLDS time series. The temporal variations in TREFHT and FLDS are mainly affected by the interannual variations in both the liquid cloud distributions and the amount of anthropogenic aerosol. The nudged sim-ulations (shown as dots) capture the changes in both variables well. This suggests that results from the relatively short nudged simulations can provide a reasonable estimate of the simulated near-surface temperature responses due to the anthropogenic aerosol forcing through the fast processes.

Conclusions
This study investigates the historical and present-day ERF aer in E3SMv1 using AMIP simulations for 1870-2014 and short nudged simulations. We analyze the relationship between key aerosol and cloud quantities, decompose the anthropogenic aerosol forcing in different ways, and evaluate the parameterization sensitivities.
Despite the overall increase in the global mean anthropogenic aerosol burden and optical depth ( τ a ) during 1870-2014, there are large differences in the historical changes in these quantities between different regions and aerosol species. For example, the anthropogenic sulfate burden and sulfate τ a over the NH polar regions and the NH midlatitude regions decrease significantly after 1970, while for carbonaceous aerosols, there is a continued overall increase in burden and τ a . Both the shortwave and longwave ERF aer estimated by E3SMv1 are outside the 2 standard deviation range of the CMIP6 RFMIP model estimates , although they are not the strongest among the CMIP6 RFMIP model results. Similar to the changes in anthropogenic aerosol burden and optical depth, we see large regional differences in the temporal evolution of ERF aer . The overall increasing or decreasing trend in ERF aer resembles the changes in sulfate burden and τ a , especially for NH midlatitude and polar regions.
The relative changes (dlnY versus dlnX) in historical global annual mean anthropogenic aerosol optical depths, CCN concentrations, and cloud droplet number concentrations show overall linear correlation. However, after around 1970, the correlations show a significant change. This is mainly caused by (1) the regional differences in the historical changes in anthropogenic aerosol burden and optical depths, as well as their impacts on the CCN and cloud droplet formation; (2) a shift from sulfate to carbonaceous aerosols in the anthropogenic aerosol composition; and (3) the "saturation effect", as the slope ln N d / ln CCN slightly decreases from 1900 to 2010 for all the inspected regions. For ice cloud properties, we find similar but weaker correlations compared to those seen in liquid cloud properties.
A forcing decomposition analysis following Ghan (2013) shows that the top-of-model ERF aer in EAMv1 is dominated by the anthropogenic aerosol effect through the interactions with clouds. There is a strong compensation between the forcings in the shortwave and longwave components. At the surface, the anthropogenic aerosol effect through the interactions with radiation contributes to more than 40 % of the total forcing. Large surface longwave ACI forcing ap-pears in the NH high-latitude and polar regions, which is caused by the increase in low liquid cloud induced by anthropogenic aerosols. The forcing analysis decomposed by aerosol species shows that anthropogenic sulfate aerosols have the largest effect (both in the shortwave and longwave) on the overall ERF aer . Anthropogenic BC aerosols exert a strong positive shortwave forcing due to absorption but a weaker negative longwave forcing by weakening vertical motions and reducing high cloud amounts. Both anthropogenic POM and SOA have significant negative forcing, but the forcing over the NH high-latitude and polar region is small.
We have also evaluated the sensitivity of simulated ERF aer to parameterization changes in EAM. Out of the four changes that were adopted in EAMv1, resetting the autoconversion tuning parameters (to the original values used in KK2000) has the largest impact on the net forcing (more negative), confirming that the autoconversion tuning parameters applied in EAMv1 indeed help to reduce the effective aerosol forcing. Two changes are made for mixed-phase clouds, i.e., the tuning of the Bergeron process rate and the change in ice nucleation scheme. Both of them significantly reduce the ERF aer in the shortwave and longwave components, but the impacts on the net forcing are relatively small. The last change, in homogeneous ice nucleation, makes the ice formation rate smaller in cirrus clouds. This also significantly reduces the aerosol forcing in both shortwave and longwave components, but the impact on net forcing is negligible.
The results summarized above suggest that in order to constrain the net ERF aer in future versions of E3SM, it is important to constrain the relative changes ( lnX = X/X) in cloud droplet number concentrations (N d ) and liquid water path (LWP). In E3SMv1, very small N d (e.g., N d < 10 cm −3 ) appears frequently under pristine conditions (to be presented in detail in a separate study). Applying a lower bound of cloud droplet number concentration can effectively increase N d under pristine conditions and reduce ln N d in some models, though this treatment is artificial and thus not desirable . One possible solution to reduce ln N d in a more physical way is to improve the representation of natural aerosols (especially for sulfate oxidized from DMS) and the associated cloud droplet formation or depletion processes in the model. Previous modeling studies have shown that the simulated aerosol radiative forcing is sensitive to DMS emission (Carslaw et al., 2013) and oxidation (Fung et al., 2022) treatments. The other way is to improve the droplet formation parameterization so that it can include some important processes (e.g., enhanced mixing induced by cloud-top radiative cooling, droplet spectral dispersion). These processes can usually only be resolved in model simulations at very high resolution (e.g., Zhang et al., 2022). How to parameterize them accurately in global models will be a challenge, and substantial changes might be needed for models like E3SM. Finally, as suggested by the sensitivity simulations, increasing the primary and/or secondary ice for- Figure 18. Annual (a) and boreal winter (b) mean changes in near-surface temperature at 2 m ( TREFHT, K, in blue) and downward longwave radiation flux ( FLDS, W m −2 , in red) over the NH polar region due to the anthropogenic aerosol effects. Both the AMIP (lines) and nudged (dots and circles) simulations are shown here. See Sect. 7 for details. mation in mixed-phase clouds could also indirectly reduce ln N d and lnLWP. There are still some missing components or processes in our model. For example, the atmospheric chemistry is underrepresented in E3SMv1, wherein we only consider simple sulfur chemistry and the oxidants are prescribed. Related to this simplified treatment, nitrate aerosols are currently not treated in E3SMv1. Previous studies (e.g., Bellouin et al., 2011) have shown that nitrate aerosols have partly replaced sulfate in some regions (e.g., Europe). The overall effective aerosol forcing could maintain negative values in the future, even through some fossil fuel emissions are expected to continuously decrease. The representation of nitrate aerosols has been developed for a newer version of E3SM (Wu et al., 2022), with similar treatment used in CAM5 (Zaveri et al., 2021) and CAM6 (Lu et al., 2021). It would be interesting to investigate the impact of nitrate aerosols on the historical and future ERF aer changes using the new model.
In response to the aerosol forcing associated with fast processes (in contrast to the slow processes involving ocean and/or sea ice changes), significant regional changes are seen in the simulated near-surface air temperature. Due to the strong shortwave anthropogenic aerosol effect through ARI and ACI, surface cooling occurs in the midlatitude and lowlatitude regions. Strong surface warming is seen in the NH high-latitude regions during boreal winter, which is mainly caused by denser liquid clouds in these regions that have a large impact on the longwave radiation. Since the warming caused by greenhouse gases is stronger over the polar regions than in lower latitudes, the high-latitude warming as well as midlatitude and low-latitude cooling caused by anthropogenic aerosols further amplify the meridional gradient in the near-surface temperature changes. The small change (close to zero) in the global mean surface temperature in the AMIP and nudged simulations suggests that the surface cooling in the second half of the 20th century simulated by the coupled model is very likely caused by the anthropogenic aerosol effects through slow processes involving ocean-sea ice responses. In the future, we plan to use both the AMIP and coupled simulations to further investigate the anthropogenic aerosol effects on the historical surface temperature change through different processes.