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 modified, 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.

1 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; Forster 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 (Textor et al., 2006) and their impacts on clouds (Zhang et al., 2016). As a result, the estimated anthropogenic aerosol effects in global climate models are quite uncertain for both the aerosol–radiation interactions (Schulz et al., 2006; 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⁻². This value is outside the likely range of −0.63 to −1.37 W m⁻² 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 multi-model estimate from Bellouin et al. (2020), who report a 68 % likelihood range of −1.6 to −0.6 W m⁻² and a 90 % likelihood range of −2.0 to −0.4 W m⁻².

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 (Ma et al., 2022).

• 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.

2 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 Rasch et al. (2019), 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.

2.1 Model description

2.1.1 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 (Golaz et al., 2002; 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 large-scale 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).

2.1.2 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, 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⁻¹ (inherited from the CAM5.4) is applied to the characteristic subgrid updraft velocity¹ 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.²

For mixed-phase clouds, E3SMv1 uses a parameterization based on classical nucleation theory (CNT). The deposition, 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<-\mathrm{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<-\mathrm{37}$ ^∘C), the homogeneous freezing of sulfate aerosols and heterogeneous immersion freezing on mineral dust are considered using the Liu and Penner (2005) 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⁻¹. 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 (Golaz et al., 2019). We discuss this parameter sensitivity in Sect. 6.

2.2 Simulations

Ma et al. (2022)Khairoutdinov and Kogan (2000)Meyers et al. (1992)(Wang et al., 2014)

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.

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 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 (Wang et al., 2018), 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., spatial 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 (Lamarque et al., 2010). 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 (interpolated) 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.

Figure 1 Temporal variations of global and annual mean aerosol burden differences between simulations with transient and pre-industrial (1850) aerosol emissions. Solid lines indicate aerosol species that are directly influenced by anthropogenic aerosol emissions (SO₄, BC, POM, SOA), and dotted lines indicate natural aerosols (DUST, SS, MOM) that are mainly driven by surface winds. Dots indicate results from the short nudged simulations. Note that for clarity, the BC Δburden is scaled by a factor of 5.

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.

3 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 (1870–2014). 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 2 Similar to Fig 1, but for anthropogenic aerosol optical depth at 0.55 µm (total and individual components). Solid lines indicate aerosol species that are directly influenced by anthropogenic aerosol emissions (SO₄, BC, POM, SOA), and dotted lines indicate natural aerosols (DUST, SS, MOM) that are mainly driven by surface winds. Dots indicate results from the short nudged simulations. The black vertical bar in panel (a) shows the estimated range of total Δτ_a by Bellouin et al. (2020). Note that the BC values are scaled by a factor 5 for easier comparison. See Sect. 3.1 for details.

3.1 Changes in anthropogenic aerosol burden and optical depth

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).

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. For example, the global biomass burning emission is largest in 1997 (van der Werf et al., 2017), an extreme El Niño year, and it contributes to a large increase in carbonaceous aerosol burden change. The impact of biomass burning emissions is more evident in cleaner areas like SH midlatitude and polar regions (Fig. 1e, f). Similar increases are also present in the NH polar region, where biomass burning emissions increase significantly after 2011 (especially over boreal North America and boreal Asia; see Fig. 9 of van der Werf et al., 2017).

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 evaluated 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 large-scale 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.

Figure 3 Contribution of sulfate, BC, POM, and SOA to anthropogenic aerosol optical depth (excluding dust, sea salt, and marine organic aerosols) in 1970, 2000, and 2010. Numbers in the panel title are the total anthropogenic aerosol optical depth (Δτ_a,ant = Δτ_a,SO4 + Δτ_a,BC + Δτ_a,POM + Δτ_a,SOA). Results are calculated using the AMIP simulations. For example, the contribution from sulfate is calculated as $\mathrm{\Delta }{\mathit{\tau }}_{\mathrm{a},\mathrm{SO}\mathrm{4}}/\mathrm{\Delta }{\mathit{\tau }}_{\mathrm{a},\mathrm{ant}}$. See Sect. 3.1 for details.

Δτ_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 midlatitude 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).

3.2 Historical changes in ERF_aer

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⁻²) first appears around 1980. In contrast, the longwave ΔF is overall positive, with the largest forcing of about 1.0 W m⁻². 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.

Compared to ΔF estimated by CMIP6 RFMIP models (Smith et al., 2020), 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.

Figure 4 Temporal variations of global and annual mean ERF_aer (W m⁻²) derived from the differences between simulations with transient and pre-industrial aerosol emissions. Solid lines indicate net ERF_aer, and dotted lines indicate the shortwave (blue) and longwave (red) components. Dots indicate results from the short nudged simulations. The cross and the thick and thin vertical bars show the mean, standard deviations, and max–min ΔF estimated by CMIP6 RFMIP models (Smith et al., 2020). See Sect. 3.2 for details.

Figure 5 Geographical distributions of annual mean net (left), shortwave (middle), and longwave (right) ERF_aer (ΔF, W m⁻²) at the top of the model (TOM) estimated for different time slices using nudged simulations. Numbers shown on the top right of each panel are global mean values. For a given time slice, the CMIP6 aerosol emissions for the corresponding year are used. Simulations are described in Table 1, and “E2010” represents the CTRL simulation. See Sect. 3.2 for details.

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 affected 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 Annually and seasonally averaged zonal mean top-of-model ΔF (W m⁻²) estimated for different time slices (emission years) indicated by different colors. Solid lines show the shortwave ΔF and dotted lines the longwave ΔF. In JJA, the shortwave ΔF for 1970 (relative to 1850) is about −14 W m⁻² near the North Pole. See Sect. 3.2 for details.

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₂ 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⁻² in 1970 to −4 W m⁻² 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.

4 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 ³ $\mathrm{\Delta }X/X_{\mathrm{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 (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 (Wan et al., 2022).

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, and 0.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.

Figure 7 Relationships between simulated anthropogenic AOD (Δτ_a) and present-day total AOD (τ_a), relative (PD-PI) changes in AOD (Δln τ_a) and in the column-integrated CCN at 0.1 % supersaturation (ΔlnCCN_0.1 %), Δln τ_a and relative changes in the column-integrated N_d (Δln N_d), and Δτ_a and changes in clear-sky shortwave radiation (ΔF_SW,clear). $\mathrm{\Delta }\ln X=\mathrm{\Delta }X/X_{\mathrm{PD}}$. Dark grey dots indicate estimated values from AMIP PD/PI simulations for the years 1870–1969, while red dots are the estimated values for the years 1970–2014. Each dot represents the ensemble-averaged global annual mean field for a certain year. Blue markers show the estimates based on nudged simulations. The orange square shows the nudged simulation results averaged from conditionally sampled high-frequency data following Ghan et al. (2016), wherein CCN (at 0.3 % supersaturation) and N_d are sampled below 1 km and at the top of liquid clouds. Grey areas in panels (a, d) indicate the multi-model estimated range from Bellouin et al. (2020) and in panels (b, c) the multi-model estimated range from Ghan et al. (2016). See Sect. 4.1 for details.

The relationship between Δln τ_a (Δτ_a/τ_a) and ΔlnCCN_0.1% is shown in Fig. 7b. They have a quasi-linear 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 $>-\mathrm{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 parameterized 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).

4.2 Susceptibility of N_d, LWP, and R_e,liq to changes in CCN

Figure 8 As Fig. 7, but for relationships between the relative changes in CCN (ΔlnCCN_0.1 %) and in liquid cloud properties including droplet number (ΔlnN_d), liquid water path (ΔlnLWP), and effective radius of cloud droplets at cloud top (Δln Re_liq). $\mathrm{\Delta }\ln X=\mathrm{\Delta }X/X_{\mathrm{PD}}$. Dark grey dots indicate estimated values from AMIP PD/PI simulations for the years 1870–1969, while red dots are the estimated values for the years 1970–2014. Each dot represents the ensemble-averaged global annual mean field for a certain year. Blue markers show the estimates based on nudged simulations. The orange square shows the nudged simulation results averaged from conditionally sampled high-frequency data following Ghan et al. (2016), wherein CCN (at 0.3 % supersaturation) and N_d are sampled below 1 km and at the top of liquid clouds. The grey area indicates the multi-model estimated range from Ghan et al. (2016). See Sect. 4.2 for details.

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).

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 hygroscopicity 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; Stevens, 2015), 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 (Δln$N_{\mathrm{d}}/\mathrm{\Delta }$lnCCN_0.1 %) calculated using regional mean fields slightly decreases from 1900 to 2010 for all the inspected regions.

Figure 9 Global and regional annual mean relative changes in CCN (ΔlnCCN_0.1 %, a) and N_d (ΔlnN_d, b) due to differences in anthropogenic aerosols between PI (1850) and PD (2010). Panel (c) shows the ratio between ΔlnN_d) and ΔlnCCN_0.1 %. See Sect. 4.2 for details.

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-to-max 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 As Fig. 7, but for relationships between ΔN_d and N_d and between relative changes in N_d and other quantities. Δln τ_liq, Δln f_liq, and ΔlnCF_SW are the relative changes in liquid cloud optical depth, liquid cloud fraction, and shortwave cloud radiative forcing, respectively. $\mathrm{\Delta }\ln X=\mathrm{\Delta }X/X_{\mathrm{PD}}$. Dark grey dots indicate estimated values from AMIP PD/PI simulations for the years 1870–1969, while red dots are the estimated values for the years 1970–2014. Each dot represents the ensemble-averaged global annual mean field for a certain year. Blue markers show the estimates based on nudged simulations. The orange square shows the nudged simulation results averaged from conditionally sampled high-frequency data following Ghan et al. (2016), wherein CCN (at 0.3 % supersaturation) and N_d are sampled below 1 km and at the top of liquid clouds. The grey area indicates the multi-model estimated range from Ghan et al. (2016). See Sect. 4.3 for details.

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.

4.4 Relationships between changes in liquid and ice cloud properties

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. $\mathrm{\Delta }\ln X=\mathrm{\Delta }X/X_{\mathrm{PD}}$. Dark grey dots indicate estimated values from AMIP PD/PI simulations for the years 1870–1969, while red dots are the estimated values for the years 1970–2014. Each dot represents the ensemble-averaged global annual mean field for a certain year. Blue markers show the estimates based on nudged simulations. See Sect. 4.4 for details.

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 column-integrated 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).

Figure 12 Annual mean column-integrated cloud droplet and ice crystal number concentrations for the pre-industrial era (E1850, panels a, b), present day (CTRL for year 2010, panels c, d), and PD–PI differences (e, f) in the nudged simulations. See Sect. 4.4 for details.

The relative change in global annual mean column-integrated 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 Liu et al. (2005) 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 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.

5 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.

5.1 Decomposition of ERF_aer due to the aerosol–radiation interaction, aerosol–cloud interaction, and surface albedo effect

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.

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.

We use the Ghan (2013) method to decompose the total ERF_aer (ΔF) into individual forcings caused by aerosol–radiation interactions (ARI, direct aerosol effect), aerosol–cloud 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): Δ(F − F_clean)

• 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 effect is much smaller than the indirect aerosol effect (Ghan et al., 2012).

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⁻²), 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.

5.2 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 anthropogenic sulfur emission has a large impact on the ΔF (Fig. 15d). The global annual mean net forcing is about −1.66 W m⁻², 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⁻²), 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⁻². 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⁻², 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 pre-industrial (1850) and present-day (2010) emissions.

Figure 15 Geographical distributions of the annual mean net (left column), shortwave (middle column), and longwave (right column) top-of-model (TOM) effective radiative forcing by individual aerosol components (i.e., compositions). For example, “E2010SU” indicates the same simulation as E1850 but with present-day (2010) sulfur emissions. See Sect. 5.2 and Table 1 for details.

Table 2 Global and annual mean net, shortwave (SW), and longwave (LW) ERF_aer (W m⁻²) at the top of the model (TOM) and surface in the nudged simulations. Model configurations are described in Table 1.

Figure 16 Global annual mean anthropogenic aerosol effects on top-of-model (TOM) radiative fluxes (ΔF and ΔCF) and hydrometeor mass and number (relative changes) in the reference and sensitivity simulations. CF_SW and CF_LW are shortwave and longwave cloud radiative forcings. LWP and IWP are the liquid and ice water path in stratiform clouds. N_d and N_i are the column-integrated droplet and ice number concentrations. $\mathrm{\Delta }\ln X=\mathrm{\Delta }X/X_{\mathrm{PD}}$. See Sect. 6 for details.

6 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 Kogan (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⁻²) 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 (Caldwell et al., 2019), 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 (Gettelman et al., 2019).

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⁻² in the reference model).

We note that the sensitivities discussed above are likely dependent on the quantity of simulated background (pre-industrial) 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 Fig. 2 in Rasch et al., 2019).

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⁻²), and downward longwave radiation flux (ΔFLDS, W m⁻²) due to the anthropogenic aerosol effects. See Sect. 7 for details.

7 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 pre-industrial 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 responses 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 simulations (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.

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⁻², 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.

8 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 (Smith et al., 2020), 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 $\mathrm{\Delta }\ln N_{\mathrm{d}}/\mathrm{\Delta }\ln \mathrm{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 appears 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 ($\mathrm{\Delta }\ln X=\mathrm{\Delta }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⁻³) 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 (Hoose et al., 2009). 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 formation 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 low-latitude 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.