Impacts of Tibetan Plateau uplift on atmospheric dynamics and associated precipitation δ 18 O

Palaeoelevation reconstructions of mountain belts have become a focus of modern science since surface elevation provides crucial information for understanding both geodynamic mechanisms of Earth’s interior and the influence of mountain growth on climate. Stable oxygen isotopes palaeoaltimetry is one of the most popular techniques nowadays, and relies on the difference between δ18O of palaeoprecipitation reconstructed using the natural archives, and modern measured values for the point of interest. Our goal is to understand where and how complex climatic changes linked with the growth of mountains affect δ18O in precipitation. For this purpose, we develop a theoretical expression for the precipitation composition based on the Rayleigh distillation and the isotope-equipped atmospheric general circulation model LMDZ-iso outputs. Experiments with reduced height over the Tibetan Plateau and the Himalayas have been designed. Our results show that the isotopic composition of precipitation is very sensitive to climate changes related to the growth of the Himalayas and Tibetan Plateau. Specifically our simulations suggest that only 40 % of sampled sites for palaeoaltimetry depict a full topographic signal, and that uplift-related changes in relative humidity (northern region) and precipitation amount (southern region) could explain absolute deviations of up to 2.5 ‰ of the isotopic signal, thereby creating biases in palaeoelevation reconstructions.


Introduction
Despite ongoing debates regarding the thermal and mechanical nature of mechanisms involved (Boos, 2015;Chen et al., 2014), the Himalayas and the Tibetan Plateau (hereafter TP) have long been considered to exert major influences on Asian atmospheric dynamics, notably by reinforcing South Asian monsoon and driving subsidence ultimately leading to onsets of deserts over central Asia (Rodwell and Hoskins, 2001;Broccoli and Manabe, 1992). Thus, reconstructing the history of Himalayas and TP uplift appears crucial to understand the long-term climate evolution of Asia. In addition, understanding the timing and scale of surface elevation growth is crucial for reconstructing the rate and style of this tectonic plates convergence (e.g. Royden et al., 2008;Tapponnier et al., 2001).
The difference between palaeoprecipitation δ 18 O detected from natural archives and modern values of the site of interest has been used to identify the effect of the surface uplift in numerous recent studies Cyr et al., 2005;Ding et al., 2014;Hoke et al., 2014;Mulch, 2016;Rowley and Currie, 2006;Rowley et al., 2001;Xu et al., 2013). In the absence of direct measurements of "palaeo" altitude-δ 18 O relationship in situ, stable isotope palaeoaltimetry is potentially hampered by the fact that the presumed constancy of altitude-δ 18 O relationships through time might not be valid. For instance for the Andes, not considering the impact of uplift on climate dynamics and related δ 18 O values has been shown to produce errors in palaeoelevation reconstruction reaching up to ±50 % (Ehlers and Poulsen, 2009;Poulsen et al., 2010). Regional climate variables and associated isotopic signal in precipitation can also be affected by global climate change (Battisti et al., 2014;Jeffery et al., 2012;Poulsen and Jeffery, 2011). Moreover, it has been suggested that climate-driven changes in surface ocean δ 18 O through the Cenozoic can also influence recorded values of precipitation δ 18 O over the continent and corrections has been applied in some studies . Over TP, mismatches between palaeoelevation estimations from palynological and stable isotope data (e.g. Sun et al., 2014) could be related to complex climatic changes and associated variations of altitude-δ 18 O relationship linked to the uplift, but still a detailed assessment of the consequences of topographic changes on precipitation δ 18 O is lacking.

Model simulations
We use an atmospheric general circulation model (GCM) developed at Laboratoire de Météorologie Dynamique, Paris, France with isotopes-tracking implement, called LMDZ-iso (Risi et al., 2010). LMDZ-iso is derived from the LMDz model (Hourdin et al., 2006) that has been used for numerous future and palaeoclimate studies Pohl et al., 2014;Sepulchre et al., 2006). Water in a condensed form and its vapour are advected by the Van Leer advection scheme (Van Leer, 1977). Isotopic processes in LMDZiso are documented in (Risi et al., 2010). Evaporation over land is assumed not to fractionate, given the simplicity of the model surface parameterisation (Risi et al., 2010). Yao et al. (2013) have provided a precise description of rainfall patterns over the TP, and showed LMDZ-iso ability to simulate atmospheric dynamics and reproduce rainfall and δ 18 O patterns consistent with data over this region.
LMDZ-iso is also equipped with water tagging capabilities, allowing us to quantify different moisture contributions from continental and oceanic evaporation sources. The advantage of this technique compared to typical backtrajectories methods is that it tracks the water rather than air masses, thus taking into account effects of phase changes. In our simulations five potential moisture sources are considered: (1) continental sources, (2) Indian Ocean, (3) Atlantic Ocean, (4) Mediterranean Sea, and (5) Pacific Ocean. We use a model configuration with 96 grid points in longitude, 72 in latitude and 19 vertical layers, with the first four layers in the first kilometre above the surface. LMDZ-iso has a stretchable grid that allows increased spatial resolution over a defined region. In our case, it gives an averaged resolution of ∼ 100 km over central Asia, which is a good trade-off between a reasonable computing time and a spatial resolution that adequately represents main features of TP topography.
Here we report results from three experiments designed to isolate the influence of Asian topography on climate and isotopic composition of precipitation. Topography is derived from a 10 min US Navy dataset and interpolated to the model grid. The control run (MOD) is a pre-industrial run, i.e. initialised with boundary conditions (insolation, greenhouse gases, sea surface temperatures (SSTs), topography) kept at pre-industrial values. For the two other experiments, we keep all boundary conditions (including albedo, rugosity, and vegetation distribution) similar to those in MOD run, except for the topography. We reduce the altitude over the area covering the Tibetan Plateau, Himalayas and a part of surrounding mountains: Tian Shan, Pamir, Kunlun and Hindu Kush to 50 % of modern elevations (intermediate, INT case) and to 250 m elevation (low, LOW case) (Fig. 1). SSTs for all runs come from the AMIP dataset (monthly SSTs averaged from 1979 to 1996; Taylor et al., 2000). Each experiment has been run for 20 years. We analyse seasonal means over the last 18 years, as the two first years are extracted for spin-up.

Theoretical framework for the precipitation composition
Our goal is to understand to what extent topography changes explain the precipitation δ 18 O signal over TP (i.e. the direct topography effect) and what part of this signal depends on other climate processes. To do so, we develop a theoretical expression for the precipitation composition.
To the first order, the δ 18 O composition of the precipitation R p follows that of the vapour R v . Deviations from the vapour composition, ε = R p −R v , are associated with a local

Initial site
Site of interest condensational or post-condensational process.
In an idealized framework of an isolated air parcel transported from an initial site at low altitude to the site of interest ( Fig. 2), the vapour composition can be predicted by Rayleigh distillation: where R v0 is the initial composition of the vapour at the initial site; α is the fractionation coefficient that depends on temperature and on the water phase (Majoube, 1971;Merlivat and Nief, 1967); and f is the residual fraction of the vapour at the site of interest relatively to the initial site of an air mass ascent. We take the initial site as characterised by a temperature and humidity T 0 and q 0 . Under these conditions, R v0 is the theoretical isotopic composition of vapour that it would have if all the vapour originated from the local evaporation over quiescent oceanic conditions. Depending on the atmospheric circulation, on deep convective and mixing processes and on the source region of water vapour, the isotopic composition of vapour may deviate from the Rayleigh distillation by dR v . The residual fractionf depends on the specific humidity q at the site of interest: The air is not always saturated near the surface, therefore: where h and T s are the relative humidity and air temperature near the surface of the site of interest. The air can be under-saturated because it can be considered as air that has been transported adiabatically from the area of minimum condensation temperature, T * (Galewsky and Hurley, 2010;Galewsky et al., 2005;Sherwood, 1996): q = q s (T * ).
The surface temperature can be predicted to first order by the adiabatic lapse rate, , and is modulated by the nonadiabatic component dT s that represents processes such as large-scale circulation or radiation: where z and z 0 are the altitudes at the site of interest and at the initial site. We use an adiabatic lapse rate equal to −5 km −1 based on the measurements of modern observed mean temperature lapse rate on the southern slope of the central Himalayas, that ranges from −4.7 to −6.1 km −1 for the monsoon season and from −4.3 to −5.5 km −1 for the rest of the year (Kattel et al., 2015). If we combine Eqs.
(1) to (5), we get that R v is a function of ε, dR v , h, dT s and z: Or in a simpler form: Parameters z 0 , q 0 , T 0 are reference values that are common to all sites of interest, all climates and geographies. Even if initial conditions for the Rayleigh distillation vary depending on the atmosphere circulation, on deep convective processes and on the site of interest, we keep the same reference values and we consider all variations in initial conditions are accommodated by dR v . This model is equivalent to that of Rowley et al. (2001) for dR v = 0 (i.e. neglecting the effects of mixing and deep convection on the initial water vapour), ε = (α − 1) ×R v (i.e. neglecting post-condensational effects), and h = 1 (i.e. assuming the site of interest is inside the precipitating cloud).

Decomposing precipitation composition differences
Our goal is to understand why R p varies from one climatic state to another. We refer to these climatic states using subscript 1 and 2 and to their difference using the notation. Differences between INT and LOW and between MOD and INT climatic states correspond to the initial and the terminate stages of the TP uplift respectively. We decompose R p = R p2 − R p1 into contribution from dR v , ε, h, dT s , and z: where R p, ε , R p, dR v , R p, h , R p, dT s , and R p, z are the contributions of dR v , ε, h, dT s , and z to R p . Non linear terms of decomposition are gathered into the residual term N . Contributions are estimated using Eq. (7) (see also Table 1): In order to decrease the sensitivity of the decomposition to the state at which it has been calculated we take z , dT s , h , dR v , and ε as centred differences: Note that ε in Eqs (10) to (13) and dR v in Eqs. (9) and (11) to (13) can be replaced by 0 without changing the result. Parameters z, dT s , h, dR v , and ε are diagnosed for the climatic states 1 and 2 from LMDZ-iso simulations (e.g. for pairs of experiments, MOD and INT cases). Parameter ε is estimated as ε = R p −R v , where R p and R v are isotopic ratios simulated by LMDZ-iso. Parameter h is the relative humidity simulated by LMDZ-iso. Altitude z is a prescribed boundary condition of the simulations. Parameter dR v is estimated by calculating the difference between the water vapour isotopic ratio simulated by LMDZ-iso (R v,LMDZ ) and that predicted by Rayleigh distillation if the initial water vapour isotopic ratio is R v0 : where q is the specific humidity simulated by LMDZ-iso and α is the isotopic fractionation as a function of the nearsurface air temperature T s simulated by LMDZ-iso. Parameter dT s is estimated from Eq. (5) by calculating the difference between the near-surface air temperature simulated by LMDZ-iso and that predicted by the adiabatic lapse rate: All the isotopic decomposition terms computed are weighted by the precipitation amount.
Total isotopic difference between state 2 and state 1 Effect of lapse rate change, associated with non-adiabatic effects, possibly due to changes in surface energy budget or in large-scale atmospheric stratification Effect of local relative humidity change, possibly due to large-scale circulation changes Effect of changes in condensational and post-condensational effects, possibly due to changes in rain reevaporation processes All other effects, including effects of deep convection, mixing, water vapour origin, continental recycling on the initial water vapour

Robustness of the decomposition
First, to check whether the linear decomposition is a good approximation of the total R p change, we estimate the nonlinear term N as a residual, i.e. for each pair of states, we calculate the deviation of R p = R p (ε 2 , dR v2 , h 2 , dT s2 , z 2 ) − R p (ε 1 , dR v1 , h 1 , dT s1 , z 1 ) from LMDZ-simulated isotopic differences between the two experiments. N represents less than 17 % of the total R p change for both stages of TP uplift. Our method to estimate the terms in Eq. (7) is equivalent to first order approximation of partial derivatives, i.e. we neglect the sensitivity of the partial derivatives to the state at which they are calculated. We tested this sensitivity by using Eqs. (9) to (13) changing z to z 1 or z 2 , dT s to dT s2 or dT s1 and so on. For example, in Eq. (13), replacing h with h 1 changes the resulting R p, z by 0.03 ‰, replacing h with h 2 has an impact of 0.09 ‰. In the same equation, replacing dT s with dT s1 and with dT s2 contributes to R p, z with 0.005 and 0.039 ‰ respectively. As it was highlighted earlier, replacing ε and dT s with ε 1 or ε 2 and dR v1 or dR v2 respectively has no impact to the resulting R p, z Thus, our method shows low sensitivity to the state.
Second, to check the influence of initial conditions R v0 , T 0 and q 0 on the decomposition, we estimate the sensitivity of the different contributions to changes in R v0 , T 0 , and q 0 , of 1 %, 1 K and 10 % respectively (Table 2). R v0 is the parameter that influences most of the decomposition terms, with a maximal sensitivity of 0.9 ‰ obtained for R p, z for a change of 1 ‰ in R v0 . Sensitivity to temperature and humidity is lower, ranging from 0 to 0.6 ‰. Overall, all the decomposition terms show a sensitivity < 1 ‰ with most (82 %) of them < 0.5 ‰, making our decomposition method robust.

Model validation in terms of simulated climate variables
LMDZ has been used for numerous present-day climate and palaeoclimate studies (Kageyama et al., 2005;Ladant et al., 2014;Sepulchre et al., 2006), including studies of monsoon region (e.g. Lee et al., 2012;Licht et al., 2014). Yao et al. (2013) showed that LMDZ-iso has the best representation of the altitudinal effect compared to similar GCM and RCM models. These authors have also provided a detailed description of rainfall patterns over the Tibetan Plateau, and showed LMDZ-iso ability to simulate atmospheric dynamics and reproduce rainfall and δ 18 O patterns consistent with data over this region. For the purpose of our experiments validation, we compare MOD experiment outputs with rainfall data from the Climate Research Unit (CRU) (New et al., 2002) ( Fig. 3a, b, c). When compared to CRU dataset, MOD annual rainfalls depict an overestimation over the high topography of the Himalayas and the southern edge of the Plateau, with a rainy season that starts too early and ends too late in the year. Over central Tibet (30-35 • N), the seasonal cycle is well captured by LMDz-iso, although monthly rainfall is Table 2. INT-LOW and MOD-INT sensitivity of the decomposition terms (in ‰) to the changes of R v0 , T 0 , q 0 , of 1 ‰, 1 K and 10 % respectively.

Northern region Southern region
always slightly overestimated (+0.5 mm day −1 ). CRU data show that the northern TP (35-40 • N) is drier with no marked rainfall season and a mean rainfall rate of 0.5 mm day −1 . In MOD experiment, this rate is overestimated (1.5 mm day −1 on annual average). Despite these model data mismatches, the ability of LMDZ-iso to represent the seasonal cycle in the south and the rainfall latitudinal gradient over the TP allows its use for the purpose of this study. Our MOD simulation is pre-industrial, consequently a comparison with modern data is expected to provide differences driven by the pre-industrial boundary conditions. Still comparing LMDZ-iso outputs with mean annual temperatures from CRU dataset (New et al., 2002) (Fig. 3d, e, f) and relative humidity from NCEP-DOE Reanalysis (Kanamitsu et al., 2002) (Fig. S1 in Supplement) shows that LMDZ-iso model captures these variables reasonably well.

Impact of TP uplift on Asian climate
Theoretically, the Tibetan Plateau has both mechanical and thermal effects on atmospheric dynamics that induce increased monsoon activity to the south and drive arid climate to the north (Broccoli and Manabe, 1992;Sato and Kimura, 2005). Thus, modifying TP height is expected to alter these large-scale atmospheric dynamics and associated climate variables (namely temperature, precipitation, relative humidity (hereafter RH), cloud cover), and in turn to affect the isotopic signature of rainfall.
In LOW experiment, strong summer heating leads to the onset of a "Thermal Low" at the latitude of maximal insolation (ca. 32 • N), similar to the present-day structure existing over the Sahara (Fig. S2). This structure is superimposed by large-scale subsidence linked to the descending branch of the Hadley cell, and both factors act to drive widespread aridity over the TP area between ca. 30 and 40 • N, associated with very low (< 40 %) RH values (Fig. S2). Subsidence also prevents the development of South Asian monsoon over the north Indian plane and favours aridity over this region. In winter, large-scale subsidence induces high surface pressures and creates an anticyclonic cell that prevents convection and humidity advection, resulting in low RH and annual rainfall amount ranging from 50 to 500 mm over the TP area (Fig. 4f).
Uplifting TP from 250 m above sea level (ASL) to half of its present-day altitude (INT case) initiates convection in the first tropospheric layers, restricting large-scale subsidence to the upper levels (Fig. 4c, e). In turn, South Asian monsoon is strengthened and associated northward moisture transport and precipitation increase south of TP (Figs. 5, 6). As a consequence the hydrological cycle over TP is more active, with higher evaporation rates (Fig. 7d). Together with colder temperatures linked to higher altitude (adiabatic effect) (Fig. 7b), the stronger hydrological cycle drives an increase in RH (Fig. 7a) and cloud cover (Fig. S3). Another consequence of increased altitude is higher snowfall rates in winter and associated rise of surface albedo (Fig. S4). When added to the increased cloud cover effect, this last process contributes to an extra cooling of air masses over the Plateau. To the north of TP, the initial stage of uplift results in increased aridity (i.e. lower RH and rainfall) over the Tarim Basin region (Fig. 6). This pattern can be explained both by a barrier effect of southern topography and by stationary waves strengthening, which results in subsidence to the north of TP. This latter mechanism is consistent with pioneer studies which showed that mountain-related activation of stationary waves prevented cyclonic activity over Central Asia and induced aridity over this region (Broccoli and Manabe, 1992).
The impact of the terminal stage of TP uplift also drives an increase in RH over the Plateau, especially during summer time, when a very active continental recycling (  amount also increases significantly (Fig. 6), driven both by increased evaporation and water recycling during summer, and intense snowfall during winter. The latter contributes to the increase in the surface albedo and associated surface cooling during winter. Conversely, the uplift to a modernlike Plateau reduces RH (down to 30 %) north of the Plateau,  and allows the onset of large arid areas. We infer that this aridification is linked to a mechanical blocking of moisture transport, both by Tian Shan topography for the winter westerlies, and the eastern flanks of TP for summer fluxes; despite changes in stationary waves structure and sensible heat (not shown), no marked shift in subsidence between INT and MOD experiments is simulated. This result is consistent with recent studies (Miao et al., 2012;Sun et al., 2009) that have suggested the potential contribution of Pamir and Tian Shan rainshadow effect to aridification in Qaidam Bassin and the creation of Taklamakan Desert.  (Fig. 8a, b). In general, the model shows a good agreement with precipitation and VSMOW-weighted modern surface waters δ 18 O, including stream, lake and spring waters (data from Bershaw et al., 2012;Hren et al., 2009;Quade et al., 2011), as testified by a Pearson coefficient of 0.86 between modelled and observed precipitation δ 18 O (Fi. 8c). This comparison shows the ability of LMDZ-iso to reproduce the decrease in δ 18 O from the Indian subcontinent to the Himalayan foothills and with minimum values over the Himalayas. Simulated increase in δ 18 O over the TP with the distance from the Himalayas is also consistent with data sampled along a southwest-northeast transect across the Plateau (Bershaw et al., 2012). However, over the northern margins of the TP, LMDZ-iso underestimates simulated δ 18 O in precipitation (Fig. 8a). This model data mismatch may result from two types of uncertainties. First despite the high resolution obtained with a zoomed grid, restricted topographic features could not be well-captured over some parts of the TP, which could lead our simulations to miss local processes affecting δ 18 O in rainfall. Second, overestimating the westerlies fluxes (see the comparison with the ERA moisture transport on Fig. 5a) could lead to underestimate δ 18 O over the northern part of the TP, through advection of depleted air masses. Nevertheless, despite our model not capturing the absolute maximal values well, the regional latitudinal gradient is correctly represented, and most observed values are within the range of simulated δ 18 O (Fig. 8b). We consider that the ability of LMDZ-iso to represent this gradient makes it reliable to carry out this study, which is focusing on sensitivity experiments with large changes in topography and associated anomalies in δ 18 O.

Simulated isotopic changes and signal decomposition
To first order, increasing topography over TP leads to more negative δ 18 O over the region (Fig. 9). In the absence of topography, precipitation δ 18 O follows a zonal pattern and undergoes a weak latitudinal depletion on the way to the continental interior, except from slight deviations over India, central China and the eastern part of the TP (Fig. 9b) (Fig. 9a). The total difference in isotopic composition of precipitation, R p , between pairs of experiments (INT-LOW, MOD-INT) is significant beyond the areas where the topography was reduced by the experimental design (Figs. 10a, 11a). Substantial differences in δ 18 O between MOD and INT experiments are simulated over the southern TP (up to 10 ‰) www.clim-past.net/12/1401/2016/  and over the Tarim Basin (up to 7 ‰). Between INT and LOW cases, the differences are over the margins of the TP, over Pamir, Tian Shan and Nan Chan. We should note that the isotopic difference becomes more important for the later stage of the plateau uplift. For clarity, we define two boxes, over the northern (from 34 to 38 • N and from 88 to 100 • E) and southern (from 27 to 33 • N and from 75 to 95 • E) part of TP (Fig. 12).

Direct topography effect on δ 18 O
The direct effect of topography change is determined as the decomposition term R p, z in Eq. (7). For the initial stage of the uplift, the altitude effect produces a decrease in precipitation δ 18 O ranging from −1 to −3 ‰ (Fig. 10b). For the terminal stage of the uplift, the isotopic decrease linked with altitude goes up to −7 ‰ (Fig. 11b). Differences between both stages are linked to the non-linear relationship between δ 18 O and elevation. Also for both stages, the difference between R p and R p, z is non-zero (Figs. 12a, 13a). These differences are particularly marked for the terminal stage, for which R p, z averages −5.5 ‰ over the northern part of TP (Fig. 13a, b), whereas the total isotopic change averages −3 ‰. Locally, the difference between R p, z and R p can reach +4 ‰. When averaged over the southern box, R p, z is less negative (−4 ‰) than R p (−4.6 ‰), with localised maximum differences reaching −4 ‰ (Table 3). Offsets between R p, z and R p are also detected for the initial stage of the uplift (Fig. 12a, b), but are lower: they reach +2 ‰ over central TP but barely reach 1 ‰ when averaged over southern and northern boxes. These offsets are related to additional effects of uplift on δ 18 O that are discussed in the following sections.  Table 1). Figure 11. (a) Total isotopic difference between MOD and INT experiments ( R p ) and spatial isotopic variations related to the following: (b) direct effect of topography changes, (c) effect of lapse rate change, associated with non-adiabatic effects, (d) effect of local relative humidity change, (e) effect of changes in post-condensational processes, (f) all other effects (see Table 1).

Non-adiabatic temperature changes impact
Besides the adiabatic temperature effects linked with the TP uplift, non-adiabatic temperature changes can be identified, in relation with surface albedo and cloud cover changes depicted in Sect. 3.2.1. The term R p, dT s in Eq. (7) ( Table 1, line 3) is associated with these non-adiabatic effects, i.e. spatial variations of the temperature lapse rate. Figures 10c and  11c show the portion of the total isotopic signal that is linked to this effect. It plays a modest role for the early phase of uplift (+1-2 ‰ locally), but is more important for the second stage. It contributes to 2-5 ‰ of total isotopic difference, with a positive sign over southeast TP interior, TP northern margins and Asia interior. Negative anomalies have a magnitude of 2-3 ‰, but are less widespread, localised over the TP interior (Fig. 11c). Positive isotopic anomalies are associated with a steeper lapse rate than expected based on adiabatic processes. Conversely, negative δ 18 O anomalies that are observed over northern TP and over Pamir are explained by a weaker lapse rate than adiabatic. Overall, these variations represent between 10 and 19 % (4-10 % for the initial stage) of the processes that are not linked to topography (Figs. 12d, e and 13d, e).

Impact of RH changes during condensation process
The term R p, h in Eq. (7) depicts the portion of total isotopic signal R p linked to local RH change during condensation process (Table 1, line 4). Over TP, R p, h is positive for both uplift phases, and RH changes act as a counterbalance to the topography effect. R p, h reaches +4 ‰ for the late stage (Fig. 11d), and maxima are located over the western part and northern part of TP for both stages of the uplift. Equation (4) shows that this positive anomaly is directly related to the increase in RH described in Sect. 3.2.1. For the initial stage, R p, h depicts also positive values (up to +3 ‰) to the southwest of TP. When averaged over northern and southern boxes, the counterbalancing effect of RH on R p ranges from 1.5 to +3 ‰, and this effect represents up to 76 % of all non-topographic processes (Figs. 12, 13). Interestingly, an opposite signal is simulated over the Tarim basin, where topography was kept constant in the three experiments. This signal is consistent with the previously-depicted decrease in RH over this region, in relation with rain-shadow effects and large-scale subsidence.

Post-condensation processes impact
The difference between δ 18 O v and δ 18 O p is linked to the post-condensation effects, mainly associated with raindrop reevaporation that can occur after initial condensation. Because lighter isotopes evaporate more easily, rain reevaporation leads to an isotopic enrichment of precipitation. Therefore, the more reevaporation, the greater the difference between δ 18 O p and δ 18 O v . We refer to the study of (Lee and Fung, 2008), where post-condensation effects are explained in detail. The contribution of such processes increases dramatically for very dry areas, where the relative humidity is less than 40 %. Estimation of term R p, ε , i.e. the change in isotopic difference between vapour and precipitation, allows to quantify the contribution of post-condensational processes to total R p signal (Figs. 10e, 11e) without appealing to the d-excess. For both stages of uplift, R p, ε is mostly negative, indicating a depletion of R p relative to R v with the uplift. Over the Plateau, contribution of post-condensational effects for the initial stage of uplift ranges from 25 % to 46 % of total non-topographic effects, whereas it represents less than 10 % for the terminal stage (Figs. 12a, 13a). The most significant signal is simulated over the northern part of the Plateau and over its western margin and adjacent areas. Post-condensational effects during the initial stage lead to up to a −5 ‰ anomaly over the western margin of TP (Fig. 12e) whereas the terminal stage creates a substantial negative anomaly only over northern TP margin and Tarim Basin (Fig. 13e).

Residual processes effect
The last term of Eq. (7), R p, dR v , corresponds to the part of the total isotopic signal that could not be explained by previously mentioned processes. These residual anomalies are rather weak for the initial stage of the uplift, explaining less than 1 ‰ of the signal over the northern plateau, and around 1 ‰ over the southern TP and adjacent parts of Asia and India (Fig. 10f). Contribution of these effects to the initial stage is 4 and 21 % to the northern and southern box respectively (Fig. 12d, e). Conversely, for the terminal stage of the TP uplift this anomaly reaches up to −4 ‰ over the southern part of the TP (Fig. 11f) and contributes to 49 % of the nontopographic processes signal (Fig. 13d, e). In the next sections we propose several mechanisms that could contribute to this residual anomaly.

Discussion
Our results suggest that TP uplift affects precipitation δ 18 O through direct topographic effect, but that a significant part of the signal is related to several other processes. These processes alter the isotopic signal not only over TP, but also over adjacent regions, where topography was kept the same by the experiment design. A second result is that despite a similar altitudinal change of TP between the two uplift stages, the topographic effect on δ 18 O is more perturbed by other processes during the terminal stage than during the initial one.
For the terminal stage, the residual effects change over the southern region dominates (49 %) the isotopic signal that is not linked to the direct topographic effect. The RH change and non-adiabatic temperature changes also have an important counterbalancing impact, together contributing to 43 % of the isotopic signal (Fig. 13e). For the northern region, the topographic effect is mainly counterbalanced by the RH change effect (2.5 ‰), ultimately leading to a 2.3 ‰ offset between R p and what is expected from topography. Here RH contributes to 76 % of the isotopic signal not linked with the topography change, while nonadiabatic temperature changes, residual effects change and post-condensational processes have an impact of 16, 7 and < 1 % respectively (Fig. 13d).

Impact of RH variations
RH alters rainfall isotopic signature through two steps, during and after condensation. As mentioned earlier, the first effect of RH, as shown in Eq. (4) and expressed as R p, h , occurs during condensation through Rayleigh distillation and induces that R p increases with increasing RH. Our model shows that RH increases over TP with the initial stage of uplift, driving precipitation δ 18 O towards less negative values. This mechanism is more efficient for the terminal stage of uplift, when RH is increased in summer as a response of a more active water cycle. South of TP, RH direct effect on δ 18 O is noticeable, as efficient moisture transport is activated with the uplift-driven strengthening of monsoon circulation (Fig. 4). Interestingly, this mechanism is not active for the second stage of the uplift, during which rainfall increases through more effective convection, not through higher advection of moisture. As a consequence, negligible RH and R p changes are simulated south of the Plateau when it reaches its full height. This suggests that an altitudinal threshold might trigger South Asian monsoon strengthening, and ultimately precipitation δ 18 O signature, a hypothesis that should be explored in further studies. Conversely, the negative values of www.clim-past.net/12/1401/2016/ Clim. Past, 12, 1401-1420, 2016 R p, h over and northeast of the Tarim basin are related to a decrease in RH during both stages. Our analysis suggests that the first uplift stage is sufficient to create both barrier effects to moisture fluxes and large-scale subsidence that ultimately drive aridity over the region.
The second effect of RH on δ 18 O concerns very dry areas (ca. < 40 %), where raindrop re-evaporation can occur after initial condensation, leading to an isotopic enrichment of precipitation compared to water vapour (Lee and Fung, 2008) (Fig. S2). Such an effect is implicitly included in the postcondensational term of our decomposition that shows opposite sign when compared to R p, h . Over the Plateau, this mechanism is effective only for the first uplift stage, where TP area transits from very low precipitation amounts and very low RH values to wetter conditions (Fig. S7).
Over TP, the opposed effects of RH almost compensate each other for the early stage of the uplift (Fig. 10d, e), but this is not the case for the final stage, since RH postcondensational effect is similar between INT and MOD experiments. Since absolute values of the impact of RH through condensation and post-condensational processes can reach 5 ‰, it is crucial to consider RH variation when inferring palaeoaltitudes from carbonates δ 18 O.

"Amount effect" and monsoon intensification
Our results also show a substantial increase in precipitation amount over northern India, the Himalayas and TP with the growth of topography for both uplift stages (Fig. 14). The inverse relation between the enrichment in heavy isotopes in precipitation and precipitation amount, named the "amount effect" (Dansgaard, 1964) is largely known for oceanic tropical conditions (Risi et al., 2008;Rozanski et al., 1993) and for Asia monsoonal areas (Lee et al., 2012;Yang et al., 2011). Over South Tibet recent studies have shown the role of deep convection in isotopic depletion (He et al., 2015). For the two stages of uplift, the residual component of the isotopic signal depicts negative values over southern TP, where annual rainfall amount is increased. Thus we infer that this anomaly can be driven, at least partly, by the amount effect that increases with growing topography.
Various climate studies have suggested that the appearance of the monsoonal system in East Asia and the onset of central Asian desertification were related to Cenozoic Himalayan-Tibetan uplift and withdrawal of the Paratethys Sea Clift et al., 2008;Guo et al., 2002Guo et al., , 2008Kutzbach et al., 1989Kutzbach et al., , 1993Ramstein et al., 1997;Raymo and Ruddiman, 1992;Ruddiman and Kutzbach, 1989;Sun and Wang, 2005;Zhang et al., 2007), although the exact timing of the monsoon onset and its intensification remains debated (Licht et al., 2014;Molnar et al., 2010). Although our experimental setup, which does not include Cenozoic palaeogeography, was not designed to assess the question of monsoon driving mechanisms nor its timing, our results suggest that uplifting the Plateau from 250 m ASL to half of its present height is enough to enhance moisture transport towards northern India and strengthen seasonal rainfall. Nevertheless, massive increase of rainfall over TP between INT and MOD experiments indicates that the second phase of uplift might be crucial to activate an efficient, modern-day-like, hydrological cycle over the Plateau. The decrease in simulated precipitation north of the Plateau also suggests that terminal phase of TP uplift triggered modern-day arid areas.

Other effects
Although precipitation amount change explains well the residual isotopic anomaly (Figs. 10f, 11f), additional processes could interplay. Continental recycling can overprint original moisture signature and shifts the isotopic ratios to higher values due to recharging of moisture by heavy isotopes from soil evaporation (Lee et al., 2012;Risi et al., 2013). In our simulation, we detect an increasing role of continental recycling in the hydrological budget of the TP (Fig. S6), especially in its central part, that likely shifts the δ 18 O to more positive values and partially compensates for the depletion linked to the "amount effect" over the central plateau. Another process frequently invoked to explain the evolution of precipitation δ 18 O patterns over TP is changes in moisture sources (Bershaw et al., 2012;Dettman et al., 2003;Quade et al., 2007;Tian et al., 2007). Except for the continentally recycled moisture, southern Himalayas precipitation moisture originates mainly from the Indian, the Atlantic and the Pacific oceans (Fig. S6). Proximate oceanic basins are known to be sources of moisture with a more positive signature than remote ones (Chen et al., 2012;Gat, 1996). Supplemental analyses with water-tagging feature of LMDZiso show that contribution of continental recycling to rainfall over TP increases with the uplift, at the expanse of Pacific and Indian sources (Fig. S6). Although we have no mean to decipher between sources and amount effect in the residual anomaly, it seems that the change of sources is not sufficient to yield a strong offset of δ 18 O values.

Relevance of palaeoelevation reconstructions based on palaeo δ 18 O
Quantitative palaeoelevation reconstructions using modern altitude-δ 18 O relationship will succeed only if R p corresponds mainly to the direct topography effect. Modern palaeoaltimetry studies cover almost all regions of the Plateau for time periods ranging from Palaeocene to Pleistocene-Quaternary (see data compilation in Caves et al., 2015). Most of these studies consider changes in δ 18 O as a direct effect of the topography uplift. Palaeoelevation studies locations (see Caves et al., 2015 for a synthesis) plotted over the anomaly maps (Figs. 12a, 13a) show for what geographical regions restored elevations should be used with an additional caution. Numerous palaeoelevation data points were located either over the northern part of the TP (from 34  Saylor et al. (2009) to 38 • N and from 88 to 100 • E) or over the southern region (from 27 to 33 • N and from 75 to 95 • E). Our model results show that when TP altitude is increased from half to full, considering topography as an exclusive controlling factor of precipitation δ 18 O over the southern (northern) region likely yields overestimations (underestimations) of surface uplift, since the topography effect is offset by RH and amount effects. Projecting our modelling results to each locality where palaeoelevation studies have been published (Table 4) reveals that topography change explains simulated total isotopic change reasonably well for only few locations (Linzhou Basin, Lunpola Basin, Kailas Basin, Huaitoutala). Indeed topography appears to be the main controlling factor for only 40 % of the sites, while 30 % is dominated by RH effects, 20 % by residual effects and 5 % by postcondensational and non-adiabatic temperature changes, respectively. Nevertheless such figures have to be taken carefully, since we ran idealized experiments testing only the impact of uplift, neglecting other factors like horizontal palaeogeography or pCO 2 variations, the latter being known to influence δ 18 O as well (Jeffery et al., 2012;Poulsen and Jeffery, 2011).
For the initial uplift stage apparent consistency occurs between the topography impact and the total isotopic composition, in relation to counteracting effects RH and postcondensational processes. For the southern region RH impact appears to be the main controlling factor for the isotopic composition of precipitation, surpassing the direct to-pography impact. Nevertheless, these processes have a different contribution for initial and terminal stages of uplift. Precipitation changes lead to overestimating altitude changes for both stages, but for the terminal stage its contribution is bigger. This effect dominates in the southern part, and more generally where the isotopic composition of precipitation strongly depends on convective activity. RH changes dominate over the western part of TP and Northern India for initial uplift stage and over the northern TP for the terminal. Differences between both stages could be partly explained by non-linearities in q s -temperature relationships, as well as in Rayleigh distillation processes (Fig. S8). Determining whether other processes contribute to this difference would be of interest, but it is out of the scope of the present study.

Conclusions
Previous studies focusing on the Andes (Ehlers and Poulsen, 2009;Poulsen et al., 2010) or North American cordillera (Sewall and Fricke, 2013) have inferred that the impact of uplift of mountain ranges on δ 18 O could be altered by the consequences of the uplift on atmospheric physics and dynamics. Our modelling results show that it is also the case for the Tibetan Plateau uplift. Additionally, we designed a decomposing analysis to quantify for the first time the different processes that can alter precipitation δ 18 O changes with uplift. As suggested for the Andes, the onset of convective rainfall plays an important role in shifting δ 18 O towards more negative values. Nevertheless this process is not the main factor, as we show that saturation of air masses, quantified by RH have two to three-time bigger effects on the final δ 18 O. We infer that increase in precipitation linked with the TP uplift would lead to overestimation of the topography uplift at sites over Himalayas and Southern TP, whereas increase in RH leads to underestimating the uplift at sites in Northern Tibet.
Our results could be applied to interpret palaeoclimate records and to reconstruct the region uplift history. Palaeoelevation reconstructions suggest the Himalayas attained their current elevation at least by the late Miocene or even earlier (Garzione et al., 2000a, b;Rowley et al., 2001;Saylor et al., 2009). Our results show overestimation of the topography impact over this region, thus the Himalayas may have attained their current elevation later than expected. In contrast, isotope-based palaeoaltimetry could underestimate surface elevation over the northern TP. This could explain why available isotope-based palaeoelevation estimates for the northern TP (Cyr et al., 2005), which estimates surface elevation at about 2 km, contradict palynological assemblages in lacustrine sediments from the Xining Basin, which show the presence of high-altitude vegetation at the same time period (Dupont-Nivet et al., 2008;Hoorn et al., 2012).
Still, our decomposition methods reveal that even if the impact of the TP uplift phases are rather straightforward (monsoon enhancement to the South, increase in continen-tal recycling over TP, moisture fluxes deflection and increased aridity to the North), the consequences in terms of δ 18 O are extremely complex, since interplays and compensation occur amongst all the processes. Limitations in our approach are related to a perfectible hydrological cycle in LMDZ-iso, and idealized boundary conditions (topography uplift scenarios, modern land-sea mask, SSTs and pCO 2 ). Model data comparison show that mean annual precipitation amount is slightly overestimated by the model for the northern TP, thus could result in underestimation of the amount effect contribution for the northern TP. On the contrary, the model overestimates the precipitation over the southern edge of the Himalayas. If it was more realistic, the contribution of the amount effect estimated by the decomposing method could be less important. Changes in vegetation cover, by altering albedo and persistence of snow cover, could affect the impact of non-adiabatic temperature changes on δ 18 O. Vegetation over Asia was shown to have a major variation through Cenozoic based on pollen (Dupont-Nivet et al., 2008;Miao et al., 2011;Song et al., 2010;Zhao and Yu, 2012) and palaeobolanical data De Franceschi et al., 2008;Kohn, 2010) and future studies would benefit to explore its impact on precipitation δ 18 O. Also it is widely known that during the Cenozoic air temperature was higher due to a higher concentration of greenhouse gases in the atmosphere (Zachos et al., 2008). Studies taking into account this feedback inferred that it could lead to even larger inaccuracy in surface uplift estimations during the Cenozoic (Poulsen and Jeffery, 2011). Thus the field of palaeoaltimetry would benefit from future studies focusing on (1) using palaeoclimate proxies to constrain specifically relative humidity, surface temperature and precipitation amount in deep time and (2) applying a decomposition method to isotope-enabled GCM simulations forced by constrained palaeogeography (land-sea mask and different scenarios for orogens) and atmospheric pCO 2 for specific geological time period. The combination of both could help refine calibration for palaeo δ 18 O-elevation relationships and refining palaeoelevation estimates.
The Supplement related to this article is available online at doi:10.5194/cp-12-1401-2016-supplement.