"Global sensitivity analysis of Indian Monsoon during the Pleistocene"

The sensitivity of Indian Monsoon to the full spectrum of climatic conditions experienced during the Pleistocene is estimated using the climate model HadCM3. The methodology follows a global sensitivity analysis based on the emulator approach of Oakley and O’Hagan (2004) implemented following a three-step strategy: (1) 5 develop an experiment plan, designed to e ﬃ ciently sample a 5-dimensional input space spanning Pleistocene astronomical conﬁgurations (3 parameters), CO 2 concentration and a Northern Hemisphere glaciation index, (2) develop, calibrate and validate an emulator of HadCM3, in order to estimate the response of the Indian Monsoon over the full input space spanned by the experiment design, and (3) estimate and interpret 10 sensitivity diagnostics, including sensitivity measures, in order to synthesize the relative importance of input factors on monsoon dynamics, estimate the phase of the monsoon intensity response with respect to that of insolation, and detect potential non-linear phenomena. Speciﬁcally, we focus on four variables: summer (JJAS) temperature and 15 precipitation over North India, and JJAS sea-surface temperature and mixed-layer depth over the north-western side of the Indian ocean. It is shown that precession controls the response of four variables: continental temperature in phase with June to July insolation, high glaciation favouring a late-phase response, sea-surface temperature in phase with May insolation, and continental precipitation in phase with 20 July insolation, and mixed-layer depth in antiphase with the latter. CO 2 variations controls temperature variance with an amplitude similar to that of precession. The e ﬀ ect


Introduction
Palaeoclimate observations and simulations with climate models consistently show a response of the Indian Monsoon to climate and forcing changes associated with 5 glacial-interglacial cycles. In particular, many authors have used general-circulationmodel (GCM) simulations to study the sensitivity of Indian Monsoon to changes in astronomical parameters (e.g. Kutzbach and Street-Perrott, 1985;Kutzbach and Guetter, 1986;Prell and Kutzbach, 1987;Anderson et al., 1988;Kutzbach and Liu, 1997;Braconnot et al., 2002Braconnot et al., , 2008Braconnot and Marti, 2003).
forcing dominates the monsoon response, which is maximum around 11-9 kyr BP. Kutzbach and Liu (1997) provide an early application of ocean-atmosphere models to palaeoclimate studies, based on an asynchronous coupling procedure. Again, the increase in the seasonal cycle insolation in the Northern Hemisphere increased the sea surface temperature in the tropical Atlantic 6000 years ago in late summer, which 15 in turn further enhanced the summer monsoon precipitation in northern Africa by more than 25 %, compared to simulations with prescribed, pre-industrial sea-surface temperatures. Braconnot and Marti (2003) used a fully coupled ocean-atmosphere model, and specifically focused on Indian monsoon and its effect on Indian and South-East 20 Asia climatology. They presented three experiments specifically chosen to study the sensitivity of the seasonality of the monsoon signal to the time at which perihelion is reached. Namely, an "early-phase" configuration (perihelion reached in April) produces a stronger monsoon, which occurs earlier in the year than a "late-phase" configuration (perihelion reached in September). 25 Intercomparison exercices, based on experiments with different climate models, allowed to further test the dependency of these conclusions to model choice. Specifically, Joussaume et al. (1999) and Braconnot et al. (2002) investigated the enhancing of the summer monsoon in the Northern Hemisphere at 6 kyr BP based Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | on a set of 17 atmosphere-only simulations, as part of the Paleoclimate Modeling Intercomparison Project (PMIP). It was found that in all simulations, the 6 kyr BP orbital configuration resulted in a continental warming over the Eurasian continent, which enhances the Indian Monsoon. They also found that the changes in monsoon precipitation depends on the magnitude of the continental warming. Later, Braconnot 5 et al. (2008), using seven ocean-atmosphere couple simulation models, compared the monsoon response to orbital parameters changes between the Eemian and Holocene period. Their study confirmed the strong relationship between increased seasonality of insolation in the Northern Hemisphere and monsoon amplification.
It is thus clear that astronomical forcing (eccentricity, longitude of perihelion, 10 and obliquity), CO 2 and ice boundary conditions determine the history of monsoon dynamics. One general difficulty is to disentangle the individial and combined effects of these five factors influences at reasonable computing cost. A classical factorial experiment, with only 2 distinct levels per factor (a minimum and a maximum range) would already require 32 experiments. Three levels would require 243 experiments. 15 Fortunately, the theory of experiment design and global sensitivity analysis with computer models provides us with strategies to address this problem. The method featured here follows Oakley and O'Hagan (2004), who further refer to the earlier contributions of Sacks et al. (1989) and Homma and Saltelli (1996). The principle rests on the combination of an adequate experiment design with a statistical model. The 20 purpose of the latter is to interpolate the GCMs outputs in order to produce appropriate numerical and visual diagnostics.
The statistical model, hence, act as a fast surrogate of the GCMs, and for this reason it is commonly named "emulator" in the literature (O'Hagan, 2006). Specifically, the term emulator refers to the following properties (O'Hagan, 2006;Petropoulos et  once the emulator is built, it is not necessary to perform any additional runs with the model.
The objective of the present study is to develop a first experience in the application of the emulator framework for a sensitivity analysis of palaeoclimates relevant for the entire period of the Pleistocene. More specifically, we address the sensitivity of the 5 Indian Monsoon to the three components of the astronomical forcing, plus carbon dioxide concentration and ice boundary conditions.
To this end, an ensemble of 61 simulations with varying input parameters are run with the climate model HadCM3 Gordon et al. (2000). The climate analysis is inspired from the earlier work of Zhao et al. (2005), and emphasis is set on benefits of the emulator 10 framework for understanding the monsoon response.
The paper is structured as follows. Section 2 provides a brief description of the emulator and the simulations used. Section 3 details the results of applying the emulator to the Indian monsoon region, focusing on the relations of the parameters under study and their impact in the climate of the monsoon. We also study the specific 15 influence of ice sheet topographic forcing. The conclusion follows in Sect. 4.

Experiment design
Five inputs factors are considered here: the three elements of astronomical forcing (eccentricity e, longitude of perihelion and obliquity ε), the concentration in carbon 20 dioxide (CO 2 ), and a variable called the ice or glaciation level, which combines ice and orography forcings associated with the presence of continental ice in the Northern Hemisphere.
The three elements of astronomical forcing are combined under the form of e sin , e cos and obliquity ε. This choice is justified by the fact that these combinations 25 produce orthogonal patterns in the season-latitude space, and generally insolation at 1614 any point and time in year is well-approximated as a linear combination of those terms (Loutre, 1993 boundary conditions representative of glacial-interglacial dynamics. Pragmatically, we sampled these boundary conditions among the series prepared by Singarayer and Valdes (2010), and supplied to us by Paul Valdes. Level 1 corresponds to presentday conditions, and levels 2 to 11 are chosen such as to represent approximately ten equally spaced top altitudes of the North American ice sheet, within the glaciation 10 phase. One limitation of this design for the present purpose is that levels 3 to 11 effectively represent similar ice sheet areas -thus similar albedo forcing -even though they sample very different ice sheet volume (see Fig. 3). The next step is to effectively select the input factors combinations used in the simulation ensemble. Theoretical considerations and experience point to the latin 15 hypercube design (McKay et al., 1979;Morris and Mitchell, 1995;Urban and Fricker, 2010) as a good starting point for computer experiments. The principle, for a latin hypercube design of n elements, is to divide the ranges covered by each input factor into n distinct categories, each experiment sampling one of the n categories without replacement. We followed the standard practice of associating the latin hypercube 20 design with additional constraints (Joseph and Hung, 2008, e.g.): namely, maximize the minimum distance between every two pair of points (specifically, the Euclidean distance in the normalized input space), and maximize the covariance matrix of the design, again expressed in the normalized input space. By maximizing the distance between points we ensure the experimental design to be space filling.

25
In the present context, two additional constraints need to be accounted for in order to avoid sampling unrealistic inputs that would be uninformative for the sensitivity analysis of climate over the Pleistocene: exclude forcings with e > 0.05, and exclude combinations of high CO 2 and high glaciation levels (and conversely), delineated by an Introduction ellipse with great and small axes as shown on Fig. 1. To satisfy these constraints, the design points generated by the latin hypercube sampling procedure and lying in the exclusion zone are geometrically projected on the allowed region. This procedure may break some of the original properties of the design (maxi-min and orthogonality), but it offers the practical advantage of enhancing the coverage of the input space near its 5 boundary.
Note that this design is in principle suitable for continuous factor ranges only. The glaciation level used for experiments is an integer obtained by rounding the value obtained by this process to the closest integer. Designs specifically adapted for input spaces mixing categorical and continuous variables could best be implemented in the 10 future (see, e.g., MacCalman, 2013, for an up-to-date review). Table 1 lists the simulations with their input parameters. The choice of 61 members is a conservative implementation of the recommendation of 10 experiments per input factors (Loeppky et al., 2009). In fact, a first 57 member design was produced using the method above, to which 4 members were added (exp. 20-23). These experiments 15 are idealised orbital changes that were performed during the first phase of this project in order to explore locally the model sensitivity to astronomical forcing.

Climate simulator
The climate model -referred to in this context as the simulator -is the General Circulation Model HadCM3 (Gordon et al., 2000) Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | ran 300 years, accounting for ice sheet topography from the beginning. Typical residual deep-ocean temperature trends are of the order of 10 −4 • C year −1 . The last 100 years of all simulations with orographic forcing were retained for analysis. Over this interval, the top-of-the-atmosphere imbalance ranges between −0.2 and −0.1 W m −2 . The last 100 years of the experiments section without orographic 5 forcing are also used for an investigation of the specific effect of the orographic forcing (cf. Sect. 3.6).

Emulator
The emulator is a statistical model calibrated on the simulator output. Its role is to predict simulator outputs for untried experiments. We follow here the Gaussian process 10 framework of Sacks et al. (1989), Kennedy and O'Hagan (2000) and Oakley and O'Hagan (2002). Let x j be the input vector associated with the j th component of the experiment design. The output of the climate model is modeled as a stochastic process combining a global response function (the regressors) with a local component. It is fully specified 15 by the mean m and a covariance V function. They have the following priors: where c(x, x ) is the Gaussian process correlation function, σ 2 its variance, h(x) 20 is a (q × 1) vector of a priori known regression functions and β is the vector of corresponding regression coefficients. Note that the () is used to denote a horizontal vector. Let y be the vector of actual outputs, obtained by running the model at the n design points. The posterior estimate of the simulator output at any input point x is given by Introduction  (Oakley and O'Hagan, 2002): The operator () T is the matrix transpose, The above expressions assume the vague prior (β, σ 2 ) ∝ σ −2 proposed by Berger et al. (2001) and used by, e.g., Oakley and O'Hagan (2002) and Bastos and O'Hagan (2009). Note that with this prior the posterior state distribution is a student-t distribution with n−q degrees of freedom, but it is close enough to being Gaussian to be considered 15 as such in the following discussion.
With this framework, the choices of the regression functions h(x) and the Gaussian process correlation function c are application dependent. This is where the user has the opportunity to inject knowledge on the expected response of the simulator.
For this application, linear regression is an adequate choice because the seasonal 20 and annual forcings are almost linear with the input factors, except possibly for glaciation level. Hence, h (x) = (1, x ). The correlation function c is the exponential decay with nugget, discussed in length in Andrianakis and Challenor (2012): a regularization antsatz to circumvent poor matrix conditioning (Pepelyshev, 2010) and mis-specification of the correlation function Gramacy and Lee (2012). In the present context it also has a more directly interpretable function. Indeed, the output of HadCM3 depends to some degree on a number of factors not explicitly accounted for here, namely initial conditions, experiment length and averaging time span. The effects of 5 these non-explicit inputs can be absorbed by the nugget.
Following Kennedy and O'Hagan (2000), hyperparameters are obtained by maximising the emulator likelihood (the expression used here is from Andrianakis and Challenor, 2012): More specifically, Andrianakis and Challenor (2012) recommend the use of a penalised likelihood, that imposes a restriction on the amplitude of the nugget: 15 where M(ν, Λ) is the Mean Squared Error between the training points and the emulator's posterior mean at the design points, and M(∞) is its asymptotic value at λ i → ∞. We use = 1.

Sensitivity measures
One of the early applications of Bayesian approach for emulators was to estimate 20 sensitivity measures that quantify the influence of a factor on the simulator output (Oakley and O'Hagan, 2004). Although originally developed to estimate uncertainties associated with uncertain inputs, the indices may be reinterpreted as indicators of variance generated by known variations in inputs. Call ρ(x), the time-wise occupation density of the input space. This occupation density along the components of the Introduction

Tables Figures
Back Close

Full Screen / Esc
Printer-friendly Version

Interactive Discussion
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | astronomical forcing can be estimated with histograms of long time series generated with known astronomical solutions, such as Berger (1978) For CO 2 and glaciation level we consider the following empirical distribution to broadly capture the observed covariance between CO 2 and glaciation level: where c , i are inputs standardised as follows: With this density at hand, it is now possible to estimate, for each factor (or combination thereof) a mean and a variance (Oakley and O'Hagan, 2004).
Define the function m p (x p ), where p refer to one or several components of x, as the expected mean of x given x p fixed, and obtained by integrating the emulator over the sub-space χ p covered by the remaining components: The conditional dependency on y introduced in Eq. (3) is dropped for readability. The notation x|x p means that the integral is made over all the remaining components of x with given x p .

20
The expected covariance of the Gaussian process given x p is given by a double integral: This quantity may further be integrated over the possible values of x p : where χ p is the domain sub-space covered by factor p. Define, consistently, 5 where the notation m 0 , V 00 implies that the space χ p that appears in Eqs. (9) and (10) is the whole input domain χ , and, finally, 10 Then, the emulator output variance associated with input factor(s) p may be measured either as -S p = Σ p − Σ 0 : this is the expected loss in output variance that would occur if factor p was fixed, all other factors varying, or -S p = Σ − Σ p : this is the expected output variance obtained by varying factor p, all 15 other factors fixed.
The distinction is especially important when ρ(x|x p ) depends on x p . Indeed, in this case the variance of the mean effect includes both a contribution from the isolated effect of the factor on the simulator output, and an implicit effect associated with the changes in the probability distribution of other factors associated with a change in x p .

20
The variance S p discards the second effect. It is therefore a better starting point to explore physical mechanisms. Following Homma and Saltelli (1996)  Note, finally, that S p may be further split into S m p + S V p , i.e., one can separate the contributions from m and V appearing in Eqs. (11)-(13). S m p is the output variance induced by expected changes of simulator output in response to changes in input p. This is thus the quantity of physical interest. S V p is an effect produced by the emulator variance, i.e., the uncertainty associated with using an emulator rather than 5 the simulator. Hence, S m p should only be considered to be significant if it is large enough compared to S V p .

Results
In order to study the Indian Monsoon, we define two regions: Northern India (NI), with coordinates 70-100 • E, 20-40 • N, and northwestern Indian Ocean (IND), with 10 coordinates 55-75 Zhao et al., 2005). The chosen regions are depicted in Fig. 2, in which the sea-level pressure and surface temperature of one of the simulations are shown. The NI region covers the Indian continent and part of the Himalayas (which is dry today), while IND covers the northwestern part of the Indian Ocean. 15 We focus specifically on four physical variables representative of the Summer Indian Monsoon process: June-July-August-September (JJAS) temperature and precipitation on the continental box, and JJAS sea surface temperature (SST) and mixed-layer depth on the Indian Ocean box. Over the experiment design, continental temperature varies between 15

Emulation validation
An emulator using all 61 experiments is calibrated using the procedure given in Sect. 2.3, with scales λ i (with i = 1, ..., 5) and nugget determined by maximization of the penalised likelihood. The performance of the emulator is then assessed following a leave-one-out cross validation approach, that is, constructing 60 emulators to predict 5 the experiment being left out. Figure 4 shows the result of this leave-one-out cross validation procedure for SST and mixed-layer depth only, the other variables being discussed later. This leads us to the following observations: 1. the optimal covariance scales are generally commensurate with the range 10 covered by the input factors. This is the ideal scenario, as it implies that the stochastic component of the Gaussian process is smooth, and may thus be suitably calibrated by the experiment design.
2. There are however some instances where the optimum covariance scales are much greater than the scale of the variables: this is observed on all output 15 variables for the response in CO 2 , and, to a lesser extent, for obliquity. A large covariance scale implies that response is linear with respect to the factor, which is indeed a realistic outcome for CO 2 , in the range considered. This is not a problem on its own. It simply informs the user that a sparser sampling of this factor would have worked as well. 3. The leave-one-out cross validation plot reveals two outliers for SST (experiments 11 and 40) and one outlier for mixed layer (experiment 40). We term "outliers" experiments generating output that is predicted with an error of more than 3 standard deviations when they are left out of the calibration procedure. The signature of these outliers is also visible on the mean effect plots (Fig. 6). These 25 plots, which will be commented on more in depth in Sect. 3.5, represent the mean effect such as given by Eq. (9) and assuming CO 2 fixed. The figure reveals convoluted contours, most notably the 26.25 and 26.5 • C isotherms on the SST plot and the 38.5 m iso-depth that conflict with our expectation of a smooth response structure.
Further inspection of the outliers reveals a clear warm/cold/warm pattern in the North Atlantic, and cooling over the rest of the ocean, exemplified here by comparing 5 experiments 11 and 15 (Fig. 11). This pattern has been seen before in HadCM3, most notably in early experiments of the Last Glacial Maximum (Hewitt, 2003). It was associated with an enhancement of the North Atlantic Overturning circulation cell, and can be annealed by addition of freshwater in the North Atlantic (Hewitt et al., 2006). Experiments 11 and 40 have, however, low to moderate glaciation levels. Hence, the 10 most reasonable explanation seems that the 100 year sampling procedure has picked up some rarely visited ocean circulation regime that affects the climate system globally. Although this regime may be of relevance for palaeoclimate analysis, it appeared here sufficiently anomalous to be considered off our focus. Consequently, the emulator was recalibrated using the remaining 59 experiments. 15 This new emulator with new scales λ i and nugget (see Table 2) then presents quite convincing validation scores (Fig. 5): 1. all emulators capture between 38 (mixed-layer) and 43 (continental temperature) of the left-one-out experiments within 1 standard deviation, and between 56 and 58 within 2 standard deviations, which roughly correspond to the 66 and 95 % 20 ratios expected for a normal distribution; 2. the errors normalised by standard deviation are compatible with a normal distribution based on the Shapiro-Wilk normality test, except for continental temperature (p = 0.03); 3. there is no error exceeding 3 standard deviations.
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | In addition, the nugget value obtained by maximizing the likelihood is overall consistent with what we know from the model variability. The nugget quantifies the uncertainty of the "observation", here related to the specific choice of one 100 year simulation sample as representative of the model mean stationary state. The residual error in the emulator is of the order ofσ 2 ν, but it can be estimated precisely by looking 5 at the posterior variance at design points. Here, the obtained nuggets induce residual uncertainties with standard deviations of 0.04 • C on continental temperature, 2.3 % on precipitation, 0.05 • C on SST, and 0.7 % on mixed-layer depth. All these values are consistent with the 100 year variances of the corresponding quantities in HadCM3. Thus, the emulator calibration procedure has been able to infer information on model 10 variability from an ensemble of simulation outputs that are all 100 years averages. This is quite remarkable. Specifically, continental summer temperature is primarily determined by precession, CO 2 and, to a lesser extent, ice volume. It shows no significant sensitivity to obliquity.

20
Continental precipitation is also mainly driven by precession and less to ice volume. Contrarily to temperature, it exhibits no sensitivity to CO 2 .
Similar to continental temperature, SST is primarily driven by precession and CO 2 and, to a lesser extent, ice volume. It also shows a larger response to obliquity. Finally, mixed-layer depth shows a pattern similar to precipitation, except that the response to obliquity is not significant compared to the sources of uncertainty induced by the emulation and sampling variance. Introduction  Figure 12 displays the effects of precession on the four variables retained for analysis. The choice here is to show the effects by fixing ice and CO 2 concentration at three distinct levels representative of the course of glaciation (from top to bottom): glaciation level 1/CO 2 = 280 ppm, glaciation level 5/CO 2 = 230 ppm and glaciation level 5 11/CO 2 = 180 ppm. Quantities are further averaged over obliquity. In order to ease the interpretation, the months representing the time at which perihelion is reached are written on the plots: June for = 90 • , September for = 180 • etc. That is, neglecting slow transient effects that could be associated with the deep ocean response, this graphical representation provides an indication of the phase lag between the climate 10 response and the precession forcing of insolation. Specifically, the temperature response is in phase with June insolation at low glaciation levels, and in phase with July insolation at mid-and high-glaciation stages. In addition, the low-glaciation response to precession shows a marked asymmetric pattern, with temperatures at June perihelion that are lower than a linear extrapolation 15 would have predicted.

Sensitivity to precession
This feature may physically be understood by considering the summer precipitation response. Precipitation enhances latent heat cooling when perihelion is around July. This effect gradually weakens as glaciation takes place and the total amount of precipitation declines, hence the drift towards a more linear response. At higher 20 glaciation levels the JJAS temperature response phase also aligns with July insolation.
The maximum precipitation is obtained when perihelion is reached in early July. Among the series of experiments shown by Braconnot et al. (2008), this is indeed the 126 000 yr BP experiment (i.e., July perihelion) experiment that shows the strongest precipitation response over India. about 6 m, consistent across different models, in 6000 yr experiments (September perihelion). Braconnot and Marti (2003) examined also two nearly-opposite precession configurations with the IPSL model, corresponding to perihelion in April and October, respectively, and they found a shoaling of the mixed-layer compared to the present-day (perihelion in January) in both cases. 5 Zhao et al. (2005) attributed the mixed-layer shoaling to a stratification effect involving the response of SST. On this point, our analysis reveals that the maximum SST response occurs when perihelion is reached in May. This is not so surprising given that the ocean thermal inertia generally imposes a lag of a few months between the forcing and the response. This response, however, induces an asymmetry between 10 perihelion in April and perihelion in October, the first one only showing anomalously high SSTs. This is consistent with the analysis of seasonal cycle response provided by Braconnot and Marti (2003).

Sensitivity to obliquity
The response of obliquity is mostly linear. The range of obliquity covered during the 15 Pleistocene induces negligible continental temperature response over the West-Indian box. It also induces a slight increase in precipitation. Regarding the Indian Ocean box, there is a somewhat larger effect on SST compared to continental temperature, but not significant. As for the mixed-layer depth, the response to obliquity is negligible.
In order to better understand obliquity effect, we considered the four ideal 20 cases experiments (simulations 20-23, see Table 1). Specifically, we discuss here experiments 22 and 23, termed OBL23 and OBL24. They use zero eccentricity, same CO 2 concentration and glaciation level, and differ by the configuration of obliquity (24 and 23 • , respectively). The temperature difference map, for JJAS, reveals the signature of obliquity-induced insolation changes, with a warming of Northern Hemisphere 25 continent, and slight cooling of significant areas of the tropical oceans (see Fig. 9

Sensitivity to CO 2 and glaciation level
The response of all variables to CO 2 is best captured by linear processes (optimal λ i largely exceeds the range covered by the experiment design). Hence, the contribution of CO 2 to the climate response may be estimated straightforwardly from the coefficientŝ β, given by Eq. (3). Namely, the continental temperature and SST responses to the  Figure 8 shows the response of continental temperature (left panel), sea surface temperature (middle panel) and mixed-layer depth (right panel) to the variations of CO 2 concentration and glaciation level. The temperature ranges covered by CO 2 and glaciation levels are of the order of 1 and 2 • C for the continent and ocean surface, respectively. The continental ice effect is mainly present between glaciation levels 1 15 and 3. With the ice sheet reconstructions used here, the ice area extent responsible for the shortwave forcing reaches almost its maximum value at glaciation level 3. Further increasing the glaciation levels affects climate predominantly through the orography forcing (cf. Sect. 3.6).
Interestingly, the figure shows that SST is insensitive to glaciation level. However, 20 Fig. 7 shows a dependency on it. What happens is that the signal is reverse for low and high glaciation levels. When integrating over obliquity, these signals are averaged, so that dependency on glaciation is masked.

Orographic effect
Finally, we consider the differences between the simulations with and without orography 25 forcing of the ice sheets. The latter is potentially important given that mountains and elevated land masses affect the atmospheric circulation and precipitation patterns, and then the whole climate system. To this end, an emulator was calibrated on the available present-day orography experiments. The net effect orography can then be seen in Fig. 10, where all four variables are plotted as a function of the glaciation level. Black solid lines show the respective variables obtained with the standard experiment design, while red solid lines show 5 the response obtained with the experiment design assuming pre-industrial orography, regardless of the presence of ice sheets. The value plotted is obtained from Eq. (9). Note that by construction this value is also implicitly a function of CO 2 concentration, which enters Eq. (9) via the factor ρ(x|x ice ). Dotted lines indicates a 1-σ deviation, in both cases, based on Eq. (10), using x p = x p .
A clear deviation is seen around glaciation level 3. This effect is due to the fact that, as explained in Sect. 2.1, levels 3-11 represent effectively similar ice sheet area, but significantly higher orography (see Fig. 3). Hence, the albedo forcing dominates over the lower range of glaciation levels (1-3), with decreasing temperatures, precipitation and mixed-layer shoaling. The orography-no orography differences appear more 15 markedly above indice 3: orography reduces the cooling trends, by as much as 1 • C on the continent at glaciation level 11, and even reverses the precipitation trend. It is known that ice orography forcing may impact monsoon precipitation regimes (Prell and Kutzbach, 1997;Yin et al., 2009), though to our knowledge the specific effect of Northern Hemisphere ice sheet orography on Indian monsoon is yet to be documented.

20
The warming signal caused by orography may be understood by considering the increase in surface potential temperature over elevated regions, similar to what is seen today over the Tibetan Plateau. Because of these high potential temperature, downsloping air is effectively warmer than it would be in absence of orography forcing, and contributes here to increase the Northern Hemisphere continental surface temperature.

25
Orographic forcing generally induces atmospheric circulation anomalies and effects on ocean circulation and stratification. Namely Fig. 10 suggests a weak positive effect on mixed-layer depth, though quite small compared the astronomical forcing effects. An in-depth analysis of these effects falls beyond the scope of the present contribution.

Conclusions
We performed a global sensitivity analysis of the climate response of the Indian Monsoon to the astronomical forcing (e sin( ), e cos( ), ε), CO 2 concentration and glaciation level. To achieve this, we make use of an emulator, that is, a statistical model calibrated 5 on the outputs of a set of well-chosen experiments, in order to explore and quantify globally the sensitivity of the model to each parameter and combinations thereof. The present study focuses on four variables: continental temperature, continental precipitation, sea surface temperature and mixed-layer depth. These variables were averaged for the JJAS season over northern India and northwestern Indian Ocean.

10
The method is divided in three steps: -Design an experiment plan. We adopted a latin-hypercube design, optimised following two constraints: maximise the minimum distance between two points in the input space, and maximise the determinant of the matrix of covariance between the input factors (orthogonality constraint). In addition, the design 15 excludes configurations with excessive eccentricity and unrealistic combinations of CO 2 and glaciation level.
-Calibrate and validate and emulator. The emulator is the Gaussian process emulator developed by Sacks et al. (1989), Kennedy and O'Hagan (2000) and Oakley and O'Hagan (2002). The validation was performed following a leave-20 one-out cross validation approach. Two experiments were excluded of the design as presenting an anomalous North-Atlantic convection pattern. The emulator calibrated on the remaining 59 experiments presents convincing validation scores.
-Quantify and visualise the individual and combined effects of the different factors on summer Indian monsoon, based on sensitivity indices and cross-section plots.

25
This analysis yielded the following conclusions: 1630 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | precession controls the response of four variables: continental temperature in phase with June-July insolation, high glaciation favouring a late-phase response, sea-surface temperature in phase with May insolation, and continental precipitation in phase with July insolation, and mixed-layer depth in antiphase with the latter.

5
-The effect CO 2 on continental temperature and SST is of similar size as that of precession on summer continental temperature and SST.
-Obliquity is a secondary effect, negligible on most variables except seasurface temperature.
-The effect of glaciation is dominated by the albedo forcing, and its effect on 10 precipitation competes with that of precession.
-The orographic forcing reduces the glacial cooling induced by the albedo forcing, and even has a positive effect on precipitation.
The originality of this study relies on the use of the emulator technique as a tool to provide reliable numerical results. We confirm that this technique has a large potential 15 for the analysis of climate model outputs. Indeed, the study of any climatic event is far from straightforward when several variables are taken into account. In addition to state-of-the-art climate models, careful statistical modelling may significantly enhance to information that can be inferred from a well-chosen set of experiments. This holds regardless of the region of focus or the climate model being considered.  Comput., 22, 713-722, doi:10.1007Comput., 22, 713-722, doi:10. /s11222-010-9224-x, 2012. 1619 Hewitt, C. D.: The effects of ocean dynamics in a coupled GCM simulation of the Last Glacial Maximum, Clim. Dynam., 20, 203-218, doi:10.1007/s00382-002-0272-6, 2003. 1624 The effect of a large freshwater perturbation on the Glacial Atlantic Ocean using a coupled general circulation model, J. Climate, 19, 4436-4447, doi:10.1175/JCLI3867.1, 2006. 1624 Fig. 7. Sensitivity analysis : shown is the standard deviation of model outputs ( √S ) of each variable, induced by variations in input factors during the Pleistocene. From left to right, top to bottom: continental precipitation, continental temperature, sea surface temperature and mixedlayer depth. The variance is partitioned between the effects associated with the response to output changes (grey) and Gaussian process variance (black), associated with using an emulator rather than direct simulator output. 33 Fig. 7. Sensitivity analysis: shown is the standard deviation of model outputs ( S) of each variable, induced by variations in input factors during the Pleistocene. From left to right, top to bottom: continental precipitation, continental temperature, sea surface temperature and mixedlayer depth. The variance is partitioned between the effects associated with the response to output changes (grey) and Gaussian process variance (black), associated with using an emulator rather than direct simulator output. Introduction 34 Fig. 8. Sensitivity to CO 2 and glaciation level. From left to right: continental temperature, sea surface temperature and mixed-layer depth. Fields were integrated over e sin( ), e cos( ) and obliquity.
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Fig. 9. Sea surface temperature difference for two ideal simulations. CO 2 concentration, glaciation level and precession remained fixed, the only difference being obliquity (23 and 24 • ). 35 Fig. 9. Sea surface temperature difference for two ideal simulations. CO 2 concentration, glaciation level and precession remained fixed, the only difference being obliquity (23 and 24 • ). Fig. 10. Orography -No orography difference. From top to bottom, left to right: effect on continental temperature, precipitation, sea-surface temperature, and mixed-layer depth, with orography forcing (black) and without (red). The dotted lines show one standard deviation of the emulator prediction. One may see a departure point from glaciation level 3 in all four fields, as this is at this point that orography forcing becomes the most significant. 36 Fig. 10. Orography-no orography difference. From top to bottom, left to right: effect on continental temperature, precipitation, sea-surface temperature, and mixed-layer depth, with orography forcing (black) and without (red). The dotted lines show one standard deviation of the emulator prediction. One may see a departure point from glaciation level 3 in all four fields, as this is at this point that orography forcing becomes the most significant. Introduction  Fig. 11. Shown is the sea surface temperature difference between simulations 11 and 15 (see Table 1). There is a clear warming pattern in the North Atlantic, which affects the mean sea surface temperature. 37 Fig. 11. Shown is the sea surface temperature difference between simulations 11 and 15 (see Table 1). There is a clear warming pattern in the North Atlantic, which affects the mean sea surface temperature. Introduction

Tables Figures
Back Close

Full Screen / Esc
Printer-friendly Version

Interactive Discussion
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Fig. 12. Sensitivity to ecos(̟) and esin(̟) for all fields. Each panel, from top to bottom, shows the four fields with a different configuration of glaciation level -CO 2 concentration. Top panels: glaciation level = 1 and CO 2 =280 ppmv. Middle panels: glaciation level = 5 and CO 2 =230. Bottom panels: glaciation level = 11 and CO 2 =180. All fields were integrated over obliquity.
38 Fig. 12. Sensitivity to e cos( ) and e sin( ) for all fields. Each panel, from top to bottom, shows the four fields with a different configuration of glaciation level -CO 2 concentration. Top panels: glaciation level = 1 and CO 2 = 280 ppmv. Middle panels: glaciation level = 5 and CO 2 = 230. Bottom panels: glaciation level = 11 and CO 2 = 180. All fields were integrated over obliquity.