Uncertainty in carbon budget estimates due to internal climate variability

Anthropogenic CO₂ emissions increase the atmospheric CO₂ burden, while land and ocean currently take up about half of the emitted anthropogenic CO₂ [14, 15]. Since the total amount of carbon in the atmosphere-ocean-land system is conserved, with the exception of input from fossil fuel emissions, total cumulative fossil fuel CO₂ emissions ($\mathop \smallint \nolimits^ {E_{FF}}$) may be partitioned into the change in atmospheric carbon burden ($\Delta {H_A}$), carbon uptake by the land ($\Delta {H_L}$) and carbon uptake by the ocean ($\Delta {H_O}$), expressed by equation (1).

$$\begin{equation}\mathop \smallint \nolimits^ {E_{FF}}\, = \,\Delta {H_A}\, + \,\Delta {H_O}\, + \,\Delta {H_L}\end{equation} \tag{ 1 }$$

The net ocean carbon uptake is predominantly driven by the increase in atmospheric CO₂ concentration. The net carbon uptake over land is the result of three primary competing processes: (1) the natural carbon cycling through the vegetation-litter-soil continuum, which responds to rising atmospheric CO₂ (through the CO₂ fertilization effect); (2) the effect of associated change in physical climate (e.g. increase in temperature, change in precipitation patterns, etc.); and (3) the carbon loss over land due to land-use and land-cover changes (particularly relevant for historical and future emission scenarios). The change in carbon pool over land can therefore be expressed as the sum of changes due to the natural land carbon pools' response to changes in atmospheric CO₂ and climate ($\Delta {H_{L,natural}}$), and due to anthropogenic land-use change (LUC) processes, where $\int {{E_{LUC}}} {\rm{ }}$ represents the cumulative LUC emissions and includes the effects of the reduced natural land carbon sink due to LUC. The carbon uptake by land ($\Delta {H_L}$) can thus be expressed as: $\Delta {H_L}\, = \,\Delta {H_{L,natural}}\, - \,\mathop \smallint \nolimits^ {E_{LUC}}$ (following Arora et al [16, 17]). Rearranging the terms of the total cumulative emissions (equation (1)) yields:

$$\begin{equation}\smallint {{E_{FF}}} {\mkern 1mu} = {\mkern 1mu} \Delta {H_A}{\rm{ }}{\mkern 1mu} + {\mkern 1mu} \Delta {H_O}{\rm{ }}{\mkern 1mu} + {\mkern 1mu} \Delta {H_{L,natural}}{\mkern 1mu} - \smallint {{E_{LUC}}} {\mkern 1mu} \end{equation} \tag{ 2 }$$

$$\begin{align}\small \smallint E \! =\! {\rm{ }}\smallint {{E_{FF}}} \! +\! {\rm{ }}\smallint {{E_{LUC}}{\rm{ }}} {\mkern 1mu} \! = \!{\mkern 1mu} \Delta {H_A}{\mkern 1mu} \! +\! {\mkern 1mu} \Delta {H_O}{\mkern 1mu} + {\mkern 1mu} \Delta {H_{L,natural}}\end{align} \tag{ 3 }$$

The term $\Delta {H_{L,natural}}$ is the terrestrial sink. Most Earth System Models (ESMs) simulate processes related to LUC interactively, and it is not possible to diagnose $\int {{E_{LUC}}} {\rm{ }}$ without another simulation that does not include anthropogenic LUC. $\mathop \smallint \nolimits^ {E_{LUC}}$ can typically be diagnosed by differencing $\Delta {H_L}$ from simulations with and without LUC [18].

In the absence of such simulations with and without LUC for the models considered here, we use an estimate of the values of CO₂ emissions from LUC, ${E_{LUC}}$, for the historical period and for the respective future SSP and RCP scenarios using estimates from [19], as was done in earlier studies [5, 20, 21]. We acknowledge that this estimate of LUC emissions differs from the actual LUC emissions that each model generates. However, since the primary purpose of this study is to quantify the role of internal variability on carbon budgets rather than to provide an exact estimate of the remaining carbon budget, the use of specified $\mathop \smallint \nolimits^ {E_{LUC}}$ is not expected to influence our conclusions. Note that LUC emissions are not relevant for estimating TCRE in CO₂-only simulations, as explained in section 2.2.

The relationship between CO₂-induced warming and cumulative CO₂ emissions has shown to be approximately linear for up to 2000 PgC emitted [1, 5, 22–24] and beyond [25, 26], and largely independent of CO₂ emission pathway for a wide range of CO₂ -only pathways [27, 28]. The linearity of this relationship breaks down only for very low or very high emission rates [29].

The TCRE [1, 5] is expressed as a ratio of CO₂-induced global mean warming ($\Delta {T_{C{O_2}}}$) to total cumulative CO₂ emissions ($\mathop \smallint \nolimits^ E$), expressed by equation (4). Cumulative CO₂ emissions are inferred as the sum of all right-hand-side terms in equation (2), since the simulations considered here are driven by CO₂-concentrations (as opposed to driven by CO₂ emissions). We calculated $\Delta {H_L}{\text{ }}$ and $\Delta {H_O}{\text{ }}$ as the time-integrated atmosphere-land and atmosphere-ocean carbon fluxes ('nbp' and 'fgco2' CMIP variable names). We use the global mean surface air temperature with full global coverage (GSAT), that is directly diagnosed from the ESM output ('tas' CMIP variable).

$$\begin{equation}TCRE{\mkern 1mu} = {\mkern 1mu} \frac{{\Delta {T_{C{O_2}}}}}{{\int {E{\rm{ }}dt} }}{\mkern 1mu} = {\mkern 1mu} \frac{{\Delta {T_{C{O_2}}}}}{{\Delta {H_A}{\mkern 1mu} + {\mkern 1mu} \Delta {H_L}{\mkern 1mu} + {\mkern 1mu} \Delta {H_O}}}\end{equation} \tag{ 4 }$$

TCRE (equation (4)) was originally defined [1] using 1pctCO₂ concentration-driven simulations in which atmospheric CO₂ increases at a constant rate of 1% per year. Concentrations of all other non-CO₂ climate forcings and land cover in the 1pctCO₂ simulation stay at their pre-industrial level and therefore, emissions from LUC are zero (${E_{LUC}}$; section 2.1).

In addition to 1pctCO₂ simulations, we also consider all-forcing simulations, which include both CO₂ and non-CO₂ forcing agents, in addition to natural forcing (solar and volcanic), and anthropogenic LUC (equation (3)). Warming as a function of cumulative emissions computed from such all-forcing simulations is referred to as the effective TCRE [30], and the presence of non-CO₂ forcing makes it non-linear and dependent on the pathway of non-CO₂ emissions [21, 31, 32]. Furthermore, emissions from LUC differ among such scenarios and are difficult to diagnose from the standard model output (see section 2.1). As noted in section 2.3, it is no longer recommended [3] to estimate remaining carbon budgets from the effective TCRE in RCP and SSP scenarios.

A remaining carbon budget can be calculated from an estimate of the anthropogenic warming to date (discussed in the next section 2.4), and either (i) the ratio of CO₂-induced warming to cumulative emissions, known as the TCRE, in addition to information on the temperature response to the future evolution of non-CO₂ emissions; or (ii) climate model scenario simulations that reach a given temperature threshold.

2.3.1. Approach (i): SR 1.5 framework for estimating remaining carbon budgets

2.3.2. Approach (ii): threshold exceedance budgets (IPCC AR5 and studies following it)

Paris Agreement temperature limits may be interpreted as levels of anthropogenic warming [5, 6], that by definition, are not subject to internal variability. To estimate the anthropogenic contribution to the observed warming, the [8] AWI can be used, by regressing out the natural influence. However, the AWI estimate is affected by internal variability in the observations. The uncertainty range due to internal in AWI can be estimated based on pre-industrial control variability, which is based on pre-industrial control simulations of CMIP5 models [8, 30]. Either method of calculating the remaining carbon budgets from the present-level day warming (either the SR1.5 approach or the previously used threshold exceedance approach) requires an estimate of the anthropogenic warming to-date, in order to know how much warming is left until the 1.5 C or 2.0 °C target, for example. Remaining carbon budgets are more accurate than the total carbon budgets, as using the present-day baseline reduces the overall uncertainty in carbon budget estimates [33].

We use simulations from climate models participating in the CMIP5 and CMIP6 intercomparisons [36, 37] that had ensemble sizes of at least five members available in all-forcing simulations, and single ensemble members from 1pctCO₂ simulations (listed in Supplementary table S1 (available online at stacks.iop.org/ERL/15/104064/mmedia)). The three models with the largest ensemble sizes are: the Canadian Earth System Model (CanESM2) with a 50-member ensemble [16, 38]; CanESM5 with a 50-member ensemble [39], and the MPI-ESM-GE grand ensemble of 100 members [40] (note that only 32 members of MPI-ESM-GE had carbon cycle output available in 1pctCO₂ simulations used to calculate TCRE). In each model's ensemble, the differences among simulations are only due to internal variability of the climate system. Such variations can be simulated by perturbing a random seed within a stochastic component of the cloud parametrisation (in the case of CanESM2) [38], and/or by choosing a different initialization year from the pre-industrial control run from which a member of the historical simulation is initialized (in the case of MPI-ESM-GE) [40]. Each model simulation was driven with prescribed atmospheric CO₂ concentration according to the historical forcing followed by the specified-concentration RCP 8.5 scenario [12] (in the case of CMIP5 models), or SSP 5–8.5 and SSP 2–4.5 scenarios [13] (in the case of CMIP6 models; Supplementary tables S1 and S2). These simulations include natural and anthropogenic forcing from CO₂ and non-CO₂ forcers [19]. We also use the 1pctCO₂ simulations from models which contributed to CMIP6. These 1pctCO₂ simulations were used to calculate TCRE values for the following models: CanESM5, CESM2, CNRM-ESM2-1, IPSL-CM6A-LR, MPI-ESM1-2-LR, NorESM2-LM, UKESM1-0-LL, ACCESS-ESM1-5, CESM2-WACCM, MIROC-ES2L (see Supplementary Material for a list of all models used for each type of simulation, and their respective ensemble sizes).

Internal variability in the climate system can be quantified as uncertainty across the individual ensemble members in a large ensemble of simulations from each model, because each ensemble member is subject to identical forced radiative forcing (see Methods section 2.5). An estimate of the forced response (with internal variability removed) can be obtained by taking a mean of individual ensemble members of the same model. However, most models only provide one ensemble member in 1pctCO₂ simulations, which is subject to internal variability. Some studies [41, 42] suggest that fitting a fourth-order polynomial to annual temperature time-series results in a reasonably good approximation of the forced response. Here we compare how different smoothing methods applied to the temperature response in 1pctCO₂ simulations affect the resulting TCRE estimates.

The TCRE estimate obtained by fitting a fourth-order polynomial to individual ensemble members of the same model results in nearly the same value as taking the 32-member ensemble mean annual mean temperature in MPI-ESM-GE (figures 1(a) and (b)) and similar results are obtained for the 10-member ensemble of CNRM-ESM 2.1 (figures 1(c) and (d)). However, in both models, the likely range (17%–83%) and 5%–95% range is much narrower when the smoothing is applied, compared with the annual time-series (figures 1(b) and (d), first yellow bar compared with the last purple bar). The remaining uncertainty in TCRE calculated using a fourth-order polynomial results from variability on longer, multi-decadal time-scales that is not accounted for by the polynomial.

Zoom In Zoom Out Reset image size

The uncertainty in TCRE due to internal variability in individual models can be as large as ±0.1 °C/1000 PgC (5%–95% range) when temporal smoothing is applied in 1pctCO₂ simulations in a single model (e.g. MPI-ESM-GE; figure 1(a) and (b), or CNRM-ESM2-1 in figure 1(c) and (d); see Supplementary table S2 for details). We note that a similarly narrow TCRE range can also be obtained by estimating TCRE as the 20-yr mean, centered on the year when 1000 PgC is emitted, which is the recommended way of obtaining the TCRE estimate [24]. The method of calculating TCRE does not influence the overall TCRE distribution in CMIP6 models appreciably, as model uncertainty dominates over the uncertainty due to internal variability and the method to calculate TCRE (figure 1(e)).

In principle, TCRE as a property of climate models based on the 1pctCO₂ simulations, and the value of TCRE should be determined by the externally forced model response (free from the unforced internal variability). In practice, there are currently limitations to estimating TCRE without the influence of unforced internal variability as illustrated here, because often only one ensemble member is available per model in 1pctCO₂ simulations. For similar analysis regarding the Transient Climate Response (TCR) estimates (i.e. warming in the year 70 of 1pctCO₂ simulations) see Supplementary Material.

In this section, we investigate the influence of internal variability on remaining carbon budgets using the threshold exceedance approach (i) (Methods section 2.3). Despite smoothing to reduce the uncertainty, we still find a significant contribution to overall uncertainty from internal variability, especially for 1.5 °C budgets. The uncertainty in warming due to internal variability as a function of cumulative emissions (figure 1) implies that even within the same model and under the same scenario, the remaining carbon budget for a given level of warming (such as 1.5 °C) is not a single number, but rather a range of numbers, even after strong smoothing is applied to isolate the forced response. This remaining uncertainty in each model response arises due to internal variability in the climate system that acts on various time scales from multiple decades to multiple centuries, and is thus still present even once the inter-annual to decadal variability is removed. This uncertainty can be reduced but not eliminated entirely given current observational and model capacities. We quantify this uncertainty in the remaining carbon budgets in figure 2.

Zoom In Zoom Out Reset image size

Here, for illustration of the effect of internal variably in all-forcing scenarios, we calculate remaining carbon budgets for 1.5 °C and 2.0 °C, with respect to the 2006–2015 baseline, directly from the effective TCRE curve in SSP2-4.5, SSP5-8.5 and RCP 8.5 scenarios (figure 2), and reflect the amount of cumulative CO₂ emissions in the year prior to which a given temperature target is exceeded (i.e. 'exceedance budgets' [11]). As noted in Methods section 2.3, this approach of calculating remaining carbon budgets is no longer recommended, as such threshold exceedance budgets calculated directly from ESM output under RCP and SSP scenarios are subject to assumptions on the particular mix of CO₂ and non-CO₂ forcings that have limited application for mitigation pathways in line with keeping warming to 1.5 °C or 2 °C. Recently, other methods (using TCRE values directly and estimates of future warming from non-CO₂ forcing contributions) have been suggested to circumvent this issue (See Methods section 2.3 and Rogelj et al [3]).

The focus of our study is to analyse the uncertainty in the remaining carbon budgets due to internal variability, rather than their nominal values, and carbon budgets calculated directly from ESM output provide such an opportunity. Note that historical and future emission scenarios (such as SSPs and RCPs) are subject to anthropogenic and natural forcing (solar forcing and volcanic forcing), in addition to the unforced climate variability. We make the assumption that the ensemble mean of the simulations in a large ensemble including both anthropogenic and natural forcing is an estimate of the anthropogenic response (i.e. the natural response is assumed to be zero). While the natural forcing response is probably relatively small in both cases, it will not be exactly zero. This may affect the threshold exceedance carbon budgets, but not the carbon budgets calculated from TCRE directly (as TCRE is obtained from 1pctCO₂ simulations; Methods section 2.2).

Internal variability affects the remaining carbon budgets for both the 1.5 °C and 2.0 °C limits by about ±30 PgC and ±40 PgC (5%–95% range for both sets of scenarios considered here; Supplementary table S3; figures 2(b) and (c)) even after estimating the forced temperature response by fitting fourth-order polynomials. Generally, differences among models contribute much more to the overall uncertainty than the differences due to internal variability. The influence of inter-model uncertainty can be minimized by adjusting the baseline period to the most recent period and calculating remaining budgets only from that point onwards and until a given temperature limit is reached [33]. In turn, the uncertainty due to internal variability is approximately constant and not affected by this change of baseline (since uncertainty in TCRE and remaining carbon budgets is dominated by the temperature component, and internal variability in temperature does not change much at different forcing levels [42]).

In figure 3 we illustrate relative contributions of internal variability and model uncertainty, for the models considered here, following Hawkins and Sutton [41, 42] approach of visualizing contributions of different uncertainty types. The relative contribution from internal variability is important and may exceed the model uncertainty on very short time-scales (as it is constant, irrespective of the baseline), but on longer time-scales, the inter-model spread dominates. Since the remaining carbon budget for 1.5 °C is close to the origin of the effective TCRE curve, with the 2006–2015 baseline (figure 2(a)), it is subject to relatively large uncertainty due to internal variability (figure 3(c)), as model uncertainty is assumed to be zero at the origin (by design due to the choice of a recent baseline). This relative contribution of internal variability to overall uncertainty diminishes for a 2.0 °C carbon budget (figure 3(d)), where model uncertainty starts to emerge and dominate.

Zoom In Zoom Out Reset image size

We note that figure 3 is for illustrative purposes only because only models for which large or medium-sized ensembles were considered. We also do not consider scenario uncertainty here, as for the time period and lead time considered there is no substantial difference in the effective TCRE between SSP5-8.5 and SSP2-4.5 (for the models considered here; figure 3(a)), and between RCP 8.5 and RCP 4.5 (Tokarska and Gillett [20]). However, the SSP and RCP scenarios represent only a small portion of all available scenarios, and for 1.5 °C and 2.0 °C targets, there are numerous combinations of non-CO₂ and CO₂ forcings to reach those targets. Therefore, ambitious mitigation targets (such as the Paris Agreement) are subject to large scenario uncertainty that would need to be explored beyond the RCP and SSP scenarios (for example, using climate model emulators and the SR1.5 scenario database [34, 35]). Such analysis is beyond the scope of this study that focuses primarily on the role of internal variability on TCRE and remaining carbon budgets, rather than an assessment of the overall uncertainty.

Estimating the remaining budget using either method -the SR1.5 approach (i) or the remaining threshold exceedance budgets approach (ii) based on AR5 requires a precise estimate of the present-day anthropogenic warming or mean warming at the time of the baseline period. The [9] AWI, is a method of assessing anthropogenic global mean warming for meeting the Paris Agreement target [3, 4] (Methods section 2.4). We qualitatively compare the uncertainty due to internal variability in AWI with the spread of the large ensemble (note that the definition of internal variability is not exactly the same between the two estimates, and AWI has the natural-only variability from solar and volcanic forcing influence removed [9]). Following Haustein et al [9], figure 4 shows the AWI calculated on a single ensemble member from the MPI-ESM-GE (black), which is treated as pseudo-observations. The resulting AWI and its overall uncertainty range is shown in orange, while the MPI-ESM-GE grand ensemble mean of 100 simulations is shown in dark blue. The spread due to internal variability in AWI (red shaded area) is similar to the spread of the 100 members of MPI-ESM-GE large ensemble (light blue lines), where individual ensemble members were smoothed using the fourth-order polynomial fits. The magnitude of the uncertainty due to internal variability (at 'present-day', i.e. by the end of these time-series, figure 4(b)) in both these methods is similar (though these uncertainty ranges are not a like-for-like comparison, as discussed above). Generally, the large ensembles approach estimate provides a reasonable way of characterising internal variability in a particular model. The uncertainty in AWI due to internal variability in the year 2019 is estimated at ±0.09 °C (figure 4 red shaded area; 5%–95% range), which corresponds to a contribution to uncertainty in remaining carbon budgets of about ±50 PgC (using a median estimate of TCRE of 1.65 °C/1000 PgC as in SR1.5 Rogelj et al [4]).

Zoom In Zoom Out Reset image size