Research Article

Revised paleoaltimetry data show low Tibetan Plateau elevation during the Eocene

See allHide authors and affiliations

Science  01 Mar 2019:
Vol. 363, Issue 6430, eaaq1436
DOI: 10.1126/science.aaq1436

Ancient height of the Tibetan Plateau

The elevation of the Tibetan Plateau has a major impact on climate, affecting the monsoons and regional weather patterns. Although some isotope proxies have suggested a roughly equivalent height for the plateau as far back as the Eocene (∼40 million years ago), other lines of evidence suggest a lower elevation in the distant past. Botsyun et al. used a model to show that several previously overlooked factors contribute to the isotopic record from the Eocene (see the Perspective by van Hinsbergen and Boschman). The results harmonize the isotopic record with other proxies and argue for a Tibetan Plateau that was about 1000 meters lower than it is today.

Science, this issue p. eaaq1436; see also p. 928

Structured Abstract

INTRODUCTION

The uplift history of the Tibetan Plateau (TP) is critical for understanding the evolution of the Asian monsoons and the geodynamic forces involved in collisional orogens. The early topographic history of the TP is uncertain, and the timing of the initiation of uplift remains controversial. The majority of studies that find evidence for an old plateau (as early as the Eocene, ~40 million years ago) rely on stable isotope paleoaltimetry. This method is based on observations and models that show depletion in heavy oxygen isotopes in rainfall during orographic ascent. Ancient carbonates of the TP record past rainfall isotopes; when available, their isotopic signature can be compared with isotopic lapse rates in order to estimate at what elevation they grew in the past, but the use of this method in deep time has many uncertainties.

RATIONALE

Applying stable-isotope paleoaltimetry in greenhouse climates makes the implicit assumption that the factors that control atmospheric distillation and rainfall oxygen isotopic composition (δ18Ow) have remained constant over millions of years. However, the impact of past climate change on δ18Ow values is unclear. In particular, for the Eocene, higher atmospheric CO2 concentration (Pco2) and markedly different Asian paleogeography, including a wide and shallow Paratethys Sea in central China and a latitudinal shift of the southern Tibet margin ~10° to the south, have been hypothesized to modify the Asian climate and regional δ18Ow values. In addition, the carbonate formation temperature is often unknown, increasing the uncertainty in the reconstructed δ18Ow.

RESULTS

We ran climate simulations with Eocene boundary conditions and varying TP elevation. We accounted for changing Pco2, land surface albedo, orbital variation, and sea-surface temperatures, which potentially cause shifts in δ18Ow, by changing the relative contribution of different air masses and the local hydrological cycle—the evaporation-to-precipitation ratio and the fractioning between convective and large-scale rainfalls. In our simulations, the south-shifted location of the entire Indian foreland induces strong convection over the southern flank of the TP and a radically different pattern of water recycling compared with that of present day. Our simulations reproduce monsoonlike seasonal precipitation over the Indian foreland, with summer convective rainfall reaching the TP. An intense anticyclonic circulation during summer months induces widespread aridity on the northern part of the Plateau. This peculiar atmospheric circulation, together with intensified water recycling and multiple moisture sources, results in a reversed isotopic lapse rate across the southern flank of the TP, with the most negative δ18Ow over northern India and increased δ18Ow northward. On the basis of Eocene experiments with varied boundary conditions, we argue that this is a robust feature of Eocene climate over the region. Furthermore, this pattern is the opposite of the present-day δ18Ow over the Himalayas, which decreases with elevation, driven by orographic rainout following a Rayleigh distillation process. Last, using our simulated temperatures and δ18Ow, we derived virtual carbonates δ18O for different elevation scenarios and compared them with the geological record. Statistical analysis shows that a low TP topography during the Eocene is the scenario that provides the best match between model and data.

CONCLUSION

Our simulations indicate that standard stable isotope paleoaltimetry methods are not applicable in Eocene Asia because of a combination of increased convective precipitation, mixture of air masses of different origin, and widespread aridity. In the Eocene, the presence of a reversed isotopic lapse rate precludes use of any of the previously developed δ18Ow-elevation relationships for estimating the Eocene elevation of the TP. Rather, a model-data comparison on the carbonates δ18O suggests that the TP reached only low to moderate (<3000 m) elevations during the Eocene, reconciling oxygen isotope data with other proxies of elevation and with geodynamic models that propose a recent (Neogene) uplift. More generally, we suggest that using climate models in conjunction with stable isotope data from the geological archives provides a powerful tool to incorporate climatic changes into the analysis of paleoelevations.

Simulated stable oxygen isotope composition of summer precipitation for the Eocene (42 million years ago) Asia.

The combination of changes in water recycling and the source of air masses induces the reversal of isotopic lapse rate. Summer winds trajectories highlight intense westerlies and easterlies at each border of the TP. Circles indicate localities of paleoaltimetry sites used in this study.

Abstract

Paleotopographic reconstructions of the Tibetan Plateau based on stable isotope paleoaltimetry methods conclude that most of the Plateau’s current elevation was already reached by the Eocene, ~40 million years ago. However, changes in atmospheric and hydrological dynamics affect oxygen stable isotopes in precipitation and may thus bias such reconstructions. We used an isotope-equipped general circulation model to assess the influence of changing Eocene paleogeography and climate on paleoelevation estimates. Our simulations indicate that stable isotope paleoaltimetry methods are not applicable in Eocene Asia because of a combination of increased convective precipitation, mixture of air masses, and widespread aridity. Rather, a model-data comparison suggests that the Tibetan Plateau only reached low to moderate (less than 3000 meters) elevations during the Eocene, reconciling oxygen isotope data with other proxies.

Because the Tibetan Plateau (TP) is the largest orogen on Earth, its uplift history is critical for understanding the evolution of the Asian monsoon systems and the geodynamic forces involved in continental collisions (1, 2). Previous reconstructions have proposed a wide range of paleogeographies for the late Eocene of Tibet, which all accommodate India-Asia plate motion through a combination of Indian plate subduction and Asian shortening and extrusion (3), but the altitude of the TP through time remains disputed (4). Numerous uplift scenarios for the TP have been suggested from a precollision, highly elevated Plateau [for example, ~100 million years (Ma) ago] (5) through stepwise uplift since initial collision (~55 Ma ago) (6) to more recent rapid uplift, particularly of the northern and eastern TP (<8 Ma ago) (7). Although any of these uplift scenarios are possible, they imply different geodynamic mechanisms and climate responses. The elevation of the TP since the Eocene as reconstructed in stable isotope paleoaltimetry studies is close to modern values (>4000 m) (810). However, other lines of evidence from paleobotany (11, 12) and paleontology (13, 14) suggest much lower elevations of the TP in the past. Improved constraints on the early topographic history of the TP are critical to better understand the geodynamic evolution of India-Eurasian collision and its consequences on the climate system—namely, the Asian monsoon onset.

Stable isotope paleoaltimetry is based on the observed and theoretically predicted decrease of precipitation water δ18O [δ18Ow, expressed in per mil (‰)] with increasing elevation across mountain ranges (15). In the Himalayas and southern Tibet, the present-day isotopic lapse rate is ~–2.8 ‰/km (16), which is in agreement with distillation models that predict decreasing δ18Ow in response to decreasing atmospheric humidity and temperature during orographic ascent (15). This altitude effect is archived in the δ18O of continental carbonates (δ18Oc), offset from δ18Ow by a fractionation factor that depends on the carbonate formation temperature (17). Stable isotope paleoaltimetry studies commonly provide δ18Oc at sites of unknown elevation; make assumptions about the temperature at the site, the near-sea level δ18Ow, and the atmospheric relative humidity; and from this, reconstruct the paleoelevation by using modern empirical or modeled lapse rates (15, 16).

However, stable isotope paleoaltimetry relies on the key assumption that the main atmospheric factors that control atmospheric distillation and air mass δ18O values have remained constant over millions of years. The impact of past climate change—including changes in lapse rates (18), sources of precipitation (9), and in monsoonal regime (19)—on δ18Ow is uncertain but has been hypothesized to affect δ18Ow over the TP (20). In addition, carbonate growth temperature is often unknown, increasing the uncertainty in the reconstructed δ18Ow (17). All these factors may importantly bias deep-time paleoelevation estimates, particularly during the Eocene Greenhouse world. Higher atmospheric CO2 concentration(Pco2) at that time (21) and markedly different Asian paleogeography, with a wide and shallow Paratethys Sea in central China (22) and a latitudinal shift of the southern Tibetan margin ~10° to the south (23), might have substantially modified the Asian hydrological cycle (24, 25).

We tested the impact of late Eocene paleogeography and topography variations on water vapor and precipitation δ18O using an isotope-enabled version of the LMDZ atmospheric circulation model and compared simulated δ18Oc values with previously published data. LMDZ-iso properly simulates atmospheric dynamics and reproduces rainfall and δ18Ow patterns consistent with present-day data over Asia, including the substantial depletion that occurs across the Himalayan flanks and the stark south-to-north positive gradient in δ18Ow across the Plateau (Fig. 1, A and B). In order to fully account for uncertainties regarding our knowledge of Asian paleogeography during the Eocene, we varied the elevation of the TP, the latitude of the TP, and the land-sea mask (Fig. S1). Additionally, we accounted for changing Pco2, land surface albedo, orbital variation, and sea-surface temperatures. These variables potentially cause shifts in δ18Ow by changing the relative contribution of different air masses and the local hydrological cycle (2628).

Fig. 1 Present-day and Eocene simulations.

(A) Simulated MJJAS precipitation-weighted present-day δ18Ow. Triangles show δ18Ow from GNIP stations (70), circles represent δ18O in streams, rivers, and springs from [(37) and references therein] compilation. (B) The δ18O south-north profile, averaged between 70°E and 90°E. Gray lines indicate minimum and maximum model δ18Ow 70°E and 90°E. (C and D) Simulated Eocene precipitation and isotopic patterns. (C) LMDZ-simulated MJJAS precipitation amount (millimeters per day), and (D) MJJAS precipitation-weighted δ18Ow (‰) and wind velocities (meters per second) for the Eocene (EOC-L) experiment.

Eocene Asian climate and water δ18O

Simulated Eocene large-scale circulation is characterized by monsoonlike seasonal precipitation over the Indian foreland (fig. S2). By contrast, over the TP and over central Asia, westerlies more intense than those of today interacted with large-scale atmospheric subsidence (figs. S2 and S3), causing widespread aridity (Fig. 1C and fig. S2). The onset of an anticyclonic circulation during summer months over the eastern part of the TP (Movie 1) induces intense westerlies and easterlies at each border of the TP. This feature is in stark contrast with the modern simulations, in which the southern TP receives a greater proportion of precipitation from air masses originating in the Bay of Bengal and propagating northward (Movie 1). These processes occur across all of our different paleogeography and paleotopography scenarios. High rainfall over the Indian foreland (Fig. 1, C and D) decreases δ18Ow from the southern tip of India (~–4‰) to a minimum at 15°N over the Indian foreland (from –8 to –11‰, depending upon TP topography) (Fig. 2 and fig. S2). The multiple sources feeding precipitation in the southern part of India explain the absence of a correlation between the rainfall and the value of δ18Ow in the region (Movie 1). Conversely, more negative δ18Ow values over the northern part of the Indian continent arise from distillation processes as specific humidity decreases as air masses progress northward across India owing to high precipitation values. In addition, extremely arid environments (<0.5 mm/day) over the northern TP and central Asia are associated with higher δ18Ow values (from –3 to –1‰) (Fig. 1, C and D). The most surprising and counterintuitive result in our Eocene simulations is an isotopic gradient across the southern flank of the TP, with δ18Ow values ranging from –8‰ at the foot of the TP to –5‰ at the top of the southern flank of the TP. Such a gradient is at odds with modern observed and simulated values that depict a decrease in δ18Ow from –8 to –18‰ as altitude increases on the southern flank of the TP (Fig. 2 and fig. S4). The combination of changes in water recycling and the source of air masses explains the reversal of the isotopic lapse rate during orographic ascent in the Eocene.

Movie 1. Model-simulated δ18Ow and precipitation amount for present-day (CTR) and Eocene (EOC) cases, monthly time slices.

Vector arrows show wind speed and directions.

Fig. 2 Precipitation-weighted MJJAS δ18Ow in simulations.

(A) EOC-L. (B) EOC-XL. (C) EOC-S. (D) EOC-Him. (E) EOC-M. (F) EOC-sea. (G) EOC-North. (H) EOC-South. Vectors show MJJAS vertically integrated moisture transport (kilograms per meter per second).

Back-trajectory analysis

Our back-trajectory analysis confirms that in the Eocene, sites across the TP are influenced by multiple air masses that originate from both the Indian Ocean and the Paratethys Sea (Fig. 3 and Movie 1). In detail, stronger westerlies during the Eocene are divided between a southern branch following the Himalayan flanks and a northern branch that deviates to the south and joins the easterlies around 95°E (Movie 1). The presence of multiple air masses over the Eocene southern TP invalidates two common assumptions (29, 30) in stable isotope paleoaltimetry: that (i) southern TP sites receive moisture from a single air mass and (ii) air masses progress over the southern margin and onto the TP, permitting the resulting δ18Ow to be treated as a Rayleigh distillation process. Instead, δ18Ow over the southern TP appears to be a mixture of multiple air masses having followed along-strike trajectories (eastward and westward winds). Therefore, the assumption of Rayleigh distillation along a single air mass trajectory cannot be used to explain δ18Ow patterns over elevated parts of the Eocene topography. Our modeling results are also supported by present-day satellite observations showing that mixing of various air masses tends to induce higher δ18Ow than a purely Rayleigh regime (31).

Fig. 3 Back trajectory analysis for present-day (CTR) and Eocene winds.

(A to F) The probability distribution of air masses along back trajectories for each of the six points for [(A) to (C)] CTR and [(D) to (F)] EOC-L. Back trajectories are calculated daily from May to September at 900 hPa level, backward for 6 days each, starting from the black circle.

Evaporation to precipitation ratio and precipitation regimes

The south-shifted location of the Eocene TP induces strong convection over its southern flank and a radically different pattern of water recycling compared with the modern pattern (Fig. 4). These features persist even in a scenario in which the TP has been latitudinally shifted (fig. S5). Our present-day simulation reproduces the large decrease in convective precipitation that exists from the Indian foreland up to the southern flank of the TP [evaporation/precipitation ratio (E/P) diminishes along the same transect] (Fig. 4A). However, during the Eocene, convective precipitation (CP) prevails over large-scale precipitation (LSP) (CP/LSP > 1), with no clear decrease in E/P over this region (Fig. 4B and fig. S5). The negligible decrease in δ18Ow in the Eocene between the Indian foreland and the southern flank of the TP together with high convective precipitation rate are therefore consistent with observational and modeling studies (3234) that show that stratiform precipitation regimes are associated with more depleted isotopic values than convective regimes. An additional simulation using a TP of the equivalent elevation as today for the Eocene does not change zonal profiles of E/P and CP/LSP, suggesting that our hypothesized mechanism to explain reverse isotopic gradients in the Eocene over the southern flank of the TP is robust.

Fig. 4 Hypothesized mechanisms that drive differences between Eocene and present-day δ18Ow.

(A and B) The ratio between convective precipitation and large-scale precipitation (CP/LSP, blue lines) and the ratio between evaporation and precipitation (E/P, red lines) along south-north TP cross sections for (A) CTR and (B) Eocene experiments. In (B), dashed lines indicate CP/LSP and E/P for the EOC-XL experiment, and the solid line indicates the same for the EOC-L experiment. Stratiform precipitation is mimicked in the model with large-scale precipitation. Gray shading indicates the topography for the corresponding simulation. Profiles are averaged between 85°E and 90°E for the CTR case and between 73°E and 78°E for the Eocene cases.

Reversed Eocene isotopic lapse rate

The increased convective precipitation rates and multiple air masses over the southern TP result in a latitudinal dipole in δ18Ow, with low δ18Ow over the India foreland and higher δ18Ow over the TP. This dipole creates an apparent positive isotopic lapse rate across the southern edge of the TP, with increased δ18Ow as elevation increases south to north. This positive isotopic lapse rate is a feature of all Eocene simulations, regardless of the prescribed topography, of the latitudinal position of the TP, and of the position of the coastline (Fig. 2). Furthermore, this pattern is opposite to the expected δ18Ow-elevation relationship derived from Rayleigh distillation and driven by orographic rainout, which predicts that δ18Ow should decrease with elevation (15). Flat isotopic lapse rates are found today in areas of high convective precipitation, such as Ethiopia and Kenya (35) because of convective instability over the region that produces enriched δ18O values at altitude. Similar negligible isotopic lapse rates have also been observed in the Andes (36). For the Eocene, this unusual isotopic lapse rate precludes use of any of the previously developed δ18Ow-elevation relationships for the Eocene TP.

Discussion and comparison with data

For the purpose of a model-data comparison, we derived virtual δ18Oc values using simulated δ18Ow and surface temperatures and compared these derived values with measured δ18Oc from previous studies (Table 1 and Figs. 5, A and B, and 6). The δ18Oc data show a reasonable agreement with the simulated south-north isotopic gradient in δ18Oc for the Eocene experiments (Figs. 5, A and B, and 6). However, simulations with relatively low Asian topography (EOC-S and EOC-M) are a better fit, according to a ranking of the model simulations based on the sum of squared residuals (Fig. 5C). The large uncertainties in the δ18Oc data preclude any precise paleoaltimetry estimates but do allow us to discriminate between different paleoelevation scenarios. In simulations with higher TP topography, δ18Oc data are systematically more negative than simulated δ18Oc values (Figs. 5A and 6). A decomposition of the effects on δ18Oc values shows that this discrepancy is mainly due to the decrease in temperature over the TP as elevation increases because of adiabatic cooling (figs. S7 to S9). Lower carbonate formation temperatures result in less negative δ18Oc values over the TP. This effect is enhanced by more intense rainfall over the Indian foreland associated with more depleted δ18Ow values (figs. S2 and S3) as TP elevation increases. These two mechanisms reinforce the positive isotopic lapse rate at the southern edge of the TP, resulting in increased δ18Oc with increasing elevation.

Table 1 Compilation of measured Eocene δ18Oc from published materials.

Corresponding sample sites, geographic locations are provided in fig. S6 (8, 10, 13, 24, 5766). δ18Oc are reported relative to V-PDB. We report the minimum measured values for lacustrine/biotic data and minimum and averaged values for paleosols. Numbers in parentheses are average values for paleosols. Dashes indicate absence of corresponding data in the original study. For δ18Oc, dashes are shown where 1 SD calculation is not applicable. PETM, Paleocene–Eocene Thermal Maximum; Min, minimum.

View this table:
Fig. 5 Model-data comparison for the Eocene.

(A and B) Maps of δ18Oc calculated from model-simulated MJJAS precipitation-weighted δ18Ow and model MJJAS mean temperatures for (A) EOC-XL and (B) EOC-S experiments. Points indicate measured δ18Oc from published studies (Materials and methods and Table 1), consisting of the most negative values for lacustrine carbonates and average values for paleosols. (C) The sum of squared residuals (SSR) for the Eocene simulations (y axis) against modeled TP elevation (Materials and methods). The best model-data fit is obtained when SSR is minimized [cases with relatively low Asian topography (EOC-S and EOC-M)].

Fig. 6 Maps of δ18Oc calculated from model simulated MJJAS precipitation-weighted δ18Ow and model MJJAS mean temperatures for Eocene cases.

(A) EOC-L. (B) EOC-Him. (C) EOC-M. (D) EOC-sea. (E) EOC-North. (F) EOC-South. For all Eocene maps, points show measured δ18Oc from published data sets, consisting of the most negative values for lacustrine carbonates and average values for paleosols (Table 1).

Despite uncertainties in our modeling approach (Materials and methods), our simulations highlight (i) the permanence of a latitudinal δ18Ow dipole in agreement with conclusions from recent studies (37, 38) and (ii) the limit of standard atmospheric distillation models based on modern calibrations. We show that most of the assumptions that underlie stable isotope paleoaltimetry—such as commonly used modern temperature lapse rates, single and invariant moisture sources, and similar convective precipitation activity—are not valid for the Eocene TP. By contrast, our simulations suggest that previously published δ18Oc are more compatible with modest Eocene TP elevations (not higher than 3000 m). These estimates reconcile stable isotopic paleoaltimetry estimates with alternative paleontological (13, 14) and paleobotanical (11, 12) data. Our results imply that substantial uplift (>2000 m) of at least some parts of the TP could have occurred after the Eocene. The complexity of atmospheric processes within greenhouse climates, combined with differing paleogeographies, highlight the necessity of using dedicated climate simulations for interpreting δ18Oc data.

Materials and methods

Numerical modeling approach

In order to assess spatial and temporal variations of climatic parameters and the isotopic composition of vapor and precipitation, we use the LMDZ-iso Atmospheric General Circulation model (AGCM) (39), a derivative of LMDZ (40). LMDZ was developed at Laboratoire de Météorologie Dynamique (LMD), Paris, France, and is the atmospheric component of the Institute Pierre Simon Laplace (IPSL) Earth System Model (ESM) (40), which participates in the Coupled Model Intercomparison Project (CMIPs) program (41). LMDZ incorporates many processes decomposed into a dynamical core, calculating the numerical solutions of general equations of atmospheric dynamics, and a physical part, calculating the details of the climate in each grid point and containing parameterizations for processes such as the effects of clouds, convection, and orography.

In the version of LMDZ used in this work, the land surface is represented as a simple bucket model, and land surface evaporation is calculated as a single flux: no distinction is made between transpiration, bare soil evaporation, or evaporation of intercepted water by the canopy. However, despite these limitations, the ability of this model to simulate atmospheric dynamics and reproduce rainfall patterns over Asia has been demonstrated in numerous studies (20, 24, 4244).

LMDZ-iso represents isotopic processes associated with all phase changes and is documented in (39). The implementation of water isotopes in the convection scheme is described in (45). Water species are transported and mixed passively by large-scale advection (Van Leer advection scheme (46)) and various air mass fluxes. Equilibrium fractionation coefficients are calculated after Merlivat and Nief (47) and Majoube (48) and kinetic effects are taken into account following Merlivat and Jouzel (49) and Jouzel and Merlivat (50). The model takes into account the evolution of the composition of both the rain and surrounding vapor as raindrops re-evaporate (45). However, the model does not account for raindrop size during evaporation as described in Lee and Fung (51).

Here, we use a model configuration with 144 grid points in longitude, 142 in latitude, 39 vertical layers, and a stretchable grid (40) that permits greater regional spatial resolution and yields an average resolution of ~50 km over central Asia (fig. S10). We force LMDZ-iso with realistic sea-surface temperatures (SSTs) simulated by FOAM (52) and by realistic albedo and roughness values, obtained with the vegetation distribution simulated by the biophysics vegetation model LPJ (fig. S11) (53). We choose this approach rather than (i) using LMDZ-iso with pre-industrial boundary conditions because it permits accounting for the greenhouse gases concentration and related changes in SSTs and (ii) use of a fully coupled GCM, because stable isotopes have not been implemented yet in the IPSL ESM, and present-day calculation capabilities do not allow performing fully coupled experiments in high resolution with the isotopes implementation. We should note, however, that FOAM and similar coarse resolution models (for example, MPIOM, T31) have been successfully applied in many other paleoclimate studies (54). In order to show the capability of the FOAM model to simulate climate, we provide here a rough comparison between the SST and the field data [see the data compilation in (55)] for the Late Eocene, from 38 to 34 Ma (fig. S12). As in similar recent climate model studies (56), we find that our simulation fits some of the data but not all and that our latitudinal gradient remains too steep when compared to data. We stress that this is a crude comparison as we include data from a long time interval (Late Eocene) during which global climate changes substantially.

Experimental design

We provide experiments depicting both present-day and Eocene climate. For the Eocene experiments, we use the Eocene paleogeography reconstruction of Licht et al. 2014 (24), which includes a southward shift of the Himalayan front, and a large Paratethys Sea to the north of the TP. To understand the role of topography on Eocene climate and the distribution of δ18Oc, we perform a series of experiments by varying the elevation of the TP. The Eocene TP is as high as 3500 m in the initial reconstruction (EOC-L simulation) (fig. S1D) (24). The elevation is increased to 5000 m in the EOC-XL simulation (fig. S1E), reduced to 1750 m in the EOC-M simulation (fig. S1C) and lowered to 200 m in the EOC-S simulation (fig. S1B). Additionally, we modify the TP geometry to obtain steep, narrow mountains analogous to the present-day Himalayas over the southern flank of the elevated topography in the EOC-Him experiment (fig. S1F). Lastly, we perform a series of complementary experiments to test different possible Eocene paleogeographies for the studied area: (i) with the TP moved to the south, EOC-South experiment (fig. S1G); (ii) to the north, EOC-North experiment (fig. S1H); and (iii) with a water body over northern India, EOC-sea experiment (fig. S1J).

For all Eocene simulations, we use Pco2 of 1120 ppm, which lies within published estimates (21), and identical SSTs to isolate the effect of topography on atmospheric dynamics and δ18Ow. We perform one additional experiment with the Eocene paleogeography with pre-industrial Pco2 (EOC-1x experiment) in order to estimate its impact on resulting δ18Oc.

We use present-day orbital parameters, which is a standard technique for deeptime paleoclimate modeling studies, given that we provide snapshots to illustrate the climatic state during the Eocene period and do not focus on representing transient climates due to orbital variations. However, in some other studies, it has been shown that such variations may have a significant impact on the hydrological cycle over Asia (24, 25). In order to address this issue, we have tested the impact of orbital parameters variations on the isotopic composition of precipitation with a series of additional paleo-climate experiments (fig. S4, B and C, EOC-CO and EOC-WO experiments). In these experiments, we test two extreme orbital configurations for the Eocene case in order to produce either (i) “warm austral” orbit (or “cold boreal” orbit), with eccentricity 0.05, obliquity 24.5, and Earth at perihelion in January (EOC-CO) and (ii) a cold austral summer (and thus a “warm boreal” summer), with eccentricity 0.05, obliquity 22.5, and Earth at perihelion in July (EOC-WO). Despite some differences in simulated MJJAS precipitation amounts for these “warm boreal” and “cold boreal” scenarios, δ18Op summer [May, June, July, August, and September (MJJAS)] patterns appear to be fairly similar for the three orbital forcings tested here. For all cases, more negative values of δ18Op are simulated over India and the southern slopes of the HTP, while an area between 20°N and 40°N is characterized by more positive values of δ18Op. This gradient is strongest for the “warm boreal” scenario case.

Model-data comparison

We compare Eocene precipitation-weighted δ18Oc derived from our simulations with published Eocene δ18Oc data (Table 1 and Figs. 5 and 6) (8, 10, 13, 24, 5766). Pedogenic carbonates on the TP commonly form during the summer season (67); thus, to derive δ18Oc from our simulated δ18Ow values, we use May–September mean surface temperatures and use the fractionation factors of (17).

With an objective of adjusting the parameters of a model function [y = f(x)] to best fit a data set (where x is δ18Oc from published data sets and y is δ18Oc simulated by LMDZ-iso), we calculate the sum of squared residuals (SSR) for 20 data points, which permits us to define the optimum model-data fit, which is when the SSR is minimized. For EOC-S and EOC-M simulations, SSR is estimated to be the smallest, 65.8 and 160 respectively, while EOC-XL and EOC-Him simulations show the highest values, 369.7 and 453.63 respectively. SSR for the EOC-sea experiment is comparable with those for the EOC-L experiment and both EOC-North and EOC-south show relatively high SSR similar to the EOC-Him case. Note that for SSR calculation for the EOC-North and EOC-South, we have adjusted the location of the data points 5° to the south and to the north respectively. The difference in the SSR between two different elevations scenarios is larger than the difference attributable to the variability in the δ18O dataset (Table 1). Consequently, we dismiss variability in the carbonate δ18O as excluding the firm conclusions that the best data-model match is with lower elevation scenarios.

Uncertainties and limits of our modeling approach

A major uncertainty in this study is due to the reconstruction of the geographical position of each site’s δ18Oc data during Eocene. Based on paleomagnetic data for the Hoh Xil Basin (modern coordinates: 34.60°N, 93°E) (68), we reconstruct its Eocene paleolatitude as approximately 21° N (as a result of northward convergence of the Hoh Xil Basin with respect to Eurasia (Siberia) since Early Eocene–Late Oligocene time). Thus, we use a constant correction of 13° for positioning all of the geologic δ18O sampling sites over our Eocene paleogeography. Such paleoaltitude estimates and corrections of about 13 ± 5° are in general consistent with the palaeolatitudes derived from a large number of high-quality palaeomagnetic studies from eastern China to Kyrgyzstan (23).

Though GCMs produce realistic Eocene temperatures over Asia and simulate the associated response of δ18Ow, numerical modeling with GCMs involves several uncertainties, including the assumption of constant seawater δ18O and the accuracy of the paleogeographic reconstruction. The latitudinal position of the Indian subcontinent in the Eocene and the height of other orogens in Asia, as well as their spatial extent, remain a point of contention (23, 69) and may impact the results of our experiments. To address this limitation, additional paleogeographic scenarios should be tested in future work. Other potential limitations of using isotope-enabled GCMs include complexity associated with our representation of soils, evapotranspiration, and vegetation. However, results of our current study permit a first-order estimate of Eocene δ18Ow and its low sensitivity to topographic change.

Decomposing δ18Oc differences

Our goal is to understand to what extent δ18Ow changes explain the δ18Oc signal over the Eocene TP and what part of this signal depends on the temperature of carbonate formation. First, based on an assumption of equilibrium between water and calcite, we relate δ18Oc to δ18Ow using the relationship of (17)Rc = Rw/1.0309 – 29.98 + 17.49*1000/T – 31.45(1)where T is the temperature of carbonate formation, Rc is the isotopic ratios in carbonate, and Rw is the isotopic ratios in paleo surface water, or in a simpler formRc = Rc (Rw, T)(2)In order to understand why Rw varies from one climatic state to another, we decompose ΔRc = Rc2Rc1 into contribution from dRw and dT terms, referring to the climatic states using subscript 1 and 2 and to their difference using the Δ notationΔRc = ΔRc,dRw + ΔRc,dT + N(3)where ΔRc,dRw and ΔRc,dT are the contributions of the part of the Rc signal related to the change of Rw and temperature accordingly. Nonlinear terms of decomposition are gathered into the residual term N.

Similar to the isotopic signal decomposition suggested in (20), contributions could be estimated asΔRc,dRw = Rc (Rw2, dT′) – Rc (Rw1, dT′)(4)andΔRc,dT = Rc (dRw′, T2) – Rc (dRw′, T1)(5)where Rw2, Rw1, T2, and T1 are δ18Ow and temperature from two climatic states, respectively, and dRw′, dT' are centered differencesdRw′ = (Rw1 + Rw2)/2(6)dT′ = (T1 + T2)/2(7)From the Eqs. 1, 3, 4, and 5, we have

ΔRc,dRw = (Rw2Rw1)/1.0309(8)ΔRc,dT = 17.49 * (1000/T2 – 1000/T1)(9)N = ΔRc – (ΔRc,dRw + ΔRc,dT)(10)

Supplementary Materials

References and Notes

Acknowledgments: Funding: This work was funded by the FP7-PEOPLE: Marie-Curie Actions [grant 316966, iTECC (interaction Tectonics-Erosion-Climate-Coupling) project]. This work was granted access to the HPC resources of IDRIS under the allocation 2017-0292 made by GENCI. S.B. was also supported by the German Ministry for Education and Research (BMBF) grant 03T0863A to T. Ehlers. J.K.C.R. is funded by an ETH Fellowship. Author contributions: S.B., P.S., and Y.D. developed the conceptual idea. S.B. designed and carried out simulations in close collaboration with P.S., C.R., and Y.D.; J.K.C.R. and S.B. did statistical analysis. S.B., A.L., and J.K.C.R. analyzed observed carbonate data. All authors discussed the results and implications and commented on the manuscript. Competing interests: None declared. Data and materials availability: All data are available in the manuscript or in the supplementary materials.

Stay Connected to Science

Navigate This Article