How Can Climate Models Be Used in Paleoelevation Reconstructions?

Paleoelevation reconstructions derived from proxy data such as stable oxygen isotope records in terrestrial archives have been determined for Cenozoic mountain ranges around the world. Recent studies have highlighted that a variety of paleoclimate processes can contribute to the isotopic composition of a measured precipitation (δ18Op) signal used in elevation reconstructions. These processes can include: regional, global, and topographic variations in paleotemperature; environmental conditions of an air mass before orographic ascent; evapotranspiration; water vapor recycling; and changes in the vapor source. In some cases, these processes can overprint the elevation signal sought in proxy data and preclude robust elevation reconstructions. Recent advances in isotope tracking climate models allow us to estimate paleoclimate changes during orogen development and associated changes in paleo δ18Op due to both climate and topographic changes. These models account for adiabatic and non-adiabatic temperature changes, relative humidity variations, changing continental evapotranspiration, vapor recycling, vapor source changes, etc. Modeling strategies using high-resolution isotopes-enabled General Circulation Models (iGCMs) together with time-specific boundary conditions and variable topography provide a powerful tool for enhancing elevation reconstructions from δ18Op proxy data. In this review, we discuss the principles, benefits and caveats of using iGCMs for interpreting isotopic records from natural archives for paleoelevation reconstructions. We also highlight future challenges for the application of iGCMs to paleoaltimetry proxy data that open up new avenues for research on tectonic-climate interactions.


INTRODUCTION
Quantitative paleoelevation techniques based on proxy data include approaches such as stable isotopes (e.g., δ 18 O, δD) measured in terrestrial deposits, paleofloral and faunal findings and their physiognomic characteristics, and clumped isotope (Δ 47 ) 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 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 δ 18 O value of precipitation (δ 18 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 δ 18 O preserved in pedogenic or lacustrine carbonates (hereafter referred to as δ 18 O c ) as a proxy for paleoelevation. Interpretation of measured δ 18 O c generally requires two steps. For the first step, the reconstruction of paleoprecipitation δ 18 O (δ 18 O pp ) using the calcite-water fractionation coefficients depend on the carbonate formation temperature as described through the equilibrium between calcite and water δ 18 O (Kim and O'Neil, 1997). Paleoaltimetry estimates based on δ 18 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 paleoconditions 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 δ 18 O pp must be attributed to a paleoelevation using a δ 18 O pp -elevation relationship (the "isotopic lapse rate") in a reference to lowelevation δ 18 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 δ 18 O p composition based on the Rayleigh distillation, can be used to predict changes in the δ 18 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 δ 18 O and the stationarity of the δ 18 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 δ 18 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 δ 18 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 timespecific temperatures, δ 18 O p and isotopic lapse rates for adequate interpretation of δ 18 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 seasurface 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 physicsbased 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 δ 18 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 δ 18 O p . Since the development of iGCMs, a diverse range of studies have been conducted documenting the sensitivity of δ 18 O p to climate change resulting from variations in pCO 2 (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 δ 18 O p and the non-stationarity of isotopic lapse rates through time. For example, in climate simulations with elevated CO 2 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 δ 18 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 preindustrial 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 δ 18 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 δ 18 O p . For example, previous work has shown that Andean Plateau surface uplift resulted in the onset of convective precipitation  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 δ 18 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 δ 18 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 δ 18 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 δ 18 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 δ 18 O p occurring over northern India and increased δ 18 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.
In recent years, iGCMs have been extensively applied for δ 18 O proxy interpretation by using model-derived: 1) δ 18 O p lapse rates for various elevation scenarios; and 2) δ 18 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 δ 18 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 δ 18 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 δ 18 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 Frontiers in Earth Science | www.frontiersin.org February 2021 | Volume 9 | Article 624542 between isotopic ratio and elevation are mutually consistent ( Figure 2).

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 δ 18 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 δ 18 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 oceanatmosphere 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 δ 18 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 timeconsuming 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 δ 18 O p through time, and comparison of model predicted δ 18 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 δ 18 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 δ 18 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 Frontiers in Earth Science | www.frontiersin.org February 2021 | Volume 9 | Article 624542 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 δ 18 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 "nonsimple" 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 δ 18 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 δ 18 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.