A simple and robust framework is proposed for the partitioning of the different components of internal variability and model uncertainty in an unbalanced multimember multimodel ensemble (MM2E) of climate projections obtained for a suite of statistical downscaling models (SDMs) and global climate models (GCMs). It is based on the quasi-ergodic assumption for transient climate simulations. Model uncertainty components are estimated from the noise-free signals of the different modeling chains using a two-way analysis of variance (ANOVA) framework. The residuals from the noise-free signals are used to estimate the large- and small-scale internal variability components associated with each considered GCM–SDM configuration. This framework makes it possible to take into account all members available from any climate ensemble of opportunity. Uncertainty is quantified as a function of lead time for projections of changes in temperature and precipitation produced for a mesoscale alpine catchment. Internal variability accounts for more than 80% of total uncertainty in the first decades. This proportion decreases to less than 10% at the end of the century for temperature but remains greater than 50% for precipitation. Small-scale internal variability is negligible for temperature; however, it is similar to the large-scale component for precipitation, whatever the projection lead time. SDM uncertainty is always greater than GCM uncertainty for precipitation. It is also greater for temperature in the middle of the century. The response-to-uncertainty ratio is very high for temperature. For precipitation, it is always less than one, indicating that even the sign of change is uncertain.
a. Uncertainty sources
A critical issue in climate change impact studies is the estimation of uncertainties associated with future projections along with the estimation of the contribution of the different uncertainty sources. As pointed out by Dobler et al. (2012), very different sources of uncertainty are involved including scenario uncertainty, model uncertainty, and model internal variability.
Scenario uncertainty is related to the poorly known future of greenhouse gas emissions. Model uncertainty, also termed response uncertainty, corresponds to the dispersion between the different climate responses classically obtained with different models for the same forcing configuration. It is due to the limitations of the model structure and parameterization used to represent geophysical processes. Model uncertainty clearly concerns global climate models (GCMs) used to simulate at the global scale the response of the climate system to atmospheric composition perturbations. It also concerns regional downscaling models (RDMs) used to derive corresponding climate scenarios at the regional scale. Model uncertainty and scenario uncertainty are classically explored via multimodel and multiscenario experiments (Solomon et al. 2007).
The internal variability of a given GCM–RDM chain supposedly represents the natural variability of regional climate at daily to multidecadal time scales (e.g., Karoly and Wu 2005). Because of the chaotic and nonlinear nature of atmospheric processes, this variability has long been observed even in a stationary climate (Madden 1976). In a nonstationary climate, this variability can remain high above the trend related to a given forcing (e.g., greenhouse gases and aerosols) (e.g., Hawkins and Sutton 2011; Deser et al. 2012; Lafaysse et al. 2014). The internal variability of a given GCM/RDM chain can actually be partitioned into its large-scale and local-scale components (Braun et al. 2012; Lafaysse et al. 2014). The large-scale component can be attributed to the chaotic variability of the climate at the global scale. Mainly produced by the GCM itself, it corresponds to the variability obtained from a given GCM experiment over a long time period for a stationary climate (e.g., Räisänen 2001; Deser et al. 2012). Very similar large-scale atmospheric circulation patterns can then lead to very different meteorological observations at the local scale. This local-scale component of natural variability is in principle mainly accounted for by the RDM. For instance, when a statistical downscaling model (SDM) is used as an RDM, it classically includes a stochastic process that produces several possible meteorological scenarios for a given large-scale pattern or a given sequence of large-scale patterns (Buishand and Brandsma 2001; Mezghani and Hingray 2009). The local-scale internal variability, sometimes referred to as the RDM’s internal variability, corresponds to the dispersion between these scenarios.
b. Partitioning and estimating uncertainties
Estimating the total uncertainty of projections from a multimember multimodel ensemble (MM2E) of climate experiments and estimating the contribution of each particular uncertainty source requires an appropriate statistical framework. Various methods have been proposed for this in the recent years. Most are based on empirical statistical analysis or more formal analysis of variance (ANOVA) of projections obtained for a particular MM2E opportunity (Hawkins and Sutton 2009; Hingray et al. 2007; Räisänen 2001; Yip et al. 2011). In all cases, a major difficulty lies in the limited number of members that have until now been classically available for most modeling chains and in particular for GCMs. When a large number of members is available for a given chain (e.g., GCM), the climate response of the chain for the considered projection lead time can be estimated from to the multimember mean of the projections. Similarly, the internal variability of the chain can be estimated from the intermember variance of the projections (e.g., Deser et al. 2012). When only a few members are available, the climate response of the chain and the effect of its internal variability are actually difficult if not impossible to separate, especially when the internal variability is nonnegligible compared to the chain’s climate response. In such a case, significantly biased estimates are likely to be obtained for the different uncertainty components, total uncertainty, and significance of projected changes (e.g., Deser et al. 2012; Sansom et al. 2013; Hingray and Saïd 2014, manuscript submitted to Geophys. Res. Lett.).
Another consequence of a small number of members is that these estimates may also strongly vary from one projection lead time to another, depending on the projection values simulated for each lead time. Nonnegligible temporal variations in GCM internal variability were for instance obtained over the next century by Yip et al. (2011) from regional warming projections based on 21 future climate experiments with two runs each. Such temporal variations are expected to be much greater for other variables presenting a nonnegligible internal variability such as precipitation (Hingray and Saïd 2014, manuscript submitted to Geophys. Res. Lett.). There is, however, obviously no reason to expect abrupt changes in internal variability in a transient climate.
In recent years, long time series have become available for the large majority of GCM experiments and in turn for a large number of GCM–RDM chains. GCM experiments from the stream2 data of the Ensemble-Based Predictions of Climate Changes and Their Impacts (ENSEMBLES) European research project cover, for instance, 140 years of the preindustrial period (1860–2000) and the 100 years of the next century (Johns et al. 2011). This provides a major opportunity to significantly improve the estimation of internal variability and model uncertainty components, total uncertainty, and in turn the significance of projected changes.
The aim of this paper is to present a theoretical framework for the partitioning of the total uncertainty of climate projections into the four following components: GCM uncertainty, SDM uncertainty, and large-scale and local-scale internal variability of the GCM–SDM modeling chains. The method in particular makes use of all data and simulation experiments available for the target region, whatever the number of runs available for instance for each GCM. For illustration, the method is applied to hydrometeorological projections based on an MM2E of experiments developed for a mesoscale alpine catchment as part of the Regional Climate, Water, Energy Resources and Uncertainties from 1960 to 2030 (RIWER2030) research project (Hingray et al. 2013).
Section 2 describes the regional climate experiments and explains some of the practical choices made for processing the data. Section 3 introduces the two-part statistical framework proposed for the analysis (basically an analysis of variance). Some technical developments are relegated to the appendixes. Results are presented in section 4. Section 5 discusses the results and presents our conclusions. The data and our code for the calculations in this application are available online, written in the Matlab computing environment (http://www.lthe.fr/RIWER2030/).
We illustrate the method by applying it to projections of the change in 20-yr means of precipitation and temperature obtained from the MM2E of opportunity derived within the RIWER2030 research project (Hingray et al. 2013; Lafaysse et al. 2014). The MM2E is extracted from an original ensemble of daily time series simulated for the years 1860–2099.
Projections were obtained for the upper Durance River catchment, a 3850 km2 wide drainage basin located in the southern French Alps. They come from Ns = 6 multivariate SDMs forced by the outputs of Ne = 11 climate simulations. The simulations correspond to Ng = 5 GCMs, an ensemble of three runs being available for three GCMs. They were obtained from the STREAM2 experiment conducted within the ENSEMBLES European project under Twentieth-Century Climate in Coupled Models (20C3M) historical forcing (with constant solar and volcanic forcing) for the 1860–2000 historical period and the Special Report on Emissions Scenarios (SRES) A1B scenario for the 2000–2100 future period (Johns et al. 2011). The GCMs, SDMs, and corresponding references are listed in Tables 1 and 2. The particular GCMs used for the analysis were chosen on the basis of data availability.
For each GCM–SDM chain, an ensemble of Nk = 100 scenarios is available, resulting from the stochastic generation process associated with each SDM (Lafaysse et al. 2014). This gives a total of Ns × Ne × Nk = 6600 times series of multivariate meteorological scenarios for 1860–2099. Figure 1 shows time series of 20-yr mean temperatures and precipitation for selected GCM–SDM configurations and members.
As proposed by Hingray et al. (2007), the uncertainty analysis is carried out on future to reference control period changes X, expressed in terms of absolute changes for temperatures and in terms of relative changes for precipitation:
where MT and MP are the 20-yr interannual mean values Y of the raw meteorological projections, where t refers to a given future or past 20-yr period centered on year t and where c refers to the 20-yr reference control period centered on year c. In the present work, the reference control period is the 1980–99 period. The choice to work on change variables rather than on raw projections follows the choice made for most impact studies. More confidence is actually usually given to changes variables assumed to allow for partially coping with biases in simulation models (e.g., GCMs and SDMs).
3. Theoretical framework
a. The quasi-ergodic assumption
In recent years, long time series of control simulations plus transient projections have become available for the large majority of GCM experiments. Time variations resulting from a GCM over the whole simulation period are a combination of the climate response of the GCM and its internal variability. Similarly, time variations of a regional variable obtained over the period for a given GCM–SDM chain are the combination of the climate response of the chain and of its total internal variability.
Let us consider that the studied variable is the interannual mean over a given n-yr subperiod of a given regional climate variable. Let us next consider one simulation member obtained with a given GCM–SDM chain for a stationary climate (e.g., a climate in the absence of any anthropogenic forcing or any external radiative forcing). Provided that the simulation period is long enough, the time average of this simulation for the studied variable corresponds to the climate response of the GCM–SDM chain for this variable in this climate context. It would also correspond to its ensemble average for any n-yr subperiod of the total simulation period, where the ensemble is constituted by all members that could be achieved with this GCM–SDM chain for this subperiod. For any member, the statistical distribution of the climate variable over time would similarly define the internal variability of the chain for this variable. It would also correspond, for any subperiod of the total simulation period, to the distribution of the variable that would be achieved from an infinite number of simulation members.
Such a dynamical system, for which the behavior averaged over time is the same as that averaged over the space of all the system’s states, is said to be ergodic. For an ergodic system, the time average of one sequence of events is the same as the ensemble average at a given time. The same applies also for higher-order statistical moments such as variance.
If long time series were available for a stationary climate context for all GCM–SDM chains, the climate response and the internal variability of each chain for a given time would be therefore easy to estimate assuming that each chain behaves as an ergodic random process. Because solar and volcanic forcings were held constant in ENSEMBLES stream 2 experiments, this would here appear to be a reasonable assumption for the 1860–2000 extended control period (or at least for a relevant part of it; see section 3b) for which GCM–SDM control experiments were produced. When 20-yr means are considered as the studied change variable and when stationarity can be assumed for the whole extended control period, the GCM–SDM climate response and its internal variability can be estimated from a sample of size 7 (140 yr/20 yr).
A sample with more data can be obtained assuming that the process is “quasi-ergodic” for any transient period over which solar and volcanic forcings are held constant. It seems reasonable to consider that if the climate response of each GCM–SDM chain varies over the transient period, this variation should be gradual and smooth, the higher-frequency variations of the time series being due to internal variability alone (Fig. 2). Another reasonable assumption is that the internal variability remains constant over the period or that it varies gradually as a linear function of the climate response of the GCM–SDM chain. The first assumption was used for instance by Hawkins and Sutton (2011) and Räisänen (2001) and in the present study for changes in surface temperature. The second assumption is used in the present study for changes in precipitation (see appendix B).
Under this quasi-ergodic assumption, extracting of the noise-free signal (NFS) from the time series makes it thus possible to estimate the climate response of the GCM–SDM chain and its possible change with time (the climate change response of the chain), as well as the noise around this response and in turn the GCM–SDM internal variability. In the present case, as simulations are provided for the 1860–2100 period, the GCM–SDM internal variability can be estimated from a sample of size 12 (240 yr/20 yr).
Finally, note also that when N different times series members are available for a given GCM–SDM chain, the quasi-ergodic assumption allows the estimation of both the GCM–SDM climate response and its internal variability from the “super sample” constituted by all subperiod data of all times series members (Fig. 2, modeling chain B). In the present case, and even if only one stochastic realization is considered, the size of this super sample is for instance 12 × 3 = 36 for GCM–SDM chains with three available GCM runs. As shown by Hingray and Saïd (2014, manuscript submitted to Geophys. Res. Lett.), the estimation of both the climate response and the internal variability of the chain is in this case expected to be much more robust than when estimated from the small number of members available for the considered lead time, even if internal variability is significant when compared to interchain dispersion.
These considerations are the base of the quasi-ergodic ANOVA framework used for partitioning the uncertainty sources for the present dataset. Let us note X(g, s, r, k, t), the simulation outputs of the studied change variable for a given run r of each GCM g, each stochastic generation k of SDM s, each year t. These outputs can be written as
where NFS(g, s, t) is the NFS of the change variable for GCM–SDM chain g–s and where η(g, s, r, k, t) are the residuals of the kth stochastic generation of SDM s for the rth run of the GCM–SDM chain g–s.
The total uncertainty of the change variable, [X(g, s, r, k, t)], corresponds to the sum of the variances of both terms of the right-hand side of Eq. (3). They respectively correspond to the model uncertainty and to the internal variability of X for the modeling chains. The following framework, summarized in Fig. 3, details how to estimate the different components of both terms. The main steps are as follows:
The noise-free signal NFS(g, s, t) in Eq. (3) is first estimated for each GCM–SDM chain g–s (top of Fig. 3). The estimation is done for absolute changes in temperatures or relative changes in precipitation (see section 3b).
In a second step, we focus on model uncertainty, associated with the noise-free signal (left-hand side of Fig. 3). The estimated noise-free signal obtained from step 1 is modeled using a two-way ANOVA framework allowing us to separate GCM uncertainty, SDM uncertainty, and residual uncertainty (see section 3c).
In a subsequent step, we focus on model internal variability, corresponding to the noise term η(g, s, r, k, t) in Eq. (3) (right-hand side of Fig. 3). The estimated counterpart of internal variability, say , comes from the estimation of NFS(g, s, t) and is modeled as the sum of two terms: δ(g, s, r, t) related to the large-scale component of internal variability and ε(g, s, r, k, t) related to the small-scale component of internal variability. The estimation of the associated variances is detailed in section 3d.
The total uncertainty is finally derived from these different model uncertainty and internal variability components (see section 3e).
b. Noise-free signal
The NFS for the change variable X (e.g., temperature changes) is estimated for each modeling chain from a trend model fitted to the raw projections Y (e.g., temperature) of the chain as
for absolute changes in temperature and
for relative changes in precipitation, where y(g, s, t) and y(g, s, c) are respectively the trend estimates of the raw projections Y for the future time period t and the reference control period c.
The specific constraints and formulations retained for the estimation of the trends are presented in appendix A. The trend models are composed from a constant value y(g, s, C) over the extended control period [1860, year0] and from a trend function for the transient period [year0 + 1, 2099] where the pivot year0 is the year where the climate is estimated to become nonstationary (Fig. 2).
A linear trend was chosen for precipitation. No satisfying adjustment was obtained with higher-order monomials or polynomials. Higher-order polynomials were observed to overfit time fluctuations of raw data and often lead to no coherent NFSs between runs available for a same GCM–SDM configuration. As reported by Räisänen (2001), a linear trend is expected to overestimate internal variability when the magnitude of this latter is small when compared to nonlinearities in the climate change response. As highlighted in section 4, this is obviously not the case for precipitation as internal variability is very large for this variable. The confidence in a linear trend is also supported by the fact that internal variability estimated for the extended control period is similar to that estimated for the detrended future period (not shown).
A linear trend is conversely not suited for changes in temperature that clearly highlight a progressive increase from the 1950s. A third-order polynomial was found to be optimal, allowing additionally for adjusting for smaller increase rates at the end of the twenty-first century. The pivot year, estimated with a sensitivity analysis, is 1950 for temperature and 1980 for precipitation.
c. Partitioning model uncertainty
It is assumed that the noise-free change response can be partitioned as
where μ(t) is the overall climate response representing the grand ensemble mean of all experiments at projection lead time t, α(g, t) and β(s, t) are respectively the mean deviations of GCM g and SDM s from the grand ensemble mean μ(t), and γ(g, s, t) is the residual.
The parameters α(g, t) and β(s, t) are called the main effects of the models, and μ(t) is the multichain mean. They are estimated using a classical two-way ANOVA framework without interactions (Berrington de González and Cox 2007; Searle 1971). Let us define the following means where the dot symbol indicates averaging over the particular index
With a least squares estimation under constraints and , the parameter estimators in the model of Eq. (6) are given by
residuals are given by
and the total variability (or, in other words, total model uncertainty) of is given by
The three model uncertainty components will hereafter be referred to as GCM uncertainty G, SDM uncertainty S, and residual uncertainty R. The theoretical variances may be estimated by their empirical counterparts, namely
The residuals γ(g, s, t) describe any GCM-dependent deviations of SDMs. These deviations may be purely random. In this case, the deviations γ(g, s, t) obtained for a given GCM–SDM chain depend on the members (run/generations) used to identify the NFS of the chains. GCM-dependent deviations of SDMs, however, may be purely a result of GCM–SDM interactions (e.g., Yip et al. 2011). In this case, the deviation term of Eq. (6) is independent of the members used to identify the NFS of the chains. For a given multimodel ensemble of NFSs, deviations come from the combined effect of randomness and interaction. Estimating the relative importance of each of these components can be obtained using a classical two-way ANOVA with interaction. However, this would require the availability of multiple members of NFSs for each GCM–SDM chain. In the present case, three runs are available for three GCMs. Three members of the NFS can therefore be obtained for each of the corresponding GCM–SDM chains. The interaction term estimated from this subset of chains was found to account for 20%–30% of the variance of the deviation term. The contribution of interactions cannot be estimated for the full RIWER2030 MM2E given that only one run is available for two GCMs; however, it is expected to be of the same order as that obtained from the subset. The uncertainty component associated with the deviation term in the ANOVA will therefore subsequently be referred to as the residual–model interaction uncertainty component (R/MI).
d. Partitioning internal variability
The noise term η(g, s, r, k, t) of Eq. (3) accounts for the internal variability (IV) of the change variable X associated to the selected run r of GCM–SDM chain g–c. It can be also partitioned in the following components:
where δ(g, s, r, t) is the deviation of the rth run of GCM g from the climate response of the GCM–SDM chain g–s and ε(g, s, r, k, t) is the deviation of the kth stochastic generation of SDM s from the rth run of the chain.
The parameter δ(g, s, r, t) is related to the large-scale component of internal variability (LSIV) resulting from the run r of GCM g used to force SDM s and ε(g, s, r, k, t) is related to the small-scale component of internal variability (SSIV) resulting from the stochastic generation k of SDM s when forced by run r of GCM g. It is assumed that δ(g, s, r, t) and ε(g, s, r, k, t) have zero mean.
Expression of LSIV and SSIV components are derived in the following for temperature. Similar expressions were derived for precipitation where relative changes are considered instead of absolute ones. They are presented in appendix B.
1) Small-scale internal variability
For a given GCM–SDM chain g–s, the SSIV of the change variable X can be derived from the SSIV of the raw data Y. When only one run is available for the chain, this latter can be estimated for each projection lead time t as the intermember variance of the 100 stochastic generations for Y. It was found to be roughly constant over the whole 1860–2099 simulation period validating the quasi-ergodic assumption for this component (see, e.g., the bottom graphs of Fig. 1 where the interpercentile distance is roughly independent of time). We could also check that for a given chain there is no correlation between the values of Y in the control and in a future climate. Then it follows that the SSIV of the change variable X = Yt − Yc, reads for any year t:
where k[Y(g, s, r, k, c)] and k[Y(g, s, r, k, t)] denote the intergeneration variance of raw data Y simulated by GCM–SDM chain g–s for the reference control period c and any other period t. To estimate k[Y(g, s, r, k, t)], define the multigeneration mean at time t obtained from the outputs of the rth run of GCM–SDM combination g–s:
Then k[Y(g, s, r, k, t)] is estimated by their empirical counterpart, that is,
It also follows that the SSIV of the change variable X for the GCM–SDM chain g–s can be estimated as the multiperiod mean of 2 times the variance for the raw data Y. The quasi-ergodic assumption also implies that when multiple GCM runs are available for the chain, the SSIV of the change variable X for this chain is the multirun mean of the SSIV estimated for each run of the chain.
Finally, the multimodel mean of this estimate is taken to be the SSIV component for the change variable X for the ensemble of climate experiments. It finally reads
where Ng,r is the number of runs available for GCM g (either 1 or 3 in the present case; see Table 1).
Note that for precipitation, for which relative changes are considered instead of absolute changes, the SSIV expression is a function of projection lead time t (see appendix B). For the sake of notation simplicity, this component will be referred to in the following as SSIV(t) for both temperature and precipitation variables.
2) Large-scale internal variability
The LSIV is the variance component of the noise that remains after having removed the noise due to the stochastic downscaling generation process of the SDMs. It is here estimated for each modeling chain from the time series obtained for each run of the GCM from the multigeneration mean of the raw data. For the rth run of the GCM–SDM chain g–s, Y(g, s, r, •, t) is actually the sum of the noise-free change response for this chain and of the noise component of the run due to LSIV. Assuming that the system is quasi-ergodic over the whole simulation period, the LSIV of the GCM–SDM chain for the raw data Y can next be estimated as the variance over time of the residuals of Y(g, s, r, •, t) from the NFS for Y. Note that for temperature changes, this assumes that LSIV is constant over the entire simulation period. Results obtained for the control and the future period independently show this is a reasonable assumption. When multiple runs are available for a chain, the LSIV for Y is estimated from the variance over time of residuals from all runs. Recall that the NFS correspond in this case to the common trend of these multiple runs.
Following Hingray et al. (2007), we assume that for a given GCM–SDM chain g–s there is no correlation between the values of raw data Y(g, s, r, •, t) obtained in the control and in a future climate. The LSIV of the change variable X is thus also 2 times that of variable Y for this chain.
The multimodel mean of the large-scale internal variability for the change variable X is finally taken to be the large-scale internal variability component. It finally reads
For precipitation, for which relative changes are considered instead of absolute changes, the LSIV expression is also a function of projection lead time t (see appendix B). For both variables, this component will be also referred to as LSIV(t) in the following.
Note finally that by construction δ(g, s, r, t) and ε(g, s, r, k, t) are expected to be uncorrelated and the total uncertainty variance resulting from internal variability is simply the sum of LSIV and SSIV variance components. The absence of correlation can be checked empirically, for each stochastic generation obtained with SDM s when forced by run r of GCM g, with the correlation coefficient obtained between time series δ(g, s, r, t) and ε(g, s, r, k, t) available over the simulation period. For our dataset, the empirical correlation coefficient is in the ]−0.4, 0.4[ interval for more than 90% of the cases.
e. Total uncertainty
Under the assumption that the different variables on the right-hand side of Eqs. (3), (6), and (18) are uncorrelated, the total variance for the change variable X for a given future period t is given by the sum of the different variance components introduced previously. It reads
The grand ensemble mean μ(t) of all experiments at t and the total variance characterize the magnitude and total uncertainty of the change variable X. Note that the total uncertainty corresponding to this total variance is the potential uncertainty that would be obtained from a large number of members for the considered set of GCM–SDM (i.e., large number of runs for each GCM and large number of generations for each GCM–SDM chain). It may be significantly different from the total variance that would be obtained from the MM2E sample of opportunity (Hingray and Saïd 2014, manuscript submitted to Geophys. Res. Lett.).
We further define the fraction of total variance explained by uncertainty source U(t) as
where U(t) is G(t), S(t), R(t), LSIV(t), or SSIV(t).
a. Total uncertainty and sources of contributions
Figure 4 presents, for 20-yr temperature and precipitation, the grand ensemble mean climate change response μ(t) and the limits of the interval , where T(t) is the total uncertainty variance.
For both variables, total uncertainty increases with lead time, especially for projections of temperature changes for which it is multiplied by a ratio of 3.5 from the beginning to the end of the century. It is conversely multiplied by a factor of only 1.4 for changes in precipitation. The main contribution for this increase is that of model uncertainty. Remember that the amplitudes of both IV components are assumed to remain constant over the whole period for temperature changes. For relative changes of precipitation, they are conversely assumed to be a linear function of the grand ensemble mean of the NFSs (see appendix B). In the present case, as the grand ensemble mean is expected to slightly decrease with time, both components of IV highlight therefore a slight decrease over the period.
The relative contributions of model uncertainty and internal variability components to the total uncertainty are given in Fig. 5. They vary significantly with projection lead time and the variable considered. For the first three decades, the contribution of internal variability is highly predominant for both variables. For the first two decades, in particular, the combined contribution of LSIV and SSIV represents nearly all the total uncertainty (more than 80% and 95% respectively for temperature and precipitation). The contribution of internal variability then decreases with lead time as the total uncertainty increases. For temperature, it rapidly drops to less than 10% by the end of the century, then becoming negligible compared to model uncertainty. For precipitation, it is still around 60% of total uncertainty in 2090 and therefore remains the greatest contribution to total uncertainty over the whole simulation period.
For both variables, the LSIV component is greater than the SSIV component. For temperature, the contribution of SSIV is negligible whatever the lead time. It is less than 5% even at the beginning of the century. This clearly reflects the fact that the small-scale variability in temperature is nearly fully explained by the large-scale variability. For precipitation, the contribution of SSIV to cumulated internal variability is around 20% (as a consequence of the quasi-ergodic assumption, this contribution is roughly constant over the simulation period). In 2090, LSIV is 50% of total uncertainty and thus remains the dominant uncertainty source. However, at this time, SSIV still represents more than 10% of total uncertainty (versus 20% for the first decades). For precipitation, SSIV is therefore far from negligible, highlighting, as already discussed, that large-scale variability cannot explain all of the variability observed for this variable (Lafaysse et al. 2014).
For both variables, the contribution of model uncertainty increases over the whole simulation period. The contribution of the error and/or interaction term is nonnegligible; however, it is far less than the contributions of GCM and SDM uncertainty. For temperature, GCM uncertainty is the main contribution for the end of the century. The SDM contribution is however nonnegligible. It reaches up to 20% for the end of the century and even tends to be higher than the GCM contribution for the middle of the century.
For precipitation, the GCM contribution is slightly smaller than the SDM contribution for the whole period. The change in the main effect of each SDM with time was estimated for the next century. Results show that the inter-SDM dispersion from a given SDM with different predictors is as large as that from different SDMs (see Fig. 8 in section 5b). Note also that the contributions of GCM and SDM uncertainty to total uncertainty are much smaller than the contribution of LSIV. They are also smaller than the contribution of SSIV except for the end of the century where the contribution of SDM uncertainty is slightly larger.
b. Significance of changes
Comparing for each projection lead time the grand ensemble mean response to the total uncertainty provides a rough idea of the significance of the estimated changes with respect to any reference level of change. As mentioned previously, the total colored area in graphs of Fig. 4 corresponds to the interval , where T(t) is total uncertainty variance. Assuming for convenience that all possible future climate projections are normally distributed, this interval corresponds to the confidence interval of possible future changes at the 90% confidence level. It allows estimating the significance at a 90% confidence level that all possible future realizations of the climate exceed a given reference level of change. In the following, we will discuss the significance for a nonzero climate change realization. A significant nonzero climate change realization is expected when the zero change value is outside the confidence interval or when the response-to-uncertainty ratio, expressed as , is outside the [−1, 1] interval. The response-to-uncertainty ratio (R/U) is presented for both variables in Fig. 6.
For temperature, a significant nonzero change is predicted from the beginning of the century. The time of emergence of a significant warming, defined here as the first future lead time for which R/U is outside the [−1, 1] interval, is found within the first decade. For precipitation, the nonzero change is not significant. The absolute value of R/U remains much smaller than 1 throughout the entire period, indicating that the sign of change is very uncertain at the 90% confidence level even at the end of the century. Precipitation that may be experienced for a given future period could be therefore higher or lower than that observed for the control period. This is mainly due to the large internal variability value for this variable.
Note that Fig. 4 also allows us to discuss the significance of a nonzero climate change response. The total area covered by the GCM uncertainty (blue), SDM uncertainty (green), and residual–model interaction uncertainty (cyan) corresponds actually to the interval , where M(t) is total model uncertainty expressed as the square root of total model uncertainty variance G(t) + S(t) + R(t). Assuming that the predicted climate change response has a normal distribution, a significant nonzero climate response at a 90% confidence level is obtained when the zero climate change response value is outside this interval (i.e., outside the model uncertainty area). For precipitation, it is here also interesting to note that, even if simulation chains roughly agree on the direction of the change response (a decrease), a nonzero climate change response is still not significant at a 90% confidence level.
c. The potential to reduce uncertainty
It will never be possible to remove uncertainties related to internal variability because they are intrinsic to the Earth system even if uncertainty in the large-scale component can be reduced for coming decades by improving the initialization of climate projections with observations (Smith et al. 2007). On the other hand, GCM and SDM uncertainties might be reduced by a better understanding of climatic and hydrological processes and resulting improvements of numerical models.
Figure 6 also presents the response-to-uncertainty ratio assuming perfect models. This highlights the potential gain that would be achieved when improving both GCMs and SDMs. For temperature, the R/U increase is important especially for the second half of the century. This traduces again the very small impact of internal variability on the significance of changes. This traduces next the large potential to narrow uncertainty by improving simulation models. The high contribution of SDM uncertainty was surprising as we thought initially that a large fraction of uncertainty should be carried by GCMs. It even tends to be higher than that of GCM for the middle of the century. This reflects actually the influence of the different large-scale temperature predictors retained for the SDMs [e.g., the extent of the large-scale predictor domain in the different versions of dsclim (see Table 2) that presents the most different effects]. This critical issue was underestimated in the development of the SDMs. SDM development and selection were based on an extensive evaluation of the model ability to reproduce, in perfect prognosis conditions [where SDMs are forced by National Center for Environmental Prediction (NCEP)–National Center for Atmospheric Research (NCAR) reanalysis; Kalnay et al. 1996], a number of key statistical characteristics of the observed regional climate for the past 50 years. As all SDMs were found to perform much higher for temperature than for the precipitation (Lafaysse 2011; Lafaysse et al. 2014), the main effort in model development was focused on precipitation. A deeper discussion is here obviously out of the scope of this paper but this clearly shows also the interest of the uncertainty analysis that allowed us to identify this likely drawback of some of the SDMs. This clearly shows also that an improved choice of the predictors would allow in this case for reducing the SDM contribution to total uncertainty in future projections. More constraining evaluation tests would have therefore to be carried out in future SDM developments for temperature. Special attention should be paid to the time transferability of the SDMs for this variable (e.g., comparison with temperature changes obtained from GCM outputs directly).
For precipitation, the potential to reduce uncertainty is also nonnegligible; however, model uncertainty accounts for a much smaller part of total uncertainty than for temperature. The R/U increases obtained when assuming perfect models are therefore low to very low, especially for the first four decades. This again reflects the dominant role of internal variability. Even if all models would agree on the ensemble mean precipitation change, this mean change would be very small compared to the internal variability (see the potential R/U in Fig. 6). This means that adapting management practices to the internal variability of precipitation should be a priority as this will likely contribute to successful adaptation to climate change.
5. Discussion and conclusions
a. Methodology and interpretation of results
We introduce a simple two-part approach for the partitioning of model uncertainty and internal variability components in multimember multimodel ensembles of climate experiments. It is based on the quasi-ergodic assumption for the outputs of climate simulations obtained over long control and transient climate period. Model uncertainty components are estimated from the noise-free responses of the different modeling chains using a classical two-way ANOVA framework. Internal variability components, produced by each modeling chain, are estimated from the time series of residuals for the multiple members available for the chain.
Our results for changes in precipitation and temperature are the following:
SDM uncertainty is of same order as that of a GCM. For the present dataset, it was even found to be greater for precipitation during the whole simulation period and for temperature during some decades in the middle of the century.
The contribution of the residual–model interaction term is much smaller than that of SDM and GCM uncertainty, but is not negligible.
Internal variability is the main component of total uncertainty, especially for the first decades. It rapidly decreases to less than 10% at the end of the century for temperature. It roughly always remains the main uncertainty component for precipitation.
The small-scale component of internal variability is negligible for temperature. It is, however, of the same order as the large-scale component for precipitation, whatever the projection lead time.
The response-to-uncertainty ratio at the 90% confidence level is very high for temperature, whatever the lead time. Its absolute value is conversely always much smaller than one for precipitation, indicating that even the sign of precipitation change is uncertain. The time of emergence of change due to global warming is found to be as soon as the first decade for regional temperature but is not expected to be within this century for precipitation.
These findings have important implications. As already suggested by numerous works, neglecting model uncertainty associated with downscaling models is expected to lead to erroneous climate change estimates (e.g., Chen et al. 2011). Depending on the studied variable, the same is likely to apply when internal variability is ignored. For precipitation, the major impact of the large-scale component of internal variability is well established (e.g., Hawkins and Sutton 2011). It is, however, often disregarded. The same applies for the small-scale component of internal variability, associated with downscaling models. Its impact on climate projections, however, is likely to be nonnegligible as highlighted in the present study for precipitation projections. Similar conclusions could be expected with dynamical downscaling models. The presence of internal variability in ensembles of nested regional climate model (RCM) simulations is actually widely acknowledged in the community working on dynamical downscaling (Nikiéma and Laprise 2013). Alexandru et al. (2007) found, for instance, that seasonal statistics of precipitation fields simulated with the Canadian Regional Climate Model are likely to be poorly estimated from a single model simulation and that a robust estimation, depending on season and region, could require a minimum number of 10 members. For climate projections, the internal variability associated to the RCM could be locally larger than the mean climate change response, as shown for instance by Braun et al. (2012). Part of this internal variability is expected to be reduced with a closer forcing of the RCM toward the driving circulation but part of it must also correspond to the freedom the system has of creating different local weather from different weather patterns. If the real contribution of this internal variability component seems to remain unclear at this time, its impact on climate projection is therefore worth more attention.
Multimodel (GCMs and downscaling models) experiments became a standard in climate change impact studies. Another standard should be to rely on multimember experiments. Impact studies based on single members of SDMs or RCMs experiments (or small ensembles) are likely to be not more relevant than those based on single runs of available GCMs (or small ensembles). When they are intended to provide information for climate change adaptation, they may lead to poor decisions. A relevant strategy would be, in the present case, to adapt to internal variability of precipitation.
b. Partitioning uncertainty components with the quasi-ergodic ANOVA framework
A strong advantage of the quasi-ergodic ANOVA (QE-ANOVA) approach presented in the present work is that all periods and experiments available in the ensemble of climate experiments opportunity can be valorized. A critical and classical problem in similar studies is actually that the number of available runs differs from one GCM to the other. As an ANOVA in such an unbalanced case is not easily tractable (Sansom et al. 2013), a unique number of runs is next usually considered for all GCMs (e.g., Hawkins and Sutton 2009; Hingray et al. 2007; Yip et al. 2011). This simplifies the ANOVA framework and allows us to treat all GCM–SDM chains equally. This leads us, however, to disregard a large number of available runs. Additionally, Deser et al. (2012) and Northrop (2013) pointed out the fact that results can depend strongly on the particular dataset chosen, so relying on a single dataset may be misleading. A major challenge is to have therefore the possibility to use all available experiments for the impact study under consideration (Sansom et al. 2013).
In the present framework, the model uncertainty analysis is carried out on the climate change responses estimated respectively for each GCM–SDM chain. The climate change response is unique for a given chain whatever the number of runs available to estimate it. In the classical ANOVA part of the method, all available runs can therefore be accounted for and all GCMs can be treated equally, whatever the number of runs accounted for. The residuals from the climate response estimated for each modeling chain allows next for estimating the different components of internal variability.
In the present work, a large number of SDM generations was made available for each run of each GCM–SDM chain. This allowed for an easy separation of both large- and small-scale components of internal variability. Such a dataset is likely to be not always available. In such a case, the partition between both components would not be possible. The method would, however, allow estimating the total internal variability resulting from both components using an equation similar to Eq. (23).
In all cases, all available data can again be used to estimate the large-scale component of internal variability. For a given modeling chain, this latter can be actually estimated from the residuals obtained for the entire control and transient simulation period from all members available for the chain.
c. Robustness of uncertainty components estimates
The more members we can rely on, the more robust the estimates of uncertainty components are obviously expected to be. The climate response for a given modeling chain, however, can be already identified from only one member. In this case, when a long control and transient period is available (e.g., 240 yr in the present study), the quasi-ergodic assumption makes it already possible to estimate the large-scale internal variability component from a quite large sample of residuals (seven in our case).
In this work, results were presented using all experiments of the RIWER2030 MM2E, and especially the three runs available for three of the five considered GCMs (Table 1). The results are roughly the same when only one run is used for each GCM. This is illustrated in Fig. 7 by the roughly unchanged partitioning of total variance obtained when only the first, second, or third run is used. The main difference is the slightly smaller contribution of the large-scale internal variability component obtained for precipitation. The robustness of the method is also illustrated by the fairly unchanged patterns of main model effects shown in Fig. 8 for such configurations. For both variables, the main effect of SDMs is similar whatever the run or multirun experiment considered. The main effect of the GCMs is somewhat more sensitive to the experiment but the general pattern of main effects is the same.
The robustness of uncertainty component estimates obviously depends on the possibility to achieve a robust estimate of each modeling chain climate response. A potentially critical problem is therefore the potential errors associated to the identification of the climate response of each modeling chain. A first difficulty concerns the choice of the type of trend used for the NFS estimation. For convenience, analytical trend functions were chosen in the present work and a same analytical model (with chain-dependent parameters) was retained for all GCM–SDM chains. For some chains, better adjustments could have been probably obtained with higher-order polynomials or with trend functions estimated with nonparametric methods. An overly adaptable trend function, however, is expected to underestimate the internal variability of the chain and to be too sensible to the sample of data used for its identification (i.e., to the member retained for the estimation).
Even when a simple trend analytical function is retained, the dependence of the trend to the data sample may be high and it could be therefore difficult to identify the true climate change response of the modeling chain. This problem is of course expected as soon as internal variability is high compared to the climate response of the chain. This is likely to explain differences in the results presented here, for precipitation especially (e.g., the rather large dispersion of the main GCM effects between one-run configurations or the smaller contribution of large-scale internal variability component to total uncertainty for one-run configurations).
From a 40-member climate experiment, Deser et al. (2012) found that the minimum number of members required for an accurate estimation of the climate change response (linear least squares trends fit to the period 2005–60) can be quite large depending on the considered variable. Because of the relative amplitudes of the climate change response and natural variability, they highlighted that one member is needed to detect a significant (at the 95% confidence level) warming in the 2050s decade compared to the 2010s at nearly all locations, compared to approximately 3–6 (>15) ensemble members for tropical and high-latitude (midlatitude) precipitation.
For a number of impact studies, the relevance and significance of the climate responses identified from available MM2E of climate projections are therefore expected to be rather problematic. For this critical identification issue, robust methods and clear consensual guidelines would be highly welcomed from the climate research community. Further methodological developments are also required to complement the QE-ANOVA framework. They should allow for a thorough analysis and quantification of the different error sources liable to impact uncertainty estimates.
The QE-ANOVA framework can be easily applied for a large number of different datasets. The main requirement is that long-term simulations are available to allow for a robust estimation of the climate response of each modeling chain and of its internal variability. In all cases, the NFS would allow for the estimation of all model uncertainty components and the residuals from this NFS would allow for the estimation of the different internal variability components.
A straightforward application is for instance to quantify the uncertainty in hydrometeorological projections obtained by simulation with a given hydrological model from a dataset similar to the one used in the present work. An application is presented by Lafaysse et al. (2014) for hydrological projections obtained for the upper Durance River basin from the present MM2E of climate projections. The same framework could be also easily extended to hydrological projections from multiple hydrological models (HMs). The NFSs resulting from all GCM–SDM–HM chains would make it possible to partition the main model effects. This could be done using a classical three-way ANOVA. As the hydrological behavior of catchments is mainly deterministic, as reflected by the large majority of hydrological models used for climate change impact studies, the HM is not expected to introduce any additional component to the internal variability of the chain. Internal variability could be partitioned, in the same way as presented here, into its large- and small-scale components related mainly to GCMs and SDMs respectively.
As long as climate experiments would be available for a number of different emission scenarios, another straightforward extension of this framework would be also to additionally estimate scenario uncertainty. Each emission scenario would actually influence the NFSs of the modeling chains. Scenario uncertainty would be thus considered as an additional “model uncertainty factor” to be estimated from the classical ANOVA part of the framework.
Note finally also that the QE-ANOVA framework presented here does not necessarily require all the model components of the simulation chain. For instance, it does not require that GCM experiments are downscaled with SDM or RCMs. It can be applied directly on GCM outputs if only experiments from different GCM are available. The partition would here only concern the uncertainty resulting from the GCM and that due to large-scale internal variability. As already mentioned, the framework does also not require all the members used in the present work (e.g., different runs of the same GCM, a large number of stochastic generations from the SDM). It could be also easily applied on the outputs of GCM–SDM experiments with a single SDM generation for each GCM–SDM chain. In such a case, the different components of internal variability could however not be partitioned and the estimate of internal variability would contain both components.
This work was initiated during the Regional Climate, Water, Energy Resources and Uncertainties from 1960 to 2030 (RIWER2030) research project (http://www.lthe.fr/RIWER2030/) funded by Electricité de France (EDF), the Centre National de la Recherche Scientifique (CNRS), and the French National Research Agency (ANR). It was finalized within the R2D2-2050 research project (https://r2d2-2050.cemagref.fr/) founded by the Ministère de l’Ecologie, du Développement Durable et de l’Energie. The author wishes to thanks Matthieu Lafaysse (Méteo-France, Grenoble), Joël Gaihlard (EDF, Grenoble), and Abdelkader Mezghani (MetNo, Oslo) for providing climate projections. We thank also Jerome Buzzi (Orsay, Paris), Douglas David Baptista De Souza (LTHE, Grenoble), and Jean-Philippe Vidal (Irstea, Lyon) for interesting discussions on this work as well as Michel Slivitzky for his enthusiastic appreciation of a former manuscript version. We finally thank the three anonymous reviewers for constructive suggestions.
Different trend models were used in previous works to fit the time evolution of future projections from climate experiments. To estimate internal variability in 17 experiments from phase 2 of the Coupled Model Intercomparison Project (CMIP2), Räisänen (2001) assumed that the climate change response for the 2000–80 period was a linear function of time, like the greenhouse run radiative forcing that was proportional to the logarithm of atmospheric CO2. Hawkins and Sutton (2009) assumed the noise-free signal to be a fourth-order polynomial of time over the 1950–2099 period.
In the present work, an analytical trend model is also used for convenience. It applies on the raw temperature or precipitation projections Y. For each GCM–SDM combination, the trend model is estimated for the whole 1860–2099 period with following constraints:
It is composed of a constant value y(g, s, C) over the extended control period [1860, year0] and of a linear or polynomial trend for the transient period [year0 + 1, 2099], where y(g, s, C) is the mean value of the variable estimated for the extended control period and constitutes the starting value of the adjustment for the transient period.
The pivot year and the polynomial degree are assumed to be the same for all GCM–SDM chains. They vary from one variable to the other. The coefficients of the linear or polynomial trend are the same for all members of a given GCM–SDM chain. They vary from one chain to the other.
The fit is made using ordinary least squares. In case of a polynomial trend, the fit insure that the derivative of the polynomial is zero for the pivot year (see equations below).
A basic assumption there is thus that the climate response of a given GCM–SDM configuration varies progressively over the 1860–2100 period with a continuous derivative for the pivot year especially. Fitting a trend without these constraints potentially leads to a noise-free climate change that overfits the raw series around the 2000 year especially, leading in turn to overestimating the model uncertainty components and underestimating the internal variability components.
Note that y(g, s, C) is not necessarily equal to the raw projection Y obtained for the reference control period c from which future changes are estimated. On the one hand, the pivot year does actually not necessarily correspond to the middle year of the reference control period c. In turn the trend estimate y(g, s, c) of the raw projection Y for the reference control period c is a priori different from the constant value y(g, s, C) estimated over the extended control period C (see Fig. 2). On the other hand, the raw value Y(c) is also expected to be significantly different from the trend estimate y(g, s, c) as soon as internal variability is nonnegligible (e.g., for precipitation in the present case). The trend function identified using y(g, s, C) as starting value for the transient climate period is therefore expected to be much more robust and thus relevant than the one that could have been obtained using the sample mean of the reference control period c instead.
The analytical expressions of the coefficients of trend function have been derived to meet the above mentioned constraints. Different trend models were tested for the raw variables considered in this work: a monomial trend of degree n (1 ≤ n ≤ 3) and a polynomial trend of degree 3. In the case of a monomial trend, the trend has the following expression:
where x0 is the constant value estimated for the extended control, τ = t − t0, and the least squares estimate of coefficient a is
where Yi is the raw projection of the studied variable for future period centered at year ti and T is the final future projection period (2080–99 here). This expression holds for every strictly positive real number. For n > 1, it insures that the derivative is zero when t = t0.
In case of a polynomial of degree 3, the trend expression reads
where the least squares estimation of both coefficients with constraints mentioned above are
The choice of the trend model was made studying the homoscedasticity of the residuals and the relevance of the trend shape over the whole simulation period. The sensitivity analysis considered different types of models and different pivot years (year0).
Internal Variability Components for Relative Changes
Uncertainty components in the case of relative changes are derived in this appendix following appendix A of Hingray et al. (2007). Let consider the simulations outputs over the 1860–2099 period of modeling chain m defined by the rth run of GCM–SDM combination g–c. Assuming that at any time the noise-free estimate of Y is a good approximation of the expected value of Y, approximations to the variance of the change variable X for this chain m can be obtained with the delta method (Stuart and Ord 1987, sections 10.5 and 10.6):
where denote the estimated total internal variability of the change variable X for modeling chain m and where and are the estimated variances of the raw data Y for the future period t and reference control period c, respectively.
Assuming that the large-scale and small-scale internal variability components are noncorrelated, both terms in the parentheses of the right-hand side of this equation can be partitioned into their large- and small-scale components.
An estimate of the small-scale variability component therefore reads as follows:
Both two terms of Eq. (B2) are equivalent to a coefficient of variation of Y with respect to the intergeneration variance. This coefficient of variation was assumed to be roughly constant for a given chain over the whole 1860–2099 simulation period. Its multiperiod mean was used for estimating the SSIV of the relative change variable X for the modeling chain m when only one run is available. It reads
When multiple runs are available for a given GCM–SDM chain, the multirun mean of these estimates was retained as SSIV for this chain. In both cases, the SSIV for chain m is a function of time via the term y(g, s, t) in the equation. The SSIV component of the relative change variable for the projection ensemble was next estimated for each future lead time t as the multimodel mean of these variances, from the following equation:
Note that Eq. (B4) has a form similar to Eq. (22). The raw variable is replaced by the raw variable normalized by its expected value and a multiplicative coefficient accounts for the relative change in expected values from the control to the future period. Note that for t ≤ c, the SSIV is 2 times the multimodel mean of the coefficient variation of Y with respect to the intergeneration variance.
The large-scale internal variability component has the same expression as that of SSIV in Eq. (B2) but, because of the limited number of runs available, the interrun variance cannot be estimated. Similarly to the assumption made for changes in temperature, we assume that the LSIV for Y with respect to the interrun dispersion is, in term of coefficient of variation, constant over the whole simulation period. It follows that for any time t
When multiple runs are available for a chain, this variance component is estimated accounted for data from all runs. The LSIV component of the relative change variable X was next estimated from the multimodel mean of the interperiod interrun variance of Y(g, s, r, •, t) as
Current affiliation: Agence Technique de l’Information sur l’Hospitalisation, Lyon, France.