How Can Climate Models Be Used in Paleoelevation Reconstructions?

Introduction

Quantitative paleoelevation techniques based on proxy data include approaches such as stable isotopes (e.g., δ¹⁸O, δD) measured in terrestrial deposits, paleofloral and faunal findings and their physiognomic characteristics, and clumped isotope (Δ₄₇) paleothermometry. These techniques (and others) have been extensively applied to various mountain ranges such as the Himalayas and Tibetan Plateau (Rowley and Currie, 2006; Gébelin et al., 2013; Ding et al., 2014), the North America Cordillera (Chamberlain et al., 2012; Gébelin et al., 2012; Cassel et al., 2014), the Andes and Andean Plateau (Mulch et al., 2010; Garzione et al., 2017) for elevation reconstructions throughout the Cenozoic, ranging from the Paleocene and Eocene (e.g., Ding et al., 2014) up to the Pleistocene (e.g., Hoke et al., 2014). Oxygen isotope paleoaltimetry is one of the most widely applied techniques and is based on a progressive decrease in the δ¹⁸O value of precipitation (δ¹⁸O_p) with elevation during the topographic (i.e., orographic) ascent of an air mass. This technique assumes Rayleigh distillation of an air mass under equilibrium conditions and without local evaporation or advection of additional vapor from outside the system (Rozanski et al., 1993; Gat, 1996).

Oxygen isotope paleoaltimetry uses δ¹⁸O preserved in pedogenic or lacustrine carbonates (hereafter referred to as δ¹⁸O_c) as a proxy for paleoelevation. Interpretation of measured δ¹⁸O_c generally requires two steps. For the first step, the reconstruction of paleoprecipitation δ¹⁸O (δ¹⁸O_pp) using the calcite-water fractionation coefficients depend on the carbonate formation temperature as described through the equilibrium between calcite and water δ¹⁸O (Kim and O’Neil, 1997). Paleoaltimetry estimates based on δ¹⁸O alone, without independent information on paleotemperatures from other proxies or a paleoclimate model, depend on assumptions made about the paleoelevation and temperature (global and regional) changes relative to the modern. Knowledge of these paleo-conditions is challenging to come by. Moreover, if a surface temperature is assumed, then the calculation of paleoelevation from the proxy involves circular reasoning. This is because the temperature used at the point of interest depends on the elevation through the temperature lapse rate and often the chosen coefficients from the equations in Kim and O’Neil (1997) already implies high elevations.

For the second step, the reconstructed δ¹⁸O_pp must be attributed to a paleoelevation using a δ¹⁸O_pp-elevation relationship (the “isotopic lapse rate”) in a reference to low-elevation δ¹⁸O. A common approach (largely due to it’s simplicity) is to assume that the modern isotopic lapse rate measured from precipitation in an elevation-transect of stations (e.g., Gonfiantini et al., 2001; Fiorella et al., 2015), or simulated using a one-dimensional thermodynamic model for δ¹⁸O_p composition based on the Rayleigh distillation, can be used to predict changes in the δ¹⁸O_p response to changes in temperature, and the relative humidity of air parcels during orographic ascent (Rowley et al., 2001; Rowley and Garzione, 2007). This approach assumes the isotopic lapse rate at the time of carbonate formation not to change over millions of years. Thus, both of the previous steps rely on several assumptions regarding the contribution (if any) of global and regional climate change to δ¹⁸O and the stationarity of the δ¹⁸O_pp-elevation relationship throughout geologic time.

However, numerous recent studies have shown that a range of atmospheric processes are important for the interpretation of δ¹⁸O_pp. For example, hemispheric-scale atmospheric circulation and associated teleconnections (Schneider and Noone, 2007; Takahashi and Battisti, 2007; Pausata et al., 2011), moisture transport changes (Ehlers and Poulsen, 2009; Poulsen et al., 2010), evapotranspiration and water vapor recycling within continental interiors (Risi et al., 2013; Chamberlain et al., 2014), and shifts in precipitation amount and ratios of convective-to-large-scale precipitation (Lee and Fung, 2008; Feng et al., 2013; Botsyun et al., 2019a) all can affect δ¹⁸O_pp. The effect of these complicating processes on isotopic lapse rates range from insignificant to major, up to the extent of flattening or inverting the paleo isotopic lapse rate relative to the modern (e.g., Moran et al., 2007). Thus, despite the intriguing potential of proxy records to reconstruct the paleoelevation history of orogens, the interpretation of the data requires several assumptions about the congruence between modern and paleo atmospheric processes and climate which may, or may not, be correct (e.g., Starke et al., 2020). Therefore, a rigorous interpretation of paleoelevations requires evaluating these assumptions on a case-by-case basis, which in practice is not always simple. One promising methods to address these concerns is an application of climate models for reconstructing time-specific temperatures, δ¹⁸O_p and isotopic lapse rates for adequate interpretation of δ¹⁸O_c signal.

Basic Principles of Isotope Tracking Climate Models… in a Nutshell

Climate (and paleoclimate) models are typically based on a set of governing primitive equations for conservation of mass, energy, and momentum. These equations are solved simultaneously by discretizing (gridding) the surface of the earth and atmosphere and applying a numerical model. If only an atmospheric General Circulation Model (GCM) is used (as opposed to a coupled ocean-atmosphere GCM) then solving these equations also requires application of boundary conditions such as sea-surface temperatures, solar radiation, and also material properties of the atmosphere and land surface such as greenhouse gas concentrations, topography, vegetation/ice cover, sea-ice extent, soil properties, etc. Given the physics-based approach of these models, they can be used for diverse purposes including weather prediction/forecasts, future climate predictions, or paleoclimate predictions so long as the appropriate boundary conditions, and land surface properties are known. Thus, paleoclimate modeling studies that are applied to paleoelevation proxy data interpretation are “simply” applying well used and tested climate models from the atmospheric sciences community, with the caveat that they require a means for predicting the isotopic composition of precipitation, discussed next.

Isotopic Modeling and (Non-)stationarity of δ¹⁸O_p-Lapse Rate

Climate modeling has been extensively applied to study the impact of mountain surface uplift on atmospheric physics and dynamics (Ruddiman and Kutzbach, 1989; Broccoli and Manabe, 1992; Sepulchre et al., 2006). The implementation of water isotope tracking in simple climate models (Dansgaard, 1964) and further into GCMs (Joussaume et al., 1984) enabled calculation of the atmospheric, land surface, and topographic processes impacting δ¹⁸O_p. Since the development of iGCMs, a diverse range of studies have been conducted documenting the sensitivity of δ¹⁸O_p to climate change resulting from variations in pCO₂ (Poulsen et al., 2007; Poulsen and Jeffery, 2011), sea surface temperatures (SSTs) (Sturm et al., 2007), sea level variations (Poulsen et al., 2007), and paleogeography (Sewall and Fricke, 2013; Roe et al., 2016; Botsyun et al., 2019a), and much more (Risi, 2009; Sturm et al., 2010 and references herein).

The previous iGCM based studies have shown, among other things, that present-day isotopic lapse rates are not always suitable for paleoelevation reconstructions. First, modeling studies have shown spatial non-stationarity of the isotopic lapse rates. For example, Galewsky (2009) uses an idealized, fully nonlinear atmospheric model to evaluate the effects of stratified atmospheric flows on orographic precipitation isotopic ratios, with the aim of improving the foundation for interpreting proxy records used for paleoelevation studies. He showed that existing Rayleigh condensation models are strictly applicable for paleoaltimetry studies only when the atmosphere is not stratified. Changes in climate or the horizontal terrain aspect ratio can change precipitation isotopic ratios at least as much as changes in surface elevation (Galewsky, 2009), suggesting that proxy records must be interpreted within a broad context of climate variability and landscape evolution.

Second, iGCMs have demonstrated the impact of global climate change on δ¹⁸O_p and the non-stationarity of isotopic lapse rates through time. For example, in climate simulations with elevated CO₂ concentrations, the middle troposphere undergoes preferential moistening and warming, which weakens vertical stratification and gives rise to shallower lapse rates (Poulsen and Jeffery, 2011). Li et al. (2016) identified spatial and temporal variations in simulated isotopic lapse rates for different high relief zones along the flanks of the Tibetan plateau during Late Quaternary climates. They found that the precipitation weighted annual mean δ¹⁸O_p lapse rate at the Himalaya is about 0.4‰/km larger during the Middle Holocene and 0.2‰/km smaller during the Last Glacial Maximum than during pre-industrial times. These changes are large enough to impact the interpretation of proxy data (Li et al., 2016). Furthermore, results from Botsyun et al. (2020) predict a similar magnitude (4–5‰ relative to present day) of change in δ¹⁸O_p resulting from the topographic development of the European Alps and paleoclimate change during the Pliocene and the Last Glacial Maximum.

Third, modeling studies show that regional (e.g. continental or subcontinental scale) climate changes result from mountain surface uplift and also have a potential impact on δ¹⁸O_p. For example, previous work has shown that Andean Plateau surface uplift resulted in the onset of convective precipitation (Poulsen et al., 2010) following a switch of the prevalent moisture source and associated transport paths from the South Pacific to the Equatorial Atlantic and initiation/strengthening of the South American Low Level Jet (Ehlers and Poulsen, 2009; Insel et al., 2010). These studies suggest that high precipitation rates enhance the isotope amount effect, leading to more negative δ¹⁸O_p at high elevations and, thus, increases in the isotopic lapse rate. In agreement with this result, experiments with reduced elevations of the Andean Plateau by Insel et al. 2012 show that use of present-day isotopic lapse rate can lead to an underestimation of surface elevation by 2,000 m.

Regional climate changes associated with surface uplift have also been shown to impact δ¹⁸O_p across the North American Cordillera (Feng et al., 2013) and Tibetan Plateau (Botsyun et al., 2016; Shen and Poulsen, 2019). For both these areas it was suggested that neither the common assumption that isotopic fractionation occurs primarily through rainout following Rayleigh distillation, nor the application of modern empirical δ¹⁸O_p lapse rates to past environments, are valid. Shifts in atmospheric processes, including shifts in local precipitation type between convective and large-scale rain and between rain and snow, intensification of low-level vapor recycling particularly on leeward slopes, development of air mass mixing and changes in wind direction and moisture source changes all contributed to the deviations of the isotopic signal.

Finally, the robustness of applying the present-day isotopic lapse rate to interpret Eocene Tibetan Plateau proxy data was recently tested in work of Botsyun et al. (2019a). In this study, an iGCM with time specific Eocene boundary conditions was applied to assess the influence of changing Eocene paleogeography on climate and δ¹⁸O_p signals. The authors found that a combination of increased convective precipitation, a mixture of air masses of different origin, widespread aridity, and an intensified water recycling resulted in a reversed isotopic lapse rate across the southern flank of the Tibetan Plateau. These processes resulted in the most negative δ¹⁸O_p occurring over northern India and increased δ¹⁸O_p northward (Figure 1). Taken together, these results indicate that standard stable isotope paleoaltimetry methods are not applicable in Eocene Asia and an alternative method, which takes into account paleo isotopic lapse rates, should be developed.

FIGURE 1. Schematic representation of the processes in the atmosphere and over land that potentially contribute to a δ¹⁸O_p signal resulting in: (A) an isotopic lapse rate according to a Rayleigh distillation model for an isolated air parcel; and (B) flattened or inversed isotopic lapse rate. Examples of proxy data interpretations based on iGCM results.

In recent years, iGCMs have been extensively applied for δ¹⁸O proxy interpretation by using model-derived: 1) δ¹⁸O_p lapse rates for various elevation scenarios; and 2) δ¹⁸O_c and soil temperatures for time-specific simulations including various paleogeographies and paleotopographies. Fan et al. (2017) and Gao and Fan (2018) used an iGCM-derived δ¹⁸O_p-elevation gradient from Feng et al. (2013) to constrain the paleoelevations of the Cordillera orogenic system (Uinta Mountains, United States). This approach allows them to consider vapor mixing in the mountain flank and the warm climate during the Paleogene. Sundell et al. (2019) interpreted observed δ¹⁸O_p across the Peruvian central Andes using different lapse rates including the thermodynamically-derived non-linear model from Rowley et al. (2001), an empirical approach from Quade et al. (2007), and reconstructed isotopic lapse rate using a climate model. Experiments with 25, 50, 75, and 100% of Andean elevations from Insel et al. (2012) and experiments with 50 and 100% Andean elevation from Poulsen et al. (2010) together with climate correction are used for the data interpretation.

Botsyun et al. (2019a) designed experiments with variable Tibetan Plateau elevations and Eocene paleogeography to compare simulated δ¹⁸O_c values with previously published data. These sensitivity tests identified that the best model-data fit corresponds to a low Tibetan Plateau in the Eocene. The authors showed that quantitative estimates of paleoaltimetry from stable isotopes depend on the outcome of a climate model, which in turn used paleoaltimetry, as input. This means that paleoaltimetry and paleoclimate models need to be iteratively tuned until the input elevation and the output relations between isotopic ratio and elevation are mutually consistent (Figure 2).

FIGURE 2. Principal scheme showing application of a paleoclimate model for paleoelevation reconstructions. t.f.—transfer function.

Limits of climate modeling technique to paleoaltimetry studies

From a stable isotope paleoaltimetry point of view the application of iGCMs have several limitations. The key caveats of climate modeling approach are linked to: 1) the limited ability of the model to simulate δ¹⁸O patterns over regions with complex topography due to the complexity of natural processes and limited scope of model parameterizations; 2) uncertainties in precipitation and δ¹⁸O changes in monsoonal areas; 3) limited computational resources for high-resolution fully-coupled ocean-atmosphere simulations; 4) uncertainties in the choice of boundary conditions (e.g. Figure 2), 5) relatively low model spatial resolution, and 6) a “social” challenge whereby the most fruitful applications involve a collaboration between proxy-oriented and model-oriented scientists.

Examples supporting the first point are as follows. Convective processes are known to have a significant imprint on the stable oxygen isotope composition of precipitation, especially in the tropics (Risi et al., 2008). However, convective processes in a iGCM strongly depend on the parameterization applied (Hourdin et al., 2006). Large uncertainties also exist in iGCMs depending on the representation of water stable isotopes parameterizations used for the isotopic exchanges between vapor and rain droplets during descent and the degree of partial reevaporation (Lee and Fung, 2008; Risi et al., 2010). Risi et al. (2010) demonstrated a sensitivity of modeling results to evaporation, especially in dry regions where this process is prevalent. In addition, the isotopic equilibration time of raindrops is drop-size dependent, with smaller equilibration time for small raindrops (Lee and Fung, 2008). Given this, if raindrop size is not taken into account an underestimation of equilibration times for heavy rains with big raindrops could be expected. Additional potential limitations of using iGCMs include complexity associated with realistic representation of soils, evapotranspiration schemes, geochemistry, vegetation, and clouds (Risi et al., 2010; Haese et al., 2013; Rio and Hourdin, 2015). These challenges are active areas of research in the climate modeling community and improvements in model parameterizations are always occurring.

The response of the hydrologic cycle to the climate forcing in monsoon areas vary depending on the model applied and, even for future climates, monsoons amplification has not been shown as a strong feature (Wang et al., 2020). Furthermore, the interpretation of water isotope variations and their link to monsoon is not straightforward even for recent climates (e.g., Middle Holocene) when there is no problem of mountain elevation change (Lin et al., 2019). Further uncertainties are associated with the application of atmosphere-only climate model with no feedbacks between the atmosphere, vegetation and the ocean. Use of modeling techniques with SSTs and sea ice extension prescribed from coupled simulations to run a separate high resolution atmospheric model suffer from not taking into account ocean variability and its feedback to the atmosphere. Moreover, vegetation can also play an important role in isotopic changes through not only evapotranspiration changes, but also through changing surface albedo feedback. Given this coupled high-resolution ocean-atmosphere and dynamic vegetation modeling approaches are preferred, but extremely challenging to conduct. However, usage of high-resolution atmosphere-only iGCMs that use previous coarser resolution coupled ocean-atmosphere models for boundary conditions have been shown to be a good trade-off between the high-resolution necessary to capture observed spatial, seasonal, and daily variations of δ¹⁸O_p and computation resources required (Botsyun et al., 2019b).

Furthermore, paleoclimate studies are sensitive to the paleogeographical reconstruction used (e.g., Baatsen et al., 2016). Producing paleogeographic reconstructions for climate models is time-consuming and elaborate. For example, paleoclimate models frequently use reconstructions where the latest state-of-the-art of plate tectonic, crustal and paleomagnetic reconstructions, paleotopography, paleobathymetry, or vegetation distribution have not yet been incorporated. The availability of realistic Cenozoic boundary conditions, such as sea-surface temperatures, is also limited, unless a time-consuming coupled ocean-atmosphere model is used. Finally, coarse model resolutions used in iGCMs to reduce simulation time often prevent accurate representation of smaller-scale topographic features such that explicit simulation of air masses ascending mountain topography are only grossly represented (Hourdin et al., 2006; Haese et al., 2013). Despite the above caveats, numerous studies have applied iGCMs to capture regional changes in δ¹⁸O_p through time, and comparison of model predicted δ¹⁸O_p to observations (e.g., Global Network of Isotopes in Precipitation (GNIP) stations data) to understand model limitations are commonly conducted when applying models to new locations.

When are state-of-the-art climate models useful and when can simpler approaches be applied?

Despite the previously mentioned limitations associated with iGCMs, applications of climate models to enhance proxy interpretations are gaining success (e.g., van Hinsbergen and Boschman, 2019). Modeling studies have shown that simpler (i.e., non-iGCM based) paleoaltimetry approaches are valid when the Rayleigh distillation process is dominant. However, iGCMs have provided several examples of when atmospheric processes lead to δ¹⁸O_p variations that are different, and sometimes opposite, of those estimated using Rayleigh distillation models of moist adiabatic condensation. For example, Shen and Poulsen (2019) used the ECHAM5-wiso model to show that Rayleigh distillation is only prevailing in the monsoonal regions of the Himalayas when the mountains are high. In contrast, when orogen topography is lower, local surface recycling and convective processes become important. This is due to weakened forced ascent causing weaker Asian monsoons. In this example, applying a standard paleoaltimetry approach is not valid. Similarly, over the North America Cordillera changes in atmospheric processes have been shown to impact δ¹⁸O_p and violate the common assumption that isotopic fractionation occurs through rainout following Rayleigh distillation (Feng et al., 2013). In this example, the atmospheric changes included shifts in local precipitation type between convective and large-scale rain and also between rain and snow, intensification of low-level vapor recycling on leeward slopes, and development of air mass mixing and changes in wind direction and moisture source. Given the uncertainty of when, or when not, a one-dimensional Rayleigh distillation model adequately represents moisture transport, the significance of lateral vapor mixing, or temporal changes in the isotopic lapse rate, the application of iGCMs is recommended (e.g., Li et al., 2016) as an initial step to characterize these processes rather than assuming they are insignificant.

There are, however, observational (non-iGCM based) techniques that have emerged to help circumvent the effects of changing atmospheric processes on paleoaltimetry data. One such technique is the δ-δ approach. This approach assumes that a baseline sets of samples collected from a low elevation stratigraphic succession provides a reference frame for interpreting higher elevation samples. The implicit assumption is that any regional or global climate change will affect stable isotopes records from both (low and high elevation) localities equally if looking at the same time interval. If true, then the difference between the two records can be interpreted as resulting from elevation change, thereby enabling a paleoelevation reconstruction (e.g., Campani et al., 2012; Gébelin et al., 2012; Mulch, 2016; Pingel et al., 2020). While the δ-δ approach holds promise and there are geographic locations where this may work well, the approach does require that a suitable baseline location has been chosen–which is difficult to do without prior knowledge. In recent work by Botsyun et al. (2020), an iGCM was used to explore the δ¹⁸O_p signal of Cenozoic topographic and paleoclimatic change in the European Alps. The results of their study suggest that the δ-δ approach is applicable for the case of the European Alps under the condition that the low-elevation reference section measured hasn’t experienced climate change, and that it is located far enough away and upwind from the area of interest to not experience climate changes associated with orogenic uplift.

In summary, when working with proxy records from “non-simple” locations where Rayleigh distillation assumptions may not be valid, the modern paleoelevation community is well equipped to make a step forward to more advanced reconstructions involving isotope enabled climate modeling. However, doing this requires a collaboration between proxy and climate modeling communities to evaluate how simple, or complex, the δ¹⁸O_p response is to changing atmospheric processes. In this review, we have highlighted how application of climate models in conjunction with geological/geochemical data provides a powerful tool to incorporate climatic change effects into the analysis of paleoaltimetry work. Moreover, paleobotanic data and stable oxygen isotope paleoaltimetry can be reconciled using climate models. Further model-data comparison studies that combine multiple proxy types, together with iGCM modeling efforts testing the impact of complex topography structures as well as atmospheric parameterizations on δ¹⁸O are paving the way for more accurate and robust paleoelevation estimates in the decades to come.

Author Contributions

SB and TAE conceived and wrote the manuscript.

Funding

This study was supported as part of the German priority research program Mountain Building Processes in 4D (MB4D) grant (EH329/18-1 and 23-1) to TAE.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

References