## Abstract

Interpretation of low-temperature thermochronological data usually relies on assumptions on the shape of isotherms. Recently, a number of thermal modelling approaches investigate and predict the theoretical influence of topography on isotherms. The application and proof of these predictions is not well confirmed by measured data. Here we present apatite fission-track (AFT) data from samples collected along the Gotthard road tunnel and its corresponding surface line to test these predictions. AFT ages broadly cluster around 6 Ma along the tunnel. No correlation of tunnel ages with superimposed topography is seen, which means that topography-induced perturbation of isotherms under given boundary conditions (topographic wavelength 12 km; relief 1.5 km; exhumation rate 0.45 km Ma^{−1}) can be neglected for the interpretation of AFT ages. Thus, in areas characterized by similar topographies and exhumation rates, apparent exhumation rates deduced from the age–elevation relationship (AER) of AFT data need no correction for topography-induced perturbation of isotherms. Three-dimensional (3D) numerical thermal modelling was carried out incorporating thermally relevant parameters and mechanisms, such as topography, geology, thermal conductivities and heat production. Modelling reveals a strong influence on the shape of isotherms caused by spatially variable thermal parameters, especially heat production rates. Therefore, not only topography has to be considered for interpreting low-temperature thermochronological data, but also other parameters like heat production rates.

**Supplementary material:** 1. Electron microprobe analyses, 2. Topography and model extend, 3. Model parameters are all available online at http://www.geolsoc.org.uk/SUP18380.

The first geological applications of low-temperature thermochronology commonly assumed that critical isotherms (e.g. 110 °C for apatite fission track) remain horizontal in relation to topography (Schaer *et al.* 1975). Thermal modelling, however, predicts that in areas with pronounced relief near-surface isotherms will be influenced by topography with compressed isotherms beneath valleys and wider-spaced isotherms beneath ridges (e.g. Stüwe *et al.* 1994). The impact for the interpretation of low-temperature thermochronological data, especially apatite fission track (AFT) and (U–Th)/He data, was investigated recently (Stüwe *et al.* 1994; Mancktelow & Grasemann 1997; House *et al.* 1998, 2001; Stüwe & Hintermüller 2000; Braun 2002; Foeken *et al.* 2007). On the other hand, bending of isotherms carries potential information to reveal minimum ages of topography and the evolution of palaeorelief (House *et al.* 1998; House *et al.* 2001; Braun 2002; Foeken *et al.* 2007).

A widely used technique to derive exhumation rates utilizes the correlation of ages of a single isotopic system, e.g. AFT or apatite (U–Th)/He, with elevation, the so-called age–elevation relationship (AER) (Wagner & Reimer 1972; Schaer *et al.* 1975; Stüwe *et al.* 1994). The great advantage, in contrast to the mineral-pair method (Wagner *et al.* 1977), is that no estimation of the geothermal gradient has to be made. However, this proxy is only valid under certain assumptions (e.g. Parrish 1983): (1) during and after passing the closure isotherm all samples followed a vertical exhumation pathway and no tectonic displacement exists between samples; (2) all samples are kinetically uniform; and (3) the closure isotherm is fixed with respect to some geographical reference horizon, e.g. sea level.

A drawback of many studies is that sample profiles are more horizontal than vertical, and thus deduced exhumation rates from such profiles are affected by topography-induced perturbation of isotherms. The magnitude of isotherm perturbation increases with topographic wavelength, relief and exhumation rate (Stüwe *et al.* 1994). Hence, in mountainous regions the assumption that AFT relevant isotherms are flat is probably not fulfilled. In this case only vertical profiles (along a borehole or steep cliff) give correct values of exhumation rates. Resulting apparent exhumation rates of AER plots from ‘near’-vertical profiles have to be corrected for the effect of isotherm perturbation due to topography. This leads to the question of the geomorphic and structural–kinematic boundary conditions that cause isotherm perturbation. Based on numerical modelling and analytical solutions (Stüwe *et al.* 1994; Mancktelow & Grasemann 1997; Stüwe & Hintermüller 2000; Braun 2002), exhumation rates of ≥0.5 km Ma^{−1} along with a topographic wavelength of ≥10 km and a relief of ≥2 km are threshold values for perturbing the 110 °C isotherm.

The aim of this study is to estimate the perturbation of isotherms under a given framework of exhumation rate, topographic wavelength, relief amplitude, rock properties and geothermal gradient by applying AFT thermochronology. AFT age variations along the tunnel transect give evidence as to whether the 110 °C isotherm was perturbed and, if so, to what extent. For this we sampled the Gotthard road tunnel as well as the corresponding surface line directly above the tunnel. The Gotthard road tunnel is located in central Switzerland and has a length of 16.3 km (Fig. 1). Present-day uplift rates range between 0.7 mm year^{−1} in the north and 1.1 mm year^{−1} in the south of the tunnel (Kahle 1997) (Figs 1b & 2a). The Gotthard transect is characterized by more or less ENE–WSW-trending ridges and valleys with a topographic wavelength of 12 km and relief of up to 1.5 km (Fig. 1b). Published exhumation rates constrained by AFT data from the Gotthard Massif range between 0.45 and 0.5 km Ma^{−1} between 10 and 6 Ma (Schaer *et al.* 1975; Wagner *et al.* 1977; Michalski & Soom 1990). These values suggest that the Gotthard transect is close to the lower limit of the proposed parameter values given for a 110 °C isotherm perturbation. Thus, this study provides a natural benchmark to verify topography-induced perturbation of isotherms predicted by existing modelling approaches.

In addition, a new three-dimensional (3D) finite-difference thermal model for predicting the positions of isotherms beneath topography is presented. Compared to previously published analytical solutions (Stüwe *et al.* 1994; Mancktelow & Grasemann 1997; Stüwe & Hintermüller 2000) our model provides the following advantages: (1) it incorporates 3D topography; and (2) it includes temporally and spatially variable parameters (e.g. thermal conductivity, internal heat production, exhumation rate and topography). Thus, the model provides information about the influence of different thermal parameters on predictions of the shape of near-surface isotherms, comparable to the program ‘Pecube’ of Braun (2003).

## Geological setting

The Gotthard road tunnel cross-cuts the entire Gotthard Massif (GM) and the southern part of the Aar Massif (AM). The GM forms an ENE–WSW-trending mountain range 80 km long and up to 12 km wide (Fig. 1a). The GM and AM belong to the external Massifs of the Alps. They consist of pre-Variscan polymetamorphic basement intruded by late Variscan granitoids and covered by Late Palaeozoic–Mesozoic sedimentary rocks (Fig. 1) (e.g. Labhart 1977, 1999; Schaltegger 1994). The boundary between GM and AM is marked by the heavily tectonized Urseren–Garvera zone, built up by steeply dipping Permo-Carboniferous and Mesozoic metasediments. The sedimentary cover of the GM and AM was mostly detached during the Alpine orogeny, forming parts of the Helvetic nappes.

Peak Alpine metamorphic conditions along the tunnel were reached at around 35–30 Ma, with greenschist-facies conditions in the north and amphibolite-facies conditions in the south (Frey & Ferreiro Mählmann 1999). The post-metamorphic evolution of the study area was mainly controlled by the northward movement of the Adriatic indenter and related thrusting in the external Alps (e.g. Schmid *et al.* 1996). Thrusting during the ‘Grindelwald stage’ (22–12 Ma) propagated towards the foreland (Schmid *et al.* 1996), resulting in thrusting of the Gotthard Massif upon the Tavetsch and Aar Massif around 20 Ma (ductile deformation) and subsequent steepening of these structures (brittle deformation) (Wyder & Mullis 1998). It is likely that brittle deformation continues up to the present, as indicated by post-glacial vertical movements that are common in the study area (Persaud & Pfiffner 2004; Ustaszewski *et al.* 2008).

## Methods

Sample sets along the Gotthard transect were collected and 31 AFT ages measured (Fig. 2a and Table 1). Minerals were separated using standard magnetic and heavy liquid techniques. Apatites were mounted in epoxy, and their internal surfaces were excavated by grinding and polishing. Afterwards the mounts were etched by 5 M HNO_{3} for 20 s at 20 °C. Irradiation was carried out at the FRM-II reactor in Garching (TU München, Germany). Mica detectors were etched to reveal induced tracks using 40% HF at 20 °C for 40 min. Fission-track counting was carried out with an optical microscope (Axioscope 1, Zeiss) under 1000× magnification using dry lenses. Samples were dated using the external detector method (Naeser 1978; Gleadow 1981) using the zeta calibration approach (Hurford & Green 1982, 1983), with a zeta of 354±7 years cm^{−2} for dosimeter glass CN5 and Durango and Fish Canyon Tuff apatite age standards. Age determination, visualization and statistics were calculated and performed with Trackkey 4.2 g (Dunkl 2000). All AFT ages are displayed as central ages and errors as ±1*σ* (Galbraith & Laslett 1993).

Measured *D*_{par} values were used as kinetic parameters (e.g. Burtner *et al.* 1994), complemented by electron microprobe analysis on selected crystals, using a JEOL Superprobe with a beam current of 30 nA, an acceleration voltage of 15 kV and a beam diameter of 10 µm.

## Inferred shape of isotherms from AFT data

The aim of this section is to investigate the perturbation of measured AFT ages under the given boundary conditions, especially with respect to topography.

Table 1 summarizes the results of AFT dating along the tunnel and the corresponding surface line. Previously published AFT data for the Gotthard area (Schaer *et al.* 1975; Wagner *et al.* 1977; Michalski & Soom 1990) are in good agreement with our data and are therefore included in our interpretations (cf. Fig. 2b).

To make sure that measured variations of AFT ages are, in fact, the result of the shape of the palaeo-isotherm, the kinetic homogeneity of our apatites was tested by etch-pit diameter (*D*_{par}) measurements. *D*_{par} values obtained for measured samples are small and relatively uniform, varying between 1.1 and 1.6 µm (Table 1). Microprobe analyses of 15 samples from different lithologies were carried out, demonstrating that all samples are close to the F-apatite end member, with Cl contents of less than 0.1 wt%. Furthermore, analysed elements such as Si, Mn, Ce and Sr do not show any significant variation, suggesting that apatites are kinetically uniform (for details see p. 1 of the supplementary material for this paper).

AFT ages of samples collected at the same elevation can be used to predict the palaeoshape of the critical isotherm (House *et al.* 1998). All samples along the Gotthard tunnel are at approximately 1.1 km elevation and yield AFT ages of 6.2±0.6 Ma, overlapping within 1σ-error limits (Fig. 2b). Individual samples characterized by low U-content show larger 1σ-error, especially those from the northern and southern ends of the transect.

Surface samples yield AFT ages ranging from 7.0 to 9.6 Ma. Surface samples collected at approximately 1.3 km elevation along the Gotthard transect (CGP 07, MRP 276, MRP 278 and MRP 279) exhibit the same ages of 7.4±0.2 Ma. Thus, AFT age patterns on two structural elevation levels (*c*. 1.1 and *c*. 1.3 km) indicate that between 6.2 and 7.4 Ma the 110 °C isotherm was flat, suggesting that topography has no visible effect on the AFT age pattern along the tunnel. Small variations in ages partly coincide with stratigraphic–tectonic boundaries (Fig. 2b), which can be explained by kinetic or physical differences (such as heat production and thermal conductivities).

Surface AFT ages show an increase with elevation (Figs 2b & 3). Slopes of AER were calculated by weighting the ages according to their errors. Slopes and uncertainties were inverted to estimate exhumation rates with defined uncertainties (Reiners *et al.* 2003).

Regression bands were plotted with a 95% confidence interval, which allows prediction of the error of the regression lines. The AER plot for samples of the northern flank of the Gotthard Massif (MRP 237, MRP 238, CGP 11 and MRP 292) yields an apparent exhumation rate of 0.35±0.1 km Ma^{−1} for the period between 10 and 5 Ma (Fig. 3a). The horizontal and vertical distance between the samples is 1.5 and 1.2 km, respectively. Samples of the southern flank (MRP 229, MRP 276, MRP 294, CGP 07 and A14 of Schaer *et al.* 1975) yield an apparent exhumation rate of 0.5±0.2 km Ma^{−1} for the period between 10 and 6 Ma (Fig. 3b). Here the horizontal and vertical distance between the samples is 4 and 1.6 km, respectively. Deduced apparent exhumation rates are in good accordance with previously published data (Schaer *et al.* 1975; Michalski & Soom 1990). The presented data suggest that critical isotherms are flat under the Gotthard transect and a correction for topography is not necessary.

## Modelling the shape of isotherms

Although we demonstrated that in the specific case of the Gotthard transect no topographic correction is required, in the following section we test the theoretical influence of topography and other thermally relevant parameters on near-surface isotherms. We present numerical 3D finite-difference thermal models developed for a crustal block of 14×26×10 km around the Gotthard road tunnel (an overview is provided on p. 2 of the supplementary material). Modelling incorporates spatial and/or temporal varying parameters, namely topography, exhumation rates, thermal conductivity and heat production (for details see p. 3 of the supplementary material). The parameter values are either taken from this study (exhumation rates) or from geophysical measurements (e.g. Busslinger & Rybach 1999). In the following, all mentioned isotherm perturbation values refer to the location of the Gotthard road tunnel transect, and are defined here as the maximum vertical deflection of the corresponding isotherm.

For simplification in the first step it is assumed that heat transfer in the crust is controlled by conduction only (Stüwe 2007). In this case Fourier's law of heat conduction can be used:
1
where κ is the thermal diffusivity with typical values of around 10^{−6} m^{2} s^{−1}, *T* is temperature, *t* is time and ∇ is the nabla operator, describing here the spatial change of temperature in 3D (e.g. Stüwe 2007). Equation (1) is modified to account for the effect of exhumation and internal heat production due to radioactive decay, resulting in:
2
where *z* is depth, *S* the heat production rate, ρ the rock density and *c*_{p} the specific heat. Initial and boundary conditions necessary for the solution of the differential equation (2) are as follows.

The lower boundary of the model was placed at −10 km from the surface and was fixed with a constant vertical geothermal gradient, for which a value of 20 °C km

^{−1}was used based on borehole and tunnel temperature measurements in or near the Gotthard road tunnel (Keller*et al.*1987). To verify that this relatively ‘shallow’ lower boundary allows the deflection of isotherms caused by thermal relevant parameters (e.g. topography and exhumation), a simple 2D sinus-shaped topography was modelled analytically using the program ‘TERRA’ (Ehlers*et al.*2005) (Fig. 4). Used parameters fit those from the Gotthard transect. Position of isotherms are unaffected by the depth of the lower boundary of the model.Topography was used as an upper-boundary condition, discretized with a rectangular grid and converted into air temperature, depending on elevation (Busslinger & Rybach 1999): 3 with

*T*_{Air}being the mean air temperature at surface,*T*_{0}the air temperature at sea level, α the atmospheric lapse rate and*h*the elevation. Kohl*et al.*(2001) showed that ground surface temperature (GST) data collected from high elevations in Alpine terrain are best described by the model of Niethammer (1910), with*T*_{0}=12 °C and α=4.6 °C km^{−1}. These values were used to calculate GST for each node based on a Shuttle Radar Topography Mission (SRTM) digital elevation model (DEM). We neglected the difference between air temperature and the corresponding GST depending on exposition (slope and orientation), vegetation, snow cover and rock surface properties (e.g. Šafanda 1999). Consequently, each node of the upper boundary is composed of an elevation and a corresponding GST.Lateral boundaries were fixed with no horizontal heat flow. For simplification a uniform distance from node to node (Δ

*x*=Δ*y*=Δ*z*) of 132.9 m was used, according to the spatial resolution of the input DEM. Thus, the total mesh was built up by 1.55×10^{6}nodes. A simplified geological model was used for the modelling, which considered five geological units: Aar Massif, Urseren Zone, Tavetsch Massif, Gotthard Massif and Schistes lustrés (Fig. 1 and Table 2). Vertical boundaries between individual geological units were assumed, consistent with high dip angles observed at the surface and during tunnel construction (Fig. 2a). The geological units are characterized by different physical and thermal properties. The foliation developed during Alpine metamorphism results in a pronounced anisotropic thermal conductivity of all geological units, with ratios of up to 1.5 between maximum (*k*_{‖}) and minimum (*k*_{⊥}) values within one unit (Kappelmeyer & Haenel 1974; Kohl*et al.*2001). Input thermal conductivities are listed in Table 2. Rock radioactivity measurements within the tunnel (Keller*et al.*1987) were used to calculate the total heat production related to the radioactive decay of uranium, thorium and potassium (Turcotte & Schubert 1982, equations 4–6). Resulting heat production values ranged from 7.54×10^{−11}to 2.56×10^{−9}W kg^{−1}for individual geological units, with 1.22×10^{−9}W kg^{−1}as an average value. Modelling was carried out either with individual values or an averaged heat production rate extrapolated to the surroundings of the Gotthard road tunnel (Table 2). Tunnelling and geophysical investigations revealed that the Aar and Gotthard Massifs continue towards depth, possibly by more than 10 km (Pfiffner & Heitzmann 1997). Therefore, heat production rates are assumed to be constant with depth. Modelling was conducted assuming both a steady-state and time-varying topography, and a constant exhumation rate of 0.45 km Ma^{−1}(mean value of the AER plots in Fig. 3). To account for a time-varying topography nodes were added/removed and temperature adjusted (according to equation 3) at the upper boundary during the modelling progress. For instance, a change in elevation from 0 to 2000 m in 2 Ma, is implemented by adding one node (132.9 m) every 0.133 Ma (reaching an elevation of 1994 m after 2 Ma). A simple 2D thermal model was numerically calculated and compared with an analytical solution calculated with the program ‘TERRA’ (Ehlers*et al.*2005) (Fig. 5), proving the correctness of the used finite-difference thermal model approach.

From several performed model runs we present five solutions, which differ by the assumed thermal conductivity, heat production and topographic evolution (Table 3). Models 1–3 were calculated until reaching a steady-state temperature field, which takes around 5 Ma depending on the parameterization. Models 4 and 5 were calculated with an initial steady-state temperature field, subsequently adapted during the modelling progress because of a changing topography.

To evaluate the influence of thermally relevant parameters (e.g. topography) on low-temperature thermochronology in the following models, closure temperatures of 110 and 60 °C were assumed for the apatite fission-track (e.g. Rahn & Grasemann 1999) and (U–Th)/He system (e.g. Ehlers & Farley 2003), respectively. AFT age differences along the tunnel are predicted, using forward modelling of AFT ages with the program HeFTy (Ketcham *et al.* 2007) based on measured kinetic parameters (*D*_{par}) and generated time–temperature (*tT*) path for selected localities. Therefore, during modelling progress the temperatures of selected localities were read out every 10 ka, yielding nearly continuous *tT* paths.

Figure 6 shows the predicted steady-state 110 °C isotherm for different model runs, and the contours of the corresponding 2D shape of the 60 and 110 °C isotherms along the Gotthard road tunnel.

The benefit of model 1 is that we can estimate the net influence of the topography on a chosen isotherm because all other input parameters are spatially uniform. The modelled 110 °C isotherm clearly demonstrates that small topographic features have no influence on the isotherm (Stüwe *et al.* 1994; Mancktelow & Grasemann 1997). Thus, according to model 1, the shape of the isotherm follows the topography in a strongly dampened fashion (Stüwe *et al.* 1994) and is clearly deformed beneath the Gotthard Massif, with a perturbation of 250 m. Resulting maximum AFT age difference along the tunnel is 0.5 Ma (modelled ages range from 6 to 5.5 Ma).

Incorporating spatially variable and anisotropic thermal conductivities (model 2) results in an overall deeper 110 °C isotherm, because average input thermal conductivities are higher than those used for model 1. The shape of the 110 °C isotherm in model 2 is similar to that from model 1, apart from a smaller isotherm perturbation of only 150 m (Fig. 6). Areas characterized by low thermal conductivities (Urseren zone and Schistes lustrés, cf. Fig. 1) accumulate heat that bulge the 110 °C isotherm. We conclude that the observed differences and anisotropy of thermal conductivities along the Gotthard road tunnel (Table 2) do only weakly influence the shape of the near-surface isotherms. For models 1 and 2, the corresponding perturbation of the 60 °C isotherm is more or less twice as high compared to the perturbation of the 110 °C isotherm. However, the resulting AFT age difference along the tunnel is only 0.2 Ma (modelled ages range from 6.7 to 6.5 Ma).

From model 3, a laterally variable heat production was added, resulting in a more complex 110 °C isotherm shape. Interestingly, the pronounced topography of the Gotthard Massif, forming an ENE–WSW-trending ridge in the southern tunnel section, is not mirrored by the 110 °C isotherm. Evidently, differences in heat production values are capable of blurring the effect of topography. As demonstrated by model run 1, the net topographic effect is small. Thus, for model 3, the perturbation of the 110 °C can be nearly completely attributed to the spatial variation of heat production, which here was assumed to persist to the depth of 10 km. Owing to the lateral decrease in total heat production rate, the shape of the 110 °C isotherm along the tunnel profile shows a steady decrease to the south, with a perturbation of 550 m. The 60 °C isotherm runs more or less horizontally up to tunnel kilometre 12 and then decreases towards the south with a total perturbation of 450 m. For this model run, variations of AFT ages in the range of 1.2 Ma (modelled ages range from 7.6 to 6.4 Ma) along the tunnel are predicted. The perturbation effect of topography and heat production at that depth becomes less pronounced in the northern part of the transect. Thus, differences in heat production are capable of compensating (northern section) and amplifying (southern section) the perturbation of isotherms owing to topography (Fig. 6).

A palaeotopographic evolution with increasing relief after 5 Ma, assuming isotropic thermal parameters, would result in flat isotherms at the time the AFT samples cross the 110 °C isotherm, following the observation of model 1. In the extreme case that the palaeotopography is flat at that time, resulting AFT ages along the tunnel are all the same. Incorporating spatial variable heat production rates (model 4) results in increased isotherm perturbation with depth, with perturbation values of the 60 and 110 °C isotherms of 250 and 500 m, respectively (Fig. 6).

To investigate the effect of a possible palaeotopographic evolution with higher relief and peak elevations in Late Miocene times (6 Ma), which decrease to the present, model 5 was carried out assuming initial maximum elevations of 5 km and a relief of 3 km. The shape of resulting isotherms is comparable to that of model 1, however, with more pronounced isotherm perturbations of 450 m (110 °C) and 750 m (60 °C). The doubling of the vertical isotherm perturbation with respect to model 1 corresponds to a doubling in the FT age difference along the tunnel.

## Influence of advective heat transport on low-temperature isotherms

To investigate the impact of heat advection due to faulting on isotherms and AFT ages simple 2D numerical modelling was carried out. A flat topography was assumed with two crustal blocks exhumed with rates of 0.5 and 1 km Ma^{−1} (Fig. 7a). Assuming steady-state conditions and an exhumation rate of 0.5 km Ma^{−1}, the 110 °C isotherm is at a depth of *c.* 3 km. Initiation of a perpendicular fault and increased exhumation (1 km Ma^{−1}) on one side of the fault results in a perturbation of isotherms depending on the duration of displacement. Short-term movements (<0.1 Ma) result in a less pronounced perturbation of the isotherm across the fault (<50 m) and a displacement of the crustal blocks of <50 m. For small displacements heat advection approximates relative rock uplift (Fig. 7b); hence advective heat flow dominates the temperature field. For long-term movements, however, the role of diffusion increases. Thus the perturbation of the 110 °C isotherm is small as compared to the displacement (e.g. for a duration of 1 Ma, a perturbation of 200 m corresponds to a relative displacement of 500 m). Infinite displacement results in a maximum possible isotherm perturbation of *c*. 320 m.

## Discussion

In the following some essential requirements for the investigation of isotherm perturbation due to topography are discussed. In addition, we discuss the measured AFT age pattern in light of the modelled predictions of isotherm perturbation.

### Advective and convective heat transport

In addition to the discussed thermally relevant parameters (e.g. heat production), near-surface temperatures are affected by advective and convective heat transport. Therefore, an important prerequisite for investigations on isotherm perturbations is that advective and convective heat transport must be negligible or quantifiable over geological timescales.

Advective heat transport may be caused by deformation and spatial differences in rock uplift rates. Studies on brittle deformation structures and hydrothermal mineral precipitation along the Gotthard road tunnel showed that the main phase of brittle faulting had ceased before the structural level of the tunnel samples cooled below 190 °C (Luetzenkirchen 2002). First-order precise levelling of the Swiss national levelling network revealed only a slight gradient of present-day vertical movements within the area of the Gotthard Massif (Fig. 1b), with rates of 0.8–1.0 mm year^{−1} (Kahle 1997). On a smaller scale, however, remote sensing analysis, field work and numerical modelling demonstrate the activity of uphill-facing scarps in the Urseren valley as a response to isostatic vertical movements after the last Ice Age (Dahinden 2001; Persaud & Pfiffner 2004; Hampel & Hetzel 2006; Ustaszewski *et al.* 2008). AFT ages along the Gotthard road tunnel show no jumps in ages (Fig. 2b), thus a differential vertical displacement of individual crustal blocks has been small and is hidden within the 1σ-error of the AFT ages. Assuming an exhumation rate of 0.5 km Ma^{−1} and an average AFT error of 0.6 Ma, only a vertical displacement of more than 300 m along a fault would be significant. We conclude that vertical displacements along fault zones were small and therefore can be neglected for the last 6 Ma. Simple 2D thermal modelling of movements along a single vertical fault show that a displacement of 300 m is capable of generating a perturbation of the 110 °C isotherm of approximately 150 m (Fig. 7). We cannot exclude the fact that mass advection along the Gotthard transect possibly led to a perturbation of isotherms of less than 150 m, which in magnitude is close to the net topography-induced isotherm perturbation (see models 1 and 2, Fig. 6).

Convective heat transport may play an important role in the total heat transport in high mountainous areas (e.g. Whipp & Ehlers 2007). During construction of the Gotthard road tunnel water inflow measurements were carried out, with maximum *in situ* measured rock temperatures of 32 °C (Fig. 2a) (Keller *et al.* 1987; Luetzenkirchen 2002). Except for an area in the central Gamsboden granite-gneiss, which is assumed to be affected by a near-surface convective hydrothermal circuit (Pastorelli *et al.* 2001; Luetzenkirchen 2002), measured rock temperatures are correlated to topography.

We conclude that advective and convective heat transport are negligible, and that conductive heat transport dominates and controls the temperature distribution in the upper crust of the study area.

### Predicted v. observed ages

A striking characteristic when comparing the modelled shape of isotherms and measured AFT ages along the Gotthard transect is the incompatibility of the more complex model 3, including spatially variable thermal conductivities and heat production rates. A possible explanation could be that measured radioactivities in the rocks of the Gotthard road tunnel and the resulting heat production rates are not representative of the entire lithological unit, and, thus, cannot be extrapolated straightforwardly. Alternatively, perturbation of isotherms may be compensated for because of the spatially variable exhumation rates increasing to the south, which is likely when the deduced exhumation rates are considered. Mean exhumation rates of between 10 and 6 Ma are higher on the southern flank of the Gotthard Massif compared to the northern flank (Fig. 3). The present-day uplift rates show the same trend (Kahle 1997). Repeating model run 3 with exhumation rates increasing from north (0.45 km Ma^{−1}) to south (0.55 km Ma^{−1}) results in an AFT age pattern similar to that measured for the surface and tunnel samples (Fig. 8). These exhumation rates are optimized to match the measured AFT ages and the AER-derived exhumation rates. Thus, small spatial differences in exhumation rates are able to compensate for and mask the perturbation of isotherms resulting from other thermal parameters. Therefore, we conclude that measured AFT age variations along the Gotthard road tunnel (6.2±0.6 Ma) do not preclude any presented model or palaeotopographic evolution. The Gotthard transect, however, is characterized by steady exhumation rates (Schaer *et al.* 1975; Wagner *et al.* 1977), and more or less vertically oriented tectonic structures and boundaries (Fig. 1b), facilitating the generation of a steady-state topography with spatially invariant ridges and valleys (Willett & Brandon 2002). But steady-state conditions are rarely achieved, especially relating to the Late Neogene climate change and the increased frequency in climatic oscillations (Whipple & Meade 2006). Most probably, positive feedback between climate change, weathering, erosion and isostatic rebound led to an increase in relief over the last million years (e.g. England & Molnar 1990; Molnar & England 1990). However, resolving the latest Neogene (<5 Ma) exhumation history of the Gotthard Massif is beyond the scope of this paper and will require the use of apatite (U–Th)/He data.

## Conclusions

Apatite fission-track data along the Gotthard road tunnel and its corresponding surface line are compared with 3D numerical thermal modelling results. Modelling allows the interplay and relative importance of 3D topography, spatially variable thermal conductivities, heat production rates and exhumation rates on the shape of low-temperature isotherms to be investigated. In addition, the possible impact of advective and convective heat transport was examined.

The following conclusions were obtained.

The AER of AFT ages yielded exhumation rates of 0.35±0.1 km Ma

^{−1}(northern flank of the Gotthard transect) and 0.5±0.2 km Ma^{−1}(southern flank), similar to previously published data (Schaer*et al.*1975; Michalski & Soom 1990).Measured AFT ages from tunnel samples were around 6.2 Ma (5.6–6.8 Ma), with no significant trend. A perturbation of AFT ages due to topography, with young ages correlating with maximum overburden, is not visible along the Gotthard road tunnel. Because AFT age errors range from 5 to 25% (0.3–1.6 Ma), a perturbation of approximately 0.6 Ma will mostly remain hidden behind the data scatter.

Modelling reveals that spatially variable heat production and exhumation rates strongly influence the shape of near-surface isotherms, and need to be considered for the interpretation of low-temperature thermochronological data. In the specific case of the Gotthard transect, the modelled perturbation of the 110 °C isotherm is in the range of 150–550 m (the corresponding difference in AFT ages along the tunnel is between 0.3 and 1.2 Ma), depending on the input parameters and boundary conditions.

This study illustrates that topography-induced perturbation of isotherms can be neglected for the interpretation of AFT ages under the given boundary conditions (topographic wavelength 12 km; relief 1.5 km; exhumation rate 0.45 km Ma

^{−1}) and petrophysical parameters of the Gotthard transect, such as thermal conductivities.

## Acknowledgments

This study was funded by the German Science Foundation (DFG), project SP 673/2. G. Höckh, D. Mühlbayer-Renner and D. Kost (Universität Tübingen) are gratefully acknowledged for the mineral separation. Thanks also to T. Wenzel (Universität Tübingen) for his support during electron microprobe analysis. M. Zattin and P. van der Beek are thanked for their constructive reviews of this manuscript.

- © Geological Society of London 2009