Gold Open Access: This paper is published under the terms of the CC-BY license.


Translating amounts and rates of rock cooling derived from low-temperature thermochronometry into denudation requires assumptions about the local geothermal gradient. The temperature gradient in the crust depends on many factors, including basal heat flow, crustal heat production, and thermal conductivity. Consequently, geothermal gradients may be variable on time scales over which rock cooling is tracked by thermochronometry. Using one-dimensional numerical modeling of heat transfer in rocks of varying thermal characteristics, we show that the geothermal gradient of the eroded layer is the most important factor for accurate estimation of denudation amounts. Using a three-dimensional numerical model (Pecube), we demonstrate the impact of crustal heat production and thermal conductivity on estimates of total denudation derived from apatite fission track data from central west Britain. We show that the regional variation in cooling ages measured in Caledonian granites can be explained by geothermal gradient variation due to the presence of a heat-producing granite batholith and removal of insulating sedimentary rocks, and does not require variable denudation. Neglecting the blanketing effect leads to twofold overestimation of the amount of denudation. The occurrence of heat-producing basement that was once covered by a sedimentary blanket is common, in particular in the core of mountain belts. Accurate determination of the amount and rate of denudation from thermochronometric studies in these situations must take into account the composition of the eroded rocks.


Low-temperature thermochronometry (LTT) is used to quantify amounts of denudation by converting paleotemperatures into paleodepths. In the absence of a clear indication of elevated basal heat flow, present-day geothermal gradients (20–30 °C/km; Turcotte and Schubert, 2002) are typically assumed. This assumption ignores the spatial and temporal variation in heat production and thermal conductivity of rocks, as well as short-lived heat-flow perturbations, such as from faulting, fluid circulation, magmatism, and rapid rock uplift (Turcotte and Schubert, 2002). Although efforts have been made to predict the effect of these factors on the upper crustal geothermal gradient (e.g., Brown et al., 1994), only recently have studies considered the role of the thermal conductivity of the rocks that are being denuded (e.g., Barbarand et al., 2013; Braun et al. 2016), while variation in crustal heat production is largely neglected.

Surface heat flow is the sum of basal heating of the lithosphere and radiogenic heat production in the crust. Typically, each contributes about equally to surface heat flow (Turcotte and Schubert, 2002). However, where the upper crust comprises high-heat-producing (>3 μW/m3) rocks, crustal heat production may contribute in excess of 80% (Neumann et al., 2000). The thermal conductivity in the uppermost crust varies from 0.5 to 7.0 W/m/K depending on temperature, pressure, porosity, quartz content, and water saturation (Eppelbaum et al., 2014). Low-conductivity rocks such as mudstone, coal, and basalt have an insulating effect on the underlying rocks (often referred to as the “blanketing effect”) that can enhance the maturity of organic material in sedimentary basins (Pollack and Cercone, 1994) and the geothermal energy potential of heat-producing basement (Majorowicz and Minea, 2012).

Acid-intermediate basement rocks are commonly sampled for LTT due to the abundance of apatite and zircon. Exhumation events commonly involve erosion of overlying sedimentary or volcanic successions. In such cases, the combined effect of the low-conductivity rocks overlying heat-producing igneous bodies on denudation estimates derived from LTT data has not yet been investigated. Here we use one-dimensional (1-D) and 3-D (Pecube) numerical models to determine the extent to which the thermal conductivity of eroded strata and heat production in the basement affect the geothermal gradient in the uppermost crust. The modeling results are tested in a case study from central west Britain, where we show that denudation can be overestimated by a factor of two if variations in the thermal structure of the crust are ignored.


Heat transfer in eroding crust follows a second-order equation: 
where T is temperature, t is time, v is velocity of rock uplift, x, y, and z are coordinates of the rock particle, κ is thermal diffusivity, and A is rate of heat production (κ = k/ρ/c and A = H/ρ/c, where k is thermal conductivity, H is heat production, ρ is density, and c is specific heat capacity). In Figure 1 we show solutions to the 1-D (vertical) form of this equation for scenarios designed to determine how the geothermal gradient changes in response to heat production in a pluton and to the thermal conductivity of an overlying sediment pile.

Figure 1A shows how geothermal gradients change in response to pluton thickness and heat production (constant thermal conductivity). For a 4–12-km-thick intrusive body, the geothermal gradient increases only by 1.2–4.4 °C/km for each 1 μW/m3 increase in H.

For a “normal” crust, where heat production is low and constant (1 μW/m3), the geothermal gradient increases rapidly when the thermal conductivity of the blanketing layer is <2 W/m/K (Fig. 1B); high-porosity–low-conductivity rocks, such as chalk (≤1.5 W/m/K; Eppelbaum et al., 2014), can produce geothermal gradients of >40 °C/km.

Figure 1C shows the combined effect on the geothermal gradient of varying radiogenic heat production and thermal conductivity. The two factors can either produce extremely low or high geothermal gradients, or almost cancel each other, leaving the gradient practically unchanged. The point to note here is that geothermal gradients in excess of 50 °C/km can readily be generated in a low-conductivity sedimentary layer overlying a high-heat-producing granite.

Figure 1D shows the thermal structure of the crust for varying heat production and with or without a sedimentary cover. For 2 km of erosion, calculating denudation using the geothermal gradient in the preserved rocks is accurate only when the sedimentary cover is absent. In all other cases, denudation is significantly overestimated, demonstrating that an accurate estimate of the geothermal gradient within the eroded rock section is crucial for accurate determination of the amount of denudation using LTT.


Rapid exhumation of onshore Britain in the early Paleogene is consistent with rock uplift in response to the arrival of the Iceland mantle plume; however, the amount and extent of denudation is controversial (Green et al., 2012; Cogné et al., 2016). Apatite fission track (AFT) data show that latest Cretaceous temperatures in the English Lake District were >110 °C and were 50–90 °C in the surrounding regions (Fig. 2B), requiring >3 km of denudation centered on the ∼60-km-wide Lake District block (assuming a geothermal gradient of 30 °C/km; Green, 1986). This amount of exhumation is at odds with stratigraphic reconstructions (Holliday, 1993) and estimates based on the thickness of magmatic underplating (Tiley et al., 2004). An elevated geothermal gradient (∼60 °C/km; Green, 2002) was introduced to reconcile these conflicting interpretations and was explained as an effect of enhanced basal heat flow from the underplating melts (Green et al., 2012). However, deep-seated magmatism has a negligible effect on the thermal structure of shallow crust (see Fig. DR2 in the GSA Data Repository1; Brown et al., 1994). Holliday (1993) proposed that the high latest Cretaceous temperatures were a consequence of high heat flow and a now-eroded low-conductivity sedimentary cover. This interpretation has remained untested.

We collected new AFT data from the Lake District, south Scotland, and north Wales to explore the effect that a low-conductivity sedimentary cover and variable heat production may have on cooling histories (Fig. 2A). New AFT ages are youngest in the Lake District (48–75 Ma) and significantly older (>150 Ma) in south Scotland and north Wales. Thermal histories were derived using a Bayesian approach of data inversion in the QTQt software (Gallagher, 2012). The models indicate that latest Cretaceous temperatures were ∼120 °C in the Lake District and ∼70 °C in south Scotland and north Wales (Fig. 2B), similar to those presented by Green et al. (1997). The full data set and paleotemperatures estimates are available in the Data Repository.

The latest Cretaceous temperature pattern roughly correlates with the present-day surface heat flow (Fig. 2C; Busby et al., 2011). Heat flow exceeds 80 mW/m2 around the large granite batholiths in north England, where radiogenic heat production is 1.9–5.2 μW/m3 (Downing and Gray, 1986). These values are two to five times the average heat production in Phanerozoic crust (0.9–1.1 μW/m3; Jaupart and Mareschal, 2005). The AFT-derived thermal histories indicate an exhumation pulse that began in the latest Cretaceous, when the granites were overlain by Mesozoic mudstones and chalk (Holliday, 1993). The low thermal conductivity of these rock types (1.2–2.0 W/m/K; Downing and Gray, 1986) implies that they may have provided thermal insulation to enhance the temperatures in the granites at the time of denudation.

Effect of Thermally Heterogeneous Crust

In order to test the effect of variable heat production and sedimentary cover on the AFT-derived thermal histories, three sets of forward models (grand total of 56) for three crust thermal structures were determined using the Pecube numerical model (Braun et al., 2012): UC, a spatially uniform crust with a geothermal gradient of 20–30 °C/km; G0, a spatially uniform crust that includes four heat-producing bodies (one underneath the Lake District and three in south Scotland); and GS, a crust structure similar to G0, but covered by a blanket of low-conductivity sedimentary rocks. Models were run with variable thermal parameters and for different values of total rock uplift. The quality of the models is evaluated based on the misfit (μ) between the observed and predicted AFT ages. All modeling parameters, misfit values obtained for particular models, and a detailed description of the modeling procedure are given in Section DR-2c of the Data Repository. Misfit values are plotted versus total rock uplift in Figure 3. For a range of total uplift consistent with regional geological constraints, low misfit values (<5) are found only by models that consider the combined effect of high-heat-producing granites and a sedimentary blanket.

Figure 4 shows, for each modeled scenario, the predicted distributions of AFT ages for total rock uplift of 2.25 km, which, for a given topography, implies exhumation of 1.25–2.25 km. Measured AFT ages are significantly younger and more variable than those that the uniform crust model (UC) generates. The AFT age distribution is better reproduced if a heat-producing body is present (G0), but the modeled ages are still older than measured values. The best fit to the data (μ < 4) is obtained when, prior to exhumation, a low-conductivity sedimentary succession covers a heat-producing intrusion (GS). In this scenario, the AFT age pattern is the result of variable thermal properties of the rocks, not variable denudation.

The model has some limitations, such as a simplified spatial distribution of thermal parameters and a lack of pre-Cenozoic uplift and subsidence events. Although taking these aspects into account would further improve our estimates of denudation, they have little effect on the regional pattern of the AFT ages on which the model results are based (see Section DR-2d in the Data Repository).

The 3-D numerical models support the results of the 1-D model by demonstrating that crustal thermal heterogeneities have to be taken into consideration when deriving the amount of uplift and denudation from LTT data. In this case, the presence of a low-thermal-conductivity layer reduces the amount of total Cenozoic denudation required in the Lake District by 50%–70%, and the AFT-derived cooling is translated into 1.25–2.25 km of denudation. These amounts overlap with the stratigraphic estimates of eroded Mesozoic cover (1.2–1.75 km; Holliday, 1993), without the need to invoke elevated basal heat flow.


This study demonstrates that not accounting for the increased geothermal gradient in a low-thermal-conductivity sedimentary layer underlain by high-heat-producing rocks may lead to significant overestimation of the amount of denudation. This issue can be particularly important in LTT studies, as granitic rocks are favored because of their content of idiomorphic apatite and zircons and they are commonly exhumed by erosion of overlying sedimentary successions. For instance, in the European Alps, the Variscan granite batholiths have been exposed by post-orogenic erosion, subsequently buried beneath Paleozoic–Paleogene sediments, and finally exhumed during Alpine orogenesis, resulting in complex LTT age patterns (e.g., Glotzbach et al., 2011) that are likely influenced by the thermal properties of the eroded crust. Similar situations can be found in other orogenic belts (Carrapa et al., 2014; Hoke et al., 2015), at rifted and passive margins (Persano et al., 2006; Wildman et al., 2015), and in intra-plate settings (Majorowicz and Minea, 2012; Barbarand et al., 2013). This study shows that quantifying the amount and regional extent of denudation should be preceded by a careful examination of the thermal properties of both sampled and eroded rocks.


The funding for this study was provided by the University of Glasgow, UK, (the College of Science and Engineering Ph.D. Scholarship and the Sir Alwyn Williams Scholarship) and Natural Environment Research Council (NERC) standard grant (GU 52562-NE/H008276/1). Kerry Gallagher is thanked for providing the code for the 1-D modeling. Andy Gleadow, Zoltán Erdős, and one anonymous reviewer are thanked for their helpful and constructive comments.

GSA Data Repository item 2017256, the apatite fission track data set and modeling parameters, is available online at or on request from