Miocene high elevation in the Central Alps

Reconstructing Oligocene–Miocene paleoelevation contributes to our understanding of the evolutionary history of the European Alps and sheds light on geodynamic and Earth surface processes involved in the development of Alpine topography. Despite being one of the most intensively explored mountain ranges worldwide, constraints on the elevation history of the European Alps remain scarce. Here we present stable and clumped isotope measurements to provide a new paleoelevation estimate for the mid-Miocene (∼14.5 Ma) European Central Alps. We apply stable isotope δ–δ paleoaltimetry to near-sea-level pedogenic carbonate oxygen isotope (δ¹⁸O) records from the Northern Alpine Foreland Basin (Swiss Molasse Basin) and high-Alpine phyllosilicate hydrogen isotope (δD) records from the Simplon Fault Zone (Swiss Alps). We further explore Miocene paleoclimate and paleoenvironmental conditions in the Swiss Molasse Basin through carbonate stable (δ¹⁸O, δ¹³C) and clumped (Δ₄₇) isotope data from three foreland basin sections in different alluvial megafan settings (proximal, mid-fan, and distal). Combined pedogenic carbonate δ¹⁸O values and Δ₄₇ temperatures (30±5 ^∘C) yield a near-sea-level precipitation δ¹⁸O_w value of $-\mathrm{5.8}\pm \mathrm{1.2}$ ‰ and, in conjunction with the high-Alpine phyllosilicate δD value of $-\mathrm{14.6}\pm \mathrm{0.3}$ ‰, suggest that the region surrounding the Simplon Fault Zone attained surface elevations of >4000 m no later than the mid-Miocene. Our near-sea-level δ¹⁸O_w estimate is supported by paleoclimate (iGCM ECHAM5-wiso) modeled δ¹⁸O values, which vary between −4.2 ‰ and −7.6 ‰ for the Northern Alpine Foreland Basin.

1 Introduction

Past elevations of mountain ranges provide insight into the coupled climatic and geodynamic processes that shape orogenic belts. The European Alps are one of the most intensively investigated mountain belts worldwide, and yet there are surprisingly few studies addressing their surface uplift history. The topographic evolution of continent–continent collision zones such as the European Alps is mainly controlled by isostatic compensation for crustal and/or lithospheric deformation caused by plate convergence (e.g., Beaumont et al., 1996; Schmid et al., 1996; Willett et al., 1993). Various mechanisms may have contributed to the increase in post-collisional Alpine surface elevation. These include horizontal shortening, thickening of the continental crust and following isostatic adjustment of the lithosphere (e.g., Pfiffner et al., 2002), mantle-scale uplift caused by slab break-off (e.g., Lippitsch, 2003), slab rollback and associated delamination within the lithosphere (Kissling and Schlunegger, 2018; Schlunegger and Kissling, 2015), and asthenospheric upwelling (e.g., Lyon-Caen and Molnar, 1989). At the same time, climate change will impact topography, e.g., through accelerated erosional unroofing forcing isostatic rebound of the unloaded lithosphere (Cederbom et al., 2004; Champagnac et al., 2009; Mey et al., 2016). Despite being a prime example of continent–continent collision following protracted convergence and subduction of oceanic lithosphere (e.g., Frisch, 1979; Froitzheim et al., 2008; Schmid et al., 1996; Stampfli et al., 1998), recent studies in the European Central Alps suggest that compressional tectonics may not be the sole driver for surface uplift and thus for the formation of Alpine topography (e.g., Kissling and Schlunegger, 2018; Schlunegger and Kissling, 2015). Convergence between the European and the Adriatic plates commenced in the late Cretaceous and led to the subsequent collision between the two continental plates at ca. 35 Ma (Schlunegger and Kissling, 2015). Post-collisional slab break-off (at ca. 32 Ma), and continuing slab rollback of the subducting lithosphere beneath the Alpine arc associated with lithospheric delamination, may have contributed to the rise of topography (Kissling and Schlunegger, 2018).

Previous studies of the Central Alps suggested a wide range of Oligocene to Miocene paleoelevation estimates. Schlunegger and Kissling (2015) inferred possible mean paleoelevations of 1900±1000 m for the post-20 Ma Alpine drainage divide from spacing and cross-sectional widths of alluvial megafan conglomerates deposited in the Northern Alpine Foreland Basin. Stable-isotope-based paleoelevation studies suggested paleoaltitudes of 2300±650 m obtained from oxygen isotope and ${}^{\mathrm{87}}\mathrm{Sr}{/}^{\mathrm{86}}\mathrm{Sr}$ compositions of shark teeth from early Miocene foreland basin deposits (20.3–17.6 Ma; Kocsis et al., 2007) and of $\mathrm{2850}+\mathrm{800}/-\mathrm{600}$ m as inferred from oxygen isotopic composition of precipitation for the middle Miocene (Campani et al., 2012). Summit (or maximum) elevations have been estimated at 2500–3000 m for large parts of the Western and Central Alps for the middle Miocene based on combined evidence from sediment budget curves, thermochronology, and sediment facies (Kuhlemann, 2007). The highest peak elevations have been provided by Sharp et al. (2005), who suggested at least 5000 m based on stable isotope compositions of minerals and fluid inclusions from fissures in the Swiss Alps, and by Jäger and Hantke (1983, 1984), who inferred >5000 m for the Swiss Bergell mountains at the Oligocene–Miocene boundary based on large erratic boulders found at great distances from their place of origin.

In this study, we complement previous work by applying stable isotope (δ–δ) paleoaltimetry (Mulch, 2016) using authigenic soil carbonates from the near-sea-level Northern Alpine Foreland Basin (NAFB) and contrast these with hydrogen isotope data from syntectonic high-Alpine fault zone silicates from the Simplon Fault Zone (SFZ). Our application of Miocene stable isotope paleoaltimetry data builds upon (and revises) former stable isotope paleoaltimetry estimates (Campani et al., 2012), as the combined analysis of well-dated mid-Miocene sediments (e.g., Kälin and Kempf, 2009; Kempf and Matter, 1999; Schlunegger et al., 1996) along and across two Miocene Alpine foreland megafans (Napf, Hörnli) with clumped isotope (Δ₄₇) thermometry data permits refining the soil carbonate-based near-sea-level δ¹⁸O estimate of meteoric water (δ¹⁸O_w).

δ¹⁸O_w values used in paleoelevation reconstructions can be affected by paleoclimate change (Botsyun et al., 2016, 2019; Botsyun and Ehlers, 2021; Ehlers and Poulsen, 2009; Insel et al., 2012; Mulch, 2016; Poulsen et al., 2010). For instance, the sensitivity of δ¹⁸O_w to regional, global, and topographic variations in paleotemperature, environmental conditions of an air mass prior to orographic ascent, evapotranspiration, water vapor recycling, and changes in vapor source has been shown to introduce uncertainties in stable-isotope-based elevation reconstructions (e.g., Mulch, 2016; Botsyun et al., 2020, Botsyun and Ehlers 2021). In particular, isotopic changes over continental Europe could be related to a variety of factors such as declining pCO₂ levels (Pagani et al., 1999), variable ocean circulation and sea surface temperatures (Flower and Kennett, 1994; Wright et al., 1992), sea-level fluctuations (Foster and Rohling, 2013), paleogeographic changes (Herold et al., 2008; Poblete et al., 2021), and other processes affecting δ¹⁸O_w (Botsyun et al., 2019; Poulsen et al., 2007; Risi et al., 2008; Roe et al., 2016; Sewall and Fricke, 2013; Sturm et al., 2010). We thus compare our newly refined near-sea-level δ¹⁸O estimate with paleoclimate simulations from the isotope-enabled ECHAM5-wiso atmospheric general circulation model (iGCM), which predicts changes in δ¹⁸O of precipitation.

Figure 1 (a) Topographic map of the Central Alps and their northern foreland basin (Swiss Molasse Basin; yellow hatched area) with the sections Fontannen (F), Aabach (A), Jona (J), and the alpine Simplon Fault Zone (SFZ) sampling locality. Dashed blue lines mark the present-day drainage divide (after Schlunegger et al., 2007). (b) Setting of the foreland basin sections within alluvial megafan deltas during the middle Miocene (15–14 Ma; drawn after Berger et al., 2005). The Fontannen section is located in the proximal part of the Napf megafan, while the Aabach and Jona sections represent a distal and a mid-fan setting, respectively, in the Hörnli megafan. The dashed blue line indicates the Oligocene–Miocene drainage divide (after Schlunegger et al., 2007).

2 Study area

2.1 The Alps and the Swiss Molasse Basin

The European Alps formed as a result of the northward drift of the Adriatic microplate and the associated formation of a southward-directed subduction zone beneath the Tethys Ocean. The late Cretaceous to Paleogene closure of the Alpine Tethys led to the collision between the Adria and Europe continental plates (Handy et al., 2010; Schmid et al., 1996; Stampfli et al., 1998). Subsequent post-collisional convergence resulted in overthrusting and stacking of nappe sheets (e.g., Schmid et al., 1996). The Northern Alpine Foreland Basin formed due to elastic down-warping of the European lithosphere resulting from subduction slab load and topographic load, and it accommodated eroded material from the northward-propagating Alpine thrust front (Fig. 1; e.g., Matter et al., 1980; Schlunegger and Kissling, 2015). The Oligo-Miocene Swiss Molasse Basin (SMB) represents the central part of the Northern Alpine Foreland Basin (Fig. 1). Deposition of several kilometer-thick sequences of basin fill in the Swiss Molasse Basin started in the early Oligocene (e.g., Pfiffner, 1986; Kempf et al., 1997, 1999) and continued until the late Miocene–early Pliocene when basin inversion led to erosion of molasse sediments (Cederbom et al., 2004, 2011; Mazurek et al., 2006). During this period deposition changed twice from shallow marine to terrestrial, resulting in two regressive shallowing, coarsening, and thickening-upward megacycles. The Swiss Molasse Basin is therefore divided into four lithostratigraphic units: the Lower Marine Molasse (“Untere Meeresmolasse”, UMM), the Lower Freshwater Molasse (“Untere Süßwassermolasse”, USM), the Upper Marine Molasse (“Obere Meeresmolasse”, OMM), and the Upper Freshwater Molasse (“Obere Süßwassermolasse”, OSM) (e.g., Matter et al., 1980; Schlunegger et al., 1996).

Advancing surface uplift of the Alpine mountain belt led to formation of extended drainage networks and alluvial megafans since ca. 32–30 Ma (Kempf et al., 1999; Kuhlemann and Kempf, 2002; Schlunegger et al., 1997) that formed large dispersal systems with cross-sectional widths of nearly 30 km and corresponding stream lengths of >150 km (Schlunegger and Kissling, 2015). A reduction in sediment flux (Kuhlemann et al., 2001) paired with ongoing basin subsidence resulted in a shift from basin overfill (Lower Freshwater Molasse) to underfill (Upper Marine Molasse) at 20 Ma (Garefalakis and Schlunegger, 2019). A subsequent increase in erosional flux (Kuhlemann et al., 2001) together with a lowering of the eustatic sea level (Garefalakis and Schlunegger, 2019) enabled propagation of the fan deltas towards the basin center and led to a shift to an overfilled basin again at ca. 17 Ma and therefore to the establishment of terrestrial OSM sedimentation by that time (e.g., Kuhlemann and Kempf, 2002).

2.2 Alluvial megafans of the OSM

The Napf and Hörnli megafans (Fig. 1b), among other fan systems, formed at the Alpine front and merged into a basin axial drainage system in the central part of the Swiss Molasse Basin. Deposition of the Napf and Hörnli megafans initiated during the OMM (Kuhlemann and Kempf, 2002) or possibly earlier (Garefalakis and Schlunegger, 2019; Schlunegger and Kissling, 2015). Persistent progradation and accumulation of debris formed subaerial deltas, and the relief of the SMB megafans may have reached several hundred meters above base level depending on distance to the apex and slope geometry (Garefalakis and Schlunegger, 2018). During the youngest depositional phase (OSM; ca. 17–11 Ma) sediments of SMB megafans were predominantly composed of amalgamated conglomerate and sandstone packages with mudstone interlayers at the basin margin and alternations of finer-grained sandstone beds, mudstones, and marls in the basin center. Sediments of the Hörnli megafan document a short-lived marine ingression at ca. 18.5–18.0 Ma (Bolliger et al., 1995; Keller, 1989), whereas deposits of the center of the Napf megafan show no evidence of marine sedimentation (Schlunegger et al., 1996). However, marine sedimentation between the megafan deltas continued for at least 1 Myr after the transition from marine OMM to terrestrial OSM within fan deltas (Schlunegger et al., 1996), resulting in lateral facies and elevation changes. As such, the sedimentary sections presented here are in close proximity to the retreating Molasse Sea; yet, depending on their fan position, they may have developed up to several hundred meters above mid-Miocene sea level. Here we present oxygen (δ¹⁸O), carbon (δ¹³C), and clumped (Δ₄₇) isotope data from three fully terrestrial SMB sections (Fontannen, Jona, and Aabach) to explore the change in environmental conditions in proximal to distal depositional settings in the Miocene Napf and Hörnli megafans (Fig. 1b), and we relate these data to published stable oxygen isotope data in the evolving Alps (Campani et al., 2012) to update the paleoelevation estimate by Campani et al. (2012) and previous authors for that time.

2.3 Sampled sections

Our study of Miocene paleosol carbonate takes advantage of numerous pedogenic soil horizons embedded in a detailed geochronologic framework determined from previous studies of paleomagnetostratigraphy (Kälin and Kempf, 2009; Kempf et al., 1997; Kempf and Matter, 1999; Schlunegger et al., 1996), radiometric age data on volcanic ash layers (Gubler et al., 1992; Schmieder et al., 2018), and mammal biostratigraphy (e.g., Bolliger, 1992; Kälin, 1997). The local paleomagnetic sections were correlated with the astronomically tuned Neogene timescale (ATNTS2012) of Hilgen et al. (2012).

The almost 1000 m thick Fontannen section is situated in the proximal part of the Napf alluvial megafan and covers an age range of 17.6–13.3 Ma (Fig. 1b). The section is composed of alternating massive conglomerates, sandstones, and silty mudstones of the OSM (Fig. 2). The base of the section (0 to ∼100 m) is composed of the Schüpferegg Conglomerate, which represents the terrestrial equivalent of the marine Luzern and St. Gallen formations of the OMM (Garefalakis and Schlunegger, 2019; Keller, 1989). The Schüpferegg Conglomerate is overlain by the Napf beds, which mainly comprise conglomerates and silty mudstones (Schlunegger et al., 1996). Massive conglomerates of up to 100 m thickness dominate the stratigraphy in this region and reveal the proximal position of the section within the alluvial fan (Figs. 2 and 3a). Well-developed paleosols with carbonate nodules occur in mud–siltstone interlayered between the conglomerate beds. Pedogenic horizons are up to 50 cm thick and occasionally show mottling in grey, purple, and yellow. Age constraints of the Fontannen section are given by deposits associated with the Ries meteorite impact (14.81±0.02 Ma; Schmieder et al., 2018) and by four mammal faunal sites in three zones (MN 4, MN 5, and MN 6; Kälin, 1997). We inferred pedogenic carbonate ages based on the local magnetostratigraphy (Schlunegger et al., 1996; Kempf et al., 1997; Figs. 2 and S1 in the Supplement). Carbonate δ¹⁸O and δ¹³C data for the Fontannen section have been previously published in Campani et al. (2012), and Δ₄₇ temperatures have been provided by Methner et al. (2020; Table S6 in the Supplement).

Figure 2 Magnetostratigraphy and pedogenic carbonate δ¹⁸O and δ¹³C values for Swiss Molasse Basin sections. In total, 141 individual carbonate nodules were analyzed for Fontannen (Campani et al., 2012) and Jona, as well as 101 for Aabach. Only the mean value for each carbonate bearing horizon is shown here (Fontannen n=50; Jona n=56; Aabach n=34). δ¹⁸O_c and δ¹³C values are plotted next to their stratigraphic position. Assignment of polarity chrons to the ATNTS2012 is according to Kempf et al. (1997), Kempf and Matter (1999), and Schlunegger et al. (1996). Also shown is the stratigraphic position of important marker beds (Küsnacht bentonite, Meilen Limestone and its equivalent for proximal regions, the Hüllistein Conglomerate) with dated ages. Stars show locations of mammal sites (⋆¹ projected Oeschgraben (MN 6), ⋆² Chatzloch (Frohberg, MN 5), ⋆³ Vermes 1 (MN 5), ⋆⁴ Tägermaustrasse (MN 4), ⋆⁵ projected Hirschengrabentunnel (uppermost MN 5), ⋆⁶ Tobelholz (middle–upper MN 5), ⋆⁷ Aabach 2 (middle MN 5), ⋆⁸ Aabach 3 (middle MN 5), ⋆⁹ Käpfnach-Rotwegstollen (middle MN5), ⋆¹⁰ Bachtel-Ornberg (uppermost MN 6), ⋆¹¹ Frohberg (MN 5), ⋆¹² Tobel Hombrechtikon (upper MN 5), ⋆¹³ Güntisberg (middle MN 5), ⋆¹⁴ Martinsbrünneli (middle MN 5), ⋆¹⁵ Hubertingen – MN 4b). Mammal sites are after Kälin (1997) and Bolliger (1992). Stratigraphic names are after Büchi (1959), Kempf et al. (1997), and Matter (1964). Colored diamonds represent obtained (Δ₄₇) temperatures for each section. Fontannen (Δ₄₇) temperatures are from Methner et al. (2020). Each stratigraphic section covers an age range of ca. 17.5 to 14 Ma (see the Supplement).

The 750 m thick Jona section covers an age range of 16.8–13.7 Ma and is located in the mid-fan to proximal part of the Hörnli alluvial megafan (Fig. 1b). It is mainly composed of alternating conglomerates and mudstones (Kälin and Kempf, 2009). Similar to the Fontannen section, the stratigraphy of the Jona section is dominated by frequent conglomerate horizons. Almost all outcrops reveal conglomerates up to several tens of meters thick overlaying grey sandstones and marls, the latter characterized by pedogenic overprint (Figs. 2 and 3d). Pedogenic features of the Jona section include carbonate nodules, calcified roots, and occasionally strong mottling with grey, yellow, and purple. In contrast to the proximal Fontannen section, paleosol horizons in the Jona section formed more frequently and are thicker. The up to 2 m thick paleosols are well developed and contain abundant carbonate nodules. The section is dated through a projection of the Küsnacht bentonite (14.91±0.09 Ma) to a level situated at ∼710 m of the Jona section and through seven mammal sites comprising faunal zones MN 4b, MN 5, and MN 6 (Bolliger, 1992). Furthermore, the Hüllistein conglomerate at ∼280 m (16 Ma; Kempf et al., 1997) was deposited during a fan-wide mass flow event and serves as an important regional marker horizon. It allows a lateral correlation with the distal counterpart Meilen limestone in the Aabach section 20 km further to the west (e.g., Bürgisser, 1984; Kempf and Matter, 1999). Our carbonate nodule age constraints are based on the magnetostratigraphic study of Kempf et al. (1997).

Figure 3 Outcrop conditions of the Swiss Molasse Basin sections and close-ups of sampled pedogenic carbonate nodules. (a–c) Proximal Fontannen section: ca. 1 m thick red paleosol profile overlain by massive sandstones and conglomerates (a). Mottled paleosol with carbonate nodules (b) and pedogenic carbonate nodules (c). (d–f) Mid-fan Jona section: ca. 2 m thick paleosol profile with intense mottling intersected by river channel sandstone (d). Ca. 20 cm long calcified root within a paleosol profile (e) and pedogenic carbonate nodules (f). (g–i) Distal Aabach section: floodplain deposits with typical strong pedogenic overprint and intense mottling in yellow and red (g). Paleosol abundant in carbonate nodules (h) and individual pedogenic carbonate nodules (i).

The 352 m thick Aabach section covers an age range of 17.3–14.8 Ma (Kempf and Matter, 1999). These sediments were deposited in the distal part of the Hörnli alluvial megafan (Fig. 1b; Kälin and Kempf, 2009) and are characterized by thick mudstones reflecting typical overbank deposits on floodplains. The lithostratigraphic units consist of alternating sequences of mudstones, marls, siltstones, and fine-grained sandstones (Fig. 2). The typical marl in this section is greyish or yellowish and shows strong evidence of pedogenic overprinting, including frequent mottling in intense purple and yellow, root traces and hackly structures resulting from bioturbation, and/or shrink–swell features (slick-and-slide structures) resulting from seasonal wetting and drying (Fig. 3g, h). Pedogenic horizons often contain abundant carbonate nodules or show caliche formation. The Aabach section reveals the highest frequency of paleosol occurrence of all three SMB sections, with pedogenic overprint in almost all mudstones and marls. Furthermore, the characteristic lack of conglomerates and the presence of lacustrine marls and the Meilen limestone at ∼170 m indicate that this section was deposited on a floodplain of the distal part of the alluvial megafan. We inferred carbonate nodule ages based on the local magnetostratigraphy of Kempf and Matter (1999). Radiometric ages of the projected Küsnacht (14.91±0.09 Ma) and Urdorf bentonites (15.27±0.12 Ma), as well as five mammal sites (all in mammal faunal zone MN 5; Bolliger, 1992), allowed for correlation with the Neogene polarity timescale (Hilgen et al., 2012).

3 Proxies and methods

3.1 Sampling strategy

Miocene Swiss Molasse Basin paleosols show typical characteristics of pedogenic overprinting including mottling, root traces, and different stages of soil carbonate development (Fig. 3). For all analyzed sections we targeted well-developed undisturbed paleosols with micritic pedogenic carbonate nodules typically 1–5 cm in diameter (Fig. 3). Carbonate-bearing horizons occurred mainly in sequences of stacked paleosols or were overlain by sandstones or conglomerates separated by a sharp erosional contact, indicating that the top of the soil profile was truncated. In total, we collected 383 pedogenic carbonate nodules from 140 layers of former B-horizons in different alluvial megafan settings (proximal, mid-fan, and distal) for analysis of δ¹⁸O, δ¹³C, and Δ₄₇. Samples were collected along small rivers and from well-exposed sections of the SMB alluvial deposits. We compiled coherent age models for all three sections with a resolution of ≤100 kyr and an error of ±80–150 kyr for individual paleosol horizons based on magnetostratigraphy studies (Kälin and Kempf, 2009; Kempf et al., 1997; Kempf and Matter, 1999; Schlunegger et al., 1996). Errors in the age models are based on uncertainties of sample placement within the stratigraphic section and were calculated for the period of the lowest sedimentation rate, thus representing the maximum error. Details on the determined age model and error calculation are given in the Supplement.

3.2 Carbonate stable isotope (δ¹⁸O, δ¹³C) analyses and Δ₄₇ paleothermometry

δ¹⁸O and δ¹³C data were obtained from pedogenic carbonate nodules in well-developed paleosol horizons. Pedogenic carbonate nodules typically form in soils of arid to sub-humid zones due to chemical precipitation from supersaturated soil water (Cerling, 1984; Cerling and Quade, 1993). Their δ¹⁸O and δ¹³C values are controlled by δ¹⁸O of meteoric (soil) water and soil CO₂, respectively, and are sensitive to changes in temperature, water availability, soil respiration rates, and the proportions of C3:C₄ biomass of local vegetation (Cerling, 1984). Carbonate clumped isotope (Δ₄₇) paleothermometry is based on the measurement of the abundance of the rare isotopes ¹³C and ¹⁸O in the carbonate mineral lattice (Eiler, 2007; Ghosh et al., 2006). The effect of “clumping” of rare isotopes is temperature-dependent and unrelated to the δ¹⁸O value of the water from which the carbonate formed. Pedogenic carbonates have been shown to reliably record primary Δ₄₇ values in the SMB sections (Methner et al., 2020). We therefore determined Δ₄₇ paleotemperatures for each of the three investigated sections (n=5 samples for Fontannen, Jona, and Aabach) to calculate δ¹⁸O_w values of meteoric waters. Both clumped isotope (Δ₄₇) and δ¹⁸O and δ¹³C analyses were performed at the Joint Goethe University–Senckenberg BiK-F Stable Isotope Facility (Frankfurt, Germany). Pedogenic carbonate oxygen and carbon isotope values are reported using standard δ-per mill notation relative to VSMOW and VPDB, respectively. Analytical details can be found in the Supplement.

3.3 Stable isotope paleoaltimetry

Stable isotope paleoaltimetry relies on the systematic decrease in δ¹⁸O and δD values of meteoric water with increasing elevation on the windward side of an orographic barrier (e.g., Currie et al., 2005; Poage and Chamberlain, 2001; Rowley and Garzione, 2007; Siegenthaler and Oeschger, 1980). Ascending air masses undergo adiabatic cooling and rainout with increasing altitude, which leads to fractionation-driven depletion of ¹⁸O in the residual water vapor. The rainfall becomes increasingly depleted in ¹⁸O with progressive rainout. Even though altitudinal lapse rates of δ¹⁸O and δD in precipitation may not necessarily be constant through time and space (Botsyun et al., 2020; Ehlers and Poulsen, 2009; Poulsen et al., 2010), averaged global oxygen isotope lapse rates show a systematic decrease in δ¹⁸O of meteoric water with increasing elevation (e.g., Poage and Chamberlain, 2001). The present-day Alpine δ¹⁸O lapse rate averages −0.20 ‰ per 100 m (Campani et al., 2012). The comparison of the δ¹⁸O and/or δD values of meteoric waters from age-equivalent low-elevation and high-elevation sites (δ–δ) strengthens paleoelevation reconstructions by reducing the impact of long-term climate change on the isotope proxy records (Mulch, 2016; Mulch and Chamberlain, 2018). In this study we evaluate different low-elevation δ¹⁸O records from foreland basin paleosols and compare them with an age-equivalent high-elevation δD record from the Simplon Fault Zone published in Campani et al. (2012). We then relate the difference in δ¹⁸O of meteoric water between these sites to the difference in paleoelevation and update the previous interpretation presented by Campani et al. (2012) (Fig. 4). Normal fault and detachment systems that served as pathways for meteoric water, which percolated along such fault zones (e.g., Mulch and Chamberlain, 2007), represent valuable archives for high-elevation rainfall δ¹⁸O (or δD) values. Hydrous silicate minerals that formed synkinematically during deformation within these fault systems undergo isotopic exchange with meteoric-derived fluids. The δD values of such synkinematic hydrous minerals hence record the δD of the infiltrating meteoric waters and can therefore be used to reconstruct δD (and δ¹⁸O) of high-elevation precipitation. High-elevation samples used here include hydrous phyllosilicates (muscovite, biotite, and chlorite) from the Zwischbergen segment of the Simplon Fault Zone, a ca. 14.5 Ma major extensional detachment in the Central Alps (Fig. 1) that has been shown to be a conduit for meteoric water (e.g., Campani et al., 2012). The formation of the Simplon Fault Zone is associated with Miocene orogen-parallel extension and coeval with orogen-perpendicular shortening (Mancel and Merle, 1987; Mancktelow, 1992). Crustal extension and development of this major extensional fault system promoted large-scale exhumation of the Lepontine metamorphic dome, which underwent rapid exhumation from 20 Ma. Slip movement along the Simplon Fault Zone persisted from 30 to 3 Ma and peaked between 18 and 15 Ma accompanied by rapid exhumation rates, accelerated footwall cooling with the highest rates around 20 Ma (Campani et al., 2010; Grasemann and Mancktelow, 1993; Mancktelow, 1992), and infiltration of meteoric fluids with δD and δ¹⁸O values of $-\mathrm{107}\pm \mathrm{2}$ ‰ and $-\mathrm{14.6}\pm \mathrm{0.3}$ ‰, respectively (Campani et al., 2012).

Figure 4 Simplified sketch of the mid-Miocene Swiss Molasse Basin (SMB, 15–14 Ma) showing the distribution of alluvial megafans along the northern Alpine flank. The sampled sections are located at different paleoaltitudes within the megafan deltas, which results in an internal elevation difference of ca. 300 m (±100 m; see text for discussion). We use the difference in precipitation δ¹⁸O values (Δδ) between the low-elevation SMB records (F, J, A; see Fig. 1) and the high-elevation (projected) Simplon Fault Zone (SFZ) to calculate the elevation difference (Δz) between these sites. The grey arrow depicts the paleo-discharge direction of the drainage system. Distribution of alluvial fans is adopted from Berger et al. (2005).

Pedogenic carbonates from the SMB sections serve as proxies for near-sea-level precipitation δ¹⁸O values. The δ¹⁸O values of pedogenic carbonate (δ¹⁸O_c) are determined by the δ¹⁸O value of soil water from which the carbonate formed, which in turn is closely linked to local meteoric water (Cerling and Quade, 1993). In combination with clumped isotope (Δ₄₇) carbonate formation temperatures we calculate δ¹⁸O_w values of meteoric water based on the δ¹⁸O_c values assuming oxygen isotope equilibrium fractionation (Kim and O'Neil, 1997; updated by Kim et al., 2007). Thus, by linking the high-elevation Simplon Fault Zone δ¹⁸O estimate with our temperature-corrected near-sea-level δ¹⁸O_w value we arrive at an updated estimation of the mid-Miocene elevation for the region surrounding the Simplon Fault Zone and the headwaters of the Napf megafan.

Soils can be affected by evaporative ¹⁸O enrichment in soil water caused by preferential loss of ¹⁶O, which can result in a bias of the soil water δ¹⁸O values (Cerling, 1984; Quade et al., 2007a). Therefore, when calculating low-elevation (near-sea-level) precipitation δ¹⁸O_w values from the δ¹⁸O carbonate record we rely on the first quartile (lowest 25 %) mean δ¹⁸O_c values to avoid any potential synkinematic bias to increased near-sea-level δ¹⁸O_w values and a resulting overestimation in our paleoelevation reconstructions. We calculate δ¹⁸O_w values for meteoric water for each SMB record for the 15.5 to 14.0 Ma time interval (instead of the entire records) to provide an age-equivalent estimate when compared to the high-elevation Simplon Fault Zone data.

3.4 Paleoclimate simulations

For regions of variable topography the application of isotope-tracking climate models allows us to estimate the impact of global paleoclimatic changes and the influence of surface uplift on regional δ¹⁸O_w values (e.g., Botsyun et al., 2019; Ehlers and Poulsen, 2009; Feng et al., 2013; Poulsen et al., 2010). Here we discuss the results from ECHAM5-wiso general circulation model (GCM) experiments, first performed in Botsyun et al. (2020), in the context of our paleoelevation reconstruction of the Central Alps. The ECHAM5-wiso is an atmospheric GCM complemented with a water isotope module that enables global simulation of δ¹⁸O patterns in precipitation (Werner et al., 2011). ECHAM5-wiso is able to reproduce both climate- and altitude-related variations of δ¹⁸O_w over the European region (Botsyun et al., 2020; Langebroek et al., 2011). Simulations were performed with a high-resolution grid (T159, corresponding to ∼80 km in latitude and longitude at the Equator) with 31 vertical levels (up to 10 hPa). Different topographic scenarios (topography experiments) were evaluated in order to quantify the signal of surface uplift preserved in δ¹⁸O_w over the Northern Alpine Foreland Basin and the Alpine orogen. In the absence of model runs with Miocene boundary conditions we rely on a pre-industrial model setup with modern land–sea distribution and ρCO₂ of 280 ppm (details on the experimental setup and boundary conditions such as orbital parameters, greenhouse gas concentrations, sea surface temperature, sea-ice cover, and vegetation cover can be found in Mutz et al., 2018, and Botsyun et al., 2020). Four experiments were conducted with Alpine topography ranging between 150 % and 0 % (topography set to 250 m) of the modern mean elevation in 50 % increments (Alps150, Alps100, Alps50, and NoAlps).

We evaluate δ¹⁸O_w in three regions: (1) Northern Alpine Foreland Basin (low-elevation region), (2) Central Alps (high-elevation region), and (3) nearshore (upwind) area of modern France (distant region), which lies on the preferential moisture trajectory to the northern Central Alps and the Northern Alpine Foreland Basin (Fig. 5a). The last was chosen following Botsyun et al. (2020), who suggested that examination of low-elevation regions that lie far enough from the high-elevation data sites enhances assessment of paleoclimate changes. Results are presented for summer months (June–July–August, JJA) because pedogenic carbonate preferentially forms during the warm season (e.g., Breecker et al., 2009), although in some settings pedogenic carbonate formation has also been reported for cool seasons (Gallagher and Sheldon, 2016; Peters et al., 2013).

Figure 5 (a) Spatial distribution of ECHAM5-wiso simulated mean JJA δ¹⁸O_w for the pre-industrial (Alps100) experiment for three regions: (1) within the blue rectangular area for elevations over 1600 m (high elev.), (2) within the green rectangular area for elevations between 200 and 600 m (low elev.), and (3) within the red rectangular area for distant regions (distant region). (b) Ranges of JJA δ¹⁸O_w for the Alps150, Alps100, and NoAlps experiments for the selected regions (Botsyun et al., 2020).

4 Results

4.1 Swiss Molasse Basin near-sea-level paleosol carbonate δ¹⁸O, δ¹³C, and Δ₄₇ data

Carbonate δ¹⁸O data in Fig. 2 are shown as a mean value per horizon (one to eight measured pedogenic nodules per carbonate-bearing layer; Figs. S2 and S3 in the Supplement for δ¹³C data). Although carbonate δ¹⁸O and δ¹³C values reveal no systematic trend throughout the investigated time interval (∼17.5 to ∼13.5 Ma), each section shows high variability. Over the entire sections δ¹⁸O_c values vary between 19.0 ‰ and 26.7 ‰ and δ¹³C values between −8.5 ‰ and 1.2 ‰, covering ranges of 7.7 ‰ and 9.7 ‰, respectively. We observe differences in δ¹⁸O_c and δ¹³C values among the three foreland basin records, which we relate to their depositional setting (proximal vs. distal) within the megafans and the varying impact of soil water evaporation as well as soil productivity and interaction with atmospheric CO₂ for δ¹⁸O_c and δ¹³C, respectively.

Fontannen δ¹⁸O_c values show a tight range and vary between 19.0 ‰ and 21.3 ‰ with a mean of 19.8±0.4 ‰ (n = 47 individual nodules; Campani et al., 2012). δ¹³C values show a higher variability and range from −6.0 ‰ to −0.6 ‰ with a mean of $-\mathrm{4.0}\pm \mathrm{1.1}$ ‰. Δ₄₇ temperatures average at 36±7 ^∘C (propagated error of the two individual samples) for two pedogenic carbonate samples (MC965b, MC961; Methner et al., 2020) covering 14.5–14.9 Ma. Reconstructed δ¹⁸O_w values for meteoric water range from −7.1 ‰ to −4.7 ‰ with a first quartile mean δ¹⁸O_w value of $-\mathrm{6.5}\pm \mathrm{0.0}$ ‰.

Jona δ¹⁸O_c values range from 20.9 ‰–26.7 ‰, with a mean of 22.9±1.5 ‰ (n=100). δ¹³C values range from −8.5 ‰–1.2 ‰ with a mean of $-\mathrm{3.8}\pm \mathrm{2.1}$ ‰. The mid-fan Jona section shows a wide range of δ¹⁸O and δ¹³C values that cover 5.8 ‰ (δ¹⁸O_c) and 9.7 ‰ (δ¹³C). Δ₄₇-derived carbonate formation temperatures were obtained for two pedogenic carbonate samples (18EK192b, 18EK172a) at 14.5 and 14.9 Ma and yield an average temperature of 30±5 ^∘C. Calculated δ¹⁸O_w values range from −6.3 ‰ to −0.7 ‰ with a first quartile mean δ¹⁸O_w value of $-\mathrm{5.8}\pm \mathrm{0.2}$ ‰.

Compared to the Jona and Fontannen sections pedogenic carbonates from the distal Aabach section yield the highest δ¹⁸O_c values ranging from 23.0 ‰ to 25.9 ‰ with an average of 24.4±0.7 ‰ (n=57). δ¹³C values vary between −6.2 ‰ and −0.1 ‰ with an average of $-\mathrm{2.3}\pm \mathrm{1.7}\mathrm{‰}$. We obtained a Δ₄₇ temperature of 33±2 ^∘C (12EK096a) at 15.6 Ma. Consequently, the reconstructed δ¹⁸O_w values range between −3.8 ‰ and 0.0 ‰ with a first quartile mean δ¹⁸O_w value of $-\mathrm{3.1}\pm \mathrm{0.3}\mathrm{‰}$.

4.2 Paleoclimate simulation of δ¹⁸O_w data

Simulated Northern Alpine Foreland Basin (low-elevation) δ¹⁸O_w values for the topography experiments range from −4.2 ‰ to −5.0 ‰ for the NoAlps experiment, from −4.4 ‰ to −5.3 ‰ (Alps50), from −5.1 ‰ to −6.5 ‰ (Alps100), and from −5.5 ‰ to −7.6 ‰ for the Alps150 experiments (Fig. 5b). For all topography experiments δ¹⁸O_w values decrease with increasing elevation. Simulated central Alpine (high-elevation) δ¹⁸O_w values range from −4.0 ‰ to −4.9 ‰ for NoAlps, from −4.7 ‰ to −6.1 ‰ (Alps50), from −6.0 ‰ to −9.4 ‰ (Alps100), and from −7.6 ‰ to −12.0 ‰ for the Alps150 experiment (Fig. 5b). For the distant region all topography simulations show slightly higher δ¹⁸O_w values when compared to the Northern Alpine Foreland Basin (low-elevation) setting and range from −3.6 ‰ to −5.0 ‰.

5 Discussion

5.1 Swiss Molasse Basin near-sea-level paleosol carbonate δ¹⁸O record

δ¹⁸O_c values of pedogenic carbonate reflect the δ¹⁸O_w values of soil water and prevalent soil temperatures during carbonate formation. Soil water δ¹⁸O_w values typically decrease with elevation but may be biased to higher δ¹⁸O_w values in response to enhanced soil water evaporation (Quade et al., 2007b, 2020). Disentangling the impacts of temperature, elevation, and (soil water) evaporation on the δ¹⁸O_c values of soil carbonates is essential for reconstructing a reliable near-sea-level δ¹⁸O_w record. In this context we want to stress that δ¹⁸O and (Δ₄₇) values of terrestrial carbonates can be affected by diagenetic alterations due to increased burial temperatures (Quade et al., 2020). Maximum erosion estimates for the OSM in the central and eastern Swiss Molasse Basin range from 350 to ∼2000 m (e.g., Cederbom et al., 2011; Schegg and Leu, 1998). Associated burial temperatures (Burkhard and Kalkreuth, 1989; Schegg, 1992) and duration of burial suggest that our estimates of Swiss Molasse Basin carbonate δ¹⁸O and Δ₄₇ compositions have not been affected by burial alterations (for more details see Sect. S4 in the Supplement).

The Swiss Molasse Basin (SMB) foreland δ¹⁸O_c records show significant differences along the slope of the SMB megafans, with the lowest mean δ¹⁸O_c value at the proximal fan position (19.8±0.4 ‰ at Fontannen), an intermediate value at the mid-fan position (22.9±1.5 ‰ at Jona), and the highest at the distal fan position (24.4±0.7 ‰ at Aabach; Fig. 2). Environmental differences between the fan positions are not surprising as the Miocene sediments were deposited in large alluvial megafans which prograded from the Alpine front towards the basin center covering a present-day downstream distance of ca. 30 km (Garefalakis and Schlunegger, 2019). Accounting for post-Miocene crustal shortening, we consider this a minimum estimate for the Miocene length of the Alpine megafans. These megafans were fed by large drainage networks that led to high sediment discharge and aggradation of several hundred meters of detrital material near the fan apex, resulting in an elevation gradient towards the orogenic front. Based on sedimentological data for late Oligocene deposits of the Rigi megafan (30–25 Ma) in central Switzerland, fan surface slopes ranged between approximately 0.2^∘ at the base and 0.9^∘ at the top (Garefalakis and Schlunegger, 2018). Using a time-averaged fan surface slope of 0.6±0.2^∘ and a fan length of 30 km the proximal part of such a fan would be at elevations of 300±100 m above the Miocene sea level. Because of this internal elevation gradient, foreland basin records from proximal locations (e.g., Fontannen; Campani et al., 2012) were deposited at higher elevations compared to sections from more distal fan sites, thus not representing near-sea-level δ¹⁸O_w values.

The distal Aabach section was hence closest to the basin axis, which was near or at sea level at that time (Kuhlemann and Kempf, 2002). The δ¹⁸O_c and δ¹³C values at Aabach, however, show a strong positive correlation (Fig. S3 in the Supplement), and they are additionally associated with high δ¹³C values. Because significant C₄ vegetation was absent at that time (Cerling and Quade, 1993; Tipple and Pagani, 2007), these values indicate low soil respiration rates during the formation of pedogenic carbonate (Quade et al., 2007b, and references therein). High δ¹³C values of soil carbonate are common during high plant water stress, low soil productivity, and low soil respiration and high soil evaporation rates (e.g., Breecker et al., 2009; Caves et al., 2016). Furthermore, the high δ¹³C values at Aabach are consistent with a shift from humid to warm and semiarid conditions after 25 Ma, inferred from paleofloral records for the Swiss Molasse Basin (Berger, 1992). Evaporative effects caused by local aridity and enhanced soil water evaporation may shift carbonate δ¹⁸O_c values to higher values (Cerling, 1984; Cerling and Quade, 1993). These processes likely account for the positive covariance of the δ¹⁸O_c and δ¹³C values, which has also been observed in modern soils (Cerling, 1984; Cerling and Quade, 1993). Consequently, we consider the rather high and highly variable Aabach pedogenic carbonate δ¹⁸O_c values to be biased by varying soil water evaporation that translates into a positive bias in reconstructed soil water δ¹⁸O_w values.

For our paleoelevation reconstruction, we thus consider the mid-fan Jona section to be a conservative estimate for near-sea-level δ¹⁸O_w values as (a) its fan position places the Jona section clearly down-slope from the proximal Fontannen section, (b) combined δ¹⁸O_c and δ¹³C values show no indication for enhanced soil water evaporation, and (c) using the mid-fan Jona section underestimates rather than overestimates reconstructed paleoelevations, as unlike Aabach, it is not located closest to the basin center. Consequently, we consider the variability in δ¹⁸O_w values of the Jona section to reflect secular changes in environmental conditions. To further reduce potential (soil) evaporation bias we only select the first quartile mean δ¹⁸O_c values (δ¹⁸O${}_{\text{c,\hspace{0.17em}lowest\hspace{0.17em}1/4}}=\mathrm{21.4}\pm \mathrm{0.2}\mathrm{‰}$ compared to mean δ¹⁸O${}_{\mathrm{c}}=\mathrm{22.7}\pm \mathrm{1.4}\mathrm{‰}$) to establish our best estimate for a mid-Miocene near-sea-level δ¹⁸O_w value.

5.2 Near-sea-level precipitation δ¹⁸O estimate

We conducted clumped isotope (Δ₄₇) thermometry on pedogenic carbonates from all three SMB sections in order to obtain estimates of carbonate formation temperatures and δ¹⁸O_w values of meteoric water. Δ₄₇ temperatures show high reproducibility and reach values of 36 ^∘C (Fontannen; n=2 samples), 30 ^∘C (Jona; n=2 samples), and 33 ^∘C (Aabach; n=1 sample), consistent with warm season pedogenic carbonate formation during the warm temperatures of the mid-Miocene Climatic Optimum (Methner et al., 2020). Using 21.4±0.2 ‰ (first quartile mean Jona section) as a near-sea-level carbonate δ¹⁸O_c value, this results in δ¹⁸O${}_{\mathrm{w}}=-\mathrm{5.8}\pm \mathrm{1.2}\mathrm{‰}$ as our best estimate for Miocene near-sea-level δ¹⁸O_w values of meteoric water. Collectively, these data revise previous estimates of near-sea-level δ¹⁸O_w based on the proximal Fontannen section (δ¹⁸O${}_{\mathrm{w}}=-\mathrm{8.9}\pm \mathrm{0.5}$ for 21 ^∘C; Campani et al., 2012), which were calculated with mean annual temperatures derived from paleobotanical analysis. The difference from the revised value is attributed to a more suitable choice of the near-sea-level record (mid-fan vs. proximal fan position) and measured (instead of estimated) carbonate formation temperatures amounting to +1.2 ‰ and +1.9 ‰, respectively.

The reconstructed near-sea-level Jona δ¹⁸O_w value is consistent with δ¹⁸O_w values from mid-Miocene SMB volcanic ash horizons (Bauer et al., 2016). Using a mineral–water fractionation temperature of 30 ^∘C (measured Δ₄₇ temperature for Jona) and a smectite–water oxygen isotope fractionation factor from Sheppard and Gilg (1996) these ash layers (Ries, Küsnacht, and Urdorf) reveal mean δ¹⁸O_w values between −6.1 ‰ and −2.9 ‰, which are equal to or even higher than our conservative near-sea-level estimate of δ¹⁸O${}_{\mathrm{w}}=-\mathrm{5.8}\pm \mathrm{1.2}\mathrm{‰}$ based on SMB Jona pedogenic carbonate. Furthermore, the near-sea-level Jona δ¹⁸O_w estimate is in good agreement with modern summer month (June–July–August) δ¹⁸O_w data from GNIP stations in Switzerland (∼250–600 m a.s.l.), which range between $\sim -\mathrm{8}$ ‰ and −2 ‰.

Figure 6 Paleoelevation calculation of the Miocene Central Alps. (a) The difference of precipitation δ¹⁸O_w values (Δ(δ¹⁸O)) between the high-elevation (Simplon Fault Zone; Campani et al., 2012) and the low-elevation sites (Fontannen, Jona, Aabach) is a measure of the respective elevation difference. We use a conservative approach to calculate the paleoelevation and only consider the lowest 25 % δ¹⁸O_w values (n=10) of the mid-fan-situated Jona section as most suitable to serve as a low-elevation reference baseline (see text for discussion). Swiss Molasse Basin precipitation δ¹⁸O_w values were calculated from pedogenic carbonate δ¹⁸O_c values using carbonate clumped isotope (Δ₄₇) formation temperatures. The figure shows foreland data from this study (Fontannen, Jona and Aabach) and from Campani et al. (2012; Fontannen, 2012). Also shown are quartz vein fluid inclusion δ¹⁸O_w values from the Lukmanier region (Swiss Alps; Sharp et al., 2005), δ¹⁸O_w values derived from volcanic ashes (Bauer et al., 2016), and the modeled range for high- and low-elevation precipitation δ¹⁸O_w values for the present-day Alps150 (ECHAM5-wiso iGCM run; Botsyun et al., 2020). (b) Jona Δ(δ¹⁸O) value and inferred elevation differences (Δz; m) (horizontal green rectangle) are based on different isotope lapse rates. Black line: −2.0 ‰ km⁻¹ (Campani et al., 2012); brown line: −2.1 ‰ km⁻¹ (Poage and Chamberlain, 2001); yellow line: −2.4 ‰ km⁻¹ (Botsyun et al., 2020); grey line: Currie et al. (2005), Rowley and Garzione (2007). The uncertainty in Δ(δ¹⁸O) is ±1.5 ‰. The present-day Alpine δ¹⁸O_w lapse rate of $-\mathrm{2.0}\pm \mathrm{0.04}\mathrm{‰}$ km⁻¹ translates the Jona Δ(δ¹⁸O) value of $-\mathrm{8.8}\pm \mathrm{1.5}\mathrm{‰}$ to a paleoelevation difference of 4420±770 m (dashed lines in black) between the foreland basin and the intra-Alpine SFZ.

5.3 Mid-Miocene stable isotope paleoaltimetry of the Central Alps

We contrast the near-sea-level ∼14 Ma δ¹⁸O_w record from the Jona section with a high-elevation δ¹⁸O_w record from the Simplon Fault Zone at 14.5 Ma (Fig. 4). Relative differences in δ¹⁸O_w between near sea level and the high-Alpine Simplon Fault Zone are expressed by Δ(δ¹⁸O_w) and calculated by subtracting the SMB δ¹⁸O_w value from the SFZ δ¹⁸O_w value.

Given δ¹⁸O_w values for near-sea elevation of $-\mathrm{5.8}\pm \mathrm{1.2}\mathrm{‰}$ and the high-Alpine Simplon Fault Zone of $-\mathrm{14.6}\pm \mathrm{0.3}\mathrm{‰}$ (Campani et al., 2012), respectively, Δ(δ¹⁸O_w) equals $-\mathrm{8.8}\pm \mathrm{1.5}\mathrm{‰}$ (Fig. 6a). For the paleoelevation calculation, we explored different oxygen isotope lapse rates (Fig. 6b), including the present-day Alpine lapse rate based on long-term meteorological station data (−2.0 ‰ km⁻¹; Campani et al., 2012), the present-day surface-water- and precipitation-based oxygen isotope lapse rate for Europe (−2.1 ‰ km⁻¹; Poage and Chamberlain, 2001), a thermodynamic-model-based lapse rate that tracks the isotopic composition of water vapor along precipitation trajectories (Currie et al., 2005; Rowley et al., 2001; Rowley and Garzione, 2007), and the output of the isotope-enabled paleoclimate model ECHAM5-wiso (−2.4 ‰ km⁻¹; Botsyun et al., 2020). Independent of the choice of these lapse rates a Δ(δ¹⁸O_w) value of $-\mathrm{8.8}\pm \mathrm{1.5}$ ‰ requires significant orographic rainout and hence elevated topography. For instance, applying the different lapse rates to a Δ(δ¹⁸O) value of −8.8 ‰ yields paleoelevations in the range of ∼3670 to 4400 m (Table S7 in the Supplement). The inferred paleoelevation estimate represents the mean elevation of the mid-Miocene paleo-catchment of precipitation in the Simplon area. To keep comparability with the previous stable isotope paleoaltimetry study of Campani et al. (2012), we use the modern Alpine oxygen isotope lapse rate of $-\mathrm{2.0}\pm \mathrm{0.04}$ ‰ km⁻¹. Based on the Δ(δ¹⁸O_w) value of $-\mathrm{8.8}\pm \mathrm{1.5}$ ‰ this is consistent with a relative elevation difference between the Swiss Molasse Basin and the Simplon Fault Zone of $\mathrm{\Delta }z\left(\mathrm{m}\right)=\mathrm{4400}\pm \mathrm{770}$ m (Fig. 6b). Assuming a much more conservative choice of isotope lapse rate (−2.4 ‰ km⁻¹; Botsyun et al., 2020) the same Δ(δ¹⁸O_w) value results in $\mathrm{\Delta }z\left(\mathrm{m}\right)=\mathrm{3670}\pm \mathrm{625}$ m. Differences in atmospheric temperature gradients in the past will inevitably affect oxygen and hydrogen isotope lapse rates, and climate change may play a decisive role when trying to assess paleoelevation of ancient mountain ranges. Especially for warmer periods of the Earth's past it is shown that warm conditions yield shallower lapse rates (Poulsen and Jeffery, 2011; Rowley and Garzione, 2007). All existing data point to warmer conditions in central Europe during the mid-Miocene Climatic Optimum when compared to today (Böhme, 2003; Methner et al., 2020; Mosbrugger et al., 2005; Pound et al., 2012). The mid-Miocene should logically have been characterized by lower isotopic lapse rates when compared to the present. Thus, application of the present-day Alpine lapse rate most likely underestimates mid-Miocene paleoelevations. Our estimated mid-Miocene (∼14.5 Ma) surface elevation of 4440 m for the region surrounding the Simplon Fault Zone is in good agreement with the proposed mid-Miocene minimum elevation of 5000 m given by Sharp et al. (2005) and the >5000 m found for the Oligocene–Miocene boundary by Jäger and Hantke (1983, 1984), and it places our estimated paleoelevation in the higher spectrum of hitherto derived paleoaltimetry estimates. For comparison, peaks in the region today reach elevations in excess of 3500 m.

5.4 Modeled low-elevation δ¹⁸O_w estimates

Collectively, the topography-modulated model experiments indicate the following. First, modeled Northern Alpine Foreland Basin δ¹⁸O_w values for the Alps50, Alps100, and Alps150 cases range from −4.4 ‰ to −7.6 ‰ and are therefore in good agreement with the reconstructed low-elevation δ¹⁸O_w estimate of −5.8 ‰ from the Swiss Molasse Basin. Furthermore, the ECHAM5-wiso simulations reveal that a Δ(δ¹⁸O_w) value of −8.8 ‰ between the high-elevation and the low-elevation records requires a mountain range with an elevation >150 % of present-day Alpine mean topography (Figs. 5b and 6b). We therefore interpret the observed difference in δ¹⁸O_w between the high-elevation and low-elevation records as resulting from a higher-than-modern mid-Miocene elevation of the region surrounding the Simplon Fault Zone. This interpretation is supported by the predicted differences in the 150 % to 100 % topographic change explored in the paleoclimate simulations.

5.5 High (and highly variable) mid-Miocene Central Alps?

Our paleoelevation estimate of >4000 m for the mid-Miocene Simplon area is surprisingly high when compared to low erosion rates between 16 and 12 Ma for the Swiss Alps, as inferred from sediment budget studies (Kuhlemann et al., 2001). At the first glance, this may seem contradictory. In the following we explain how the formation of a complex Alpine topography can resolve this contradiction. We propose a scenario in which the pre- to mid-Miocene Central Alps changed from a rather simple, and most likely cylindrical, to a more diverse landscape with a complex topography at ∼15 Ma at the latest. This change may have been mainly driven by rapid exhumation of deep-seated core complexes, followed by a rearrangement of the drainage system.

The history of central Alpine topography is considered to have started no later than the late Oligocene–early Miocene after slab break-off of subducted oceanic lithosphere ∼32–30 Myr ago (Schlunegger and Kissling, 2015). Increased sediment discharge in the Swiss Alps has been documented for several periods (30–28, 24–21, and 18–17 Ma) and linked to exhumation and surface uplift events as well as an increase in topographic relief (Kuhlemann et al., 2001). A transition from a poorly dissected, plateau-like orogen to a rugged mountain range with steep relief is proposed for the Alps at ∼27 Ma as a possible response to slab break-off (Garefalakis and Schlunegger, 2018), promoting the buildup of high mountain peaks and establishing a N–S-oriented drainage network with a drainage divide situated close to the area of inferred slab break-off (Kuhlemann et al., 2001). The situation then changed from ca. 20 Ma onward, when basin- and orogen-related reorganizations are likely to have affected the topographic evolution on a regional scale since at least the early Miocene. These include exhumation in the Lepontine and Aar regions (Herwegh et al., 2017; Schlunegger and Willett, 1999), a hypothesized reversal in slab polarity (Kissling et al., 2006; Lippitsch, 2003) or slab delamination (Handy et al., 2021) beneath the Eastern Alps, a switch in the regional tilt of the Swiss Molasse Basin with an associated change in the basin-axial discharge direction at ∼17 Ma to a more complex pattern thereafter (Berger et al., 2005; Kuhlemann and Kempf, 2002; Kühni and Pfiffner, 2002; Garefalakis and Schlunegger, 2019), a change in sediment provenance in the SMB deposits (e.g., Anfinson et al., 2020; Von Eynatten et al., 1999), and the beginning of the reorganization of the drainage network into an orogen-parallel E–W-oriented system (Bernard et al., 2021; Kühni and Pfiffner, 2002; Schlunegger et al., 2001). Four observations point to high spatial topographic heterogeneity in the mid-Miocene Central Alps: (1) the rapid and nearly vertical rise of the Aar Massif at ∼20 Ma is associated with a rearrangement of the drainage network, leading to a shift of the drainage divide towards the uplifted crystalline block (Bernard et al., 2021; Kühni and Pfiffner, 2001; Schlunegger et al., 2001), which is made up of basement rocks of the European plate comprising granites and granodiorites (Herwegh et al., 2020) with the lowest erodibility in the Central Alps (Kühni and Pfiffner, 2002). A shift of the drainage divide towards rocks with low erodibility implies a decrease in overall erosion and denudation rates by reducing the erosional potential for several million years. This potentially resulted in increased surface uplift rates and thus in a local rise of topography (Bernard et al., 2021). (2) Associated rerouting of streams by rise of such a crystalline block in the center of an orogen is further consistent with the initiation of an orogen-parallel-oriented drainage system similar to the change from an orogen-normal N–S-directed to an orogen-parallel E–W-oriented drainage pattern that occurred during the early Miocene. Such a change will inevitably affect the overall relief structure. We need to add, however, that an orogen-parallel pattern was not yet fully established in the mid-Miocene, as studies show that the Napf megafan was still receiving deposits from the Aar region at this time (Stutenbecker et al., 2019). (3) Interestingly, at the same time, rapid exhumation of the nearby Lepontine area through slip along the Simplon Fault Zone was associated with a reduction in sediment supply to the Swiss Molasse Basin (Kuhlemann et al., 2001). Normal fault-induced rapid exhumation of the Lepontine dome footwall, with the highest rates at ∼20 Ma, renders it likely that mean elevation decreased in the Lepontine area in response to rapid tectonic exhumation and associated unroofing. (4) The stable isotope paleoaltimetry data from the Simplon Pass region, however, indicate rather high >4000 m coeval surface elevations only 45 km to the west of the Lepontine Dome.

The co-existence of regions with different elevations on a small spatial scale within the Miocene Central Alps points to a regional landscape characterized by significant topographic complexity that was sensitive to alterations of the local drainage network and that was accompanied by laterally variable exhumation rates. Therefore, on a regional scale the Miocene Alps were most probably characterized by a heterogeneous and spatially transient topography with high elevations locally exceeding 4000 m. A heterogeneous, non-cylindrical Alpine tectonic structure with corresponding topography has implications, e.g., for the spatially variably impact of slab dynamics on Earth surface process by provoking perturbations of surface loads as well as for assumptions made in orogen-scale denudation and landscape development models.

6 Conclusions

Our revised stable isotope paleoaltimetry estimates indicate high (>4000 m) paleoelevations for the mid-Miocene (15.5–14.0 Ma) central Alpine region surrounding the Simplon Fault Zone. This result is based on (a) the pedogenic carbonate record best representing the near-sea-level δ¹⁸O_w values in the Swiss Molasse Basin, (b) quantification of carbonate formation temperatures by Δ₄₇ thermometry (rather than assuming those from other terrestrial temperature proxies), and (c) ECHAM5-wiso modeled Northern Alpine Foreland Basin δ¹⁸O_w values that cover the same range as reconstructed near-sea-level mean δ¹⁸O_w values. Based on coupled δ¹⁸O_c values and clumped isotope (Δ₄₇) paleothermometry on SMB pedogenic carbonates (present-day elevation of ∼500–800 m a.s.l) we estimate a mean δ¹⁸O_w value for near-sea-level precipitation of $-\mathrm{5.8}\pm \mathrm{1.2}\mathrm{‰}$. Our conservative paleoaltimetry estimation yields a difference in paleoelevation of 4400±770 m between the Swiss Molasse Basin and the region surrounding the Simplon Fault Zone. We propose a change from a uniform pre- to mid-Miocene landscape to a more complex one with highly variable topography at the latest in the middle Miocene. We link this change with the exhumation of the Aar Massif at ∼20 Ma and the associated reorganization of the Alpine drainage network. Subsequent paleoaltimetry studies will have to focus on changes in surface elevation through space and time. The combination of stable isotope paleoaltimetry and analysis of regional exhumation patterns in conjunction with paleoclimate modeling can contribute to establishing not only absolute elevations, but also to reconstructing past elevation difference along and across the orogen by linking long-wavelength subsurface with short-wavelength surface processes which affect the Alpine topographic evolution. Further experiments with isotope-tracking global atmospheric circulation models that test time-specific boundary conditions of the middle Miocene will represent key elements in verifying observed changes in δ¹⁸O_w and will provide important information on links and feedbacks between surface uplift and associated climate change.