Quantifying Multiple Erosion Events in the Distal Sector of the Northern Alpine Foreland Basin (North-Eastern Switzerland), by Combining Basin Thermal Modelling with Vitrinite Reﬂectance and Apatite Fission Track Data

: This work quantiﬁes the amount of erosion associated with the Cretaceous and Miocene erosional unconformities recognised in the distal part of the Northern Alpine Foreland Basin (NAFB), north-eastern Switzerland. To achieve this goal, the basin thermal modelling approach is applied, calibrated by two different sets of data collected in previous studies: vitrinite reﬂectance (%R o ) and the temperature estimated from apatite ﬁssion tracks (AFT) data modelling. The novelty of this approach is the possibility to constrain the timing and magnitude of multiple erosion events by integrating thermal modelling with thermochronologic data. Combining these two methods allows the erosional events to be separated which would not be possible using only irreversible paleothermometers, such as vitrinite reﬂectance data. Two scenarios were tested, based on the data of two published thermochronology studies. For the Cretaceous unconformity, similar results are obtained for the two scenarios, both indicating that the deposition and the subsequent complete erosion of Lower Cretaceous deposits, in the order of 500–1300 m, depending on the area, are necessary, in order to attain the temperatures estimated by the thermal history modelling of AFT data. Thus, a depositional hiatus for this period is not likely. For the Miocene-Quaternary unconformity, the magnitude of erosion calculated for the two scenarios differs by 300–1400 m, depending on the AFT data considered. The two scenarios lead to a different evaluation of the subsidence and uplift rate of the study area, thus to a different interpretation of the tectono-stratigraphic evolution of this distal sector of the NAFB.


Introduction
This work focuses on the area located in the distal part of the Northern Alpine Foreland Basin (NAFB) (Figure 1a), in the north-eastern sector of Switzerland and extended throughout the Cantons of Aargau, Zürich and Thurgau (Figure 1b). Several studies were conducted in this area, for characterising the shallow and deep subsurface, initially for coal and hydrocarbon exploration, and later for nuclear waste storage [1]. The interest in this area has been recently renewed for both shallow and deep geothermal energy exploration and geological storage of CO 2 [2]. For proper and safe use of the subsurface, it is necessary to understand the present and past temperatures attained by the rocks, as well as to understand the tectono-stratigraphic evolution of the area. For these purposes a fundamental step is to estimate the thickness of the deposits eroded during the evolution of the basin. These data from the stratigraphic record enable one to (i) define the original  Several studies have been done to estimate the magnitude of the erosion in the study area. Several methods have been used such as analysis of vitrinite reflectance data (%Ro) [3][4][5][6], stratigraphic and porosity extrapolation [7]; subsidence analysis [8], clay mineral transformation [9], seismic interval analysis [7], and apatite fission tracks [6,[10][11][12][13][14][15].
As yet there is no consensus for the amount and timing of erosion. One of the major difficulties is to correctly estimate the amount of missing section for each of the main erosional events recognised in the stratigraphic record of the NAFB, represented by Permian-Triassic, Cretaceous-Eocene and Miocene-Quaternary unconformities.
In this work, we use a basin thermal modelling approach to estimate the erosion in the study area [16][17][18][19][20]. The advantage of basin thermal modelling compared to other methods is that it integrates by a software multiple basin components, such as (i) the sediment decompaction, (ii) the burial/uplift timing, (iii) the lithological conductivity of the basin infill, (iv) the circulation of fluids, and (v) the basin basal heat flow. To calibrate the basin thermal modelling results, the theoretical thermal conditions estimated by the software are fitted to paleo-thermal data measured in the rocks, which indicate the maximum paleotemperature attained by the stratigraphic record at a given depth. However, when the basin is affected by multiple burial and erosional events, as is the case of the Northern Alpine Foreland Basin, to provide reliable estimate of the maximum temperature attained at each burial event and the magnitude of the following erosion is challenging. This is because the thermal transformation of the proxies commonly used for calibration (clay minerals, vitrinite reflectance, biomarker composition, etc.) are irreversible, which means that if the basin is affected by several thermal pulses, only the hottest is registered. Thus, it is not possible to reconstruct whether the thermal transformation/maturation of the sedimentary record occurred gradually or primarily reflects punctual strong burial and/or thermal events. As a result, there will be multiple basin evolutions that can match the observed paleothermometers. To overcome this limitation, thermochronometers able to estimate the timing of the main relevant heating events, such as apatite fission tracks data (AFT), were combined with the paleothemometers commonly used to calibrate the thermal modelling results (vitrinite reflectance data, %Ro). Using this approach, the magnitude of the erosion for each erosional unconformity recognised in the NAFB was quantified.

Geological Setting
The studied area is part of the Swiss Northern Alpine Foreland Basin (NAFB), which formed in three major geodynamic phases (Figure 2), represented by three stratigraphic mega-sequences [21]: (i) the basement, (ii) the Mesozoic substratum and (iii) the infill of the Alpine foreland basin. The basement formed during the Variscan and pre-Variscan orogenies, is locally segmented in elongated grabens. The graben formed in the Late Carboniferous to Early Permian extensive tectonic regime that accompanied the collapse of the Variscan orogenic belt, and filled with syn-tectonic continental sediments [22,23]. An intense magmatic and hydrothermal activity characterises this period, which involved significant changes in the thermal regime [23,24]. The stratigraphic record of the Mesozoic substratum represents a passive margin sequence related to the opening of the Tethys Ocean, which was located between the Eurasian and African plates [25]. The Triassic and Jurassic sequences are uniformly represented throughout the basin, whereas the Cretaceous deposits are represented only in the western sector of the NAFB, limited to an Early Cretaceous age and are completely absent in the eastern sector [4,8,22,26], where the study area is located. The third and uppermost mega-sequence is the foreland infill, which gives the wedge shape to the basin. It formed in response to the Alpine collision, as a consequence of the flexural bending of the European plate under the thrust and crustal thickening of the Alpine orogeny [21,[27][28][29][30][31]. The stratigraphic record of this sequence is composed of shallowing-upward cycles of shallow marine to fluvial-lacustrine deposits, Rupelian to Serravalian in age [28,[32][33][34][35][36]. In the Miocene, an important displacement of the Alpine thrust front, 30-40 km north-westward, along the Triassic evaporitic detachment units, formed the Jura fold-andthrust belt chain [21,[37][38][39]. By consequence, the most external part the foreland basin was gradually exhumed and eroded. The rest of the foreland basin remains relatively undeformed, whereas subsidence and sedimentation most likely continued up to around 5 Ma [40], when the entire NAFB was uplifted and partially eroded [12].

Basin Erosional Unconformities
In the Northern Alpine Foreland Basin (NAFB) three main erosional unconformities are described. The first one is the pre-Mesozoic unconformity, which separates the Mesozoic sequence from the previous Variscan and post-Variscan deposits ( Figure 2). The In the Miocene, an important displacement of the Alpine thrust front, 30-40 km north-westward, along the Triassic evaporitic detachment units, formed the Jura fold-andthrust belt chain [21,[37][38][39]. By consequence, the most external part the foreland basin was gradually exhumed and eroded. The rest of the foreland basin remains relatively undeformed, whereas subsidence and sedimentation most likely continued up to around 5 Ma [40], when the entire NAFB was uplifted and partially eroded [12].

Basin Erosional Unconformities
In the Northern Alpine Foreland Basin (NAFB) three main erosional unconformities are described. The first one is the pre-Mesozoic unconformity, which separates the Mesozoic sequence from the previous Variscan and post-Variscan deposits ( Figure 2). The Paleozoic rocks are composed of pre-Carboniferous metamorphic rocks and Carboniferous and Permian sediments. Uplift and erosion events have been suggested to have occurred in the Early Permian (Saalian phase), as testified by the large amount of siliciclastic conglomeratic and breccia alluvial debris series deposited at that time in the grabens (e.g., Rotliegend unit), indicating an active relief between troughs and shoulders areas [41]. The second unconformity is located at the top of the Mesozoic sequence, where the Oligocene-Miocene foreland units lie unconformably (Figure 2) on the Mesozoic sequence. This is observed throughout the entire NAFB, extending to the western part of the German Molasse Basin [8,26,33]. The stratigraphic gap of this unconformity increases from west to east. The Lower Cretaceous deposits are preserved in the western sector of the Swiss Foreland Basin, whereas they are completely absent in the eastern part [8,26]. In the Zürich High area (extending from northeast Switzerland to the southern Rhine valley) the Mesozoic series was eroded down to the Oxfordian [33]. The third unconformity surface is located at the top of the Cenozoic foreland sequence (Molasse deposits, Oligocene to Miocene in age), separating these from the Quaternary units ( Figure 2). This last unconformity registers the uplift and exhumation of the foreland basin during the last Alpine thrusts pulses and was accompanied by erosion of the uppermost units [4,42,43]. This last unconformity is diachronous throughout the entire NAFB basin, as the Quaternary deposits lies on Molasse deposits which get older going westward [28]. Toward the Jura fold-and-thrust belt domain the unconformity time gap is greater, with Quaternary deposits often lying directly on top of Jurassic or Triassic rocks. In these areas, the pre-Eocene and pre-Quaternary unconformities are hardly distinguishable.

Apatite Fission Tracks (AFT) Investigations in Northern Switzerland
Several AFT studies have been carried out in the study area [6,10,14,15] thanks to investigations of deep crystalline basement rocks, drilled by the Swiss National Cooperative for the disposal of radioactive waste (Nagra), to investigate the properties of the basement [44]. One of the most referenced studies in the area is that of Muzerek et al. [6] which included a series of AFT data obtained by Nagra [45] on a large number of deep wells (Benken-1, Herdern-1, Weiach-1, Böttstein-1, Schafisheim-1 wells, Figure 1b). These authors found two successive stages of heating. The first one was in the Early Cretaceous, associated to burial and to a heat flow peak in this period the second was in the Miocene, also associated with burial with the formation of the Northern Alpine foreland basin and the deposition of the Molasse foreland units. Temperatures attained by the Miocene event exceeded those of the Cretaceous in the more proximal part of the foreland but not in the distal part (Jura region), where burial depths were lower. Timar-Geng et al. [10] suggest for the same area (Kaisten-1, Riniken-1, Leuggern-1 wells, Figure 1b) two important cooling events in the Early Cretaceous and a heating event in the Eocene, followed by cooling to present-day temperatures. The Eocene heating event is mostly ascribed to deep hydrothermal fluid circulation, associated to the Upper Rhine Graben rift and volcanic activity. In Kuhlemann and Rhan [14] a Late Jurassic/Early Cretaceous cooling event followed by isothermal conditions throughout most of the Cretaceous-early Cenozoic is discerned from thermal modelling AFT data in the Benken-1 well. A second heating event is recognised in the Oligocene with the deposition of the Molasse foreland units. Finally the most recent AFT study carried out on the basement and Paleozoic sediments of the area (Böttstein-1, Weiach-1, Riniken-1, and Leuggern-1 wells), Villagomez et al. [15], shows a heating event in the Early Cretaceous followed by a slow cooling in the Late Cretaceous continuing into/through the Pliocene.
The consensus data resulting from these works is a cooling event in the Cretaceous, interpreted as the effect of the lithosphere thermal relaxation after an important heating event in the Late Jurassic/Early Cretaceous. However, it is not clear if these heating and cooling events are associated to the deposition and subsequent erosion of Cretaceous deposits [6] or caused by regional-scale hydrothermal fluid migration from the nearby Black Forest [46], which could justify a depositional hiatus interpretation for the Jurassic-Eocene unconformity [5,21,26,34]. The main discrepancy among these studies is related to the Oligo-Miocene heating event, which is recognised only in [6,14] and associated to the important burial (>1000 m) of the foreland basin followed by a correspondent exhumation and erosion event. In [10,15] the Miocene heating peak is not identified, inferring a less significant foreland burial. In [10] this difference is explained by the palaeogeographically more distal settings of the studied area, located slightly northward with respect to the wells studied in [6]. However, if we observe in the map the distribution of the wells studied by these authors (Figure 2), this statement is not a constraint; the Riniken-1, analysed by [10], is an internal position with respect to the Böttstein-1 well, at a similar position to the Weiach-1 well, both of which were studied in [6]. Furthermore, the Leuggern-1 well (studied by [10]) is very close to the Böttstein-1 well (studied by [6]) (3.5 km apart). The difference among the results and interpretations provided by these studies can be related to different AFT data acquisition and thermal modelling methods, which have evolved over time [15].

Methods and Data
To estimate the amount of the erosion in the NAFB, the basin thermal modelling approach was applied [16][17][18][19][20]. Modelling results are calibrated with paleothermometers (vitrinite reflectance) and thermochronometers (apatite fission track data). Thermal modelling was performed by using the Petromod ® software provided by Schlumberger. The 1D models were reconstructed based on the data from six wells drilled between 1983 and 2006 [6,44] (Figure 3). Wells were chosen by considering the availability of calibration data, which have to include both vitrinite reflectance and apatite fission track data ( Figure 1b). For each well the model inputs consist in the age, thickness and lithology of the stratigraphic units. The facies characterising each stratigraphic unit were defined by considering literature data and well reports [47][48][49][50][51][52] and expressed in lithologies percentage ( Figure 3 and Table 1).
To constrain the model, upper and lower thermal boundary conditions were defined. The upper boundary is the sediment-water interface temperature (SWIT). Paleotemperature distribution maps are automatically calculated by the modelling software, which defines the evolution temperature at sea level considering variations of global mean surface temperature and latitudinal variation of the study area through time [53]. The lower boundary condition is the heat flow at the bottom of the basin. It is related to the geodynamic setting where the basin forms and evolves, which controls, among other processes, the original lithosphere thickness, the stretching forces, the magmatic activity, and the circulation of deep fluids [17,[54][55][56][57].

Basal Heat Flow
The basal heat flow trend through time was estimated by considering the tectonic settings proposed in the literature for formation of the Swiss Foreland Basin and its Mesozoic substratum [4,6,22,23,25,26,54,[58][59][60] (Figure 4). Furthermore, for the extensional events, an estimate of the heat flow values, attained as a consequence of the lithosphere stretching, was deduced by the magnitude of tectonic subsidence registered in the basin [55,61]. For further detail see Omodeo-Salé et al. [62].   The basal heat flow trend through time was estimated by considering the tecto settings proposed in the literature for formation of the Swiss Foreland Basin and its M ozoic substratum [4,6,22,23,25,26,54,[58][59][60] (Figure 4). Furthermore, for the extensio events, an estimate of the heat flow values, attained as a consequence of the lithosph stretching, was deduced by the magnitude of tectonic subsidence registered in the ba [55,61]. For further detail see Omodeo-Salé et al. [62]. An important thermal surge, with heat flow up to 100-120 mW/m 2 is placed in Early Permian, expressing the thinning and delamination of the lithosphere and the m matic activity characterising this period [6,23]. This is followed by a gradual decrease heat flow associated with thermal subsidence, reaching values of 50-60 mW/m 2 throu the Triassic. In the Jurassic the heat flow gradually increased as a consequence of li sophere stretching which affected the European margin during the opening of the Teth Ocean. A maximum heat flow peak of 65-75 mW/m 2 is placed at the Late Jurassic. additional peak in the Early Cretaceous has been proposed by Mazurek et al. [6], wh could be related to further extension of the lithosphere and the formation of oceanic cr in the western Tethys [22,23,25,43,[63][64][65][66]. After the Jurassic-Early Cretaceous therm anomaly, it is assumed that heat flow decreased gradually to values of 50-60 mW/m 2 to the Pliocene, as the basin evolved toward a convergent geodynamic setting, with gradual thickening of the lithosphere, thus a cooling of the geothermal regime. At 2 M an important heat flow peak up to 75-105 mW/m 2 is suggested, in order to be compati with the present-day high geothermal gradient measured in most of the wells of the a [6,67,68]. These high values, compared to what is expected for a foreland basin, are int preted as a consequence of advective heat transfer in the deep crystalline basement, w upflows circulating along the Permo-Carboniferous troughs and faults [6,68].
Preliminary heat flow values were assigned in the model with a range of +/− mW/m 2 around the values reported above. These were then iteratively changed to obt a best fit with the calibration data (vitrinite reflectance and temperature modelled fro AFT data). The best results are indicated by the bold continuous line shown in Figure   3.

Calibration Process and Data
In the first step of the calibration process, the thermal conditions simulated by th mal modelling were validated by vitrinite reflectance data (%Ro) as reported in literatu An important thermal surge, with heat flow up to 100-120 mW/m 2 is placed in the Early Permian, expressing the thinning and delamination of the lithosphere and the magmatic activity characterising this period [6,23]. This is followed by a gradual decrease of heat flow associated with thermal subsidence, reaching values of 50-60 mW/m 2 through the Triassic. In the Jurassic the heat flow gradually increased as a consequence of lithsophere stretching which affected the European margin during the opening of the Tethys Ocean. A maximum heat flow peak of 65-75 mW/m 2 is placed at the Late Jurassic. An additional peak in the Early Cretaceous has been proposed by Mazurek et al. [6], which could be related to further extension of the lithosphere and the formation of oceanic crust in the western Tethys [22,23,25,43,[63][64][65][66]. After the Jurassic-Early Cretaceous thermal anomaly, it is assumed that heat flow decreased gradually to values of 50-60 mW/m 2 up to the Pliocene, as the basin evolved toward a convergent geodynamic setting, with a gradual thickening of the lithosphere, thus a cooling of the geothermal regime. At 2 Ma, an important heat flow peak up to 75-105 mW/m 2 is suggested, in order to be compatible with the present-day high geothermal gradient measured in most of the wells of the area [6,67,68]. These high values, compared to what is expected for a foreland basin, are interpreted as a consequence of advective heat transfer in the deep crystalline basement, with upflows circulating along the Permo-Carboniferous troughs and faults [6,68].
Preliminary heat flow values were assigned in the model with a range of +/−20 mW/m 2 around the values reported above. These were then iteratively changed to obtain a best fit with the calibration data (vitrinite reflectance and temperature modelled from AFT data). The best results are indicated by the bold continuous line shown in Figure 4.

Calibration Process and Data
In the first step of the calibration process, the thermal conditions simulated by thermal modelling were validated by vitrinite reflectance data (%Ro) as reported in literature for the six wells considered (Weiach-1, Böttstein-1, Benken-1; Riniken-1, Schafisheim-1, and Herdern-1, Figure 1b and Annex S1) [6,[69][70][71][72]. For most of the wells, %Ro data show a linear increasing trend with depth and age of the stratigraphic units, as expected in sedimentary basins [54,73,74]. However, in some cases, scattered %Ro values are present, which could be related to altered and/or resedimented vitrinite particles. Unfortunately, a better data interpretation and quality evaluation is not possible, as in the literature the minimum and maximum %Ro values, deviation standard and number of particles measured are not uniformly available for all the wells. To calibrate the thermal conditions simulated by the thermal models, the theoretical vitrinite values computed by the software, applying the Sweeney and Burnham [75] kinetic, were fitted to the measured values, by iteratively changing the variables that influence the temperature in the basin, namely the basal heat flow and the erosion magnitude ( Figure 5). for the six wells considered (Weiach-1, Böttstein-1, Benken-1; Riniken-1, Schafisheim-1, and Herdern-1, Figure 1b and Annex S1) [6,[69][70][71][72]. For most of the wells, %Ro data show a linear increasing trend with depth and age of the stratigraphic units, as expected in sedimentary basins [54,73,74]. However, in some cases, scattered %Ro values are present, which could be related to altered and/or resedimented vitrinite particles. Unfortunately, a better data interpretation and quality evaluation is not possible, as in the literature the minimum and maximum %Ro values, deviation standard and number of particles measured are not uniformly available for all the wells. To calibrate the thermal conditions simulated by the thermal models, the theoretical vitrinite values computed by the software, applying the Sweeney and Burnham [75] kinetic, were fitted to the measured values, by iteratively changing the variables that influence the temperature in the basin, namely the basal heat flow and the erosion magnitude ( Figure 5).  In a second step, the modelling results were calibrated with respect to the temperature reconstructed by the thermal history modelling of the AFT data, performed in the area by previous studies [6,15]. In these works, by means of appropriate software that integrates the age and confined track length, a range of temperatures and timings for the most relevant heating and cooling events of the basin are provided ( Table 2). This results in different interpretations proposed by these authors (see paragraph above). Two scenarios were tested: the first scenario considers a first heating event in the Late Jurassic/Early Cretaceous, followed by cooling until the Eocene-Oligocene, and a second heating event in the Oligo-Miocene, followed by cooling until the present day, as indicated by [6]. The second scenario considers only the Late Jurassic-Early Cretaceous heating event, followed by a gradual cooling until the present day, as indicated by [15]. A description of the main differences between the methods applied by these authors is presented in the Annex S2. For more detail see the methods described in [6,15]. Only the AFT data collected in those wells where %Ro data was also available were considered (Figure 1b). To compare our basin thermal modelling results with the temperature obtained by these two authors for the main heating events (Table 2), temperature versus time curves were extracted for the units located at the same depth as the AFT data considered.

Results
One dimension basin thermal modelling was carried out for the six wells considered (Weiach-1; Böttstein-1; Riniken-1; Benken-1; Schafischeim-1; Herdern-1) and calibrated with both vitrinite reflectance data and the temperature estimated from the thermal modelling of the apatite fission tracks data in [6,15]. To calibrate the model results the amount of erosion and the basal heat flow values first were changed to fit %Ro data ( Figure 5). Second, for each well two scenarios were obtained ( Figure 6) which validate the thermal histories respectively estimated in [6,15]. For each of the two scenarios an estimate of the stratigraphic thickness eroded during the main uplift and erosional phases recognised in the NAFB are provided ( Table 3). The Permian-Triassic erosion event could be evaluated only for the Weiach-1 well, where the Permian and Carboniferous strata have been completely penetrated, and thus sampled throughout for vitrinite reflectance and AFT measurement. In the Riniken-1 well, despite AFT data coming only from the Permian unit [15], paleothermal data along depth are not available for other portions of the Paleozoic sequence. This makes it difficult to recognise a clear geothermal trend for that time. Therefore, an estimate of the magnitude of the Permian-Triassic erosion throughout the entire studied area is not provided.  In the case that, for a single well AFT data were not measured by both [6,15] studies, temperature values were extrapolated from the nearest well where AFT data by the same author and at similar depths were available. As a result, the erosion estimated in those cases is less reliable and should be considered as only an indicative value (indicated by an asterisk in Table 3). Due to the several modelling uncertainties and assumptions (heat flow, calibration data uncertainties, software limitations, etc.), an error bar of hundreds of meters in the results obtained herein can be envisaged (see Section 5.3).

Weiach-1 Well
The Weiach-1 well is located in the distal part of the NAFB, between the Jura fold thrust Belt and the Tabular Jura mountain (Figure 1b). In this well the entire stratigraphic record of the NAFB is represented, from the Quaternary to the Variscan crystalline basement ( Figure 3). The foreland units are only 150 m thick, lying directly on the Upper Jurassic deposits. %Ro data versus depth show an increase of the paleothermal gradient in the Permian and Carboniferous time with respect to the uppermost sequence ( Figure 6). To calibrate this trend, the best result was obtained by assigning a very high heat-flow value (150 mW/m 2 ) during the Permo-Carboniferous time, and a low erosion thickness (200 m) during the Permo-Triassic unconformity. Assigning a large amount of erosion (>1000 m) and a lower heat flow does not provide a satisfactory calibration. Since the Permo-Carboniferous stratigraphic record is not represented in other wells, it is not possible to determine if this high heat flow is the result of a regional paleo-thermal anomaly, or if it is the effect of recent deep hot fluids circulation in the specific Weiach-1 well area. The current value of the geothermal gradient measured in this well (75-105 • C/km) points to the presence of advective heat transport along basement faults [6,67,68].
In Scenario 1, in order to calibrate the %Ro data and to approximate the temperatures estimated by thermal modelling of AFT data in [6] (Table 2), an erosion of, respectively, 700 and 1100 m, for the Cretaceous and Miocene unconformities, is necessary. The assignment of a heat flow peak in the Early Cretaceous gives a better fit to the measured data ( Figure 6). In Scenario 2 the best calibration was obtained assuming a constant heat flow of 60 mW/m 2 along the entire Mesozoic section. The Early Cretaceous peak gives thermal conditions and temperatures in the Cretaceous that are too high with respect to the measured thermal data. In this scenario, to calibrate the %Ro data and to approximate the temperature estimated by thermal modelling of AFT data in [15] ( Table 2) an erosion of 800 and 300 m for the Cretaceous and Miocene unconformities, respectively, is necessary.
In order to understand the effect of the deposition and erosion of Cretaceous deposits on the thermal history and to compare it with the available calibration data, an additional hypothesis was tested for this well (Figure 6). A unique erosion event in the Miocene was assumed, with the removal of 1700 m of deposits, whereas it was assumed no erosion in the Cretaceous and a constant heat flow of 60 mW/m 2 during that period. Even with that adjustment, this hypothesis also provides a good fit with the measured %Ro data; the temperature attained in this scenario in the Cretaceous and in the Miocene are, respectively, lower and higher than the temperature estimated by thermal modelling of AFT data in both the studies considered [6,15].

Böttstein-1 Well
The Böttstein-1 well is located in the Tabular Jura area (Figure 1b). In this well, only the Mesozoic sequence has been preserved, consisting of Triassic deposits (Figure 3). To calibrate the %Ro data nearly 2000 m of erosion is necessary. Assuming that this erosion occurred entirely during the Miocene-Pliocene Jura Mountain uplift is not likely, as it would not be coherent with other observations in the area, where the Jurassic deposits are overlain by the Eocene unit. This indicates that the erosion of the Cretaceous occurred before the initiation of the foreland basin. Furthermore, a unique erosion event in the Miocene-Pliocene would result in the record of temperature in the Miocene and in the Cretaceous, being higher and lower, respectively, than those estimated by the AFT data in both scenarios. This is similar to what is shown in the Weiach-1 well ( Figure 6).
In Scenario 1, in order to calibrate the %Ro data and to approximate the temperatures estimated by thermal modelling of AFT data in [6] (Table 2), an erosion of 1000 m for both the Cretaceous and Miocene unconformities is necessary. The assignment of a heat flow peak in the Early Cretaceous performs a better fit with the measured data. In Scenario 2 the amount of erosion in the Cretaceous necessary to calibrate %Ro and the temperature estimated by thermal modelling of AFT data in [15], for a rock located at nearly 1000 m of depth (Table 2), varies between 1000 and 1300 m, depending whether it is considered a constant heat flow of 60 mW/m 2 or a heat flow peak of 75 mW/m 2 in the Early Cretaceous. Both options correctly fit the calibration data. In the Miocene, 600 m of erosion are necessary to validate the temperature estimated by thermal modelling of the AFT data ( Table 2).

Riniken-1 Well
The Riniken-1 well is located between the Fold-and-Thrust Belt Jura and the Tabular Jura mountain (Figure 1b). The Quaternary deposits lie directly on the Upper Jurassic (Malm) unit, whereas the entire Cenozoic section is missing (Figure 3). The Riniken-1 well penetrates a very thick Permian sequence (≈1000 m). As opposed to the Weiach-1 well, the crystalline basement in Riniken-1 was not reached. By seismic data interpretation, it is assumed that, below this well, more than 3-4 km of Permo-Carboniferous deposits can be still present [76]. Thus, this area could potentially represent the depocenter of the Paleozoic trough.
For the Riniken-1 well only AFT data from [15] are available. Therefore, to calibrate Scenario 1, the temperature obtained by [6] in the Böttstein-1 and Weiach-1 wells, at nearly 1400 m of depth, were considered. To validate these and the %Ro data, an erosion of 500 and 1100 m, for the Cretaceous and Miocene unconformities, respectively, is necessary. In Scenario 2, to calibrate %Ro data and to approximate the temperature estimated by thermal modelling of AFT data in [15], directly measured in the Riniken-1 well at nearly 1800 m of depth (Table 2), an erosion of 600 and 400 m, for the Cretaceous and Miocene unconformities respectively, is necessary. In Scenario 2 a constant heat flow of 60 mW/m 2 in the Early Cretaceous, provides a better calibration fit than assuming a Cretaceous heat flow peak (Figure 6).

Benken-1 Well
In this well, the Triassic and Jurassic sequence is directly overlain by nearly 300 m of Molasse deposits, similar to the neighbouring Weiach-1 well (Figure 3). However, here in Benken-1 the Mesozoic deposits lie directly on the crystalline basement, and the Permian and Carboniferous sequences are missing. AFT data are available only from [6]. To calibrate scenario 2, temperature data were extrapolated from the AFT data obtained in [15] for the Weiach-1 well.
In Scenario 1, to calibrate %Ro data and to approximate the temperature estimated by thermal modelling of AFT data in [6] (Table 2), an erosion of 900 and 1200 m for the Cretaceous and Miocene unconformities, respectively, is necessary. These values are slightly higher than those estimated by [6] with the same thermal modelling approach and using the same input data. Results similar to what were calculated by this author, can be obtained by increasing the heat flow peak in the Early Cretaceous to 80-100 mW/m 2 , which is not likely. Furthermore, if the heat flow is increased the theoretical %Ro curve does not fit the measured %Ro data satisfactorily. The difference between our results and what was proposed in [6] can be due to the use of a new version of the thermal modelling software, because the embedded calculations have been greatly improved in recent decades.
In Scenario 2, to calibrate %Ro data and to approximate the temperature estimated by AFT data modelling by [15] at nearly 1940 m of depth in the Weiach-1 well (Table 2), 1200 and 500 m of erosion are necessary for the Cretaceous and Miocene unconformities, respectively. A constant heat flow of 60 mW/m 2 results in the best calibration result. Assuming a heat flow peak of 75 mW/m 2 in the Early Cretaceous is not likely, because the temperature computed is too high with respect to that estimated by thermal modelling of AFT data in [15] in the Weiach-1 well at the same depth. Furthermore, the %Ro data are less well calibrated than with a constant heat flow of 60 mW/m 2 ( Figure 6).

Schafisheim-1 Well
The Schafisheim-1 well is located close to the Fold-and-Thrust Belt Jura (Figure 1). The stratigraphic sequence penetrated by the well is similar to that of the Benken-1 and Weiach-1 wells, with nearly 300 m of Molasse deposits on top of the Jurassic unit ( Figure 3). Similar to Benken-1, in Schafisheim-1 the Mesozoic deposits lie directly on the crystalline basement. For Scenario 1 only temperature estimated from AFT data for the Miocene heating event are available [6]. Therefore, to calibrate the Cretaceous heating, temperatures are extrapolated from the data obtained in the Weiach-1 well in [6]. In [15], no data are available for this well. Therefore, for Scenario 2 the temperatures estimated from the thermal modelling of AFT data obtained in [15] in the Riniken-1 and Weiach-1 wells were used.
In Scenario 1, to calibrate %Ro data and the temperature estimated by AFT data modelling ( Table 2) 1400 m of erosion is necessary in the Miocene. For Cretaceous time, it is necessary to approximate the temperature estimated by thermal modelling of AFT data in [6] in the Weiach-1 and Böttstein-1 wells at around 1400-1500 m of depth (Table 2), which implies 500 m of erosion. In Scenario 2 it was necessary to use the data from 1840 m in the Weiach-1 and Riniken-1 wells (Table 2) to calibrate the %Ro data and approximate the temperature estimates based on thermal modelling of the AFT data in [15]. This presumes an erosion of 700 and 500 m, respectively, for the Creataceous and Miocene unconformities. Assigning a heat flow peak of 75 mW/m 2 during the Early Cretaceous corresponding to the lower Cretaceous erosion (500 m) is required to calibrate the data. This option is plausible, as the temperatures estimated by the model over these times are still consistent with the AFT data.

Herdern-1 Well
The Herdern-1 well is located in the southernmost part of the studied area (Figure 1b). A thick Molasse sequence (1292 m) is recorded in this area, lying directly on the Upper Jurassic deposits (Figure 3). In this well the crystalline basement is reached at 2130 m, and, as in the Benken-1 well, it is directly overlain by the Mesozoic sequence. For this well only AFT data for the Miocene heating are available, measured in [6]. Therefore, for Scenario 2, the AFT data measured in the Weiach-1 well was used.
In Scenario 1, to calibrate %Ro data and to approximate the temperature estimated by thermal modelling of AFT data in [6] in the Molasse unit at Miocene time (Table 2), an erosion of 1700 m of is necessary. A similar calibration is performed assigning a Cretaceous erosion of between 300 and 800 m. As for this well we do not dispose of AFT data for the Cretaceous heating, a more precise value cannot be constrained.
In Scenario 2, by considering the temperature estimated by thermal modelling of AFT data in the Weiach-1 well where the crystalline basement is reached at similar depth as in Herdern-1 well (2130-2150 m) the Cretaceous should have attained temperatures in the range of 90-120 • C (Table 2). To reach these temperatures, a Cretaceous erosion of 700-800 m is necessary, depending on the heat flow assumed in the Early Cretaceous (60 and 75 mW/m 2 , respectively).

Discussion
The amount of erosion calculated by thermal modelling based on the %Ro and AFT data for the two scenarios is summarised in Table 3 and illustrated in Figure 7. For the Triassic/Jurassic-Eocene unconformity the erosion estimated in the entire studied area is of a similar order of magnitude for both scenarios, ranging between 500 and 1000 m for Scenario 1 and between 700 and 1300 m for the Scenario 2. The greatest erosion is recorded in the north and it decreases towards the south and south-west (Figure 7). For the Miocene-Quaternary unconformity, the erosion estimate differs remarkably for the two scenarios: between 1100 and 1700 m for Scenario 1 and between 300 and 600 m for Scenario 2. In both scenarios, the erosion amount increases towards the south (Figure 7). The greater erosion values calculated for Scenario 1 are necessary to attain the Miocene heating peak reconstructed by thermal modelling of the AFT data in [6], which is not observed in [15]. Thus, in Scenario 2 a lower magnitude of erosion in the Miocene is necessary. In the following paragraphs, the obtained results are discussed with respect to the geodynamic evolution of the area.

The Jurassic/Cretaceous-Eocene Unconformity
The extent of erosion represented by the Jurassic-Cretaceous to Eocene unconformity has engendered much debate. The hypothesis of very little or no sedimentation in the distal part of the NAFB during the Late Jurassic to Eocene period, has been based on i) the very reduced Lower Cretaceous deposits thicknesses preserved in the Jura area, ii) a Cre-

The Jurassic/Cretaceous-Eocene Unconformity
The extent of erosion represented by the Jurassic-Cretaceous to Eocene unconformity has engendered much debate. The hypothesis of very little or no sedimentation in the distal part of the NAFB during the Late Jurassic to Eocene period, has been based on i) the very reduced Lower Cretaceous deposits thicknesses preserved in the Jura area, ii) a Cretaceous sea level fall and regression phase which could strongly limit the deposition at that time, and iii) the end of the subsidence trend at the Late Jurassic [5,21,26,34]. To support this hypothesis, and in order to calibrate the organic paleothermal data measured in the area, a unique burial event in the Oligocene-Miocene has to be postulated, resulting in the higher thermal maturity indicated by the %Ro data. However, a hypothesis of a unique maximum heating event in the Oligocene-Miocene is contraindicated by apatite fission track data [6,10,14,15] which indicate an important heating event also in the Cretaceous. Our basin thermal modelling results show that the deposition of 500-1000 m of Lower Cretaceous deposits, depending on the area and on the scenario assumed, is the best explanation for the temperatures obtained by thermal history modelling of the AFT data for this interval of time (see calibration results for this scenario in the Weiach-1 well, Figure 6). Furthermore, in the case that no deposition is assumed in the Cretaceous, a greater thickness of eroded material has to be assigned to the Miocene erosion event, which would result in a higher temperature than that estimated by the thermal history modelling of AFT data ( Figure 6, Weiach-1 well). Thus, our results support the Cretaceous heating event registered by AFT data as most likely caused by burial in the Early Cretaceous. The sedimentary sequence deposited at that time was likely eroded in the Late Cretaceous to Paleocene time, when the basin was uplifted and affected by a gradual cooling trend [6,10,14,15].
To justify the accommodation of several hundreds of meters of sediments in the Early Cretaceous there must have been a subsidence acceleration at that time. The subsidence acceleration can be related to extension pulses affecting the European margin during the spreading of the Tethys in the Late Jurassic-Early Cretaceous. Extension tectonics in the Early Cretaceous have been observed in much of western Europe (from the Bay of Biscay through the Aquitanian Basin and Provence into the Western Alps) [77][78][79]. Normal faulting in the Helvetian shelf, Cretaceous in age, has been also reported [65,80], which could be also active northward, along the European margin, where the studied area is located. At the end of the Early Cretaceous, the rift evolved to drift in the Bay of Biscay with the formation of oceanic crust and important thermal anomalies (e.g., [81]). The regional extension affecting the Tethys area could explain the acceleration of the subsidence rate in the European margin. This was necessary to accommodate a thick sedimentary sequence (0.5-1 km) and enabled the rocks to attain the temperature estimated by thermal modelling of AFT data in this area. A heat flow peak could be associated to the lithosphere extension, as proposed by Mazurek et al. [6]. In fact, in Scenario 1, a better fit of the calibration data is generally obtained by introducing a heat flow peak of 70-80 mW/m 2 in the Late Jurassic-Early Cretaceous, higher than the assumed constant heat flow of 60 mW/m 2 . On the other hand, in Scenario 2 in most of the wells %Ro and AFT data are better calibrated by using a constant heat flow in the Early Cretaceous, although in some cases both options are plausible. Therefore, with the data available at present, the heat flow magnitude during the Cretaceous period it is hard to constrain.
To explain the erosion of the entire Cretaceous and part of the Jurassic deposits an exhumation of the area must be envisioned before the deposition of the Oligocene-Miocene foreland sequence. The subaerial emersion of the Mesozoic deposits, since at least the Paleocene-Eocene time, is supported by the findings of Mesozoic karstified structures infilled by younger deposits [79] of continental origin and dated as Eocene (Sidherolitic unit, [82,83]). The aerial exposure and erosion of the Mesozoic sequence can be related to the uplift of the Alpine foreland forebulge. This led to the progressive erosion of older units towards the northwest and to the northward unconformable onlap of the foreland sequence on the Mesozoic passive margin deposits [84][85][86].
An earlier uplift of the Mesozoic sequence starting in the Late Cretaceous is also possible. This is suggested by the cooling trend reconstructed at that time by [15]. In the Late Cretaceous the European plate started to subduct as a consequence of the African plate drift direction change from east to north, with a subduction zone forming between the Piedmont Ocean and the Adriatic continent in the first Cretaceous Alpine orogeny [22,87,88]. This initial subsidence phase could have led to the formation of a first forebulge structure in the European plate, thus impeding the deposition of the Upper Cretaceous record and initiating the erosion of the Lower Cretaceous deposits. Important inversion features in the Cretaceous have been observed in the German Molasse Basin, resulting in northwestward truncation of Cretaceous and Upper Jurassic strata [33,89]. The same authors have described a similar structure in the Swiss part of the Zürich high, formed during a Cretaceous inversion phase [33]. Therefore, an uplift of the basin since that time can be supposed, which would have led to the erosion of 0.5-1 km of Cretaceous and Jurassic deposits. The erosion rate increased northward, where the forebulge uplift was more accentuated.

The Miocene-Quaternary Molasse Erosion
Previous studies have estimated a significant variation of the eroded thickness throughout the NAFB, from an average of 2500 m in the southwestern sector to 350 m in the eastern sector [4,8,9,12]. In the north-eastern sector of the distal part of the NAFB, the amount of erosion estimated by previous authors is variable. Depending on the method used, up to 4 km from rock density and compaction studies [90]; up to 2500 m based on porosity and velocity evaluations in Middle Jurassic rocks [7]; from 0 to 700 m by applying subsidence analysis, clays transformation and extrapolation of coalification trends [4] and references herein); nearly 1000-1100 m by means of thermal modelling and apatite fission tracks [6]; from 500 to 1000 m by recently acquired apatite fission tracks data [15]. For this sector of the NAFB, our thermal modelling results calculates, for the Miocene-Quaternary unconformity, erosion values varying between 1100 and 1700 m for Scenario 1 and between 300 and 600 m for Scenario 2. The high values estimated in Scenario 1 are necessary to calibrate the thermal heating event identified in [6] in the Miocene. This scenario implies that, at that time, the external part of the foreland basin was affected by important subsidence, with the deposition of a sedimentary infill thickness in the same range as the middle-internal part (from 1000 to 1500 m, [12]). Another possibility that could explain the Miocene heating event rather than significant burial is to consider a heat flow peak during the foreland infill time, high enough to attain the paleo-temperature estimated by %Ro and AFT data. This hypothesis has been proposed by [5,10], which supports an anomalous heating trend in Eocene-Miocene time ascribed to deep hydrothermal fluid circulation associated with the Upper Rhine Graben rift and volcanic activity.
Results from more recent AFT studies do not recognise major heating events during the Cenozoic [14,15] which suggests a reduced burial of the foreland basin in this area. This requires only deposition and erosion of a low Molasse unit thickness (300-600 m, Scenario 2). This scenario supposes that the northern distal area of the NAFB was subsiding less than the rest of the basin, a typical trend observed in a foreland basin. Reduced subsidence could have been accentuated by Alpine foreland forebulge uplift, occurring until Miocene time, and/or to interference with the exhumation of the northern crystalline Black Forest Massif taking place during the Late Oligocene to Early Miocene times [91,92]. The erosion is greater toward the south of the studied area (Figure 7), most likely affected by the northward Alpine thrusting in the Miocene to Pliocene period. Due to the good quality of both of the data sets considered, we think the results obtained for both scenarios are reliable. To solve the discrepancies highlighted further work will be needed; specifically increasing the number of samples and wells analysed and/or integrating other thermochronometers such as (U-Th-Sm)/He or clumped isotopes, among others.

Modelling Results Uncertainties
There are uncertainties in the erosion values estimated herein by thermal modelling due to several assumptions about the data input (e.g., timing of the erosional unconformities), on the boundary conditions (e.g., heat flow, paleo-water depth), and potential calibration data errors (e.g., poor constrain on the %Ro and AFT results due to scarce material, uncertainties in AFT data modelling, etc.). Furthermore, software limitations have to be also considered. Our experience demonstrates that the erosion magnitude calculated by 1D thermal modelling generally overestimates what is obtained by 2D and 3D modelling by some hundreds of meters. This is because the 2D and especially the 3D model can take into account the lateral convection heat transfer, which is different from the 1D model where heat transfer can be simulated only by conduction in a vertical direction [93]. This could result in lower values than we have presented here. Thus, further 2D and 3D basin modelling of the area could constrain our results.
An additional uncertainty must be considered for the Miocene-Quaternary erosion. Greater values than proposed in this contribution could be assumed if we consider a sequential repetition of deposition and erosion events, which in a tectonically active foreland basin might be more realistic than a continuous depositional event followed by a unique great erosion event. The alternance of depositional and erosional events would make it difficult to attain significant burial depths, and thus relevant heating, but despite that, cumulatively, a greater sediment thickness would be deposited and eroded, with respect to what is estimated herein. However, because thermal modelling does not simulate very short-term sequential deposition and erosion events this hypothesis cannot be properly tested.

Conclusions
The results presented in this work provide an estimate of erosion during the main unconformities recognised in the north-eastern sector of the distal part of the NAFB. Unlike work we used a thermal modelling approach. Several variables such as the sediment decompaction, thermal conductivity, paleo-water depth, surface sediment-water interface temperature, the heat flow variation, and the timing of burial were integrated. The novelty of our modelling approach is that we did not restrict our modelling calibration to data from irreversible paleothermomenters, such as vitrinite reflectance. Thermochronometers, such as AFT data, were used, giving a timing and a temperature of the main heating and cooling events. By integrating these data with vitrinite reflectance data the amount of erosion corresponding to each unconformity was quantified. Different tectono-stratigraphic evolution scenarios were reconstructed, based on two representative thermochronology studies proposed in the literature.
For the Cretaceous unconformity, the deposition of a sequence of between 0.5 and 1 km thick is necessary to justify the heating event recorded by AFT data in the Early Cretaceous. These deposits would be subsequently eroded in the Cretaceous-Eocene. Our thermal modelling results calculate similar results for two scenarios considered. The option of a heat flow peak during the Early Cretaceous, instead of assuming a constant trend, is the less likely, although it cannot be excluded. The formation of accommodation space at that time could be triggered by repeated extensional pulses affecting the European margin during the spreading of the Tethys in the Late Jurassic-Early Cretaceous. The subsequent erosion of the Lower Cretaceous deposits, as well as part of the uppermost Jurassic sequence, can be related to the uplift of the Alpine foreland forebulge, occurring from the Late Cretaceous to the Eocene time. This would explain the northward increasing erosion trend. The gradual Late Cretaceous cooling, observed in the AFT data, can be related to uplift and thickening of the lithosphere as a consequence of the transition toward a compressional setting.
For the Miocene-Quaternary unconformity, the two scenarios considered provide quite different results. One assumes important deposition and subsequent erosion of thicknesses during the formation and relative uplift of the Alpine foreland basin of 1100-1700 m whereas the second postulates considerably lower values of 300-600 m. This second scenario infers lower subsidence in the distal sector of the foreland with respect to the more internal part most likely limited by the forebulge uplift and/or by the exhumation of the Black Forest Massif.
Supplementary Materials: The following are available online at https://www.mdpi.com/2076-326 3/11/2/62/s1, Annex S1: Vitrinite reflectance data collected from literature. Annex S2: Apatite fission tracks data collected from two previous literature studies [6,15] and description of the acquisition methods used by these authors. The set of data presented by these authors differ in the acquisition method (to obtain the AFT age) and in the thermal modelling procedure used to obtain the timetemperature histories. Mazurek et al. (2006) [6] methods: for most samples 238U was determined via 235U-induced fission tracks resulting from a nuclear thermal irradiation on an external detector (Ufree mica) as a proxy. A couple of samples were determined by the "population method" (multigrain method) mostly used in the 70s and now less employed ("population method" is less precise than the external detector method [94]). Data processed by forward model is obtained using MonteTrax [95]: the forward procedure tests how the analytical data fits with an assumed annealing model (in this case [96]. The modelling provides a single solution (single T-t path). Use of the annealing model of [96]. This annealing model is well known for providing late cooling events, which might lead to an interpretation of a kilometre or more of unroofing. This is usually considered a modelling artefact with no geological significance (see discussion in [97]). Villagomez et al. (2020) [15] method: -238U is directly determined by LA-ICP-MS. Data are processed by inverse model using the HeFTy ® software (version 1.9.3; August 2017). The inverse modelling procedure [98] predicts parameters for various thermal history paths and compares them with the observed analytical data. The controlled random search procedure identifies those thermal histories that most closely match the analytical data. The software provides a huge range of possible solutions. The modelling takes into account annealing variability due to c-axis projection [98]. It used the Dpar (mean maximum diameter of fission track etch figures parallel to the crystallographic c-axis, [99]) as a proxy for the apatite chemistry required for the latest multicompositional annealing algorithms. Funding: This study has been funded by the Swiss Federal Office of Energy (SFOE) and by the Swiss Federal Geological Survey Office (Swisstopo) and is part of the UNCONGEO project carried out by the GE-RGBA team at the University of Geneva (https://www.unige.ch/ge-rgba).

Data Availability Statement:
The study did not report any data. All the data used are in the manuscript and in the annexes.