Using historical temperature to constrain the climate sensitivity, the transient climate response, and aerosol-induced cooling

The most recent generation of climate models that has informed the Sixth Assessment Report (AR6) of the Intergovernmental Panel on Climate Change (IPCC) is characterized by the presence of several models with larger equilibrium climate sensitivities (ECSs) and transient climate responses (TCRs) than exhibited by the previous generation. Partly as a result, AR6 did not use any direct quantifications of ECSs and TCRs based on the 4×CO2 and 1pctCO2 simulations and relied on other evidence when assessing the Earth's actual ECS and TCR. Here I use historical observed global-mean temperature and simulations produced under the Detection and Attribution Model Intercomparison Project to constrain the ECS, TCR, and historical aerosol-related cooling. I introduce additivity criteria that disqualify 8 of the participating 16 models from consideration in multi-model averaging calculations. Based on the remaining eight models, I obtain an average adjusted ECS of 3.5 ± 0.4 K and a TCR of 1.8 ± 0.3 K (both at 68 % confidence). Both are consistent with the AR6 estimates but with substantially reduced uncertainties. Furthermore, importantly I find that the optimal cooling due to short-lived climate forcers consistent with the observed temperature record should, on average, be about 47 % ± 39 % of what these models simulate in their aerosol-only simulations, yielding a multi-model mean, global-mean, and annual-mean cooling due to near-term climate forcers for 2000–2014, relative to 1850–1899, of 0.24 ± 0.11 K (at 68 % confidence). This is consistent with but at the lower end of the very likely uncertainty range of the IPCC's AR6. There is a correlation between the models' ECSs and their aerosol-related cooling, whereby large-ECS models tend to be associated also with strong aerosol-related cooling. The results imply that a reduction in the aerosol-related cooling, along with a more moderate adjustment of the greenhouse-gas-related warming for most models, would bring the historical global-mean temperature simulated by these models into better agreement with observations.

1 Introduction

The equilibrium climate sensitivity (ECS) is a well-established (Arrhenius, 1896) yet, despite progress, poorly constrained property of the climate system (Knutti et al., 2017; Forster et al., 2021; Smith and Forster, 2021). For a hypothetical doubling of the atmospheric CO₂ content above preindustrial levels, it states the associated surface temperature increase at equilibrium. Similarly, the transient climate response (TCR) measures the warming simulated in simulations with CO₂ increasing at 1 % per annum (a⁻¹) above its preindustrial abundance at the time of CO₂ doubling. For both quantities, disagreement amongst climate models, particularly in the most recent generation (Meehl et al., 2020) and persisting despite ever-improving model physics and resolution, is an impediment to narrowing their associated long-standing uncertainties. The large spreads in ECSs and TCRs characterizing the present generation of climate models are contributing to some substantial inter-model spread in simulated end-of-century warming in future-scenario simulations (Lee et al., 2021). It is therefore desirable to reduce these model disagreements to more confidently project future climate under any climate scenario.

A standard optimal detection framework assumes that (a) changes in an observed or simulated quantity such as temperature are the sum of changes driven by individual climate forcers (referred to as “additivity”) in the presence of climatological noise and (b) model imperfections can be captured by introducing scaling or correction factors to those model responses to external climate forcers. This amounts, in a regression analysis, to optimally reproducing the observed variations and thereby constraining the TCR and ECS (Schurer et al., 2018, and references therein). The presence of climatological noise means that attribution benefits from using as many models and simulations as possible. Approaches of this kind have been used in many studies using the Coupled Model Intercomparison Project Phase 5 (CMIP5) and earlier model generations (e.g. Hasselmann, 1993; Hegerl et al., 1997; Allen and Tett, 1999; Schurer et al., 2018) and again for the most recent generation, CMIP6 (Gillett et al., 2021). When separating the greenhouse gas (GHG) and aerosol influences, the main problem is that GHGs and aerosols cause warming and cooling effects on climate that are similar in terms of overall temporal developments and therefore can be difficult to distinguish statistically in observational or simulated records of temperature (see below). Furthermore, especially before CMIP6, single-model ensembles used to be small (the largest single-model ensemble used by Schurer et al., 2018, had six ensemble members). These two problems, along with general model disagreements regarding the sizes of these effects, combine to yield substantial uncertainties characterizing successive quantifications of TCR, ECS, and aerosol-induced cooling, such as 1.2 to 1.9 K for the TCR and −0.7 to −0.1 K for the aerosol-induced cooling in CMIP6 (Gillett et al., 2021). Both uncertainties contribute to uncertain projections of future global warming, e.g. 2.1 to 3.5 K of warming in 2081–2100 relative to 1850–1900 in the middle-of-the-road Shared Socioeconomic Pathway (SSP) 2-4.5 (Lee et al., 2021).

In this work, I explore what “historical” all-forcing experiments and single-forcing experiments conducted for the Detection and Attribution Model Intercomparison Project (DAMIP; Gillett et al., 2016) imply for the ECS, TCR, and the aerosol-driven cooling, which partially offsets global warming.

Decreases in future aerosol loading are thought to contribute to projected warming (Andreae et al., 2005), but the size of this effect is highly uncertain (Forster et al., 2021; Watson-Parris and Smith, 2022). In individual CMIP6 models, mismatches between observed and simulated historical global-mean surface temperatures have been associated with a misrepresentation of the aerosol-induced cooling (Andrews et al., 2020; Smith and Forster, 2021; Golaz et al., 2022), such that despite the increases in the mean ECS and TCR characterizing the Coupled Model Intercomparison Project Phase 6 (CMIP6) ensemble of models relative to CMIP5, the simulated historical warming in CMIP6 is actually smaller than in CMIP5 (Flynn and Mauritsen, 2020; Smith and Forster, 2021; Flynn et al., 2023). Smith and Forster (2021) conclude that differences between CMIP5 and CMIP6 historical simulations are due to an increased ensemble-mean climate sensitivity in CMIP6 versus CMIP5, compensated by a marginally increased aerosol radiative forcing and associated cooling. This causes lower temperatures during 1960–1990 and a larger post-1990 warming trend in the CMIP6 models (Flynn et al., 2023).

The CMIP6 generation of climate model simulations differs from previous generations in some important respects.

• The CMIP6 ensemble now contains several single-model “large” ensembles (i.e. consisting of at least 10 ensemble members, Table 1), including for some single-forcing experiments.

• Studies targeting CMIP5 generally only used single-forcing simulations for GHG and natural influences, inferring the influence of aerosols from these and the historical, all-forcing simulations (Schurer et al., 2018). Here I will exploit single-forcing simulations explicitly targeting anthropogenic aerosols (as have Gillett et al., 2021). Using these “hist-aer” simulations means that I will explicitly test for and not simply assume additivity, as has underlain these earlier studies.

• Since the 1990s, aerosol forcing has been on a declining trend (Szopa et al., 2021; Hodnebrog et al., 2024). This trend reversal means that the aerosol influence is becoming more easily distinguishable from the dominant influence of the ever-increasing GHG-induced warming. The additional years in the CMIP6 DAMIP simulations since this turnaround might critically improve the detectability of the aerosol influence.

Meehl et al. (2020)Dix et al. (2019a, b, 2020a, b, c)Meehl et al. (2020)Ziehn et al. (2019a, b, 2020a, b, c)Meehl et al. (2020)Wu et al. (2018, 2019a, b, c)Xin et al. (2019)Meehl et al. (2020)Swart et al. (2019a, b, c, d, e)Meehl et al. (2020)Danabasoglu (2019a, b, c, 2020)Meehl et al. (2020)Voldoire (2018, 2019a, b, c, d)E3SM (2022e)E3SM (2022a, b, c, d)Scafetta (2023)Li (2019, 2020a, b, c, d)Meehl et al. (2020)Krasting et al. (2018)Horowitz et al. (2018a, b, c)John et al. (2018)Meehl et al. (2020)NASA/GISS (2018a, b, c, d)Meehl et al. (2020)Ridley et al. (2019); Jones (2019a, b, c); Good (2019)Meehl et al. (2020)Boucher et al. (2018a, b, c, d, 2019)Meehl et al. (2020)Tatebe and Watanabe (2018)Shiogama (2019a, b, c)Shiogama et al. (2019)Meehl et al. (2020)Wieners et al. (2019)Müller et al. (2019a, b, c, d)Meehl et al. (2020)Yukimoto et al. (2019a, b, c, d, e)Meehl et al. (2020)Seland et al. (2019a, b, c, d, e)

Table 1 DAMIP models, ensemble sizes of the experiments, ECS values of the models from the literature, and references for the ECS values and the model data. MPI-ESM1-2-LR actually has not completed hist-nat. For this model, the hist-nat temperature change is inferred from its “hist-sol” and “hist-volc” ensembles (Gillett et al., 2016).

A fundamental problem remains: that nature has produced only one realization that might deviate in unknown ways from what the average of a large ensemble of hypothetical realizations would indicate. This problem probably requires a probabilistic approach to tackle, which would be a further development of the method laid out below. I will return to this problem in Sect. 4.

Below I develop a regression approach to evaluate, for all models participating in DAMIP, the fidelity of the simulation of both GHG-driven surface warming and aerosol-related cooling in the CMIP6 ensemble.

2 Data and method

2.1 DAMIP models and experiments

Simulations produced for DAMIP, historical simulations (Eyring et al., 2016) extended to 2020 using the SSP2-4.5 simulations (Gidden et al., 2019), and published ECS and TCR values that are based on 4×CO2 and 1pctCO2 simulations (Eyring et al., 2016) form the basis of this analysis. I use all 16 models for which “hist-GHG”, hist-aer, and “hist-nat” simulations are available (Table 1). The hist-GHG, hist-aer, and hist-nat experiments are identical to the historical coupled all-forcing simulations, except that forcings other than the GHGs, aerosols or their precursors, and natural (solar or volcanic) influences, respectively, are held at their 1850 values (Gillett et al., 2016).

The models that have completed these simulations span large ranges of ECSs (between 2.5 and 5.6 K) and TCRs (between 1.5 and 2.7 K). For 12 models, the simulation period for DAMIP simulations is 1850–2020, but 4 models (CESM2, E3SM-2-0, GISS-E2-1-G, and MPI-ESM1-2-LR) end their DAMIP simulations in 2014.

Furthermore as an observational reference, I use the HadCRUT5 global-mean temperature climatology (Morice et al., 2021). The choice of this reference is not crucial; other reference datasets yield very similar results to those shown here. I form ensemble-, global-, and annual-mean temperature time series from the available simulations.

HadCRUT5 is an amalgamation of sea-surface temperatures and near-surface air temperatures over land and sea ice (Morice et al., 2021). For a rigorous comparison with the temperature fields provided by the CMIP6 models, I thus form, individually for every model, similar amalgamates of “surface temperature” (ts) and “surface air temperature” (tas) using

$$\begin{array}{}\text{(1)}& T=\stackrel{\mathrm{‾}}{(\mathrm{1}-f-s)T_{\mathrm{s}}+(f+s)T_{\text{sa}}},\end{array}$$

where T is the global-, annual-, and ensemble-mean temperature used in the analysis below; T_s is the monthly-mean latitude–longitude resolved surface temperature; T_sa is the monthly-mean surface air temperature; f is the land fraction at every grid point; s is the monthly-mean sea ice fraction; and the overline marks global and annual averaging. For sea ice, I use sea ice concentration on the ocean grid (siconc), regridded to the atmosphere grid using nearest-neighbour interpolation. Only for the GISS-E2-1-G model do I use sea ice concentration on the atmosphere grid (siconca) because of the unavailability of the siconc variable.

Furthermore, to reduce the influence of interannual variations, I apply a 15-year boxcar filter to the hist-aer and hist-GHG ensemble means but not to the hist-nat, historical, and HadCRUT5 datasets. This choice is based on the understanding that interannual variations in hist-aer and hist-GHG reflect random climate variability, which I do not expect to correlate with the climate variability in the historical ensemble or in HadCRUT5. Volcanic eruptions produce forced variations on the annual timescale in the hist-nat and historical ensembles and in HadCRUT5, which I want to preserve, hence the asymmetric application of the filter.

Figure 1 summarizes the behaviour of surface/surface air temperature in the 16 models over the historical period. Panel (a) indicates that there is an approximate proportionality between the simulated warming attributable to GHGs (as taken from the hist-GHG ensembles) and the models' tabulated ECS and TCR values. Furthermore, panel (b) indicates that with some notable exceptions, the models are mostly additive in the sense that the sum of the warming simulated in the three DAMIP experiments is mostly quite similar to the warming simulated in the models' historical ensembles.

Figure 1 (a) Simulated global- and annual-mean warming for 2000–2014 in the hist-GHG ensembles of the DAMIP models relative to the 1850–1899 mean versus their ECSs (black) and TCRs (violet, all in K). The width of the horizontal lines corresponds to $\text{SD}\left(T_i\right)/\sqrt{\mathrm{15}}$, where the T_i values are the annual-mean temperatures for 2000–2014. Solid lines: best-estimate proportional fits. The models' names are abbreviated to three characters. AC1 is ACCESS-ESM1-5 and AC2 is ACCESS-CM2. Also stated are the best-fit proportionality constants and correlation coefficients. (b) Simulated global-mean warming between 1850–1899 and 2000–2014 in the historical ensembles versus the sum of the warming simulated in the respective hist-GHG, hist-aer, and hist-nat ensembles. The solid line marks the diagonal and dashed lines the 80 % and 120 % lines. The lengths of the bars in both directions correspond to the statistical uncertainties at 68 % confidence. Models marked in violet are excluded from the multi-model mean calculations because they do not satisfy the additivity constraints (Eqs. 4 and 5).

2.2 Regression models

Using a linear regression approach, I derive rescaling factors α₁, β₁, γ₁, α₂, β₂, and γ₂ for the temperature responses to GHG, aerosol, and natural forcings, such that the resultant sums of the rescaled simulated temperature anomalies minimize the root-mean-squared deviations ϵ₁ and ϵ₂ versus the ensemble-, global-, and annual-mean historical temperature anomaly T_hist and the HadCRUT5 temperature anomaly record T_obs, respectively, over the period 1850–2020 or 2014 (171 or 165 years):

$$\begin{array}{}\text{(2)}& {\displaystyle }& {\displaystyle }{T}_{\text{obs}}={\mathit{\alpha }}_{\mathrm{1}}{T}_{\text{hGHG}}+{\mathit{\beta }}_{\mathrm{1}}{T}_{\text{haer}}+{\mathit{\gamma }}_{\mathrm{1}}{T}_{\text{hnat}}+{\mathit{\delta }}_{\mathrm{1}}+{\mathit{ϵ}}_{\mathrm{1}}\text{ and}\text{(3)}& {\displaystyle }& {\displaystyle }{T}_{\text{hist}}={\mathit{\alpha }}_{\mathrm{2}}{T}_{\text{hGHG}}+{\mathit{\beta }}_{\mathrm{2}}{T}_{\text{haer}}+{\mathit{\gamma }}_{\mathrm{2}}{T}_{\text{hnat}}+{\mathit{\delta }}_{\mathrm{2}}+{\mathit{ϵ}}_{\mathrm{2}}.\end{array}$$

T_hGHG, T_haer, and T_hnat are all normalized relative to their 1850–1899 averages. δ₁ and δ₂ are intercepts that account for uncertainties (due to natural variability and other factors) in this normalization process; i.e. the other regression coefficients are insensitive to the normalization. ϵ₁ and ϵ₂ are the regression residual time series. This approach is similar to Gillett et al. (2021).

In the absence of climatological noise and if the regression model was complete, for a model which is perfectly additive in the anthropogenic (GHG and aerosol) and natural forcings, the regression coefficients α₂, β₂, and γ₂ would be 1. However, omitted here is the influence of ozone. Since 1900, the temperature changes due to ozone have been around −30 % of those of the aerosols in the Intergovernmental Panel on Climate Change's (IPCC's) best estimate (Fig. 7.8 of Forster et al., 2021). Therefore, the terms β₁T_haer and β₂T_haer will capture the contributions of both aerosol and ozone changes (i.e. the “near-term climate forcers”, NTCFs) to the evolutions of T_obs and T_hist. Because of the offsetting role of ozone, I thus expect β₂<1 for the perfect model, whereas for the GHG influence I expect α₂≈1.

With this understanding, I introduce thresholds

$$\begin{array}{}\text{(4)}& \mathrm{0.8}\le {\mathit{\alpha }}_{\mathrm{2}}\le \mathrm{1.2}\end{array}$$

and

$$\begin{array}{}\text{(5)}& \mathrm{0.6}\le {\mathit{\beta }}_{\mathrm{2}}\le \mathrm{1.1}\end{array}$$

to identify and define models that satisfy additivity and use only those models satisfying both criteria in an “emergent-constraint” approach to calculate best estimates and uncertainty ranges for the ECS, TCR, and aerosol-induced cooling (see below). This is the main difference with respect to Gillett et al. (2021), who use all models without considering their additivity properties. Essentially these two conditions remove models from multi-model emergent-constraint calculations if their regression coefficients versus T_hist deviate substantially from expectations. I will discuss the role of ozone separately below.

I note that there is considerable joint uncertainty resulting from substantial anticorrelations between T_hGHG and T_haer but much smaller correlations (and practically no joint uncertainty) between T_hGHG and T_nat and between T_haer and T_nat. This allows me to simplify the analysis and focus in the following only on the GHG and aerosol influences. To study their joint uncertainty, I calculate, as a function of α, β, and time, the regression error time series

$$\begin{array}{}\text{(6)}& r_{\text{obs}}(\mathit{\alpha },\mathit{\beta })=T_{\text{obs}}-\mathit{\alpha }{T}_{\text{GHG}}-\mathit{\beta }{T}_{\text{haer}}-{\mathit{\gamma }}_{\mathrm{1}}{T}_{\text{hnat}}-{\mathit{\delta }}_{\mathrm{1}}.\end{array}$$

I thus define an error function,

$$\begin{array}{}\text{(7)}& E_{\text{obs}}(\mathit{\alpha },\mathit{\beta })\equiv \sqrt{{\displaystyle \frac{\stackrel{\mathrm{‾}}{r_{\text{obs}}^{\mathrm{2}}(\mathit{\alpha },\mathit{\beta })}}{\stackrel{\mathrm{‾}}{{\mathit{ϵ}}_{\mathrm{1}}^{\mathrm{2}}}}}},\end{array}$$

where the overbar denotes the 171- and 165-year means over 1850–2020 and 1850–2014, respectively.

I define that a regression fit differs from the optimal fit of Eq. (2) if E_obs>E_max, where E_max is a value of the error function, to be determined below, where the fits associated with such an rms residual differ significantly from the optimal fit.

Substituting $Q_{\text{obs}}(\mathit{\alpha },\mathit{\beta })=\stackrel{\mathrm{‾}}{r_{\text{obs}}^{\mathrm{2}}(\mathit{\alpha },\mathit{\beta })}$ with $\mathrm{\Delta }\mathit{\alpha }=\mathit{\alpha }-{\mathit{\alpha }}_{\mathrm{1}}$ and $\mathrm{\Delta }\mathit{\beta }=\mathit{\beta }-{\mathit{\beta }}_{\mathrm{1}}$, I express Q_obs as a quadratic form expanded around the minimum:

$$\begin{array}{}\text{(8)}& Q_{\text{obs}}(\mathit{\alpha },\mathit{\beta })=\stackrel{\mathrm{‾}}{{\mathit{ϵ}}_{\mathrm{1}}^{\mathrm{2}}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}(\mathrm{\Delta }\mathit{\alpha },\mathrm{\Delta }\mathit{\beta })\mathbf{\text{M}}\left(\begin{array}{c}\mathrm{\Delta }\mathit{\alpha }\\ \mathrm{\Delta }\mathit{\beta }\end{array}\right),\end{array}$$

where terms linear in (Δα,Δβ) are 0 because the expansion is around the minimum of Q_obs and

$$\begin{array}{}\text{(9)}& {\displaystyle }\begin{array}{rl}{\displaystyle }\mathbf{\text{M}}& {\displaystyle }=\left(\begin{array}{cc}\frac{{\partial }^{\mathrm{2}}}{\partial {\mathit{\alpha }}^{\mathrm{2}}}& \frac{{\partial }^{\mathrm{2}}}{\partial \mathit{\alpha }\partial \mathit{\beta }}\\ \frac{{\partial }^{\mathrm{2}}}{\partial \mathit{\alpha }\partial \mathit{\beta }}& \frac{{\partial }^{\mathrm{2}}}{\partial {\mathit{\beta }}^{\mathrm{2}}}\end{array}\right)Q_{\text{obs}}\\ {\displaystyle }& {\displaystyle }=\mathrm{2}\left(\begin{array}{cc}\stackrel{\mathrm{‾}}{T_{\text{hGHG}}^{\mathrm{2}}}& \stackrel{\mathrm{‾}}{T_{\text{hGHG}}{T}_{\text{haer}}}\\ \stackrel{\mathrm{‾}}{T_{\text{hGHG}}{T}_{\text{haer}}}& \stackrel{\mathrm{‾}}{T_{\text{haer}}^{\mathrm{2}}}\end{array}\right).\end{array}\end{array}$$

The Hesse or curvature matrix M is characterized by its two positive eigenvalues, λ₁ and λ₂, and associated eigenvectors, e₁ and e₂, with λ₁<λ₂. The extreme case of T_haer∼T_hGHG, i.e. M is degenerate, would imply λ₁=0. This is not actually the case for any of the models considered here, but the two regressors are similar enough that λ₁ is close to 0. The analysis implies that the error function E_obs forms ellipses around the minimum with two orthogonal axes that point in the directions of the eigenvectors, with curvatures in these directions proportional to $\sqrt{{\mathit{\lambda }}_{\mathrm{1}}}$ and $\sqrt{{\mathit{\lambda }}_{\mathrm{2}}}$.

I interpret the eigenvectors e₁ and e₂ as the directions in (α,β) parameter space that correspond to optimal cancellation (for e₁) and optimal reinforcement (for e₂) of the warming effects due to GHGs and NTCFs. For the case of optimal cancellation, GHG warming and NTCF-induced cooling are statistically distinguishable in the observed temperature record because of a trend reversal in SO₂ precursor emissions in the late 20th century (Szopa et al., 2021), causing anthropogenic aerosol-induced cooling to be on a declining trend (in absolute terms) since then, in contrast to the monotonically increasing warming since 1850 associated with GHGs. This means that variations in the contributions of both processes in the direction of optimal cancellation cause detectable variations in the temperature trend of the final 20 years of the regression fit (2001–2020 or 1995–2014) that I will relate to the trend uncertainty in the observed temperature record. This analysis will define bounds E_max on the cost function E_obs and consequently the regression parameters (α₁,β₁). The analysis implies that regression parameters outside these bounds yield significantly and detectably inferior regression fits.

Variations in the direction of optimal reinforcement (e₂), by contrast, produce shifts in the regression fits away from the optimum in either direction. By comparing these shifts to the uncertainty in the mean of the detrended 2001–2020 (or 1995–2014) global temperature record (∼0.03 K in HadCRUT5), I find bounds on the cost function that are substantially more restrictive than the bounds associated with variations in the direction of optimal cancellation discussed above. I will therefore only present an analysis of variations in the direction of cancellation, e₁, which yields wider, more conservative error bounds.

Analogously, I define an error function E_hist(α,β) for the regressions to the historical ensemble means:

$$\begin{array}{}\text{(10)}& E_{\text{hist}}(\mathit{\alpha },\mathit{\beta })\equiv \sqrt{{\displaystyle \frac{\stackrel{\mathrm{‾}}{r_{\text{hist}}^{\mathrm{2}}(\mathit{\alpha },\mathit{\beta })}}{\stackrel{\mathrm{‾}}{{\mathit{ϵ}}_{\mathrm{2}}^{\mathrm{2}}}}}},\end{array}$$

where

$$\begin{array}{}\text{(11)}& r_{\text{hist}}(\mathit{\alpha },\mathit{\beta })=T_{\text{hist}}-\mathit{\alpha }{T}_{\text{GHG}}-\mathit{\beta }{T}_{\text{haer}}-{\mathit{\gamma }}_{\mathrm{2}}{T}_{\text{hnat}}-{\mathit{\delta }}_{\mathrm{2}}.\end{array}$$

The ellipses spanned by E_hist have the same orientations and aspect ratios of the two main axes as those spanned by E_obs.

3 Evaluation

3.1 Example model calculations

As an example, Fig. 2 shows the result of the analysis for the HadGEM3-GC31-LL model. GHGs drive a warming of around 2 K in this model over 1850–2020 (light-green line), offset by aerosol-driven cooling of around −0.9 K by 2020 (dark-green line). Natural influences explain the temporary features associated with volcanic eruptions and solar forcing (blue line). Optimal regression parameters α₂, β₂, and γ₂ for T_hist are close to 1 (i.e. HadGEM3-GC31-LL is nearly additive; violet line). However, the regression against T_obs requires substantial reductions in the parameters describing both the GHG and the aerosol influences (α₁ and β₁), to the point that the aerosol cooling would need to be reduced by around 75 % and the GHG influence by 30 % to match the observed record (orange line).

Figure 2 Ensemble-, global-, and annual-mean temperature anomalies relative to the 1850–1899 average for HadGEM3-GC31-LL. Black symbols and dark-green, light-green, and blue lines: the DAMIP and historical/SSP2-4.5 ensemble means as indicated. The hist-GHG and hist-aer temperature evolutions have been smoothed using a 15-year boxcar filter. Violet: optimal regression fits to T_hist following Eq. (3). Orange: optimal regression fits to T_obs, the HadCRUT5 reconstruction, following Eq. (2). Red symbols: HadCRUT5 (Morice et al., 2021). The regression coefficients α₁, β₁, γ₁, δ₁, α₂, β₂, γ₂, and δ₂ that are stated in orange and violet are as defined in Eqs. (2) and (3).

Figure A1 contains equivalent plots for the remaining 15 models.

3.1.1 When are two regression fits statistically indistinguishable?

I note that the observational record exhibits a nearly linear warming trend towards the end of the record (Fig. 2). Furthermore, during much of the 20th century the aerosol-induced cooling is directly opposed to the GHG warming, but in the 1990s its trend changes sign in the HadGEM3-GC31-LL hist-aer ensemble. This trend reversal is a major reason that the GHG and aerosol influences are at all statistically distinguishable in the historical record. I thus define two regression fits to be significantly different if their 20-year trends for 2001–2020 (or 1995–2014 for CESM2, E3SM-2-0, GISS-E2-1-G, and MPI-ESM1-2-LR) differ by more than the observational uncertainties at 95 % confidence in these trends (κ=6.7 and 6.3 mK a⁻¹, respectively).

I evaluate the regression fits in regression parameter space (α,β) along the lines that correspond to optimal cancellation of the warming and cooling impacts of GHGs and aerosols, respectively. This is the line spanned by the eigenvector corresponding to the smaller eigenvalue λ₁, e₁; i.e.

$$\begin{array}{}\text{(12)}& (\mathit{\alpha },\mathit{\beta })=({\mathit{\alpha }}_{\mathrm{1}},{\mathit{\beta }}_{\mathrm{1}})+c\cdot {\mathit{e}}_{\mathbf{1}}.\end{array}$$

This line marks the direction of maximum joint uncertainty in the regression parameters.

I plot the error functions E_obs against the 2001–2020 (or 1995–2014) trends in the associated fits, $T(\mathit{\alpha },\mathit{\beta })=\mathit{\alpha }{T}_{\text{hGHG}}+\mathit{\beta }{T}_{\text{haer}}+{\mathit{\gamma }}_{\mathrm{1}}{T}_{\text{hnat}}+{\mathit{\delta }}_{\mathrm{1}}$, evaluated along the line described by Eq. (12) (Fig. 3).

Figure 3 Error function E_obs as a function of the 2001–2020 (for CESM2, E3SM-2-0, and GISS-E2-1-G: 1995–2014) temperature trends in the regression fits to T_obs (Eq. 2) for the DAMIP models. The trends are evaluated along the line in (α,β) parameter space described by Eq. (12). Solid vertical line: trend in the optimal fit (that minimizes E_obs). Dashed vertical lines: trends in the sub-optimal fits that differ from the optimal trend by the observational trend uncertainty κ. Horizontal line: value of the error function corresponding to this trend uncertainty. The number in the titles is the value of the error function at these points.

By evaluating E_obs at the two trend values that differ from the trend in the optimal solution by κ, I find model-dependent values for E_max. For all but three of the models, E_max≤1.13. These models (BCC-CSM2-MR, FGOALS-g3, and NorESM2-LM) do not simulate the trend change in the aerosol-induced cooling characterizing the other models (Fig. A1). The other 12 models yielding smaller, regular values for E_max all simulate a substantial change in the rate of cooling, such that their hist-aer temperature time series become statistically independent from their hist-GHG temperatures, and almost all exhibit warming trends during the final 2 decades of their hist-aer ensembles (Figs. 2 and A1).

3.2 Joint uncertainty analysis of the GHG and aerosol influences for all models

Figure 4 illustrates firstly that additivity does not extend to all models; i.e. the centres of many open ellipses are outside the “additivity rectangles”. However, a subgroup of eight models does satisfy Eqs. (4) and (5) (ACCESS-CM2, ACCESS-ESM1-5, CESM2, E3SM-2-0, GISS-E2-1-G, HadGEM3-GC31-LL, MIROC6, and MPI-ESM1-2-LR). I note that this criterion disqualifies four models which also do not satisfy additivity in their simulated warming as expressed in Fig. 1b, i.e. CanESM5, GFDL-ESM4, MRI-ESM2-0, and NorESM2-LM.

Figure 4 Error functions E_obs (colours; Eq. 7) and E_hist (contours; Eq. 10) where these functions are smaller than $max(E_{\text{max}},\mathrm{1.1})$ for the 16 DAMIP models. Rectangle: window of additivity defined by Eqs. (4) and (5). “Aer. cooling” is the global-, ensemble-, and annual-mean cooling for 2000–2014 relative to 1850–1899, as simulated in the models' hist-aer ensembles. “ECS” is as in Table 1.

In all cases except for three models (BCC-CSM2-MR, FGOALS-g3, and MPI-ESM1-2-LR), there are no substantial overlaps between the regression uncertainty ellipses for the fits to T_obs and T_hist. This means that 13 of the models have systematic differences between the simulated historical and the observed temperature evolutions large enough to show in this lack of overlap in the regression parameters. This implies irreconcilable scaling factors α and β for the GHG or NTCF influences, or for both.

Specifically regarding the GHG scaling factor α (plotted on the horizontal axes in Fig. 4), for some large-ECS models including ACCESS-CM2, CanESM5, CESM2, HadGEM3-GC31-LL, and IPSL-CM6A-LR, the analysis suggests that α₁<1; i.e. to better match the HadCRUT5 time series, the GHG influences in these models need to be scaled down (Gillett et al., 2021).

Lastly, for almost all models the analysis suggests that the NTCF scaling factor β₁<0.6; i.e. the filled ellipses are centred below the rectangles of additivity in Fig. 4. Exceptions are the MIROC6, MPI-ESM1-2-LR, and MRI-ESM2-0 models where the regression yields an NTCF influence consistent with no rescaling. Exaggeration of the NTCF influence is large and unambiguous for ACCESS-CM2, ACCESS-ESM1-5, CanESM5, E3SM-2-0, HadGEM3-GC31, and NorESM2-LM; these models all simulate at least 0.66 K of cooling in their hist-aer ensembles (Gillett et al., 2021).

3.3 Emergent constraints for the ECS and the aerosol cooling influence

3.3.1 The GHG influence

Figure 5 shows the result of the regression analysis (Sect. 2.2) for all models. The ensemble includes three models (CanESM5, CESM2, and HadGEM3-GC31-LL) that have ECSs exceeding the “very likely” range given by AR6 (2–5 K; Forster et al., 2021; Fig. 5a). For these three models, the GHG correction factors are α₁<1 (i.e. reductions in the GHG influences would bring their historical simulations into better agreement with the observations). At the other end of the spectrum, MIROC6, with an ECS of 2.6 K, requires an increase in the GHG-induced warming by 23 % to bring its historical evolution into agreement with HadCRUT5. In general, the distribution (panel a) can be approximated by $\text{ECS}\sim {\mathit{\alpha }}_{\mathrm{1}}^{-\mathrm{1}}$. Equivalently, panel (b) shows the ECS×α₁ versus α₁. The thus “adjusted” ECSs (i.e. ECS×α₁) are now within the AR6 “likely” range (2.5–4 K) for all eight models that satisfy additivity, and the multi-model spread of these models' ECSs is smaller than the AR6 uncertainty range. I obtain a multi-model-mean adjusted ECS of 3.5 ± 0.4 K at 68 % confidence (Fig. 5b). Replacing the ECSs with TCRs in the above analysis gives essentially the same result. Two models (CanESM5 and HadGEM3-GC31-LL) have TCRs outside the very likely AR6 range (1.2–2.4 K). Multiplying the TCR by α₁ yields adjusted TCRs that are now almost all within the likely range of AR6 (1.4–2.2 K) for models satisfying additivity. The multi-model mean adjusted TCR of 1.8 ± 0.3 K compares very well to the AR6 estimate.

Figure 5 (a) The models' ECSs (K) versus α₁ (Eq. 2). The lengths of the horizontal lines depict the regression uncertainties at 68 % confidence. Solid line: regression fit assuming $\text{ECS}\cdot {\mathit{\alpha }}_{\mathrm{1}}=\text{const.}$ Dashed and dotted lines: uncertainty ranges at 68 % and 95 % confidence. The yellow and green regions are the likely (i.e. 66 % confidence) and very likely (90 %) ECS intervals assessed by AR6 (Forster et al., 2021). (b) Same as (a) but for ECS⋅α₁ versus α₁. The solid, dashed, and dotted lines are the mean, 68 %, and 95 % confidence intervals. (c) and (d): same as (a) and (b) but for the TCR. (e) The tabulated ECSs versus cooling simulated in the hist-aer ensembles for 2000–2014 relative to 1850–1899. Solid line: best linear fit. (f) The aerosol-induced cooling taken from hist-aer times β₁ (K) versus the correction factors to aerosol cooling β₁ (Eq. 2). Solid and dashed blue lines denote the means and 68 % confidence limits of both quantities. The black line and the green box are the AR6 best estimate and the 5 % to 95 % uncertainty range for the temperature change due to NTCFs (Fig. 7.8 of Forster et al., 2021). The numbers in black denote the means and standard deviations of β₁ and of aerosol-induced cooling times β₁ (K). The models marked in violet are excluded from this averaging.

3.3.2 The aerosol influence

The second part of this analysis concerns the aerosol-induced cooling. There is an anticorrelation (with a correlation coefficient of −0.49) between the cooling attributable to anthropogenic aerosol increases, as discerned from hist-aer, and the ECSs of the 16 models (Fig. 5e). This means that large-ECS models tend to compensate for some of their GHG-induced warming by simulating relatively strong cooling due to aerosols. In other words, the biases in both properties are coupled.

Focussing here only on the eight models that satisfy additivity, their NTCF rescaling factors, β₁, are in the range of 0.1 to 1.25 with a mean (standard deviation) of 0.47 (0.39). In this group, the MPI-ESM1-2-LR and MIROC6 models are the only models that are consistent with no rescaling of their NTCF-induced cooling. Average adjusted aerosol-induced cooling between 1850–1899 and 2000–2014 in this group amounts to 0.24 ± 0.11 K, when their average unadjusted cooling is 0.67 ± 0.31 K (both at 68 % confidence). For comparison, this places my analysis at the lower end of the AR6 estimate of 0.31 (0.15 to 0.57) K (at 5 % to 95 % confidence) of cooling due to aerosols and ozone combined. The AR6 range is inferred from the data accompanying Fig. 7.8 of Forster et al. (2021). If the models were perfect and the ozone influence was proportional to that of the aerosols, offsetting around 30 % of the aerosol-induced cooling, a value of β₁≈0.7 would be expected, but only MIROC6 and GISS-E2-1-G get close to this. The fact that five of the other models have rescaling factors, β₁, in the range 0.1 to 0.4 makes it implausible that ozone is the explanation here. Four of these models (MPI-ESM1-2-LR, MIROC6, GISS-ES-1-G, and ACCESS-ESM1-5) agree within their 68 % uncertainty ranges with the AR6 best estimate (0.31 K). The other four (CESM2, E3SM-2-0, ACCESS-CM2, and HadGEM3-GC31-LL) require a smaller NTCF-induced cooling than that. Flynn et al. (2023) find that models that better reproduce the observed temperature evolution all simulate relatively small aerosol-induced cooling, in agreement with this study, and Gillett et al. (2021) also find small rescaling factors for the aerosol influence for several of the same models used here.

4 Discussion and conclusions

Mismatches between the observed global-mean surface temperature and CMIP6 simulated historical temperature have been documented before and attributed to, in some cases, a deficient simulation of aerosol-related cooling (Andrews et al., 2020; Flynn and Mauritsen, 2020; Smith and Forster, 2021; Golaz et al., 2022; Flynn et al., 2023). Here I exploit these mismatches to derive scaling factors for the GHG-induced warming and the aerosol-related cooling that in a hypothetical model would bring the simulated historical temperature into optimal agreement with the HadCRUT5 climatology. I then relate these scaling factors to the warming attributable to GHGs, on the one hand, and on the other hand to the cooling attributable to anthropogenic NTCFs. The GHG scaling factors very approximately follow an inverse relationship to the ECSs of the models, such that the products of the ECSs and the scaling factors are in better agreement with the AR6 evaluation of the planetary ECS than the modelled ECSs themselves. Particularly for three large-ECS models with ECSs outside the AR6 very likely range (Forster et al., 2021), this adjustment brings these ECSs into agreement with the AR6 estimate. Essentially the same holds for the TCR.

These results are consistent with quantifications of the aerosol and GHG influences based on energy balance calculations (e.g. Storelvmo et al., 2016; Smith and Forster, 2021; Smith et al., 2021) but arrived at using an independent approach. Smith and Forster (2021) use energy budget constraints to find that despite reductions in historical aerosol and GHG forcing from CMIP5 to CMIP6, stronger climate feedback in CMIP6 models, which is reflected in the increased ECSs in CMIP6 models, causes both stronger aerosol cooling and, from 1990, increased GHG warming. This analysis generally confirms their results but shows that a tendency to overestimate aerosol-induced cooling pertains not just to the high-ECS models but to many of the moderate-ECS models as well. Only for one model (MPI-ESM1-2-LR) did I find a rescaling factor larger than 1 (Gillett et al., 2021). At 0.26 K, the mean aerosol cooling simulated by MPI-ESM1-2-LR is the second-smallest in the ensemble. Storelvmo et al. (2016) employ a purely observations-based approach to find that aerosols masked approximately one-third of the GHG-induced warming since the 1960s, leaving a TCR of 2 ± 0.8 K. Again, my results are consistent with but largely independent of their results. Smith et al. (2021) use an energy balance approach constrained by CMIP6 input data to quantify the ECS and the TCR at 3.1 and 1.8 K, respectively, in agreement with my results, and give a very wide uncertainty range for the aerosol-induced radiative forcing of −1.8 to −0.5 W m⁻² since 1750. Applying a conversion factor of 0.5 K W⁻¹ m² (Forster et al., 2021) yields 0.25 to 0.9 K of cooling due to aerosol. Further adjusting for some anthropogenic aerosol increase between 1750 and 1850–1899 and accounting for some offset by ozone, my result is consistent but near the low end of this range.

Hodnebrog et al. (2024), in a recent study, attribute much of the observed recent increase in the Earth's energy imbalance to a reduction in the anthropogenic aerosol forcing. While the methodology used by these authors differs from my approach, this finding is likely inconsistent with the conclusion arrived at here that aerosols, in the adjusted multi-model mean, exert a substantially smaller influence on climate than the ever-increasing GHGs. It also appears inconsistent with the plain results of the hist-aer and hist-GHG simulations (Figs. 2 and A1), where the warming attributable to aerosols over 2001–2019 is generally smaller than that simulated in the hist-GHG ensembles over the same period. This is even before this warming is scaled down, as is laid out above. More research is required to reconcile these findings.

I find an anticorrelation between the total simulated aerosol-induced cooling since 1850 and the ECSs of the models, suggesting that both quantities are coupled and that there could be a compensation of errors between these two processes. The eight models considered here that satisfy additivity span a considerable range in simulated aerosol-induced cooling for 1850–1899 to 2000–2014, namely 0.3 to 1.2 K, with a mean (standard deviation) of 0.67 K (0.31 K). Applying the correction factor β₁ reduces the range to 0.04 to 0.36 K, with a mean (standard deviation) of 0.24 (0.11) K. This is now interpreted as the cooling due to aerosols between 1850–1899 and 2000–2014 offset by warming due to ozone in the same period.

There are several limitations to the analysis presented here. The first is that the hist-GHG experiment quantifies the responses of the climate models to all GHGs in combination, whereas the ECSs and TCRs are expressions of the sensitivity of climate to CO₂ increases only. I have shown that there are near-perfect proportionalities and high degrees of correlation (0.80 and 0.87, respectively) between the warming simulated in hist-GHG and the ECSs and TCRs in the 16 models used here (Fig. 1), suggesting that the substantial model diversities that exist for these quantities are due to the same processes, i.e. climate feedback (e.g. due to cloud adjustments) that is not sensitive to the detailed properties of the driving GHGs.

A further, more fundamental limitation is that the models do not respond perfectly additively to GHG and aerosol forcing. This is expressed in deviations from 1 for the α₂, β₂, and γ₂ parameters in Eq. (3) (Fig. 4). The non-additivity is due to a variety of reasons, including forcings not included in the analysis (such as land use and ozone changes). I have tested the sensitivity of the results to including the impact of warming due to ozone changes using the “hist-totalO3” simulations that were produced under DAMIP and by correspondingly expanding Eqs. (2) and (3) by a term for the temperature anomalies simulated under this experiment. Five models (CanESM5, GISS-E2-1-G, HadGEM3-GC31-LL, MIROC6, and MPI-ESM1-2-LR) have completed this experiment, but two of these (CanESM5 and MPI-ESM1-2-LR) are “non-additive” in the expanded regression model. Ozone changes are the most important anthropogenic radiative forcing agent after those considered here (GHGs and aerosols; Forster et al., 2021). Repeating the analysis based on just GISS-E2-1-G, HadGEM3-GC31-LL, and MIROC6 but with a term added to the regression models for ozone-induced warming, I obtain quite similar values for the ECS (3.2 ± 1 K), TCR (1.8 ± 0.4 K), and aerosol-induced cooling (−0.29 ± 0.11 K). This suggests that ozone forcing is the leading explanation for neither the non-additivity nor the weak aerosol-induced cooling found here. More models completing the hist-totalO3 simulations would help. Furthermore, a dedicated experiment for the NTCFs (i.e. focusing on aerosol and ozone forcings combined) would also be very helpful.

A further reason for non-additivity could be a substantial influence of random variations on the regression. I note that the model with the joint-smallest hist-aer ensemble size (GFDL-ESM4) exhibits too substantial a non-additivity to be included in the multi-model averages (Sect. 3.3), and the models with the best additivity (HadGEM3-GC31-LL, MIROC6, and MPI-ESM1-2-LR) have all contributed large ensembles. However, there is not a consistent association of additivity with ensemble sizes. For example, CanESM5 has large ensembles for all experiments but still exhibits substantial non-additivity. Quite a few models have ensemble sizes of three for their hist-aer or hist-GHG ensembles, but this group includes some models with good and with poor additivity. In some cases, the non-additivity can be traced to an extremely small eigenvalue, λ₁, of the covariance matrix M. For HadGEM3-GC31-LL, as an example of an additive model, I find λ₁=0.02. Examples of poor additivity include FGOALS-g3 (${\mathit{\lambda }}_{\mathrm{1}}=\mathrm{8}\times {\mathrm{10}}^{-\mathrm{4}}$), IPSL-CM6A-LR (${\mathit{\lambda }}_{\mathrm{1}}=\mathrm{3}\times {\mathrm{10}}^{-\mathrm{3}}$), and NorESM2-LM (${\mathit{\lambda }}_{\mathrm{1}}=\mathrm{6}\times {\mathrm{10}}^{-\mathrm{3}}$). However, the group of additive models also includes MIROC6 ($\mathrm{2}\times {\mathrm{10}}^{-\mathrm{3}}$). Thus, near degeneracy of M also does not completely explain why some models exhibit non-additivity. Other factors may come into play, including the fact that some models simply respond non-linearly to the applied forcings or even that errors exist in the experimental setups. It is beyond the scope of this paper to fully diagnose these occurrences of non-additivity. However, removing such models from the emergent-constraint calculations of Sect. 3.3 substantially improves model consensus.

Several models indicating that large reductions in aerosol cooling would be beneficial for bringing the simulated historical temperature record into better agreement with observations, including ACCESS-ESM1-5, E3SM-2-0, HadGEM3-GC31-LL, and MRI-ESM2-0, all have β₂>0.6; i.e. these models behave relatively additively, and the inference that exaggerated historical aerosol-induced cooling contributes substantially to errors in the simulations of global-mean temperature by these models is quite well founded.

Bellouin et al. (2020) review the constraints posed by observed temperature on the effective radiative forcing of aerosols. Applying a conversion factor of 0.5 K W⁻¹ m² (Forster et al., 2021), I translate their effective radiative forcing estimate into a cooling influence over the historical period. In their discussion, which considers both direct and indirect aerosol effects on radiation, the best-estimate uncertainty interval of aerosol radiative forcing (−1.6 to −0.35 W m⁻², translated into approximately 0.18 to 0.8 K of aerosol-induced cooling) includes most models considered here for both unadjusted and adjusted aerosol-induced cooling, although the best-estimate multi-model mean cooling due to NTCFs inferred here (0.24 K) is near the low end of this range. Accounting for the influence of ozone (Forster et al., 2021), this translates into a best estimate of 0.38 K, which compares better to the Bellouin et al. (2020) headline estimate. As noted above, my estimate for the cooling due to NTCFs skews slightly smaller than the one in AR6 but with overlapping uncertainty ranges.

The results qualitatively confirm Smith and Forster (2021), who find that excessive cooling due to aerosols in 1960–1990 causes cold biases in this period in many CMIP6 historical simulations. My analysis leaves open a question regarding whether the excessive response to aerosol forcing in most models is due to too much aerosol being produced in these models (i.e. a problem with the CMIP6 forcing data) or whether the internal model physics of aerosol–radiation and aerosol–cloud interactions is flawed. The fact that this behaviour is common to most models suggests that the former is a likely factor.

Figure 3 shows that the regression coefficients versus the HadCRUT5 temperature, α₂ and β₂, are subject to substantially larger uncertainties than those versus the ensemble-mean simulated temperatures, α₁ and β₁. This is a reflection of the greater noise associated with observed temperature compared to an ensemble-mean temperature. Despite this larger uncertainty, the analysis does indicate, for most models, differences between these two parameter pairs that are irreconcilable within the statistical uncertainties. In a follow-on paper, I will investigate this aspect further using a probabilistic approach.

In summary, I have used a reconstruction of global-mean merged surface/surface air temperature, DAMIP, and historical simulations by 16 contemporary climate models to derive constraints for the GHG-induced warming and the aerosol-induced cooling, by far the leading influences driving global warming. Using an emergent-constraint approach, I derive a “corrected” ensemble-mean equilibrium climate sensitivity of about 3.5 ± 0.4 K and a corrected TCR of 1.8 ± 0.3 K (both at 68 % confidence), in excellent agreement with the AR6 estimates but with reduced uncertainties (Forster et al., 2021). For the eight models with relatively good additivity, I find that reductions in the NTCF-induced cooling, along with some reductions in the GHG-induced warming for models with large ECSs, would bring their historical simulations into better agreement with the observational record. The results presented here highlight ongoing difficulties in correctly simulating climate feedback in global models. Substantial, systematic, and nearly community-wide issues in representing historical global surface temperature reduce confidence in quantitative projections of global warming by models affected by these problems. Interestingly, at least some CMIP3 models were consistent with observations without any need for rescaling of the aerosol and GHG signatures (Stone et al., 2007b), including the precursor of CESM2 (Stone et al., 2007a). This may suggest that at least for this model, development occurring in the intervening time has introduced this problem.

The analysis is limited by the substantial anticorrelation between the GHG and the aerosol global-mean warming signatures. I anticipate that as anthropogenic aerosol production continues, as projected, to decline in the future (Lee et al., 2021), the anticorrelation between GHG-induced warming and aerosol-induced cooling will reduce, allowing for a more confident attribution of their respective roles in driving global warming.

Appendix A: Regression fits for the remaining models

For models not represented in Fig. 2, the regression fits are given in Fig. A1.

Figure A1 Same as Fig. 2 but for the ACCESS-CM2, ACCESS-ESM1-5, BCC-CSM2-MR, CanESM5, CESM2, CNRM-CM6-1, E3SM-2-0, FGOALS-g3, GFDL-ESM4, GISS-E2-1-G, IPSL-CM6A-LR, MIROC6, MPI-ESM1-2-LR, MRI-ESM2-0, and NorESM2-LM models.