Comparison of CMIP6 historical climate simulations and future projected warming to an empirical model of global climate

The sixth phase of the Coupled Model Intercomparison Project (CMIP6) is the latest modeling effort for general circulation models to simulate and project various aspects of climate change. Many of the general circulation models (GCMs) participating in CMIP6 provide archived output that can be used to calculate effective climate sensitivity (ECS) and forecast future temperature change based on emissions scenarios from several Shared Socioeconomic Pathways (SSPs). Here we use our multiple linear regression energy balance model, the Empirical Model of Global Climate (EM-GC), to simulate and project changes in global mean surface temperature (GMST), calculate ECS, and compare to results from the CMIP6 multi-model ensemble. An important aspect of our study is a comprehensive analysis of uncertainties due to radiative forcing of climate from tropospheric aerosols (AER RF) in the EM-GC framework. We quantify the attributable anthropogenic warming rate (AAWR) from the climate record using the EM-GC and use AAWR as a metric to determine how well CMIP6 GCMs replicate human-driven global warming over the last 40 years. The CMIP6 multi-model ensemble indicates a median value of AAWR over 1975–2014 of 0.221 C per decade (range of 0.151 to 0.299 C per decade; all ranges given here are for 5th and 95th confidence intervals), which is notably faster warming than our median estimate for AAWR of 0.157 C per decade (range of 0.120 to 0.195 C per decade) inferred from the analysis of the Hadley Centre Climatic Research Unit version 5 data record for GMST. Estimates of ECS found using the EM-GC assuming that climate feedback does not vary over time (best estimate 2.33 C; range of 1.40 to 3.57 C) are generally consistent with the range of ECS of 1.5 to 4.5 C given by the IPCC’s Fifth Assessment Report. The CMIP6 multi-model ensemble exhibits considerably larger values of ECS (median 3.74 C; range of 2.19 to 5.65 C). Our best estimate of ECS increases to 3.08 C (range of 2.23 to 5.53 C) if we allow climate feedback to vary over time. The dominant factor in the uncertainty for our empirical determinations of AAWR and ECS is imprecise knowledge of AER RF for the contemporary atmosphere, though the uncertainty due to time-dependent climate feedback is also important for estimates of ECS. We calculate the likelihood of achieving the Paris Agreement target (1.5 C) and upper limit (2.0 C) of global warming relative to pre-industrial for seven of the SSPs using both the EM-GC and the CMIP6 multi-model ensemble. In our model framework, SSP1-2.6 has a 53 % probability of limiting warming at or below the Paris target by the end of the century, and SSP4-3.4 has a 64 % probability of achieving the Paris upper limit. These estimates are based on the assumptions that climate feedback has been and will remain constant over time since the prior temperature record can be fit so well assuming constant climate feedback. In addition, we quantify the sensitivity of future warming to the curbing of the current rapid growth of atmospheric methane and show that major near-term limits Published by Copernicus Publications on behalf of the European Geosciences Union. 546 L. A. McBride et al.: Comparison of CMIP6 warming to an empirical model of global climate on the future growth of methane are especially important for achievement of the 1.5 C goal of future warming. We also quantify warming scenarios assuming climate feedback will rise over time, a feature common among many CMIP6 GCMs; under this assumption, it becomes more difficult to achieve any specific warming target. Finally, we assess warming projections in terms of future anthropogenic emissions of atmospheric carbon. In our model framework, humans can emit only another 150± 79 Gt C after 2019 to have a 66 % likelihood of limiting warming to 1.5 C and another 400±104 Gt C to have the same probability of limiting warming to 2.0 C. Given the estimated emission of 11.7 Gt C per year for 2019 due to combustion of fossil fuels and deforestation, our EM-GC simulations suggest that the 1.5 C warming target of the Paris Agreement will not be achieved unless carbon and methane emissions are severely curtailed in the next 10 years.

larger values of ECS (median 3.74°C; range of 2.19-5.65°C). The dominant factor in the uncertainty for our empirical determinations of AAWR and ECS is imprecise knowledge of AER RF for the contemporary atmosphere. We calculate the likelihood of achieving the Paris Agreement target (1.5°C) 30 and upper limit (2.0°C) of global warming relative to pre-industrial for seven of the SSPs using both the EM-GC and the CMIP6 multi-model ensemble. In our model framework, SSP1-2.6 is the 1.5°C pathway with a 64.8% probability of limiting warming at this level by the end of century and SSP4-3.4 is the 2.0°C pathway, with a 74.0% probability of achieving the Paris upper limit. In addition, we quantify the sensitivity of future warming to the curbing of the current rapid growth of atmospheric methane and show 35 major near-term limits on the future growth of methane are especially important for achievement of the 1.5°C goal of future warming. Finally, we assess warming projections in terms of future anthropogenic emissions of atmospheric carbon. In our model framework, humans can emit only another 268 Gt C after 2019 to have a 66% likelihood of limiting warming to 1.5°C, and another 565 Gt C to have the same probability of limiting warming to 2.0°C. Given the current emission of 11.7 Gt C per year due to 40 combustion of fossil fuels and deforestation, our EM-GC simulations suggest the 1.5°C warming target of the Paris Agreement will not be achieved unless carbon and methane emissions are severely curtailed in the next two decades.

Introduction
The goals of the Paris Agreement, negotiated in December of 2015, are to keep global warming below 45 2.0°C relative to the start of the Industrial Era and pursue efforts to limit global warming to 1.5°C. General circulation models (GCMs) project future temperature change using various evolutions of greenhouse gases and determine the likelihood of achieving the goals of the agreement. Many GCMs are participating in the sixth phase of the Coupled Model Intercomparison Project (CMIP6) to quantify how the models represent different aspects of climate change . Having accurate projections of future 50 temperature is critical for achieving the goals of the Paris Agreement. Chapter 11 of IPCC's Fifth Assessment Report shows that some of the previous generations of these models participating in phase 5 of the Coupled Model Intercomparison Project (CMIP5) (Taylor et al., 2012) tended to overestimate the https://doi.org/10.5194/esd-2020-67 Preprint. Discussion started: 8 September 2020 c Author(s) 2020. CC BY 4.0 License. increase in global mean surface temperature (GMST) for the 21 st century (Kirtman et al., 2013). In this analysis we use a multiple linear regression energy balance model to quantify the change in GMST from 55 1850-2019, project future changes in GMST, compare to the CMIP6 multi-model ensemble, and determine the likelihood of achieving the goals of the Paris Agreement.
Several prior studies have used a multiple linear regression approach to model the GMST anomaly in order to quantify the impact of anthropogenic and natural factors on climate (Foster and Rahmstorf, 2011;Rind, 2008, 2009;Zhou and Tung, 2013). Typically, total solar irradiance, volcanoes, 60 and El Niño southern oscillation (ENSO) are the natural components represented in the multiple linear regression, and greenhouse gases and aerosols are the anthropogenic factors. We use multiple linear regression, in connection with a dynamic ocean module that accounts for the export of heat from the atmosphere to the ocean, to represent the natural and anthropogenic components of the climate system.
In addition to the typical natural factors listed above, we include the Atlantic meridional overturning 65 circulation (AMOC), Pacific decadal oscillation (PDO), and Indian Ocean dipole (IOD) to provide a robust representation of the natural climate system Hope et al., 2017). Our anthropogenic components also include the effect of land use change (i.e., deforestation) on Earth's albedo and the export of heat from the atmosphere to the ocean as the atmosphere warms.
Our analysis builds on the work of Canty et al. (2013) and Hope et al. (2017) and includes several 70 key updates. One is the extension back in time of our analysis to 1850. The Hadley Center Climatic Research Unit (Morice et al., 2012) and Berkley Earth Group (Rohde et al., 2013) provide GMST records starting in 1850, which now allows for a simulation of GMST that covers 170 years. The second update is the use of the Shared Socioeconomic Pathways (SSPs) (O'Neill et al., 2017) as our climate scenarios to designate future evolution of greenhouse gas and aerosol abundances. The third is the adoption of an 75 active upper ocean to our model, formulated in a manner that matches the equations of Bony et al. (2006) and Schwartz (2012). A description of the model, the various input parameters used, and the updates listed above is given in Sect. 2. Section 3 provides results of CMIP6 comparing to the historical climate record, equilibrium climate sensitivity (ECS), as well as comparisons of our model and CMIP6 projections of future GMST change. Discussion of these results is provided in Sect. 4, along with concluding remarks. 80 https://doi.org/10.5194/esd-2020-67 Preprint. Discussion started: 8 September 2020 c Author(s) 2020. CC BY 4.0 License.

Empirical model of global climate
In this analysis we use the empirical model of global climate (EM-GC), which provides a multiple linear regression, energy balance simulation of GMST. As detailed in the following paragraphs, the EM-GC 85 solves for ocean heat uptake efficiency (κ) and six regression coefficients to minimize the cost function in Eq. (1). (1) In this equation, ΔTOBS represents a time series of observed monthly GMST anomalies, ΔTMDL is the modeled monthly change in GMST, σOBS is the 1-sigma uncertainty associated with each temperature 90 observation, i is the index for each month, and NMONTHs is the total number of months used in the analysis.
The observed GMST anomalies in this analysis are blended near surface air and sea surface temperature differences relative to the GMST anomaly over 1850-1900, which is assumed to represent pre-industrial conditions.
We consider several anthropogenic and natural factors as components of ΔTMDL. The radiative 95 forcing (RF) due to greenhouse gases (GHGs), anthropogenic aerosols (AER), land use change (LUC), and the export of heat from the atmosphere to the world's oceans are the anthropogenic components of ΔTMDL. The influence on GMST from total solar irradiance (TSI), El Niño southern oscillation (ENSO), the Atlantic meridional overturning circulation (AMOC), volcanic eruptions that reach the stratosphere and enhance stratospheric aerosol optical depth (SAOD), the Pacific decadal oscillation, (PDO) and the 100 Indian Ocean dipole (IOD) are the natural components of ΔTMDL. Equation (2) In Eq.
of 6 months for SAOD is representative of the time needed for the atmosphere to respond to a change in the aerosol loading due to a volcanic eruption (Douglass and Knox, 2005). This lag is the same as used 110 by Lean and Rind (2008) and Foster and Rahmstorf (2011). The 1 month delay for TSI yields the maximum value of C2, the solar irradiance regression coefficient. Lean and Rind (2008) and Foster and Rahmstorf (2011) also use a 1 month lag for TSI in their analyses. The 2 month delay for the response of GMST to ENSO is the lag needed to obtain the largest value of the correlation coefficient of the Multivariate ENSO Index version 2 (MEI.v2) (Wolter and Timlin, 1993;Zhang et al., 2019) versus the 115 value of TENSO calculated by Thompson et al. (2009). In Thompson et al. (2009), TENSO is the simulated response of GMST to variability induced by ENSO, taking into consideration the effective heat capacity of the atmospheric-ocean mixed layer. Lean and Rind (2008) used a 4-month lag for ENSO.
The term AMOCi represents the influence of the change in the strength of the thermohaline circulation on GMST (Knight et al., 2005;Medhaug and Furevik, 2011;Stouffer et al., 2006;Zhang and 120 Delworth, 2007). We use the Atlantic multidecadal variability, based on the area weighted monthly mean sea surface temperature (SST) in the Atlantic Ocean between the equator and 60°N (Schlesinger and Ramankutty, 1994), as a proxy for the strength of AMOC. A strong AMOC is characterized by northward flow of energy that would otherwise be radiated to space, which occurs in both the ocean and atmosphere and leads to particularly warm summers in Europe (Kavvada et al., 2013) as well as a number of other 125 well documented influences in other climatic regions (Nigam et al., 2011). The total anthropogenic RF of climate is used to detrend the AMOC signal because this method provides a more realistic approach to infer the changes in the strength of AMOC and its effect on GMST than other detrending options .
Our model explicitly accounts for the export of heat from the atmosphere to the world's oceans 140 (i.e., ocean heat export or OHE). The quantity QOCEAN in Eq.
(2) to be comparable to the formulation for this term used by Bony et al. (2006) and Schwartz (2012 The κ term is the ocean heat uptake efficiency (W m −2 °C −1 ) and is based on the definition used in Raper 150 et al. (2002), where κ is the ratio between the atmosphere and ocean temperature difference that best fits observed OHC data (Sect. 2.2.8 describes the OHC data records used in our analysis). The value of κ is determined based upon the best fit (described below) between QOCEAN and the observed OHC record. The term ΔTOCEAN,HUMAN represents the temperature response of the well-mixed, top 100 m of the ocean due to the total anthropogenically driven rise in OHC. This formulation of ΔTOCEAN,HUMAN allows the model 155 ocean to warm in response to an atmospheric warming. We use a 6 year lag (72 months) for QOCEAN to account for the time needed for the energy leaving the atmosphere to heat the upper ocean and penetrate to depth, based on Schwartz (2012). Our analysis of modeled GMST is insensitive to whether this 6 year lag or the 10 year lag from Lean and Rind (2009) is used. The tSTART and tEND limits on the integral in Eq.
(5) are the start and end years, associated with each OHC record. The start and end years vary between 160 the 5 OHC records (see supplement for the different start and end years). The constant f0 term in Eq. (5) is a combination of the heat capacity of ocean water, the fraction of total ocean volume in the surface layer, and the fraction of total QOCEAN that warms the surface layer, and is equal to 8.76×10 -5 °C m 2 W −1 .
We represent the global ocean as being 1 km deep for 10% of the ocean area (representing the continental shelves) and 4 km deep for the remaining area, which approximates the average depth of the actual world's 165 oceans to within 3%; 3.7 km compared to 3.682-3.814 km from Charette and Smith (2010). Based upon our analysis of decadal ocean warming as a function of depth extracted from CMIP5 GCMs, we have determined that 13.7% of the rise in total OHC occurs in the well mixed, upper 100 m of the ocean, the term represented by ΔTOCEAN,HUMAN in equation (4). The bottom rung of Fig. 1 compares our modeled OHC to the observed OHC record based upon the average of five data sets; the value of κ resulting in the 170 best simulation of observed OHC is shown.
We use the reduced chi-squared (χ 2 ) metric to define the goodness of fit between the modeled and measured GMST anomaly for the atmosphere and also between simulated and observed OHC. Equation (6) and Eq. (7) show the calculations for χ 2 for the atmosphere, and Eq. (8) shows the calculation for χ 2 for the ocean. As noted above, minimization of the difference between the measured and modeled GMST 175 anomaly results in the EM-GC being able to replicate the observed rise in temperature over the past 170 years quite well, as shown in Fig. 1. We have added two additional new features to the model to assure accurate representation of the rise in OHC as well as the rise in GMST since 1940. The first new feature, Eq. (7), was added because of a change in the specification of the uncertainty of the GMST anomaly (σOBSi in Eq. (2)) given by the Hadley Center Climatic Research Unit (CRU). A recent update resulted in 180 much larger uncertainties being ascribed to the GMST anomaly for the entire data record, which caused some solutions to yield visually poor simulations of the rise in GMST over the past 4 to 5 decades. The second new feature, Eq. (8), was added because for some selections of the radiative forcing due to tropospheric aerosols (AER ΔRFi in Eq. (2)), the original model formulation was converging but producing simulations of OHC that seemed physically improper, based on visual inspection of observed 185 and modeled OHC. As a result of these two issues, all calculations shown here are subject to three goodness-of-fit constraints, described by Eq. (6) to (8): https://doi.org/10.5194/esd-2020-67 Preprint. Discussion started: 8 September 2020 c Author(s) 2020. CC BY 4.0 License. 190 Here, <ΔTOBS>, <ΔTMDL>, and <σOBS> in Eq. (6) and Eq. (7)  The calculation of χ 2 RECENT shown in Eq. (7) is used to constrain the model to match the observed changes in GMST over the time frame 1940-2019, a total of 80 years (NYEARS,REC equals 80). This time frame was chosen to include a full cycle of AMOC, as the strength of the thermohaline circulation tends to vary on a period of 60-80 years (Chen and Tung, 2018;Kushnir, 1994;Schlesinger and Ramankutty, 1994 rise in GMST over the past 80 years, which results in modeled GMST that replicates observed GMST for 215 the entire time series. Figure 1 shows the observed and modeled GMST anomaly from 1850-2019, and the various anthropogenic and natural components that constitute modeled GMST. Figure 1a shows the value of climate feedback, 1.38 W m −2 °C −1 , that is needed to achieve a best fit to the climate record for this simulation, resulting in values of χ 2 ATM = 0.71 and χ 2 OCEAN = 0.32. Figure 1b is the total contribution of 220 Figure 1. Measured and modeled GMST anomaly (ΔT) relative to a pre-industrial (1850-1900) baseline. (a) Observed (black) and modeled (red) ΔT from 1850-2019. This panel also displays the values of λΣ and χ 2 ATM (see text) for this best-fit simulation. (b) Contributions from total human activity. This panel also denotes the numerical value of the attributable anthropogenic warming rate from 1975-2014 (black dashed) as well as the 2σ uncertainty in the slope. (c) Solar irradiance (light blue) and major volcanoes (purple). (d) Influences from ENSO on ΔT. (e) Contributions from AMOC to ΔT and to observed warming from 1975-2014. (f) Influences from PDO (blue) and IOD (pink) on ΔT. (g) Measured (black) and modeled (red) ocean heat content (OHC) as a function of time for the average of five data sets (see text), the value of χ 2 OCEAN for this run, as well as the ocean heat uptake efficiency, κ, needed to provide the best-fit to the OHC record. The error bars (blue) denote the uncertainty in OHC used in this analysis (see Sect. 2.2.8).
human activity to variations in GMST, which includes GHGs, AER, LUC, and the export of heat from the atmosphere to the ocean. For the simulation shown, the aerosol radiative forcing is −0.9 W m −2 , the best estimate given by IPCC 2013 . This panel also notes the time rate of change of GMST attributed to humans from 1975-2014, or the attributable anthropogenic warming rate (AAWR (see Sect. 2.3)). Figure 1c illustrates the contribution to the GMST anomaly from TSI (Solar) and SAOD Sect. 2.2.8) and the modeled value of OHC (red line). For this simulation, the OHC data is best fit for a value of κ equal to 1.28 W m −2 °C −1 , which falls within the range of empirical estimates for this parameter given by Raper et al. (2002). The sum of the contributions of human activity, TSI, SAOD, ENSO, AMOC, as well as the PDO and the IOD to temporal variations in the GMST anomaly shown in Fig. 1b to 1f plus the value of C0 equals the modeled GMST anomaly, shown by the red line in Fig. 1a.

Temperature Data
We use global mean surface temperature anomalies from four data centers, the Hadley Centre Climatic Research Unit (CRU, (Morice et al., 2012)) from 1850-2019, National Centers for Environmental 240 Information (NCEI, (Smith et al., 2008;Zhang et al., 2019)) from 1880-2019, NASA Goddard Institute of Space Studies (GISS, (Hansen et al., 2010)) from 1880-2019, and Berkley Earth Group (BEG, (Rohde et al., 2013)) from 1850-2019. Our analysis primarily uses the CRU GMST data set, but in some sections, results are shown for all four data sets. All temperature anomalies are with respect to a pre-industrial baseline . To transform each data record so that the temperature anomaly is relative to a pre-245 industrial baseline, we adjust all data sets relative to the CRU baseline of 1961-1990 because we primarily use the CRU data record in this analysis. We then transform each data set to the pre-industrial baseline, as described in the methods section of Hope et al. (2017). https://doi.org/10.5194/esd-2020-67 Preprint. Discussion started: 8 September 2020 c Author(s) 2020. CC BY 4.0 License.

250
For this analysis, we use the estimates of the future abundances of greenhouse gases and aerosols provided by the SSPs. There are twenty-six scenarios, five baseline pathways and twenty-one mitigation scenarios.
The baseline pathways follow specific narratives for factors such as population, education, economic growth, and technological developments of sources of renewable energy Fricko et al., 2017;Fujimori et al., 2017;Kriegler et al., 2017;van Vuuren et al., 2017) to represent several possible 255 futures spanning different challenges for adaptation and mitigation to climate change (O'Neill et al., 2014). The twenty-one mitigation scenarios follow one of the baseline pathways but include specific climate policy to reach a designated radiative forcing at the end of the century.
As part of CMIP6, the ScenarioMIP experiment (O'Neill et al., 2016) includes eight SSPs (SSP1-1.9, SSP1-2.6, SSP4-3.4, SSP2-4.5, SSP4-6.0, SSP3-7.0, SSP5-8.5, and SSP5-3.4-OS) that GCMs use to 260 project future GMST. The first number is the reference pathway that the scenario follows (i.e. SSP1 follows the first SSP narrative) and the numbers after the dash are the target radiative forcing at the end of the century (i.e. SSP1-2.6 reaches around 2.6 W m −2 in 2100). The ScenarioMIP experiment designates Tier 1 and Tier 2 scenarios. The Tier 1 scenarios, SSP1-2.6, SSP2-4.5, SSP3-7.0, and SSP5-8.5 are considered high priority and required for all of the GCMs participating in ScenarioMIP (O'Neill et al.,265 2016). The Tier 2 scenarios are SSP1-1.9, SSP4-3.4, SSP4-6.0, and SSP5-3.4-OS (an overshoot pathway that follows SSP5-8.5 until around 2040, where carbon dioxide emissions drastically decrease and become negative in 2065). The Tier 2 scenarios are not required for modeling groups to run in order to participate in ScenarioMIP (O'Neill et al., 2016). Our analysis includes seven of the eight ScenarioMIP SSPs: all but the overshoot pathway. We highlight four in the main paper: two Tier 1 (SSP1-2.6 and 270 SSP2-4.5) and two Tier 2 (SSP1-1.9 and SSP4-3.4) scenarios. Analysis of the other three SSPs is included in supplement. Figure 2 shows the time evolution of the atmospheric abundance of the three major anthropogenic GHGs (carbon dioxide, methane, and nitrous oxide) for each of the seven SSPs we consider https://doi.org/10.5194/esd-2020-67 Preprint. Discussion started: 8 September 2020 c Author(s) 2020. CC BY 4.0 License. as well as observations of the global mean atmospheric abundance for these gases to the end of 2019 (Dlugokencky, 2020;Dlugokencky and Tans, 2020).

Greenhouse gases
The historical values of GHG mixing ratios were provided by Meinshausen et al. (2017) from 1850-2014.
We used the equations from Myhre (1998) and IPCC 2013 to calculate the change in RF due to carbon dioxide (CO2), methane (CH4), nitrous oxide (N2O), ozone depleting substances (ODS), 280 hydrofluorocarbons, perfluorocarbons, and sulfur hexafluoride relative to RF in year 1850. The radiative forcing of CH4 also includes the 15% enhancement from the increase in stratospheric water vapor due to rising atmospheric CH4 (Myhre et al., 2007). Values of GHG mixing ratios, other than ODSs, from 2015-2100 were from the SSP Database Fricko et al., 2017;Fujimori et al., 2017;Kriegler et al., 2017;Rogelj et al., 2018;van Vuuren et al., 2017) and are provided on a decadal basis. The mixing 285 ratios were interpolated onto a monthly time scale. We used the estimates of future ODS abundances provided in Table 6-4 of the 2018 Ozone Assessment Report (Carpenter et al., 2018), because the SSP database did not provide these estimates. We also include tropospheric ozone (O3 TROP ) as a GHG, because tropospheric ozone rivals N2O as the third most important anthropogenic GHG (Fig 8.15 of ). The RF due to O3 TROP from the RCPs provided by the Potsdam Institute for Climate Impact 290 Research (Meinshausen et al., 2011) is used, because the SSP database does not provide estimates. Values of RF due to O3 TROP from RCP2.6, RCP4.5, RCP6.0, and RCP8.5 are substituted in for SSP1-2.6, SSP2-  4.5, SSP4-6.0, and SSP5-8.5, respectively. We created new time series for the RF due to O3 TROP for SSP4-3.4 and SSP3-7.0 using linear combinations of RF time series from RCP2.6 and RCP8.5, with weights based upon the end of century total RF value due to all GHGs of the respective time series. Finally, the 295 RF time series for O3 TROP from RCP2.6 was also used for SSP1-1.9. Figure S2 shows the ozone RF time series used in this analysis and the supplement provides more information about the creation of the time series for the RF due to O3 TROP .

300
The value of the change in total aerosol radiative forcing in 2011 relative to pre-industrial (AER RF2011) is highly uncertain. Chapter 8 of the IPCC 2013 report gives a best estimate of AER RF2011 as −0.9 W m −2 , a likely range between −0.4 and −1.5 W m −2 , and a 5 th to 95 th percent confidence interval between −0.1 and -1.9 W m −2 . This substantial range in AER RF2011 results in a large spread in future projections of global GMST. Figure 3 shows the effect of varying the value of AER RF2011 on 305 projections of GMST in our EM-GC framework, for the same SSP4-3.4 GHG scenario. The middle panel on Figs. 3a, 3b, and 3c shows the contribution to GMST of GHGs, LUC, AER, as well as net human activities. As the value of AER RF2011 decreases and aerosols cool more strongly, the value of climate feedback (model parameter λΣ) rises, and the net contribution of human impact on GMST by the end of the century increases. Depending on which value of AER RF2011 is used, the rise in GMST by year 2100 310 for the SSP4-3.4 pathway could range from 1.3°C (Fig. 3a) to 2.6°C ( Fig. 3c) relative to pre-industrial.
Strong aerosol cooling offsets a substantial fraction of GHG-induced warming, and a large value of climate feedback (λΣ = 2.32 W m −2 °C −1 ) is needed to fit the historical climate record (Fig. 3c). In this case, future warming is large, well above the goals of the Paris Agreement by the end of the century.
Conversely, weak aerosol cooling offsets only a small fraction of GHG-induced warming, resulting in a 315 small value of climate feedback (λΣ = 0.73 W m −2 °C −1 ) needed to fit the observed GMST record (Fig.   3a). The use of any of the values of AER RF2011 in Fig. 3 can result in a very good fit to the climate record (i.e., χ 2 ATM ≤ 2, χ 2 RECENT ≤ 2, and χ 2 OCEAN ≤ 2).
We use the total aerosol RF time series provided by the SSP database for each SSP scenario. The database provides AER RF from 2005-2100, with values for all SSPs nearly identical until about 2015 320 (Riahi et al., 2017;Rogelj et al., 2018 RF2011 nearly equal to −1.0 W m −2 , which we take as the SSP-based best estimate of the change in total aerosol radiative forcing in 2011 relative to pre-industrial. Next, the single continuous time series is 330 scaled, again by a constant multiplicative factor, to match the IPCC 2013 best estimate and range of uncertainty for AER RF2011 . This procedure results in five additional time series of AER RF. Six time series of AER RF are thus created for each SSP, having values of AER RF2011 equal to −0.1, −0.4, −0.9, −1.0, −1.5, and −1.9 W m −2 . Figure S3 shows these six AER RF time series for SSP1-2.6 and SSP4-3.4. In the EM-GC framework, we further scale these six time series to create a total of 400 335 AER RF time series to fully analyze the range of AER RF2011 given by Myhre et al. (2013).

Total solar irradiance and stratospheric aerosol optical depth
We use the TSI time series provided for the CMIP6 models from 1850-2014 (Matthes et al., 2017)  (2) are differences of monthly mean values minus the long-term average (i.e., TSI anomalies). Consistent with prior studies (e.g., Lean and Rind (2008) and Foster and Rahmstorf (2011)) variations in solar irradiance due to the 11-year solar cycle have a small but noticeable effect on the EM-GC simulation of the GMST anomaly (Fig. 1c). For projections of future warming, we set the term TSIi in Eq. (2) equal to zero from the start of 2020 until 2100 (i.e., we

El Niño southern oscillation, Pacific decadal oscillation, and Indian Ocean dipole
We use the MEI.v2 (Wolter and Timlin, 1993;Zhang et al., 2019) to characterize the influence of ENSO on GMST. The MEI.v2 provides two month averages of empirical orthogonal functions of five different climatic variables from 1979 to present (Zhang et al., 2019). To have the ENSO index extend back to 365 1850, we compute differences in SST anomalies over the tropical Pacific basin as defined by the MEI.v2 from 1850-1870 using HadSST3 (Kennedy et al., 2011). Our internal computation of this surrogate for the MEI index is then appended to the MEI.ext of Wolter and Timlin (2011), which extends from 1871-1978, and the MEI.v2 index of (Zhang et al., 2019(Zhang et al., ) (1979(Zhang et al., -2019. This full time series provides a representation of ENSO that covers from 1850 to present. Consistent with prior regression-based 370 approaches (Foster and Rahmstorf, 2011;Lean and Rind, 2008), we find a significant portion of the monthly and at times annual variation in GMST is well explained by ENSO (Fig. 1d). As for the other natural terms, we assume ENSOi in Eq. (2) is zero for 2020-2100.
The Pacific decadal oscillation is the leading principal component of North Pacific monthly SST variability poleward of 20°N (Barnett et al., 1999). The PDO index maintained by the University of 375 Washington provides monthly values from 1900-2018. The PDO varies on a multidecadal time scale and affects climate in the North Pacific and North America, and has secondary effects in the tropics (Barnett et al., 1999).
The Indian Ocean dipole is based upon the difference in the anomalous sea surface temperatures (SST) between the western equatorial Indian Ocean (50°-70° E and 10° S-10° N) and the south eastern 380 equatorial Indian Ocean (90° E-110° E & 10° S-0° N) as defined in Saji et al. (1999). We use 1°  1° SSTs from the Centennial in situ Observation-Based Estimate (COBE) (Ishii et al., 2005) to create an IOD index from 1850-2019. As noted above and shown on Fig. 1f, the regression coefficients for PDO and IOD are quite small. We find little influence of either PDO or IOD in the CRU time series of GMST, but these terms are retained for completeness. We assume PDOi and IODi in Eq. (2) are zero after the start of 385 2019 and 2020, respectively. https://doi.org/10.5194/esd-2020-67 Preprint. Discussion started: 8 September 2020 c Author(s) 2020. CC BY 4.0 License.

Atlantic meridional overturning circulation
We use the Atlantic multidecadal variability (AMV) index as the area weighted, monthly mean SST from HadSST3 (Kennedy et al., 2011), between the equator and 60° N in the Atlantic Ocean (Schlesinger and 390 Ramankutty, 1994) to characterize the influence of variations in the strength of the AMOC on GMST.
The AMV index is detrended using the RF anomaly due to anthropogenic activity over the historical time frame of the analysis, as discussed in Sect. 3.2.3 of Canty et al. (2013), because this detrending option removes the influence of long-term global warming on the AMV index. The detrended AMV index serves as a proxy for variations in the strength of the AMOC (Knight et al., 2005;Medhaug and Furevik, 2011;395 Zhang and Delworth, 2007), which has particularly noticeable effects on climate in the Northern Hemisphere (Jackson et al., 2015;Kavvada et al., 2013;Nigam et al., 2011). For this analysis, the index has been Fourier filtered to remove frequencies above 9 yr −1 to retain only the low frequency, high amplitude component of the thermohaline circulation ). As noted above and shown in considerable debate about the validity of the use of a proxy such as the AMV index as a surrogate for the climatic effects of AMOC that is centered mainly around how much of the variability of the index is either internal (i.e., natural variability) or externally forced (i.e., driven by anthropogenic factors) (Haustein et al., 2019;Knight et al., 2005;Medhaug and Furevik, 2011;Stouffer et al., 2006). We stress, as explained 405 below, none of our scientific conclusions are altered if we neglect AMV as a regression variable.

Ocean heat content records
Ocean heat content data records from five recent and independent papers are used in this study. We utilize OHC data from Balmaseda et al. (2013), Carton et al. (2018), Cheng et al. (2017), Ishii et al. (2017), and 410 Levitus et al. (2012), as well as the average of the records to model the export of heat (OHE) from the atmosphere to the ocean. Figure S4 shows these five OHC records as well as the multi-measurement average. While most of these data sets have a common origin, they differ in how extensive temporal and spatial gaps in the coverage of ocean temperatures have been handled, ranging from data assimilation https://doi.org/10.5194/esd-2020-67 Preprint. Discussion started: 8 September 2020 c Author(s) 2020. CC BY 4.0 License. (Carton et al., 2018) to an iterative radius of influence mapping method (Cheng et al., 2017). The five 415 data sets have been normalized (by applying a constant additive offset) to the same value of zero in 1986, which is the midpoint of the multi-measurement time series, for visual comparison. Since OHE, in units of W m −2 , is based upon the slope of each OHC data set, this offset has no impact on the computation of OHE from OHC that is central to our study. For the computation of OHE from OHC, we use a value of the surface area of the world's oceans equal to 3.3  10 14 m 2 . The OHC records 420 we analyze are for the upper 700 m of the ocean. To calculate the OHE for the whole ocean, we multiply the OHE by 1/0.7 to account for the fact that the upper 700 m of the ocean holds 70% of the heat (Sect. 5.2.2.1 (Solomon, 2007)). When we subtract the amount of heat going into the ocean in Eq. 2 (QOCEAN), we also must account for the difference in surface area between the global atmosphere and the world's oceans. Since the QOCEAN term is computed for the surface area of the ocean, but the forcing is applied to 425 the whole atmosphere, we multiply the QOCEAN term by the ratio of the surface area of the ocean to the surface area of the atmosphere, which is 0.67.
As noted above, the calculation of χ 2 OCEAN shown in Eq. (8) is used to constrain our model representation of the temporal rise in OHC. Only model runs that provide a good fit to the observed OHC record are shown below. For these five OHC data sets, uncertainty estimates are not always provided. 430 Furthermore, some studies that do provide uncertainties give estimates that seem unreasonably small (see Fig S5 and supplement). Because of the discrepancy in uncertainties between OHC records, we create a new uncertainty time series using both the standard deviation of the mean and the uncertainties from Cheng et al. (2017). We create this new uncertainty from 1955-2019 by a monthly time step and use either the standard deviation of the mean or the uncertainties from Cheng et al. (2017), whichever is larger, for 435 that month. We use the Cheng et al. (2017) uncertainties because these estimates are the largest of the five data sets. Additionally, the standard deviation from the mean is very low in the 1980s, which is an artifact of our normalization treatment, not inherent to any of the records. This combined uncertainty estimate is substituted in for each individual data set and the average, resulting in our use of the same time varying uncertainty in OHC for all data sets. Figure  warming rate and equilibrium climate sensitivity. The attributable anthropogenic warming rate, or AAWR, is the time rate of change of GMST due to humans from 1975-2014. We use AAWR as a metric in the EM-GC framework to quantify the human 455 influence on global warming over the past few decades, and most importantly to also assess how well the CMIP6 GCMs can replicate this quantity. This analysis is motivated by the study of Foster and Rahmstorf (2011), who examined the human influence on the time rate of change of GMST from 1979-2010 using a residual method. We extend the end year of our analysis to 2014 because this is the last year of the CMIP6 Historical simulation. We pushed the start year back to 1975 so that our analysis covers a forty- 460 year period, over which the effect of human activity on GMST rose nearly linear with respect to time ( Fig. 1b).

Attributable anthropogenic warming rate
We calculate AAWR using the EM-GC based upon a linear fit to the ΔTHUMAN,ATM term shown below:  Table   S1), as well as whether or not the AMOC, PDO, or IOD terms are included in the regression framework 475 Hope et al., 2017). We are able to fit the climate record better (i.e. smaller values of χ 2 in Eqs. (6), (7), and (8)

Equilibrium Climate Sensitivity
The equilibrium climate sensitivity (ECS), which represents the warming that would occur after climate 515 has equilibrated with atmospheric CO2 at the 2×pre-industrial level (Kiehl, 2007;Otto et al., 2013;Schwartz, 2012) is also used to compare results of our EM-GC to CMIP6 multi-model output. To calculate ECS from the EM-GC, we use the following equation: That represents the rise in GMST for a doubling of CO2, assuming no other perturbations as well as 520 equilibrium in other components of the climate system (i.e., QOCEAN = 0) . The expression for the radiative forcing of CO2 is from Myhre (1998). The quantity γ in Eq. (10), which represents the sensitivity of the GMST to feedbacks within the climate system, is the only variable component of ECS. We only use values of γ that result in good fits (χ 2 ≤ 2 for Eq. (6) to (8)) between modeled and observed GMST and modeled and observed OHC. 525 For the estimate of ECS from the CMIP6 multi-model ensemble, we use the method described by Gregory et al. (2004). Near surface air temperature output from the Abrupt 4CO2 and piControl simulations, as well as net downward radiative flux output from the Abrupt 4CO2 simulation is used to calculate ECS. At the time of this analysis, 28 models released the necessary output to the CMIP6 archive (see Table S4 for the list of models and individual values of ECS of ΔT to a quadrupling of CO2. This equilibrium response is then divided by two (Jones et al., 2019) to 535 arrive at the equilibrium climate sensitivity (Fig. S6).

Aerosol Weighting Method
Probabilistic forecasts of the future rise in GMST for various SSPs are an important part of our analysis.
Probabilities of AAWR and ECS are computed by considering the uncertainty in AER RF2011. We also 540 provide probabilistic estimates of AAWR and ECS. All of these quantities are computed by incorporating the uncertainty in the radiative forcing of climate due to tropospheric aerosols within results of our EM-GC simulations. We use an asymmetric Gaussian to assign weights to the value of GMST, AAWR or ECS found for various time series of radiative forcing by aerosols associated with particular values of AER RF2011. Figure 5a shows the asymmetric Gaussian function we use to maximize the values of AAWR  how well the CMIP6 multi-model ensemble (see Table S5 for a list of CMIP6 GCMs analyzed in this study) is able to simulate the human influence on global warming over the past several decades. The EM-570 GC results in Fig. 6 have been constrained by blended near surface air temperature (TAS) and the temperature at the interface of the atmosphere and the upper boundary of the ocean (TOS) (Griffies et al., 2016 these GCMs, and calculation of AAWR for 1975-2014 as described in Sect. 2.3, we find that AAWR based upon blended CMIP6 temperature is only 3.5% lower than AAWR found when using only TAS. Tokarska et al. (2020) similarly report a much smaller influence upon the use of blended CMIP6 temperature output than suggested by Cowtan et al. (2015). Since this difference is so small and does not affect any of our conclusions, we use TAS output from the CMIP6 multi-model archive because this 585 choice allows the behavior of many more GCMs to be examined. The box and whisker symbol labeled CMIP6 in Fig. 6   represents the 25 th , 50 th , and 75 th percentiles, the whiskers denote the 5 th and 95 th percentiles, and the stars show the minimum and maximum values of AAWR from the EM-GC based upon the aerosol weighting method described in Sect. 2.5. The red box labeled "CMIP6" shows the 25 th , 50 th , and 75 th percentiles, the whiskers represent the 5 th and 95 th percentiles, and the stars denote the minimum and maximum values of AAWR from the 50 member CMIP6 multi-model ensemble.
60% larger than the 50 th percentile value of AAWR found using the CRU data set for GMST noted above.
The 5   considerably between the Third Assessment Report (TAR) and the Fourth Assessment Report (AR4). The change in temperature over time for the TAR and AR4 only span 17 and 10 years, respectively (Hausfather et al., 2020). In Fig. 6, we examine the ability of the GCMs to simulate the rise in GMST attributed to humans over a 40 year time period, which provides a better measure of how well the models simulate the observations than when using a shorter time period. The temperature change over time for the TAR and

ECS
Equilibrium climate sensitivity (ECS) is a metric often used to compare the sensitivity of warming among GCMs, as well as with warming inferred from the historical climate record. Figure 7 shows values of ECS inferred from the climate record using our EM-GC, four GMST data sets, and five OHC records (as indicated). As for AAWR, the largest variation in ECS is driven by uncertainty in AER RF2011. The Earth's climate will exhibit a larger rise in GMST to reach equilibrium than if the value of QOCEAN inferred https://doi.org/10.5194/esd-2020-67 Preprint.  Table S4 for ECS values for specific models). In contrast, the average   Meehl et al., 2020;Sherwood et al., 2020;Zelinka et al., 2020). The IPCC 2013 report gives a likely range of 1.5°C to 4.5°C for ECS (Stocker et al., 2013), yet some of the CMIP6 GCMs have values of ECS more than 1°C above this range. However, some in the climate community seem to currently 690 doubt whether the very large values of ECS are representative of the real world (Forster et al., 2020;Lewis and Curry, 2018;Tokarska et al., 2020). Gettelman et al. (2019)  portion of Fig 5c) are assigned low weights due to the assessment by Myhre et al. (2013) that the large 715 AER RF2011 associated with these ECS values is unlikely.
Four empirical determinations of ECS (our study plus Lewis and Grünwald (2018), Otto et al. (2013), and Skeie et al. (2018)) and the CMIP6-based estimate of Nijsse et al. (2020) are in slight contrast with the 2.3-4.7°C range for ECS (5 th and 95 th confidence interval) published recently by Sherwood et al.
(2020) (Fig 8). As noted above, Sherwood et al. (2020) use paleoclimate data to rule out the high range 720 of ECS. They rely on a determination that the feedback due to clouds is moderately to strongly positive to rule out the low range of ECS found by our analysis and the four studies noted above. We caution that knowledge of the cloud feedback from observations is generally limited to databases such as the International Satellite Cloud Climatology Project (ISCCP) (Schiffer and Rossow, 1983) and Pathfinder Atmospheres Extended (PATMOS-x) (Foster and Heidinger, 2013) that, while monumental in terms of 725 complexity and scope, cover only a fairly short (i.e., about 36 years) part of the century and a half climate record Sherwood et al., 2020). Most assessments of total cloud feedback rely on some combination of observations such as ISCCP, PATMOS-x, or other satellite records together with the results of regression analysis, GCM projections, and large eddy simulations that are able to resolve some of the convective processes involved in cloud formation Sherwood et al., 2020). The  Fig. 3). This fact, combined with the recent study by Weaver et al. (2020) who report no long term statistically significant trend in global cloud reflectivity at 340 nm averaged between 45° S and 45° N based on analysis of data collected by a variety of NOAA and NASA satellite instruments, causes us to suggest the true value of ECS may lie below the 2.3°C lower limit (5 th percent confidence interval) given by Sherwood et al. (2020). 740 In our model framework, the largest uncertainty in ECS is driven by imprecise knowledge of the radiative forcing of climate by tropospheric aerosols. As shown in Fig. 5c can be inferred from the century and a half long climate record. We stress that each value of ECS shown in Fig. 5c is based upon a simulation for which χ 2 ATM, χ 2 RECENT, and χ 2 OCEAN are all less than or equal to 2. Better knowledge of AER RF for the contemporary atmosphere would lead to a reduction in the 745 uncertainty of ECS. Numerous studies of the climate record, including our century and a half simulations, infer the possibility of lower values of ECS than was given by a recent analysis of studies that involve examination of data from compendiums such as ISCCP and PATMOS-x (Sherwood et al., 2020).
However, Sherwood et al. (2020) did not examine consistency of the inferred value of ECS with the ability of models to accurately simulate the GMST anomaly between 1850 and present.

CMIP6
The CMIP6 multi-model archive provides future projections of the GMST anomaly relative to preindustrial (ΔT) using the ScenarioMIP Shared Socioeconomic Pathways (SSPs). Figure 9 shows the 755 CMIP6 multi-model ensemble projections of ΔT for the four SSPs (SSP1-1.9, SSP1-2.6, SSP4-3.4, and SSP2-4.5) highlighted in our analysis. Each SSP scenario has varying amounts of gridded, monthly mean TAS projections submitted to the CMIP6 archive by GCMs (indicated on each plot). Global, monthly ΔT was created by averaging the TAS output over the globe with a cosine latitude weighting. The global, monthly ΔT time series for all of the runs for each CMIP6 GCM were averaged together to obtain one 760 time series of ΔT. The varying amount of GCM output available for each SSP scenario is due to the fact that: a) SSP1-2.6 and SSP2-4.5 are Tier 1 scenarios (O'Neill et al., 2016) and are designated as priority over the other SSPs (as described in Sect. 2.2.2), and b) not all GCMs have provided results to the CMIP6 archive at the time of the analysis. More CMIP6 multi-model output will likely become available as modeling groups who have not submitted output to the CMIP6 archive finalize their results. However, we 765 do not expect additional GCM simulations will affect our conclusions unless the GCM output is significantly different than that currently available.
https://doi.org/10.5194/esd-2020-67 Preprint. Discussion started: 8 September 2020 c Author(s) 2020. CC BY 4.0 License. Fig. 9 labeled as the IPCC 2013 likely range is the same trapezoid as that displayed on Figure 11.25b from chapter 11 of the IPCC 2013 report (Kirtman et al., 2013). All of the  GCMs having future values of ΔT that are considerably higher than others. This divergence for GCM projections of ΔT is especially evident in Fig. 9a, c, and d. The two CMIP6 GCMs that have the highest values of ΔT across the four SSPs are CanESM5 (Swart et al., 2019) and UKESM1 (Sellar et al., 2020). 785 The CanESM5 and UKESM1 GCMs have the highest values of AAWR (0.349°C/decade and 0.294°C/decade, respectively), large values of ECS (5.72°C and 5.44°C, respectively), and exceed observed ΔT reported by CRU for the past few decades (apparent in Fig. 9).

790
The EM-GC is also used to project future changes in ΔT using the SSPs. Figure 10 shows the GMST anomaly in 2100 from pre-industrial (ΔT2100) as a function of the climate feedback parameter and AER RF2011, for the four SSPs highlighted throughout. Only model runs from the EM-GC that achieved a good fit to the climate record (χ 2 ATM ≤ 2, χ 2 RECENT ≤ 2, χ 2 OCEAN ≤ 2) are shown. The EM-GC runs that satisfy these three χ 2 constraints but fall outside of the IPCC 2013 range for AER RF2011  are 795 shaded grey (left hand side of each panel). We do not consider the EM-GC projections that lie outside of the IPCC 2013 range for AER RF2011 in our projections of ΔT, yet these results are shown to illustrate that the EM-GC can fit the climate record with estimates of the RF due to tropospheric aerosols that lie below (i.e., less cooling) of the 5 th and 95 th confidence interval of −0.1 to −1.9 W m −2 for AER RF2011 given by IPCC 2013. We cannot establish any good fits of the climate record for AER RF2011 with a 800 cooling stronger than about −1.6 W m −2 . The range of ΔT2100 we compute using the EM-GC for SSP1-1.9, SSP1-2.6, SSP4-3.4, and SSP2-4.5 are 0.65-2.16°C, 0.82-2.78°C, 1.00-3.28°C, and 1.21-3.78°C, https://doi.org/10.5194/esd-2020-67 Preprint. Discussion started: 8 September 2020 c Author(s) 2020. CC BY 4.0 License. respectively. Results for SSP4-6.0, SSP3-7.0, and SSP5-8.5 are shown in Fig. S8: ΔT2100 ranges are 1.41-4.47°C, 1.84-5.56°C, and 2.13-6.75°C for these three scenarios. 805 The large range of T2100 found for any given SSP scenario (i.e., a factor of 3.1 difference between the smallest and largest end of century warming for SSP2-4.5) is caused by the fact that the climate record can be fit nearly equally well by a considerably large combination of the climate feedback parameter (our ) and scenarios for radiative forcing due to tropospheric aerosols. The more aerosols have cooled, offsetting the relatively well-known warming due to GHGs, the larger  must be to fit the climate record. 810 Since the RF of aerosols will diminish in the future due largely to public health concerns (Lelieveld et al., 2015;Shindell et al., 2016;Smith and Bond, 2014), the part of our model ensemble requiring relatively

Comparing CMIP6 and EM-GC
Time series of future projections of ΔT from the EM-GC can be illustrated as probabilistic forecasts. Figure 11 shows the change in future ΔT for SSP1-1.9, SSP1-2.6, SSP4-3.4, and SSP2-4.5 colored by the probability of reaching at least that rise in ΔT by the end of the century. The EM-GC probabilities are 820 computed from ensemble members for model runs constrained by the CRU data records for GMST and the average of 5 OHC data records (Fig. S4) based upon the aerosol weighting method, described in Sect.
2.5. The trapezoid from chapter 11 of IPCC 2013 (Kirtman et al., 2013) is shown on Fig. 11 in black to highlight that the EM-GC projections of the future rise in ΔT lie within the IPCC 2013 likely range of warming. The Paris Agreement target and upper limit are included to compare the EM-GC projections of 825 future ΔT to the Paris Agreement goals. The white shaded region is the EM-GC's median estimate of future ΔT for each SSP scenario. The median estimate for ΔT2100 for simulations using SSP1-1.9 and SSP1-2.6 falls below the Paris Agreement target at 1.0°C and 1.3°C, respectively. The median estimate of ΔT2100 from the EM-GC for SSP4-3.4 is between the Paris Agreement target and upper limit at 1.6°C.
For SSP2-4.5 the median estimate of ΔT2100 is just below the Paris Agreement upper limit at 1.9°C. The 830 CMIP6 minimum, multi-model mean, and maximum projections of ΔT, identical to those in Fig. 9, are also shown in Fig. 11. The CMIP6 minimum projection of the rise in ΔT falls near the EM-GC median estimate of ΔT for each SSP scenario. The CMIP6 multi-model mean value of the future change in ΔT falls below the EM-GC maximum value of ΔT, while the CMIP6 maximum value is far above the maximum projections of the future rise in ΔT using the EM-GC. Results for SSP4-6.0, SSP3-7.0, and 835 SSP5-8.5 are provided in Fig. S9.
https://doi.org/10.5194/esd-2020-67 Preprint. Discussion started: 8 September 2020 c Author(s) 2020. CC BY 4.0 License. Figure 12 compares probability distribution functions (PDFs) for the projection of ΔT2100 for the EM-GC and the CMIP6 multi-model ensemble. For the CMIP6 multi-model results, we compute the probabilities of achieving the Paris Agreement target of 1.5°C and upper limit of 2.0°C (by the end of the 840 century) by calculating how many of the GCMs participating in each scenario have projections of ΔT2100 below the target or upper limit. In contrast, the probabilities for the projections of ΔT2100 using our EM-GC are computed using the aerosol weighting method, described in Sect. 2.5. The height of each histogram represents the probability that a particular range of ΔT2100, defined by the width of each line segment, will occur. The left-hand y-axis displays the probability of ΔT2100 using the EM-GC, while the 845 right-hand y-axis represents the probability of ΔT2100 using the CMIP6 multi-model simulations. The values on the CMIP6 multi-model ensemble y-axis are double the values on the EM-GC y-axis, for visual Numerical values of probabilities for staying at or below the Paris Agreement target or upper limit for all seven SSP scenarios are given in Table 1. The CMIP6 multi-model projections exhibit lower probabilities of achieving the goals of the Paris Agreement than the projections using the EM-GC. In the creation of ScenarioMIP, SSP1-2.6 was designed to be a scenario that achieved the Paris Agreement goals and likely (greater than 66% probability (Stocker et al., 2013)) limited warming below 2.0°C, and was 855 expected to produce a future rise in ΔT2100 of 1.7°C (O'Neill et al., 2016). The CMIP6 multi-model probability of SSP1-2.6 to stay at or below 2.0°C is 51.5%, as shown in Table 1. Based on our analysis, Figure 12. Probability density functions (PDF) for ΔT2100 found using the EM-GC and CMIP6 multi-model results.
the CMIP6 multi-model ensemble does not indicate SSP1-2.6 as being a 2.0°C pathway, because it will only provide about a 50:50 likelihood of limiting warming below 2.0°C. 860 Table 1. List of SSP scenarios analyzed in this study and the probabilities of achieving the Paris Agreement target or upper limit based on the EM-GC and the CMIP6 multi-model ensemble. The probabilities using the EM-GC are computed using the aerosol weighting method. The probabilities using the CMIP6 models are computed by calculating how many of the models for that scenario are below the temperature limits compared to the total number of models. Projections of ΔT2100 based on the EM-GC provide more optimism for achieving the Paris Agreement goals than the CMIP6 multi-model ensemble. The SSP1-1.9 scenario results in an 84.1% probability of ΔT2100 staying at or below 1.5°C, while SSP1-2.6 gives a 64.8% likelihood of global warming staying at or below 1.5°C by end of century (Table 1). The SSP1-1.9 scenario involves extreme 870 climate mitigation that is unlikely to happen in the next few years with atmospheric CO2 peaking close to present day values (Fig. 2a). The SSP1-2.6 scenario requires less climate mitigation than SSP1-1.9 (though still requires net negative emissions towards the end of the century) and provides greater than a 50% likelihood of staying at or below the Paris Agreement target, thus we designate SSP1-2.6 as the 1.5°C pathway in our model framework instead of SSP1-1.9. This result is supported by Tokarska et al. 875 (2020), who show that the SSP1-2.6 scenario has a likely range of warming at 1.33-1.99°C above https://doi.org/10.5194/esd-2020-67 Preprint. preindustrial by end of century, based upon filtering CMIP6 GCM output on the level of agreement with the observed climate record. Previous studies suggested that a 2.6 W m −2 scenario was in line with the 2.0°C goal O'Neill et al., 2016;Riahi et al., 2015). However, our analysis suggests the 2.6 W m −2 scenario provides about a 66% probability of achieving the more stringent 1.5°C 880 target, and that a 3.4 W m −2 scenario (i.e., ) is in line with the 2.0°C goal and has about a 74% probability of limiting warming to 2.0°C (Table 1). We therefore designate SSP4-3.4 as the 2.0°C pathway. Significant climate mitigation efforts will be required to keep the growth of CO2, CH4, and N2O below the trajectories shown for SSP1-2.6 (1.5°C pathway in our model framework) and SSP4-3.4 (2.0°C pathway) (Fig. 2). 885

Transient climate response and carbon budgets
The transient climate response to cumulative emissions (TCRE) relates the rise in ΔT to the cumulative amount of carbon released to the atmosphere by human activities. We illustrate TCRE from the EM-GC as probabilistic forecasts, as shown in Fig. 13, to analyze future projections of ΔT. Figure  intersection of these horizontal lines with the dark grey and the black curves are the 66% and 50% probabilities, respectively of the Paris Agreement target or upper limit being attained. The SSP5-8.5 https://doi.org/10.5194/esd-2020-67 Preprint.  Carbon Budget project (Friedlingstein et al., 2019) Paris Agreement thresholds. The range of years given represents when the Paris Agreement thresholds will be passed based upon the rate of emissions from SSP5-8.5 or continuing the current rate of emissions of 11.7 Gt C yr −1 .  (2051)(2052)(2053)(2054)(2055)(2056)(2057)(2058)(2059)(2060)(2061)(2062)(2063)(2064)(2065)(2066)(2067) 685 Gt C (2056-2077) a Year the 1.5°C target or 2.0°C upper limit will be exceeded assuming the rate of emission inferred from SSP5-8.5 b Year the 1.5°C target or 2.0°C upper limit will be exceeded assuming current rate of emission of 11.7 G C yr -1 An analysis by van Vuuren et al. (2020)  achieving the Paris Agreement target of limiting the rise in ΔT below 1.5°C in 2100. They base this estimate on an analysis of climate sensitivity and carbon cycle components, including an adjustment to TCRE for the tendency of CMIP5 GCMs to warm too quickly that had been suggested by Millar et al. (2017). In our model framework, we find a 66% probability of limiting warming to 1.5°C upon the release 935 of 369 Gt C between 2010 and 2100. It is not surprising our analysis provides somewhat more latitude for the probabilistic forecasts of limiting warming to 1.5°C compared to estimates based on analyses of GCM output, given the tendency of CMIP5 GCMs (Hope et al., 2017) and CMIP6 GCMs (Sect. 3.1) to warm so much faster than the observed climate system. Regardless, between 2010 and 2019, about 101 Gt C has been released to the atmosphere (Friedlingstein et al., 2019), so the remaining budget after 2019 940 for limiting warming to 1.5°C is about 127 Gt C according to van Vuuren et al. (2020). At the present pace of emissions of 11.7 Gt C yr −1 , society will cross this threshold in about a decade. Our model framework suggests a remaining budget of 268 Gt C (Table 2). Society has at most about 20 years to severely limit carbon emissions to have a 66% probability to achieve the target of the Paris Agreement.

Blended methane
Atmospheric abundances of methane will likely continue to increase as society expands natural gas production and agriculture, making it important to analyze the impact of various methane scenarios on the rise of GMST. It is unlikely future atmospheric methane abundances will progress as indicated by SSP1-2.6 (see Fig. 2), a low radiative forcing scenario. Current observations shown in Fig. 2 illustrate 950 that the methane mixing ratio is following SSP2-4.5 and has missed the initial decline needed to follow the SSP1-2.6 pathway. To analyze the effect varying future methane abundance pathways will have on GMST, we have generated linear interpolations of the SSP1-2.6 and SSP3-7.0 methane future abundances and created four alternate scenarios (see Fig. S11), which we call blended methane scenarios. We can substitute one of the blended methane scenarios into the EM-GC instead of using the projection of 955 methane specified by the SSP database to quantify the sensitivity of future warming to various evolutions of methane on the rise in GMST. Figure 14 shows the probability of staying at or below the Paris Agreement target (gold colors) or upper limit (purple colors) for SSP1-2.6 (solid) and SSP4-3.4 (dotted) as a function of the methane mixing https://doi.org/10.5194/esd-2020-67 Preprint. Discussion started: 8 September 2020 c Author(s) 2020. CC BY 4.0 License. ratio in 2100. The lowest atmospheric methane mixing ratio value in 2100 of 1.15 ppm is from the SSP1-960 2.6 methane pathway, the highest mixing ratio in 2100 of 3.20 ppm is from the SSP3-7.0 methane pathway, and the four in between are the blended methane scenarios. As the atmospheric methane abundance increases, the likelihood of achieving the goals in the Paris Agreement decreases. For SSP1-2.6, the probability of limiting the rise in GMST below the 1.5°C target begins at 65% using the SSP1-2.6 designated methane pathway and decreases as the blended scenarios are considered. The probability 965 of achieving the Paris Agreement target declines to just under 50% if methane reaches 2.4 ppm in 2100 and to 34% if methane increases to 3.2 ppm in 2100. Even though we have labeled SSP1-2.6 the 1.5°C pathway in our model framework, limiting future warming to this challenging amount can likely only be achieved by strict limits on both emissions of carbon and methane.
In Sect. 3.3.3, we showed that if all GHGs follow the SSP4-3.4 scenario there would be a 74% 970 probability of limiting warming to 2.0°C. If the methane pathway instead follows SSP1-2.6, which has an Figure 14. Probability of staying at or below the Paris Agreement target and upper limit for SSP1-2.6 and SSP4-3.4 as a function of varying methane scenarios using the EM-GC. The atmospheric methane scenarios are calculated using linear combinations of methane abundances from SSP1-2.6 and SSP3-7.0 to span the range of future methane abundances.
https://doi.org/10.5194/esd-2020-67 Preprint. Discussion started: 8 September 2020 c Author(s) 2020. CC BY 4.0 License. end of century mixing ratio of only 1.15 ppm, then the probability of achieving the Paris Agreement goal rises to 82%. However, if the methane pathway follows SSP3-7.0 and the end of century mixing ratio increases to 3.2 ppm, then the probability of achieving the Paris Agreement goal declines to 65%.
Reducing the future anthropogenic emissions of methane might be more challenging then 975 controlling future emissions of carbon dioxide, simply because methane has such a wide variety of sources related to energy, agriculture, and ruminants (Kirschke et al., 2013). Given the current widespread use of methane as a source of energy in the United States and parts of Europe (Saunois et al., 2020), combined with the continued growth in the global number of ruminants (Wolf et al., 2017), it seems unrealistic for atmospheric methane to follow the peak and sharp decline starting in 2025 of the SSP1-2.6 pathway ( Fig.   980 3b). Our analysis suggests failure to limit methane to the SSP1-2.6 trajectory will have a larger impact on the achievement of the 1.5°C Paris goal compared to the 2.0°C upper limit. Figure 14 is designed to provide some perspective on the importance of future controls on limiting the growth of methane on projections of end of century warming.

Conclusions
In this paper we use a multiple linear regression energy balance model (EM-GC), to analyze and project changes in the future rise in global mean surface temperature (GMST), calculate the attributable anthropogenic warming rate (AAWR, the component of the rise in GMST caused by human activities) over the past four decades, and compute the equilibrium climate sensitivity (ECS, the rise in GMST that 990 would occur after climate has equilibrated with atmospheric CO2 at the 2×pre-industrial level).
Projections of the rise in GMST (T) are conducted for seven of the Shared Socioeconomic Pathway (SSP) projections of GHGs (O'Neill et al., 2017). We compare computations of AAWR, ECS, and projections of T to values for each quantity computed from archived output provided by GCMs as part of CMIP6 . A critical component of our study is comprehensive analysis of 995 uncertainties in AAWR, ECS, and projections of T in our EM-GC framework, due to the rather large uncertainty in radiative forcing of climate from tropospheric aerosols (AER RF).
The best estimate of values of AAWR from 1975-2014 computed using our EM-GC constrained by the century and a half long record for GMST provided by the CRU data record (Morice et al., 2012) https://doi.org/10.5194/esd-2020-67 Preprint. Discussion started: 8 September 2020 c Author(s) 2020. CC BY 4.0 License. is 0.135°C/decade and the 5 th , 25 th , 75 th , and 95 th percentiles are 0.097, 0.115, 0.160, and 0.195°C/decade, 1000 respectively. The median value of AAWR from the CMIP6 multi-model ensemble is 0.217°C/decade and the 5 th , 25 th , 75 th , and 95 th percentiles are 0.150, 0.191, 0.242, and 0.294°C/decade, respectively. We show that the component of GMST attributed to human activity within the CMIP6 multi-model ensemble warms considerable faster than observations over the past four decades, a result that is consistent with a recent analysis of output from the CMIP6 multi-model ensemble  as well as output 1005 from CMIP5 GCMs assessed in AR5 (i.e, Fig. 11.25b of Kirtman et al. (2013)). This finding differs from the conclusion of Hausfather et al. (2020), who showed fairly good agreement between projections of global warming from GCMs and observed T. As detailed in Sect. 3.1, this paper examined GCMs that proceeded CMIP5 and examined T for a time period that ends in 2017, a time when global temperature was temporarily warming by a strong ENSO event. The majority of the uncertainty in our EM-GC based 1010 estimate of AAWR is due to imprecise knowledge of the true value of AER RF.
In our model framework the best estimate of ECS is 2.01°C and the 5 th , 25 th , 75 th , and 95 th percentiles are 1.12, 1.49, 2.50, and 4.12°C, respectively. The median value of ECS from the CMIP6 multi-model ensemble is 3.74°C, which is almost double the value of ECS inferred from the observed climate record. The 5 th , 25 th , 75 th , and 95 th percentiles of ECS from the CMIP6 multi-model ensemble are 1015 2.19, 2.83, 4.93, and 5.65°C, respectively. We obtain a wide range of ECS values using the EM-GC because of the uncertainty in AER RF. With an AER RF2011 equal to −1.6 W m −2 , the EM-GC calculates a value of ECS similar to the maximum value of ECS from the CMIP6 multi-model mean. We cannot rule out the very high value of ECS, but we assign a low probability based on the IPCC 2013 low likelihood for the needed value of AER RF2011. Our empirically based determination of ECS is in good 1020 overall agreement with the recent empirical determinations of Lewis and Grünwald (2018) (1.87°C, range of 1.1-4.05°C) and Skeie et al. (2018) (2.0°C, range of 1.2-3.1°C) and the slightly older empirically determination reported by Otto et al. (2013) (2.0°C, range of 1.2-3.9°C) (all range values are for the 5 th and 95 th percent confidence interval). Nijsse et al. (2020) reported a median value of ECS of 2.6°C (range of 1.52-4.05°C) based upon analysis of CMIP6 output tied to the actual climate record. A recent review the range of 2.3 to 4.7°C at the 5 th to 95 th percent confidence intervals; their lower bound for ECS is quite a bit higher than the lower bound found in our analysis, as well as by Grünwald (2018), Nijsse et al. (2020), Otto et al. (2013), and Skeie et al. (2018).
We also examined the probability of limiting the future rise in GMST below the Paris Agreement 1030 target of 1.5°C and upper limit of 2.0°C. Our probabilistic forecasts of projections of T include a comprehensive treatment of the uncertainty in AER RF, a capability outside the scope of the GCM intercomparisons conducted for CMIP6. Our analysis indicates that the SSP1-2.6 scenario is the 1.5°C pathway, providing a 64.8% likelihood of keeping the end of century rise in T below the Paris Agreement target of 1.5°C (relative to pre-industrial). We find that the SSP4-3.4 scenario is the 2.0°C 1035 pathway, as this scenario provides a 74.0% likelihood of limiting global warming to below the Paris Agreement upper limit of 2.0°C by end of century. In contrast, the CMIP6 multi-model mean only suggests a 15.2% probability of achieving the Paris Agreement target for SSP1-2.6 and a 16.7% probability of attaining the Paris Agreement goal for SSP4-3.4. This result is not surprising, given the tendency of most CMIP6 GCMs to warm faster than has been observed over the past four decades. Our 1040 projections of T using a physically based model tied to observations of ocean heat content, quantification of natural as well as anthropogenic drivers of variations in GMST, and consideration of uncertainty in AER RF are shown to be remarkably similar to the expert assessment of the future rise in GMST that was sketched out in Fig. 11.25b of AR5 (Kirtman et al., 2013), and the emprically-based filtering of CMIP6 model output recently published by Tokarska et al. (2020). 1045 We also quantify the sensitivity of the probability of achieving the Paris Agreement target (1.5°C) or upper limit (2.0°C) to future atmospheric abundances of methane. The end of century mixing ratio of methane in the SSP1-2.6 scenario is 1.15 ppm, considerably less than the contemporary abundance of 1.88 ppm. The likelihood of attaining the 1.5°C target for SSP1-2.6 decreases as future methane emissions increase, declines to just under 50% if methane reaches 2.4 ppm in 2100 and to 34% if methane increases 1050 to 3.2 ppm at end of century. Our analysis described in Sect. 3.3.5 demonstrates that major near-term limits on the future growth of methane are especially important for achievement of the 1.5°C limit to future warming that constitutes the goal of the Paris Agreement. https://doi.org/10.5194/esd-2020-67 Preprint. Discussion started: 8 September 2020 c Author(s) 2020. CC BY 4.0 License.
Finally, we have also quantified in the EM-GC framework the remaining budgets of carbon (i.e, CO2) emissions that can occur while attaining either the goal or upper limit of the Paris Agreement. We 1055 find that after 2019, society can only emit another 108, 268, or 336 Gt C to have either a 95%, 66%, or 50% chance of limiting warming to 1.5°C. These future emissions estimates rise to 295, 565, and 685 Gt C to have a 95%, 66%, or 50% chance of limiting warming to 2.0°C. Given that current anthropogenic emissions of carbon due to combustion of fossil fuels, cement production, gas flaring, and land use change are about 11.7 Gt per year (Friedlingstein et al., 2019), our study indicates that the target (1.5°C warming) 1060 of the Paris Agreement will be achieved unless carbon emissions are severely curtailed in the next two decades.
We concluded by noting that the CMIP6 multi-model ensemble provides many useful parameters such as sea level rise, sea ice decline, and precipitation changes, that provide a great societal understanding of the impact of climate change. We do not mean to undermine the importance of the 1065 CMIP6 GCMs by this analysis. Rather, we hope that studies such as this, along with other recent evaluations of CMIP6 multi-model output such as Nijsse et al. (2020) and Tokarska et al. (2020) will provide improved use of the CMIP6 multi-model ensemble for policy decisions. Our EM-GC was built to specifically simulate and project changes in GMST; we do not examine numerous other components of the climate system that affect society. Our study indicates that unless society can implement steep 1070 reductions in the emissions of carbon (CO2) and methane (CH4) rather soon, the Paris Agreement will fail to be achieved. We suggest there is slighly more time to achieve these steep reductions than indicated by a literal interpretation of the CMIP6 multi-model mean. The consequences for society of 1.5°C warming, 2.0°C warming, or even larger rises in GMST rely entirely on the incredibly valuable output of the CMIP6 GCMs.

Competing Interests
The authors declare that they have no conflict of interest.

Acknowledgements
We would like to acknowledge the World Climate Research Programme for coordinating and promoting CMIP6 through its Working Group on Coupled Modelling. We thank the climate modeling groups participating in CMIP6 for producing and making their model results available, the Earth System Grid Federation (ESGF) for archiving the data and providing access, and the several funding agencies 1135 who support ESGF and CMIP6. This project could not occur without the results from CMIP6. We appreciate very much financial support from the NASA Climate Indicators and Data Products for Future https://doi.org/10.5194/esd-2020-67 Preprint. Discussion started: 8 September 2020 c Author(s) 2020. CC BY 4.0 License.