Global sensitivity analysis of the climate – vegetation system to astronomical forcing : an emulator-based approach

A global sensitivity analysis is performed to describe the effects of astronomical forcing on the climate–vegetation system simulated by the model of intermediate complexity LOVECLIM in interglacial conditions. The methodology relies on the estimation of sensitivity measures, using a Gaussian process emulator as a fast surrogate of the climate model, calibrated on a set of well-chosen experiments. The outputs considered are the annual mean temperature and precipitation and the growing degree days (GDD). The experiments were run on two distinct land surface schemes to estimate the importance of vegetation feedbacks on climate variance. This analysis provides a spatial description of the variance due to the factors and their combinations, in the form of “fingerprints” obtained from the covariance indices. The results are broadly consistent with the current understanding of Earth’s climate response to the astronomical forcing. In particular, precession and obliquity are found to contribute in LOVECLIM equally to GDD in the Northern Hemisphere, and the effect of obliquity on the response of Southern Hemisphere temperature dominates precession effects. Precession dominates precipitation changes in subtropical areas. Compared to standard approaches based on a small number of simulations, the methodology presented here allows us to identify more systematically regions susceptible to experiencing rapid climate change in response to the smooth astronomical forcing change. In particular, we find that using interactive vegetation significantly enhances the expected rates of climate change, specifically in the Sahel (up to 50 % precipitation change in 1000 years) and in the Canadian Arctic region (up to 3 in 1000 years). None of the tested astronomical configurations were found to induce multiple steady states, but, at low obliquity, we observed the development of an oscillatory pattern that has already been reported in LOVECLIM. Although the mathematics of the analysis are fairly straightforward, the emulation approach still requires considerable care in its implementation. We discuss the effect of the choice of length scales and the type of emulator, and estimate uncertainties associated with specific computational aspects, to conclude that the principal component emulator is a good option for this kind of application. Published by Copernicus Publications on behalf of the European Geosciences Union. 206 N. Bounceur et al.: Climate–vegetation response to astronomical forcing


Introduction
The seasonal and spatial distribution of incoming solar radiation (insolation) at the top of the atmosphere is determined by three astronomical parameters: Earth's eccentricity, the longitude of the perihelion, and Earth's obliquity.The variation in these parameters causes sufficient changes in insolation to significantly affect climate, such as the distribution of surface temperature, vegetation cover, monsoon rainfall, Arctic sea ice, etc.These changes can be simulated and studied by means of experiments with global climate models.One classical approach consists in identifying two epochs in the past for which sufficient data are available, running the climate model (with the implicit assumption that simulated climate is quasi-stationary with respect to the astronomical forcing), and then comparing the two resulting simulated climates.This is the approach followed, for example, by the Paleoclimate Modelling Intercomparison Project (Braconnot et al., 2007).
The difference between the two epochs may then be further decomposed into contributions of several factors.Suppose we want to compare the beginning and the end of the last interglacial.These two periods, distant by 11 000 years, are characterised by significant climatic differences (e.g.Sanchez-Goñi et al., 1999).The difference in insolation forcing is caused by a change in the position of the perihelion on Earth's orbit (this is the precession process) and a decrease in obliquity.Based on a series of transient experiments, Crucifix and Loutre (2002) suggested that the sum of individual effects of precession and obliquity during the Eemian are less than their combined effects.In the model experiments discussed in that article, this is caused by feedbacks associated with the responses of vegetation and sea ice.Tuenter et al. (2003) investigated the difference between individual and combined effects more generically, without reference to specific epochs of the Quaternary.They performed a set of seven experiments with different combinations of precession and obliquity spanning the range of values reached during the Quaternary.They found that "the amplitude of the [North African] precipitation response to obliquity depends on precession, while the precipitation response to precession is independent of obliquity".The fact that the response to obliquity depends on precession may be interpreted as a second-order effect.
One difficulty appears when one points to the possibility that second-order effects may be significant during certain epochs, or at certain critical locations only.A local analysis, focused on specific periods, such as the study by Crucifix and Loutre (2002), may thus by chance over-or underemphasise second-order effects.The approach of Tuenter et al. (2003) is more global in the sense that it does not refer to a specific epoch, but the coverage of astronomical forcing configurations may be too coarse to a detect second-order effects reliably enough.
To see how the problem may be formulated statistically, let us use vector x to express the astronomical forcing.Then, let us denote the output of a model run at configuration x by f (x).The vector x varied in the past, and we can estimate its distribution by reference to existing astronomical solutions.If we assume that the climate model represents reality and that the climate is quasi-stationary with respect to the astronomical forcing, then f (x) reflects the past evolution of climate.We can in particular enquire about the variance in f (x) caused by the distribution of x, decompose this variance into contributions from individual factors (precession and obliquity in this example), and then check for secondorder effects.This is the principle of global sensitivity analysis (Homma and Saltelli, 1996;Saltelli et al., 2004), for which there exist methods specifically tailored for complex numerical simulators.
We therefore propose to carry out a global sensitivity analysis of the climate model of intermediate complexity LOVE-CLIM (Goosse et al., 2010) to the astronomical forcing.We provide geographic distributions of the contributions of obliquity and precession on precipitation and temperature.We also attempt to detect regions where fast climate responses may occur in response to the slow changes in astronomical forcing.The objectives of this work are twofold.The first is climatic: we would like to determine the respective roles of astronomical factors on interglacial climate change and contribute to the discussion on the mechanisms of glacial inception and interglacial duration with a focus on climatevegetation interactions, adding to the discussion in de Noblet et al. (1996), Claussen et al. (1999), Crucifix and Loutre (2002), Kageyama et al. (2004), and Meissner et al. (2003).Our second objective is methodological: while global sensitivity analysis is well established in statistics, it has only recently been applied to climate problems.Lee et al. (2011), for example, performs a global sensitivity analysis of a global atmospheric model associated with a complex aerosol model in order to decompose model output uncertainty into contributions from eight uncertain parameters (see also Lee et al., 2013;Carslaw et al., 2013).Although our scientific objective is different (we want to decompose forced climate variances induced by variances in forcing factors), the methodological approach that will be followed here is broadly similar to that of Lee et al. (2011) and Oakley and O'Hagan (2004): (a) choose the experimental design to efficiently fill the input space (here, the space of astronomical forcings); (b) run LOVECLIM in these experiments; (c) train and validate an emulator, that is, a stochastic statistical model used to predict the function f (x) at any input point, based on the output of the experiments actually run; and (d) use the emulator to estimate sensitivity indices.Emulation has been increasingly used in climate science in recent years as a tool to explore input spaces with the aim of calibrating the model on observations and estimating climate sensitivity (Rougier et al., 2009;Holden et al., 2010;Schmittner et al., 2011).We explore the potential of this methodology for our specific application, with a discussion of its possible advantages, challenges, and drawbacks compared to more classical approaches.

Choice of input factors
In palaeoclimatology, it is common to refer to the time of the year using the true solar longitude (λ), that is, the geocentric angle between the vernal equinox and the position occupied by the Earth at any point during the year.For example, the June solstice corresponds to λ = 90 • , the September equinox to λ = 180 • , etc.For the purpose of computing insolation at a given time of year, we need the true solar longitude at perihelion, that is, the true solar longitude corresponding to the shortest Earth-Sun distance.This quantity is denoted .The shape of the Earth's orbit is elliptic and characterised by eccentricity e.Finally, the angle between the ecliptic and the equator is called the obliquity and denoted ε.It may then be shown that the secular evolution of the top-ofthe-atmosphere incoming solar radiation at any latitude and any true solar longitude is reasonably well approximated by a linear combination of i 1 = e sin , i 2 = e cos , and i 3 = ε (Loutre, 1993).The quantity i 1 is often referred to as the climatic precession parameter.As i 1 , i 2 , and i 3 are not correlated, the three inputs can be viewed as a canonical form of the astronomical forcing parameters.In particular, their signature on the season-latitude distribution of incoming solar radiation is characteristic: i 1 and i 2 control the Earth-Sun distance at any true solar longitude with very little effect on annual mean insolation1 , and i 3 controls the seasonal contrast as well as the annual mean insolation.
The time evolution of the astronomical parameters over the Pleistocene is well known (Berger, 1978b;Laskar et al., 2004).Note that, 99 % of the time, eccentricity e < 0.05 and the inner 99th percentile of ε is 22.3 to 24.3 • , these boundaries differing by less than 0.1 • if the Berger (1978b) or Laskar et al. (2004) solution is used as the reference.

Experiment design
As specified above, the climate model we use is the oceanatmosphere-vegetation model of intermediate complexity, "LOVECLIM" (Goosse et al., 2010).The choice of inputs for the experiments are encoded in an input factor matrix denoted X, which has three columns corresponding to the inputs i 1,2,3 , and which has as many rows as ensemble members (27 experiments in our case).There is a rich literature on experiment design sampling techniques, and, for example, the monograph by Santner et al. (2003) is specifically dedicated to computer experiments.As a general rule, given a fixed number of experiments to be run, one objective is to maximise the information that can be inferred from the experiment set.If we use a Gaussian process (GP) emulator (as here), this criterion corresponds to minimising the posterior variance in the GP predictions (Sacks et al., 1989).Another objective is to minimise the bias in the quantities to be estimated.Finding an adequate design is thus an optimisation problem, and generally not a straightforward one, because the GP hyper-parameters are not known a priori, and so the design must trade off the need to learn these parameters with the need to minimise the prediction variance.In practice, an effective approach consists in using Latin hypercube designs (McKay et al., 1979;Morris and Mitchell, 1995;Urban and Fricker, 2010) with specific properties, such as maxi-min properties (where we maximise the minimum Euclidean distance between any pair of design points) or orthogonality (in this case, maximising X X).Hybrid designs combine several such properties and can be theoretically justified as good solutions for Gaussian process emulation (Joseph and Hung, 2008).Out-of-the-box packages can be used to produce Latin hypercube designs (e.g. the lhs package; Carnell, 2012), but they are unsuited to our needs as we face an additional constraint: we want to avoid wasting computational effort sampling eccentricities > 0.05, which translates into a constraint on the quantity i 2 1 + i 2 2 .(Vernon et al., 2010) and (Draguljić et al., 2012) provide examples of designs for non-rectangular input spaces.The algorithm used here (Appendix A) is similar to that of (Vernon et al., 2010): it satisfies the constraint on eccentricity while aiming at near orthogonality and max-min properties.The choice of n = 27 experiments arises from the general recommendation of (Loeppky et al., 2009) to perform 10 experiments per input dimension.We used 27 members (and not 30) in order to compare with a factorial design of 3 3 members (Bounceur, 2015).
The resulting design is shown in Fig. 1.It is executed three times.Two ensembles use the standard version of LOVE-CLIM with the VECODE vegetation model (Brovkin et al., 1997) but with two distinct sets of initial conditions.The first ensemble uses the pre-industrial conditions provided by default in the LOVECLIM model package.The second uses the final state of member 2 of the first ensemble.This is a so-called warm orbit (high obliquity, high eccentricity, and 90 • ) that produces extensive vegetation cover in the Northern Hemisphere and in the Sahara region.Every experiment is run for 2000 years, and the data used for the following analysis are obtained by averaging the last 500 years.
The purpose of using two distinct initial conditions is to detect potential co-existence of multiple steady-state solutions.Brovkin et al. (1998), for example, found two stable equilibria for certain orbital configurations when using the same vegetation model as here, but a different atmosphere-ocean system.If such multiple states were to co-exist, then the simulator output would need to be emulated using a multi-modal process rather than a (unimodal) Gaussian process as used here.Our initial intention was to build two distinct emulators: one using the warm orbit as initial conditions, and one using the cold orbit.In practice, the two emulators ended up being redundant (Sect.3.1), and thus a third of our computational budget was thus spent on replicating essentially identical runs.However, this could not have been easily foreseen, and we feel this was worth the cost so as to avoid the risk of missing the existence of multiple steady states.The third ensemble uses the same initial conditions as the first but considers the original "ECBILT" surface scheme with fixed vegetation (Opsteegh et al., 1998) and not the VECODE model.We chose to build a third independent emulator in this case, but we note that "multi-level emulation" approaches, which associate similar but different simulators (Cumming and Goldstein, 2009), could potentially have increased our predictive accuracy by allowing us to make better use of the limited computational resource.We leave this option as a possible subject for further investigation.

Global variance measures
In the present application the inputs (the astronomical forcing) are known, and we want to estimate the simulator output variance induced by separate and/or combined variations of the different forcing factors.Recall that we defined f (x) to be the simulator response at input vector x, which contains the values of the three astronomical forcing terms.Let ρ(x) be the time-wise occupation density of the input space during the Pleistocene, which can be estimated from standard astronomical solutions (Berger, 1978b;Laskar et al., 2004).The total variance in the output associated with the variations of factors may then be expressed as where the variance operator means that we are sampling (varying) x following the distribution ρ(x).As f (x) is a vector, V is a matrix with elements giving covariances between any output pair.Suppose we are now interested in the variance in the simulator output caused by the variation of a subset of the input factors only.Let x i denote a subset of the components in x, and let x −i be the remaining components.For example, if i is obliquity, −i will be the indices associated with eccentricity and precession.
We define the following quantities (e.g.Saltelli et al., 2004): the main effect, η(x i ), is the expected output conditional on the value of x i , i.e. η(x i ) = E[f (x)|x i ]; the main effect variance is the variance in the main effect with respect to x i : On the other hand, the output variance associated with factors x i varying while the factors x −i are fixed is denoted Var(f (x)|x i ).This is a function of the value of the fixed factors x −i .If we further average this variance over possible values of x −i , that is, we take the expectation of this quantity with respect to x −i , we obtain the total effect variance associated with factor x i , denoted T i : (2) The law of total variance implies V i := V − T −i , so that we can see that the main effect variance may be interpreted as the expected loss of output variance resulting from fixing (knowing) the value of x i .The standardised quantity V i /V is commonly referred to as the main effect index, and T i /V as the total effect index (cf.Homma andSaltelli, 1996, andSaltelli et al., 2004).
To estimate these sensitivity indices, we extend the framework established by Oakley and O'Hagan (2004) for onedimensional scalar simulator outputs to multi-dimensional vector outputs.Let X , X i , and X −i be the domains of the input factors x, x i , and x −i , respectively.The total variance can be explicitly written as (3) The main effect is computed as follows: Earth Syst.Dynam., 6, 205-224, 2015 www.earth-syst-dynam.net/6/205/2015/where ρ(x −i |x i ) is the conditional density of x −i given x i .
Because our knowledge of f (x) is limited to the ensemble of model runs, we are uncertain about the value of η(x i ) for all values of x i .The approach described in Oakley and O'Hagan (2004) is to build an emulator of f (x) that can predict its value for any input configuration x.They use a Gaussian process (GP) model for the emulator, with mean function m(x) and covariance function between f (x) and f (x * ), which we will denote (x, x * ).Because the emulator is a Gaussian process, the main effect η(x i ), which is a linear transformation of f (x), is also a Gaussian process, and has mean and variance functions It follows that where E f denotes expectation with respect to the emulator, and The expectation of the total variance V of the simulator output is From now on, it will be implicit that the V, V i , and T i are estimated with the emulator and the symbol E f will be dropped.Specifically, we refer to T {e } as the total effect variance associated with precession, and T ε as that associated with obliquity.Finally, we refer to the quantity V−T {e } −T ε = T ε −T ε as the synergy term.More generally, T i − V i is a measure of how much the factor i is involved in interactions with any other input variable.The word "synergy" is commonly used in the climate literature to express the difference between the model response to different factors varied together, and the sum of the responses to factors individually.Unlike the global sensitivity indices discussed here, synergy terms are classically estimated on the basis of a reference experiment, and comparing changes one at a time with combined changes in the different factors, following a method called factor separation analysis (Stein and Alpert, 1993;Alpert and Sholokhman, 2011;see, for example, Crucifix and Loutre, 2002;Ganopolski et al., 1998;Claussen, 2009;Braconnot et al., 1999;Berger, 1999;Henrot et al., 2010;Wohlfahrt et al., 2004).

Motivation
As the outputs of our simulator are spatially resolved climate quantities, we need to build an emulator capable of modelling multivariate outputs.A simple pragmatic solution is to train independent emulators for each grid point, which is done by Lee et al. (2011).There are, however, several other possibilities for multivariate emulation (see, for example, Rougier, 2008).The main challenge is defining a covariance function for generating the covariance matrix (x, x * ) in order to produce a valid covariance function for the Gaussian process.See Fricker et al. (2013) for discussion.Here, we propose the principal component (PC) emulator (Higdon et al., 2008;Wilkinson, 2010) as a cost-effective and statistically reasonable alternative to point-wise emulation.The advantages of this choice are commented on in Sect.2.6.The derivation of sensitivity indices with PC emulation is largely based on published material, which we elaborate on here in order to introduce the notation leading to the final Eq. ( 18).To our knowledge, this has not been given elsewhere.The reader mainly interested in the emulator performance and climatological analysis may immediately jump to Sect. 3.

Principal component decomposition
Let Y denote the matrix in which each column represents the output of one experiment, i.e.Y = [f (x (1) ), . .., f (x (n) )], where x (j ) is the input of experiment j .For example, if we want to emulate annual mean surface temperature, f is a vector of p = 2048 components (the number of grid points).Denote the p vector of row averages of Y by Y .We now define the centred matrix Y = Y − Y 1 p , with 1 p a vector of length p with all components equal to 1, and consider the singular value decomposition (SVD) Y = UDV , where D is a diagonal matrix and U and V are square orthonormal matrices.The columns of U represent the basis vectors, {u k }, and the projection coefficients are given by VD.For the j th experiment, the coefficient for the kth basis vector is a k (x (j ) ) = V j,k D kk .Wilkinson (2010) keeps the first eigenvectors only (ordered by decreasing eigenvalues).This is sometimes known as "hard-thresholding regularisation" (Silverman, 1996), and results in a reduced-order model: f (x) ≈ k=1 a k (x)u k .
Thus, for any input x, the emulator output prediction will lie in the space spanned by the {u 1 . ..u }.The error associated with hard-thresholding is accounted for in the covariance of the estimator (see Sect. 2.5).

Emulation of PC scores
We then need to predict the a k (x) based on the experiment output a k (x i ).This is done by considering Gaussian process (GP) models (Rasmussen and Williams, 2005).We suppose the prior mean of each GP is h(x) β k , where h(x) is a vector of q known regression functions and β k is the vector of corresponding regression coefficients.This prior mean is conditional on parameters β k , which will then need to be estimated.For the present application, h(x) is simply defined as (1, x 1 , x 2 , x 3 ) (linear regressors), so that h(x) We then define the following: -H, the matrix which has row i equal to the regressors h(x i ); , the vector of PC scores we wish to emulate; )), the correlation function for a k (•).This is typically a monotonically decreasing function of the distance between the two points.We let A k be the n×n Gram matrix, with If each a k (•) is modelled as a Gaussian process, then we have a Gaussian prior predictive distribution for y conditional on The interpretation of this model is that the simulator response is the sum of a mean response function, expressed as the linear combination of regressions, and a stochastic component (a zero-mean Gaussian process) that absorbs deviations from the mean response.
In order to estimate σ 2 k and β k , we assume a conjugate non-informative prior distribution, i.e. π (β k , σ 2 k ) ∝ σ −2 k (Berger et al., 2001).This allows σ 2 k and β k to be marginalised out of the analysis, resulting in a posterior distribution of the simulator output that is a Student t distribution with n − q degrees of freedom, with mean and variance functions where We primarily use the squared exponential covariance function with a nugget term for each c k (as discussed at length in Andrianakis and Challenor, 2012): where k is a scaling matrix chosen to be diagonal and ν k is a "nugget", which we will discuss shortly.The diagonal elements of k are commonly called the length scales.The parameters k and ν k are then found by optimising a penalised likelihood L (pen) k of the Gaussian process associated with the principal component k, as a function of k and ν k .The expression of L k is given in (Andrianakis and Challenor, 2012): The role of the penalty is to guarantee smaller Gaussian process variances than would be obtained by least-squares regression.
The nugget term, ν k I i=j , was originally introduced to account for measurement errors in geospatial data analysis (Cressie, 1993).In emulators of deterministic systems, the nugget may also be justified as a regularisation ansatz to avoid poor matrix conditioning (Pepelychev, 2010) or as a way to account for the mis-specification in the correlation function (Gramacy and Lee, 2012).In climate model applications, the nugget may also be justified as a way to account for "internal variability".Indeed, the chaotic dynamics of the simulator are such that a particular climate average over a given time window can be viewed as a stochastic quantity, even though the simulator is deterministic.In the climate modelling parlance, the effect is referred to as uncertainty associated with the internal simulator variability.For example, in Araya-Melo et al. (2015), we found that our estimate of the nugget is consistent with the assumption that this term represents the uncertainty due to simulator variability.This is also the interpretation adopted by Williamson et al. (2014) (both studies use the climate model HadCM3).In the present application, we use rather long climate averages (500 years) and we anticipate that internal variability will be of the same order as the error associated with the truncation of the principal components.It may thus not be appropriate to interpret ν k as an indicator of internal variability, and we therefore chose not to do so.
Finally, two options will be considered to estimate the different k and ν k associated with the different principal components: (1) maximise L k independently for each k or (2) use the same parameters for all k, i.e. k = , and ν k = ν and optimise k=1 L k .Although option (2) does not maximise the likelihood of the emulator taken as a whole, it presents computational benefits in the context of global sensitivity analysis as we show next.Whichever option is used, the k and ν k are, once chosen, considered as known.In particular, the means and variances given in Eqs. ( 12) and ( 13) are conditional on the values of these parameters.

Recombination of PC scores
The emulator posterior distribution for predictions of LOVE-CLIM's outputs f (x) follows a distribution with mean and covariances functions given by .
The covariance matrix of the emulator for LOVECLIM is thus of dimension p×p and provides information on the joint uncertainty of any two simulator outputs.
Variance indices may now be obtained by plugging m k and boldsymbol k (Eqs.12 and 13) into Eqs.( 17) and ( 18) to obtain m and and then using these expressions in Eq. ( 7).Although these operations can be performed numerically, there is a computational advantage in processing the equations symbolically.Details are given in Appendices B and C.

Short discussion of possible advantages over the independent emulator approach
Now that the notation and relevant concepts have been introduced, the potential advantages and drawbacks of the PC emulator and the independent emulator approaches may be briefly summarised as follows: -The PC emulation is based on the calibration of Gaussian process models per output field (temperature, precipitation, GDD).We use = 10 (justified below).
We note that computational cost may be saved by using the same length scales for all Gaussian processes, though in practice computational costs remain affordable even when using independently optimised length scales.Consequently, the impact of the same lengthscale assumption may be assessed more easily than it would be if we used 2048 independent emulators (i.e. the number of grid points) for each output.
-As the PC emulation approach requires fewer emulators, more time can be spent individually validating each of them.-In some regions there may be only small variation in the simulated output as the input parameters change.If independent emulators are used for each grid cell, estimating the hyper-parameters for these cells can be difficult without applying some sort of parameter regularisation, and, moreover, the computational effort of building the emulator is unnecessary (as the output is constant).The global principal component emulator is therefore preferable in these situations, as these constant regions are automatically accounted for in the principal component variance decomposition.
-Finally, the PC emulator provides covariance indices between any two simulator outputs.It therefore allows us to analyse the spatial structure of the simulator response to individual and combined factors.

Sensitivity to initial conditions
For all ensemble members but two (experiments 20 and 27), the runs with distinct initial conditions converged to the same output, modulo small variations that can be attributed to sampling variability (see Supplement).Experiment 27 shows a higher-amplitude variability pattern, but clearly oscillates around one mean value and, as we will shortly see, this mean is correctly captured by the emulator.It is therefore kept without further discussion for all subsequent analyses.Experiment 20 used the lowest configuration of obliquity ( 22• ).This is lower than any obliquity that occurred during the Pleistocene (22.07 • following Laskar et al., 2004).In this configuration, LOVECLIM develops a slow oscillation pattern that may be reminiscent of Dansgaard-Oeschger oscillations: millennial transitions between a warm and a cold North Atlantic phase, with fast warming and slow cooling (Fig. 2).The phenomenon is a known feature of LOVECLIM (Goosse et al., 2002;Loutre et al., 2014).It can be described as the apparition of a cold North Atlantic phase that is being visited stochastically and increasingly frequently as obliquity decreases.This cold phase is being visited shortly once during an entire (additional) experiment at obliquity of 22.5 • (lowest 7th percentile of Pleistocene obliquities).The oscillation itself could be of physical relevance for past climate variability, but the limits of the phase space region in which the oscillation occurs are also likely to depend on many other uncertain parameters of the model.One possible action would be to use a sequential design strategy to delineate the region of occurrence of the phenomenon and develop an emulator aimed at characterising this oscillation.In particular, historymatching theory provides adequate concepts and methods to this end (Williamson et al., 2013).Given the likely sensitivity of the oscillation on model parameters, the significance of this enterprise for palaeoclimate interpretation is unclear.
We choose to ignore the experiment for the time being (the following diagnostics ignore experiment 20), but we discuss the possible consequences of this decision in the final discussion.In statistical terms, we provisionally condition the analysis on the hypothesis that these oscillations do not occur in the phase space.

Validation and choice of PC emulator
We concentrate on three output fields: annual precipitation, growing degree days (GDD), and annual mean temperature.GDD is defined here as the annual sum of daily temperatures (in degrees Celsius) exceeding 0 • C. It is used as a calendarindependent indicator of summer intensity and length in extra-tropical regions, and the vegetation model VECODE, used in LOVECLIM, uses GDD and annual precipitation to predict the dynamics of vegetation (Brovkin et al., 1997).For simplicity, GDD is estimated from monthly means, that is, 30 times the sum of monthly mean temperatures for which this temperature is above zero.In equatorial and subtropical regions this information is equivalent to annual mean temperature.Furthermore, we use the logarithm of annual precipitation rather than precipitation, as the former is closer to being Gaussian distributed than the latter.The decomposition in principal components is effective, with 99 % of the variance on average over all grid points captured by the first four (annual temperature) to eight principal components (Fig. 3).As discussed in the methodology section, two options are considered for the estimation of the length scales and nugget variance and ν.We first attempt to use different correlation parameters obtained by maximising the penalised likelihood for each principal component emulator, independently of the others.It is then observed that the Gaussian process likelihood decreases with the index of the PC (Fig. 4).This is a natural result if we think of the fact that, as the index of the PC increases, its spatial pattern becomes more noisy and dependent on idiosyncrasies of the analysis such as the specific experiment design, experiment length, and initial conditions.They are thus less informative about the model itself, and scores are more difficult to predict with a smooth Gaussian process.The likelihood stabilises around PC 10 to a minimum value that indicates that the calibrated GP is not more informative than assuming independence of outputs on inputs.We therefore use = 10.The alternative approach consists in using the same correlation parameters for all PCs, in which case they are found by maximising the product of Gaussian process likelihoods.
Our evaluation strategy is based upon the leave-one-out cross-validation approach: for each member of the experiment design, a PC emulator is trained using the remaining design runs (using the correlation parameters estimated above).The means and standard deviations of the resulting emulator are then found for the design member left out.Figure 5 shows (bars) the number of grid points correctly predicted within the central 66th, 95th, and 99th credible intervals.A wellcalibrated emulator would, on average, correctly predict 66, 95, and 99 % of the points, respectively, in each category.Some remarks are in order: 1. Based on this diagnostic only, using constant rather than PC-specific correlation parameters does not significantly affect the overall performance.This is explicitly shown for GDD, but is also true of the other fields.
2. However, all fields exhibit an excessive number of predictions outside the central 99th credible interval.Annual mean precipitation is, in this respect, less well predicted than the others, perhaps not surprisingly given that precipitation responds less straightforwardly to insolation changes than temperature.
3. There are an excessive number of predictions within the 66th central credible interval.

Astronomically forced variance vs. other effects
The total variance V resulting from the astronomical forcing can be estimated from Eq. ( 11).Here we compare the estimated total variance with other effects that may broadly be described as sources of uncertainties on this quantity (Fig. 6).
For the assessment specific to this subsection, we considered the uniform distribution ρ(x) = 1 over the cube in order to be able to provide analytical integrals, and hence isolate the effects associated with Monte Carlo sampling.The following observations can be made: -The variance associated with the input factors largely dominates other sources of uncertainty.
-The error caused by the Monte Carlo approximation of the integrals, estimated by comparing these Monte Earth Syst.Dynam., 6, 205-224, 2015 www.earth-syst-dynam.net/6/205/2015/q q q q q q q q q q q q q q q q q q q q q q q q 0 5 10 15 20 q q q q q q q q q q q q q q q q q q q q q q q q Annual precipitation # of PCs q q q q q q q q q q q q q q q q q q q q q q q q 0 5 10 15 20 1e−04 1e−03 1e−02 1e−01 1e+00 PC q q q q q q q q q q q q q q q q q q q q q q q q Growing Degree Days # of PCs q q q q q q q q q q q q q q q q q q q q q q q q 0 5 10 15 20 1e−04 1e−03 1e−02 1e−01 1e+00 PC q q q q q q q q q q q q q q q q q q q q q q q q Annual temperature # of PCs Shown are the cumulated normalised variances resolved by the principal components (red) and the remaining variance (blue), which is modelled as white noise in Eq. ( 18), as a function of .Quantities are grid-cell averages.
q q q q q q q q q q 2 4 6 8 10 0.00 0.05 0.10 0.15 λ 1 (no unit) PC PANN q q q q q q q q 2 4 6 8 10 0.00 0.05 0.10 0.15 λ 1 (no unit) PC GDD q q q q q q q q q q 2 4 6 8 10 0.00 0.05 0.10 0.15 λ 1 (no unit) PC TANN PC q q q q q q q q q q 2 4 6 8 10 λ 3 (degrees) PC q q q q q q q q q q 2 4 6 8 10 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 λ 3 (degrees) PC q q q q q q q q q q 2 4 6 8 10 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 λ 3 (degrees) PC PC q q q q q q q q q 2 4 6 8 10 0.00 0.05 0.10 0.15 ν (no unit) PC q q q q q q q q q 2 4 6 8 10 0.00 0.05 0.10 0.15 ν (no unit) PC q q q q q q q q q q 2 4 6 8 10 0.00 0.05 0.10 0.15 ν (no unit) PC PC q q q q q q q q q q 2 4 6 8 10 Log likelihood PC q q q q q q q q q q q q q q q q q q q q 2 4 6 8 10 Log likelihood PC q q q q q q q q q q q q q q q q q q q q 2 4 6 8 10 0 20 40 60 80 100 120 Log likelihood PC q q q q q q q q q q PC Figure 4. Gaussian process parameters maximising the penalised log likelihood for three variables: GDD, log(annual precipitation), and annual temperature, either (black) optimised for each PC independently or (blue) optimised based on the product of the likelihoods of the first 10 PCs, assuming that the same correlation parameters are used on all PCs.The maximised log likelihood associated with each PC is given for reference.
Carlo integrals with the analytical solution in the particular case of a uniform distributions for x (ρ(x) ∝ 1), is of the order of 0.5 % of the total variance.
-Slightly different variance estimates are obtained depending on whether we use the same length scales for all PC or not, or whether we use independent emulators or the PC emulator.The difference over a grid point is on average 2.5 % of the mean grid-point variance.Different estimates will also be obtained with the independent emulators over all grid boxes depending on the length scales being used; only one length scale was tested here.The Supplement further shows that the patterns of the global sensitivity indices are similar regard- experiment index q q q q q q q q q q q q q q q q q q q q q q q q q q 0 5 10 experiment index q q q q q q q q q q q q q q q q q q q q q q q q q q 0 5 10 15 20 GDD with constant hyperparams experiment index q q q q q q q q q q q q q q q q q q q q q q q q q q 0.000 experiment index q q q q q q q q q q q q q q q q q q q q q q q q q q 0.00 less which emulator is used, so that this choice is of no consequence for the scientific interpretation.
-The term "GP var" explicitly refers to S tot defined in Eq. ( 11): shown here is the mean of the diagonal of this matrix, equal to the mean of the Gaussian process variance over all grid points, averaged over the input space.Again this is a small term, which is of the same order of magnitude as the grid-box mean uncertainty associated with the choice of initial conditions.
-The absolute value of the synergy term, measured as the difference between the total variance and the sum of the mean sensitivity indices, is also of the same order as the different sources of uncertainty just discussed.

Grid-point-wise variance analysis
The variance indices over the different grid points are the diagonal elements of the matrices V, T {e } , and T ε .They pro-vide essentially the same information as could be obtained using emulators independently calibrated over all grid boxes, as follows (see Fig. 7): 1. Precipitation is mainly controlled by precession over western Africa and Australia.This is expected given the known control of precession on monsoon dynamics (e.g.Zhao et al., 2007).The absence of large variance patterns in South East Asia and South America, otherwise expected, may be a consequence to the limitations of LOVECLIM in simulating tropical weather systems.Note also the significant influence of obliquity in the most western part of North Africa.
3. Annual mean temperature has the highest variance near the poles.Again, Northern Hemisphere temperature is equally controlled by precession and obliquity, while obliquity dominates the variance in the Southern Hemisphere.

Fingerprint analysis
The PC emulator, however, allows us to go one step further than the independent emulator strategy.As the full covariance matrices V and T i are available, linear fingerprints may be obtained by performing an SVD of these covariance matrices (Fig. 8).Specifically, the eigenvectors T {e } and T ε are called the "fingerprints" of precession and obliquity, respectively.The first fingerprint of obliquity explains more than 90 % of the variance in all three variables considered here.Precession aggregates two inputs (e sin and e cos ).It can therefore be expected to have at least a second significant fingerprint.This is the case, but this second component represents less than 30 % of the variance.We come to that shortly.Compared to the point-wise variance analysis above, the main advantage of the fingerprint analysis is that it provides information on the in-phase or anti-phase relationships between climate variables, as follows: -Obliquity produces in-phase effects on monsoon-related precipitation both in the Northern and in the Southern Hemisphere (compare Africa and northern Australia, while precession causes opposite-phase effects.This pattern is easily explained by reference to insolation.Indeed, obliquity causes in-phase responses of summer insolation in both hemispheres, while precession causes anti-phased responses. -Obliquity produces an equator-pole see-saw response of annual mean temperature, with, however, a weakeramplitude response at the equator than in the extratropical regions.Again, this consistently reflects the pattern of annual mean insolation (Loutre et al., 2004).
Note also that fingerprints of precession and obliquity are not orthogonal and thus cannot be readily recovered by principal component analysis of the model outputs.For reference, we estimated the variance decomposition of the PC scores associated with precession and obliquity (decomposition and maps of principal components available in the Supplement).The variance analysis reveals a mixture of precession and obliquity effects on each principal component.
Coming back to precession, we expect the simulated response phase to differ from place to place.To illustrate this point, we plot the emulated precipitation as a function of the longitude of the perihelion for three points along the African monsoon flow.We assume obliquity and eccentricity typical of the Holocene (Fig. 9), and indicatively denote longitudes of perihelion corresponding to the time of the beginning of the Holocene (11 000 years ago), as well as 6000 years ago, a reference period used for model intercomparison exercises (Braconnot et al., 2007).The timing of the maximum response is gradually shifted towards a late phase response as one travels northwards.This observation can be explained by considering the seasonal development of monsoon dynamics, along with the course of the zenithal sun.According to this analysis, the most favourable epoch for a "Green Sahara" experiment with a global climate model would therefore be around 9000 years, corresponding to the choice of early modelling experiments on this subject (Street-Perrott et al., 1990).

Detection of fast changes
Assuming that the climate system responds fast enough to changes in astronomical forcing, the time evolution of the climate system may in principle be simulated by forcing the emulator with a realistic evolution of astronomical forcing.Unfortunately, the output cannot be readily compared to observations because we are neglecting here other significant forcing elements, such as changes in land ice cover and greenhouse gas concentrations (Araya-Melo et al., 2015, as well as Bounceur, 2015).This exercise may, however, help us to detect regions where, potentially, the climate system may respond with steep gradients to the smooth astronomical forcing changes in interglacial conditions.Specifically, we estimate the maximum rate of change of a climate variable, expressed in terms of units per thousands years (Fig. 10).In order to enhance the palaeoclimatological interest of this discussion, we compare the experiments with interactive vegetation to those with fixed land-surface properties, hereafter referred to as VOFF.As a reading guide, a climate variable responding linearly and exclusively to precession with standard deviation 1 would show a maximum rate of change of 0.82 ky −1 ; a variable responding linearly and exclusively to obliquity with 1 standard deviation would show a maximum rate of change of 0.35 ky −1 .On this basis, comparing Fig. 10 with 7 allows us to detect two regions of potentially rapid changes: The western Sahara: rate of changes expressed on the log scale are of the order of 0.4 per thousand years, that is, about 50 % precipitation change in 1000 years.This is a well-known feature explained by feedbacks between vegetation and climate (Brovkin et al., 1998), discussed specifically in LOVECLIM by Renssen et al. (2003).
The North American sector of the Arctic, including both northern Canada and sea-ice-covered regions: rates of temperature changes are of the order of 3 • C per thousand years.
These hotspots almost disappear with the fixed-land-cover scheme, underlining the role of vegetation response.However, the fact that the North American hotspot region extends over the Arctic Ocean also suggests a role for sea-ice cover, which further amplifies vegetation-induced effects.Such interactions between Arctic vegetation and sea ice have already been suggested to have played a role in Arctic climate change during the Holocene (Ganopolski et al., 1998).

Discussion
Let us first review and comment upon the results of palaeoclimate significance presented here.Naturally, they are conditional on the use of the specific simulator considered here (LOVECLIM) and must be considered critically given that LOVECLIM is an imperfect representation of reality: 1. Precession and obliquity both contribute to annual temperature.Precession generally has a greater effect in the Northern Hemisphere and tropical regions, and obliquity is the dominant forcing in the Southern Hemisphere.  Figure 9. Emulated annual precipitation for three locations at the North Atlantic/African sector as a function of the longitude of the perihelion.For ease of interpretation, the longitude of the perihelion is denoted as the time of the year at which perihelion (closest point to the Sun) is reached.
insolation, but it may affect the annual mean climate by acting on the seasonal cycle of albedo associated with sea-ice, snow, and vegetation feedbacks.The latter two naturally depend on the presence of continental masses, which occupy a larger fraction of the Northern than Southern Hemisphere.This dichotomy between Southern and Northern Hemisphere responses was previously noted by Yin and Berger (2012), based on LOVECLIM simulations of previous interglacial periods, who then referred to it as one of the elements needed to explain the occurrence of the "Mid-Brünhes Event" (Yin, 2013).
It is also consistent with observations, namely -the prominence of obliquity signals in Southern Hemisphere records, and more particularly Antarctic cores, be it CO 2 concentration (Petit et al., 1999;Siegenthaler et al., 2005;Luethi et al., 2008) or deuterium excess (Vimeux et al., 2002); -the contrasting dynamics between southern records and northern continental records, such as Baïkal's, during isotopic stage 11 (Prokopenko et al., 2002).
2. GDD is used here as a measure of summer length and intensity.We considered it because it is used in VE-CODE as a predictor for vegetation changes.This quantity is also mathematically equivalent to the positivedegree days (PDD) index used as a predictor of net snow accumulation balance over ice sheets (e.g.Pollard and DeConto, 2005).We see here that GDD is, in the Northern Hemisphere, approximately equally sensitive to precession and obliquity.Crucifix (2011), based on Berger (1978a), noted that the Milankovitch's caloric season insolation is also equally sensitive to precession and obliquity.Hence, this result is consistent with the proposal by Ruddiman (2007) to use caloric season insolation as a predictor for ice age inception.
3. We find a fairly strong obliquity effect on North African precipitation.This has been noted before, in particular by Tuenter et al. (2003)   to changes in large-scale atmospheric circulation patterns.
4. Vegetation feedbacks substantially increase the climate change rate in the Arctic and in the Sahel.The response is slightly non-linear (best seen, for example, in Fig. 9), but not to the point of generating multiple steady states in response to a same astronomical configuration.In other words, all 27 experiments of the design converge to the same state (the particular case of experiment 20 is discussed in one of the next points).This is consistent with previous transient simulations with this model (Renssen et al., 2003), but we note that Brovkin et al. (1998) report multi-stability of the western Sahara in the early Holocene, leading to a bifurcation associated with the abrupt desertification of the Sahara during the mid-Holocene (Claussen et al., 1999).It remains unclear whether possible multiple stable states persist when more sophisticated land-surface-vegetation schemes associated with finer spatial resolution are used (see, for example, Kleidon et al., 2007;Dekker et al., 2010).On the other hand, model intercomparisons support the existence of a single stable state in the high latitudes (e.g.Brovkin et al., 2003).
5. Our methodological approach allows us to document response phases of precipitation patterns associated with the African monsoon.In particular, we found that the northward penetration of African monsoon is at a maximum when the perihelion is reached in August.Indeed, in this configuration, average levels of spring insolation can prevent excessive warming of the ocean surface during this season, while high positive June-July-August insolation anomalies effectively enhance the warming of the subtropical continent.This combination maximises the contrast between ocean and continental temperature during the monsoon season.Again, the result needs to be qualified with the usual cautionary remarks about the simplification of tropical dynamics in a model like LOVECLIM.
6. Instabilities of the North Atlantic circulation develop at very low obliquities.The effects of these instabilities have not been taken into account in the variance diagnostics discussed here.We note in particular that obliquities as low as those necessary to trigger the oscillations would, in the real world, also be associated with the development of ice sheets.The latter could further complicate the dynamics in the North Atlantic region.We therefore limit ourselves to observe that such oscillation dynamics, if they were to occur in specific interglacial configurations, would dominate the astronomical sources of variance examined here.
A perhaps more surprising result of our analysis is the smallness of the synergy terms.Crucifix and Loutre (2002) outlined the significance of a synergy between precession and obliquity during the last interglacial period.They noted that effects of precession and obliquity may combine and produce non-linear effects associated with the taigatundra transition in the Arctic area, not incompatible with what we find here about rates of changes in the Canadian Arctic.We explain this paradox by observing that metrics provided by global sensitivity analysis are aimed at determining whether the total effect variances, assessed over the whole Pleistocene, add up linearly or not.In LOVECLIM, total variances indeed add up roughly linearly.This does not exclude that non-linear effects may episodically dominate at critical periods.

Conclusions
We have presented a global sensitivity analysis of the effects of astronomical forcing on the climate model LOVECLIM in interglacial conditions.The work relies on the methodology of PC-Gaussian process emulation to explore the input space and deliver spatially resolved variance indices.In particular, we introduce the fingerprints as the eigenvectors of the covariance indices obtained from global sensitivity analysis to provide a spatial description of the effects of the individual factors.From a palaeoclimatological perspective, the results shown here are broadly consistent with the current understanding of Earth's climate response to the astronomical forcing.Compared to standard approaches based on a small number of simulations for well-defined past epochs, the methodology presented here allows us to identify more systematically regions susceptible to experience rapid climate change in response to the smooth astronomical forcing changes, and examine the response phase of climate change to precession.We do not have to rely on transient experiments, so the methodology is readily applicable to more complex climate models, but we have to rely on the assumption of quasi-stationarity of the climate response to the astronomical forcing.Although the mathematics are fairly straightforward, the emulation approach requires considerable care in its implementation.We discussed the effect of the choice of length scales and the type of emulator and estimated uncertainties associated with specific computational aspects such as the Monte Carlo estimates of integrals, and have concluded that the PC emulator is a reasonable option.We therefore recommend its use for future applications.

Figure 1 .
Figure 1.Two-dimensional projections of the experimental design.The experiment marked in red (experiment 20) was discarded from the analysis.

Figure 2 .
Figure 2. Slow oscillations developing in experiment 20 (e = 0.040, = 334.6 • , ε = 22 • ).The surface annual temperature over one of the North Atlantic grid points is shown (inset, in • C) along with the geographic distribution of the difference between the warm and the cold phases.The horizontal red line in the inset is the emulator prediction, calibrated on the 26 remaining experiments.

Figure 3 .
Figure 3. PCA decomposition of the logarithm of annual precipitation, GDD, and annual temperature fields based on the experiment design.Shown are the cumulated normalised variances resolved by the principal components (red) and the remaining variance (blue), which is modelled as white noise in Eq. (18), as a function of .Quantities are grid-cell averages.

Figure 5 .
Figure 5. Evaluation of the PC emulators.The bars give the fraction of grid points for which the emulation correctly predicts the value of the experiment left out of the training set, within the 66th, 95th, and 99th inner quantiles of the distribution.The horizontal lines indicate the theoretical position of the percentiles for well-calibrated emulators.Dots provide root mean squares of the differences between predicted and actual values.The graphical layout is adapted from the recommendation of the "Modelling Uncertainty in Computer Model" project, http://mucm.aston.ac.uk/MUCM/MUCMToolkit/index.php?page=ExamMultipleOutputsPCA.html.

Figure 6 .Figure 7 .
Figure6.Total variances associated with inputs considered grid-point-wise, and then averaged over all grid points.These variances are then compared to different elements of variances associated with the PC emulator, namely variance in the discarded principal components, error estimated from Monte Carlo sampling (estimated by comparison with analytical integrals), Gaussian process variance ( tot ), difference in variance associated with the choice of correlation parameters (constant for all PC, or PC-dependent), or the choice of emulator (PC emulator vs. independent emulators), difference due to the choice of initial conditions (experiment 20 excluded), and synergy term defined as V − T {e } − T ε .

Figure 8 .
Figure8.First eigenvectors of T {e } and obliquity T ε for annual precipitation (log) and annual surface temperature, along with the fraction of the explained variance by these eigenvectors.These eigenvectors are referred to here as the fingerprints.

Figure 10 .
Figure10.Maximum rate of change, in units per thousand years, estimated using the PC emulator and assuming a quasi-stationary climate and interglacial conditions.