Climate variability from annual to multi-decadal timescales in a two-layer stochastic energy balance model: analytic solutions and implications for general circulation models

Abstract A low-order stochastically forced two-layer global energy balance model (EBM) admitting an analytical solution is developed for studying natural inter-annual, decadal and multi-decadal climate variability, and ultimately to better understand forced climate change. The EBM comprises upper and lower oceanic layers with a diffusive coupling, a radiative damping term including feedbacks and stochastic atmospheric forcing. The EBM is used to analyse the influence of radiative forcing, feedbacks and climate system inertia on the global mean surface temperature variance (climate variability) and to understand why Coupled Model Intercomparison Project, Phase 5 (CMIP5) models exhibit such a wide range in the level of variability in globally averaged surface air temperature. We examine the influence of the model parameters on the climate variability on different timescales. To this end, we derive the Fokker–Planck equation for the EBM and then obtain the analytical expression that quantifies the sensitivity coefficients for all model parameters. For all timescales, the most influential factors are as follows: (1) the magnitude of the stochastic forcing, (2) the feedback mechanisms, (3) the upper layer depth, (4) the diffusion parameter and (5) the lower ocean depth. Results from the EBM imply that the range of stochastic forcing in the CMIP5 climate models is around twice as important as the strength of radiative feedback or upper layer depth in causing the model-to-model spread in the magnitude of globally averaged climate model variability.


Introduction
Global temperatures are characterised by large globalscale 'natural variability' on timescales from months to decades (Power et al., 2006;Hawkins and Sutton, 2009;Deser et al., 2012). The importance of variability has been particularly prominent in recent years during the debate of the role of decadal contributions to a purported 'hiatus' in global warming (Collins et al., 2013). Although the precise nature, cause and importance of the 'hiatus' are still the subject of research (e.g. Kosaka and Xie, 2013), what is clear is that natural variability on timescales from years to multiple decades can either substantially reduce or exacerbate the effects of global warming, and it plays a critical role in the detection and attribution of climate change.
It has long been known that climate change projections under a range of projected greenhouse gas increases vary by a factor of around three in coupled general circulation models (CGCMs) (e.g. Bony et al., 2006), and this range has not decreased over a long period of model development (Flato et al., 2013). Perhaps less appreciated is that the range of natural variability found in GCMs is equally stark. Colman and Power (2018) find the standard deviation of global mean temperature at decadal timescales varies by a factor of more than four across Coupled Model Intercomparison Project, Phase 5 (CMIP5) models, with a similar range apparent in an earlier generation (CMIP3). In contrast to climate change, however, much less effort has gone into understanding the reasons behind this large spread.
Part of the challenge is the complexity of processes underlying natural variability, including features such as El Niño-Southern Oscillation (ENSO) and the Interdecadal Pacific Oscillation (Power et al., 2006), and Southern Ocean processes (Van Westen and Dijkstra, 2017). However, recent findings suggest that these and other coupled dynamic processes, important as they may be for regional-scale variability (such as ocean basin temperatures or regional rainfall), may be less prominent at global scales (Liu, 2012). The review paper of Liu (2012) concludes that the 'significant advancement of the last two decades is the recognition of the (atmospheric) stochastic forcing as the dominant generation mechanism for almost all interdecadal variability' (our emphasis). Consistent with this, Middlemas and Clement (2016) found that a model with only mixed layer physics showed decadal variability only 5% weaker than a version of the same GCM with full ocean dynamics. These two factors together suggest, prima facie, that much can be learned exploring the role of simpler 'Hasselmann' (1976)-type responses to stochastic radiative forcing in understanding variability in models.
A major advantage of the use of simplified, for example, energy balance models (EBMs) is that the hugely complex 'multi-dimensionality' of CGCM response can potentially be reduced to a few simple parameters, allowing straightforward exploration and understanding of model sensitivities. Under climate change scenarios, simple EBMs have been extensively used and it has been found that they can describe much of the multi-decadal temperature evolution under very different idealised external forcing, such as instantaneously quadrupled CO 2 , or 1% compounded increases (Geoffroy et al., 2013a(Geoffroy et al., , 2013bGregory et al., 2015) or even under more complex scenarios (Watterson, 2000;IPCC, 2013). Indeed, they have shown surprising robustness and utility in representing global temperature changes out to 2300 and at extreme values of up to 12 K (Palmer et al., 2018). Although much less work has been directed at understanding CGCM variability using such models, one recent study has attempted to use simple models to link inter-annual variability with climate change in CMIP5 CGCMs (Cox et al., 2018).
A key point of interest in the present study is to explore the relative roles of radiative climate feedbacks versus oceanic heat uptake in setting the level of climate variability and the magnitude of the annual to multi-decadal response to external forcing. There is now considerable evidence supporting a significant role for radiative feedbacks in influencing inter-annual and decadal variability (Minschwaner and Dessler, 2004;Hall and Qu, 2006;Dessler, 2010Dessler, , 2013Qu and Hall, 2014;Andrews et al., 2015;Brown et al., 2015;Gregory and Andrews, 2016;Ying and Huang, 2016;Zhou et al., 2016;Colman and Hanson, 2017). Furthermore, radiative feedbacks also determine equilibrium climate sensitivity (ECS) (Bony et al., 2006) and, an associated measure of climate sensitivity, the transient climate response (TCR) (Randall et al., 2007), which is a function of both feedbacks and oceanic heat uptake. If links could be established between the magnitude of inter-annual/decadal/multi-decadal climate variability and ECS/TCR, then this has the exciting potential to provide constraints on future climate change based on climate variations derived from the instrumental record or palaeo-reconstructions (Colman and Power, 2018;Cox et al., 2018).
With this ultimate goal in mind, and following the discussion above, we develop a simple two-layer EBM of the climate system, and then we examine the 'natural variability' evident in the EBM that arises in response to stochastic radiative forcing. We will obtain analytical solutions for the second-moment statistics of the model (i.e. the global mean surface temperature variance) under stochastic forcing. We will then explore the sensitivity of the EBM variability on different timescales to the parameter settings. This will help us to understand the relative importance of the different physical processes included in the EBM, such as radiative feedbacks, mixed layer depth and deep oceanic mixing. We will then use the EBM to try and better understand the reasons why current CGCMs exhibit such a wide range in globally averaged surface air temperature (GSAT) and variability.
In summary, we emphasise that we are not seeking a fully detailed representation of variability in the real world or in GCMs, but rather seeking to understand the simplest possible model that can represent some of the key features of global-scale variability (and response to external forcing) and understand the characteristics of that model. Section 2 describes the EBM, derivation of analytic solutions and sensitivity analysis methods. Section 3 will present the sensitivity analysis, discussion and comparison with CMIP5 results, and Section 4 will present the conclusions.

Two-layer energy balance model
To represent the earth's climate system, we formulate a two-layer EBM that consists of a pair of coupled linear subsystems (Gregory and Mitchell, 1997;Gregory, 2000;Held et al., 2010;Geoffroy et al., 2013aGeoffroy et al., , 2013b. The first subsystem represents the upper portion of the climate system: the atmosphere, the land surface and the ocean mixed layer combined, while the second subsystem represents the deep ocean. These subsystems are characterised by specific effective heat capacities C and C D such that C ( C D . Their states are defined by the temperatures T and T D , which are globally averaged anomalies with respect to the initial state and which satisfy the following equations: where k is a climate radiative feedback parameter (in Wm À2 K À1 ), c is a parameter that describes the coupling strength between the two subsystems (in Wm À2 K À1 ), while F d and F s are deterministic and stochastic radiative forcing terms, respectively (in Wm À2 ).
A key motivation of using such a model is that similar two-layer models have been considered and analysed in a number of papers considering the response to forced climate change. For example, Geoffroy et al. (2013aGeoffroy et al. ( , 2013b obtained and discussed the general analytical solutions of the two-layer model for different hypothetical climate forcing scenarios, and suggested the approach of calibrating the model parameters to imitate the time response of CGCMs from CMIP5 to radiative forcing. Gregory et al. (2015) analysed the two-layer model and approximations to (1) and (2), and discussed the TCR, the global mean surface air temperature change T under two scenarios, one with a step forcing (the abrupt 4xCO 2 experiment) and one with the 1pctCO 2 scenario (atmospheric CO 2 increasing at 1% per year). Despite the fact that the two-layer model is one of the simplest tools to mimic climate dynamics under external radiative forcing, it was able to simulate the evolution of average global surface temperature over time in response to both abrupt and time-dependent forcing reasonably accurately (Geoffroy et al., 2013b).
Indeed, the degree of complexity of the EBM is chosen in the present study very deliberately. The advantages, quite apart from the close links to EBMs used in climate change, include the simplicity, and the evidence that stochastic forcing dominates global decadal variability (as discussed in the Introduction). Further refinements could of course be added, for example as discussed by Held et al. (2010), the EBM could be modified to include features such as a timevarying value for the 'efficacy' factor for deep ocean heat uptake (see e.g. Winton et al., 2010). Geoffroy et al. (2013b), however, found impacts of including such uptake 'efficacy' were small over decadal to multi-decadal timescales, and they will not be examined further here. Note also that since we wish to keep the forcing and response at its simplest (again consistent with the way the model is used for studying secular climate change)we do not consider the annual cycle in forcing or temperature variation.
In accordance with the theory of sensitivity in dynamical systems (e.g. Eslami, 1994;Rozenvasser and Yusupov, 2000;Cacuci, 2003), the inverse of the climate feedback parameter a ¼ 1=k is referred to here as the 'climate sensitivity parameter' or 'climate sensitivity coefficient' (in Wm À2 K À1 ). Meanwhile, following the Intergovernmental Panel on Climate Change (IPCC, 2013) terminology, 'climate sensitivity' is defined as the equilibrium response of the global mean surface temperature to a doubling of the atmospheric CO 2 concentration. To avoid confusion, herein we will refer to this as the 'ECS'. To understand the role of climate feedbacks, it is useful to introduce a 'reference climate sensitivity parameter', a 0 with only the 'Planck' response from surface temperature operating (Roe, 2009). The Plank response is the radiative cooling that would be obtained if nothing changed in the warming atmosphere, apart from uniform heating through the depth of the atmosphere (of the magnitude of the surface warming), without changes to water vapour amounts or clouds or to surface albedo (Bony et al., 2006). We can write: where e ¼ 0:62 is the emissivity of the Earth (Karper and Engler, 2013), r ¼ 5:67 Â 10 À8 kg s À3 K À4 the Stephan-Boltzmann constant, T 0 ¼ 287K a 'reference' global mean surface temperature. The corresponding k 0 ¼ 1=a 0 % 3:13K W m À2 ð Þ À1 . We introduce the dimensionless feedback factor f , which represents the fraction of the Planck cooling that is offset by the radiative feedbacks in the system, specifically those due to changes in water vapour, lapse rate, clouds and surface albedo (Schlesinger, 1985;Bony et al., 2006). Not all of these individual feedbacks are positive (e.g. lapse rate, Bony et al., 2006). Taken together, however, these radiative feedbacks are 'positive' in that they offset the Planck cooling and therefore amplify the response of the system to perturbation (Randall et al., 2007). To reiterate, however, the net radiative response of the system (the Planck cooling plus the positive feedbacks) must be a dampening one, otherwise the system is unstable to variations. Since the Planck response is indeed large, this is the case both in the real world and in climate models (Flato et al., 2013). We can represent the parameter k as follows: (Schlesinger, 1985).

Deterministic and stochastic radiative forcing approximation
We assume that the deterministic radiative forcing F d is due to changes in greenhouse gas (and specifically CO 2 ) atmospheric concentrations. For the classic 1% compounding CO 2 experiment, we assume F d is linear function of time (i.e. assuming a logarithmic relationship between F d and CO 2 concentration), which gives (Geoffroy et al., 2013b): Bony et al., 2006). Function (4) provides the 1% yr À1 growth in CO 2 concentration until To describe the natural variability, we include an additive stochastic radiative forcing F s ðtÞ (Hasselmann, 1976) on the right-hand side of Equation (1). Note that the random term represents additive noise if it does not depend on the model state variables. Ideally, we would have a formulation that describes the time dependence of stochastic forcing standard deviation. Considering the 'multiplicative' noise, we would need to introduce additional parameters, which affect the amplitude of stochastic forcing. The choice of extra parameters will generate additional complexity in our simple model but not guarantee better results, so consideration of 'multiplicative' noise is left to a separate study. Random fluctuations of the atmospheric heat flux are modelled by a Gaussian white noise with zero mean hF s ðtÞi ¼ 0 and correlation function given by: where r 2 s is the variance of the noise, d is the Dirac delta function, and angular brackets denote ensemble averaging.
Since our objective is to study climate variability on annual, decadal and multi-decadal (30 years) timescales, the normally distributed stochastic forcing is assumed to be uncorrelated on monthly timescales. Estimates of stochastic forcing from the CMIP5 CGCMs were calculated by: (1) first detrending pre-industrial (PI) experiment annual mean temperatures and top of atmosphere (TOA) and surface radiation; (2) removing the mean seasonal cycles; and (3) removing TOA and surface radiation fluctuations correlated with global mean temperature. This latter removal was carried out on the assumption that radiative perturbations on monthly timescales contain a component (particularly in the long wave) which is in response to the temperature perturbations and is therefore not 'forcing' the system. This component was removed by linearly regressing radiation change against temperature, then subtracting the temperature-dependent fraction. The fraction removed was typically very small, at around 5% of the total standard deviation of net radiation on average across the CGCMs.
The remaining variation is taken to represent the stochastic forcing in the models and is found to be dominated by shortwave variations (not shown). As the result, we have obtained the value of 0.60 Wm À2 at the TOA and $0.46 W m À2 at the surface for the monthly timescale r s . The latter will be used in the calculations that follow. Observational estimates based on CERES (Clouds and the Earth's Radiant Energy System) satellite data indicate that global-scale total TOA variability has a standard deviation of around 0.62 (0.28) W m À2 on monthly (annual) timescales (Trenberth et al., 2014) a value comparable to the multi-model mean. In parallel, analysing detrended surface temperature time series in the CGCMs, we calculated the values of 0.00701, 0.00212 and 0.00071 K 2 for PI inter-annual, decadal and multi-decadal (30 years) variances of the global mean surface temperature, respectively. These values will be used in the analysis of two-layer model discussed in this article.

Model parameters
The two-layer model has two state time-dependent variables T and T D , and six free parameters: C; C D ; f ; c; g and r s . Since the values of parameters g and r s are already discussed above, here we will pay attention to the remaining four parameters.
Note that we would expect parameter values to determine the 'timescale' of the temperature response to forcing. Specifically, there are two timescales set by the relaxation times s A and s O for the equilibration of upper and deep ocean temperatures, respectively (Geoffroy et al., 2013a). Critical parameters that define s A and s O are, correspondingly, the effective heat capacities C and C D . For inter-annual variability, we assume that C % 7:5 W yr m À2 K À1 (Geoffroy et al., 2013a). The corresponding depth of the upper ocean layer h is about 60 m since C ¼ c 0 h, where c 0 ¼ 4 Â 10 6 J m À3 K À1 is the volumetric heat capacity (Gregory, 2000). For decadal variability, we assume that C % 20W yr m À2 K À1 (Gregory, 2000;Schwartz, 2007), which gives the estimate h % 150m. Note that these differences in the upper ocean reflect differing physical processes operating on these timescales. On shorter timescales, relatively shallow mixed layer processes dominate the effective depth of the ocean heat exchange, whereas on longer timescales overturning circulations, vertical advection and convective mixing at high latitudes cause an exchange in heat at much deeper levels in the ocean (Power and Hirst, 1997;Rhein et al., 2013).
The question arises on possible variations of these parameters on multi-decadal timescales. In a stationary climate, there is no a priori reason to expect time dependence of the parameters, and we are assuming a stable PI climate (compared with changes under secular climate change). This assumption of unchanging variables is also consistent with the assumptions elsewhere, for example under the modelling of Geoffroy et al. (2013a) and Palmer et al. (2018). Investigation of introducing an 'efficacy' factor by Geoffroy et al. (2013b) for deep ocean processes (relevant for multi-decadal timescales) found little impact on global temperature changes at less than centennial timescales, so would not likely impact the results here. Nevertheless, there is an assumption here of stationarity. In fact, true climate stationarity is unlikely to exist in nature (Masson-Delmotte et al., 2013). Furthermore, even in a hypothetical 'stationary' climate, variances are likely to themselves vary on very long timescales due to the internal dynamics of the climate system. Climate models show temporal variability in variance and other climate statistics, from processes such as inter-decadal to inter-centennial variability of ENSO (Wittenberg, 2009). Such temporal variations are beyond the scope of the present study, however, so are not investigated further here.
The values of C D and c are taken in accordance with values consistent with the CMIP5 multi-model mean under climate change derived by Geoffroy et al. (2013a): C D % 100 W yr m À2 K À1 (this value corresponds to an equivalent deep ocean layer depth h D % 800 m) and c % 0:7 Wm À2 K À1 .
The feedback factor, f, plays an influential role in the behaviour of the coupled dynamical system. We assume that 0< f <1, meaning that the feedback amplifies the system response to external forcing (Colman, 2003;Roe, 2009), and that the system is stable to external forcing. In practice, in CGCMs the net feedback factor is the result of contributions from temperature-dependent changes to water vapour, lapse rate, surface albedo and clouds (Colman, 2003;Bony et al., 2006). Using the multi-model mean of the climate radiative feedback parameter k ¼ 1:13 W m À2 K À1 (Geoffroy et al., 2013a) and applying the relationship between f and k (see Section 2.1), we obtain f % 0.64.

Mean and second-moment equations
Equations (1) and (2) can be represented as where nðtÞ is a Gaussian white noise with unit variance, The associated Fokker-Planck or forward Kolmogorov equation that describes the time evolution of the probability density function P P T; T D ; t ð Þfor these equations is of the form: Multiplying both sides of Equation (9) by T, and then integrating over all possible values of T and T D , one can derive the mean equation for variable T (Nicolis, 1998;Heinz, 2011): The partial derivative by t is replaced here by the regular derivative because the expectation values hTi and hT D i of variables T and T D are only functions of t. Similarly, we can also obtain the second equation, which describes the time evolution of the mean value hT D i The second-moment equation for the variable T can be derived by multiplying both sides of the Fokker-Planck Equation (9) by T 2 and then integrating by parts: where hdT 2 i ¼ hT 2 iÀhTi 2 is the surface temperature variance, dT ¼ TÀhTi and dT D ¼ T D ÀhT D i. Using the same procedure, we can obtain the following second-moment Equations: For the steady-state probability distribution, the Equations (12)-(14) take the form: hdT 2 i is then given by: It follows from Equation (18) that in the EBM, the surface temperature variance (variability) is independent of the deterministic radiative forcing F d or, in other words, on concentrations of atmospheric greenhouse gases such as CO 2 . Therefore, the parameter g, which defines deterministic radiative forcing by (4), can be excluded from our consideration.
Note that the study of the variance of the model is confined to that in the PI climate (i.e. without external forcing). A full analysis of the evolution of variability under climate change (e.g. using the approach in Majda and Gershgorin (2010)) would be a useful extension in a further study.

Sensitivity analysis method
Model parameters, denoted by the vector p ¼ C; C D ; f ; c; r s ð Þ T , influence the model output and, in particular, the surface temperature variance defined by Equation (18). Most of these parameters are empirical and therefore not uniquely defined. However, the CGCM and observational results provide plausible ranges (e.g. on radiative feedbacks). The impact of varying a parameter on the EBM output itself varies from parameter to parameter. For exploring the influence of the uncertainty in the model parameters on the climate variability, we employ sensitivity coefficients defined as partial derivatives of surface temperature variance dT 2 h i with respect to parameters (e.g. Eslami, 1994;Rozenvasser and Yusupov, 2000): where p i is any particular parameter of interest and p Ã i its base (nominal) value around which the sensitivity is estimated. Sensitivity coefficients can be readily determined by using (19) with respect a given parameter, keeping other parameters unchanged.
Sensitivity coefficients for different model parameters, however, can naturally be expected to differ in magnitude, simply because of differences in their units. For assessing, and ranking, the relative influence of parameter uncertainties on the climate variability, we use the relative sensitivity coefficients defined as The relative sensitivity, then, describes how temperature variance is affected by fractional changes in a particular variablefor example, how a 10% change in mixed layer depth changes the variance, which can then be compared, say with a 10% change in the ocean deep layer.

Results and discussion
To evaluate the two-layer model and in particular to rank the impact of model parameters on the surface temperature variance dT 2 h i , we use the so-called monothetic analysis, or one-variable-at-a-time (OVAT), method, changing each parameter over its range, and keeping others at their 'base' valuesthat is, changing only one parameter at a time (e.g. Daniel, 1973;Murphy et al., 2004). The downside of this approach is its inability to explore interactions among the parameters. Since in our simple model we do not explicitly account for interactions between parameters and the number of parameters is small, the use of OVAT method is a good starting point. The OVAT approach requires maximum and minimum values, and base values for all model parameters. Geoffroy et al. (2013a) derived ranges of parameter values for the two-layer model from analysis of the CMIP5 models under simplified forcing (instantaneous CO 2 quadrupling). For reference, we display the three values for each parameter used by Geoffroy et al. (2013a) in Table 1. Note that the range of radiative feedback parameter k (Table 1) corresponds to the range of f max % 0:8 and f min % 0:47 for the dimensionless feedback factor f, which is used in our model instead of k. We also want to explore sensitivities of our model to parameters outside this range, both for reference purposes, and also acknowledging that credible values in the climate system may lie outside the values found by Geoffroy et al. (2013a). Table 2 lists the maximum, Table 1. Maximum, minimum and multi-model mean values of radiative feedback parameter k and climate system inertia parameters C, C D and c obtained from analysis of the CMIP5 models (Geoffroy et al., 2013a)  minimum and base values of all two-layer model parameters used in calculations. Figures 1 and 2 give a general representation on the impact of all model parameters on the surface temperature variance dT 2 h i . Before proceeding to the analysis of these figures, it is important to note that both the magnitude of calculated surface temperature variance and the estimated impact of one or another model parameter on the behaviour of dT 2 h i in the parameter space strongly depend on the 'climate timescale' considered.
To explore how the model reproduces climate variability from annual to multi-decadal timescales under the influence of stochastic radiative forcing, we calculated the surface temperature variance as a function of the upper layer heat capacity. This was chosen, as we would expect this to set the timescale of response, over short time periods at least to vary sensitively with this variable (Gregory, 2000). Figure 1a (red curve) presents the dependence of dT 2 h i on the parameter C, assuming that  all other parameters hold their base values. This figure shows that the variance dT 2 h i is subject to decrease when the parameter C is growing, that is when the timescale of climate processes increases. (Equation (18) reveals that decreasing rate is in fact proportional to 1/C 2 , since kC D ) cC.) Thus, the model behaves qualitatively as we would intuitively expect: (1) the dependence of the surface temperature variability on the climate timescale is nonlinear and (2) the magnitude of dT 2 h i decreases when the climate timescale increases (see Manabe and Stouffer (1996) and Section 2.2). Recall that in our model the timescale is defined by the effective heat capacity of the upper layer (e.g. Schwartz, 2007;Knutti et al., 2008). Figure 1a also shows four additional curves that correspond to the values of 0.2, 0.4, 0.75 and 0.8 for the dimensionless feedback factor f . For a given 'climate timescale' (i.e. value of C), the surface temperature variance increases when the feedback factor f increases. Thus, analogous to the case for forced climate change (e.g. Bony et al., 2006), the uncertainty present in the feedback factor f (e.g. Colman and Hanson, 2017) has a marked effect on the modelled climate variability, particularly for high values of f. Figure 1b shows that the impact of the parameter f on the surface temperature variance depends on the climate timescale considered: the shorter the timescale, the greater the surface temperature variance (at a given value of the parameter f ). This seems to suggest that uncertainties in radiative feedbacks are more important for driving differences in short timescale variability than long. This matter will be discussed further below, when we analyse relative sensitivity coefficients.
The heat exchange parameter, c, also noticeably affects the climate variability, as shown in Fig. 1c, and this influence increases, as expected, as C decreases (i.e. as the depth of the upper layer decreases, or the timescale of the upper layer increases). For increasing c, there is a decrease in surface temperature variance for all timescales. This is because increasing c effectively increases the depth of water responding to the stochastic forcing.
A less significant effect on the climate variability of various climate timescales is provided by the parameter C D (Fig. 1d). For a particular climate timescale, the influence of change in C D on dT 2 h i is insignificant, except perhaps for very low values of C D (which are unrealistic in the models, given it represents the lower ocean). We might expect greater sensitivity to C D for larger values of c. Figure 2 displays the influence of radiative stochastic forcing r s on the surface temperature variance for different timescales. Since the dependence of dT 2 h i on r s is quadratic (see Equation (18)), the uncertainty in stochastic forcing can produce a significant uncertainty in the model output.
We now compare the EBM variability with variability in the CGCMs. To do this, we first set the variance of the stochastic forcing to the mean of the CGCM values. Verifying the calculated temperature variance against the PI inter-annual, decadal and multi-decadal variances (see Section 2.2), we conclude that the model can quite realistically simulate the inter-annual and decadal climate variability (see Table 3) using standard settings of all parameters, with changes to C alone. However, for these 'base' parameter values, the multi-decadal temperature variance is twice the value evident in the CMIP5 models. By changing parameters c, C D and r s , it is possible to fit the multi-decadal variability to the surface temperature variance obtained from analysis of the CMIP5 models to some extent. Table 3 includes some results calculated for different values of c, C D and r s , which, however, belong to the range of values obtained by Geoffroy et al. (2013a). This illustrates two things. First, the EBM can be 'tuned' to represent variability simulated by the CGCMs across a wide range of timescales. This is analogous to the tuning that was carried out by Geoffroy et al. (2013a) for forced climate change. Second, the model temperature variance is sensitive to the choice of parameters, but that sensitivity depends in turn on the baseline selection of those parameters. As the aim is to use the EBM to understand the relationships between variability at different timescales and response to forced climate change in CGCMs, we need to systematically explore our model sensitivity to parameter settings. That sensitivity can be estimated via 'sensitivity coefficients' (19). If dp i is a small variation in the model parameter p i , the change in the surface temperature variance D dT 2 h i induced by dp i can be estimated to a first-order accuracy in the following way: Calculated sensitivity coefficients with respect to the model parameters are illustrated in Figs. 3 and 4 for different climate timescales. As shown in these figures,  coefficients S f and S rs are positive, while all others are negative, which indicates that the infinitesimal perturbations in these parameters lead to positive (negative) changes in the surface temperature variance dT 2 h i .
Furthermore, all sensitivity coefficients are nonlinear functions of the corresponding parameters with the exception of r s (the standard deviation of the stochastic forcing). Importantly, the effects of variations in parameters C; f ; c and r s all decrease in magnitude with the climate timescale considered, with only C D increasing. This can be understood as the absolute variance decreases with timescale, so that given parameter changes are less effective in changing temperature variance. For C D , this is not the case because the deep ocean is seen less at shorter timescales. Each sensitivity coefficient, as shown in Figs. 3 and 4, varies noticeably in the range of the corresponding parameter. Therefore, the same parameter changes induce different changes in dT 2 h i at different parts of the parameter range. For example, if dC ¼ 0:3750 W yr m À2 K À1 , which is equal to 5% change with respect to the base value, then for inter-annual timescales we obtain D dT 2 h i % 3:78 Â 10 À4 K 2 , while for decadal -D dT 2 h i ¼ 1:45 Â 10 À4 K 2 since the sensitivity coefficient S C decreases in absolute value from 1:01 Â 10 À3 to 1:45 Â 10 À4 , that is decreases by almost an order of Table 4. Changes in the absolute value of surface temperature variance D dT 2 h i caused by variations in model parameters whose absolute values dp i deviate 5% from the base values p Ã i . Parameter Variation in absolute value of parameter dp i 0.3750 0.0320 0.0350 5 0.0230 Change in D dT 2 h i 3.87 Â 10 À4 4.50 Â 10 À4 1.40 Â 10 À4 6.54 Â 10 À6 7.86 Â 10 À4 magnitude. To some degree, a similar pattern emerges for the remaining parameters. For instance, when the feedback factor increases from 0.4 to 0.8, then it can be seen the increase in D dT 2 h i is almost a factor of 10from 1:36 Â 10 À4 to 1:18 Â 10 À3 K 2 because of the growth of the sensitivity coefficient from 0:68 Â 10 À2 to 2:95 Â 10 À2 . Sensitivity coefficients allow us to estimate the influence of model parameter uncertainty on the uncertainty of calculated surface temperature variance, but the absolute effects depend on the scale of the parameter being varied. To assess the relative effects of parametric uncertainty, we shall assume that the parameter variations are 5% of their base values, that is dp i ¼ 0:05 Â p Ã i . Table 4 shows how the absolute value of surface temperature variance changes when model parameters deviate 5% from their base values. It follows from this table that the changes in D dT 2 h i due to changes in the parameters C; f ; c and r s are quantities of the same order; however, the effect of variations in r s is almost twice as large as the others. The influence of variations in the parameter C D on changes in dT 2 h i is about two orders of magnitude less than affects caused by other parameters. Along with the absolute sensitivity coefficients analysed above, we also calculated relative (normalised) sensitivity coefficients S R pi which allows one to quantitatively assess the relative influence of a parameter change on the surface temperature variance and, consequently, rank the parameters in accordance with their degree of influence on D dT 2 h i . The graphs of relative sensitivity coefficients with respect to parameters C; f ; c and C D are plotted in Fig. 5. There is inherent nonlinearity in all coefficients and the coefficients themselves vary significantly in the range of their values. In contrast, the relative sensitivity coefficient S R rs is independent of the parameter r s : S R rs ¼ 2 for all r s . Table 5 presents values of relative sensitivity coefficients calculated for different climate timescales. This table shows that the stochastic forcing variability has the largest rank for all timescales. Feedback parameter and upper layer effective heat capacity rank second and third, respectively, followed by heat exchange coefficient c and lower layer effective heat capacity C D . With increasing climate timescale, the coefficients S R C and S R c decrease, while the coefficients S R f and S R CD , in contrast, increase. This confirms then that for this model the radiative feedback is more important for setting variability for long timescales than for smaller. The upper layer ocean depth also becomes less important, which is understandable as the relative importance of the deep layer increases with increasing timescale. Radiative feedback parameter remained important at all timescales.
The results discussed above are for the base values of model parameters. However, a key question is how the range of model parameters found from fitting GCM results (Geoffroy et al., 2013a) contributes to the range of variability found in the GCMs. Table 6 shows the spread in global mean surface temperature variance dT 2 h i taking the variables C; C D ; f , c and r s at their base values, then varying each parameter in turn across the range found by Geoffroy et al. (2013a) (or as derived here for r s Þ. The largest impact comes from r s (0.01 K 2 ) followed by the upper layer heat capacity C and feedback parameter k at around half the value (0.0056 K 2 and 0.0053 K 2 respectively). The impact of the diffusion parameter range is around half the value again. Finally, the range in the deep ocean heat capacity is the least significant, at 0.0002 K 2 . The values for C, c, k and C D are obtained from analysis of the CMIP5 models (Geoffroy et al., 2013a), and the values for r s are from this study.

Concluding remarks
Global mean surface temperature variances simulated by CMIP5 CGCMs deviate from each other substantially, on inter-annual, decadal and multi-decadal timescales (e.g. Colman and Power, 2018). The understanding of 'natural' causes of climate variability and the role of various physical mechanisms in the formation of climate variability therefore represent key challenges of climate science. In this article, we evaluate the effect of random radiative forcing and the relative roles of feedbacks, ocean heat uptake and climate system thermal inertia on the formation of 'PI' annual, decadal and multi-decadal climate variability. The main tool of this study is a simple two-layer energy balance model (EBM, e.g. Gregory, 2000), in which climate variability is associated with the global mean surface temperature variance, and the climate timescale is determined by the model parameters.
The EBM we use includes stochastic radiative forcing, an upper ocean layer, deep ocean heat uptake via a diffusive parameter, a deep ocean layer and radiative feedbacks. We derived the Fokker-Planck equation corresponding to the forced EBM, and we derived the analytical expression for the surface temperature variance in terms of these model parameters. We then compared the global mean surface temperature variance on inter-annual, decadal and multi-decadal timescales, with their counterparts in CMIP5 CGCMs forced under 'PI' conditions. We found that the EBM can represent the magnitude of the surface temperature variance with appropriate parameter settings, although an adjustment of the upper layer depth was important for setting the timescale of the variability.
A key aim of the study was to understand the relative importance of parameters uncertainty, and in particular, uncertainty in the radiative feedback, to the range of variability implied by the EBM. We found that for all but extremely long timescales, the most important parameters were in descending order: the magnitude of the stochastic forcing, the radiative feedback parameter, the depth of the upper ocean layer, the diffusion parameter and the depth of the deep ocean layer. However, the relative importance of these parameters varied strongly across timescales, as the timescale lengthened fractional changes in the feedback parameter and depth of the deep ocean became more important, and the depth of the upper ocean and diffusion parameter less so.
Finally, we compared the range of temperature variance that is implied by the range in model parameters obtained from the study of Geoffroy et al. (2013a) from fits of a similar two-layer model to climate change experiments (and using the range of stochastic radiative forcing derived here from the same subset of CGCMs). We found that the range of stochastic forcing implied approximately twice the range of temperature variance as implied by each of the feedback and upper ocean layer depths (which had a roughly equal magnitude). The variation in the diffusion to the deeper ocean was found to be around half as important again, and the deep ocean depth unimportant. This implies that radiative feedbacks play an important, but not dominant, role in setting the range of temperature variation across the CMIP5 CGCMs. It underscores that more research is therefore needed to understand the relationship between tropical temperature variance and climate sensitivity found by Colman and Power (2018). If this relationship is found to be robust, then radiative feedbacks may be critical in setting the magnitude of both climate sensitivity and tropical temperature variance.
Of course, there are important caveats to our conclusions. The most important is of course that the EBM is highly simplified and while we can tune the EBM to replicate CGCM variability, the EBM might be simulating the right level of variability for the wrong reasons. Another is that despite the similarity in overall magnitude of radiative feedbacks under climate change and variability (Colman and Hanson, 2017), they may differ sufficiently in pattern that they drive very different global temperature response to variability and climate change (Xie et al., 2016). Also of course, the effect of the parameter changes on variance will depend on multi-parameter combinations, and not on the changes to single parameters undertaken here. Nevertheless, the results here indicate that further research is needed to understand the relationship between climate sensitivity, climate variability and the time rate of temperature change during transient climate change (e.g. as measured by the TCR). A further study will address these issues. responsible for CMIP, and we thank the climate modelling groups for producing and making available their model output. For CMIP the US Department of Energy's Program for Climate Model Diagnosis and Intercomparison provides coordinating support and led development of software infrastructure in partnership with the Global Organization for Earth System Science Portals.

Disclosure statement
No potential conflict of interest was reported by the authors.