Climate-driven chemistry and aerosol feedbacks in CMIP6 Earth system models

Feedbacks play a fundamental role in determining the magnitude of the response of the climate system to external forcing, such as from anthropogenic emissions. The latest generation of Earth system models includes aerosol and chemistry components that interact with each other and with the biosphere. These interactions introduce a complex web of feedbacks that is important to understand and quantify.

This paper addresses multiple pathways for aerosol and chemical feedbacks in Earth system models. These focus on changes in natural emissions (dust, sea salt, dimethyl sulfide, biogenic volatile organic compounds (BVOCs) and lightning) and changes in reaction rates for methane and ozone chemistry. The feedback terms are then given by the sensitivity of a pathway to climate change multiplied by the radiative effect of the change.

We find that the overall climate feedback through chemistry and aerosols is negative in the sixth Coupled Model Intercomparison Project (CMIP6) Earth system models due to increased negative forcing from aerosols in a climate with warmer surface temperatures following a quadrupling of CO₂ concentrations. This is principally due to increased emissions of sea salt and BVOCs which are sensitive to climate change and cause strong negative radiative forcings. Increased chemical loss of ozone and methane also contributes to a negative feedback. However, overall methane lifetime is expected to increase in a warmer climate due to increased BVOCs. Increased emissions of methane from wetlands would also offset some of the negative feedbacks. The CMIP6 experimental design did not allow the methane lifetime or methane emission changes to affect climate, so we found a robust negative contribution from interactive aerosols and chemistry to climate sensitivity in CMIP6 Earth system models.

1 Introduction

Climate feedback quantifies the change in the Earth's radiation budget as the surface temperature varies. Overall, this feedback must be negative for a stable climate; i.e. the net radiation budget must decrease as surface temperature increases. The dominant negative feedback comes from increased longwave emissions from a warmer surface (Planck response). Warmer surface temperatures lead to changes in the physical climate system (water vapour, lapse rate, surface albedo, clouds) that further modify the radiation budget, contributing additional positive and negative feedbacks (Sherwood et al., 2020). Earth system models (ESMs) extend the complexity of physical climate models by coupling land and ocean biospheres, atmospheric chemistry and aerosols to the physical climate. Within these models, natural processes, chemical reactions and biological transformations respond to changes in climate, and these processes in turn affect the climate. Therefore, the physical climate system and the biogeochemical cycles are coupled, leading to climate feedbacks that may act to further amplify or dampen the climate response to a climate forcing (Arneth et al., 2010; Ciais et al., 2013; Heinze et al., 2019). The importance of biogeochemical feedbacks has long been recognised for the longer timescales involved in palaeoclimate studies, but the realisation of their relevance in the context of anthropogenic climate change is more recent. A multitude of biogeochemical feedbacks have been identified, but the evaluation of their importance for future climate change remains very limited. A recent review of Earth system feedbacks (Heinze et al., 2019) examined the extensive range of feedbacks possible in an Earth system framework. The largest biogeochemical feedback contribution comes from the carbon cycle (Friedlingstein, 2015). Arneth et al. (2010) considered a range of terrestrial biogeochemical feedbacks interacting with the carbon cycle. O'Connor et al. (2010) reviewed potential feedbacks involving methane. Carslaw et al. (2010) reviewed climate feedbacks involving natural and anthropogenic aerosols. Climate change can impact both the source strength of natural aerosols such as sea salt, dust, biomass burning aerosols or their precursors (dimethyl sulfide (DMS), biogenic volatile organic compounds) and the lifetime of natural and anthropogenic aerosols through changes in transport and dry and wet deposition (Bellouin et al., 2011; Raes et al., 2010). Here, we choose to focus especially on those feedbacks that are mediated through changes in the abundances of reactive gases and aerosols, using data from CMIP6 (Coupled Model Intercomparison Project phase 6) (Eyring et al., 2016) Earth system models that conducted the AerChemMIP (Aerosols and Chemistry Model Intercomparison Project) simulations (Collins et al., 2017).

Note that in this paper we use change in global mean surface temperature as our measure of climate change and for simplicity assume changes in other climate variables are proportional to this. For many of the forcing agents considered here, the forcing pattern varies strongly on regional scales and would be expected to cause larger regional temperature changes than represented by the global mean.

In Sect. 2, we describe the theoretical framework used to diagnose the feedbacks. In Sect. 3, we describe how the different Earth system models implement the biogeochemical processes. Section 4 quantifies the feedbacks as implemented in the models and compares these results with previous modelling and theoretical studies. Section 5 concludes the paper. The Supplement contains further details of the models used and additional figures to support the analysis in Sect. 4.

2 Theoretical framework to analyse biogeochemical feedbacks

2.1 Theory

In order to compare climate feedbacks, we need to compare them on a common scale of the change in the top-of-atmosphere radiation balance following a unit warming (in W m⁻² K⁻¹) (e.g. Gregory et al., 2009). Following Gregory et al. (2004), the radiative imbalance ΔN from an imposed forcing ΔF is given by

$$\mathrm{\Delta }N=\mathrm{\Delta }F+\mathit{\alpha }\mathrm{\Delta }T,$$

where ΔT is the global mean change in surface temperature and α is the climate feedback parameter ($=\frac{\mathrm{d}\mathrm{\Delta }N}{\mathrm{d}\mathrm{\Delta }T}$). The total derivative $\frac{\mathrm{d}\mathrm{\Delta }N}{\mathrm{d}\mathrm{\Delta }T}$ can be split into a set of partial derivatives:

$${\displaystyle \frac{\mathrm{d}\mathrm{\Delta }N}{\mathrm{d}\mathrm{\Delta }T}}={\sum }_i{\displaystyle \frac{\partial \mathrm{\Delta }N}{\partial \mathrm{\Delta }{C}_i}}{\displaystyle \frac{\partial \mathrm{\Delta }{C}_i}{\partial \mathrm{\Delta }T}}={\sum }_i{\mathit{\alpha }}_i,$$

where the α_i indicates the individual feedback terms due to a change in a climate variable C_i. For feedbacks involving changes in composition, the ΔC_i can represent changes in reactive gas or aerosol burdens or emissions. ${\mathit{\alpha }}_i=\frac{\partial \mathrm{\Delta }N}{\partial \mathrm{\Delta }{C}_i}\frac{\partial \mathrm{\Delta }{C}_i}{\partial \mathrm{\Delta }T}$ can then be expressed as ϕ_iγ_i, where ϕ_i is the radiative efficiency of the species per burden (W m⁻² Tg⁻¹) or per emission (W m⁻² (Tg yr⁻¹)⁻¹), and γ_i is the change in species burden or emission with climate (Tg K⁻¹ or Tg yr⁻¹ K⁻¹). The radiative efficiencies are based on effective radiative forcing (ERF) (Myhre et al., 2013a) to include rapid adjustments to changes in composition. Since climate change can also affect the atmospheric lifetime of a species, $\frac{\partial \mathrm{\Delta }{\text{Burden}}_i}{\partial \mathrm{\Delta }T}$ does not necessarily scale with $\frac{\partial \mathrm{\Delta }{\text{Emission}}_i}{\partial \mathrm{\Delta }T}$.

2.2 Applying the theory to Earth system models

With Earth system models, the ϕ_i and γ_i coefficients can be diagnosed from idealised simulations in which only climate or composition are changed. Here, we use the set of simulations specified under the CMIP6 project (Eyring et al., 2016).

The γ_i values are diagnosed from a pair of idealised climate change scenarios: a control climate (piControl) where composition is maintained at a level representative of 1850 conditions, and a warmer climate (abrupt-4xCO2) where temperatures have increased following an abrupt quadrupling of CO₂. To quantify the sensitivities to this temperature change, we take the 30-year time means from years 121–150 of these simulations for both the surface temperature change and the burden or emission changes. The global mean surface temperature changes are therefore not the same as the model-based equilibrium climate sensitivity (ECS) calculations but are temperatures consistent with the averaging period for the burden or emissions. The γ_i are calculated from the change in emission or burden divided by the temperature change. For the dust and sea salt (these are the aerosols with single sources), rather than the burden, we diagnose the aerosol optical depth (AOD) change (K⁻¹), where available, as being the quantity most closely related to the radiative forcing (Myhre et al., 2013b). For DMS and organic aerosol emissions, we use the emission change (Tg yr⁻¹ K⁻¹), as changes in aerosol lifetime will also affect AODs from other sources of sulfate and organic aerosol (OA) that we do not have ERF calculations for. For reactive gases, both emissions-based and concentration-based calculations are used. CO₂ can have climate effects beyond its global warming; for instance, CO₂ directly cools the stratosphere and can affect vegetation with implications for dust and biogenic volatile organic compound (BVOC) emissions. With the AerChemMIP setup, it is not possible to distinguish these adjustments to CO₂ concentration from the impacts of surface temperature increase.

The ϕ_i coefficients for changes in emissions are derived from pairs of the AerChemMIP simulations defined in Collins et al. (2017): piClim-control where composition and climate are maintained at a level representative of 1850 conditions, and experiments piClim-2x (Table 1) in which individual natural emission fluxes are doubled. The climate change in these simulations is restricted by using fixed sea surface temperatures and sea ice cover (Collins et al., 2017) for a 30-year mean of the piControl simulation. The ERFs are determined by the mean difference in top-of-atmosphere radiative fluxes between the piClim-2x and the piClim-control experiments over a 30-year period. The ϕ_i values are calculated from the ERF divided by either the change in AOD or change in emissions, depending on the units of γ_i above. The specific simulation variant numbers are listed in Table S2 in the Supplement.

Table 1 List of simulations for diagnosing ERFs of natural emitted species. The specified natural emission fluxes are doubled compared to the 1850 control.

The theoretical framework in Sect. 2.1 is inherently linear, whereas the Earth system may well not be. The climate changes used to diagnose γ_i are of the order 4–7 K (Table 5), which are much larger than the remaining ∼0.5–1 K goals of the Paris Agreement. The doubled natural emission changes used to diagnose ϕ_i are larger than the changes found in the 4xCO2 experiments and larger still than those expected from a climate following the Paris goals.

For ${\mathit{\varphi }}_{{\mathrm{O}}_{\mathrm{3}}}$, the ozone radiative forcing (Tables 10 and 11) is diagnosed from the changes in the 3-D ozone distributions multiplied by a 3-D kernel of ozone radiative efficiencies from Skeie et al. (2020). The uncertainty in radiative transfer modelling was estimated to be only 10 % in Stevenson et al. (2013), but we increase that to 15 % as a conservative estimate comparable to the 14 % radiative modelling uncertainty for methane (Etminan et al., 2016). Radiative modelling uncertainties are negligible compared to the other uncertainties in Sect. 4.

The ESM setups here, even with tropospheric chemistry, still constrain methane to specified concentrations at the surface. This means that any feedbacks mediated through changes in oxidising capacity have a negligible effect on methane. It is, however, possible to diagnose the change in methane that would be expected, if it were not constrained, from the change in its lifetime:

$${\displaystyle \frac{\mathrm{\Delta }C}{C}}={\left({\displaystyle \frac{\mathrm{\Delta }\mathit{\tau }}{\mathit{\tau }}}+\mathrm{1}\right)}^{\mathrm{f}}-\mathrm{1}\approx f{\displaystyle \frac{\mathrm{\Delta }\mathit{\tau }}{\mathit{\tau }}},$$

where C is the methane concentration, τ is the total methane lifetime (including loss to soils), and f is the feedback of methane on its own lifetime (Fiore et al., 2009). The effective radiative forcing from the change in concentration is $\mathrm{7.0}\times {\mathrm{10}}^{-\mathrm{4}}$ W m⁻² ppb⁻¹, calculated using the formula from Etminan et al. (2016) from a methane baseline of 802 ppb representative of 1850 (Myhre et al., 2013a); this is scaled by 1.52 to account for the additional chemical production of ozone (0.4) and stratospheric water vapour (0.12). These values are reduced from 0.5 and 0.15 in Myhre et al. (2013a) (their Sect. 8.SM.11.3.2) as the 25 % increase in radiative efficiency from Etminan et al. (2016) does not affect the ozone or water vapour. This gives 1.11 W m⁻² per fractional change in methane lifetime or 0.011 W m⁻² %⁻¹. Changes in methane concentration due to changes in emissions ΔE are given by

$$\mathrm{\Delta }C=\mathrm{\Delta }E\mathit{\tau }f\left({\displaystyle \frac{m_{\mathrm{air}}}{m_{{\mathrm{CH}}_{\mathrm{4}}}}}\right)/M_{\mathrm{atm}},$$

where τ=9.1 years (Prather et al., 2012), and f=1.34 (Myhre et al., 2013a). m_air and $m_{{\mathrm{CH}}_{\mathrm{4}}}$ are the relative molecular masses of air and methane (28.97 and 16.0).

3 Model descriptions

3.1 Model implementation of aerosols, tropospheric and stratospheric chemistry

We use results from seven Earth system models that contributed simulations under the AerChemMIP piClim-2x experimental setup. All seven models have interactive aerosol schemes and five have interactive stratospheric chemistry of which four also have interactive tropospheric chemistry (Table 2). The level of sophistication of the chemistry can affect the modelled responses to the emissions of reactive gases. For instance, in models without interactive tropospheric chemistry, changes in BVOCs affect only organic aerosols, whereas in models with interactive tropospheric chemistry, they also affect ozone, methane lifetime and potentially the oxidation of other aerosol precursors. For each model, one ensemble member was run for each experiment.

Table 2 Sophistication of gas-phase chemistry used in the Earth system models (for further details, see Thornhill et al., 2020).

3.2 Model implementation of natural emissions of aerosols and ozone precursors

3.2.1 Land

The land-based natural emissions analysed here are dust, BVOCs and wetland methane (Table 3).

Table 3 Levels of complexity of vegetation included in the land-based emissions schemes of dust and BVOCs for the ESMs, including dependence on photosynthetically active radiation (PAR) and leaf area index (LAI).

“N/A” signifies that the given diagnostic was not available from that model.

Dust emissions are parameterised as a function of surface wind speeds or wind stress, and account for the amount of bare soil, soil type and aridity (Ackerley et al., 2012; Collins et al., 2011; Evan et al., 2014; Fiedler et al., 2016; Huneeus et al., 2011; Shao et al., 2011; Zender et al., 2004). There is a variation between the models in the sizes considered, whether binned or modal, and the optical properties of the dust particles (Kok et al., 2018; Xie et al., 2018). Table S1 lists the parameterisations for desert-dust aerosol for the contributing models and the simulated dust–aerosol sizes.

BVOC emissions are parameterised as a function of vegetation type and cover, and also temperature and photosynthesis rates (gross primary productivity) (Guenther et al., 1995; Pacifico et al., 2011; Sporre et al., 2019; Unger, 2014). Some parameterisations also include dependence on CO₂ concentrations (Pacifico et al., 2012). Models differ in the speciation of the VOCs emitted but typically include isoprene and monoterpenes, with different emission parameterisations for different species. The ability of VOCs to form secondary organic aerosol is typically parameterised as a fixed yield (Mulcahy et al., 2020). For further details, see Table S1 in the Supplement and references therein.

3.2.2 Marine

The ocean emissions analysed here are sea salt, DMS and primary organic aerosols (Table 4).

Table 4 Levels of complexity of marine emissions in the ESMs.

The air–sea exchange processes for these emissions are parameterised as a function of wind speed and sometimes temperature (Gong, 2003; Jaeglé et al., 2011).

Changes in DMS emissions can be initiated by various factors such as changes in temperature, insolation, depth of the ocean-mixed layer, sea ice extent, wind strength, nutrient recycling or shift in marine ecosystems (Heinze et al., 2019). The DMS fluxes into the atmosphere are prescribed in some models (CNRM-ESM2-1, GFDL-ESM4, MIROC6, CESM2-WACCM) and calculated interactively from ocean biogeochemistry in others (UKESM1, NorESM2). Further details on the current generation of marine biogeochemical models, including the representation of DMS emission scheme, can be found in Séférian et al. (2020). Oceanic organic aerosol emissions are also wind-speed dependent and in addition depend on chlorophyll concentrations generated either from interactive biogeochemistry or observation-based chlorophyll concentrations in models without ocean biogeochemistry components.

3.2.3 Lightning

The models with tropospheric chemistry (UKESM1, GFDL-ESM4, CESM2-WACCM, GISS-E2-1) all include parameterisations of the emission of nitrogen oxides (NO_x) from lightning, related to the height of the convective cloud top (Price et al., 1997; Price and Rind, 1992). The lightning frequency depends strongly on the convective cloud-top height, and the ratio of cloud-to-cloud vs. cloud-to-ground lightning depends on the cold cloud thickness (from 0 ^∘C to the cloud top). The precise implementation of lighting emissions and their height profile varies between the models.

4 Quantification of feedbacks

The feedbacks in this section are all derived from the difference between the piControl and abrupt-4xCO2 CMIP6 experiments. The Earth system models all respond with different levels of climate change, so all climate feedbacks are normalised to the change in global mean surface temperature between abrupt-4xCO2 and piControl for the 30-year period (years 121–150; Table 5) to derive the γ_i (Sect. 2.1). There is a factor of nearly 2 between the temperature responses of the models. Since this time frame is not long enough for the models to have reached equilibrium (which may take many centuries), these temperatures are not the same as ECS.

Table 5 Change in global mean surface temperature following an abrupt quadrupling of CO₂ concentrations. The difference between abrupt-4xCO2 and piControl is averaged over the years 121–150. Uncertainties refer to the standard deviation of the interannual variability.

4.1 Aerosol species

4.1.1 Desert dust

The 2xdust perturbation is applied by scaling the parameterisation in the emission scheme. Since changing dust emissions will affect the boundary layer meteorology, the net effect is not an exact doubling of the emissions (Table 6). Four of the six models in AerChemMIP have a negative radiative forcing for doubled dust (Figs. 1a, S2–S4, Table 6). The models all agree on a negative ERF over the oceans close to the source regions. They differ in the sign of the ERF over the deserts themselves, with most (four out of six) showing a positive longwave ERF (Fig. S4). The shortwave ERF is more variable (Fig. S3) and is also affected by any changes in low cloud amount. For CNRM-ESM2-1 and UKESM1, this positive ERF over the deserts outweighs the oceanic negative ERF. The ERF for GFDL-ESM4 is not significantly different from zero. UKESM1 has by far the largest dust emissions (and change from doubling) because it includes particles that are emitted and deposited in the same time step. CNRM-ESM2-1 also includes large particles (up to 20 µm). These models, however, have similar changes in dust AOD compared to the other models, and hence the magnitude of the forcing efficiency per change in AOD (Table 6) is not out of line with the others. MIROC6 has the strongest forcing even with the lowest emissions and smallest change in AOD, thus giving it the largest forcing efficiency per AOD.

Table 6 Dust radiative efficiencies by emission and AOD from 2xdust experiments. Changes in emissions and AOD are from abrupt-4xCO2. “Scaled” refers to scaling the 2xdust relations between AOD and emissions by the 4xCO2 changes in emissions. The α values are calculated assuming ERF is proportional to emissions or AOD. Uncertainties for each model are errors in the mean based on interannual variability. Uncertainties in the multi-model results are standard deviation across the models. The multi-model α terms are the average of the individual model α rather than the product of the multi-model ϕ and γ. Multi-model means are not shown for the emissions as some models include coarse particles, whereas others do not.

“N/A” signifies that the given diagnostic was not available from that model.

Figure 1 Multi-model mean (a) ERF from piClim-2xdust vs. piClim-control, (b) change in dust emissions for abrupt-4xCO2 vs. piControl. Stippling shows areas where the mean changes by more than the standard deviation across models.

The response of dust aerosols to abrupt-4xCO2 (Figs. 1b, S1) is substantially different across the model ensemble. Four models (CNRM-ESM2-1, MIROC6, GFDL-ESM4 and GISS-E2) show an increase in dust emission in a 4xCO2 climate due to increased aridity and near-surface wind speeds, whereas UKESM1 has a decrease in dust emissions with more CO₂ due to increased fertilisation of the vegetation (hence less bare soil) paired with decreased near-surface winds. NorESM2 shows near-zero change. The spatial pattern of the opposing response of dust emission to 4xCO2 in the two most extreme models, UKESM1 and CNRM-ESM2-1, is consistent with the responses in near-surface wind speed to 4xCO2 (Fig. S5). These reflect larger increases in mean winds over regions where the mean emission amount is larger for 4xCO2 compared to the pre-industrial climatology. The increase or decrease in winds is also likely to be affected by changes in vegetation in semi-arid regions, e.g., the Sahel.

As well as affecting the emissions, changing climate can also affect the removal of dust through changes in both dry and wet deposition. In all models except UKESM1, the lifetime of dust increases (Table 6). The effect of an increase in lifetime can be seen by comparing the change in AOD. The modelled changes in dust AOD in the abrupt-4xCO2 experiment are a factor of 1.5–2 larger (for those models where lifetime increases) as would be expected assuming a linear scaling with emissions across all size ranges (“scaled AOD” in Table 6).

The climate feedback parameter for dust (α) is given by the product of the radiative efficiencies (ϕ) with the sensitivities to climate (γ). These vary from −0.012 to +0.0020 W m⁻² K⁻¹ with a multi-model mean of $-\mathrm{0.0026}\pm \mathrm{0.0048}$ W m⁻² K⁻¹, i.e. consistent with zero. Scaling with AOD change rather than emission change gives a slightly larger magnitude, with a range of −0.016 to +0.0048 W m⁻² K⁻¹ and a multi-model mean of $-\mathrm{0.0040}\pm \mathrm{0.0072}$ W m⁻² K⁻¹. Although some models obtain similar feedback terms, this is not necessarily for the same reason. For instance, GFDL-ESM4 and NorESM2 have small feedback terms. NorESM2-LM has a large ERF for doubled dust emissions but a small change in dust emission for 4xCO2, whereas GFDL-ESM4 has a large change in emissions but a small ERF.

Dust–aerosol feedback assessments are a relatively new area of research due to the large uncertainties of climate models in simulating dust aerosols with changes in atmospheric composition. For instance, the spread in model estimates for dust aerosol changes in the 21st century is the largest among wildfires, biogenic SOA and DMS sulfate (Carslaw et al., 2010). Predictions for future dust emission range from an increase (Woodward et al., 2005) to a decrease (Mahowald and Luo, 2003). The modelled feedbacks in Table 6 are smaller in magnitude compared to the theoretical model estimates of −0.04 to +0.02 W m⁻² K⁻¹ by Kok et al. (2018).

The model ranges in dust forcing and feedbacks are not surprising in light of past studies that highlight model differences in dust-emitting winds and dust–aerosol parameterisations that contribute to the model diversity in the dust–aerosol loading, optical properties and radiative effects (Ackerley et al., 2012; Evan et al., 2014; Huneeus et al., 2011; Shao et al., 2011; Zender et al., 2004). For instance, the parameterisation of the planetary boundary layer plays an important role in determining the dust loading, forcing and regional feedbacks on winds (Alizadeh Choobari et al., 2012). Influencing factors for regional differences in the dust radiative effects are the surface albedo and aerosol size distribution (Kok et al., 2018; Xie et al., 2018), whereas feedbacks on winds depend also on meteorological factors (Heinold et al., 2008). The substantial model differences in the dust emission response to 4xCO2 paired with corresponding differences in mean 10 m wind speed in this study suggest that also the dust feedback parameter critically relies on accurately simulating atmospheric dynamics. Modelling atmospheric circulation has been identified as a grand challenge in climate research (Bony et al., 2015). Currently, we have no estimate which of the dust feedbacks shown are the most plausible, because convective dust storms are missing in such models, but this dust storm type is believed to be important for north African dust emissions (Heinold et al., 2013). Moreover, natural aerosol–climate feedbacks are thought to depend on the anthropogenic aerosol burden and might therefore be both time-dependent and underestimated in the present-day polluted atmosphere (Spracklen and Rap, 2013). Taken together, we have low confidence in the feedback estimates for dust aerosols to increases in atmospheric concentrations of greenhouse gases.

4.1.2 Sea salt

All models show a strong negative forcing to double sea salt emissions (Figs. 2a, S7, Table 7), although the ERF for MIROC6 is considerably smaller than for the others. The emissions and mass loading for the CNRM-ESM2-1 model are approximately 20 times those of the other models, largely due to including a size bin up to 20 µm. This coarse bin contains a large mass but a lower number of particles, so the AOD change is similar to other models. All models show a similar forcing efficiency per AOD change. All models show an increase in sea salt emissions in the Southern Ocean in 4xCO2 (Figs. 2b, S6) due to increased wind speeds, with a general tendency for decreases elsewhere due to rising temperatures (Jaeglé et al., 2011). The global mean change in emissions is positive in all models except MIROC6 and GISS-E2-1 (where the lower-latitude decreases outweigh the high-latitude increases). For models showing an increased sea salt lifetime in a 4xCO2 climate, the modelled increase in AOD is larger than that expected from scaling the emissions change (“scaled AOD” in Table 7). Although emissions (and the mass burdens) of sea salt decrease in MIROC6 and GISS-E2-1, the AODs increase. The mean feedback is $-\mathrm{0.027}\pm \mathrm{0.032}$ W m⁻² K⁻¹ based on emissions and $-\mathrm{0.049}\pm \mathrm{0.050}$ W m⁻² K⁻¹ based on the increase in AOD. The signs are consistently negative, except for the emission-based feedbacks for MIROC6 and GISS-E2-1.

Table 7 Radiative efficiencies by emission and AOD from 2xss (sea salt). Changes in emissions and AOD are from 4xCO2. “Scaled” refers to scaling the 2xss relations between AOD to emissions by the 4xCO2 changes in emissions. α values are calculated assuming ERF is proportional to emissions or AOD. Uncertainties for each model are errors in the mean based on interannual variability. Uncertainties in the multi-model results are the standard deviation across the models. The multi-model α terms are the average of the individual model α rather than the product of the multi-model ϕ and γ. Multi-model means are not shown for the emissions as these are so variable.

“N/A” signifies that the given diagnostic was not available from that model.

Figure 2 Multi-model mean (a) ERF from piClim-2xss vs. piClim-control, (b) change in sea salt emissions for abrupt-4xCO2 vs. piControl. CNRM-ESM2-1 emissions are excluded from the multi-model emissions in panel (b) as they include a coarse bin which dominates. Stippling shows areas where the mean changes by more than the standard deviation across models.

4.1.3 DMS

Four models ran the 2xDMS experiment. Interactive biogeochemistry or interactive DMS emissions are not a prerequisite for the 2xDMS experiment. However, interactive emissions are required to calculate a feedback α; hence, we exclude CNRM-ESM2-1 from Table 8. Two models include interactive ocean biogeochemistry (UKESM1 and NorESM2). The ERF for 2xDMS is negative for all three models that ran this experiment (Figs. 3a, S9, Table 8), though less strongly so for GISS-E2-1. UKESM1 and NorESM2 show a decrease in sulfur emissions in 4xCO2, where the tropical decrease more than compensates for the increase along the edge of the sea ice retreat, whereas GISS-E2-1 shows an increase in overall sulfur emissions. The multi-model mean is shown in Fig. 3b and the individual models in Fig. S8. The multi-model mean emission-based α is slightly positive but consistent with zero. In spite of decreased DMS emissions in UKESM1 and NorESM2, there is an increased sulfur mass in all models in the 4xCO2 simulation due to an increase in the sulfate lifetime of around 2 % K⁻¹. Since this lifetime change applies to all sulfate, not just that from DMS, the radiative efficiency from 2xDMS will not necessarily apply, and we do not calculate an AOD or mass-based feedback, but note that it would be negative.

Table 8 Radiative efficiencies by emissions from 2xDMS. Changes in emissions are from the 4xCO2 experiment. α values are calculated assuming ERF is proportional to emissions. Uncertainties for each model are errors in the mean based on interannual variability. Uncertainties in the multi-model results are the standard deviation across the models. The multi-model α terms are the average of the individual model α rather than the product of the multi-model ϕ and γ.

Figure 3 Multi-model mean (a) ERF from piClim-2xDMS vs. piClim-control, (b) change in DMS emissions (in g(S)) for abrupt-4xCO2 vs. piControl. Stippling shows areas where the mean changes by more than the standard deviation across models.

DMS is produced by marine biological activity in the ocean, and it is assumed to be the largest natural source of sulfur to the atmosphere. Up to now, there has been no comprehensive model effort to include all the important effects, and therefore the DMS emission strength change under climate change is still uncertain. The slightly positive mean here is in contrast to the −0.02 W m⁻² K⁻¹ feedback from AR5 (Ciais et al., 2013) based on results from only one model (HadGEM2-ES).

Modelling studies including ocean biogeochemistry have shown that under climate change, an increased stratification of the ocean at low and midlatitudes leads to a reduction in nutrients supply into the surface ocean and thus a reduction in DMS emissions, whereas at high latitudes, retreat of sea ice can lead to increased biological activity and increase in DMS production (Kloster et al., 2007). Previous models which include ocean biogeochemistry have shown a slight increase in DMS production and emission to the atmosphere in a warming climate (Bopp et al., 2004; Gabric et al., 2004; Gunson et al., 2006; Vallina et al., 2007).

Some more recent studies have included the impact of ocean acidification on ocean DMS production (Schwinger et al., 2017; Six et al., 2013). Both studies used a very similar description of the ocean biogeochemistry and extended it with an observationally based relation between ocean alkalinity and ocean DMS production. Assuming a medium sensitivity of the DMS production on pH, Six et al. (2013) found a global DMS emission decrease by 18 % in 2100 under the Special Report on Emissions Scenarios (SRES) A1B scenario, and Schwinger et al. (2017) found an emission reduction by 31 % in 2200 under the Representative Concentration Pathways (RCP) 8.5 scenario. In addition, recent work has provided evidence for major pathways in the oxidation of DMS in the atmosphere which are not included in any of these ESMs (Berndt et al., 2019; Wu et al., 2015).

4.1.4 Organic aerosol

Biogenic VOC emissions lead to both organic aerosol and ozone production (in those models with tropospheric chemistry). It is therefore necessary to distinguish the two in the ERFs in these models. The ozone stratospheric-temperature-adjusted radiative forcing (SARF) from the ozone changes is diagnosed offline (Sect. 2.1). This is subtracted from the ERF to give the ERF due to aerosols only as shown in Table 9 (ozone is the only non-aerosol forcing agent that varies). For NorESM2, there is no ozone change. The ERF before subtracting the ozone SARF is shown in Fig. 4. These estimated aerosol forcing changes are large (up to −0.69 W m⁻²). All the ERF-SARF${}_{{\mathrm{O}}_{\mathrm{3}}}$ values are negative, apart from UKESM1, which has a large positive forcing from cloud changes (diagnosed from comparing all-sky and clear-sky diagnostics; not shown).

Table 9 Non-O₃ ERF (subtracting the O₃ SARF from Table 10; for NorESM2, there is no O₃ change). Radiative efficiencies by emission of BVOC from 2xVOC. Changes in emissions of BVOC are from the 4xCO2 experiment. α values are calculated assuming ERF is proportional to emissions. Uncertainties for each model are errors in the mean based on interannual variability. Uncertainties in the multi-model results are the standard deviation across the models. The multi-model α terms are the average of the individual model α rather than the product of the multi-model ϕ and γ.

Figure 4 Multi-model mean (a) ERF from piClim-2xVOC vs. piClim-control, (b) change in BVOC emissions for abrupt-4xCO2 vs. piControl. Stippling shows areas where the mean changes by more than the standard deviation across models.

In terms of aerosol, there is an increase in OA mass and expected increase in AOD with a very similar spatial pattern when the emission of BVOCs is doubled. The patterns of BVOC increase for the 4xCO2 experiments are much more similar between models (Fig. S10) in terms of pattern and sign than for the previous species (dust, sea salt, DMS), although the magnitude is considerably lower for UKESM1. In the 4xCO2 experiments, these models also simulate an increase in primary organic aerosol emissions from the ocean which adds to the OA mass on top of the effect of BVOC emissions. The feedback factors are negative, apart from UKESM1, and are very strong in some models (NorESM2 with −0.28 W m⁻² K⁻¹).

4.2 Ozone and methane feedbacks

4.2.1 Biogenic VOCs

The ozone SARF is diagnosed offline (Sect. 2.1) and shown in Table 10. For all except UKESM1, the magnitude of the ozone forcing is smaller than that for aerosols, leading to a net negative ERF from BVOCs. For UKESM1, the non-ozone forcing is positive (Sect. 4.1.4) and the ozone adds to this. The ozone SARF per Tg VOC emission is similar between the models with CESM2-WACCM slightly lower. The overall feedback is therefore dominated by the variation in the sensitivity of BVOC emissions to climate. This ranges from 0.005 W m⁻² K⁻¹ for UKESM1, which has the lowest BVOC increase with climate, to 0.014 W m⁻² K⁻¹ for CESM2-WACCM and GISS-ES-1, which have the strongest BVOC response to climate.

Table 10 Ozone SARF and radiative efficiencies for 2xVOC emissions. Changes in emissions are from the 4xCO2 experiment. α values are calculated using the ozone SARF. Uncertainties for each model are errors in the mean based on interannual variability and assuming a 15 % uncertainty in the ozone radiative efficiency (Sect. 2.2). Uncertainties in the multi-model results are the standard deviation across the models.

Table 11 ERF and ozone SARF radiative efficiencies for 2xNOX lightning NO_x emissions. Changes in emissions are from the 4xCO2 experiment. α values are calculated assuming ERF or ozone SARF. Uncertainties for each model are errors in the mean based on interannual variability and assuming a 15 % uncertainty in the ozone radiative efficiency (Sect. 2.2). Uncertainties in the multi-model results are the standard deviation across the models.

At the multi-model mean level, the cooling associated with an increase in organic aerosol ($-\mathrm{0.04}\pm \mathrm{0.04}$ W m⁻² K⁻¹ – for the four models with chemistry) dominates over the warming associated with an increase in O₃ (0.011±0.004 W m⁻² K⁻¹), leaving an overall negative feedback.

Using multi-annual simulations of global aerosol, Scott et al. (2018) diagnosed a feedback from biogenic secondary organic aerosol of −0.06 W m⁻² K⁻¹ globally and −0.03 W m⁻² K⁻¹ when considering only extratropical regions. This global feedback value was composed of a direct aerosol radiative feedback of −0.048 W m⁻² K⁻¹ and an indirect aerosol (i.e., cloud albedo) feedback of −0.013 W m⁻² K⁻¹. Using observations from 11 sites, Paasonen et al. (2013) estimated an indirect aerosol feedback of −0.01 W m⁻² K⁻¹ due to biogenic secondary organic aerosol. The ability of models to account for changes in vegetation has a large effect on the feedback. Sporre et al. (2019) found that interactive vegetation enhanced BVOC emissions by 63 % relative to prescribed vegetation, producing more organic aerosol and an increase in (negative) aerosol forcing.

The level of compensation between increased aerosol forcing and increased ozone is dependent on the model (here positive feedback for GFDL-ESM4, negative for UKESM1 and CESM2-WACCM). Unger (2014) found a positive feedback in NASA GISS ModelE2, whereas Scott et al. (2014) found a negative feedback in HadGEM2-ES.

4.2.2 Lightning NO_x

Lightning NO_x leads to ozone production and changes in methane lifetime. As for BVOCs (Sect. 4.2.1), ozone radiative kernels are used to quantify the ozone SARF. The ERF and SARF${}_{{\mathrm{O}}_{\mathrm{3}}}$ agree for all models except UKESM1 (Table 11), suggesting that there is little effect on aerosols in these models. In UKESM1, NO_x is known to increase the formation of new sulfate particles (O'Connor et al., 2020), partially offsetting the positive ozone forcing. The SARF${}_{{\mathrm{O}}_{\mathrm{3}}}$ per Tg emission varies by a factor of 2 (0.023 to 0.048 W m² (Tg(N) yr⁻¹)⁻¹) between the highest and lowest values.

The changes in lightning NO_x emissions vary widely across the models, with three showing increases (UKESM1, CESM2-WACCM, GISS-E2-1) but a slight decrease in GFDL-ESM4. Although they all use variations on the cloud-top-height schemes (Sect. 3.2.3), the differences in how this is implemented and how the modelled clouds vary with climate change all affect the emission response. The feedback is positive for the three models with increased lightning (0.009 to 0.016 W m⁻² K⁻¹), based on the ozone changes, but slightly negative for GFDL-ESM4 (−0.001 W m⁻² K⁻¹). Including the aerosol response to lightning for UKESM1 would reduce its feedback to 0.005 W m⁻² K⁻¹, but this seems to be particular to this model.

The ESMs used in CMIP6 all use a cloud-top-height parameterisation of lightning. Such schemes have previously been found to increase lightning production in warmer climates, whereas more sophisticated schemes based on convective updraft mass flux or ice flux show decreases in lightning with temperature (Clark et al., 2017; Finney et al., 2016b, 2018). The result from the Atmospheric Chemistry and Climate Model Intercomparison (ACCMIP) of 0.44 Tg(N) yr⁻¹ K⁻¹ (Finney et al., 2016a) lies within the range of the models with increased lightning under 4xCO2 (0.27 to 0.61 Tg(N) yr⁻¹ K⁻¹).

4.2.3 Methane lifetimes

BVOC and NO_x emissions also affect the methane lifetime. Methane does not change in the AerChemMIP experimental setup, but the methane changes that would be expected if methane were allowed to evolve freely can be diagnosed from the change in methane lifetime. The methane lifetime to OH (troposphere and stratosphere) is diagnosed in the models. The losses to chlorine oxidation and soil uptake are assumed to be 11 and 30 Tg yr⁻¹, respectively (Saunois et al., 2020). All models show an increase in methane lifetime with BVOC emissions (0.018–0.035 % (Tg(VOC)yr⁻¹)⁻¹ and a decrease due to lightning NO_x emissions (−2.4 % to −6.8 % (Tg(VOC)yr⁻¹)⁻¹ (Table 12). From these, the expected lifetime changes with climate can be deduced from the changes in emissions with temperature. These lifetime changes are then converted to feedbacks using the radiative efficiency (including impacts on ozone and stratospheric water vapour) for methane lifetime changes in Sect. 2.2 (0.011 W m⁻² %⁻¹). The feedbacks range from 0.012 to 0.061 W m⁻² K⁻¹ for BVOCs and −0.042 to +0.001 W m⁻² K⁻¹ for lightning NO_x, where the variability is mostly due to the different sensitivities of BVOC or lightning emissions to climate in the models. For BVOC, the methane lifetime feedback is larger than that due to ozone production, thus increasing the overall feedback. For lightning NO_x, the methane lifetime feedback is of opposite sign to that from ozone production, with approximate compensation for UKESM1 and GFDL-ESM4 (net 0.002 and 0.000 W m⁻² K⁻¹, respectively) and an overall negative lightning feedback from CESM2-WACCM and GISS-E2-1 (−0.009 and −0.028 W m⁻² K⁻¹, respectively). For UKESM1, a feedback of −0.004 W m⁻² K⁻¹ could be added to the total lightning feedback to account for the increase in sulfate.

Table 12 Percentage change in methane lifetime for BVOC and lightning NO_x emissions. Estimated change in lifetime following changes in BVOC and NO_x emissions from the 4xCO2 experiment. α values are calculated assuming a radiative efficiency of 0.015 W m⁻² %⁻¹. Uncertainties for each model assume a 14 % uncertainty in the methane radiative efficiency (Etminan et al., 2016). Uncertainties in the multi-model results are the standard deviation across the models.

4.2.4 Wetland emissions

Two models diagnosed changes in wetland emissions due to 4xCO2. Although the wetland emissions do not directly affect methane concentrations in the model, changes in emissions can be converted to concentration changes (Sect. 2.2). UKESM1 and CESM2-WACCM, both of which are models with interactive wetland emissions, show strong responses to climate change (Table 13), leading to a feedback of 0.16±0.03 W m⁻² K⁻¹.

Wetland emissions are more strongly sensitive to CO₂ concentrations than to temperature or precipitation (Melton et al., 2013), so the values presented here are more likely to be “adjustments” to the CO₂ rather than feedbacks and hence could be considered part of the CO₂ ERF. We find emission increases following quadrupled levels of CO₂ of 130 %–160 %. This compares with results from the Wetland and Wetland CH4 Inter-comparison of Models Project (WETCHIMP) of 20 %–160 % following an increase in CO₂ of a factor of 2.8 (Melton et al., 2013). The CMIP6 simulation specifications do not include free-running methane concentrations; therefore, the effects of these increased wetland emissions will not be realised in any of the CMIP6 experiments. Outside CMIP6, ESMs are starting to include free-running methane (Ocko et al., 2018), so for these it will be important to understand the effects of changing CO₂ and meteorology on wetland emissions.

Table 13 Sensitivity of wetland emissions to 4xCO2 in two models. Feedback parameter assuming pre-industrial conditions. Uncertainties for each model assume a 14 % uncertainty in the methane radiative efficiency (Etminan et al., 2016). Uncertainties in the multi-model results are the standard deviation across the models.

4.2.5 Meteorological drivers

As well as through changes in natural emissions, climate change can affect ozone burden and methane lifetime directly, as the production and loss reactions are sensitive to temperature and water vapour (Johnson et al., 2001). Here, we add the expected changes in ozone SARF and methane lifetime due to changes in BVOCs and lightning NO_x from Sect. 4.2.1 and 4.2.2 and compare those to the changes diagnosed from the 4xCO2 experiments (Table 14). Since lightning NO_x and BVOCs are the dominant climate-sensitive emissions of (non-methane) species affecting ozone and methane, the residual is then the direct effect of climate. UKESM1, GFDL-ESM4 and GISS-E2-1 all diagnosed ozone changes for the abrupt-4xCO2 experiment (Fig. S12). All three showed decreased tropospheric ozone and increased stratospheric ozone (apart from the tropical lower stratosphere) in the 4xCO2 climate. The ozone SARF (calculated using radiative kernels) is negative, whereas the expected change from lightning NO_x and BVOCs would be positive; hence, the residual attributed to meteorological changes is negative.

Table 14 Comparison of expected changes in ozone SARF and methane lifetime with that diagnosed from 4xCO2. The residual is given by the difference and is converted to a feedback using radiative efficiencies for methane lifetime.

For UKESM1, GFDL-ESM4 and GISS-E2-1, the meteorological changes decrease methane lifetime, leading to an overall decrease in lifetime for the 4xCO2. In CESM2-WACCM, the meteorological changes increase methane lifetime, adding to the strong increase from BVOC emissions. This is surprising since there is no known mechanism whereby temperature and humidity increases can increase the methane lifetime. This could be due to non-linearity, whereby the effect of increased VOCs on methane lifetime is larger than expected from scaling the 2xVOC experiment.

Combining the results from ozone and methane lifetime changes leads to overall feedbacks from temperature of −0.15, −0.14 and −0.08 for UKESM1, GFDL-ESM4 and GISS-E2-1.

The three models showing decreased methane lifetime are in approximate agreement with ACCMIP, which found a sensitivity of $-\mathrm{3.4}\pm \mathrm{1.4}$ % K⁻¹ (Naik et al., 2013; Voulgarakis et al., 2013). ACCMIP found a variation in sign of the ozone feedback amongst models $-\mathrm{0.024}\pm \mathrm{0.027}$ W m⁻² for a 1850–2000 change in climate. The ACCMIP models generally did not include stratospheric chemistry, so they either explicitly prescribed the cross-tropopause flux of ozone or imposed a climatology of ozone above the tropopause. The four CMIP6 models here all treat the chemistry seamlessly across the troposphere and stratosphere, so the impact of changes in stratosphere–troposphere exchange (STE) of ozone on the tropospheric column is likely to be different from ACCMIP.

Changes in the stratospheric ozone following a quadrupling of CO₂ are driven by cooling temperatures in the stratosphere. This is likely to be due to temperature adjustments to the stratospheric CO₂ concentrations and so part of the ERF for CO₂ rather than a feedback. Feedbacks and adjustments cannot be distinguished with this experimental setup.

4.3 Overall feedback

The multi-model mean feedbacks are summarised in Table 15 and Fig. 5. The totals assume that feedbacks are additive, which is the basis of the framework in Sect. 2.1. The subsets of model used to generate the multi-model means are different for each process, so the total feedback is a mixture of these different subsets. The largest individual feedbacks are due to the generation of aerosols by BVOCs ($-\mathrm{0.090}\pm \mathrm{0.099}$ W m⁻² K⁻¹) and the emission of methane from wetlands (0.16±0.03 W m⁻² K⁻¹). The overall uncertainty is calculated by adding the inter-model uncertainty on each feedback component in quadrature. This is dominated by the uncertainty in the aerosol response to BVOC emissions. Nearly all the feedbacks are negative, mostly because they come from an increase in aerosol emissions with temperature and increased ozone and methane removal with temperature and humidity. For BVOC emissions, the increase in aerosols outweighs the increases in ozone and methane. For lightning NO_x, the decrease in methane lifetime outweighs the ozone increase. For wetland, we have attributed all the methane emission changes to temperature, whereas a significant proportion are likely to be an adjustment to CO₂ concentrations rather than a feedback (Sect. 4.2.3).

Table 15 Feedback parameters of all the aerosol and chemical processes addressed in this study. Uncertainties are the inter-model standard deviations.

Figure 5 Feedback parameters of all the aerosol and chemical processes in Table 15. Multi-model mean and individual models are shown. Uncertainties are the inter-model standard deviations. BVOC and lightning are the sum of aerosol, ozone and methane lifetime effects (points are only shown for models that include all effects). Ozone and CH₄ lifetime are the chemical effects (i.e. excluding BVOC and lightning emissions). Non-CH₄ is the sum and excludes methane lifetime effects and wetland feedback.

There will be additional systematic uncertainties in the overall feedback term. As described above, the use of a CO₂ perturbation to generate the climate change may lead to different feedback sensitivities compared to climate change caused by other forcing agents. There will also be an uncertainty caused by using a pre-industrial baseline atmosphere rather than the present day. We are unable to quantify the likely magnitudes of these systematic uncertainties.

The ESMs that use the abrupt-4xCO2 experiment to quantify the climate sensitivity do not allow methane to vary, so we also quantify the non-methane feedbacks that will be contributing to the diagnosed climate sensitivity in these models. This feedback is significantly negative ($-\mathrm{0.183}\pm \mathrm{0.111}$ W m⁻² K⁻¹), suggesting the climate sensitivity of ESMs might be expected to be lower than for their physical-only counterparts. This analysis (and climate sensitivity in general) is focused on the global mean, but it should be noted that the cooling effects of increased aerosols will be heterogenous and some regions will experience less warming than a global climate sensitivity might suggest.

5 Conclusions

Earth system models include more processes than physical-only climate models. These models will inherently include additional climate feedbacks and so have a different overall climate feedback (and climate sensitivity) to their physical counterparts. In this study, we consider seven Earth system models (CNRM-ESM2-1, UKESM1, MIROC6, NorESM2, GFDL-ESM4, CESM2-WACCM and GISS-E2-1). In total, six of these (CNRM-ESM2-1, UKESM1, MIROC6, NorESM2, GFDL-ESM4 and GISS-E2-1) participated in the aerosol-related feedback experiments and four (UKESM1, GFDL-ESM4, CESM2-WACCM and GISS-E2-1) in the ozone- and methane-related feedback experiments.

We focus in this study on the responses to an abrupt forcing of quadrupled CO₂ concentrations as that is the usual method to diagnose climate feedbacks. By convention, the feedbacks are quantified as a response to temperature (in W m⁻² K⁻¹), but they may not necessarily be applicable to drivers of climate change other than CO₂ as some of the “feedbacks” may instead be adjustments to CO₂ concentrations. It should also be noted that abrupt-4xCO2 feedbacks are based on atmospheric conditions representative of the 1850s and thus may not be applicable to future responses starting from present-day conditions. For many of the forcing agents considered here, the forcing pattern varies strongly on regional scales and would be expected to cause larger regional temperature changes than represented by the global mean. Thus, aerosol-mediated feedbacks may alter the pattern of climate response as well as the magnitude.

Here, we find that the dominant feedbacks are negative; i.e. they act to dampen the response to an imposed forcing. The total feedback, excluding inferred changes in methane, is $-\mathrm{0.183}\pm \mathrm{0.111}$ W m⁻² K⁻¹. The increase in organic aerosols from increased emission of VOCs from vegetation makes the largest contribution to both the magnitude of the feedback and its uncertainty ($-\mathrm{0.09}\pm \mathrm{0.10}$ W m⁻² K⁻¹), with increases in sea salt and DMS emissions also contributing.

Contributions from increases in ozone production from biogenic VOCs and lightning NO_x are offset by decreased tropospheric ozone lifetime in a warmer climate, leading to an overall negative feedback through ozone. Diagnoses of changes in wetland emissions of methane indicate that if ESMs did allow methane to vary interactively the combined aerosol and chemical feedbacks would be substantially less negative and consistent with zero.

The aerosol and chemistry feedbacks listed here contribute up to the order of −0.2 W m⁻² K⁻¹. This is smaller in magnitude than the carbon cycle response to climate (of the order of 0.5 W m⁻² K⁻¹; Ciais et al., 2013) or the physical climate feedbacks (of the order of 1–2 W m⁻² K⁻¹; Sherwood et al., 2020).