Long-range persistence in global surface temperatures explained by linear multi-box energy balance models

TL;DR

This paper shows that long-range persistence in global surface temperatures can be explained by simple linear multi-box energy balance models, providing a physical basis for observed scale invariance without nonlinear effects.

Contribution

It demonstrates that scale-invariant temperature fluctuations can be modeled with linear multi-box energy balance models, explaining long-range persistence from months to centuries.

Findings

Linear models reproduce observed temperature scaling.

Two-box models explain land and sea temperature differences.

Parameters fitted from data match observed scale invariance.

Abstract

The temporal fluctuations in global mean surface temperature is an example of a geophysical quantity which can be described using the notions of long-range persistence and scale invariance/scaling, but this description has suffered from lack of a generally accepted physical explanation. Processes with these statistical signatures can arise from non-linear effects, for instance through cascade-like energy transfer in turbulent fluids, but they can also be produced by linear models with scale-invariant impulse-response functions. This paper demonstrates that on time scales from months to centuries, the scale-invariant impulse-response function of global surface temperature can be explained from simple linear multi-box energy balance models. This explanation describes both the scale invariance of the internal variability and the lack of a characteristic time scale of the response to external forcings. With parameters estimated from observational data, the climate response is approximately scaling in these models, even if the response function is not chosen to be scaling a priori. It is also demonstrated that the differences in scaling exponents for temperatures over land and for sea-surface temperatures can be reproduced by a version of the multi-box energy balance model with two distinct surface boxes.

Tables

Table 1. Table 1: Parameters estimated from data. The time scales τ k subscript 𝜏 𝑘 \tau_{k} and their weights b k , k = 1 , 2 , 3 formulae-sequence subscript 𝑏 𝑘 𝑘 1 2 3 b_{k},\;k=1,2,3 uniquely determine the estimated response functions. The parameter σ monthly subscript 𝜎 monthly \sigma_{\text{monthly}} is the estimated standard deviation of the monthly resolved stochastic forcing. T 0 , instr subscript 𝑇 0 instr T_{0,\text{instr}} and T 0 , moberg subscript 𝑇 0 moberg T_{0,\text{moberg}} are the initial temperature anomalies estimated for global instrumental temperature and the Moberg reconstruction. We note that there is some remnant seasonal variability that we were not able to remove by subtracting a mean seasonal cycle. This is not so apparent in a double-logarithmic plot of the PSD, but it causes our estimates of σ monthly subscript 𝜎 monthly \sigma_{\text{monthly}} to be slightly inconsistent with the actual residual variability. We also point out that our estimates of σ monthly subscript 𝜎 monthly \sigma_{\text{monthly}} are only valid for monthly resolved temperatures, or monthly temperatures sampled at a different time resolution.

	Unit
$({\tau }_1,{\tau }_2,{\tau }_3)$	yr	(0.5, 5, 50)	(1, 10, 100)	(1, 20, 400)
$b_1$	${\text{K m}}^2{\text{ W}}^{-1}{\text{ yr}}^{-1}$	0.198	0.165	0.187
$b_2$	${\text{K m}}^2{\text{ W}}^{-1}{\text{ yr}}^{-1}$	0.033	0.022	0.018
$b_3$	${\text{K m}}^2{\text{ W}}^{-1}{\text{ yr}}^{-1}$	0.011	0.005	0.001
${\sigma }_{\text{monthly}}$	${\text{W m}}^{-2}{\text{yr}}^{-1/2}$	0.73	0.70	0.68
$T_{0,\text{instr}}$	K	-0.13	-0.14	-0.14
$T_{0,\text{moberg}}$	K	-0.18	-0.17	-0.15
$S_{\text{eq}}$	${\text{K m}}^2{\text{ W}}^{-1}$	0.79	0.85	1.07
$C_1$	${\text{W yr m}}^{-2}{\text{ K}}^{-1}$	4.15	5.09	4.86
$C_2$	${\text{W yr m}}^{-2}{\text{ K}}^{-1}$	15.4	22.0	31.7
$C_3$	${\text{W yr m}}^{-2}{\text{ K}}^{-1}$	24.1	41.8	149
${\kappa }_1$	${\text{W m}}^{-2}{\text{ K}}^{-1}$	1.26	1.18	0.94
${\kappa }_2$	${\text{W m}}^{-2}{\text{ K}}^{-1}$	5.61	3.30	3.50
${\kappa }_3$	${\text{W m}}^{-2}{\text{ K}}^{-1}$	1.74	1.2	0.91

Table 2. Table 2: Parameters estimated from global sea and land surface temperatures. The fixed time scales are chosen to be τ 1 = 0.5 subscript 𝜏 1 0.5 \tau_{1}=0.5 yrs, τ 2 = 5 subscript 𝜏 2 5 \tau_{2}=5 yrs and τ 3 = 50 subscript 𝜏 3 50 \tau_{3}=50 yrs.

sea:			land:
$b_1$	0.092	${\text{K m}}^2{\text{ W}}^{-1}{\text{ yr}}^{-1}$	$r_1$	0.10	${\text{K m}}^2{\text{W}}^{-1}$
$b_2$	0.035	${\text{K m}}^2{\text{ W}}^{-1}{\text{ yr}}^{-1}$	$r_2$	1.40
$b_3$	0.009	${\text{K m}}^2{\text{ W}}^{-1}{\text{ yr}}^{-1}$
${\sigma }_{\text{monthly,sea}}$	0.76	${\text{W m}}^{-2}{\text{yr}}^{-1/2}$	${\sigma }_{\text{monthly,land}}$	0.2	${\text{W m}}^{-2}{\text{yr}}^{-1/2}$
$T_{0,\text{sea}}$	-0.13	K	$T_{0,\text{land}}$	-0.15	K
$S_{\text{eq,sea}}$	0.69	${\text{K m}}^2{\text{ W}}^{-1}$	$S_{\text{eq,land}}$	1.07	${\text{K m}}^2{\text{ W}}^{-1}$

Equations

$$C\frac{d\Delta T}{dt}=-\frac{1}{S_{\text{eq}}}\Delta T+\Delta F(t).$$

$$S_{\text{eq}}=\frac{\Delta T}{\Delta F}.$$

$$\Delta T(t)={\int }_{-\infty }^{t}R(t-s)\Delta F(s)ds,$$

$$R(t)=\frac{1}{C}{e}^{-t/\tau }.$$

$$C_1\frac{d\Delta T_1}{dt}$$

$$C_2\frac{d\Delta T_2}{dt}$$

$$\mathbf{C}\frac{d\Delta \mathbf{T}}{dt}=\mathbf{K}\Delta \mathbf{T}+\Delta \mathbf{F}(\mathbf{t})$$

$$\mathbf{C}=\left(\begin{array}{cc}{C}_1& 0\\ 0& C_2\end{array}\right),\Delta \mathbf{T}=\left(\begin{array}{c}\Delta T_1\\ \Delta T_2\end{array}\right),\Delta \mathbf{F}(\mathbf{t})=\left(\begin{array}{c}\Delta F(t)\\ 0\end{array}\right)$$

$$\mathbf{K}=\left(\begin{array}{cc}-({\kappa }_1+{\kappa }_2)& {\kappa }_2\\ {\kappa }_2& -{\kappa }_2\end{array}\right).$$

$$\Delta \mathbf{T}(t)={\int }_{-\infty }^t{e}^{(t-s)\mathbf{A}}\Delta \mathbf{F}(s)ds,\text{ with }\mathbf{A}={\mathbf{C}}^{-1}\mathbf{K},$$

$$R(t)={\left(e^{t\mathbf{A}}\right)}_{11}=\sum _{k=1}^N{b}_k{e}^{-t/{\tau }_k}.$$

$$C_{\text{L}}\frac{d\Delta T_{\text{L}}}{dt}$$

$$C_1\frac{d\Delta T_1}{dt}$$

$$\Delta T_{\text{L}}=\frac{F_{\text{L}}(t)+k\mu \alpha \Delta T_1/f_{\text{L}}}{{\lambda }_{\text{L}}+k/f_{\text{L}}}.$$

$$\Delta T_{\text{global}}=f_{\text{L}}\Delta T_{\text{L}}+f_{\text{O}}\Delta T_1=\left(f_{\text{O}}+\frac{k\mu \alpha }{{\lambda }_{\text{L}}+k/f_{\text{L}}}\right)\Delta T_1+\frac{f_{\text{L}}}{{\lambda }_{\text{L}}+k/f_{\text{L}}}{F}_{\text{L}}(t).$$

$$dx(t)=-\theta x(t)dt+\sigma dB(t),$$

$$x(t)={\int }_{-\infty }^{t}R(t-s)dB(s)\text{ with }R(t)=\sigma e^{-\theta t}.$$

$$⟨x(t)x(t+\tau )⟩={\int }_0^{\infty }R(t^{\prime })R(t^{\prime }+\tau )d{t}^{\prime }=\frac{{\sigma }^2}{2\theta }{e}^{-\theta t},$$

$$\Delta T_{1}(t)=\int_{-\infty}^{t}\Big{(}\sum_{k=1}^{N}b_{k}e^{\lambda_{k}(t-s)}\Big{)}d\Delta F(s)=\sum_{k=1}^{N}b_{k}\Delta T_{1,k}(t),$$

$$\Delta T_{1,k}(t)={\int }_{-\infty }^t{e}^{{\lambda }_k(t-s)}d\Delta F(s).$$

$$d\Delta F(t)=\Delta F_{\text{det}}(t)dt+\sigma dB(t).$$

$$\Delta T_{1}(t)=\Delta T_{1,\text{det}}(t)+\sigma\int_{-\infty}^{t}\Big{(}\sum_{k=1}^{N}b_{k}e^{\lambda_{k}(t-s)}\Big{)}dB(s)=\Delta T_{1,\text{det}}(t)+\sigma\sum_{k=1}^{N}b_{k}x_{k}(t),$$

$$x_k(t)={\int }_{-\infty }^t{e}^{{\lambda }_k(t-s)}dB(s)$$

$$\Delta T_{1,k}(\omega )=\frac{\Delta F(\omega )}{i\omega +1/{\tau }_k},$$

$$S_1(\omega )=\underset{T\to \infty }{\lim }\frac{1}{T}⟨|\Delta T_1(\omega ){|}^2⟩=\underset{T\to \infty }{\lim }\frac{1}{T}⟨|\Delta F(\omega ){|}^2⟩\left[S^{(0)}(\omega )+S^{(\text{cr})}(\omega )\right],$$

$$S^{(0)}(\omega )=\sum _k\frac{b_k^2}{\left({\omega }_k^2+{\omega }^2\right)}\text{ , }{S}^{(\text{cr})}(\omega )=\sum _k\sum _{j<k}\frac{2b_k{b}_j\left({\omega }_k{\omega }_j+{\omega }^2\right)}{\left({\omega }_j^2+{\omega }^2\right)\left({\omega }_k^2+{\omega }^2\right)},\text{ and }{\omega }_k=1/{\tau }_k.$$

$$\mathbf{C}=\left(\begin{array}{ccc}{C}_1& 0& 0\\ 0& C_2& 0\\ 0& 0& C_3\end{array}\right),\Delta \mathbf{T}=\left(\begin{array}{c}\Delta T_1\\ \Delta T_2\\ \Delta T_3\end{array}\right),\Delta \mathbf{F}(\mathbf{t})=\left(\begin{array}{c}\Delta F(t)\\ 0\\ 0\end{array}\right)$$

$$\mathbf{K}=\left(\begin{array}{ccc}-({\kappa }_1+{\kappa }_2)& {\kappa }_2& 0\\ {\kappa }_2& -({\kappa }_2+{\kappa }_3)& {\kappa }_3\\ 0& {\kappa }_3& -{\kappa }_3\end{array}\right).$$

$$\Delta T_{1,k}(t)={\int }_{t_0}^t{e}^{-(t-s)/{\tau }_k}\Delta F(s)ds\text{ for }k=1,2\text{ and }3,$$

$$\Delta T_1(t)=b_1\Delta T_{1,1}(t)+b_2\Delta T_{1,2}(t)+b_3\Delta T_{1,3}(t),$$

$$C_{\text{1}}\frac{d}{dt}\Delta T_{\text{1}}=-{\kappa }_1^{\prime }\Delta T_{\text{1}}+F_{\text{N}}+F_{\text{O}}(t)+\frac{k}{f_{\text{O}}\left({\lambda }_{\text{L}}+k/f_{\text{L}}\right)}{F}_{\text{L}}(t)$$

Full text

1 Introduction

Instrumental measurements and proxy reconstructions of Earth’s surface temperatures show temporal variability on a range of different time scales (Lovejoy 2015; Huybers and Curry 2006). For the global mean surface temperature (GMST) the variability can be parsimoniously described as scale invariant, since the estimated power spectral densities (PSDs) are well approximated by power-laws $S(f)\sim 1/f^{\beta}$ from monthly to centennial scales (Rypdal et al. 2013). The typical scaling exponent is $\beta\approx 1$, and the signals are well described as a so-called $1/f$-noise, or pink noise. Some of the low-frequency variability in the temperature records can be accounted for by the variability in the radiative forcing of the planet, but even the residual fluctuations are well described as a scaling stochastic process, with a slightly lower exponent $\beta$. This suggests that scale-invariant dynamics is an intrinsic property of the climate system, a claim that is supported by the observation of scaling PSDs in unforced control runs of general circulation models (GCMs), on time scales from months to centuries (Fredriksen and Rypdal 2016; Rybski et al. 2008; Fraedrich and Blender 2003).

A signal with power-law PSD can be modelled as a stochastic process with long-range dependence (LRD), and examples of such processes are the fractional Gaussian noises (fGns) and the fractional autoregressive integrated moving average (FARIMA) models. Stochastic processes that exhibit LRD provide more accurate descriptions of the unforced GMST variability compared to the traditional red-noise models, such as the Ornstein-Uhlenbeck (OU) processes and the autoregressive processes of order 1 (AR(1)) (Rypdal and Rypdal 2014). The latter are characterised by a single time scale, and are incapable of describing the multi-scale nature of the climate fluctuations. Despite this, the LRD processes are largely ignored by many climate scientists, and some consider LRD to be an exotic and redundant notion in climate science (Mann 2011).

One of the aims of this paper is therefore to de-mystify the notion of LRD in the climate system by demonstrating that the observed phenomena can be produced by simple multi-box energy balance models (EBMs). With this, we demonstrate that the exotic physics may be no more than vertical heat conduction in the ocean, and that it is reasonable to think of LRD as an approximation to the linear response of EBMs with multiple characteristic time scales. Only a few boxes are needed to obtain power-law PSDs on scales from months to centuries. We also demonstrate how we can construct box models that are consistent with the observation that the exponent $\beta$ is lower for land temperatures than for sea surface temperatures (SSTs) (Fredriksen and Rypdal 2016).

Only a few of the studies that analyze LRD in surface temperatures, focus on the mechanisms behind the phenomenon (e.g. Fraedrich and Blender (2003), Fraedrich et al. (2004), Blender et al. (2006), Franzke et al. (2015)). Most treat LRD-processes merely as statistical models that fit well with data (Vyushin et al. 2012; Rybski et al. 2006; Franzke 2010). Statistical inference for LRD processes requires special care to avoid the fallacy of circular reasoning, that is, falsely attributing trends in the forced signal to natural variability (Benestad et al. 2016). Incautious trend-significance testing using LRD null models (Cohn and Lins 2005) have led some climate scientists to view LRD processes as exotic mathematical objects that somehow fit with the “climate denier agenda” (Mann 2011; Benestad et al. 2016). This is paradoxical, since climate response models that exhibit LRD actually display more “heating in the pipeline”, and compared with other response models they predict that emissions of greenhouse gases must be reduced earlier and more drastically to avoid dangerous anthropogenic influence (K. Rypdal 2016; Rypdal and Rypdal 2014).

Other climate scientists consider scaling to be an inherent property of atmospheric turbulent flows and certain types of regime switching dynamics (Lovejoy and Schertzer 2013; Franzke et al. 2015), and as such a signature of the nonlinearity of the underlying dynamics. In fact, Huybers and Curry (2006) hypothesize that the persistent scaling of surface temperatures observed on decadal to multicentennial scales is due to a nonlinear cascade driven by the seasonal forcing. They present a bicoherence spectrum in favor of this hypothesis, but the phase correlations that give rise to high bicoherence do not imply an effective nonlinear energy transfer between the seasonal and the multidecadal scales. We have also had problems in reproducing the bicoherence spectra reported in this paper. In a forthcoming paper we will examine this hypothesis in depth.

Since the ocean has a large heat capacity compared to the atmosphere, the observation that ocean temperatures are more persistent than land temperatures (Fraedrich and Blender 2003; Fredriksen and Rypdal 2016), is an indication that the observed persistence in global temperature to a larger extent must be attributed to ocean heat content and ocean dynamics, and to a lesser extent to nonlinear processes in the atmosphere. This hypothesis is further strengthened by the results of Fraedrich and Blender (2003), who find that only models with full ocean circulation show persistence on scales longer than about a decade. In the present paper we model the slowly responding components of the climate system by including “boxes” that exchange heat with the more rapidly responding mixed layer. This is clearly an oversimplification of the ocean dynamics, but reproduces the multi-scale characteristics of the surface temperature response.

The paper is structured as follows. Sec. 2 discusses the construction of multi-box EBMs and their corresponding response functions, and in Sec. 3 we demonstrate how the superposition of different response times can be used to approximate an LRD response. Furthermore, we estimate parameters and explore how the response of sea surface temperatures differ from the response of land temperatures. Sec. 4 presents some concluding remarks.

2 Multibox EBMs

The simplest climate model we can imagine is the so-called one-box EBM for the global temperature,

[TABLE]

In this equation, $C$ denotes the average heat capacity per square meter of the surface, $\Delta T$ the temperature anomaly relative to an equilibrium state, $S_{\text{eq}}$ the equilibrium climate sensitivity, and $\Delta F(t)$ is the forcing, i.e., the perturbation of effective radiative forcing from the initial equilibrium state $\Delta T=0$. As a response to a constant perturbation $\Delta F$ the temperature will reach a new equilibrium $\Delta T$, and the change in equilibrium temperature relative to the change in radiative forcing is equal to the equilibrium climate sensitivity, i.e.,

[TABLE]

For a time-dependent forcing $\Delta F(t)$, the temperature $\Delta T(t)$ is given by a convolution integral

[TABLE]

where the impulse response function is an exponentially decaying function with a characteristic time scale $\tau=CS_{\text{eq}}$:

[TABLE]

In the one-box model there is no heat exchange with the deep ocean, but this can be included by extending the model to also include a box with a larger heat capacity $C_{2}$. If the energy exchange between the upper and lower box is proportional to the temperature difference between the two boxes, we obtain what is known as the two-box EBM (Geoffroy et al. 2013; Held et al. 2010; Rypdal 2012; Caldeira and Myhrvold 2013):

[TABLE]

The equations can be written on matrix-form:

[TABLE]

where we introduce the notation

[TABLE]

and

[TABLE]

For convenience we denote $\kappa_{1}=1/S_{\text{eq}}$, but it should be noted that the physical meaning of $\kappa_{1}$ is different from $\kappa_{2}$. While $\kappa_{2}$ is a coefficient of heat transfer between two ocean layers, $\kappa_{1}$ is determined by the response of the outgoing long-wave radiation (OLR) to changes in the surface temperature. It includes all the atmospheric feedbacks and is sometimes referred to as the equilibrium global climate feedback (Armour et al. 2013), or simply the feedback parameter.

The natural generalization of the two-box model is to consider $N$ vertically distributed boxes. The model is formulated as in Eq. (6), with $\Delta\mathbf{T}$ and $\Delta\mathbf{F}$ being $N$-vectors, and $\mathbf{C}$ and $\mathbf{K}$ being $N\times N$ matrices. $\mathbf{C}$ will be a diagonal matrix with the heat capacities of each box along the diagonal, and $\mathbf{K}$ will be a tridiagonal matrix. We note that when $N$ is large, this model set-up can approximate a vertical diffusion model.

In an $N$-box EBM the temperature $N$-vector can be written using matrix-exponential notation:

[TABLE]

and it follows that the surface temperature is given by a convolution integral similar to the one in Eq. (2), but where the impulse response function is now a weighted sum of $N$ exponentially decaying functions:

[TABLE]

The characteristic time scales are defined as $\tau_{k}=-1/\lambda_{k},\;k=1,...,N$, where $\lambda_{k}$ are the eigenvalues of the matrix $\mathbf{C}^{-1}\mathbf{K}$. Since $-\mathbf{K}$ is symmetric and positive definite, the eigenvalues $\lambda_{k}$ are real and negative.

The model defined by Eqs. (6) is meant to describe vertically distributed boxes, but boxes can also be aligned horizontally. This can be useful in order to include the atmosphere over land in the model. In principle we can have interactions between all boxes, making the matrix $\mathbf{K}$ less sparse. The mathematical form of the response function remains the same though, but the characteristic time scales and the weights $b_{k}$ are changed. Several horizontally distributed boxes could also be useful for modeling a space-dependent depth of the mixed layer.

To make separate boxes for the upper ocean layer and atmosphere over land we adopt the asymmetric heat exchange between land and sea used by Meinshausen et al. (2011) to obtain the equations

[TABLE]

Here it is assumed that the temperature in the atmosphere over oceans $\Delta T_{\text{O,atmos}}$ is proportional to the temperature in the mixed layer, i.e. $\Delta T_{\text{O,atmos}}=\alpha\Delta T_{1}$, where the factor $\alpha>1$ describes the effect of changing sea ice cover (Raper et al. 2001). The parameter $\mu>1$ quantifies the asymmetry in the heat transport between the atmosphere over the ocean and the atmosphere over land. The parameter $f_{\text{L}}=0.29$ is the proportion of Earth’s surface that is covered by land, and $f_{\text{O}}=1-f_{L}$. $F_{\text{N}}$ represents the heat transport into the deep oceans, and $F_{\text{L}}$ and $F_{\text{O}}$ are the forcing terms over land and ocean respectively. From the Coupled Model Intercomparison Project Phase 3 (CMIP3) models one finds that the typical values of $\mu$ are in the range 1 to 1.4 (Meinshausen et al. 2011). This implies that, when a new equilibrium is reached after a perturbation of the forcing, the land temperature will have changed more than the SST.

In the limit $C_{\text{L}}\to 0$, Eq. (9) becomes

[TABLE]

Hence, land temperature appears as a weighted sum of the SST and an instantaneous response to the forcing over land. The GMST anomaly is given by

[TABLE]

3 Approximate scale invariance from aggregation of OU processes

The Ornstein-Uhlenbeck stochastic process is defined via the stochastic differential equation

[TABLE]

where $dB(t)$ is the white noise measure. The equation has a stationary solution on the form

[TABLE]

The parameter $\sigma$ is called the scale parameter, and $\theta$ is the damping rate. Since $dB(t)$ is a white noise it follows that

[TABLE]

hence, the characteristic correlation time of an OU process is $\tau=1/\theta$. In the multi-box EBM with $N$ vertically distributed boxes, the temperature $\Delta T_{1}(t)$ is given by

[TABLE]

where

[TABLE]

If we consider the perturbations of the radiative forcing caused by volcanoes, solar variability, and anthropogenic activity as “deterministic”, and the perturbations from the chaotic atmospheric dynamics as random, then it is natural to model the forcing as a superposition of a deterministic component and a white-noise random process;111In this paper we follow Rypdal and Rypdal (2014) and model the random component of the forcing as white noise. However, the models can easily be modified to other stochastic models for the random forcing.

[TABLE]

Since the $N$-box models we consider are linear, the decomposition of the forcing yields a straightforward decomposition of the temperature response;

[TABLE]

where the processes

[TABLE]

are dependent OU processes with characteristic time scales given by the eigenvalues of the matrix ${\bf C}^{-1}{\bf K}$ via the relations $\tau_{k}=-1/\lambda_{k}$. Taking the Fourier transform of Eq. (14) yields

[TABLE]

and the PSD of $\Delta T_{1}(t)$ becomes

[TABLE]

where

[TABLE]

Here $T$ is the length of the time series $\Delta T_{1}(t)$, $S^{(0)}(\omega)$ is the PSD of an independent superposition of the processes $\Delta T_{1,k}(t)$, and $S^{(\text{cr})}(\omega)$ is the contribution to the PSD from the cross terms, which can not be neglected since the processes $\Delta T_{1,k}(t)$ are driven by the same forcing $\Delta F(t)$. For the stochastic component of the process Eq. (15) we can replace the forcing by a white noise process, such that $\lim_{T\rightarrow\infty}\frac{1}{T}\langle|\Delta F(\omega)|^{2}\rangle$ is a constant in $\omega$.

The PSD in Eq. (16) can easily be made to approximate a power law. For instance we can pick time scales $\tau_{k}$ such that $\tau_{k+1}=a\tau_{k}$ and weights $b_{k}$ such that $b_{k+1}=\sqrt{a^{\beta-2}}b_{k}$. Then we have the approximate relations $S^{(0)}(\omega)\sim\omega^{\beta}$, $S^{(\text{cr})}(\omega)\sim\omega^{\beta}$ and hence $S_{1}(\omega)\sim\omega^{\beta}$. In Fig. 1 this is demonstrated for a superposition of $N=5$ terms. The sum $S^{(0)}(\omega)$ is shown by the red line in (a) and the total sum $S^{(0)}(\omega)+S^{(\text{cr})}(\omega)$ by the blue line in (b). The idea that a long memory process can be produced by aggregating OU processes is the same as presented by e.g. Granger (1980) and M. Rypdal (2016). In this paper the time scales $\tau_{k}$ and their weights $b_{k}$ are estimated from data without a priori assumptions of scale invariance.

3.1 Example 1. The two-box model

For the classical two-box model, Geoffroy et al. (2013) estimated parameters by fitting the GMST-response in the two-box model to the corresponding response in CMIP5 models. The forcing scenario used for fitting were abrupt quadrupling of atmospheric CO2 concentration. They find multi-model mean parameter estimates $\hat{C}_{1}=7.3$ W yr m*-2* K*-1*, $\hat{C}_{2}=106$ W yr m*-2* K*-1*, $\hat{\kappa}_{1}=1.13$ Wm*-2K-1*, $\hat{\kappa}_{2}=0.73$ Wm*-2K-1*, which correspond to characteristic time scales of $3.88$ years and $242$ years. The two-box model provides a good fit to CMIP5 abrupt 4$\times$CO2 experiments and 1$\%$ per year CO2 increase experiments over 140 years, but if forced with white noise, the PSD is $S_{1}(\omega)\propto S^{(0)}(\omega)+S^{(\text{cr})}(\omega)$. Using the parameters estimated by Geoffroy et al. (2013) for three different CMIP5 models, we have plotted this expression in Fig. 2(a). As seen from the figure, the PSD can be approximated by two different power laws; one in the high-frequency range and another in the low-frequency range. For frequencies corresponding to time scales greater than a few decades the PSD can be approximated by $S(f)\propto 1/f^{\beta}$ with $\beta=0.3$, and for the higher frequencies it can be approximated by $S(f)\propto 1/f^{2}$. This result is inconsistent with the PSDs estimated from CMIP5 control run temperatures, which are well approximated by one power-law over the entire range of frequencies from months to centuries (Fredriksen and Rypdal 2016).

The inability of the two-box model to simultaneously describe the the average GMST response to certain forcing scenarios in GCMs, as well as the PSD of the background fluctuations, serves as a motivation to analyze more general $N$-box models, and in the next example we consider the EBM with three vertically distributed boxes. We will demonstrate that the three-box model provides accurate descriptions of both the “deterministic” response to historic radiative forcing and the statistical properties of the response to random forcing.

3.2 Example 2. The three-box model

The three-box model is given by Eq. (6) with

[TABLE]

and

[TABLE]

To estimate the parameters in the model we will make use of the HadCRUT4 data set for the GMST since 1850 (Morice et al. 2012), and the global effective forcing data, both with annual resolution. The forcing data is an updated version of Hansen et al. (2011), available at http://www.columbia.edu/~mhs119/Forcings/. We also use the Moberg Northern Hemisphere reconstruction (Moberg et al. 2005), and the Crowley forcing data (Crowley 2000) for the years 1000-1979. We fix a set of three well separated time scales ($\tau_{1}$, $\tau_{2}$ and $\tau_{3}$) and compute the responses

[TABLE]

to the historical forcing data $\Delta F(t)$, where the integral is estimated by a sum. As in Eq. (13), the GMST response is a linear combination of the responses $\Delta T_{1,k}(t)$:

[TABLE]

and our approach is to estimate the parameters $b_{1}$, $b_{2}$ and $b_{3}$ from historical data of GMST and forcing. We will subsequently demonstrate that for the range of time scales we consider in this paper, the results are largely insensitive to the choice of time scales ($\tau_{1}$, $\tau_{2}$ and $\tau_{3}$), as long as these are sufficiently separated. We do not only require that the “deterministic” response to radiative forcing fits well with observations, but also that the PSD of the stochastic component of the response is consistent with the estimated PSD of the residual observational signal (the difference between observed temperature record and the model response to the “deterministic” forcing). Rypdal and Rypdal (2014) employ a maximum-likelihood method to estimate the variance of the temperature residual $\sigma_{T}^{2}$ and the spectral exponent $\beta$ of a long-memory model. This method is inadequate here, since we have more free parameters and the method favorizes a good fit for the smallest time scales.

We employ instead an iterative routine, which in addition to weighting all the time scales in the estimation, allows us to fit to a composite spectrum. We first compute the residual temperature by guessing the model parameters, computing the deterministic responses to the Crowley and Hansen forcing time series, and then subtracting the deterministic responses from the observed Moberg and HadCRUT4 temperature records. The PSD of the residual is estimated using the periodogram, and subsequently log-binned. The theoretical PSD of the response to white-noise forcing in three-box model is then fitted to the composite residual PSD derived from the Moberg and HadCRUT4 temperatures by minimizing the square distance between the theoretical and estimated residual PSDs. From this we obtain new estimates of the relative size of the weights $b_{k}$ at each time scale, and using these we perform a new regression analysis with HadCRUT4 to determine the initial temperature $T_{0}$ and the absolute strength of the response. The procedure is repeated with these new model parameters as a starting point. The parameter estimation converges rapidly; only a few iterations are needed to obtain our estimates.

In Fig. 3(a) and Fig. 3(c) we show the deterministic responses to historical forcing with estimated parameters. The three different colors represent different choices of the time scales $\tau_{1}$, $\tau_{2}$ and $\tau_{3}$, and the parameter estimates are presented in Table 1, where we also present the corresponding values of the parameters $\kappa_{k}$ and $C_{k}$, as well as the equilibrium climate sensitivity of the model. We note that the colored curves in Fig. 3(a) and Fig. 3(c) are almost indistinguishable, and they all closely follow the HadCRUT4 and Moberg records.

Fig. 4(a) shows the theoretical PSD of the stochastic component (the response to white noise forcing) in the model. The estimated parameters are shown in Table 1, and the choice of time scales are $\tau_{1}=1$ yr, $\tau_{2}=10$ yr and $\tau_{3}=100$ yr. The PSD of the model fits well with the PSD estimated from observational data in the time-scale range from months to centuries, and in this range it is close to a power-law with an exponent $\beta=0.65$.

In Fig. 4(b) we show the three-box model responses to a unit step-forcing scenario. The different colors correspond to different choices of the time scales $\tau_{1}$, $\tau_{2}$ and $\tau_{3}$, and the gray curves are the corresponding two-box model responses with parameters estimated by Geoffroy et al. (2013) by fitting to abrupt 4$\times$CO2 experiments in CMIP5 models. From the response to a unit-step forcing we can also derive that the equilibrium climate sensitivity for an $N$-box model is given by $S_{\text{eq}}=\sum_{k=1}^{N}b_{k}\tau_{k}$.

The difference we observe in the step-forcing responses is not a result of a difference between the two- and three-box models, but rather a difference in estimation strategy. In our estimation procedure, we use only observational and proxy data, and have hence chosen the parameters that best reproduce both the residual spectra and responses to historical forcing. The parameters estimated by this method are not the same as the parameters that best describe the 4$\times$CO2 runs. This discrepancy can have several explanations; perhaps is the random forcing not accurately modeled as a white noise, or the forcing in the 4$\times$CO2 experiments is too strong for a linear approximation to be valid. It is likely that the feedback parameter or parameters related to ocean mixing can change during the strong and abrupt climate change following a quadrupling in CO2 concentration. The large differences between the different CMIP5 step-responses also reflect the large uncertainty associated with these model runs.

3.3 Example 3. Separate boxes for land and ocean

With two surface boxes, one for land and one for ocean, the equation for the ocean surface temperature is

[TABLE]

where

[TABLE]

The response function can hence be estimated in the same way as for global temperature, and the results are given in Table 2. For the estimation we have used the global SST dataset HadSST3 (Kennedy et al. 2011a, b).

Simplifying the constants in Eq. (10), and separating the stochastic and deterministic parts results in

[TABLE]

where $\sigma_{\text{L}}dB_{\text{L}}(t)$ is the direct stochastic forcing of the land surface temperature. For the ocean response we use the previously estimated parameters, given in Table 2. The remaining parameters are chosen such that the deterministic response

[TABLE]

is similar to global land surface temperature (LST), and such that the PSD of the residual temperature obtained after subtracting the deterministic response is similar to the PSD expected for the stochastic part of Eq. (20):

[TABLE]

Here we assume that the instantaneous response in ocean temperature to changes in the direct forcing of the land temperature is small compared to the land temperature response to this forcing. For the global LST we use the CRUTEM4v dataset (Jones et al. 2012). The observed SST, LST and the response to deterministic forcing with the parameters listed in Table 2 are shown in Fig. 5. This figure also shows the estimated and theoretical PSD of the response to stochastic forcing with the same parameters. The global temperature response $\Delta T_{\text{G}}=f_{\text{L}}\Delta T_{\text{L}}+f_{\text{O}}\Delta T_{1}$ is similar to the three-box temperature response estimated directly from global temperature.

The theoretical PSD of $\Delta T_{1}(t)$ fits well with the estimated PSD of the SST residual, and both are well-approximated by a power law with $\beta_{\text{O}}\approx 1$. In the same way, the theoretical PSD of $\Delta T_{\text{L}}(t)$ fits well with the estimated PSD of the LST, and both can be approximated by a power law with $\beta_{\text{L}}\approx 0.5$. These results are similar to the estimated PSDs of the linearly de-trended global LST and global SST analyzed in Fredriksen and Rypdal (2016). With this model, the only reason global LST shows persistence is because of the influence by global SST, but the persistence is weaker for land than for sea due to the component responding instantly to forcing.

We note that in the model presented in this paper, the relation between the equilibrium climate sensitivities for land temperatures and ocean temperatures is

[TABLE]

If $r_{1}/r_{2}\ll S_{\text{eq,sea}}$, there will be a near constant ratio between land and ocean temperature change, consistent with the findings of Lambert et al. (2011).

4 Concluding remarks

Simple climate models can be divided in two classes: EBMs and tuned impulse-response (IR) linear statistical models (Good et al. 2011). The power-law response model proposed by Rypdal and Rypdal (2014) is an example of an IR model that reproduces observed temperature variability quite well, but in this mathematical idealization conservation of energy is lost. Some may also find it problematic that it is not derived directly from physics. In this paper we demonstrate that such a model is closely approximated by the response of a multi-box EBM. This shows that LRD models can be seen as a compact mathematical descriptions of the effect of simple and concrete physical mechanisms.

Using only three boxes, we can produce a spectrum of the internal temperature variability that is practically indistinguishable from a power-law spectrum over a frequency range spanning more than two orders of magnitude. The response closely approximates a power-law response, even if we do not assume this a priori. The results of this paper are consistent with the findings of Fraedrich et al. (2004) and Lemke (1977), who demonstrate $1/f$-noise characteristics in linear diffusion models.

As suggested by (Ragone et al. 2015), we belive that the global temperature response in the Holocene can be well approximated as linear. Extending the linear response to local and regional temperatures is more problematic, especially for the temperature responses in the Arctic, where strong non-linear effects such as the sea-ice albedo feedback are present. However, the multi-box models can easily be extended to include non-linear terms, for instance to describe the rapid sea ice loss in the Arctic. They can also be extended to include different types of tipping points, and this allows ut to study critical transitions in systems that exhibit LRD. The effect of LRD on the so-called early-warning indicators associated with critical transitions in multi-box EMBs is a topic that will be pursued in future work.

Acknowledgements.

The authors thank K. Rypdal and O. Løvsletten for useful discussions. The paper was supported by the Norwegian Research Council (KLIMAFORSK programme) under grant no. 229754.