fair-calibrate v1.4.1: calibration, constraining, and validation of the FaIR simple climate model for reliable future climate projections

Simple climate models (also known as emulators) have re-emerged as critical tools for the analysis of climate policy. Emulators are efficient and highly parameterised, where the parameters are tunable to produce a diversity of global mean surface temperature (GMST) response pathways to a given emission scenario. Only a small fraction of possible parameter combinations will produce historically consistent climate hindcasts, a necessary condition for trust in future projections. Alongside historical GMST, additional observed (e.g. ocean heat content) and emergent climate metrics (such as the equilibrium climate sensitivity) can be used as constraints upon the parameter sets used for climate projections. This paper describes a multi-variable constraining package for the Finite-amplitude Impulse Response (FaIR) simple climate model (FaIR versions 2.1.0 onwards) using a Bayesian framework. The steps are, first, to generate prior distributions of parameters for FaIR based on the Coupled Model Intercomparison Project (CMIP6) Earth system models or Intergovernmental Panel on Climate Change (IPCC)-assessed ranges; second, to generate a large Monte Carlo prior ensemble of parameters to run FaIR with; and, third, to produce a posterior set of parameters constrained on several observable and assessed climate metrics. Different calibrations can be produced for different emission datasets or observed climate constraints, allowing version-controlled and continually updated calibrations to be produced. We show that two very different future projections to a given emission scenario can be obtained using emissions from the IPCC Sixth Assessment Report (AR6) (fair-calibrate v1.4.0) and from updated emission datasets through 2022 (fair-calibrate v1.4.1) for similar climate constraints in both cases. fair-calibrate can be reconfigured for different source emission datasets or target climate distributions, and new versions will be produced upon availability of new climate system data.

1 Introduction

Simple climate models (also known as emulators) are designed to replicate the large-scale behaviour of more complex Earth system models. Emulators can be statistically based, such as Gaussian process emulators, or physically based, where the equations of the model can be written analytically, and relationships are based on physical understanding, where possible. The Finite-amplitude Impulse Response (FaIR) model (Millar et al., 2017; Smith et al., 2018; Leach et al., 2021) and many other reduced complexity climate models (Nicholls et al., 2020, 2021) are of the latter type. Emulators project mean temperatures for the whole globe or a few aggregated regions on a monthly or annual time step, rather than replicating a full 3D atmosphere and ocean at sub-hourly time steps such as in Earth system models (ESMs). What emulators lack in spatial, temporal, and physical detail is made up for in efficiency and flexibility. Some emulators may only report global mean surface temperature (GMST) as a climatic output. However, several regional climate variables (Mathison et al., 2024; Wells et al., 2023) and climate impacts (Shiogama et al., 2022) are shown to scale with GMST, and GMST is often used as a proxy for impacts and damages in climate policy discussions (e.g. the 1.5 and 2 °C warming levels of the Paris Agreement) and economic models (Howard and Sterner, 2017). Emulators are efficient and may run at tens, hundreds, or thousands of model years per wall clock second, compared to the model years per wall clock day yardstick for Earth system models. Simple climate models are also flexible and highly parameterised, meaning that a wide range of climate behaviour can be explored by varying parameter choices.

These two features of efficiency and flexibility make it possible to run large probabilistic ensembles using emulators to explore the range of climate uncertainty to a given emission scenario. While a number of ESMs exist, allowing us to explore differences in model responses to forcing, their relatively small number represent an ensemble of opportunity (Tebaldi and Knutti, 2007), meaning that projections using ESMs alone likely under-explore the uncertainty space. It has also been well-publicised that several Coupled Model Intercomparison Project (CMIP6) models have equilibrium climate sensitivity (ECS) outside of the very likely (nominal 5 %–95 %) range assessed by the Intergovernmental Panel on Climate Change (IPCC) Sixth Assessment Report (AR6) (Forster et al., 2021), with other expert assessments coming to similar conclusions about the range of ECS (Sherwood et al., 2020). Many CMIP6 models show a poor reconstruction of historical temperatures (Smith and Forster, 2021), with future climate projections run with only a small number of Shared Socioeconomic Pathway (SSP) scenarios (O'Neill et al., 2016) that start in 2015. These simulations are therefore rapidly becoming outdated, which means that unadjusted GMST projections from CMIP6 models are often not appropriate for understanding climate change responses to anthropogenic emissions and assessing impacts of climate policy, particularly on the short timescales that policymakers need.

Flexibility can be a double-edged sword. Emulators are only useful if the climate projections they provide are reliable. It is therefore critical that emulators are calibrated to reproduce, at the very least, the time series of historical GMST to a satisfactory standard. The IPCC AR6 Working Group 1 (WG1) provided a rigorous calibration of four emulators (MAGICC v7.5.3, FaIR v1.6.2, CICERO-SCM, and OSCAR v3.1.1) against historical observations of GMST and ocean heat content (OHC) change and IPCC-assessed distributions of ECS, transient climate response (TCR), transient climate response to cumulative CO₂ emissions (TCRE), present-day aerosol forcing, and future projections of warming under SSP scenarios, including their uncertainties. Three of the emulators, including FaIR, were assessed to be suitable to be taken forward for use by the IPCC AR6 Working Group 3 (WG3) to produce warming projections from emission pathways derived from integrated assessment models (IAMs) (Riahi et al., 2022). Over 1800 scenarios were assessed by WG3, rendering this task impossible for ESMs and necessitating the existence of reliable, well-calibrated emulators.

In this paper, we develop and formalise the calibration code for FaIR, developed originally as part of the IPCC AR6 WG1–WG3 handshake over the course of 2021 and 2022 (Kikstra et al., 2022). The fair-calibrate package is available as an open-source Python and R library that builds upon the IPCC AR6 WG1 calibration process for the FaIR model and is designed to work with FaIR model versions starting at v2.1.0, with a future backport to v2.0.0 planned. The versions of fair-calibrate described in this paper are run with FaIR v2.1.3. fair-calibrate is designed to be flexible, easy to update, and has a clearly defined version control strategy. We aim to provide updated constrained probabilistic projections of near-term and 21st century warming using FaIR at least annually to coincide with the Indicators of Global Climate Change (IGCC) project (Forster et al., 2023) as new emissions and data for updating observational constraints become available. The headline calibration version in this paper, v1.4.1, is the first example of this, with emissions and observational constraints updated through 2022. For comparison, we also provide an updated IPCC AR6 calibration (v1.4.0), using historical emissions up to 2014 and projections thereafter, showing the significant impact of using different historical emission datasets for projections.

Section 2 discusses the code requirements and version control strategy. Section 3 describes the process chain for calibrating FaIR, focusing on fair-calibrate v1.4.1. Section 4 shows the results of the calibrations v1.4.1 and v1.4.0 compared to IPCC-assessed climate indicators and their updates. Section 5 concludes.

2 Calibration requirements, versions, and versioning strategy

2.1 Requirements and reproduction

fair-calibrate is a collection of Python and R scripts and is developed on GitHub, with each version's source code, intermediate data, and final output released with digital object identifiers (DOIs) on Zenodo (Smith, 2024). Required dependencies are Python version 3.8 or later and R≥4.1.1. The fair-calibrate package requirements are managed through the Anaconda Python and R package manager, which is also required. fair-calibrate sits independently of the FaIR source code, which is deliberately kept clean.

Each calibration release contains one or more comma-separated value (CSV) files of parameters and model configuration settings that allow for the reproducibility of the calibration of any emission scenario run in FaIR and a larger ZIP file containing all results, source files, and intermediate output data produced by the calibration code so that users can inspect and quickly perform their own analysis on the prior ensemble generated without having to re-run the calibration. The ZIP files also contain diagnostic plots generated by the code, many of which are included in this paper. Intermediate output files and plots are not part of the GitHub repository, owing to their file sizes.

2.2 Version control strategy

fair-calibrate does not strictly adhere to semantic versioning, but sequential version control allows for exact reproducibility and easy comparison of calibrations. As with semantic versioning, the version string is of the form vX.Y.Z. Any change in calibration strategy that represents a departure from previous logic would increment the major version X, congruent with a “breaking change” in semantic versioning parlance. If an update to an existing calibration or constraining process would change previously submitted results if they were to be re-run with the same emissions and constraints, then this is a minor version Y increment. Examples of minor version updates include bug fixes and changes in some of the prior distribution ranges used for sampling (Sect. 3.2). The micro-version Z pertains to either the constraint set or the historical emission data used. This allows different sets of emissions or constraints to be run with the same overall calibration strategy for easy comparison. Unlike in semantic versioning, an increment of Z does not necessarily imply a bug fix or that a more recent version is in some way superior than an older version or any parallels in the Z value between different vX.Y since calibrations are developed and released whenever a new use case arises. It is not always possible for different Z micro-versions to be exactly directly comparable, but the overall sentiment should be to change as little as possible, other than emissions and/or constraints.

2.3 Calibration versions in the v1.4 series

The most recent minor version 1.4 is the focus of this paper. While the methods and results presented here are specific to v1.4, this paper is designed to serve as an overall reference to the fair-calibrate method and is intended to be a valid guidance document for many future versions.

2.3.1 v1.4.1: best-estimate historical emissions 1750–2022

fair-calibrate v1.4.1 uses up-to-date historical emissions as far as possible, and the emissions are as follows:

• CO₂ emissions for fossil fuel and industrial (FFI) and agriculture, forestry, and other land use (AFOLU) CO₂ are from the Global Carbon Project 2023 v1.0 (Friedlingstein et al., 2023).

• CH₄ and N₂O from non-biomass-burning sources, plus SF₆, NF₃, and aggregated hydrofluorocarbons (HFCs) and perfluorocarbons (PFCs), are from PRIMAP-Hist v2.5 (Gütschow and Pflüger, 2023; Gütschow et al., 2016), prioritising third-party (TP) data sources over country reported emissions.

• Short-lived climate forcers, comprising black carbon (BC), organic carbon (OC), sulfur dioxide (SO₂), nitrogen oxides (NO_x), ammonia (NH₃), carbon monoxide (CO), and volatile organic compounds (VOCs) from fossil, industrial, and agricultural sources, are from the Community Emissions Data System (CEDS) v2021.04.06 (O'Rourke et al., 2021; Hoesly et al., 2018).

• Biomass burning emissions of CH₄, N₂O, and short-lived climate forcers (SLCFs) are taken from the Global Fire Emissions Database (GFED) (van der Werf et al., 2017) v4.1, which includes the BB4CMIP dataset prepared for CMIP6 historical simulations (van Marle et al., 2017).

• Emissions of Montreal Protocol greenhouse gases (CFCs, HCFCs, halons, and chlorinated and brominated gases), along with SO₂F₂, are estimated using inverse greenhouse gas concentrations that have been prepared for the IGCC (Forster et al., 2023), as no inventories of these emission datasets are available to our knowledge.

All emission datasets are produced for 1750–2022, except CEDS, which has a 2019 end-date. To extend SLCFs from CEDS to 2022, we use the “2-year blip” scenario that estimates the decline and recovery from emissions due to COVID-19 from Forster et al. (2020) and is extended by Lamboll et al. (2021), based on proxy activity data. We take the ratios of SLCF emission species over 2020–2022 to 2019 in the 2-year blip scenario and apply them as a scaling factor to CEDS emissions in 2019. Such a version-controlled strategy allows for the calibration to be updated as newer emission data become available. Emission data prepared to the end of 2023 will be available over the course of 2024, and an anticipated update to CEDS should also bring non-biomass-burning SLCFs until at least the end of 2022 (Hoesly et al., 2023). This demonstrates that “operational” calibrations are often a moving target.

We use the “third-party” emissions from PRIMAP-Hist rather than country-reported values, based on the assumption that we expect solely country-reported values to be an underestimate of true emissions. We demonstrate that third-party emissions still appear to be an underestimate for many species, based on best-estimate greenhouse gas lifetimes and concentration estimates.

2.3.2 v1.4.0: RCMIP historical emissions prepared for AR6 (1750–2014)

For consistency and comparison with the FaIR projections used in the IPCC AR6, we produce a calibration using historical emissions from RCMIP (Nicholls et al., 2020, 2021) using v5.1.0 of the Reduced Complexity Model Intercomparison Project (RCMIP) emission dataset available from Nicholls and Lewis (2021). The RCMIP emissions contain global annual total emissions of CO₂ and SLCFs that were prepared for running CMIP6 models. Emissions of non-CO₂ greenhouse gases were back-calculated to reproduce the CMIP6 best-estimate historical concentrations (Meinshausen et al., 2017). These concentrations time series were also used to drive CMIP6 models.

For SSP scenarios, emissions from 2015 to 2100 were produced using IAMs, which were then extended to 2500 using simplified assumptions (Meinshausen et al., 2020). We use the same climate constraints on GMST, CO₂ concentration, and OHC as for v1.4.1 (Sect. 3.3) datasets, which run to 2022. For the bridging period 2015–2022 between the end of the CMIP6 historical and the observational climate data, we use emissions from SSP2-4.5, expected to be the closest Tier 1 SSP to current policies (Hausfather and Peters, 2020) and, as shown later, the closest Tier 1 scenario to post-2015 emissions.

One adjustment is made to the RCMIP emissions to correct NO_x. For accounting purposes, we express NO_x in units of Tg NO₂ yr⁻¹. The source datasets for RCMIP were earlier versions of CEDS, which reports emissions in Tg NO₂ yr⁻¹ for fossil fuel and agricultural emissions, and GFED, which reports emissions in Tg NO yr⁻¹ for biomass burning. The conversion for GFED emission data was not made in RCMIP v5.1.0.

Neither v1.4.1 nor v1.4.0 of fair-calibrate includes forcing from aviation contrails. Forcing from contrails and its temperature impact were assessed in the IPCC AR6 WG1 (Forster et al., 2021), with best-estimate contributions to present-day forcing of 0.06 W m⁻² and warming of 0.02 °C, and were included in the WG1 calibration of FaIR. However, contrail forcing was excluded from the WG3 IAM emission projections, rendering the WG1 and WG3 projection sets slightly inconsistent. To project contrail forcing into the future requires estimates of aviation activity. FaIR can accept a time series of contrails forcing directly or estimate it, using a linear combination of emission species. By default, FaIR uses NO_x emissions from the aviation sector to estimate contrail forcing (Smith et al., 2018). Neither aviation activity nor NO_x emissions from aviation are provided in IAM scenarios in general, so contrail forcing could not be assessed in WG3. Aviation NO_x emissions are provided in the RCMIP historical and SSP future emissions and could be included in fair-calibrate v1.4.0. However, in order to apply the calibrations consistently to as many scenarios as possible, we calibrate without them.

3 Process

The set of output FaIR parameters is produced in three steps: (1) calibration, (2) sampling, and (3) constraining. The description and results in this section apply generally to all calibration versions to date. We focus on calibration v1.4.1 and describe methods pertinent to v1.4.0 where they differ. Figure 1 details the general process chain of fair-calibrate v1.4.1.

Figure 1 Schematic of the process chain in fair-calibrate v1.4.1. Square brackets detail sources of data, and round brackets detail section numbers in which processes are described in more detail. Dashed borders are optional processes which are not required to calibrate the history only.

3.1 Calibration

3.1.1 Climate response

The climate response module of FaIR v2.1.3 is an impulse response formulation of the three-layer stochastic energy balance model of Cummins et al. (2020). We calibrate this model using 150-year 4×CO₂ experiments from 49 CMIP6 models and using GMST (ΔT₁) and the top-of-atmosphere energy imbalance (ΔN) as anomalies relative to each model's pre-industrial control run, subtracting a linear trend from the appropriate branch point of each model's control to account for any residual drift. This calibration is performed using the maximum likelihood method of Cummins et al. (2020), and the EBM R package that accompanies Cummins et al. (2020) is used in the fair-calibrate process chain (Cummins, 2021).

The three-layer stochastic energy balance model is written as

$$\begin{array}{}\text{(1)}& {\displaystyle }{C}_{\mathrm{1}}{\displaystyle \frac{\mathrm{d}{T}_{\mathrm{1}}\left(t\right)}{\mathrm{d}t}}& {\displaystyle }=F\left(t\right)-{\mathit{\kappa }}_{\mathrm{1}}{T}_{\mathrm{1}}\left(t\right)-{\mathit{\kappa }}_{\mathrm{2}}\left(T_{\mathrm{1}}\right(t)-T_{\mathrm{2}}(t\left)\right)+\mathit{\xi }\left(t\right),\text{(2)}& {\displaystyle }{C}_{\mathrm{2}}{\displaystyle \frac{\mathrm{d}{T}_{\mathrm{2}}\left(t\right)}{\mathrm{d}t}}& {\displaystyle }={\mathit{\kappa }}_{\mathrm{2}}\left(T_{\mathrm{1}}\right(t)-T_{\mathrm{2}}(t\left)\right)-\mathit{\epsilon }{\mathit{\kappa }}_{\mathrm{3}}\left(T_{\mathrm{2}}\right(t)-T_{\mathrm{3}}(t\left)\right),\text{(3)}& {\displaystyle }{C}_{\mathrm{3}}{\displaystyle \frac{\mathrm{d}{T}_{\mathrm{3}}\left(t\right)}{\mathrm{d}t}}& {\displaystyle }={\mathit{\kappa }}_{\mathrm{3}}\left(T_{\mathrm{2}}\right(t)-T_{\mathrm{3}}(t\left)\right).\end{array}$$

In Eqs. (1)–(3), T₁, T₂, and T₃ are the temperature anomalies of the three ocean layers (starting from the surface); C₁, C₂, and C₃ are their heat capacities; κ_j represents the heat transfer coefficients between layers j−1 and j for j≥2; −κ₁ is the climate feedback parameter (often denoted λ); ε is the deep-ocean efficacy parameter (Held et al., 2010; Winton et al., 2010; Geoffroy et al., 2013); ξ is a stochastic disturbance term in the temperature response that does not affect the top-of-atmosphere energy imbalance; and F is the effective radiative forcing (ERF).

The effective radiative forcing is the sum of a deterministic and stochastic component $F=F_{\text{det}}+\mathit{\zeta }$. The stochastic forcing component ζ is modelled as a continuous-time red noise process

$$\begin{array}{}\text{(4)}& {\displaystyle \frac{\mathrm{d}\mathit{\zeta }}{\mathrm{d}t}}=-\mathit{\gamma }\mathit{\zeta }+\mathit{\eta },\end{array}$$

where η is white noise, and γ>0 controls the strength of temporal auto-correlation (Cummins et al., 2020). In FaIR, the stochastic behaviour can be switched off, and Eqs. (1)–(4) reduce to a deterministic energy balance model when $\mathit{\xi }=\mathit{\eta }=\mathrm{0}$ (Geoffroy et al., 2013; Leach et al., 2021).

The top-of-atmosphere energy imbalance N is given as

$$\begin{array}{}\text{(5)}& N\left(t\right)=F\left(t\right)-{\mathit{\kappa }}_{\mathrm{1}}{T}_{\mathrm{1}}\left(t\right)+(\mathrm{1}-\mathit{\epsilon }){\mathit{\kappa }}_{\mathrm{3}}\left(T_{\mathrm{2}}\right(t)-T_{\mathrm{3}}(t\left)\right),\end{array}$$

and the Earth's energy uptake, used as a model constraint, is the time integral of N.

For each of the 49 CMIP6 models, we obtain a set of 11 parameters $\mathit{\{}{C}_{\mathrm{1}},C_{\mathrm{2}},C_{\mathrm{3}},{\mathit{\kappa }}_{\mathrm{1}},{\mathit{\kappa }}_{\mathrm{2}},{\mathit{\kappa }}_{\mathrm{3}},\mathit{\epsilon },\mathit{\gamma },{\mathit{\sigma }}_{\mathit{\xi }},{\mathit{\sigma }}_{\mathit{\eta }},F_{\mathrm{4}\times {\mathrm{CO}}_{\mathrm{2}}}\mathit{\}}$ that describes the magnitude and rate of warming to a 4×CO₂ forcing and the behaviour of internal variability, where σ_ξ and σ_η are the standard deviations of ξ and η around the zero mean. $F_{\mathrm{4}\times {\mathrm{CO}}_{\mathrm{2}}}$ is the effective radiative forcing from a quadrupling of pre-industrial CO₂ concentrations. The comparison of one stochastic realisation of each model's energy balance model calibration (black) compared to the actual CMIP6 model (red) for the temperature response to an abrupt 4×CO₂ forcing is shown in Fig. 2. In almost all cases, the FaIR calibration is an excellent representation of the underlying CMIP6 model. The calibrated parameters are shown in Table S1.

Figure 2 Comparison of temperature projections from abrupt 4×CO₂ simulations as calibrated in FaIR (black) to the original CMIP6 model results (red) for 49 CMIP6 models. For FaIR, we show one realisation with stochastic internal variability included; different random seeds would produce different internal variability profiles.

The energy balance model parameters can be written as a matrix equation that describes the time evolution of each temperature layer (Cummins et al., 2020; Leach et al., 2021). The impulse response form of the temperature evolution in each layer can be calculated from the eigenvalues and eigenvectors of the energy balance matrix. From this, the ECS and “theoretical” TCR for each model calibration can be directly estimated from the impulse response coefficients as described in Leach et al. (2021, Sect. 2.4). The ECS calculated here is a true equilibrium value rather than as a regression over a 150-year simulation as usually performed from ESM output (the so-called effective sensitivity, EffCS). The theoretical TCR is not precisely what each model would predict after 70 years of a 1 % compound increase in atmospheric CO₂ concentrations but is usually close and has the advantage that model simulations do not need to be run to determine this value (Fig. S1b in the Supplement).

3.1.2 Minor greenhouse gas emissions

This section describes the emission adjustment procedure in fair-calibrate v1.4.1 for emissions of minor greenhouse gases. In this context, “minor” means any species that is not CO₂ or CH₄. This includes N₂O, hydrofluorocarbons (HFCs), perfluorocarbons (PFCs), SF₆, and NF₃. This emission adjustment is not required in v1.4.0, where emissions from all species are provided by the RCMIP emission datasets (Sect. 2.3.2).

HFCs and PFCs are provided in PRIMAP-Hist as aggregate values reported in CO₂-equivalent (AR6 GWP₁₀₀) emissions. We disaggregate these emissions by scaling the annual historical emission totals in CO₂-equivalent emissions from RCMIP historical + SSP2-4.5 for 1750–2022 to the PRIMAP-Hist reported values and then by multiplying this scaling by the RCMIP individual species emission value in each year. Table S2 details the HFC and PFC gases included in the disaggregation.

The following step calculates atmospheric concentrations when run forward using a single time-constant decay model with the PRIMAP-Hist emission and time constants equal to atmospheric lifetimes assessed in IPCC AR6 (Smith et al., 2021b). The calculated concentration time series is compared to the best-estimate historical concentrations from Forster et al. (2023), which is an update of the AR6 concentrations in IPCC (2021) to 2022 using recent AGAGE and NOAA station data. In many cases, the calculated and observed concentrations differ substantially, and the calculated concentrations are usually lower than the observed. This implies that either the reported emissions in PRIMAP-Hist do not capture all true emissions or that the reported atmospheric lifetimes are too short (a third, less likely, possibility is that the reported concentrations are too high). A correction can be obtained by either lengthening the lifetimes or scaling up the emissions. We choose to adjust the emissions on the basis that countries under-reporting due to incomplete data is plausible, and scaling the emissions brings some species much closer to RCMIP estimates which are derived from inverting atmospheric concentrations. The scaling is performed in order to match the projected concentrations to the historical best estimates in 2019. In many cases the scaling is mild (for N₂O, emissions are scaled up by a factor of 1.08; Fig. 3a) but can be large (NF₃ is scaled by a factor of 7.5; Fig. S2). This implies that countries are severely under-reporting emissions of some greenhouse gases (GHGs) compared to the increasing stock of these gases observed in the atmosphere.

PRIMAP-Hist does not provide emissions of SO₂F₂ or of Montreal Protocol GHGs. We estimate their emissions by inverting the concentrations time series in Forster et al. (2023).

For future projections, we harmonise to 2022 (Gidden et al., 2018) the eight Tier 1 and Tier 2 SSP scenarios to our scaled calculated historical emissions. This produces SSPs that take into account the recent past. We can then compare the harmonised adjusted future concentration projections to those created for the SSP scenarios that used MAGICC6 (Meinshausen et al., 2020). Figure 3b shows recreated historical and future N₂O concentration projections to 2100 under eight SSP scenarios using the harmonised scaled emissions (thick lines) in FaIR and their comparison to the SSP concentrations time series (thin lines) from Meinshausen et al. (2017, 2020). Note that the historical concentrations differ between Fig. 3a and b as the dataset sources differ. For N₂O, the correspondence between FaIR and CMIP6 is very good for all eight SSPs for future projections.

Figure 3 (a) Comparison of best-estimate historical N₂O emissions (black), the concentration projected from emissions in PRIMAP-Hist + GFED (dotted grey), and the concentrations after scaling up the emissions by a factor of 1.08 to get correct recent historical concentrations (solid grey). Note that a single lifetime cannot accurately reproduce best-estimate historical concentrations between 1850 and 1950. (b) Harmonised SSP projections using the scaled historical emissions (thick lines) compared to the SSP historical + future projections (thin lines) from Meinshausen et al. (2017, 2020).

3.1.3 Methane lifetime

A new feature of FaIR introduced in v2.1.0 is a variable methane lifetime that depends on burdens of chemically reactive species and climate. This is an update from v2.0.0 that used a methane lifetime self-feedback (methane concentrations and temperature affect climate) and previous versions that did not modify the lifetime of methane at all.

A methane lifetime scaling factor ${\mathit{\alpha }}_{{\mathrm{CH}}_{\mathrm{4}}}$ is applied to the base lifetime ${\mathit{\tau }}_{{\mathrm{CH}}_{\mathrm{4}},\mathrm{base}}$ calculated as

$$\begin{array}{}\text{(6)}& \log {\mathit{\alpha }}_{{\mathrm{CH}}_{\mathrm{4}}}=\log \left(\mathrm{1}+S_{\mathrm{T}}\mathrm{\Delta }{T}_{\mathrm{1}}\right)+\sum _i\log \left(\mathrm{1}+S_i\mathrm{\Delta }{A}_i\right).\end{array}$$

In Eq. (6), S_i denotes a sensitivity to species i or GMST anomaly (ΔT₁), and ΔA_i represents abundances of species i (emissions rate for SLCFs and concentrations for GHGs) of chemically reactive species. If the anomalies in temperature and abundances are relative to the pre-industrial period, ${\mathit{\alpha }}_{{\mathrm{CH}}_{\mathrm{4}}}=\mathrm{1}$ in pre-industrial conditions and ${\mathit{\tau }}_{{\mathrm{CH}}_{\mathrm{4}},\mathrm{base}}$ is the pre-industrial lifetime.

Unlike for minor GHGs, emissions are not scaled for CH₄ in fair-calibrate v1.4.1, and we instead calibrate the atmospheric chemical lifetime. Owing to dependence of the lifetime of several simultaneously changing emission species, as well as climate, there is not a unique invertible concentration to emission pathway for methane.

The UKESM1.0-LL, GFDL-ESM4, GISS-E2.1-G, and MRI-ESM2.0 Earth system models provide a complete set of results from the Aerosol Chemistry Model Intercomparison Project (AerChemMIP) single-forcing experiments that enable the estimation of the sensitivity in methane lifetime to climate (Thornhill et al., 2021a) and chemically reactive species (Thornhill et al., 2021b). We use results reported in Thornhill et al. (2021b) and Thornhill et al. (2021a) for methane lifetime in 1850 and its relative sensitivity to changes in CH₄, N₂O, and equivalent effective stratospheric chlorine (EESC) concentration; emissions of NO_x and VOCs; and global mean surface temperature between 1850 and 2014 in each of the four models. For each atmospheric species, the fractional change in lifetime in 2014 relative to 1850 is normalised by the burden change to provide lifetime changes in each model in terms of parts per billion (ppb) concentration change or Mt yr⁻¹ emissions. The four models that provide data are used as minimum and maximum ranges of a parameter search (in v1.4.1, we expand the search range by a factor of 2, since the PRIMAP-Hist methane emissions are again likely to be an underestimate and do not find suitable parameters within the model range) to minimise the difference between observed CH₄ concentrations from Forster et al. (2023) and those calculated from Eq. (6). The 1750 emissions are subtracted from the time series when performing the lifetime calibration, as it is assumed that pre-industrial concentrations of methane are in approximate equilibrium with pre-industrial emissions.

The historical best-estimate calibrations are shown in Table 1. It can be seen that the methane lifetime in fair-calibrate v1.4.1 is nearly 17 years in the pre-industrial period, which is much longer than typically determined from ESMs. The best-estimate lifetime in FaIR from historical emissions is shown in Fig. 4a (grey line) and is indeed longer than that calculated from the sensitivities in each CMIP6 model across most of the historical period, though close to the AR6 value in the present day. In Fig. 4b, the historical concentrations from Forster et al. (2023) (black) are compared to the best estimate from FaIR using the lifetime calculated in Fig. 4a and run forward with best-estimate historical emissions. In Fig. 4c, the SSP methane concentrations are projected with the harmonised emissions, starting in 2022, and compared to the SSP concentrations time series (Meinshausen et al., 2017, 2020). In general, the harmonised methane concentration projections from fair-calibrate v1.4.1 are lower than in CMIP6 for high-methane emission scenarios and higher for low-emission futures. This is due in part to the nearly 10 years of additional historical emissions in the best-estimate time series compared to the SSPs, which started to diverge from a common history in 2015. For these projections, we use the best-estimate GMST anomalies from the SSPs derived in Lee et al. (2021).

Figure 4 Methane lifetime calibration (v1.4.1). (a) Methane lifetime in the historical + SSP3-7.0 scenario for four ESMs (colours) and the lifetime from the FaIR calibration (grey). (b) Methane concentration calculated from historical methane emissions from PRIMAP-Hist + biomass burning emissions using the lifetime in panel (a), using FaIR (grey), and using the observed atmospheric concentrations (black) for 1750–2022 from IGCC (Forster et al., 2023). (c) Methane concentrations calculated from methane emissions for the eight main SSP scenarios using the harmonised future emission projections (thick lines) compared to the SSP scenarios (thin lines) (Meinshausen et al., 2017, 2020).

Table 1 Baseline CH₄ lifetime and sensitives (S_i) in lifetime due to changes in greenhouse gas concentrations, short-lived climate forcer emissions, and temperature in calibrations v1.4.1 and v1.4.0. Note that ppt stands for parts per trillion.

The lifetimes, historical concentrations, and future concentrations for the RCMIP emissions (calibration v1.4.0) are shown in Fig. S3, where it is observed that lifetimes and concentration projections are much closer to AR6 and CMIP6. This demonstrates that, first, the calibration is plausible (CMIP6 emissions give CMIP6 concentrations) and, second, that the methane lifetime calibration is very sensitive to the historical emission time series used. In Fig. S4b, we compare the methane emissions from the v1.4.0 and v1.4.1 calibrations. As 1750 emissions are subtracted from the total to report changes away from a pre-industrial equilibrium, the change in emissions (1750–2022) in v1.4.1 from PRIMAP-Hist is smaller than in v1.4.0, leading to longer atmospheric lifetimes necessary to reproduce concentrations.

Unlike in versions of FaIR prior to 2.0.0, we do not assume any natural methane emissions. In v1.3 of FaIR, for example, natural emissions were back-calculated with the assumption of a constant methane lifetime and held constant for future projections (Smith et al., 2018). It is well-known that wetlands emit large quantities of methane, and it is very likely that this effect is climate-dependent (Zhang et al., 2017). As the climate continues to warm, biogenic methane will be released from permafrost soils and clathrates – sources that most ESMs do not include at present. Including these natural sources is a development priority for future versions of FaIR.

It should be noted that the methane lifetimes derived are the best fits to observed concentrations across the 1750 to present-day period for each emission pathway and may not necessarily maintain an equilibrium concentration in 1750 with 1750 emissions. In v1.4.1, methane emissions in 1750 were around 38 Mt CH₄ and around 19 Mt CH₄ in v1.4.0, though v1.4.0 has a shorter lifetime for the same concentration. Methane emissions were not in equilibrium in 1750 and have steadily climbed over the last 2000 years (Meinshausen et al., 2017), with substantial variations due to agricultural and natural influences before then (Singarayer et al., 2011). Methane's relatively short lifetime and reactive nature make its calibration more difficult than longer-lived greenhouse gases such as CO₂ and N₂O, and the calibration strategy of the methane cycle depends on the goal of the user. In most cases using FaIR, this will be historical and future anthropogenic influences on climate for which the calibration that ensures historical emissions reproduce historical concentrations is most appropriate. Other use cases may require different calibration strategies.

3.1.4 Carbon cycle feedbacks

The carbon cycle is parameterised as a simple atmospheric decay model with four time constants, based on the impulse response functions of Joos et al. (2013). The time constants are scaled by a lifetime scaling factor that mimics the influence of carbon cycle feedbacks. This treatment is unchanged since the work of Leach et al. (2021, Sect. 2.1). A positive carbon cycle feedback reduces the efficacy of carbon sinks, thus effectively lengthening the atmospheric lifetime of CO₂.

The lifetime scaling factor is a function of the time-integrated airborne fraction of a CO₂ pulse over 100 years I₁₀₀ (Millar et al., 2017). I₁₀₀ is modified as

$$\begin{array}{}\text{(7)}& I_{\mathrm{100}}=r_{\mathrm{0}}+r_{\mathrm{U}}\mathrm{\Delta }{C}_{\mathrm{U}}+r_{\mathrm{T}}\mathrm{\Delta }T+r_{\mathrm{A}}\mathrm{\Delta }{C}_{\mathrm{A}},\end{array}$$

where r₀, r_U, r_T, and r_A are the pre-industrial time-integrated airborne fraction and its sensitivity to cumulative carbon uptake in land and ocean sinks ΔC_U, surface temperature anomaly ΔT, and airborne carbon ΔC_A respectively. Total cumulative emissions since pre-industrial is ΔC_A+ΔC_U.

The process for calibrating the carbon cycle feedbacks to 11 CMIP6 ESMs containing interactive carbon cycles is described in Leach et al. (2021, Sect. 3.2). The same coefficients derived in Leach et al. (2021) for the 11 ESMs are used in all calibrations to date.

3.1.5 Aerosol–cloud interactions

The effective radiative forcing due to aerosol–cloud interactions ERF_aci has been generalised:

$$\begin{array}{}\text{(8)}& {\mathrm{ERF}}_{\mathrm{aci}}=\mathit{\beta }\left[\log \left(\mathrm{1}+\sum _i{s}_i{A}_i\right)-\log \left(\mathrm{1}+\sum _i{s}_i{A}_{i,\mathrm{base}}\right)\right],\end{array}$$

where A_i is the emissions or concentration of a species, and the base subscript denotes its reference (usually pre-industrial) abundance. β is a scale factor, and s_i describes how sensitive a species is in contributing to ERF_aci. The generalisation allows for inclusion of more species that affect ERF_aci in addition to SO₂, BC, and OC that was modelled previously. The generalisation is useful as there is evidence of a large ERF_aci response to CH₄ in UKESM1-0-LL through methane's effect on competing for atmospheric oxidants, including OH, affecting the rate of new particle formation (O'Connor et al., 2022). As with earlier versions of FaIR, the form of Eq. (8) is inspired by Stevens (2015) but without any physical significance attached to the sensitivities s_i, allowing near-linear global mean responses in ERF_aci to changes in precursor abundances as postulated by some authors (Booth et al., 2018; Kretzschmar et al., 2017) and exhibited in some models (Smith et al., 2021a).

A total of 13 CMIP6 models provided results from transient aerosol experiments in AerChemMIP and the Radiative Forcing Model Intercomparison Project (RFMIP) (Table 2) that allow calculation of aerosol ERF. The breakdown of shortwave aerosol ERF into aerosol–radiation interactions (ERF_ari) and ERF_aci is performed using the approximate partial radiative perturbation (APRP) method (Taylor et al., 2007), following the logic of Zelinka et al. (2014, 2023). Longwave contributions to ERF_aci are estimated from the cloud radiative effect, with ERF_ari estimated as the difference between the longwave components of ERF and ERF_aci.

From the diagnosed ERF_aci in each model, a least squares curve fit of ERF_aci to historical emissions by fitting $s_{{\mathrm{SO}}_{\mathrm{2}}}$, s_BC, s_OC, and β is found (Table 2) using Eq. (8). The comparison of model-derived ERF_aci to the best fit from Eq. (8) is shown in Fig. 5.

Using Eq. (8), a wide range of ERF_aci trajectories are possible, and parameter estimates for β and individual species sensitivities span orders of magnitude. Where one or two of $s_{{\mathrm{SO}}_{\mathrm{2}}}$, s_BC, and s_OC are close to zero (CanESM5 and UKESM1-0-LL), this indicates that the species has little influence on ERF_aci in that model (e.g. UKESM1-0-LL's ERF_aci response is purely driven by sulfate in aerosol-only forcing experiments). Where all three of $s_{{\mathrm{SO}}_{\mathrm{2}}}$, s_BC, and s_OC are close to zero, and β has large magnitude (the two Geophysical Fluid Dynamics Laboratory (GFDL) models and NorESM2-LM), this indicates that ERF_aci behaves linearly in emissions from the Taylor expansion of log (1+x) for small x (Smith et al., 2021a). In the case of NorESM2-LM, the coefficient for BC is so small that it is effectively zero, with the ERF_aci response being linear with sulfate and OC.

Table 2 Models used to calibrate forcing from aerosol–cloud radiation interactions and their parameter best fit values from Eq. (8).

Figure 5 Calibrations of the ERF_aci relationship in FaIR (Eq. 8; coloured lines) to the derived ERF_aci from 13 CMIP6 models (grey lines). Extrapolation back to 1750 is shown in all cases, and extrapolation forward to 2100 is shown under SSP2-4.5 emissions where model simulations were not extended beyond 2014.

3.1.6 Ozone

The best-estimate historical ozone ERF time series from Skeie et al. (2020) is used to calibrate the role of ozone precursors to ozone forcing. As in AR6, tropospheric and stratospheric ozone are not considered separately. Again following the AR6 methodology, we select six models from the 12 coupled historical CMIP6 models analysed in Skeie et al. (2020) that are relatively independent from each other, have full stratospheric and tropospheric chemistry enabled, and reproduce expected behaviour for the overall time history of ozone ERF. The six models used are BCC-ESM1, CESM2(WACCM6), GFDL-ESM4, GISS-E2-1-H, MRI-ESM2-0, and OsloCTM3. Skeie et al. (2020) provides historical ozone forcing for 1850–2010 in these models, and following Skeie et al. (2020), we add +0.03 W m⁻² to the time series to represent the change from 1750 to 1850. The Oslo-CTM3 model provided results under SSP2-4.5 to 2020, which was also used in calibration.

As ozone ERF includes a contribution from temperature change and is calibrated from coupled historical runs, historical warming is backed out using a temperature feedback of −0.037 W m⁻² K⁻¹ (Thornhill et al., 2021a) and historical GMST from Forster et al. (2023). For this “no-feedback” ERF time series, we find a least squares fit to the change in emissions of NO_x, VOC, and CO and concentrations of CH₄, N₂O, and EESC (Fig. 6). The lower and upper bounds of the search ranges for the parameter fits are the very likely range for each precursor in Thornhill et al. (2021b), which is also scaled up to account for the difference in best-estimate ozone forcing between models participating in AerChemMIP in Thornhill et al. (2021b) and the six-model subset in Skeie et al. (2020).

Figure 6 Comparison of the ozone ERF time series from Skeie et al. (2020) (black) to the estimate from emissions and concentration precursors (grey). The estimated impact of temperature on ozone forcing has been backed out of the time series from Skeie et al. (2020) and is not included in the model fit.

Similar to the methane lifetime calibration, we derive a coefficient for each precursor species relating emissions or concentrations of each to the ozone ERF. Uncertainty sampling for the prior distribution is described in Sect. 3.2.5.

3.2 Sampling

We produce a 1.6 million member prior ensemble of FaIR projections, with parameter choices drawn from probability distributions that are informed by CMIP6 model calibrations (Sect. 3.1) or AR6-assessed ranges. Different components of FaIR are sampled independently, but within each component (e.g. climate response), the correlation structure between parameters is maintained to ensure internally consistent parameter choices. In many cases, probability distributions for parameters are constructed from a Gaussian kernel density estimate, which is a non-parametric method that attempts to estimate the underlying probability density function from a finite sample size, and can be used to preserve correlation structure in multi-variate cases (Scott, 1992).

Kernel density estimates to sample parameters are used since several parameters do not have many CMIP6 models to calibrate to (a data-sparsity issue), parameter values can span several orders of magnitude, and correlations between parameters that arise from the calibration can be included. In each case, we use the scipy.stats.gaussian_kde implementation of the multivariate kernel density estimate (https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.gaussian_kde.html, last access: 30 June 2024). Including the correlation between parameters reduces (though does not eliminate) the likelihood of physically implausible combinations being sampled, and using kernel density estimates rather than parametric multivariate distributions allows for variability in the distribution shapes of each parameter, such as admitting left-skewed and multi-modal shapes. Kernel density methods have drawbacks, such as being sensitive to outliers. However, parametric distributions assume some prior knowledge about the dataset, and selecting one model per parameter does not fully sample the potential space of plausible climate models.

In this section, prior distributions that are not sampled from kernel density estimate calibrations to CMIP6 models are shown in individual tables.

In total, 45 parameters are sampled. In the processing chain, fixed random seeds are used to ensure reproducibility. Internal variability is switched on, and again each parameter set has a random seed associated with it in order to reproduce the same pattern, and key climate metrics are saved out of the prior ensemble.

3.2.1 Climate response

An 11-dimensional kernel density estimate is generated from the energy balance model parameters that were calibrated on 49 CMIP6 models (Fig. S5). $F_{\mathrm{4}\times {\text{CO}}_{\mathrm{2}}}$ is not used in the climate response of FaIR but is used in the theoretical calculation of ECS and TCR. All parameters of the energy balance model are strictly positive, so parameter sets containing negative values are discarded and redrawn until the 1.6 million threshold is reached. We also discard and redraw instances of κ₁<0.3 W m⁻² K⁻¹, C₁<1.8 W yr m⁻² K⁻¹, C₃<C₂, C₂<C₁, and γ<0.5. The κ₁ threshold puts an upper bound on the ECS prior of around 13 °C, and the other limits ensure model stability.

3.2.2 Aerosol–cloud interactions

Similar to the climate response, we draw correlated kernel density estimates for $\log (s_{{\mathrm{SO}}_{\mathrm{2}}}$), log (s_BC), and log (s_OC). We calculate an unscaled ERF_aci for the 2005–2014 mean relative to 1750 for each parameter set. The unscaled ERF_aci is then scaled to reproduce a draw from a trapezoid distribution with limits at −2.2 and +0.2 W m⁻² and plateau from −1.6 to −0.4 W m⁻² to represent the ERF_aci for 2005–2014 relative to 1750, which selects the β value to use for that parameter set. This process is similar to that of both Smith et al. (2021a) and AR6 (Forster et al., 2021). The prior distribution is chosen to give a wide but plausible range around the ERF_aci distribution for the present day assessed by the IPCC (Forster et al., 2021), which was −1.0 W m⁻² for a nominal 2014 date relative to 1750.

3.2.3 Aerosol–radiation interactions

The ERF_ari contributions are not sampled directly from CMIP6 models, though much of the basis of this assessment is rooted in AerChemMIP (Thornhill et al., 2021b). AR6 assessed that several species (CH₄, N₂O, halogenated compounds, sulfate, BC, OC, nitrate, and VOCs) contribute directly or indirectly to ERF_ari, though only sulfate, BC, OC, and NH₃ are significant. We use the contributions to ERF_ari assessed in AR6 with the relative uncertainty from each precursor (Szopa et al., 2021) as prior distributions (Table 3) and scale both the best-estimate and uncertainty range of the ERF_ari from each precursor to reproduce the IPCC AR6 distribution of $-\mathrm{0.3}\pm \mathrm{0.3}$ W m⁻² (Forster et al., 2021). All ranges quoted are for 5th to 95th percentile, unless otherwise stated.

Table 3 Distributions of the contributions to the direct aerosol ERF sampled in fair-calibrate v1.4.1. Uncertainty ranges are shown as 90 % ranges and sampled from a Gaussian distribution.

3.2.4 Carbon cycle and initial CO₂ concentration

A four-dimensional kernel density estimate is drawn from the r₀, r_U, r_T, and r_A parameters from the 11 models calibrated in Leach et al. (2021). As part of the carbon cycle sampling, we draw CO₂ concentration values in 1750 using the IPCC AR6 best estimate and uncertainty of the 278.3±2.9 ppm (5 %–95 %) range (Gulev et al., 2021), using a Gaussian distribution.

3.2.5 Ozone

The coefficients relating emissions or concentrations of chemically relevant precursors to ozone ERF take their mean value from the bounded least square fit derived in Sect. 3.1.6, and their uncertainty values are sampled by applying the scaled 5 %–95 % uncertainty range from Thornhill et al. (2021b) to this best-estimate value. This means that some precursor ranges are outside the range of that described by Thornhill et al. (2021b), though only seven models (fewer for some precursors) provided the necessary experiments in Thornhill et al. (2021b), and thus AerChemMIP represents a small ensemble of opportunity.

Table 4 Distributions of the contributions to the ozone ERF sampled in fair-calibrate v1.4.1. Uncertainty ranges are shown as 90 % ranges and sampled from a Gaussian distribution.

3.2.6 ERF scalings

Forcing uncertainties in ERF_ari, ERF_aci, and ozone are sampled from the contribution to total forcing from their precursor species, as described in previous sections. For other major categories of forcings, we use the IPCC AR6 ranges (Forster et al., 2021) as relative uncertainty factors to scale the ERF (Table 5).

For CO₂, we use the sampled $F_{\mathrm{4}\times {\text{CO}}_{\mathrm{2}}}$ value from the climate response calibration and perform a quantile mapping to derive a scaling factor for CO₂ forcing that is Gaussian. While this does not preserve the shape of the $F_{\mathrm{4}\times {\text{CO}}_{\mathrm{2}}}$ distribution kernel, it does map low 4×CO₂ forcings to low CO₂ scalings, and vice versa.

Table 5 Forcing scaling factors used to translate the raw best estimate from FaIR to IPCC-assessed uncertainty ranges (Forster et al., 2021). Scaling uncertainty ranges are 5 %–95 %. Except for the solar trend, median distribution values are 1.

^* Contrail forcing is not used in v1.4.0 and v1.4.1 but is included in other versions.

3.3 Constraining

The 1.6 million member prior ensemble of FaIR climate projections is compared to historical observations and assessments of climate metrics from either the IPCC AR6 (Forster et al., 2021) or their updates, based on more recent data (Forster et al., 2023).

3.3.1 Step 1: root mean squared difference with respect to historical

The root mean squared (rms) difference in each ensemble member's GMST anomaly projection compared to the historical values for 1850–2022 is used as a simple pass/fail criterion for ruling out parameter sets that are inconsistent with historical observed warming. Ensemble members that have an rms difference that is greater than 0.17 °C are rejected. The mean of four GMST datasets (HadCRUT5, Berkeley Earth, NOAAGlobalTemp, and Kadow) from Forster et al. (2023) is used as the historical GMST dataset for comparison. The choice of 0.17 °C is somewhat arbitrary, which balances sufficient variability in the historical record to allow for observational uncertainty with the need for projections that are true to observations. By design, this threshold roughly reproduces the uncertainty range in present-day GMST relative to the pre-industrial range assessed by the IPCC (Gulev et al., 2021), whereas a more stringent threshold may over-constrain both the historical observational uncertainty and scope for future climate projection uncertainty (Fig. 7). Internal variability is switched on for this historical comparison to allow for the possibility that the historical record can be well-simulated by chance in mean state climate configurations that would be warmer or cooler than expected (e.g. a strong pattern effect; Andrews et al., 2018). This step reduces the ensemble size from 1.6 million to 224 342, ruling out around 86 % of the original ensemble.

Figure 7 compares the 10 ensemble members with the lowest RMSE relative to observations (blue; RMSE ≈0.10 °C) with the 10 largest RMSE members that still meet the RMSE constraint (red; RMSE ≈0.17 °C). Figure 7 shows that runs with low internal variability tend to result in the closest correspondence with historical observed temperature, and therefore, the final ensemble could be biased towards ensemble members with smaller variability. A formal analysis of the internal variability characteristics in relation to observations is not performed in this version of fair-calibrate, though it could be added to the constraining criteria in the future.

Alongside or instead of RMSE, a correlation metric could be used to evaluate goodness of fit between the observations and the model. However, RMSE encapsulates goodness of fit into a single number and is sensitive to model runs that overall warm too quickly or too slowly. Correlation coefficients would not differentiate simulations that had the right shape of historical warming but warmed too quickly or slowly.

Figure 7 Comparison of the 10 ensemble members with the smallest RMSE error (blue) compared to the historical best-estimate GMST from the Indicators of Global Climate Change 2022 (Forster et al., 2023) (black) with 10 of the largest RMSEs (red) that passed this first historical constraining step.

3.3.2 Step 2: reweighting based on observed and assessed climate metrics

The second constraining step takes the ensemble members that passed the RMSE threshold and simultaneously fits the projections to eight target distributions (Fig. 8). For each target distribution, either a Gaussian (if symmetric) or skew normal (if asymmetric) continuous probability distribution is constructed from the 5th, 50th, and 95th percentiles of the variable's uncertainty range. As a three-parameter distribution, a skew normal can uniquely fit three specified quantiles. For symmetric distributions, the number of degrees of freedom is reduced to two (by imposition of symmetry), and the Gaussian is a natural choice, as well as being a general form of the skew normal. The percentiles of the target distributions are shown in the first eight rows of Table 6. Emergent parameters (ECS, TCR, and aerosol forcing ranges) are taken from the IPCC AR6 WG1 Chap. 7 (Forster et al., 2021), and updated climate observations (GMST, OHC, and CO₂ concentrations) are taken from the Indicators of Global Climate Change 2022 (Forster et al., 2023).

Figure 8 Comparison of distributions of key climate metrics (Table 6) in each step of the constraining process. The prior distributions from the 1.6 million member prior ensemble are in blue. The first constraining step using the RMSE comparison to historical temperature is in yellow. The second constraining step that reweights each distribution to its target is in red. The target distribution is in black. The goal is for the red distribution to be as close as possible to the black across all metrics.

The ensemble size in the final reweighted constrained distribution is a user choice. Typically, ensemble sizes of a few hundred to a few thousand are used for projections using reduced-complexity models (Nicholls et al., 2021), which allows for full exploration of the uncertainty space while keeping the number of simulations small enough to allow for efficient computation. For the final posterior distribution in calibrations v1.4.0 and v1.4.1, we select 841 ensemble members from an effective ensemble size of 4356. Moreover, 841 is 1 more than a highly composite number and allows many quantiles of the full distribution to correspond to a single-ensemble member at each point in time.

The posterior ensemble size being one more than a highly composite number is simply an author preference; it is more important to ensure that the posterior is (1) large enough to provide a dense coverage of posterior constraint distributions and (2) small enough that it can provide an unbiased sample size after likelihood weighting. Condition (1) generally imposes a lower bound of around 500 ensemble members, and condition (2) suggests that the effective sample size should be around 5 or more times larger than the target posterior size. If both conditions cannot be simultaneously met, a larger or differently sampled prior or a relaxation of one or more constraints is required.

The evolution of GMST projections from the prior ensemble to the historical RMSE constraint, and finally the reweighted constrained ensemble, is shown in Fig. 9. The prior ensemble allows for a wide range of projections, the majority of which are clearly incompatible with historical GMST (Fig. 9a). The RMSE threshold step, alongside producing historically reasonable projections, substantially reduces the range in projected future warming (Fig. 9b). However, low and particularly very high future warmings pass the historical RMSE constraint. The reweighting step provides a narrower band on historical warming, as well as reducing the spread in future warming further (Fig. 9c). The 5 %–95 % ranges of future warming are similar between the RMSE constraint and the reweighted posterior, but the latter distribution constrains out much of the warm and cool tails of the distribution that passes the RMSE constraint.

Figure 9 Progression of projections using the historical + harmonised SSP2-4.5 emissions for (a) all prior ensemble members, (b) the RMSE <0.17 °C first constraining step, and (c) the final reweighted and constrained posterior. In each plot, progressively darker shaded regions correspond to the minimum–maximum, 5 %–95 %, and 16 %–84 % ranges. The black line is the ensemble median, and the blue line is a historical best-estimate GMST from the Indicators of Global Climate Change 2022 (Forster et al., 2023).

Figure 10 shows the distributions of the 45 parameters used to construct the prior samples (blue histograms) and the reweighted posterior (red histograms). Table S3 lists the parameters and the part of the model that is being affected, as well as its location within the paper. For some distributions, the constraining steps create posteriors that are differently shaped to the priors. Sometimes this is by design. For example, κ₁, the climate feedback parameter, is inversely related to ECS, and the IPCC constraint downweights the likelihood of “hot” combinations (noting that the prior distribution is constructed from CMIP6 models, many of which have higher climate sensitivity than the 95th percentile of 5 °C assessed in IPCC AR6). Occasionally, distributions are multi-modal, such as the parameters that define the ERF_aci shape, due to the model calibrations themselves spanning several orders of magnitude.

Figure 10 Prior (blue) and reweighted posterior (red) distributions of the 45 parameters sampled. For a description of what the parameters correspond to, refer to Table S3.

4 Characteristics of calibrations v1.4.1 and v1.4.0

As a demonstrative case, we show GMST projections for the eight Tier 1 and Tier 2 SSPs using the harmonised emission scenarios in Fig. 11 using calibration v1.4.1. Alongside SSP projections, we use the posterior parameter sets and run concentration-driven runs with a compound 1 % per year CO₂ concentration increase for 140 years. This allows the determination of the airborne fraction of CO₂ at the time of doubling (70 years) and quadrupling (140 years), an estimate of the TCRE obtained at the point of crossing 1000 Gt C of emissions, and a CMIP-consistent approach to calculating TCR (Fig. S1).

For the emission-driven SSP scenarios, the large-scale warming behaviour is in line with expectations, with high-emission scenarios such as SSP5-8.5 and SSP3-7.0 showing several degrees of warming over the next 2 centuries, and lower-emission scenarios warming less. Scenarios where CO₂ emissions turn net negative (SSP1-1.9, SSP1-2.6, and SSP5-3.4 overshoot) show peak and decline behaviour in the ensemble median, though some extreme high-ensemble members continue to warm beyond net zero, owing to a positive zero-emission commitment (Palazzo Corner et al., 2023).

Figure 11 Projections using the weighted posterior for the eight main SSP scenarios. Shaded ranges are (from dark to light) minimum to maximum, 5 %–95 %, and 16 %–84 % of the distribution. Solid lines are distribution medians, and black lines are best-estimate historical warming.

For a more rigorous comparison, we compare the reweighted constrained posterior from fair-calibrate v1.4.1 to the assessed ranges in the AR6 WG1 assessments in Table 6 (see Cross Chapter Box 7.1 in Forster et al., 2021, and Smith et al., 2021b). The first eight rows of the table are the distributions used to reweight the posterior. By design, the fit to the target distribution in these eight cases is very good (in most cases, the “Relative difference” columns in Table 6 are not in bold type). The slight disagreement with the lower bound of the transient climate response is due to the IPCC assessment of the lower end of the very likely range of TCR being lower than the lowest TCR in any of the CMIP6 models which are used to create the prior distribution sample. A better fit to the IPCC-assessed range could be achieved by increasing the samples in the prior TCR distribution at the lower end. The disagreement in the upper bound of ERF_aci is large in relative terms but small in absolute terms. Similarly, no comparison for the upper bound of ERF_ari is provided to avoid division by zero.

Table 6 Comparison of IPCC AR6 WG1 (Forster et al., 2021; Lee et al., 2021; Gulev et al., 2021) or updated (Forster et al., 2023) observational and assessed distributions (“Target” columns), the distributions of the posterior from calibration v1.4.1 (“Reweighted posterior”), and the relative percentage difference. Distributions denoted with an asterisk were assessed as likely ranges in IPCC AR6 WG1, interpreted as ±1 SD (standard deviation), and have been converted to 5 %–95 % ranges here for consistency with other values. Metrics with “Yes” in the “Fit?” column are part of the multiple constraining described in Sect. 3.3.2. Bold text in the “Relative difference” column shows where metrics are more than 5 % from the target for the central estimate and more than 10 % from the target for the upper and lower ranges.

Figure 12 (a) CO₂ emissions in calibration v1.4.1 (GCP 2023 v1.0 up to 2022, harmonised SSP projections after) in solid lines, calibration v1.4.0 (RCMIP v5.1.0) in dotted lines. (b) Median CO₂ concentration projections from v1.4.1, v1.4.0, and CMIP6 (thin lines). The range of 5 %–95 % from v1.4.1 is shown in shaded regions. (c) Median global mean surface temperature projections from calibration v1.4.1 (solid lines), v1.4.0 (dotted lines), and v1.4.0 calibration with historical emissions extended to 2022 under SSP2-4.5 and future scenarios harmonised from 2022 (dashed lines).

The remaining assessed ranges in Table 6 are used for validation- and sense-checking. FaIR under-predicts and provides a narrower range of airborne fraction at 2×CO₂ and 4×CO₂ and TCRE. However, the sensitivities of the carbon cycle feedbacks in FaIR are already well-constrained when comparing the 1750 to 2022 CO₂ emissions with observed concentrations, which places a tight bound on the historical cumulative airborne fraction. The IPCC assessment of airborne fraction is taken from CMIP6 idealised 1pctCO2 runs and is entirely based on the CMIP6 model (Arora et al., 2020), and emission-driven CMIP6 ESMs do not reproduce present-day CO₂ concentrations as tightly as our observational constraint (Lee et al., 2021). In idealised frameworks, TCRE is proportional to the product of airborne fraction and TCR (Jones and Friedlingstein, 2020). The IPCC TCRE assessment is wider than the product of the TCR and airborne fraction individual assessments in quadrature, and as such, distribution fitting to the AR6-assessed ranges of TCR, airborne fraction, and TCRE simultaneously is not possible.

We also compare the emission-driven SSP temperature projections in FaIR to the assessed ranges from the IPCC AR6 WG1 (Lee et al., 2021). For the strong mitigation scenarios SSP1-1.9 and SSP1-2.6, the SSP warming is above the IPCC-assessed ranges, particularly at the 95th percentile. We suggest three reasons. First, concentration (not emission)-driven runs were used to derive the IPCC warming ranges, which excludes the impact of carbon cycle sensitivity uncertainty in a future spread in CO₂ concentrations and thus over-constraining the uncertainty range. In addition, no other line of evidence used by the IPCC for ranges for temperature projections from SSP scenarios included uncertainties in the CO₂ concentrations due to differing carbon cycle feedbacks. Second, the spread in aerosol forcing in our calibration is larger than in CMIP6 (Smith et al., 2020) and the constrained emulator used in the IPCC (Forster et al., 2021). Third and most importantly, the starting point for the future scenario is now 2023 rather than 2015, and emissions have been higher in reality over the last 8 years than in the original SSP1-1.9 and SSP1-2.6 scenarios. The influences of the first and third effects can be visualised by comparing the emissions and projected concentrations of CO₂, and the projected global mean surface temperature anomalies, between v1.4.0 (dotted lines) and v1.4.1 (dashed lines; Fig. 12). Figure 12a also confirms that CO₂ emissions in the recent past can be well-approximated with the SSP2-4.5 scenario.

Conversely, the high-emission SSP3-7.0 and SSP5-8.5 scenarios are projected to warm less in fair-calibrate v1.4.1 compared to the assessments in AR6 WG1 (Fig. 12c). As for the low-emission scenarios, the high-emission scenarios have started to diverge from recent history for CO₂ (Fig. 12a). The emission-driven projections from FaIR tend to result in lower CO₂ concentrations than in the equivalent CMIP6 scenarios (derived using MAGICC6), likely due to the carbon cycle sensitivities being higher in the CMIP6 calibration of MAGICC6 (Fig. 12b). We can also test the influence of different emissions with the same calibration. Figure 12c shows median warming projections from the five main SSPs for the v1.4.0 calibration but with historical emissions updated to 2022 under SSP2-4.5 and other SSPs harmonised from a 2022 start date (dashed lines). Comparing dashed and dotted lines, it can be seen that the higher-emission scenarios are projected to warm less, and lower-emission scenarios warm more, for a 2022 harmonisation compared to SSPs that started in 2015, showing the influence of updating historical simulations for future projections.

We show the comparison to the AR6-assessed ranges for fair-calibrate v1.4.0 in Table S4. In general, these are closer to the IPCC assessments than for v1.4.1, particularly for SSP warming projections, noting that the SSP emissions start in 2015. One reason for the “narrowing” of projections in v1.4.1 (lower scenarios are warmer; higher scenarios are cooler) is the additional 8 years of near-constant CO₂ emissions for the 2015–2022 period in the harmonised scenarios used, reducing the range of climate outcomes in 2100 that are possible with SSP scenarios that satisfy recent historical constraints. One important corollary of this is that median peak warming in the updated harmonised SSP1-1.9 scenario is 1.69 °C in calibration v1.4.1 compared to 1.57 °C in v1.4.0, meaning that is now very unlikely that any realistic mitigation scenario could limit warming to 1.5 °C with no or low overshoot (Dvorak et al., 2022).

5 Conclusions

This paper describes a package, fair-calibrate, that calibrates the responses of the FaIR simple climate model to complex Earth system models, generates a large Monte Carlo ensemble sample, and constrains the results to observations and expert assessments. We claim that a rigorous calibration process that produces ensemble results that are consistent with historically observed climate is a necessary (though not sufficient) condition for trustworthy future climate projections using simple climate models.

We demonstrate two calibrations in this paper: v1.4.1, based on the most up-to-date estimates of all emitted greenhouse gases and short-lived climate forcers, and v1.4.0, which uses emission time series prepared for CMIP6 and AR6 (but are now becoming increasingly outdated). The two different versions presented in this paper produce notably different future projections. The choice of calibration to use depends on user application, and care should be taken to ensure the correct calibration is used for the supplied emissions. Additional calibrations using alternative emission time series and/or constraints can be generated under similar procedures to that described in the paper and accompanying code. Furthermore, the calibration mechanism could be extended to account for different constraints, for example, on TCRE, the zero-emission commitment, warming rates, or future scenario warming. Addition of further constraints should be done with care to ensure internal consistency, particularly when correlated with other constraints, and would likely require a larger prior ensemble size or alternative sampling strategy.

We intend to produce operational updates to fair-calibrate on at least an annual basis. A calibration could be updated based on new climate constraints such as the anticipated yearly updates to Indicators of Global Climate Change (Forster et al., 2023), new source emissions (such as an expected update to CEDS, which will update SLCF emissions to 2022), or new future emission scenarios (such as those from Network for Greening the Financial System). Operationally updated calibrations of emulators and scenarios that reflect the latest scientific knowledge, from which near-future warming can be assessed, will be a beneficial tool in tracking progress towards Paris Agreement aims.