Statistical estimation of global surface temperature response to forcing under the assumption of temporal scaling

Reliable quantification of the global mean surface temperature (GMST) response to radiative forcing is essential for assessing the risk of dangerous anthropogenic climate change. We present the statistical foundations for an observation-based approach using a stochastic linear response model that is consistent with the long-range temporal dependence observed in global temperature variability. We have incorporated the model in a latent Gaussian modeling framework, which allows for the use of integrated nested Laplace approximations (INLAs) to perform full Bayesian analysis. As examples of applications, we estimate the GMST response to forcing from historical data and compute temperature trajectories under the Representative Concentration Pathways (RCPs) for future greenhouse gas forcing. For historic runs in the Model Intercomparison Project Phase 5 (CMIP5) ensemble, we estimate response functions and demonstrate that one can infer the transient climate response (TCR) from the instrumental temperature record. We illustrate the effect of long-range dependence by comparing the results with those obtained from one-box and two-box energy balance models. The software developed to perform the given analyses is publicly available as the R package INLA.climate.

1 Introduction

Despite decades of research and development regarding global circulation models (GCMs) and Earth system models (ESMs), the discrepancies between models remain substantial, even as we describe physical processes with increasing accuracy and resolution. Part of the model spread is associated with a lack of understanding of the shortwave cloud feedback (Qu et al., 2018). However, there are several other modeling choices and compromises that contribute to the uncertainty (Flato, 2011). As a consequence, several studies have focused on constraining model results on climate sensitivity on observational data; see, e.g., the work of Annan and Hargreaves (2006) or the more recent studies of Cox et al. (2018) and Rypdal et al. (2018b, a). These studies focus on the equilibrium climate sensitivity (ECS) as an essential metric of the climate response, as have numerous paleoclimate studies (Hansen et al., 2013; von der Heydt and Ashwin, 2017; Köhler et al., 2017).

A simpler approach is to adopt a linear approximation and to apply statistical methods to extract information on the climate response from data on global surface temperature and radiative forcing in the instrumental era. Under the assumption of a linear and stationary response, the global surface temperature anomaly ΔT can be expressed as a filtering of the global radiative forcing F, mathematically expressed as

$$\begin{array}{}\text{(1)}& \mathrm{\Delta }T\left(t\right)=\underset{-\mathrm{\infty }}{\overset{t}{\int }}G(t-s)\left(F\left(s\right)\mathrm{d}s+\mathit{\sigma }\mathrm{d}B\left(s\right)\right),\end{array}$$

where σdB(t) represents a white-noise forcing that gives rise to internal climate variability, and G is the response function, or Green's function, characterizing the relation between forcing and the temperature anomaly. A model of this form can arise from the simplest energy balance model, i.e., the equations

$$\begin{array}{}\text{(2)}& {\displaystyle \frac{\mathrm{d}\mathrm{\Delta }Q}{\mathrm{d}t}}=-\mathit{\lambda }\mathrm{\Delta }T+F\end{array}$$

and ΔQ=CΔT, where ΔQ is the change in the system's heat content corresponding to a temperature change ΔT, and C is heat capacity (Rypdal, 2012). If a white-noise forcing term is included on the right-hand side of Eq. (2) it becomes a stochastic differential equation with a stationary solution to the form of Eq. (1), with $G\left(t\right)=C^{-\mathrm{1}}{e}^{-t/\mathit{\tau }}$ and $\mathit{\tau }=C/\mathit{\lambda }$. The process has a natural decomposition into the response to the known forcing,

$$\begin{array}{}\text{(3)}& \mathrm{\Delta }{T}_{\text{det}}\left(t\right)={\displaystyle \frac{\mathrm{1}}{C}}\underset{-\mathrm{\infty }}{\overset{t}{\int }}{e}^{-(t-s)/\mathit{\tau }}F\left(s\right),\mathrm{d}s\end{array}$$

and a stochastic term,

$$\begin{array}{}\text{(4)}& X\left(t\right)={\displaystyle \frac{\mathit{\sigma }}{C}}\underset{-\mathrm{\infty }}{\overset{t}{\int }}{e}^{-(t-s)/\mathit{\tau }}\mathrm{d}B\left(s\right),\end{array}$$

which for this particular model is an Ornstein–Uhlenbeck process. Rypdal and Rypdal (2014) show how the parameters of the two terms can be estimated simultaneously from time series of forcing and the GMST using the maximum likelihood (ML) method. They also demonstrate that the resulting process is inconsistent with observations. The stochastic term X(t) does not exhibit the strong positive decadal-scale serial correlations that are observed in the GMST in the instrumental era, and the model's response to reconstructed forcing for the last millennium does not show sufficient low-frequency variability compared to Northern Hemisphere temperature reconstructions.

The inconsistency of the simple energy balance model is due to the slow climate response associated with the energy exchange with the deep ocean. One can easily incorporate this effect within the framework of Eq. (1) by generalizing the zero-dimensional (one-box) model to a two-box model that includes a layer representing the deep ocean (Geoffroy et al., 2013; Held et al., 2010; Caldeira and Myhrvold, 2013) or the more general m-box model discussed by Fredriksen and Rypdal (2017).

The generalization from the zero-dimensional (one-box) energy balance model to the two-box model, or the m-box models, means that the number of free parameters increases. Concerning statistical inference, this could be problematic due to potential over-fitting. Mathematically, the generalization of Eq. (2) is of the form

$$\begin{array}{}\text{(5)}& \mathbf{C}{\displaystyle \frac{\text{d}\mathrm{\Delta }\mathit{T}\left(t\right)}{\text{d}t}}=\mathbf{K}\mathrm{\Delta }\mathit{T}\left(t\right)+\mathit{F}\left(t\right),\end{array}$$

where the diagonal elements of C are the heat capacities of each box, and the matrix K contains coefficients describing heat exchange between boxes and the feedback parameter λ.

The system in Eq. (5) is solved by bringing the matrix C⁻¹K to diagonal form, and the surface temperature anomaly can be written as in Eq. (1), where G(t) is the weighted sum of m exponential functions (Fredriksen and Rypdal, 2017):

$$\begin{array}{}\text{(6)}& G\left(t\right)=\left\{\begin{array}{ll}{\sum }_{k=\mathrm{1}}^m{w}_k{e}^{-t/{\mathit{\tau }}_k},& t\ge \mathrm{0}\\ \mathrm{0},& t<\mathrm{0}\end{array}\right..\end{array}$$

The characteristic timescales ${\mathit{\tau }}_k=-\mathrm{1}/{\mathit{\mu }}_k$ are defined from the eigenvalues μ_k of C⁻¹K, and w_k denotes the weight of the kth exponential function.

On the other hand, global temperature variability exhibits a scaling symmetry. For instance, both the forced and the unforced global temperature variability have power spectral densities (PSDs) that are approximate power laws,

$$\begin{array}{}\text{(7)}& S\left(f\right)\sim f^{-\mathit{\beta }},\end{array}$$

for frequencies corresponding to timescales ranging from months to centuries (Rypdal and Rypdal, 2016; Rybski et al., 2006; Lovejoy and Schertzer, 2013; Huybers and Curry, 2005; Franzke, 2010; Fredriksen and Rypdal, 2016). The global temperature fluctuations are consistent with a fractional Gaussian noise (fGn), which can formally be defined by the integral analogous to Eq. (4), but with the exponential response function replaced with a scale-invariant response function:

$$\begin{array}{}\text{(8)}& G\left(t\right)={\left({\displaystyle \frac{t}{\mathit{\gamma }}}\right)}^{\mathit{\beta }/\mathrm{2}-\mathrm{1}}\mathit{\xi }.\end{array}$$

Here, γ is a scale parameter with the dimension of time, and ξ is a variable needed in order for G(t) to have the correct physical dimensions. The scaling exponent β (defined from the PSD in Eq. 7) relates to the so-called Hurst exponent of the fGn via the formula $\mathit{\beta }=\mathrm{2}H-\mathrm{1}$. Based on this Rypdal and Rypdal (2014) proposed a fractional linear response model in the form of Eq. (1), in which the parsimonious expression in Eq. (8) replaces the linear combination of exponential functions in Eq. (6). The cost of the reduction in model complexity is that the fractional linear response model does not have a well-defined ECS, and in general, we cannot write the model as a system of differential equations as in Eq. (5). But on timescales up to approximately 10³ years, the model provides an accurate description of both forced and unforced surface temperature fluctuations (Rypdal and Rypdal, 2014; Rypdal et al., 2015), and the millennial-scale climate sensitivity in the estimated fractional linear response model correlates strongly with ECS over the ensemble of models in the Coupled Model Intercomparison Project Phase 5 (CMIP5) (Rypdal et al., 2018a).

Temporal-scale invariance in global temperature fluctuations is an empirical observation and we cannot deduce the parameters in the fractional linear response model from physical principles. This paper presents a statistical methodology that makes it possible to fit the model to observational data and estimate all model parameters. Parameter estimation is done within a Bayesian framework by making use of the methodology of integrated nested Laplace approximation (INLA) for latent Gaussian models introduced in Rue et al. (2009). Barboza et al. (2019) use this framework to investigate model formulations and forcing components in paleoclimate reconstructions.

The INLA methodology and inference for our statistical model, assuming the scale-invariant response function in Eq. (8), are described further in Sect. 2. This section explains how to compute the marginal posterior distributions of the model parameters. As the model has a nonstandard form, this includes certain modifications of the INLA methodology to ensure computational efficiency. We discuss applications in Sect. 3.

In Sect. 3.1 we fit the model to the temperature and forcing dataset generated by the GISS-E2-R (Schmidt et al., 2014) ESM. Here we show how to extract the GMST response to the known forcing using a Monte Carlo sampling approach. In Sect. 3.2, the model is used for temperature forecasting wherein the Representative Concentration Pathway (RCP) trajectories describe the future CO₂ forcing. Section 3.3 describes how the transient climate response (TCR) can be estimated using our model. We obtain estimates for 19 temperature series and their associated adjusted forcing series.

We compare the resulting estimates of TCRs with the TCRs obtained directly from the respective ESMs and with TCR estimates from the historical HadCRUT4 temperature dataset using different forcing data. The applications are incorporated in the R package INLA.climate. This package also includes the option of using the exponential response functions as defined by Eqs. (3) and (4) for the one-box model and Eq. (6) for m-box models. A discussion and final conclusions are given in Sect. 4.

2 Discrete-time modeling and statistical inference

Rypdal and Rypdal (2014) use an ML estimator to estimate the model parameters from the observational yearly time series of $\mathbf{\Delta }\mathit{T}=(\mathrm{\Delta }{T}_{\mathrm{1}},\mathrm{\dots },\mathrm{\Delta }{T}_n)$ of GMST and the corresponding vector of radiative forcing $\mathit{F}=(F_{\mathrm{1}},\mathrm{\dots },F_n)$. Here, we estimate parameters by adopting a Bayesian framework, making use of the INLA methodology (Rue et al., 2009, 2017). This approach implies that parameters are treated as stochastic variables and assigned prior distributions. The information given by the priors is then combined with the likelihood of the observations and updated to give posterior distributions using Bayes' theorem.

In a discrete-time model, we assume that ΔT_t has a Gaussian distribution with a random mean expressed by the linear predictor

$$\begin{array}{}\text{(9)}& {\mathit{\eta }}_t={\mathit{\sigma }}_{\mathrm{f}}\sum _{s=\mathrm{1}}^t{G}_{ts}\left(H\right)\left(F_s+F_{\mathrm{0}}\right)+{\mathit{\epsilon }}_t,t=\mathrm{1},\mathrm{\dots },n,\end{array}$$

where ${\mathit{\sigma }}_{\mathrm{f}}={\mathit{\gamma }}^{-\mathit{\beta }/\mathrm{2}+\mathrm{1}}$, while F₀ denotes a shift parameter that gives the initial forcing value. G_ts denotes a discretely indexed element of the function,

$$\begin{array}{}\text{(10)}& G_{ts}\left(H\right)=\left\{\begin{array}{ll}{\left(t-s+\frac{\mathrm{1}}{\mathrm{2}}\right)}^{H-\frac{\mathrm{3}}{\mathrm{2}}},& \mathrm{1}\le s\le t\le n\\ \mathrm{0},& \text{otherwise}\end{array}\right..\end{array}$$

Further, the vector $\mathit{\epsilon }=({\mathit{\epsilon }}_{\mathrm{1}},\mathrm{\dots },{\mathit{\epsilon }}_n)$ denotes a zero-mean fGn process, implying that the covariance between ε_t and ε_s is

$$\begin{array}{}\text{(11)}& \begin{array}{rl}{\mathrm{\Sigma }}_{ts}& ={\displaystyle \frac{{\mathit{\sigma }}_{\mathit{\epsilon }}^{\mathrm{2}}}{\mathrm{2}}}\left(|t-s+\mathrm{1}{|}^{\mathrm{2}H}+|t-s-\mathrm{1}{|}^{\mathrm{2}H}\right.\\ & \left.-\mathrm{2}|t-s{|}^{\mathrm{2}H}\right),t,s=\mathrm{1},\mathrm{\dots },n,\end{array}\end{array}$$

where σ_ε=σσ_f. In vector notation, the predictor is then given by

$$\begin{array}{}\text{(12)}& \mathit{\eta }=\mathit{\mu }+\mathit{\epsilon },\end{array}$$

where

$$\begin{array}{}\text{(13)}& \mathit{\mu }=\mathit{\mu }(H,{\mathit{\sigma }}_{\mathrm{f}},F_{\mathrm{0}})={\mathit{\sigma }}_{\mathrm{f}}\mathbf{G}\left(H\right)\left(\mathit{F}+F_{\mathrm{0}}\right).\end{array}$$

The covariance matrix of the predictor is $\mathbf{\Sigma }=\mathbf{\Sigma }(H,{\mathit{\sigma }}_{\mathit{ϵ}})$ with the elements in Eq. (11). Notice that the matrix G(H) is lower triangular with elements given by Eq. (10). The given formulation implies that the vector μ represents the GMST response to the known forcing F, while ε is the GMST response to the random forcing, i.e., the unforced climate variability.

The statistical regression formulation in Eq. (9) has a hierarchical structure in which the expected temperature anomalies are modeled in terms of the random predictor η with elements specified by Eq. (9). The predictor depends on additional model parameters $\mathit{\theta }=(H,{\mathit{\sigma }}_{\mathit{\epsilon }},{\mathit{\sigma }}_{\mathrm{f}},F_{\mathrm{0}})$. This setup implies that we need to assign priors to both the predictor and the model parameters. By assigning a Gaussian prior to η, the resulting model becomes a latent Gaussian model, which can be analyzed using the INLA methodology. In general, this class of models introduces a latent Gaussian field x, which contains all the random components of a linear predictor, including the predictor itself. In our case, the latent field is equal to the linear predictor, $\mathit{x}=\mathit{\eta }=\mathit{\mu }+\mathit{\epsilon }.$ However, inference for this model is not straightforward as the model parameter H appears in both of the terms μ and ε. We choose to circumvent this problem by considering the sum μ+ε to be a single model component, i.e., as a fractional Gaussian noise process with mean vector μ and covariance matrix Σ. The dependence between two components implies that we will not get separate posterior estimates for μ and ε directly.

Using p(⋅) as a generic notation for probability density functions, we can summarize the three-stage hierarchical structure of latent Gaussian models, including distributional assumptions, as follows.

• The first stage specifies the likelihood of the model. The observed temperature anomaly ΔT_t is assigned a Gaussian distribution with negligible fixed variance and mean η_t. The observations are assumed to be conditionally independent given the latent field x and parameters θ, i.e.,

  $$p\left(\mathrm{\Delta }\mathit{T}\right|\mathit{x},\mathit{\theta })=\prod _{t=\mathrm{1}}^{n}p(\mathrm{\Delta }{T}_t|x_t,\mathit{\theta }).$$

• The second stage specifies the prior distribution for the latent field. Given the parameters θ, the latent field x is assigned a Gaussian prior distribution with mean vector ${\mathit{\mu }}_x=\text{E}\left[\mathit{x}\right|\mathit{\theta }]$ and precision matrix, $\mathbf{Q}=\mathbf{Q}(H,{\mathit{\sigma }}_{\mathit{\epsilon }})$, defined as the inverse covariance matrix, i.e.,

  $$p\left(\mathit{x}\right|\mathit{\theta })=\sqrt{{\displaystyle \frac{det\mathbf{Q}}{(\mathrm{2}\mathit{\pi }{)}^n}}}\exp \left(-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\left(\mathit{x}-{\mathit{\mu }}_x\right)}^T\mathbf{Q}\left(\mathit{x}-{\mathit{\mu }}_x\right)\right).$$

• The third stage specifies independent priors for the parameters:

  $$p\left(\mathit{\theta }\right)=p\left(H\right)p\left({\mathit{\sigma }}_{\mathit{\epsilon }}\right)p\left({\mathit{\sigma }}_{\mathrm{f}}\right)p\left(F_{\mathrm{0}}\right).$$

  The shift parameter F₀ is assigned a zero-mean Gaussian prior, while the other parameters are assigned penalized complexity (PC) priors (Simpson et al., 2017). The class of PC priors represents a recently developed framework to compute priors based on specific principles, including support for Occam's razor. The PC prior of the two scaling parameters σ_f and σ_ϵ can be computed to equal the exponential distribution, while the PC prior of H is computed numerically (Sørbye and Rue, 2018).

The joint posterior for all components of the latent field and all of the model parameters is then summarized by

$$p(\mathit{x},\mathit{\theta }|\mathrm{\Delta }\mathit{T})\propto \prod _{t=\mathrm{1}}^{n}p(\mathrm{\Delta }{T}_t|x_t,\mathit{\theta })p\left(\mathit{x}\right|\mathit{\theta }\left)p\right(\mathit{\theta }).$$

Our main objective is to estimate the marginal posterior distribution for all components of the latent field

$$\begin{array}{}\text{(14)}& \begin{array}{rl}& p\left(x_t\right|\mathrm{\Delta }\mathit{T})=\int p(x_t|\mathit{\theta },\mathrm{\Delta }\mathit{T})p\left(\mathit{\theta }\right|\mathrm{\Delta }\mathit{T})\text{d}\mathit{\theta },\\ & t=\mathrm{1},\mathrm{\dots },n\end{array}\end{array}$$

and the marginal posteriors for all the model parameters

$$\begin{array}{}\text{(15)}& p\left({\mathit{\theta }}_j\right|\mathrm{\Delta }\mathit{T})=\int p(\mathit{\theta }\left|\mathrm{\Delta }\mathit{T}\right)\text{d}{\mathit{\theta }}_{-j},j=\mathrm{1},\mathrm{\dots },\mathrm{4}.\end{array}$$

Here, the notation θ_−j is used to denote the vector θ excluding the jth parameter. The posterior distributions provide a complete description of the latent field components and the parameters in our model. From the marginals in Eqs. (14)–(15) we can extract summary statistics such as the mean, variance, quantiles, and credible intervals.

Traditionally, marginal posterior distributions have been approximated using Markov chain Monte Carlo methods (Robert and Casella, 1999). Such methods are simulation-based and can potentially be very time-consuming for hierarchical models. The INLA methodology represents a computationally superior, but still accurate, alternative and is available using the R package R-INLA. This package can be downloaded for free at http://www.r-inla.org (last access: 10 December 2019). INLA provides a deterministic approach by approximating the posterior distributions in Eqs. (14)–(15) using numerical optimization techniques, interpolations, and numerical integration. Among other features, this includes the use of the Laplace approximation (Tierney and Kadane, 1986), which is an old technique to compute high-dimensional integrals. Specifically, the joint posterior distribution for the model parameters in Eq. (15) is approximated by employing a Laplace approximation evaluated at the mode x^∗(θ):

$$\begin{array}{}\text{(16)}& p\left(\mathit{\theta }\right|\mathrm{\Delta }\mathit{T})\approx {\left.{\displaystyle \frac{p(\mathit{x},\mathit{\theta },\mathrm{\Delta }\mathit{T})}{p_G\left(\mathit{x}\right|\mathit{\theta },\mathrm{\Delta }\mathit{T})}}\right|}_{\mathit{x}={\mathit{x}}^{\ast }\left(\mathit{\theta }\right)},\end{array}$$

where $p_G\left(\mathit{x}\right|\mathit{\theta },\mathrm{\Delta }\mathit{T})$ is a Gaussian approximation of

$$p\left(\mathit{x}\right|\mathit{\theta },\mathrm{\Delta }\mathit{T})\propto p(\mathit{x}\left|\mathit{\theta }\right)p\left(\mathrm{\Delta }\mathit{T}\right|\mathit{x},\mathit{\theta }).$$

This approximation is usually very accurate as we know that p(x|θ) is already Gaussian. The marginal for each model parameter is then obtained by assuming a normal distribution modified to allow for skewness,

$$p\left({\mathit{\theta }}_j\right|\mathrm{\Delta }\mathit{T})\approx \left\{\begin{array}{ll}\mathcal{N}\left(\mathrm{0},{\mathit{\sigma }}_{j+}^{\mathrm{2}}\right),& {\mathit{\theta }}_j>\mathrm{0}\\ \mathcal{N}\left(\mathrm{0},{\mathit{\sigma }}_{j-}^{\mathrm{2}}\right),& {\mathit{\theta }}_j\le \mathrm{0}\end{array}\right..$$

The scaling parameters σ_j+ and σ_j− are found using the approximate joint posterior distribution of Eq. (16); see Martins et al. (2013) for details. To compute Eq. (14), the Laplace approximation in Eq. (16) is combined with a simplified and computationally faster version of the Laplace approximation of $p\left(x_t\right|\mathit{\theta },\mathrm{\Delta }\mathit{T})$. Finally, the integrand of Eq. (14) is evaluated for values of θ in a grid efficiently covering the parameter space for θ; see Rue et al. (2009, 2017) for details.

A key assumption for the numerical approximations to be computationally efficient is that the latent Gaussian field x has Markov properties; i.e., x needs to be a Gaussian Markov random field having a sparse precision matrix Q (Rue and Held, 2005). This is not the case for fGn as the long-range dependency structure of this process gives a dense precision matrix. We resolve this problem by approximating ε as a weighted sum of m independent first-order autoregressive (AR(1)) processes, i.e.,

$$\begin{array}{}\text{(17)}& \stackrel{\mathrm{̃}}{\mathit{\epsilon }}=\sum _{i=\mathrm{1}}^m\sqrt{w_i}{\stackrel{\mathrm{̃}}{\mathit{x}}}_i.\end{array}$$

To capture the correlation structure between $\stackrel{\mathrm{̃}}{\mathit{\epsilon }}$ and each of the AR(1) processes ${\stackrel{\mathrm{̃}}{\mathit{x}}}_i$, the latent field must be extended to also include the underlying AR(1) processes, i.e., $\mathit{x}=(\mathit{\eta },\mathit{\mu }+\stackrel{\mathrm{̃}}{\mathit{\epsilon }},{\stackrel{\mathrm{̃}}{\mathit{x}}}_{\mathrm{1}},\mathrm{\dots },{\stackrel{\mathrm{̃}}{\mathit{x}}}_m)$. The weights $\mathit{\{}{w}_i{\mathit{\}}}_{i=\mathrm{1}}^m$ and the first-lag autocorrelation coefficients of the AR(1) processes are selected such that the resulting autocorrelation function of $\stackrel{\mathrm{̃}}{\mathit{\epsilon }}$ best approximates that of fGn. In addition to ensuring computational efficiency, this approximation also proves to be remarkably accurate. For further details about this approximation, see Sørbye et al. (2019), who also provide a discussion from a statistical perspective. For a physical interpretation of this approximation we refer to Fredriksen and Rypdal (2017).

Currently, there are no built-in model components in R-INLA that suit our specifications. This means that we have to construct one manually using rgeneric, a modeling tool that allows generic model components to be defined for INLA. To make this accessible to applied scientists we have developed an R package called INLA.climate, which includes functions that take care of the technical part of the fitting procedure and presents important information and summary statistics in a readable format. This package contains a versatile and user-friendly interface to fit the model in Eq. (9) and includes functions to replicate all results presented in this paper. The package is available at the GitHub repository at https://github.com/eirikmn/INLA.climate (last access: 11 February 2020). A detailed description of the package and its features is available in its accompanying documentation.

3 Applications

3.1 Estimating the forced temperature response

As explained in Sect. 2, our model formulation implies that the sum μ+ε is viewed as one model component. Consequently, INLA will give an estimate for the posterior distribution of the sum, and not the marginal posterior distributions for each of the terms μ and ε.

Figure 1 The marginal posterior distributions of the parameters obtained using INLA.climate to fit our model to the GISS-E2-R temperature and forcing dataset. The vertical lines show the mean and 95 % credible intervals.

In this example, we illustrate how we can approximate the marginal posterior distribution for the temperature response attributed to the known forcing, p(μ_i|ΔT) by combining INLA with Monte Carlo sampling. We first fit our model to the GISS-E2-R temperature and the corresponding forcing data using INLA. This gives the estimated marginal posterior distributions for each of the model parameters $\mathit{\theta }=(H,{\mathit{\sigma }}_{\mathit{\epsilon }},{\mathit{\sigma }}_{\mathrm{f}},F_{\mathrm{0}})$, as shown in Fig. 1.

Next, we use the inla.hyperpar.sample function from the R-INLA package to draw 100 000 samples from the approximate joint posterior distribution of H,σ_f, and F₀. For each of these samples, we compute μ according to Eq. (13). The resulting samples give approximate marginal posterior distributions for each μ_i, which can then be used to calculate summary statistics.

Figure 2 The marginal posterior mean and 95 % credible intervals of the temperature response to known forcing μ obtained by 100 000 Monte Carlo simulations compared to GISS-E2-R temperature data. Panel (a) shows the results under the scale-invariant assumption, panel (b) shows the results using a single-exponential response, and panel (c) shows the results using a mixture of two-exponential response functions.

For comparison, we apply the same approach to estimate the given marginal posterior distributions under the assumption of an exponential response function in Eq. (3). In this case, the discretized unforced response described in Eq. (4) is an AR(1) process. More generally, the discretized response functions corresponding to Eq. (6) are a weighted sum of m AR(1) processes. Notice that a mixture of a few AR(1) processes will in general have short-range dependence properties, while the approximation in Eq. (17) is constructed to exhibit the long-range dependency structure of fGn. Using the scale-invariant response function or the exponential response functions corresponding to the one- and two-box models, we can compute the marginal posterior means and 95 % credible intervals for each μ_i. The results are shown in Fig. 2. The marginal posterior means are very similar. However, we observe significantly wider credible intervals for the model using the single-exponential response function. The larger uncertainty suggests that a smaller portion of the variance is explained by the unforced climate variability, leaving more of the variation to be explained by the response to the known forcing. Using the INLA.climate package, we obtain full inference in seconds on a personal computer. The code to run the example is as follows.

3.2 Temperature predictions for Representative Concentration Pathway trajectories

Once trained on historical temperature and forcing data, the response model can easily be used to obtain temperature predictions for different future forcing scenarios. Here, we present global temperature predictions for the years 2016 to 2100 based on the HadCRUT4 temperature data and the greenhouse gas component of the Hansen forcing data for 1850 to 2015. For future forcing, we use the Representative Concentration Pathways (RCPs) RCP2.6 (van Vuuren et al., 2007), RCP4.5 (Clarke et al., 2007; Smith and Wigley, 2006; Wise et al., 2009), RCP6 (Fujino et al., 2006; Hijioka et al., 2008), and RCP8.5 (Riahi et al., 2007). These trajectories were first published in 2000 and cover the years 2000 to 2100. In our analyses, we use the RCPs for the year 2016 to the year 2100 and adjust each of them so that the forcing in 2015 equals the greenhouse gas forcing in the Hansen data in 2015. The forcing scenarios are shown in Fig. 3.

Figure 3 The greenhouse gas component of the Hansen forcing (black) followed by the RCP2.6 (blue), RCP4.5 (orange), RCP6 (red), and RCP8.5 (dark red) forcing scenarios.

Prediction is carried out using INLA.climate by appending the future scenario to the forcing of the past $\mathit{F}=({\mathit{F}}_{\text{past}},{\mathit{F}}_{\text{future}})$. The package automatically replaces missing observations with NA values and give predictions for these based on the model fitted to the observed data.

Figure 4 Panels (a)–(d) describe the marginal posterior means and 95 % credible intervals of the predicted temperature response to future forcing according to the RCP2.6, RCP4.5, RCP6, and RCP8.5 trajectories, respectively, using a scaling response function. These are compared to the AR5 projections (black boxes).

As in the previous example, we compare the results using the scale-invariant response function versus the exponential response functions corresponding to the one- and two-box models. Training and predictions only take seconds to carry out on a personal computer. Figure 4 shows the marginal posterior means and 95 % credible intervals using the scale-invariant response, while Figs. 5–6 show the corresponding results using the exponential response functions. The figures also show comparisons with the AR5 projections listed in table SPM.2 in IPCC (2013b). We observe that both of the exponential response models seem to fail in describing the persistence in the temperature response and underestimate the global warming increase projections. The predictions obtained using the scale-invariant response function give higher future temperatures estimates, slightly overestimating the AR5 projections.

Figure 5 Panels (a)–(d) describe the marginal posterior means and 95 % credible intervals of the predicted temperature response to future forcing according to the RCP2.6, RCP4.5, RCP6, and RCP8.5 trajectories, respectively, using a single-exponential response function. These are compared to the AR5 projections (black boxes).

Figure 6 Panels (a)–(d) describe the marginal posterior means and 95 % credible intervals of the predicted temperature response to future forcing according to the RCP2.6, RCP4.5, RCP6, and RCP8.5 trajectories, respectively, using a mixture of two-exponential response functions. These are compared to the AR5 projections (black boxes).

Figure 7 Panels (a)–(d) describe the deviations between the predicted marginal posterior means and the predictions obtained from each ESM.

In Fig. 7 we compare the scale-invariant model's temperature projections for RCP scenarios to ESM temperature projections. In this analysis, we tune the statistical model to historical runs of different ESMs in the CMIP5 ensemble. As in the other analyses, we use model-specific forcing for each of the ESMs. Using the same method as for historical forcing, we also estimate model-specific RCP forcing for each of the ESMs and each RCP scenario. From this RCP forcing, and using the estimated parameters, we project GMST under the RCP scenarios and compare with the corresponding projections in the ESMs. Figure 7 shows the differences between the two types of projections. As seen from the figure, there is a negative trend for almost all models and all scenarios, indicating that the predictions made using the statistical model slightly overestimate the temperature increase in the ESMs. This overestimation is not a statistical bias. Instead, it shows that scale invariance is too crude an approximation for several ESMs. One can interpret the differences as a weakness of the scale-free model. However, not all climate models have scaling properties consistent with temperature reconstructions, and it could also reflect the weaknesses of the ESMs.

3.3 Estimating the transient climate response

As a final application, we describe how our suggested method using a scale-invariant response function can be used to estimate the TCR. The TCR is defined as the average temperature response between 60 and 80 years following a gradual CO₂ doubling, assuming a 1 % annual increase. In this scenario, the forcing increases linearly according to

$$f\left(s\right)={\displaystyle \frac{Q_{\mathrm{2}\times {\mathrm{CO}}_{\mathrm{2}}}}{\mathrm{70}\text{ years}}}(s+F_{\mathrm{0}}),\text{for }s=\mathrm{1},\mathrm{\dots },\mathrm{80}\text{ years}.$$

Here, $Q_{\mathrm{2}\times {\mathrm{CO}}_{\mathrm{2}}}$ is a model-specific coefficient describing the forcing corresponding to a CO₂ doubling. We obtain these coefficients from Forster et al. (2013) for all ESMs analyzed in this paper.

The computation of TCR is carried out by inserting the forcing into Eq. (1) and performing the matrix multiplication

$$\mathit{v}={\mathit{\sigma }}_{\mathrm{f}}\mathbf{G}\left(H\right)\mathit{f},$$

where G(H) is the 80×80 matrix with elements defined as in Eq. (10) and $\mathit{f}={\left(f\left(\mathrm{1}\right),\mathrm{\dots },f\left(\mathrm{80}\right)\right)}^{\top }$. This implies that TCR is computed by

$$\begin{array}{}\text{(18)}& \text{TCR}={\displaystyle \frac{\mathrm{1}}{\mathrm{20}\text{ years}}}\sum _{t=\mathrm{60}\text{ years}}^{\mathrm{80}\text{ years}}{v}_t.\end{array}$$

As in Sect. 3.1, the approximate marginal posterior distribution for TCR is obtained by combining INLA with Monte Carlo sampling. We first generate samples from the joint posterior distribution of the model parameters p(θ|ΔT). For each of these samples, we calculate TCR, which then gives the approximate posterior distribution for TCR.

Table 1 The average absolute value of the bias, the root mean square error, and the correlation obtained when comparing the marginal posterior mean estimates of the TCR with the TCR obtained directly from the ESMs.

Hansen et al. (2010)Schmidt et al. (2014)Collins et al. (2011)The HadGEM2 Development Team (2011)Dufresne et al. (2013)Bentsen et al. (2013)Iversen et al. (2013)Bi et al. (2012)Watanabe et al. (2011)Watanabe et al. (2010)Chylek et al. (2011)Gent et al. (2011)Voldoire et al. (2013)Donner et al. (2011)Dunne et al. (2012)Dunne et al. (2013)Rotstayn et al. (2012)Jeffrey et al. (2013)Wu et al. (2014)Wu et al. (2014)Dunne et al. (2012, 2013)Volodin et al. (2010)Giorgetta et al. (2013)Yukimoto et al. (2012)

Table 2 The Earth system models used in this paper. The table includes ID number, references, time interval, TCR obtained directly from the ESM, and slope coefficient for the forcing corresponding to a CO₂ doubling.

Table 3 This table contains marginal posterior means and 95 % credible intervals for the model parameters and the transient climate response obtained from fitting the scale-invariant response model to temperature data from the first nine ESMs.

Table 4 This table contains marginal posterior means and 95 % credible intervals for the model parameters and the transient climate response obtained from fitting the scale-invariant response model to temperature data from the last 10 ESMs.

For our analyses, we use temperature datasets generated from 19 ESMs in the Coupled Model Intercomparison Project Phase 5 (CMIP5) ensemble; see Table 2. We obtain the forcing by combining the forcing data from Forster et al. (2013) and Hansen et al. (2010) such that the 18-year moving averages of the two are equal. We use the instrumental HadCRUT dataset (Morice et al., 2012), which combines the land temperatures of the CRU dataset (Jones et al., 2012) with the sea surface temperatures of HadSST3 (Kennedy et al., 2011).

To assess the accuracy of the TCR estimations from Eq. (12) we compare the estimates from each of the 19 ESMs with the TCR obtained from the ESMs directly (Forster et al., 2013). Inference is obtained by producing 100 000 Monte Carlo simulations of the TCR. Summary statistics for our model are shown in Tables 3–4, which include the marginal posterior means and 95 % credible intervals for the TCR and the model parameters used to compute it.

Figure 8 Scatter plot of the TCR obtained directly from the 19 ESMs against the corresponding marginal posterior mean estimates when using a scale-invariant response function (a), a single-exponential response function (b), and a two-exponential response function (c).

To assess the approach using the scale-invariant versus the exponential response functions, we compare the posterior mean estimates with the values obtained directly from the ESMs. Specifically, we calculate the average absolute value of the bias, the root mean square error (RMSE), and the correlation between the posterior mean estimates and the TCR values from the ESMs; see Table 1. We observe that the method using the single-exponential response model clearly performs worse than the other two approaches, having higher error and lower correlation with the TCR estimates from the ESMs. The two-exponential response function and the scale-invariant approach perform similarly, with the latter approach showing slightly higher correlation. These two methods have approximately the same average absolute value of the bias, but the two-exponential approach slightly underestimates TCR, while the scale-invariant approach slightly overestimates TCR. This is illustrated by the scatter plots in Fig. 8. Using INLA.climate, we obtained, for a typical analysis, inference in around 13 s using a scale-invariant response and 35 s using the exponential response functions.

Figure 9 The posterior mean estimates and 95 % credible intervals of the historical TCR for all forcing datasets using a scale-invariant (blue), a single-exponential (red), and a two-exponential (green) response function. ID number 0 denotes the Hansen forcing.

To obtain estimates for the TCR of the HadCRUT4 dataset we use 19 different forcing data points associated with the ESMs listed in Table 2 and the Hansen radiative forcing, which we will assign ID number 0. The Monte Carlo simulations are carried out separately for each forcing dataset using 100 000 samples and a forcing slope coefficient $Q_{\mathrm{2}\times {\mathrm{CO}}_{\mathrm{2}}}=\mathrm{3.8}$ W m⁻² (IPCC, 2013a). This is again performed using both the scale-invariant and the exponential response functions. For the scale-invariant response, the posterior means and credibility intervals of the TCR and the parameters used to compute the TCR for each ESM are shown in Tables 5–6. The marginal posterior mean estimates and 95 % credible intervals for the TCR using all approaches are illustrated in Fig. 9 where we observe wider credible intervals when using a single-exponential response function. We obtain an estimated posterior distribution for the TCR across all models by aggregating all TCR samples obtained from each analysis, totaling 2 million simulations. The posterior density is obtained from the Monte Carlo samples using the density function in R. The resulting density function is depicted in Fig. 10, where it is compared with a histogram describing the TCRs obtained directly from the ESMs. We observe a mean of 1.43, 1.35, and 1.61 K and standard deviation of 0.19, 0.40, and 0.40 K when using a scale-invariant, a single-exponential, and a two-exponential response function, respectively. The given estimates mainly fall in the lower half of the range of TCRs and fail to capture the mode of the histogram. However, since the histogram is generated from only 19 different values of TCR, its form is influenced by the size of the bins. Using bins of width 0.5 K (resulting in three total bins) would describe a more unimodal distribution with a mode in the 1.5–2.0 interval. The posterior distribution obtained from using either a scale-invariant or the exponential response functions is still on the lower side of this.

Table 5 This table contains marginal posterior means and 95 % credible intervals for the model parameters and the transient climate response obtained from fitting the scale-invariant response model to the HadCRUT dataset using forcing data from Hansen et al. (2010) (denoted by ID 0) and from the first nine ESMs.

Table 6 This table contains marginal posterior means and 95 % credible intervals for the model parameters and the transient climate response obtained from fitting the scale-invariant response model to the HadCRUT dataset using forcing data from the last 10 ESMs.

Figure 10 Histogram of the TCR obtained from the different ESMs with the posterior density function estimated from the Monte Carlo simulations using a scale-invariant (solid), a single-exponential (dashed), and a two-exponential (dotted) response function. The vertical lines describe the means of the three approaches.

4 Conclusions

This paper presents a Bayesian formulation to analyze a linear temperature response model to radiative forcing, incorporating long-range temporal dependence using a scale-invariant response function. Computational efficiency is achieved by incorporating the model within the R-INLA framework and adopting the approximation introduced in Sørbye et al. (2019). The benefits of this methodology are threefold. First, the model is both accessible and adaptable to more advanced models that require more trends and effects. Second, the approximations ensure low costs in both computational complexity and memory, even for long time series. Third, the method yields full Bayesian inference, giving a full description of the behavior of the time variables and model parameters.

In addition to providing parameter estimates, the model has been used to produce temperature predictions as responses to the four RCP forcing trajectories used to describe future radiative forcing. For comparison, we have also included prediction results using the one-box and two-box energy balance models, giving single-exponential or a mixture of two-exponential response functions, respectively. We observe that the exponential response models underestimate the predicted temperature compared to the projections made by the IPCC. On the other hand, the scale-invariant response models tend to overestimate future temperatures but are overall more accurate than using the exponential response functions.

We further demonstrate that the model can be used to estimate the transient climate response in instrumental data. Our best estimate is that the TCR is 1.43 K with a standard deviation of 0.29 K. This estimate falls right in the middle of the range put forward in the IPCC report (0.8–2.5 K), and the accuracy is consistent with the TCR obtained directly from the ESMs. The presented model has also given coherent estimates for the equilibrium climate sensitivity compared with running ESMs (Rypdal et al., 2018a).

Accurate linear response models for global temperature are essential alternatives to ESMs in studies in which one needs to explore a large number of emission scenarios, and the modeling framework presented here can easily be included in integrated assessment models. Moreover, since the models are invertible, they can efficiently compute forcing scenarios corresponding to given future scenarios for global temperatures. Hence, in combination with linear models for the CO₂ response to emissions, they can be used to obtain observation-based estimates of the remaining carbon budget in scenarios in which we reach the goals of the Paris Agreement.

In combination with dedicated ESM experiments, the methods presented in this paper can also be used to estimate global and regional climate sensitivity as a function of background state and timescale. One can use such estimates to study the effect of nonlinear responses across timescales and to obtain insight into how sensitivities and fluctuations change in the vicinity of climate tipping points.