Precipitation d18O on the Himalaya-Tibet orogeny and its relationship to surface elevation

Abstract. The elevation history of the Himalaya–Tibet orogen is central to
understanding the evolution and dynamics of both the India–Asia collision
and the Asian monsoons. The surface elevation history of the region is
largely deduced from stable isotope (δ18O, δD)
paleoaltimetry. This method is based on the observed relationship between the
isotopic composition of meteoric waters (δ18Op,
δDp) and surface elevation, and the assumption that
precipitation undergoes Rayleigh distillation under forced ascent. Here we
evaluate how elevation-induced climate change influences the
δ18Op–elevation relationship and whether Rayleigh
distillation is the dominant process affecting
δ18Op. We use an isotope-enabled climate model,
ECHAM-wiso, to show that the Rayleigh distillation process is only dominant
in the monsoonal regions of the Himalayas when the mountains are high. When
the orogen is lowered, local surface recycling and convective processes
become important, as forced ascent is weakened due to weaker Asian monsoons.
As a result, the δ18Op lapse rate in the Himalayas
increases from around −3 to above −0.1 ‰ km−1, and has little
relationship with elevation. On the Tibetan Plateau, the meridional gradient
of δ18O decreases from ∼1 to
∼0.3 ‰ ∘−1 with reduced elevation, primarily due to
enhanced sub-cloud reevaporation under lower relative humidity. Overall, we
report that using δ18Op or
δDp to deduce surface elevation change in the
Himalayan–Tibetan region has severe limitations and demonstrate that the
processes that control annual-mean precipitation-weighted
δ18Op vary by region and with surface elevation.
In summary, we determine that the application of
δ18O paleoaltimetry is only appropriate for 7 of the 50 sites
from which δ18O records have been used to infer past elevations.



Introduction
Surface elevation is a fundamental characteristic of the Earth surface, directly affecting atmospheric circulation patterns and surface temperatures, precipitation and surface hydrology, erosion and sediment transport, and the distribution and diversity of life (Aron and Poulsen, 2018). The evolution of surface elevation resulting from India-Asia convergence and creation of the Himalaya-Tibet Orogen was the defining event of the Asian continent in the Cenozoic and had global environmental 20 implications. Surface uplift of the orogen has been implicated in the onset and strengthening of the southeast Asian monsoon system (Boos and Kuang, 2010;Zhang et al., 2015) and the position of atmospheric stationary waves (Kutzbach et al., 1989); the evolutionary diversification and biogeographic distribution of fauna and flora in central Asia (Zhao et al., 2016;Yang et al., 2009;Antonelli et al., 2018); and, the intensification of chemical weathering of exposed rocks and transport of nutrients to the ocean that contributed to global atmospheric CO2 drawdown (Galy et al., 2007;Maffre et al., 2018). Surface elevation 25 also lends a first-order constraint on the crustal and upper mantle dynamics that create topography ). Surface elevation estimates for the Himalaya-Tibet orogen, and specifically evidence for high elevations of the orogen since the late Eocene, have been instrumental in supporting geodynamical models of Tibetan Plateau growth through early deformation and crustal thickening (Rowley and Currie, 2006;Rohrmann et al., 2012;Hoke et al., 2014).
The paramount importance of surface elevation to our understanding of Cenozoic environmental and tectonic evolution has 30 led to a proliferation of studies that infer past Himalaya-Tibet surface elevations from ancient proxy materials (Cyr et al., 2005;Rowley and Currie, 2006;Li et al., 2015). Stable isotope paleoaltimetry, one of the few quantitative methods to reconstruct past surface elevations, relies on the water isotopic composition of ancient materials, including pedogenic and lacustrine carbonates, authigenic clay and hydrated volcanic glass, that were formed in contact with ancient surface waters.
The method is predicated on the observed, modern decrease in water stable isotopic compositions (d 18 O, dD) with elevation gain (the isotopic lapse rate) (Chamberlain and Poage, 2000) a relationship that is commonly attributed to rainout of heavy 5 isotopologues during stably forced ascent of a saturated air over high elevation and is modeled as a Rayleigh distillation process (e.g. Rowley and Garzione, 2007). Stable isotope paleoaltimetry has been used to reconstruct past surface elevation of many of the world's major mountain belts, including the North America Cordillera (e.g. Poage and Chamberlain, 2002;Fan et al., 2014), the Andes (e.g. Garzione et al., 2008), and the Himalayan-Tibetan orogen (e.g. Rowley and Currie, 2006), due to the robustness of the isotopic lapse rate in modern orogenic regions and the ubiquity of proxy materials. 10 The interpretation of stable isotopes in ancient materials to infer past surface elevation is complicated by factors related both to the mineralization of proxy materials and to the isotope-elevation relationship of meteoric waters from which the proxies form (Poage and Chamberlain, 2001). With regard to the latter factor, studies using global climate models have demonstrated that the isotopic lapse rate can be dependent on a mountain range's elevation due to processes that are not described by Rayleigh distillation Poulsen et al., 2010;Feng et al., 2013;Botsyun et al., 2016). Indeed, Feng et 15 al. (2013) showed that the simulated d 18 O of precipitation (d 18 Op) during uplift of the Eocene North America Cordillera was substantially influenced by changes in vapor mixing, surface recycling, moisture source change, and precipitation type.
Similarly, Botsyun et al. (2016) investigated d 18 Op across the Himalaya-Tibet region in response to surface uplift and showed that direct topographic effects only partially accounted for total d 18 Op changes.
On the Himalayan slope, d 18 O in surface waters and precipitation has been widely observed to decrease with elevation at a 20 rate of ~3 ‰/km (e.g. Rowley et al., 2001). The d 18 O-elevation relationship has been attributed to orographic rainout and modeled as a Rayleigh distillation process (e.g. Rowley and Currie, 2006). On the Tibetan plateau, d 18 O in surface water and precipitation increases linearly with latitude by ~1 ‰/° over nearly uniform elevation (e.g. Bershaw et al., 2012). The source of d 18 O variations on the plateau and whether d 18 O can be used for paleoaltimetry on the high Tibetan Plateau has received little attention. Quade et al., (2011) proposed that paleoelevations could be inferred from proxy d 18 O after removing the 25 meridional d 18 O gradient on the Tibetan Plateau. This south-to-north gradient in d 18 O has been reported to have existed since the early Eocene (Caves Rugenstein and . However, little is known about either the processes that contribute to the meridional d 18 O gradient or how the meridional d 18 O gradient varied when the plateau was lower. The goal of this study is to identify and quantify the processes that control d 18 Op variations across the Himalayas and the Tibetan Plateau and to evaluate the utility of d 18 Op as a paleoaltimeter in these regions. To do this, we use an isotope-enabled 30 global climate model, ECHAM5-wiso, with prescribed elevation scenarios and compare the annual-mean precipitationweighted d 18 Op calculated by the climate model with that expected due to Rayleigh distillation alone. Botsyun et al. (2016) used the LMDZ-iso model to decompose the influence of adiabatic elevation changes on d 18 Op from other influences due to non-adiabatic temperature changes, local changes in relative humidity and post-condensational processes. Building on Botsyun et al. (2016), we take a more process-oriented approach to quantify the isotopic fluxes attributed to specific mechanisms and demonstrate that the contributions from these processes vary spatially and in response to elevation change.
Finally, we discuss the implications of our results for reconstructing paleoaltimetry of the Himalayas and Tibetan Plateau.

Model and experimental design
In this study, we employ ECHAM5-wiso, a water isotope-enabled atmospheric global climate model (AGCM). The model has been widely used both for modern and past climate and isotope simulations. For instance, Feng et al. (2013) and employed ECHAM5-wiso to explore climate and isotopic responses to Cenozoic surface uplift and climate change in western North America. ECHAM has also been shown to simulate many aspects of Asian climate (e.g. 10 Battisti et al., 2014) and isotopic compositions (d 18 Op) simulated by ECHAM5-wiso generally agree well with observed modern stream and precipitation δ 18 O across the Tibetan Plateau as shown by Li et al. (2016) (see their Fig. 11).
We use a model configuration with 19 vertical levels, and a spectral triangular truncation of 106 horizontal waves, approximately equivalent to a 100-km grid spacing. This horizontal resolution, though still relatively coarse, is about twice that of recent simulations used to validate the simulation of water isotopes over the Tibetan plateau  and 15 recent paleoclimate simulations of the Tibetan Plateau (e.g. Roe et al., 2016). The AGCM is coupled to a slab-ocean model, In ECHAM5-wiso, water isotopologues are included as independent tracers in the atmosphere. When water evaporates from 20 the sea, both equilibrium and non-equilibrium distillation processes occur as a function of sea surface temperature, wind speed, relative humidity, isotope composition in seawater and vapor above the ocean surface (Hoffmann et al., 1998).
Convective rains are assumed to have larger raindrops that reach only partial (50%) isotopic equilibration with surrounding vapor, while large-scale precipitation has smaller raindrops and attains almost complete (90%) equilibration with the environment (Hoffmann et al., 1998). For kinetic fractionation of raindrops during partial evaporation in ECHAM5-wiso, the 25 fractionation factor is formulated to depend on the sub-cloud relative humidity of the entire grid box (Hoffmann et al., 1998).
Only large lakes, the size of at least one-half grid cell, are resolved in the model, and fractionation from the land surface is not included since its impact on precipitation d 18 Op is negligible (Haese et al., 2013).
We conducted two sets of sensitivity experiments ( Fig. 1) in addition to a control simulation with modern conditions (CNTL). In the first set of sensitivity experiments, topography is uniformly lowered to 80% (TOPO80), 60% (TOPO60), 30 40% (TOPO40) and 20% (TOPO20) of its modern elevation over a domain that includes the Himalaya and Tibetan plateau.
In the second set of sensitivity experiments, we conducted two experiments with non-uniform elevation modifications over this domain. The first experiment includes a high Himalayan front with a Tibetan Plateau reduced to 20% of its modern elevation (TOPO20a). This experiment is inspired by the widely accepted notion that southeast Asia had an Andean type mountain belt before the collision of the Eurasia and the Indian plate (Royden et al., 2008). TOPO20a also serves as a test of the Himalaya on the regional climate. In the second experiment, the outer edge of the Tibetan Plateau remains, but the inside 5 is lowered to 20% of its modern height (TOPO20b). The second experiment is a sensitivity test to investigate the role of plateau heating on regional climate and isotopic compositions. In both sets of experiments, we tapered the topography along the borders of the domains to avoid any abrupt topography boundaries. Except for topography, all other boundary conditions are kept the same among all experiments. Each experiment was run for 20 years, with the last 15 years used for analysis.
Summertime (June-July-August) climate variables and annual-mean precipitation-weighted δ 18 O are analysed and presented, 10 reflecting that total precipitation over most of the region is dominated by summer precipitation and that carbonates preferentially form in summer when precipitation peaks (Peters et al., 2013).
Mean climate conditions today vary across the Himalaya; the western Himalaya is characterized by peak precipitation in winter and early spring, while the central Himalaya is dominated by the Indian Summer Monsoon (IM) and the eastern Himalaya by the East Asia Summer Monsoon (EASM) (Yao et al., 2013). Because of this heterogeneity, we separate the 15 Himalayas into four distinct regions for analysis purposes, including the western Himalaya; a transitional area between western and central Himalaya; central Himalaya; and eastern Himalaya (transitional areas between central and eastern Himalaya are excluded because climate and isotopic signals are similar to IM and EASM regions). The strength and pattern of precipitation and wind, and isotopic compositions in the two eastern-most transitional areas are similar to those in the IM and EASM in most cases. Thus, in the following we present results only for the western, transitional, IM and EASM regions 20 as shown in Fig. 2.

Rayleigh distillation
We developed an open-system, one-dimensional, altitude-dependent Rayleigh distillation model (RDM) in order to estimate decreases in d 18 Op due to Rayleigh distillation during ascent. The RDM tracks the isotopic composition of an air parcel as it ascends adiabatically from low to high altitude, becomes saturated, and loses condensate through precipitation. In the RDM, 25 an air parcel cools at the dry adiabatic lapse rate before condensation and at the moist adiabatic lapse rate upon saturation (Rowley and Garzione, 2007). The RDM is run using terrain-following coordinates and is initialized with three different moisture sources: (1) fixed air temperature (T=20°C) and relative humidity (RH=80%), (2) local, summertime low-level T and RH from ECHAM5, and (3) fixed T=20°C and summertime RH from ECHAM5. In this way, we are able to quantify the influence due to total moisture source change ((1) minus (2)), and further decompose this influence into the changes in T ((3) 30 minus (2)) or RH ((1) minus (3)).
In order to estimate how much of the mass flux of 18 O in total precipitation is contributed by the Rayleigh distillation process, we assumed that all large-scale precipitation, Pl, forms in response to stable upslope ascent and participates in Rayleigh distillation. We then estimated the isotopic flux of water undergoing Rayleigh distillation as: where RD has units of g m −2 h −1 ; 5:;<= (in g m −3 ) is the density of water; Pl (in m h −1 ) is the large-scale precipitation rate; 5 d 18 ORDM is the isotopic composition simulated by the RDM at the same elevation as the grid points in 2.3 where the mass flux in total precipitation is estimated. Note that this estimation of RD stands for the upper limit of the contribution of RD, since not all large-scale precipitation is triggered by the Rayleigh distillation process.
To compare the relative importance of Rayleigh distillation with other isotopic fractionation processes in ECHAM, we quantify the change in upslope d 18 Op attributable to Rayleigh distillation. We do this by comparing the rate of upslope d 18 Op 10 change (hereafter referred to as the d 18 Op lapse rates) in ECHAM with that estimated using the RDM. We define the d 18 Op lapse rate as the slope of the linear regression equation of precipitation-weighted d 18 Op regressed on elevation and use the coefficient of determination (R 2 ) of this regression to evaluate the robustness of the d 18 Op-elevation relationship. We use a ratio of the d 18 Op lapse rate in ECHAM to that in the RDM (p_percent) to approximate the contribution of Rayleigh distillation to the ECHAM d 18 Op lapse rate. When R 2 and p_percent are both above 0.5 for a particular domain and elevation 15 scenario, we consider the ECHAM d 18 Op lapse rate to agree with the RDM d 18 Op lapse rate. d 18 Op decreases as an air parcel travels inland. This continental effect may contribute to a significant portion of the decrease in d 18 Op over elevated regions in the low elevation scenarios. To evaluate whether a strong d 18 Op-elevation relationship exists in low-elevation scenarios, we examine the ratio of the slope of d 18 Op regressed against latitude on the subcontinent to the south of the Himalayas to that on the Himalayan slope. This method accounts for the fact that d 18 Op changes with both 20 latitude and elevation. A significant d 18 Op-elevation relationship is signified by a large ratio, indicating that Dd 18 Op with elevation is greater than that with latitude.
Note that the centered-finite-difference method is used in discretizing the derivatives in Eq. (4). This method could potentially bring errors in comparison to the spectral method used in the dynamical core of ECHAM5.
Recycling of surface water vapor transports 18 O to the atmosphere from lower boundary. To estimate this contribution to 18 O, 10 the recycled mass flux (in g m −2 h −1 ) is calculated as: where d 18 O is the isotopic composition of the evaporated water and E (m h −1 ) is the surface evaporation rate.
The mass flux of 18 O in total precipitation is estimated following (3) as: where d 18 O is the isotopic composition of precipitation and P (m h −1 ) is the precipitation rate.
Note that this method does not allow us to isolate within-column processes, including vertical mass exchanges through convective updrafts and downdrafts and through phase changes. We encourage future studies to isolate these within-column processes and how they evolve with topography. Plateau in winter and spring. Simulated d 18 Op is 4 ‰ greater than stream sample values along a cross section extending westward from 85° E and centered on 30° N. Li et al. (2016) attributed this mismatch to local factors, systematic model bias, and the influence of freshwater discharge from higher altitudes in the watershed.
We further compare ECHAM5 CNTL d 18 Op with the modern surface water isotope dataset reported in Li and Garzione (2017) (Fig. 5). Co-located ECHAM and water sample d 18 Op agree within 2 ‰ at 49.2% of sites, and within 3 ‰ at 73.3% of sites. Several discrepancies between simulated and observed d 18 Op exist: Firstly, ECHAM5 d 18 Op is lower than sampled 5 d 18 Op over northwestern Tibet (80° E-85° E, 35° N-37° N). This mismatch could be associated with the higher relative humidity in ECHAM (Fig. S1) than that in observations. This high relative humidity results in weaker evaporation both from land surface and below cloud-base, lowering d 18 O in surface waters. Secondly, ECHAM5 d 18 Op is more depleted, by 2-5 ‰, The LMDZ-iso model shows a similar mismatch of 1-4 ‰ (Botsyun et al., 2016). In other regions of Tibet, simulated d 18 Op in LMDZ-iso are very similar to those in ECHAM5, within 2‰ (Yao et 10 al., 2013). There are two potential explanations for the mismatch from modern surface water over east-central Tibet. One possibility is that the sampled d 18 Op does not reflect mean climatic conditions. We note that the water samples from the eastcentral and northern Tibet regions were collected over a short two-year span. This possibility is supported by the large interannual variability in d 18 Op (up to 9 ‰) in both ECHAM  and in precipitation sample spanning 1986 to 1992 from GNIP (AEA/WMO, 2017). The other source for the mismatch over east-central Tibet is from the higher than observed 15 precipitation simulated by ECHAM5, which can result in lower d 18 Op through the amount effect (as shown in 3.6). The simulated summertime precipitation rate over east-central Tibet (89° E-102° E, 32° N-35° N) is 3.9 mm/day, higher than the CMAP value of 2.8 mm/day. This high precipitation is systematic and is also seen in other models, such as the LMDZ-iso (Botsyun et al., 2016) and ECHAM4 (Battisti et al., 2014).

Climate response to Tibetan-Himalayan surface elevation 20
Numerous modeling studies have found that lowering the height of the Tibetan Plateau influences regional temperature, wind, precipitation, and relative humidity (e.g, Kitoh, 2004;Jiang et al., 2008). In this section, we describe regional climate changes on the western Himalayan slope, the Tibetan Plateau, the IM region and the EASM region, that occur as high elevations are lowered.
Under modern conditions, near-surface temperatures in the Tibetan-Himalayan region vary with elevation following a moist 25 adiabatic lapse rate (~5 °C km −1 ) (Fig. S2) and range from >30 °C on the Indian subcontinent at the foot of the Himalaya to <10 °C across the Tibetan Plateau (Fig. S3a). Lowering elevations in our experiments causes near-surface temperatures across the Tibet-Himalayan region to increase ( Fig. S3b-g) at approximately the lapse rate for CNTL. As a result, temperature lapse rates vary little among elevation scenarios, for instance, ranging from 4.9-5.4 °C km −1 in the IM region ( Fig. 1). 30 Wind patterns, precipitation and RH respond dramatically to reductions in elevation. These changes vary across regions. On the western Himalayan slope, wind directions nearly reverse with southerly winds in the high-elevation scenarios switching to northwesterly winds in low-elevation scenarios (Fig. 6). This wind reversal results in transport of arid air from the north, lowering total-column relative humidity by ~40%. With this substantial decrease in RH, summer precipitation decreases from ~3 mm day −1 in CNTL to ~0.1 mm day −1 on the western Himalayan slope in TOPO20. A similar decrease in RH, by ~20% from CNTL to TOPO20, is simulated on the Tibetan Plateau. Accompanying this reduction in RH, precipitation on the Tibetan Plateau decreases from ~4 mm day −1 in CNTL to ~1 mm day −1 in TOPO20. This reduction in Tibetan Plateau 5 precipitation is linked to a weakening of the Asian monsoonal systems and moisture delivery through monsoonal winds.
The responses in the monsoonal regions are somewhat different from those on the western slope and Tibetan Plateau. A reduction in surface elevation (from CNTL to TOPO20) ( Fig. 1) leads to weakening of the IM, as indicated by slowing of summer southwesterly winds over the Bay of Bengal and the Arabian Sea and a decrease in summer precipitation (by more than 20 mm day −1 in ECHAM) along the central Himalayas (Fig. 6). IM weakening is also demonstrated by the WSI1 10 monsoon index (Fig. 7, following Wang and Fan (1999)), which is defined by the vertical wind shear between the lower (850 hPa) and upper (200 hPa) troposphere in the region (5-20° N, 40-80° E) during summer. Monthly WSI1 index values (not shown) indicate that the Indian monsoon persists through all elevation scenarios, although summertime WSI1 (Fig. 7) shows that this Indian Monsoon weakens abruptly once elevations are reduced to between 40-60% of modern values, a threshold reported in previous modelling studies (e.g. Abe et al., 2003). As a result of IM weakening, total-column average RH 15 decreases from >90% in CNTL to 70% in TOPO20 and summer precipitation decreases from >16 mm day −1 in CNTL to ~6 mm day −1 in TOPO20.
ECHAM5 captures a similar threshold behavior in monsoon activity in the EASM region. With lowering of the Himalayan front to 40% of its modern elevations, the broad humid belt that characterizes central China in high elevation scenarios (CNTL, TOPO80, TOPO60, TOPO20a and TOPO20b) shifts southward, resulting in an expanded arid belt in the region (in 20 TOPO40 and TOPO20). This shift in precipitation is associated with a southward retreat of the southwesterly monsoonal winds that penetrate much of eastern Asia (Fig. 6d, e).

Climate response to Himalayan surface-elevation
ECHAM5 experiments TOPO20a and TOPO20b isolate the influence of Himalayan elevations on the regional climate.
These experiments generally indicate that the Himalayas, rather than the Tibetan Plateau, govern regional precipitation and 25 circulation patterns. The IM and EASM are strong in both TOPO20a and TOPO20b. The IM is shown by strong low-level wind over the Arabian sea and heavy precipitation across both the Indian subcontinent and the central Himalayas in TOPO20a and TOPO20b ( Fig. 6f-g). Likewise, the EASM in these simulations is similar to that in the CNTL as indicated by heavy precipitation to the north of Yangtze river and southerly winds penetrating central China ( Fig. 6f-g).
To further elucidate the contribution of the Himalaya to monsoonal dynamics, we calculated the equivalent potential 30 temperature (Fig. S4), which is commonly used to denote the location of monsoonal heating for the Indian Monsoon. In TOPO20a and TOPO20b, the equivalent potential temperature maxima are reduced but in a similar location as in the CNTL case, supporting our conclusion that the Himalayas are the dominant driver of the IM. When the Himalayas are lowered, monsoonal heating decreases and the locus shifts southeastward as cold, dry extratropical air moves southward and mixes with warm, humid subcontinental air, consistent with the results in Boos (2015), though the extent of the shift is smaller in ECHAM5.
In contrast to this well-established mechanism for the Indian Monsoon, the mechanism for the southward shift of the EASM is not well understood. Uplift of the Himalaya-Tibet orogen (Guo et al., 2008;Liu et al., 2017), retreat of the Paratethys 5 (Guo et al., 2008), and forcing by atmophere pCO2 (Licht et al., 2014;Caves Rugenstein and Chamberlain, 2018) have been proposed as possible factors triggering this southward shift, though the timing for this shift is highly debated from Eocene to early Miocene. Our simulation of a strong EASM over central China in both TOPO20a and TOPO20b suggests that uplift of the central and western Himalaya would have been capable of forcing this southward shift.

Moisture source influence on RDM d 18 Op 10
Under lower elevation scenarios, d 18 Op values increase relative to the modern elevation scenario both on the Himalayas and the Tibetan Plateau (Fig. 4). The rate of change of d 18 Op with elevation and latitude decreases substantially in the monsoonal regions as Tibetan-Himalayan elevations are reduced (Fig. 4a, b).
The oxygen isotope compositions of precipitation on mountain slopes are traditionally assumed to systematically decrease in response to adiabatic cooling, condensation, and rainout of ascending air parcels, a process described by Rayleigh distillation 15 and the basis for the application of d 18 Op-paleoaltimetry. In this section, we evaluate the degree to which Rayleigh distillation accounts for up-slope decreases in d 18 Op by comparing RDM to ECHAM d 18 Op in four separate regions (Fig. 2).
Note that the northern Tibet slope is excluded here because most of the northerly air is diverted rather than being forced to ascent (Fig. 6).
For each of the four regions, we initialize the RDM with both fixed moisture sources (with initial T=20°C and RH=80%) as 20 in paleoaltimetry studies and with ECHAM moisture sources (with T and RH that varies by region and case). The d 18 Op predicted by the RDM using each moisture source is shown in Fig. 8-9 (red diamonds vs. black circles). With ECHAMderived T and RH, the RDM d 18 Op values are close to those simulated by ECHAM5 in all elevation scenarios ( Fig. 8- fact that the RH of the initial air parcel is lower in ECHAM than in the prescribed case (Fig. S1). The impact of moisture source differences on d 18 Op is further decomposed in Table 1 to estimate the contributions adiabatic temperature versus RH changes. As seen in Table 3, the effect of adiabatic temperature changes is consistently small (~−1 ‰) across all elevation scenarios, reflecting the fact that temperature lapse rates vary little among elevation scenarios (Fig. S2). In contrast, as the initial RH decreases with lowering of elevation, ECHAM-sourced d 18 Op is lowered by as much as 3.5 ‰ (Table 1, last  30 column). In sum, these results demonstrate that elevation-related changes in moisture source characteristics substantially impact RDM d 18 Op estimates.

Performance of ECHAM-sourced Rayleigh distillation in the Himalayans
Our estimates using an RDM implicitly assume that Rayleigh distillation is the dominant process controlling the d 18 Opelevation relationship. To test this assumption, we compare ECHAM-sourced RDM d 18 Op and ECHAM5 d 18 Op.
Decreases in d 18 Op with elevation in ECHAM and the RDM generally agree under high-Himalaya scenarios and are consistent with modern observations (Table 2). Under modern topographic scenarios for the western Himalaya (Fig. S5) and 5 the monsoonal regions (Fig. 8, 9), R 2 and p_percent values are greater than 0.77 and 0.51, indicating that RDM and ECHAM d 18 Op match well. R 2 and p_percent values are similarly high, above 0.68 and 0.53, for the western Himalaya and the monsoonal regions in high-Himalaya scenarios (TOPO80, TOPO60, TOPO20a, and TOPO20b in the IM region; TOPO80, TOPO20b in the EASM region) again indicating a good match between ECHAM and RDM d 18 Op and suggesting that Rayleigh distillation drives isotopic compositions in these regions. 10 In other regions and under low-elevation scenarios, however, the comparison between RDM and ECHAM d 18 Op is poor (Table 2). For instance, on the western Himalayas, in the TOPO80 case (Fig. S5), the lapse rate of ECHAM d 18 Op is much smaller than the lapse rate predicted by the RDM (Fig. S5), as shown by p_percent values of less than 0.15. Under even lower elevation scenarios (TOPO60, TOPO40, TOPO20), orographic precipitation is not triggered over the Himalayas (Fig.   S4c-e), making the RDM an unsuitable representation of precipitation processes. In the transitional region, the d 18 Op can 15 vary by more than 5‰ at a specific elevation under all elevation scenarios (Fig. S6). This large spread is represented by a low R 2 value of 0.12 for the CNTL case, and even lower values (less than 0.10) for other topographic scenarios. In TOPO40 and TOPO20, ECHAM d 18 Op shows little relationship with elevation (Fig. 4d).
In the monsoonal regions, the relationship between ECHAM5 d 18 Op and elevation is weak in low elevation scenarios (Table   2 and Fig. 4) and compares poorly with the RDM. In the IM region, p_percent values for the TOPO40 and TOPO20 are less 20 than 0.29. In the EASM region, ECHAM5 d 18 Op is higher than RDM d 18 Op (Fig. 8, 9) in the TOPO60 scenario and the agreement is low (with p_percent value of 0.31). Under even lower elevation scenarios (TOPO40 and TOPO20), ECHAM d 18 Op shows no relationship with elevation (Table 2 and Fig. 4b) and also compares poorly with the RDM (with p_percent less than 0.29).
Among the high-Himalaya-low-Tibet cases (TOPO20a and TOPO20b), RDM and ECHAM d 18 Op match well in the IM 25 region. The match is similarly good for TOPO20b in the EASM region. However, in the case with a low eastern flank (TOPO20a), the comparison is poor with an R 2 of 0.006 and a p_percent of −0.05.

Factors influencing d 18 Op-elevation relationship on the Himalayan slope
As shown in 3.5, Rayleigh distillation cannot explain d 18 Op variations with elevation for most regions in the reduced elevation (TOPO40, TOPO20) scenarios and in many regions under higher elevation scenarios (TOPO80, TOPO60). To 30 understand the factors influencing d 18 Op, we quantify the mass fluxes of 18 O (Fig. 10), distinguishing the processes that increase the mass flux of 18 O of the column (mixing and surface recycling) from those that decrease it (Rayleigh distillation and convective rainfall). Note that this method of taking the vertical column as a whole does not isolate the processes occurring within the air column (e.g. sub-cloud re-evaporation, vertical advection and mixing); this limitation does not impact our ability to identify the contributions of Rayleigh distillation and local processes. The results from this method yield very different contributions on the western slope from those in other regions, thus, the western Himalayas are reported 5 separately.
On the western Himalayas, local surface recycling and convective rainfall contributes substantially to the total mass flux of 18 Op in the highest elevation scenarios (Fig. 10a). The contributions from these processes account for the poor match between RDM and ECHAM d 18 Op in this region. In the TOPO20a and TOPO20b scenarios, enhanced transport of enriched vapor from the south (see  (Fig. 10c, d). This change in relative importance of the two sources represents an increase in the importance of local versus remote sources as monsoon strength weakens under low-elevation scenarios.
In high elevation scenarios, Rayleigh distillation acts as the dominant sink of 18 O in the monsoonal regions causing d 18 O to decrease markedly with elevation ( Fig. 10c-d). Under reduced elevation scenarios, the absolute mass flux from Rayleigh distillation decreases, as large-scale precipitation due to stable upslope ascent decreases and convective precipitation 20 increases (Fig. 11). With reduced elevation, the percentage of large-scale precipitation to total precipitation falls from 86% to 18% in the transitional region, from 93% to 18% in the IM region and from 80% to 15% in the EASM region, mirroring the decrease in the mass contribution of Rayleigh distillation, which falls from 85% to 11% in the transitional region, from 93% to 22% in the IM region and from 80% to 18% in the EASM region. As a result of the reduction in (large-scale) precipitation by stable upslope ascent, convective precipitation becomes the largest 18 O sink in TOPO20 and TOPO40 in 25 both the transitional region and the IM region, and in TOPO60, TOPO40, TOPO20 and TOPO20a in the EASM region.
These are also the scenarios that exhibited a poor match between RDM d 18 Op and ECHAM d 18 Op (Sect. 3.5). The increase in convective rainfall in these cases leads to greater kinetic fractionation through sub-cloud evaporation of falling rain, which is only partially equilibrated with the surrounding vapor (see 2.1), and an enrichment in the isotopic composition of rain. The RDM does not capture this enrichment because it does not include sub-cloud evaporation. A similar reduction in large-scale 30 precipitation with reduction in elevation is also captured in the Andes region (Insel et al., 2009) and is associated with a decrease in the rate of change of d 18 Op with elevation.
Note that, although Rayleigh distillation is the primary sink under high-elevation scenarios, RDM d 18 Op does not match ECHAM d 18 Op well in the transitional region because of the large spread in ECHAM d 18 Op (Fig. S6). This spread is due to the bifurcated sources from northwest India (relatively enriched) and the Bay of Bengal (relatively depleted), since air parcels follow separate trajectories before mixing at the peak (Fig. S7).
In sum, the mismatch between RDM and ECHAM d 18 Op is caused by a weakening of Rayleigh distillation under low-5 elevation scenarios, triggered by a reduction in large-scale precipitation.

d 18 Op-latitude relationship on the Tibetan Plateau
In contrast to d 18 Op on the Himalayas, d 18 Op on the Tibetan Plateau increases linearly with latitude. It has been proposed that  Fig.12a, b) and variations are larger in high-elevation cases (e.g. from 1.90 ‰/° to 0.20 ‰/° in CNTL). Secondly, a linear fit is generally good with a high coefficient of determination (R 2 > 0.8) east of 85°W for all cases but TOPO20a and TOPO20b (Fig. S8). The high goodness of fit indicates a robust d 18 Op-latitude relationship for the cases with uniform reductions in topography. The poor fit in TOPO20a and TOPO20b is due to the larger variations in elevation across the Tibetan Plateau. In these two cases, dry conditions (Fig. S1) on the steep leeward side of the Himalaya favor strong below 20 cloud-base re-evaporation, increasing local latitudinal d 18 Op gradients. Lastly, in cases with uniform lowering of topography, the median meridional gradient (Fig. 12b) decreases almost linearly with reductions from 100% to 60% of modern topography, changes more abruptly between 60% and 40%, and varies little once the topography is lowered below the monsoon threshold in TOPO40 and TOPO20.
The modern latitudinal d 18 Op gradient on the Tibetan Plateau has been attributed to both surface recycling and sub-cloud re-25 evaporation. Surface recycling increases d 18 Op by adding more enriched vapor from the land surface, while sub-cloud evaporation enriches 18 O in precipitation by partial evaporation of falling raindrops in an unsaturated air column. To understand the cause of the decline in the meridional d 18 Op gradient, we quantified both the sources of 18 O (surface recycling and vapor mixing) and the isotopic contribution due to sub-cloud re-evaporation. Sub-cloud evaporation is quantified separately because this effect is contained in but cannot be isolated from the sources of 18 O. 30 Surface recycling is the primary source of 18 O on the Tibetan Plateau with largely consistent contributions under all elevation scenarios (Fig. 13). The mass flux of 18 O due to vapor mixing is small and becomes a sink in low-elevation scenarios, since there are increasingly more depleted sources from the south in these scenarios (e.g., 17 out of 42 in CNTL vs. 30 out of 42 in TOPO20 for one location (33° N, 90° E)). To further quantify the contribution of surface recycling to total precipitation, we where d 18 Op is the isotopic composition in precipitation; d 18 Os is the condensate of recycled vapor; E (m h −1 ) is the surface evaporation rate and P (m h −1 ) is the total precipitation rate. Results from Eq. (7) show that d 18 O of surface recycling contributes less than 0.1 ‰ in CNTL and slightly more (−2.67 ‰) in TOPO20 (Table 3). The overall small contribution from surface recycling in ECHAM5 is due to the fact that d 18 O values in soil are similar to those in precipitation (Fig. S9).
This small contribution grows in low-elevation scenarios due to the increased fraction of evaporated vapor to total 10 precipitation (Fig. S10). Nonetheless, surface recycling plays a secondary role in decreasing the meridional d 18 Op gradients.
Sub-cloud re-evaporation occurs within an unsaturated air column as falling raindrops undergo kinetic fractionation and become isotopically enriched (Stewart, 1975). Enrichment is reduced in heavier rain where the relative humidity is high, resulting in the observed anti-correlation between d 18 Op and precipitation rate referred to as the amount effect. To quantify this enrichment due to sub-cloud re-evaporation, we show d 18 Op against daily precipitation (Fig. 14). The slope of this d 18 Op-15 precipitation relationship denotes the strength of the amount effect. In low-elevation scenarios, the d 18 Op-precipitation slope is shallower (Fig. 14), indicating a weaker amount effect and stronger sub-cloud re-evaporation even at high precipitation rates.
This shallow slope of sub-cloud evaporation results in the shallower meridional d 18 Op gradient on the Tibetan Plateau in lowelevation scenarios. In each individual elevation scenario, the enrichment of d 18 Op due to sub-cloud evaporation is stronger at 20 higher latitudes since there are more events of lower precipitation rates at higher latitudes than at lower latitudes (Fig. 15).
As a result of this different distribution of precipitation, rainfall at higher latitudes is more enriched than that at lower latitudes. To further compare this contribution of sub-cloud re-evaporation with surface recycling, the excess enrichment of d 18 Op due to sub-cloud re-evaporation is quantified by the d 18 Op difference at high and low precipitation rates, weighted by precipitation rates (Table 3). The excess enrichment due to sub-cloud re-evaporation is much larger than that due to surface 25 recycling and decreases with reduced elevation.
To explain this stronger sub-cloud re-evaporation in low-elevation scenarios, we refer to the kinetic fractionation process in ECHAM5 during partial evaporation of raindrops. As shown in Hoffmann et al. (1998), kinetic fractionation in ECHAM5 is formulated as: where a is the fractionation factor; × b d, represents the ratio of diffusivities between 16 O and 18 O and has the constant value of 0.9727; RH is the effective relative humidity of the grid box. To estimate how this kinetic fractionation changes between different elevation scenarios, we approximated the effective relative humidity to be the total-column-averaged RH.
As seen in Fig. S11, the total-column-averaged RH at any given precipitation rate is higher in high-elevation scenarios than that in low-elevation scenarios. Specifically, when the precipitation rate is very high at 40mm day −1 , the RH is at ~100% in CNTL, suggesting very little kinetic fractionation and weak sub-cloud re-evaporation. In comparison, in TOPO20 the RH is much lower at ~85%, indicating the presence of sub-cloud re-evaporation even at high precipitation rates. 5 In sum, meridional d 18 Op gradients on the Tibetan Plateau decrease with lower elevation, and this reduction is due to stronger sub-cloud evaporation in low-elevation scenarios.

Processes impacting paleoaltimetry
The use of d 18 Op as a paleoaltimeter is based on the principle that atmospheric vapor becomes saturated and condenses under 10 forced ascent in orographic regions leading to preferential rainout of 18 O. Observations of meteoric d 18 O in orographic regions support the use of stable isotope paleoaltimetry and show a robust and significant decreasing relationship in d 18 Op with surface elevation (e.g. Poage and Chamberlain, 2001;Fiorella et al., 2015;Li and Garzione, 2017) that can be modeled by Rayleigh distillation (Rowley and Garzione, 2007). However, there is no a priori reason that modern d 18 O-elevation relationships should hold in the past when surface elevations and associated atmospheric conditions were different (e.g. 15 Ehlers and Poulsen, 2009;Poulsen et al., 2010;Feng et al., 2013;Botsyun et al., 2016).
Our ECHAM5 results confirm that the processes that govern d 18 Op vary spatially and change in response to changes in surface elevation. Across the Himalayan front, d 18 Op generally decreases with elevation, particularly in scenarios with high Himalayan elevations. However, Rayleigh distillation often does a poor job of explaining the d 18 Op-elevation relationship in the Himalaya-Tibetan region. There are two primary reasons. Firstly, local processes including convection and surface 20 recycling can dominate the land-surface exchange of vapor (Fig. 10). These local processes are especially strong under lowelevation scenarios in monsoonal regions, where forced ascent is weak, and in the western Himalayas. Secondly, in regions with multiple moisture sources, mixing of air masses with different d 18 O can cause the d 18 Op-elevation signal to deviate from that expected from Rayleigh distillation. The transitional region on the Himalaya, which receives moisture from both the Arabian Sea and the Bay of Bengal, is such a case (Fig. S7). On the Tibetan Plateau, the meridional d 18 Op gradient decreases 25 in response to reduced elevation. This meridional d 18 Op gradient is primarily controlled by sub-cloud re-evaporation under modern elevations. Sub-cloud re-evaporation increases due to wetter conditions in low-elevation scenarios, reducing the meridional d 18 Op gradient. Our conclusions are similar to those of Feng et al. (2013) for the North American Cordillera, in which it was shown that the isotopic fractionation of precipitation was not primarily due to Rayleigh distillation.
Our results also highlight the large influence that the choice of moisture source characteristics has on RDM d 18 Op and are 30 consistent with those of Botsyun et al. (2016). When we use fixed T and RH, ECHAM5 and RDM d 18 Op agreement is poor with considerable enrichment in RDM. When we use ECHAM5 T and RH, the agreement is considerably improved. This change in moisture sources represents elevation-induced climate change unrelated to rainout and not captured by Rayleigh distillation. The overall impact from temperature and RH (red diamonds in Fig. 8-9) is a slight underestimation in highelevation scenarios and severe overestimation in low-elevation scenarios (e.g. by ~100% in TOPO20 in the IM region).

Implications for d 18 Op paleoaltimetry 5
As discussed above, d 18 Op paleoaltimetry is only appropriate for monsoonal regions in high elevation scenarios. Nonetheless, proxy d 18 Op from sites across the Himalaya-Tibet region have been used to infer paleoaltimetry throughout the Cenozoic (Fig. 16, Table S1). We classify these sites into five types indicating their utility for d 18 Op paleoaltimetry.
As shown in Fig. 16 (black dots), the modern d 18 Op-elevation relationship is well captured by ECHAM5 and well represented by the Rayleigh distillation process under high-elevation scenarios for 7 of 50 sites located on the high 10 Himalayas. However, for these seven sites, the d 18 Op-elevation relationship breaks down as height of the Himalaya-Tibet orogen is reduced. For 13 of 50 sites in central Tibet (green dots), there is no direct d 18 Op-elevation relationship, since elevation is largely uniform and d 18 Op increases linearly with latitude due to sub-cloud re-evaporation and recycling (as shown by 3.7). Although it has been proposed that d 18 Op-paleoaltimetry is still possible after removing this meridional d 18 Op gradient from d 18 Op (Bershaw et al., 2012), the meridional gradient varies by as much as 70% with elevation and is largely 15 unknown in the past (as shown by 3.7); thus, d 18 Op-paleoaltimetry is not appropriate for these sites. For 7 of 50 sites (red dots), elevations are too low so that d 18 Op has no relationship with elevation and changes very little throughout different elevation scenarios (Fig. 4). As a result, d 18 Op is not indicative of elevation. For sites in the western and transitional regions (6 of 50), d 18 Op either relates very little to elevation or exhibits a large range at the same elevation (Fig. S6), and is thus not appropriate for d 18 Op-paleoaltimetry. Lastly, at sites in northern Tibet (yellow dots, 15 of 50), moisture is mostly diverted by 20 the Tibetan Plateau, Rayleigh distillation is not triggered, and d 18 Op does not vary with elevation (Fig. 4).
Although d 18 Op-elevation relationship at sites in monsoonal regions (Fig. 16, black dots) are well explained by Rayleigh distillation in our experiments, d 18 Op-paleoaltimetry should still be applied with care for these sites even under high elevation conditions due to the influence of climate change. Proxy d 18 Op values as low as −18-−16 ‰ are reported for the early Eocene  and around −14 ‰ for the late Eocene (Rowley and Currie, 2006). These values almost 25 certainly reflect high elevations. Nonetheless, elevation-independent factors, including atmospheric pCO2 (Poulsen and Jeffrey, 2011) and paleogeography (Roe et al., 2016), add substantial uncertainty to the quantification of past surface elevations.

Caveats
Like all models, ECHAM has limitations. Most pertinent to this study, ECHAM simulates higher precipitation along the steepest slope of large mountain than indicated by satellite observations (Roe et al., 2016 and reference there in), which is a common problem in most GCMs. The model tends to overestimate the modern precipitation amount and RH on the western Tibetan plateau. This overestimation also exist in earlier version of ECHAM (Roe et al., 2016) and in LMDZ-iso (Zhang and 5 Li, 2016). The overestimation of total-column-averaged RH might weaken sub-cloud re-evaporation process on the western plateau, and this weaker re-evaporation could potentially lower the meridional d 18 Op gradient in this region. Another limitation of the model is that sub-grid scale lakes are not included. Although it has been proposed that the lakes provide Another limitation of our modelling strategy is our use of a slab ocean model, which does not account for ocean circulation changes that would result from the changes in topography that we prescribe. To the best of our knowledge, no study has specifically investigated the response to a reduction in the elevations of the Himalayas and the Tibetan Plateau. It is not hard 15 to imagine regional sea surface changes that might influence inland precipitation. For example, we speculate that under lower elevations and weaker monsoon winds, ocean upwelling along the western coast of Bay of Bengal and the Arabian Sea would be reduced, leading to higher sea-surface temperatures (SSTs). In a study of the East Asian response to historical SST warming, higher SSTs led to greater precipitation over the Indian Ocean and Pacific Ocean due to enhanced local convection and less precipitation along the Himalayan front, and further weakening of the Asian monsoon (Li et al., 2010) In this study, we use a series of idealized simulations to investigate the response of water isotopes to mountain uplift and to understand the mechanisms that control d 18 Op variations on the Himalayas and the Tibetan Plateau. We do not compare and evaluate our model results directly against proxies, since the simulations are not meant to represent specific time slices in the 30 geologic past, and do not include changes in time-specific boundary conditions such as greenhouse gas composition, vegetation, orbital parameters (Tabor et al., 2018), glacial boundary conditions, tectonic displacement of the Indian subcontinent, or inclusion of the Paratethys. Changes in these boundary conditions, consistent with those that occurred through the Cenozoic, are likely to have a substantially influence on simulated d 18 Op. For instance, high early Cenozoic atmospheric CO2 may have increased d 18 Op over high elevation regions (Jeffery et al., 2012;Poulsen and Jeffery, 2011) and on the Tibetan Plateau by as much as 8‰ (Poulsen and Jeffery, 2011). In addition, Earth's orbital variations have been shown to contribute as much as 7 ‰ to oxygen isotope changes on the Tibetan Plateau (Battisti et al., 2014), which is comparable to the isotope difference from CNTL to TOPO20 (about -11 ‰) in Fig. 4. These large variations associated 5 with orbital fluctuations are not observed in long records spanning 10 6 yrs (Deng and Ding, 2015;Kent-Corson et al., 2009), presumably because the terrestrial proxy archives of d 18 O integrate over orbital time periods.

Conclusion
The isotopic composition of ancient meteoric waters archived in terrestrial proxies is often used as a paleoaltimeter under the assumption that rainout during stable air parcel ascent over topography leads to a systemic isotopic depletion through 10 Rayleigh distillation. We use an isotope-enabled GCM, ECHAM5-wiso, to evaluate the extent to which oxygen isotopes can be used as a paleoaltimeter for the Himalaya-Tibet region and to explore the processes that control the d 18 O-elevation relationship. Overall, our study highlights the myriad processes that influence d 18 Op in the Himalaya-Tibet region now and during its uplift.
We find that Rayleigh distillation describes most of the d 18 Op variation with elevation in the monsoonal regions under high-15 topography scenarios. In contrast, Rayleigh distillation does a poor job of describing d 18 Op variation with elevation under high-topography scenarios in the western Himalayas due to the dominance of local convection and surface recycling in that region. When the Himalaya-Tibet elevations are reduced to below one-half of their modern heights, d 18 Op exhibits no relationship with elevation. At these reduced elevations, d 18 O fractionation occurs primarily through local convection and surface recycling. On the Tibetan Plateau under modern elevations conditions, d 18 Op linearly increases with latitude 20 primarily due to sub-cloud re-evaporation. The d 18 Op gradient decreases as the plateau is lowered due primarily to stronger sub-cloud re-evaporation under drier conditions and secondarily to increased moisture sources from surface recycling.
Because of these elevation-independent processes, we conclude that only 7 out of the 50 paleoaltimetry sites are appropriate for d 18 O-paleoaltimetry. Taken together, these results indicate that stable isotope paleoaltimetry in the Himalaya-Tibet region, as in other orogenic regions, is at best a blunt instrument for inferring past surface elevations.               Percentage of occurence CNTL low lattitue CNTL high lattitue TOPO20 low lattitue TOPO20 high lattitue