How local crustal thermal properties influence the amount of denudation derived from low-temperature thermochronometry

Katarzyna Łuszczak, Cristina Persano, Jean Braun, and Finlay M. Stuart Institute of Geological Sciences, Polish Academy of Sciences, Twarda 51/55, 00-818 Warsaw, Poland School of Geographical and Earth Sciences, University of Glasgow, Glasgow G12 8QQ, UK Helmholtz Centre Potsdam, German Research Center for Geosciences (GFZ), Telegrafenberg, 14473 Potsdam, Germany Isotope Geoscience Unit, Scottish Universities Environmental Research Centre, East Kilbride G75 0QF, UK  In the Łuszczak et al. (2017) paper, we present the results of modeling the impact of spatially variable crustal heat production and thermal conductivity on apatite fission-track thermochronometry (AFTT). We welcome this opportunity to clarify the reasoning and conclusions of our study in response to the Comment of Westaway (2018). Firstly, we want to emphasize that the main aim of our study was not to provide a detailed estimate of Cenozoic denudation of Britain. Centralwest Britain was used as a case study. We clearly state that our model has simplifications, and that overcoming them will improve denudation estimates. In the first paragraph of his Comment, Westaway misrepresents our work by stating that our main conclusion is that previous studies of Green (2002) and Green et al. (2012) have overestimated denudation. We do not make such a statement anywhere in the paper. These papers suggest the need for a high paleo-geothermal gradient ( TP) —as we point out; the issue our paper addresses is the cause of the elevated TP. Westaway argues that our criticism of the Green et al. (2012) explanation of high TP is invalid and suggests that our assumption of instantaneous underplating lowers the impact of heating. However, instantaneous emplacement of a thick underplating layer influences the thermal structure of the upper crust more than when it is emplaced as successive thinner layers (Clift and Turner, 1998). To calculate the thermal perturbations caused by underplating, Green et al. (2012) applied the equation originally used to quantify thermal perturbations due to the removal of the entire lithospheric mantle. In this case, the crust is exposed to abnormally high temperatures for a long period of time; however, when underplating occurs, the thermal anomaly decays faster because the magmatic layer is within the lithosphere and is surrounded by cooler material; more so if the underplated material is emplaced as short-lived pulses. Setting the lower boundary of the crust to 1100 C for 10 m.y. is clearly incorrect, as it assumes that underplating (1) has the same thermal effect as the asthenosphere, and (2) lasted for 10 m.y., whereas the bulk of the British Paleogene Igneous Province was emplaced within a few million years (e.g. Wilkinson et al., 2017). In the third paragraph, Westaway points to a discordance between the late Cenozoic uplift rates used in our Pecube model and those proposed elsewhere. Previous studies of the Neogene and Quaternary uplift of Britain indicate several hundreds of meters of post mid-Pliocene denudation. However, we stress that AFTT data are insensitive to this, unless it occurs within a crust characterized by an extremely high geothermal gradient. The constant post-57 Ma uplift rate of 0.0088 mm/yr used in our study cannot be unequivocally constrained by our data, but was taken as the simplest solution to account for the remaining 500 m of Cenozoic denudation, as predicted by other studies. The rate at which this more-recent denudation occurred does not affect our model; more importantly, it lies beyond the resolution of AFTT data and thus beyond the scope of our paper. Later, Westaway states that we did not present any result for <2000 m of uplift. This is because such scenarios would significantly overestimate AFTT ages for all reasonable values of thermal conductivity. Subsequently, he links our results to the amount of Cenozoic denudation in Britain inferred from the isostatic model of Brodie and White (1994) and calls this model over-simplistic. The two cited papers, however, do not discuss the Brodie and White model and we are, therefore, not able to respond to this point. Westaway criticizes our work for a lack of data from high elevations, and we acknowledge the usefulness of vertical profiles in studying the impact of basal heat flow changes on TP. However, the main conclusion derived from our 1-D model is that the low conductivity of the sedimentary ‘blanket’ raises the overall temperature of the underlying basement, without affecting its geothermal gradient; reconstructing the thermal structure of the shallow crust from a present vertical profile, without considering the eroded section, will actually be misleading. The estimate of TP presented by Green (2002) lacks a verifiable explanation for the cause of the TP increase. Also, the robustness of the paleotemperature estimates, based on which the TP is calculated, is equivocal. The data were modeled individually rather than as a quasi-vertical profile, and for some samples, the AFTT ages were based on a statistically low number of counted grains and measured track lengths. Given that Green (2002) does not provide the uncertainties associated with his models, nor an explanation of how the model works, assessing his conclusions is impossible. The wide range of geothermal gradients that could fit the data of Green (2002) is not irrefutable evidence for high TP within the preserved rock section. Moreover, nowhere in our paper do we state that the heat flow, and thus TP in the preserved rocks, was not elevated. This is a possibility, and we emphasize that the denudation estimates presented in our paper are maxima. Models that consider both variable thermal properties of rocks and the possibility of a slightly elevated basal heat flow would be surely interesting. Nature is always much more complicated than the models we create to explain it; clearly stating these simplifications and the uncertainties associated with the data used to test the modeled scenarios, is, we feel, a strong way forward to ascertain the robustness of our conclusions.


INTRODUCTION
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/ m 3 ) 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.

QUANTIFYING GEOTHERMAL GRADIENTS
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/m 3 increase in H.
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.

CASE STUDY-CENTRAL WEST BRITAIN
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 Repository 1 ; Brown et al., 1994).Holliday (1993) proposed that the 1 GSA Data Repository item 2017256, the apatite fission track data set and modeling parameters, is available online at http://www.geosociety.org/datarepository /2017/ or on request from editing@geosociety.org.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/m 2 around the large granite batholiths in north England, where radiogenic heat production is 1.9-5.2μW/m 3 (Downing and Gray, 1986).These values are two to five times the average heat production in Phanerozoic crust (0.9-1.1 μW/m 3 ; 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.0W/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 AFTderived 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.

IMPLICATIONS
This study demonstrates that not accounting for the increased geothermal gradient in a lowthermal-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 postorogenic 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.

Figure 2 .
Figure 1.A-C: Plots showing effect of heat production (H) and thermal conductivity (k) on geothermal gradient in uppermost crust made of: heat-producing body of variable H and thickness, k = 2.5 W/m/K (constant), intruded in "normal" crust (A); 1-km-thick sedimentary layer with variable k and H = 0 µW/m 3 (constant), overlying "normal" crust (B); or 1-km-thick sedimentary layer with variable k (values in figure) and H = 0 µW/m 3 overlying 12-km-thick body (H varies, k = 2.5 W/m/K) (C).D: Change of temperature with depth (black lines) for different crustal compositions (GS-granite with sedimentary cover; CS-"normal" crust with sedimentary cover; UC-uniform "normal" crust; G0-granite only; z-elevation relative to paleosurface).Granite is 12 km thick with H = 5 µW/m 3 and k = 2.5 W/m/K.Blanket layer is 1 km thick with k = 1.5 W/m/K and H = 0 µW/m 3 .Shaded area represents rocks preserved after 2 km of denudation.Gray lines indicate predicted denudation if blanket layer is not accounted for.In all models, "normal" crust has k = 2.5 W/m/K and H = 1 µW/m 3 .

Figure 3 .
Figure 3. Misfit of Pecube forward models versus amount of total uplift for three analyzed scenarios: UC-uniform crust with geothermal gradient of 20-30 °C/km; G0-crust with heatproducing granite; GS-crust with granite and sedimentary cover.Yellow shaded area shows range of uplift values acceptable by regional geological constraints.

Figure 4 .
Figure 4. Pecube-derived forward models of distribution of apatite fission track (AFT) ages for total uplift of 2.25 km and different crustal compositions (UC-uniform crust; G0-crust with heat-producing granite; GS-crust with granite and sedimentary cover) for study area in Great Britain (see Fig. 2 for location).Lower graphs show predicted versus observed ages for sampled localities.Gray line is 1:1 reference line.Panel on left (Topo) shows digital elevation model, location of granite batholiths (dashed lines), and range of thermal parameters used for batholiths and sedimentary blanket: r-radius of batholith; A-heat production rate (LD-Lake District granite; Scot-Scottish granites); K-thermal diffusivity.