Opportunities and Challenges for Paleoaltimetry in “Small” Orogens: Insights From the European Alps

Many stable isotope paleoaltimetry studies have focused on paleoelevation reconstructions of orogenic plateaus such as the Tibetan or Andean Plateaus. We address the opportunities and challenges of applying stable isotope paleoaltimetry to “smaller” orogens. We do this using a high‐resolution isotope tracking general circulation model (ECHAM5‐wiso) and explore the precipitation δ18O (δ18Op) signal of Cenozoic paleoclimate and topographic change in the European Alps. Results predict a maximum δ18Op change of 4–5‰ (relative to present day) during topographic development of the Alps. This signal of topographic change has the same magnitude as changes in δ18Op values resulting from Pliocene and Last Glacial Maximum global climatic change. Despite the similar magnitude of the isotopic signals resulting from topographic and paleoclimate changes, their spatial patterns across central Europe differ. Our results suggest that an integration of paleoclimate modeling, multiproxy approaches, and low‐elevation reference proxy records distal from an orogen improve topographic reconstructions.

Oxygen isotope paleoaltimetry is based on the relationship between the oxygen isotopic ratio of rainfall (δ 18 O p ) and elevation (Rowley et al., 2001). Relationships between δ 18 O p values and elevation have been determined both empirically (e.g., Poage & Chamberlain, 2001) and theoretically using thermodynamic principles (Rowley et al., 2001). However, despite extensive application, this method has several uncertainties due to the sensitivity of both δ 18 O p values and the δ 18 O p -elevation gradient (isotopic lapse rate (ILR)) to regional and global climate change and topographic changes that can deflect air masses, induce adiabatic and nonadiabatic temperature changes, relative humidity variations, thereby changing continental evapotranspiration and vapor recycling (Botsyun et al., 2016(Botsyun et al., , 2019Ehlers & Poulsen, 2009;Insel et al., 2012;Mulch, 2016;Poulsen et al., 2010).
A growing amount of quantitative stable isotope paleoaltimetry studies are available for orogens smaller than the Andes and Tibet, including studies of the Pyrenees (Huyghe et al., 2012), the Tauride Mountains of Anatolia (Meijers et al., 2018), the Southern Alps of New Zealand (Chamberlain et al., 1999;Zhuang et al., 2015), the Cascade Mountains (Methner et al., 2016;Takeuchi & Larson, 2005), and the European Alps (Campani et al., 2012;Kocsis et al., 2014;Sharp, 2005). Here we test the hypothesis that if the topography of the Alps changed over the Cenozoic then a δ 18 O p signal of this elevation change can be recovered from paleoaltimetry studies, independent of global climate change effects. While ongoing debate surrounds the timing of surface uplift in the European Alps (see below), resolving this uncertainty is beyond the scope of this study. Instead, we focus on the first step needed to address this debate and evaluate how much signal of topographic change from δ 18 O p could exist relative to global climate-driven changes in δ 18 O p , and where future proxy samples of δ 18 O p should be collected.

Timing and Rate of Surface Uplift in the European Alps
Geodynamic reconstructions show that the European Alps experienced orogen-parallel extension during the Late Oligocene-Early Miocene, promoting rapid exhumation, at least in the Central (Steck, 2008, and references inside) and Eastern segments (Behrmann, 1988), that ultimately resulted in accelerated erosion in the Alpine orogenic wedge during the Oligo-Miocene times (Schlunegger et al., 2007;Willett et al., 2006). Based on sedimentological records, Schlunegger and Kissling (2015) infer high alpine elevations starting at~25 Ma after slab breakoff. However, these estimates contradict, for example, the study of Hergarten et al. (2010) who use results from a geomorphological sediment budget analysis to show that the present-day elevation of the Alps was acquired at~5-6 Ma while the Miocene topography was much lower. This potential tectonic event also correlates with an increase in sedimentation in basins surrounding the Alps starting at~5 Ma (Kuhlemann, 2007). Legrain et al. (2015) based on 10 Be-deriven erosion rates infer an important surface uplift event in the western segment of the Alps at~5 Ma potentially linked to the inversion in the Pannonian Basin. Fox et al. (2015) also report an increase in erosion rate over the past 2 Ma, which was controlled by rock uplift due to the progressive detachment of the European slab under the Western Alps.
Previous paleoaltimetry estimates based on the comparison of hydrogen isotopic ratios (δD) in silicate minerals at high-elevation points of the Central Alps with meteoric water compositions deduced from carbonate-bearing paleosols of the North-Alpine foreland basin, suggest that the central European Alps attained mean elevations similar to those of today no later than the middle Miocene (14.5 Ma; Campani et al., 2012). This result is consistent with Early Miocene paleoelevation estimates of 2,300 + −650 m based on δ 18 O from fossilized shark teeth (Kocsis et al., 2014), but conflicts with estimates of~5,000-m elevation in the Miocene based on fluid inclusions in quartz veins (Sharp, 2005). Coupling high-and low-elevation δ 18 O and δD records from the same time interval (δ-δ approach; Mulch, 2016) yields more reliable paleoaltimetry estimates by disentangling the competing effects of topography and climate on δ 18 O and δD in precipitation (Campani et al., 2012;Gébelin et al., 2013). However, the impact of global climate change, including glacial/interglacial variations, on paleoclimate proxy data is often enigmatic. Depending on the duration and timing (seasonal/cyclic or protracted over longer periods of time during, e.g., slow growth or recrystallization) short-term climatic biases may affect stable isotope paleoclimate and paleoelevation proxy data.
Here we aim at identifying potential sources for this controversy in paleoelevation estimates by (a) quantifying the signal of topographic change that could be preserved in δ 18 O p and (b) evaluating the magnitude of impact of paleotopography and long-term (paleo)climate change have on δ 18 O p .

Materials and Methods
We apply the ECHAM5-wiso atmospheric general circulation model. This model includes isotopologue (H 2 16 O, H 2 18 O, and HDO) tracking . The ability of ECHAM5-wiso to reproduce modern and paleoclimate conditions together with water isotopologue distributions on both global and regional scales has been shown in multiple studies Mutz et al., 2016). We perform ECHAM5-wiso simulations at a T159 spectral resolution (equivalent to a grid spacing of~0.75°or~80 km in latitude and longitude) with a vertical resolution of L31 (31 levels up to 10 hPa), which allows representation of mean alpine topography reasonably well. However, it does not capture high mountain peaks.
First, we provide a control Atmospheric Model Intercomparison Project-style present-day (PD) simulation, with present-day boundary conditions including the Atmospheric Model Intercomparison Project version 2 sea surface temperature and sea ice data from 1957 to 2000 and observed greenhouse gas concentrations for the same period (Nakicenovic et al., 2000). The simulation was conducted for >40 model years. A climatological reference period of 30 years was established for the analysis presented here using the simulation years 1970-1999. This simulation is used to validate the ability of the model to reproduce present-day observed isotopic and climatic patterns.
We present two sets of experiments. The first set includes topographic change sensitivity experiments (Alps150, Alps100, Alps50, and Alps0) conducted to test the climatic and isotopic response to an orogen at 150%, 100%, 50%, and~0% of the present-day alpine elevation ( Figure S1). These simulations isolate the influence of alpine topography on regional climate and δ 18 O p . The control simulation (modern topography, Alps100) was conducted with preindustrial boundary conditions (e.g., insolation, greenhouse gases, sea surface temperatures) and with current (100%) topography of the Alps. Sea surface temperatures and sea ice concentrations are derived from a transient coupled ocean-atmosphere simulation (Lorenz & Lohmann, 2004) with pCO 2 set to 280 ppm. When discussing changes in climate or δ 18 O p values, we refer to deviations from this control experiment. For the three other experiments (Alps150, Alps50, Alps0), we keep all boundary conditions, including albedo, rugosity, and the vegetation distribution identical to those in the Alps100 LGM − Alps100. Isolines on subplot (a) show topography for Alps100 experiment; isolines from 500 m, with a 500-m contour interval. simulation, with the exception of different topographies used in each simulation. We reduce the elevation over the area covering the Alps and alpine forelands to 50% of modern elevations (Alps50) and to 250-m elevation (Alps0). We also conducted an additional simulation with the topography increased to 150% of modern elevations (Alps150).
The second set of experiments includes two simulations, representing climate conditions of the Pliocene (PLIO;~3 Ma) and the Last Glacial Maximum (LGM;~21 ka). These simulations were conducted to estimate the magnitude of global climate change over Europe and to compare these changes with the climate impact of the first set (topographic change) of experiments. PLIO and LGM experimental setups and boundary conditions are identical to those of Mutz et al. (2018), with the addition of isotope prediction enabled. In our paleosimulations we account for changing pCO 2 , land surface conditions, including vegetation change and land ice, albedo, orbital variation, and sea surface temperatures, which potentially cause shifts in δ 18 O p values.

Geophysical Research Letters
Each experiment from both sets were run for 18 years, including 3 years necessary for model spin-up. In both sets of experiments, we discuss summer (June-July-August (JJA)) climate variables as numerous paleoclimate/paleoaltimetry studies (e.g., Garzione, 2008;Quade et al., 2007) or for the Alps (Campani et al., 2012) rely on pedogenic (soil) carbonate, a proxy material that forms preferentially during the warm and dry season (Breecker et al., 2009). In order to overcome the uncertainty linked to the exact season of pedogenic carbonate formation in the Alps, we provide an additional analysis of mean annual ECHAM5wiso outputs (see Figures S2-S5) and conclude that our results are completely transferable for annual mean values.

Climate Model Results
Simulated δ 18 O p values for the control simulation with PD boundary conditions are in good agreement (within 1-2‰ for most locations) with the Global Network of Isotopes in Precipitation (accessible at  Campani et al., 2012). The PD simulation also successfully reproduces observed precipitation ( Figure S3) and temperature ( Figure S4) patterns on both annual and seasonal scales. However, the model tends to underestimate present-day precipitation amount over the Alps, especially for the summer. The simulation with preindustrial boundary conditions (Alps100) shows~0.5°C lower annual mean temperatures over Europe and slightly higher rainfall than for the present day. The differences in δ 18 O p between the PD and Alps100 experiments over Europe are within 1‰ (not shown).
In the Alps100 experiment δ 18 O p values decrease with higher elevation over the Alps (down to −10‰) in JJA (Figure 1a), following the adiabatic temperature decrease with higher elevation (Figure 2a) and high rainfall amounts (Figure 3a). On the annual scale, δ 18 O p values decrease over the Alps down to −12-13‰ ( Figure S5). For two experiments with modified topography over the Alps (Alps0 and Alps150) the magnitude of difference in δ 18 O p values is up to 5‰ relative to the Alps100, with less negative values than the Alps100 for the Alps0 simulation and more negative values for the Alps150 simulation (Figures 1b and  1d). Note that the maximum change in δ 18 O p occurs over the highest alpine topography, while over the forelands the magnitude of difference between experiments with modern and modified topography is less and within 1-2‰ (Figures 1b, 1d, and 4).
Shifts in δ 18 O p values between pairs of simulations with different topography (e.g., Figures 1b-1d) are mainly controlled by temperature and precipitation effects. More specifically, spatial variations in the δ 18 O p signal follow the adiabatic temperature decrease to the highest elevation ( Figure 2) and are amplified by an increase in orographic precipitation associated with surface uplift (Figure 3). Furthermore, humidity and the evaporation-precipitation ratio (e/p), which describe the local hydrological cycle, are subjected to only minor variations within different experiments ( Figure S6), thereby differing from previous findings for the much larger Himalaya-Tibet orogen (Botsyun et al., 2016). Estimated ILR for the Alps150, Alps100, and Alps50 simulations are~2.9,~2.8, and 2.5‰/km, respectively (Figure 4), and follow presentday modeled (Rowley et al., 2001) as well as present-day empirical global (Poage & Chamberlain, 2001) and regional (Swiss Alps; Campani et al., 2012) elevation-δ 18 O relationships.
Experiments with LGM and Pliocene boundary conditions show pronounced changes in climate over the study area (Figures 2e, 2f, 3e, and 3f). Both 2-m air temperature and precipitation simulated by ECHAM5-wiso are consistent with the Paleoclimate Modeling Intercomparison Project version 3 results ( Figures S7-S9). The simulated LGM climate is characterized by lower than modern JJA temperatures, with >10°C of difference over Scandinavia and~2-8°C of difference over central Europe (Figure 2f). This follows a global decrease in temperature due to the presence of ice sheets in the Northern Hemisphere. The LGM JJA precipitation over Europe is also less than modern precipitation (Figure 3f), especially over low-elevation regions of Europe (~20-80 mm/month of difference) and the alpine region (up to 100 mm/month of difference) due to (i) a reduction of precipitable water in the atmosphere and (ii) a decrease in the convective precipitation contribution to the hydrological cycle over the Alps.
LGM δ 18 O p values over Europe are also strictly different from the present day. In general, LGM δ 18 O p values follow the pattern of a decrease in temperature. Maximum isotopic changes over Scandinavia are >8‰ (Figures 1f and S10). Over high-elevation alpine regions, simulated LGM δ 18 O p values are comparable to modern values, while low-elevation δ 18 O p values are~4‰ more negative (Figures 1 and 4). Over the Alps, the LGM ILR reduction results from the collapse of orographic precipitation due to generally cooler temperatures over the European continent. Simulated Pliocene warming over Europe is~2°C, that is, consistent with the results of the PlioMip project. This overland warming, together reduced orographic precipitation over the Alps contributes to less negative δ 18 O p values than under modern conditions. Pliocene atmospheric temperature lapse rate reduction directly influences the ILR, which is reduced to nearly 0 in the Pliocene Alps compared to modern (Figure 4).

Discussion
The results of our experiments with varied alpine elevations suggest that surface uplift of the Alps affects δ 18 O p values by about 4-5‰ within the high-Alpine region (e.g., if Alps100 and Alps0 simulations are compared; Figures 1 and 4). For regions adjacent to the Alps, where topography was kept constant, the δ 18 O p change is within 1‰ (Figure 4). Our paleoclimate simulations with realistic paleoenvironmental boundary conditions (PLIO and LGM) show that global climate change may induce a change in δ 18 O p values on the order of 5‰ over central Europe (Figures 1e, 1f, and 4). This is comparable in magnitude to the spatially focused signal from the topographic uplift experiments, but the global climate change effect produces a spatially different pattern across central Europe.
Our modeling results hence show that observed changes in δ 18 O p values due to surface uplift are in line with assumptions commonly being made for stable isotope paleoaltimetry in the region. Further uplift-related changes in δ 18 O p are expected to be of similar magnitude as changes in δ 18 O p values between glacialinterglacial time intervals, yet occur on different spatial (and most likely temporal) scales. This allowing regional patterns in δ 18 O p values to be used to discriminate between signals of surface uplift and signals of (short-term) climate change. Given that soil carbonate as a common paleoaltimetry proxy material forms over hundreds to tens of thousands of years, such proxy δ 18 O p records are prone to contain a signal of climate change. This result highlights the importance of quantifying the (potentially equally large) climatically induced signal preserved in proxy data prior to interpreting proxy records for surface uplift (e.g., Ehlers & Poulsen, 2009;Mulch, 2016).
On that account, low-elevation δ 18 O estimates, which provide a reference/low-elevation baseline for paleoelevation reconstructions, are collected in many studies to complement high-elevation samples (δ-δ approach; e.g., Campani et al., 2012;Gébelin et al., 2012;Mulch, 2016). The δ-δ approach assumes minimal surface elevation change at low elevations and provides a means for deducing the amount of surface elevation difference between foreland and high-elevation orogen interiors from time-equivalent high-elevation samples. Our modeling results show that low-elevation areas surrounding the Alps experience δ 18 O p changes (1-2‰), following mountain uplift (Figure 4; Alps0, Alps100, Alps150 experiments), whereas the same areas experience larger δ 18 O p changes (up to 4‰) in response to changes in global climate (Figure 4; LGM and PLIO experiments compared to Alps100). This simulated difference in the isotopic signal at low-elevation and high-elevation locations shows that the basic assumption that low-elevation δ 18 O p records at the foothills of the orogen are largely unaffected by downwind surface elevation changes is consistent with model results. Our results further show that at least half of the low-elevation signal can be exclusively attributed to global climatic change and not to topographic 10.1029/2019GL086046

Geophysical Research Letters
change. Since low-and high-elevation δ 18 O records are both affected by the impact of climate change on δ 18 O p , tracking differences between δ 18 O p at low and high elevations over time (Δ(δ 18 O p )) is very likely to strongly reduce climate change bias on stable isotope paleoaltimetry reconstructions (e.g., Mulch, 2016).
Interestingly, JJA δ 18 O p in the Alps150/Alps100 and LGM experiments attains similar values at the highest elevations, but shows significant differences west of the alpine crest, in the upwind direction of major winds across central Europe (westerlies; Figure 4a). Whereas in the LGM simulation, the δ 18 O p values gradually decrease (starting at 4°W, at or very close to the North Atlantic moisture source) toward the highest elevations (at~10°E), the δ 18 O p values in the Alps150/Alps100 simulations rapidly decrease close to the mountain range (driven by orographic precipitation changes). This result is not overly surprising as in both cases moist air is exposed to cooler atmospheric temperatures either during continental transport and/or increasing elevation. The interpretation of proxy data may overcome the consequences of short-term global climate changes if the analyzed proxy material for paleoaltimetry was sampled at both low and high elevations and the sample material integrates over long enough time (e.g., kiloyears) to minimize short-term climate fluctuations. However, in the case when short-term proxies (e.g., clam shells, stalagmites, varied sediments, fossil teeth) transient warm-cold phases interpretation of climate versus surface uplift signals may be compromised.
Our starting hypothesis was that if topography of the Alps changed over the Cenozoic then a δ 18 O p signal of this elevation change could be recovered through paleoaltimetry studies, independent of global climate change effects. In evaluating this hypothesis, we find that oxygen isotope values from multiple lowelevation sites with different distances (up to 4°upwind; Figure 4) from the mountain range, in combination with high-elevation sample sites, may enable differentiation between a primarily topographic and paleoclimate induced δ 18 O p signal. However, the spatial variability of δ 18 O p values (1-2‰) in lowelevation regions adjacent to mountains (Figure 1) needs to be addressed by using (i) multiple lowelevation geological sections in different geographic areas, including those remote from the mountain belt regions, and (ii) multiple proxy approaches (e.g., Δ 47 paleothermometry, paleobotany-related methods). For example, if measurements of soil carbonate δ 18 O values are paired with carbonate clumped isotope temperatures (T(Δ 47 )) to more reliably reconstruct low-elevation δ 18 O p values, then uncertainties in changing temperature or precipitation in a region can be minimized. Combining these results with sedimentological, paleofloral, and other isotopic studies (e.g., on fossil teeth) with high-elevation isotopic records (such as hydrous silicates) to reconstruct Δδ 18 O p allows improved estimates of paleoelevation changes (e.g., Mulch, 2016).
Additionally, our model results support previous work that documents decreases in ILR during warm climate conditions that mask the real elevation of mountain ranges (Poulsen & Jeffery, 2011) and hence lead to an underestimation of the true surface elevation. Changes in the ILR during topographic development have been hypothesized in several recent studies (Botsyun et al., 2019;Insel et al., 2012;Poulsen et al., 2010;Poulsen & Jeffery, 2011). These studies show that using the present-day δ 18 O p -elevation relationship (e.g., empirical global mean isotopic lapse rate as defined in Poage and Chamberlain (2001) or theoretical local isotopic lapse rate as in Rowley et al. (2001)) could imply a significant error in paleoelevation estimates during warmer past climates. Here we observe that surface uplift of the Alps has a minor impact on the regional ILR for this region, while global climate change has an imprint up to 2.5‰/km (e.g., for the Pliocene case). As the Alps have a relatively small size, a warming of~2°C (PLIO run) makes this mountain range "invisible" for a "standard" paleoaltimetry approach that uses the present-day ILP (Figure 4). This strong climate impact on Pliocene ILRs helps to evaluate similar effects under warmer climate conditions and also suggests that paleoclimate simulations should be applied on a case-by-case basis for reconstructing any temporal variations in the ILR.
This study evaluated the potential of δ 18 O p proxy data to recover surface elevation changes relative to global climate change. We have done this through an example using the European Alps. However, the results presented may have broader implications for similar "small" orogens elsewhere. The European Alps are (for example) similar in size and relief to the Caucasus, Atlas, Southern Alps (New Zealand), and Washington Cascades. The results presented here document that the Alps produce a 3-4‰ lower isotopic signal of surface uplift compared to larger orogens such as the Andes and the Himalaya ( et al., 2010). The smaller δ 18 O p signal predicted for the Alps results from their small size and decreased ability to block atmospheric winds, as well as the lower magnitude of elevation change compared to larger orogens. Other small orogens may experience a diminished isotopic signal of surface uplift which would require a similar modeling analysis as we present here conducted on a case-by-case basis to identify what multiproxy methods and far-field baseline sample locations are maximize the impact of the δ-δ approach. Future work in the European Alps to reconstruct the timing and magnitude of surface uplift will require time-specific climate simulations with appropriate paleogeography, as well as additional far-field low-elevation samples along upwind trajectories for application of the δ-δ approach.