Supplementary Information for FaIRv2.0.0: a generalised impulse-response model for climate uncertainty and future scenario exploration

Abstract. Here we present a Generalised Impulse Response (GIR) model for use in probabilistic future climate and scenario exploration, integrated assessment, policy analysis and teaching. This model is based on a set of only six equations, which correspond to the standard Impulse Response model used for greenhouse gas metric calculations by the IPCC, plus one physically-motivated additional equation to represent state-dependent feedbacks on the response timescales of each greenhouse gas cycle. These six equations are simple and transparent enough to be easily understood and implemented in other models without reliance on the original source code, but flexible enough to reproduce observed well-mixed greenhouse gas (GHG) concentrations and atmospheric lifetimes, best-estimate effective radiative forcing, and temperature response. We describe the assumptions and methods used in selecting the default parameters, but emphasize that other methods would be equally valid: our focus here is on identifying a minimum level of structural complexity. The tunable nature of the model lends it to use as a fully transparent emulator of complex Earth System Models, such as those participating in CMIP6, while also reproducing the behaviour of other simple climate models. We argue that this GIR model is adequate to reproduce the global temperature response to global emissions and effective radiative forcing, and that it should be used as a lowest-common denominator to provide consistency and continuity between different climate assessments. The model design is such that it can be written in tabular data analysis software, such as Excel, increasing the potential user base considerably.


Concentrations ppm ppb ppb ------ppb S1 FaIRv2.0.0 parameter defaults  -    (Meinshausen et al., 2017(Meinshausen et al., ) up to 2014, and fixed at the 2014 level thereafter. These are calculated using the total change in ERF arising from a 1 t emission pulse of each emission type in 2015. S2 Radiative forcing of carbon dioxide, methane and nitrous oxide Here we compare the concentration-forcing relationships of CO 2 , CH 4 and N 2 O used in FaIRv2.0.0, which exclude the interaction terms between gases, to the standard simple formulae detailed in Etminan et al. (2016). Figure Figure S1). This is due to the lack of interaction terms, and results in a maximum absolute error of 0.115 W m −2 at concentrations of CO 2 = 2000 ppm, CH 4 = 3500 ppm, and N 2 O= 525 ppm (the green triangle in the far right-hand 10 side subplot in Figure S1). We believe that in the context of other uncertainties associated with such a high concentration scenario this error is defensible. We note (as was done in Etminan et al.) that the absolute uncertainty in the OLBL calculation is estimated to be 10% for CO 2 and N 2 O and 14% for CH 4 .  (Etminan et al., 2016). We show absolute differences from OLBL approach, grouped by CO2 concentration. Marker size indicates source (Etminan et al. or  Here we illustrate the impact of including an interactive methane lifetime in FaIRv2.0.0. Running RCP8.5  through FaIRv2.0.0-alpha, the evolution of methane lifetime over history and through 2100 is shown in Figure S2. The simulated methane lifetime in 2010 is 8.90 years, compared to 9.15 years in Holmes et al. (2013) from the state-dependent adjustment factor, α(t)). All other non-linearities are due to forcing agents with non-zero f 1 or f 3 parameters. Here we attempt to understand and characterise the non-linearities within FaIRv2.0.0 using a set of pulse-response experiments. We use pulse-response experiments as even in a simple model such as FaIRv2.0.0, the existing gas-cycle and forcing non-linearities can result in relatively complex behaviour. Using SSP2-45 as the reference scenario, we add in emission 25 pulses at the present-day (2019) over several orders of magnitude (0.01 to 1000 GtCO 2 -eq) to the input carbon dioxide and methane emission timeseries in turn. The results of these experiments are shown in Figure S3. In the next paragraph, we step through how to understand these results, but a key message to emphasize is that these non-linearities are not very significant on the whole, and especially on longer (centennial) timescales.
In terms of characterisation, we will step through each variable in turn, explaining the processes behind their behaviour.

30
Throughout this section, we refer to the difference between the experiments and the reference as the anomaly. In Figure   S3, these are normalised by 1/(pulse size). We refer to the difference between the experimental anomalies relative to the smallest pulse experiment, normalised by 1/(pulse size) as the non-linearity anomaly; if the model were linear, the non-linearity anomalies would be zero. For the CO 2 pulse experiments: α CH4 (and hence the CH 4 lifetime) is reduced by a CO 2 pulse due to the overall increase in temperature. The non-linearity 35 anomalies increase as the pulse size increases, mirorring the non-linearity anomalies in temperature.
α CO2 is increased by a CO 2 pulse due to the increased CO 2 uptake and temperature. The non-linearity anomalies increase as the pulse size increases.
-CH 4 concentrations are reduced by a CO 2 pulse due to the reduction in CH 4 lifetime. The non-linearity anomalies increase (ie. the overall anomaly becomes less negative) as the pulse size increases.

40
-CO 2 concentrations are increased by a CO 2 pulse due to the increase in α CO2 . The non-linearity anomalies increase as the pulse size increases, hence CO 2 concentrations are super-linear with emissions. However, the magnitude of the differences due to non-linearities are still small: < 0.5 ppm for pulses lower than 100 GtCO 2 .
-CH 4 forcing anomalies and non-linearity anomalies follow the same behaviour as CH 4 concentrations.
-CO 2 forcing anomalies follow the same behaviour as CH 4 concentration anomalies. However, the non-linearity anoma-45 lies display different behaviour. Since CO 2 forcing is sub-linear with concentrations, the non-linearity anomaly is lower the larger the pulse size (the concentration-forcing sub-linearity dominates the emission-concentration super-linearity).
-Temperature anomalies are, as expected, positive for a CO 2 pulse. However, temperature non-linearity anomalies follow the same behaviour as CO 2 forcing, and are therefore sub-linear with CO 2 emissions. We note that the magnitude of the differences arising due to model non-linearity are < 0.005 K on all timescales, and < 0.001 K on centennial timescales 50 for pulse sizes of < 100 GtCO 2 .
The CH 2 pulse experiments can be analysed in an identical fashion. An additional feature of the behaviour during these experiments is that the sub-linearity of CH 4 forcing with concentrations dominates on very short timescales (meaning that the initial temperature non-linearity anomaly is negative), but for timescales of > 10 years, the super-linearity of CH 4 concentrations with emissions dominates; resulting in an overall long-term super-linear temperature response with CH 4 emissions. Note that 55 even for emission pulses of CH 4 , centennial non-linearities are driven by the carbon-cycle due to the short atmospheric lifetime of CH 4 . However, we again stress that the overall differences arising due to model non-linearity are < 0.05 K on all timescales, and < 0.001 K on centennial timescales for pulse sizes of < 100 GtCO 2 -eq.
Overall, the non-linearities appear insignificant for pulses of less than 10 GtCO 2 -eq. It is important to note that for realstic emission scenarios these non-linearities would grow over time, so over long time periods non-linearities may be important 60 even if the emission rates are significantly lower than 10 GtCO 2 -eq. show anomalies with respect to the reference (no pulse) experiment, scaled by 1/(pulse size). Rows 2 and 4 (marked as "relative") show the scaled anomalies relative to the anomaly for a pulse of size 0.01 GtCO2-eq. The ratio of rows 1 and 2 (or 3 and 4) gives the fractional contribution of the nonlinearity relative to the size of the anomaly for each experiment. Non-linearities are not visible at these scales for pulses of less than 1 GtCO2-eq. S5 CMIP6 data pre-processing CMIP6 data throughout this study is pre-processed as described in Nicholls et al. (2021), using a 21-year-running mean in the piControl experiment to normalise the data and calculate anomalies. The number of ensemble members per model for each experiment is given below in Table S4. in Section 3.1; and in the SSPs displayed in Figure 10.

S6 Energy balance model parameters
Using the method described in Cummins et al. (2020), we tune parameters to CMIP6 models for a 3-box version of the energy balance model described in Geoffroy et al. (2013), with the ocean heat uptake efficacy factor as Winton et al. (2010) included.

S7 Global Warming Index calculation
The Global Warming Index follows the methodology in Haustein et al. (2016), but updates several components. We carry out the calculation using 6 observational warming products (Lenssen et al., 2019;Cowtan and Way, 2014;Vose et al., 2012; 70 Morice et al., 2011;Rohde et al., 2013;Morice et al., 2020), incorporating observational uncertainty either through the provided ensemble product, or if none exists, through the HadCRUT5 ensemble errors. We then generate 5000 realisations of historical ERF as follows. For all forcings excluding aerosol forcing, we take the best-estimate historical timeseries from Smith (2020), and scale them by factors drawn from the distributions detailed in table 6. For aerosol forcing, we take the best-estimate historical timeseries of ERFaci and ERFari, and scale them by scaling factors drawn from a skew-normal distribution that 75 matches the quantiles of the constrained ERFaci and ERFari distributions stated in table 4 of Smith et al. (2020). We consider 18 different response model parameterisations, spanning the ranges of possible realised warming fraction (Millar et al., 2015) and response timescale (Geoffroy et al., 2013). Finally, we include uncertainty due to internal variability through timeseries from the piControl experiment of different CMIP6 models, rejecting models with a too large drift (here "too large" is > ±0.15 K / century); resulting in 102 different representations of internal variability (two random samples per model). Combining these 80 sources of uncertainty gives a 1,836,000,000 member ensemble of the global warming index. For each observational product, we then randomly subsample 500,000,000 members which are used to estimate the distribution of the current level and rate of anthropogenic warming, and thus the FULL ensemble selection probabilities. In this section, we outline the details of how we have implemented the development version of the FaIRv2.0.0 model (FaIRv2.0.0alpha) in python. We note that these notes may not be relevant for all programming languages, but hope that they are useful regardless, especially where a forward difference timestepping scheme is used.
In our implementation, we use a forward difference scheme, in which output variables are assumed to be intra-timestep averages, and auxiliary variables (such as R or G a ) are calculated instantaneously at the start of each timestep. Within each 90 timestep, we carry out the following sequence of computations: 1. Calculate α(t) using equation 3 with variables values from the previous timestep.
2. Calculate C(t) by (a) calculating R i (t) = E(t) * a i * α(t) * τ i * [1 − exp(∆t/(α(t) * τ i ))]/∆t + R(t − 1) * exp(∆t/(α(t) * τ i )) For clarity, we also provide a schematic, Figure S4, indicating at which point in each timestep each variable is assumed to reside in our implementation. In addition to our timestepping scheme, the other key feature of our implementation (and important property of FaIRv2.0.0) is that it is fully vectorised. Variables are represented by -and computations run over -105 multi-dimensional arrays. In our implementation, these arrays have dimensions for (where required): emission/ forcing scenario, gas/forcing parameter set, climate response parameter set, gas species, time, and any additional required dimensions (reservoirs for the carbon cycle or thermal boxes for the climate response). For example, the size of the array containing modelled concentrations in an emission driven run would be: S×P×C×G×t, where S is the number of scenarios, P the number of gas parameter sets, C the number of climate response sets, G the number of gas species and t the total number of timesteps. This 110 ability of FaIRv2.0.0 to be vectorised significantly reduces computation time for array-focused languages such as Python or MatLab. However, an important point to note is that this "complete" vectorisation results in the stored variable arrays becoming large very rapidly (eg. a tenfold increase in each of the S, P and C dimensions would result in a thousandfold increase in the size of the stored arrays, and corresponding memory required). It is therefore important to check that the total size of the stored