1. Introduction

Global temperatures are characterised by large global-scale ‘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₂, or 1% compounded increases (Geoffroy et al., 2013a, 2013b; Gregory 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, 2010, 2013; Qu 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.

2. Methods

3. Results and discussion

To evaluate the two-layer model and in particular to rank the impact of model parameters on the surface temperature variance $〈{\delta T}^2〉$, 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’ values – that 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₂ 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 $\lambda$ (Table 1) corresponds to the range of $f_{max}\approx 0.8$ and $f_{min}\approx 0.47$ for the dimensionless feedback factor f, which is used in our model instead of $\lambda$. 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, 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 $〈\delta T^2〉$. 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 $〈\delta T^2〉$ 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 $〈\delta T^2〉$ on the parameter C, assuming that all other parameters hold their base values. This figure shows that the variance $〈\delta T^2〉$ 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², since λC_D ≫ γC.) 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 $〈\delta T^2〉$ 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, γ, 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 γ, there is a decrease in surface temperature variance for all timescales. This is because increasing γ 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 $〈\delta T^2〉$ 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 γ.

Figure 2 displays the influence of radiative stochastic forcing ${\sigma }_s$ on the surface temperature variance for different timescales. Since the dependence of $〈\delta T^2〉$ on ${\sigma }_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_D and ${\sigma }_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_D and ${\sigma }_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 $\delta p_i$ is a small variation in the model parameter $p_i$, the change in the surface temperature variance $\Delta 〈\delta T^2〉$ induced by $\delta p_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_{{\sigma }_s}$ 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 $〈\delta T^2〉$. Furthermore, all sensitivity coefficients are nonlinear functions of the corresponding parameters with the exception of ${\sigma }_s$ (the standard deviation of the stochastic forcing). Importantly, the effects of variations in parameters $C, f, \gamma$ and ${\sigma }_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 $〈\delta T^2〉$ at different parts of the parameter range. For example, if $\delta C=0.3750\hspace{0.17em}\mathrm{W}$ yr m⁻² K⁻¹, which is equal to 5% change with respect to the base value, then for inter-annual timescales we obtain $\left|\Delta 〈\delta T^2〉\right|\approx 3.78\times {10}^{-4}$ K², while for decadal – $\left|\Delta 〈\delta T^2〉\right|=1.45\times {10}^{-4}$ K² since the sensitivity coefficient $S_C$ decreases in absolute value from $1.01\times {10}^{-3}$ to $1.45\times {10}^{-4}$, that is decreases by almost an order of 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 $\left|\Delta 〈\delta T^2〉\right|$ is almost a factor of 10 – from $1.36\times {10}^{-4}$ to $1.18\times {10}^{-3}$ K² because of the growth of the sensitivity coefficient from $0.68\times {10}^{-2}$ to $2.95\times {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 $\delta p_i=0.05\times 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 $\left|\Delta 〈\delta T^2〉\right|$ due to changes in the parameters $C, f, \gamma$ and ${\sigma }_s$ are quantities of the same order; however, the effect of variations in ${\sigma }_s$ is almost twice as large as the others. The influence of variations in the parameter C_D on changes in $〈\delta T^2〉$ 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_{p_i}^R$ 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 $\Delta 〈\delta T^2〉$. The graphs of relative sensitivity coefficients with respect to parameters $C, f, \gamma$ 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_{{\sigma }_s}^R$ is independent of the parameter ${\sigma }_s$: $S_{{\sigma }_s}^R=2$ for all ${\sigma }_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 $\gamma$ and lower layer effective heat capacity $C_D$. With increasing climate timescale, the coefficients $S_C^R$ and $S_{\gamma }^R$ decrease, while the coefficients $S_f^R$ and $S_{C_D}^R$, 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 $〈\delta T^2〉$ taking the variables $C, C_D, f$, $\gamma$ and ${\sigma }_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 ${\sigma }_s)$. The largest impact comes from ${\sigma }_s$ (0.01 K²) followed by the upper layer heat capacity $C$ and feedback parameter $\lambda$ at around half the value (0.0056 K² and 0.0053 K² 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².

4. 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.