the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Short communication: age2exhume – a MATLAB/Python script to calculate steadystate vertical exhumation rates from thermochronometric ages and application to the Himalaya
Peter van der Beek
Taylor F. Schildgen
Interpreting cooling ages from multiple thermochronometric systems and/or from steep elevation transects with the help of a thermal model can provide unique insights into the spatial and temporal patterns of rock exhumation. Although several wellestablished thermal models allow for a detailed exploration of how cooling or exhumation rates evolved in a limited area or along a transect, integrating large, regional datasets in such models remains challenging. Here, we present age2exhume, a thermal model in the form of a MATLAB or Python script, which can be used to rapidly obtain a synoptic overview of exhumation rates from large, regional thermochronometric datasets. The model incorporates surface temperature based on a defined lapse rate and a local relief correction that is dependent on the thermochronometric system of interest. Other inputs include sample cooling age, uncertainty, and an initial (unperturbed) geothermal gradient. The model is simplified in that it assumes steady, vertical rock uplift and unchanging topography when calculating exhumation rates. For this reason, it does not replace more powerful and versatile thermal–kinematic models, but it has the advantage of simple implementation and rapidly calculated results. We also provide plots of predicted exhumation rates as a function of thermochronometric age and the local relief correction, which can be used to simply look up a firstorder estimate of exhumation rate. In our example dataset, we show exhumation rates calculated from 1785 cooling ages from the Himalaya associated with five different thermochronometric systems. Despite the synoptic nature of the results, they reflect known segmentation patterns and changing exhumation rates in areas that have undergone structural reorganization. Moreover, the rapid calculations enable an exploration of the sensitivity of the results to various input parameters and an illustration of the importance of explicit modeling of thermal fields when calculating exhumation rates from thermochronometric data.
 Article
(6884 KB)  Fulltext XML

Supplement
(452 KB)  BibTeX
 EndNote
The steady accumulation of thermochronometric data from around the world provides an opportunity to constrain spatial patterns of longterm (millionyear timescale) exhumation with high granularity over vast swaths of the Earth's surface. This information can, in turn, provide clues to the driving mechanisms of orogen development and landscape evolution. Several wellestablished thermal models can be used to extract detailed cooling histories or exhumation rates from input cooling ages spread over a limited area or along an elevation transect. However, integrating information from large datasets, comprising cooling ages from multiple thermochronometers spread over a wide region, remains challenging due to the lack of easytouse tools that will handle such vast, multisystem datasets.
The most advanced modeling tools in common use by the thermochronology community include Pecube (Braun et al., 2012), HeFTy (Ketcham, 2005), QTQt (Gallagher, 2012), and GLIDE (Fox et al., 2014). Pecube is unique in its ability to handle forward and inverse thermal–kinematic modeling of spatially distributed data, including the options for timevarying topography as well as spatially and temporally variable rock uplift driven by defined fault geometries and kinematics. This complexity, however, entails substantial setup requirements and relatively high computational demands, which tend to limit the spatial extent of modeled datasets to ∼10^{2}–10^{3} km^{2}. HeFTy and QTQt, in contrast, model thermal histories only, for individual samples or samples that are assumed to fall into a pseudovertical alignment. GLIDE (Fox et al., 2014) was developed with the aim of extracting exhumation histories from regional datasets. While powerful, the temporally and spatially continuous coverage of calculated exhumation rates that the model produces requires interpolations that can be challenging to interpret without careful consideration of the spatial and temporal distribution of the input data (Fox et al., 2014; Schildgen et al., 2018).
Here we present a simple thermal model, age2exhume, which is optimized to provide a synoptic overview of exhumation rates from large regional datasets. This model, inspired by the original age2edot code (Brandon et al., 1998), takes the form of a MATLAB or Python script that solves for steadystate exhumation rates from input thermochronometric ages, assuming vertical exhumation pathways and unchanging topography. A key difference between age2edot and age2exhume is that the former (despite its name) solves for ages given input exhumation rates, whereas our new model solves for exhumation rates given input ages. This difference makes age2exhume more suitable for calculating exhumation rates from regional datasets, since individual sample characteristics (e.g., an elevationdependent surface temperature and local relief correction), included together with age in an input file, can be used to calculate an exhumation rate for each sample. A preliminary version of this code was used to visualize regional thermochronometric datasets in Schildgen et al. (2018); here, we provide more detailed background to the model and incorporate the individual sample characteristics mentioned above into the revised model.
The regional (constant) inputs to the model include crustal thermal properties that can be approximated or derived from the literature (an initial, unperturbed geothermal gradient, thermal model thickness, and thermal diffusivity) and kinetic parameters for the relevant thermochronometric systems, for which default values are provided. Samplespecific inputs include a local relief factor that can be extracted using standard GIS functions from a digital elevation model, elevation, thermochronometric system, age, and age uncertainty. From our example dataset of 1785 cooling ages derived from five different thermochronometric systems in the Himalaya, steadystate, vertical exhumation rates with their uncertainties can be calculated within seconds on a standard laptop computer. Despite the synoptic nature of the results, we show how they reflect several fundamental features of the mountain belt, including strong regional differences that reflect known segmentation patterns and changing exhumation rates in areas that have undergone recent structural reorganization.
2.1 Existing thermal models: their applications and limitations
Brandon et al. (1998) presented a simple, firstorder approach to predict thermochronologic ages from input exhumation rates, in the form of a code called age2edot. Age2edot calculates a steadystate conductive–advective geotherm and uses the approach of Dodson (1973) to predict the coolingratedependent closure temperature of a given thermochronometric system. It then combines the predicted closure temperature and the steadystate geotherm to find the closure depth and subsequently calculates a thermochronometric age by dividing the closure depth by the input exhumation rate. Kinetic parameters required for the Dodson (1973) calculation of closure temperature (see Sect. 2.2 below) are derived from diffusion experiments for noblegasbased systems (i.e., (U–Th)$/$He and ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$) and from fitting an Arrhenius relation to experimental annealing data for fissiontrack systems (see Reiners and Brandon, 2006, for more detail). Simplifying assumptions in the age2edot approach include (1) thermal steady state, (2) vertical exhumation paths, (3) unchanging topography, and (4) constant exhumation rates over the modeled time span. The most recent version of the age2edot code was released more than 15 years ago (Ehlers et al., 2005) and, because it was distributed as a Microsoft Windows executable, it is now obsolete.
Willett and Brandon (2013) published a modification to the age2edot approach, in which the steadystate geotherm solution was replaced by an (inherently transient) halfspace solution, a correction for the sample elevation with respect to the regionally averaged elevation was introduced, and a bestfit exhumation rate is predicted from an input age and a modern (i.e., final) geothermal gradient. The code was provided as a MATLAB script. Although it is computationally efficient, two aspects of this model limit its use for modeling large regional datasets in our view; one is of a practical nature, whereas the other is more fundamental. The practical limitation lies in the need to provide a value (or bounding values) for the modern geotherm for each prediction. Although this requirement makes conceptual sense, since only the modern geotherm can potentially be measured, it is of limited practical use because geothermal gradients are generally not known at more than very coarse spatial resolution, particularly in mountain belts. Moreover, the requirement is impractical when dealing with large datasets of widely varying ages, as geothermal gradients vary strongly in regions of variable exhumation rates. If the estimated bounding geotherms are poorly estimated (e.g., too low or high for a given thermochronometric age), no exhumation rate is returned. The more fundamental issue lies in the choice of a thermal halfspace model, which leads to a strong sensitivity of the geotherm to exhumation rate and the persistence of transient thermal conditions even after several tens of millions of years of steady exhumation (Willett and Brandon, 2013). One type of data that allows assessing if, and how rapidly, thermal steady state might be achieved in mountain belts is detrital thermochronology from sedimentary sequences in foreland basins. Several such datasets show constant lag times (i.e., thermochronometric age minus depositional age), interpreted as recording establishment of thermal steady state in the source area after only a few million years, including in the western European Alps (Bernet et al., 2001, 2009), the central and eastern Himalaya (Bernet et al., 2006; Chirouze et al., 2013), the eastern Himalayan syntaxis (Bracciali et al., 2016; Lang et al., 2016; Govin et al., 2020), Taiwan (Kirstein et al., 2010), and the Southern Alps of New Zealand (Lang et al., 2020). As argued by Bracciali et al. (2016), modeling these constant lag times using a thermal halfspace model would require decreasing exhumation rates through time, with a rate of decrease that exactly offsets the transient upward advection of the geotherm, in all the above cases. More probably, these data indicate that the thermal halfspace model is not ideal for representing orogenic geotherms.
A completely different approach is taken by the thermalhistory modeling codes HeFTy (Ketcham, 2005) and QTQt (Gallagher, 2012). These codes aim at predicting a thermal history from thermochronometric ages and additional measurements (in particular fissiontrack length distributions, but also kinetic indicators) for single samples, although the most recent versions of these codes allow modeling suites of vertically offset samples. The output of these models, when run in inverse mode, is an optimal time–temperature history and its uncertainty. These thermal history results require assumptions about the past geothermal gradient to be translated to a burial/exhumation history. Gallagher and Brown (1999) and Kohn et al. (2002) spatially interpolated thermal histories derived from large numbers of individual samples, using a precursor of the QTQt code, and combined them with heatflow maps to derive regional to continentalscale images of denudation over geological time. This laborintensive approach requires multiple thermochronometric systems and/or tracklength data for each included sample in order to resolve meaningful thermal histories.
Pecube (Braun et al., 2012) is a threedimensional thermal–kinematic code that predicts thermochronometric ages for various userdefined tectonic and geomorphic scenarios, taking into account the spatial and temporal perturbation of the geotherm by rock advection and transient topography. Pecube allows modeling both vertical and nonvertical exhumation paths, the latter controlled by a simple faultkinematic model, and can be coupled to the neighborhood algorithm (Sambridge, 1999a, b) to run in inverse mode. The code has been used in a wide variety of tectonic and geomorphic settings (see Braun et al., 2012, for an overview), including at the scale of a small orogen (Curry et al., 2021). However, the fairly high computational demands of the code, particularly when run in inverse mode, make it best suited for models of more limited spatial extent (i.e., not exceeding several tens of km in length and width), where simple fault kinematics and/or spatially uniform rock uplift can reasonably represent the tectonic deformation patterns.
GLIDE (Fox et al., 2014) comprises a linear inverse method to infer spatial and temporal variations in exhumation rate from spatially distributed thermochronometer datasets. GLIDE uses a numerical thermal model with a flux boundary condition at the base. The inversion assumes vertical exhumation and a smooth spatial variation in exhumation rates that can be described by a spatial correlation function. In this way, it uses exhumation constraints from one sample to help constrain exhumation in nearby regions, producing exhumation histories that are continuous in space and time. However, it has been argued that the code translates abrupt spatial variations in thermochronological ages, such as across faults, into temporal increases in exhumation rates (Schildgen et al., 2018), unless the faults (or other features) are explicitly included in the correlation structure (Fox et al., 2014; Ballato et al., 2015). Willett et al. (2021) argued that such issues occur mainly in areas of insufficient data coverage without, however, quantifying this term; Schildgen et al. (2018) argued that most sampled regions on Earth with sharp spatial variations in exhumation have insufficient data coverage for unbiased prediction of exhumationrate histories using GLIDE, if those variations are not taken into account.
From the above abbreviated review, we conclude that a simple, firstorder method to assess large regional datasets in a consistent manner is currently lacking from the thermochronology toolbox. We aim to provide such a simple method with the age2exhume code.
2.2 age2exhume method
Figure 1 shows a sketch outline and flowchart for the age2exhume model. Input parameters for the model include the sealevel temperature T_{0}, atmospheric lapse rate H, the initial, unperturbed geothermal gradient G_{init}, thermal diffusivity κ, and model thickness L. The latter can represent the crustal thickness or, more appropriately, the maximum depth from which rocks have been exhumed, such as the depth to a regional detachment horizon. Input data for each sample include a thermochronometric age and its uncertainty at locations x and y, sample elevation h, and local relief correction Δh. Kinetic parameters for the main low to intermediatetemperature thermochronometric systems (apatite and zircon (U–Th)/He and fission track, mica ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$) are included as default values but can be modified if desired.
When calculating exhumation rates from thermochronometric ages, a local relief correction (Δh) is needed to account for the difference in elevation of a sample (h) relative to an averageelevation (h_{avg}) surface that mimics the shape of the closure isotherm (Stüwe et al., 1994; Braun, 2002). We follow the procedure of Willett and Brandon (2013) in estimating the shape of that surface by averaging surface topography over a circle with a radius of π×z_{c}, where z_{c} is an estimated closure depth for the relevant thermochronometric system. The local relief correction Δh is then calculated for each sample as
A brief guide for how to calculate this correction using a digital elevation model in ESRI ArcMap or in QGIS is provided in Appendix A. To predict a steadystate exhumation rate from a thermochronometric age, surface temperature, and the local relief correction, the model starts with an initial guess of the closure depth (z_{c}) and exhumation rate ($\dot{e}$) from an initial, unperturbed linear geothermal gradient (G_{init}), a nominal closure temperature (T_{c}), and a surface temperature (T_{s}):
T_{s} is calculated from an input sealevel temperature (T_{0}), the surfacetemperature lapse rate (H), and the sample elevation at the position of h_{avg}: ${T}_{s\left(h\right)}={T}_{\mathrm{0}}H{h}_{\text{avg}}$. We use h_{avg}, rather than the actual sample elevation for this surfacetemperature correction to simulate how surface temperature affects the thermal field at depth. For highertemperature thermochronometers with deeper closure depths, h_{avg} becomes more smoothed, and the associated impact of surface temperature on z_{c} is reduced. Note that the initial unperturbed geothermal gradient (G_{init}) is only used to calculate an appropriate basal temperature and to provide an initial estimate of the exhumation rate using Eqs. (2) and (3).
The model then iteratively adapts T_{c}, z_{c}, and $\dot{e}$ until convergence to a steadystate solution. Importantly, Δh is not recalculated after the initial estimate. Given the generally low sensitivity of Δh to moderate variations in z_{c}, we believe this simplification is worthwhile, considering the consequent reduced computational demands. At each iterative step, first the advective perturbation of the geotherm due to exhumation is calculated following Mancktelow and Grasemann (1997):
where T_{(z)} is the temperature at depth z, T_{L} is the temperature at the base (z=L) of the model (${T}_{L}={T}_{\text{avg}}+{G}_{\text{init}}L$, where T_{avg} is the temperature at the average elevation of the whole dataset), and κ is the thermal diffusivity. Equation (4) can be solved for the closure depth z_{c}:
Next, the closure temperature is reestimated as a function of the cooling rate at the closure depth. First, the depth derivative of Eq. (4) is used to estimate the geothermal gradient:
Equation (6) is evaluated at the closure depth z_{c}. Because $\dot{e}=\mathrm{d}z/\mathrm{d}t$, the cooling rate ($\dot{T}$) is
The model then uses the Dodson (1973) equation to relate closure temperature to cooling rate:
where E_{a} (activation energy), D_{0} (diffusivity at infinite temperature), and a (diffusion domain size) are experimentally determined kinetic parameters for each thermochronological system, A is a geometry factor, and τ (characteristic time) is
Once a new estimate for T_{c} is obtained, z_{c} is updated using Eq. (5) and a new estimate for the exhumation rate is obtained with Eq. (3). The model steps through Eqs. (3)–(9) iteratively (Fig. 1b) until the change in exhumation rate between successive steps ($\mathrm{\Delta}\dot{e}$) is smaller than a threshold value; here we use $\mathrm{\Delta}\dot{e}/\dot{e}<{\mathrm{10}}^{\mathrm{3}}$. To ensure smooth convergence, the exhumation rate used in each successive step is the average between the previous and the newly calculated rate.
3.1 General model predictions
Figures 2 and 3 show contours of predicted exhumation rates for different combinations of age and Δh; Fig. 2 shows results for moderate exhumation rates (<2 km Myr^{−1}) and thermochronometric ages up to 30 Ma, whereas Fig. 3 zooms in on the youngest ages (<5 Ma) and shows results for exhumation rates up to 5 km Myr^{−1}. Input parameters for these models are as in Table 1, except that a constant surface temperature (T_{s}) of 10 ^{∘}C was used, because absolute sample elevation is not included in these generic models. Kinetic parameters for the apatite (U–Th)$/$He (AHe) system are derived from Farley (2000), for the zircon (U–Th)$/$He (ZHe) system from Reiners et al. (2004), and for the apatite (AFT) and zircon (ZFT) fissiontrack systems from Reiners and Brandon (2006). These results can be thought of conceptually as showing age–elevation profiles for different constant exhumation rates, with elevation measured relative to an average regional elevation as defined in Sect. 2.2. They can also be used as a plotted lookup table for rapidly inferring exhumation rate from a given age, Δh combination.
3.2 Results from a Himalayan example dataset
Our example dataset from the Himalaya comprises 1785 thermochronologic ages compiled from papers published through July 2022; data sources are provided in the Supplement. We have excluded some reported ages from the Siwaliks (subHimalayan foldthrust belt), as that sedimentary unit commonly yields unreset ages. We have also excluded the western and eastern syntaxis regions, where extremely rapid exhumation is driven by processes that are different from those in the main part of the Himalaya (Zeitler et al., 2014; Butler, 2019). Finally, we exclude any preHimalayan ages (>60 Ma), as these are not directly linked to exhumation during Himalayan mountain building. Our dataset comprises 345 white mica ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ (MAr) ages, 236 ZFT ages, 783 AFT ages, 281 ZHe ages, and 140 AHe ages. All ages and sample details are included in a single Excel file, with columns that include a sample ID number, latitude, longitude, elevation, Δh value, age, 1σ age uncertainty, and a numeric code for the thermochronologic system (Schildgen and van der Beek, 2022b). Table 1 shows the parameters we assume for the surface temperature (T_{0}, H) and the thermal model (L, G_{init}, κ). Kinetic parameters used for the AHe, AFT, ZHe, and ZFT systems are the same as for the general model predictions presented in Sect. 2.1 above; we used the parameters from Hames and Bowring (1994) for the MAr system.
A map of the calculated exhumation rates for the Himalaya (Fig. 4) shows exhumation plotted such that rates derived from lowertemperature systems plot on top of those from highertemperature systems. When the symbol for a lowertemperature system is darker than the symbol of a highertemperature system plotted below it, this implies that exhumation rates have slowed through time. Conversely, a lighter color for the lowertemperature system plotted over a highertemperature system implies exhumation rates have increased through time. The map reveals patterns in space and time that reflect wellknown structural patterns of the range. In general, a band of rapid exhumation rates occurs at the topographic front of the high Himalaya, with slower rates recorded to the north and south. Within this band, the highest rates are generally recorded by the lowertemperature AHe and AFT thermochronometers, suggesting increasing exhumation rates with time. Note that such variable exhumation rates recorded by different colocated thermochronometers formally violate the assumption of constant exhumation rates through time implicit in the model. The rates inferred from the highertemperature thermochronometers should therefore be considered rough estimates only; they will generally be overestimated in the case of increasing rates through time, and the corresponding rate change will therefore be underestimated. The focused rapid rates at the foot of the high Himalaya together with an increase in exhumation rates for lowertemperature systems are consistent with exhumation being driven by thrusting over a largescale ramp in the Main Himalayan Thrust (MHT), the interface between the underthrusting Indian continent and the overlying Himalayan units, often associated with duplex development (e.g., Robert et al., 2009; Herman et al., 2010; Coutand et al., 2014; Dal Zilio et al., 2021; van der Beek et al., 2023).
The highest exhumation rates (>2 km Myr^{−1}) outside of the Himalayan syntaxes occur in central Nepal (∼84^{∘} E), Sikkim (∼88^{∘} E), the Kumaun Himalaya (∼80^{∘} E), and the Sutlej Valley (∼78^{∘} E). High rates (between 1 and 2 km Myr^{−1}) are recorded along the high Himalayan front throughout northwest India (∼76–80^{∘} E) and more sporadically in eastern Nepal (∼87^{∘} E) and western Bhutan (∼89^{∘} E). The lowest exhumation rates along the high Himalayan topographic front (<0.8 km Myr^{−1}) are found in Kashmir (west of ∼75^{∘} E), western Nepal (∼81^{∘} E), and from western Bhutan (∼90^{∘} E) to the east. These lateral variations in exhumation rates have been interpreted as reflecting lateral variations in the presence/absence and geometry (location, height, and dip) of the midcrustal ramp in the MHT, together with duplex formation and local outofsequence thrusting (Hubbard et al., 2021; Dal Zilio et al., 2021; van der Beek et al., 2023). In some of the more slowly exhuming regions, in particular in Bhutan, exhumation rates appear to be decreasing through time, with lowertemperature systems recording lower exhumation rates than highertemperature systems. Decreasing exhumation rates in Bhutan can be linked to slowing convergence across the Bhutan Himalaya due to transfer of deformation to the Shillong Plateau to the south (Clark and Bilham, 2008; Coutand et al., 2014, 2016). Similar to the caveats described above concerning increasing exhumation rates, in areas of decreasing exhumation rates, the change in rates through time recorded by different systems will also be underestimated.
The above example illustrates how this method can rapidly provide internally consistent estimates of exhumation rates from multiple thermochronometers from different elevations over a large region. Inferred patterns of exhumation rates can be linked to structural and geophysical observations of orogen segmentation, as above, or to orogenwide topographic measures for assessing firstorder linkages between exhumation rates and morphology (e.g., Clubb et al., 2022).
4.1 Importance, uncertainties, and sensitivity
An advantage of the rapid calculations performed by age2exhume is that it is easy to explore the sensitivity of the calculated exhumation rates to different input parameters (i.e., samplespecific information and crustal/thermal properties), in addition to evaluating how the iterative method compares to simpler estimates of exhumation rates. Regarding the latter, we can compare calculated exhumation rates from age2exhume to those that would be obtained by assuming a simple linear geotherm and fixed nominal closure temperature, T_{c}. Figure 5a compares “initial” exhumation rates, calculated using Eqs. (2) and (3) (hence, a linear geotherm and fixed T_{c}), with the final exhumation rates predicted by age2exhume, which incorporate perturbations to the geotherm and T_{c}. Initial exhumation rates are calculated using the same thermal parameters of Table 1 and nominal closure temperatures of 70 ^{∘}C for the AHe system, 120 ^{∘}C for the AFT system, 180 ^{∘}C for the ZHe system, 220 ^{∘}C for the ZFT system, and 350 ^{∘}C for the MAr system. The comparison shows that for exhumation rates up to ∼0.5 km Myr^{−1}, there is little difference between the two methods (Fig. 5a). At higher exhumation rates, the methods deviate substantially, with the initial estimate systematically overestimating the exhumation rate. For example, at exhumation rates ≥2 km Myr^{−1}, overestimates mostly fall between 100 % and 300 %. These findings can be explained by considering the relative importance of two competing influences on the closure depth z_{c} (Fig. 1a), which directly determines the exhumation rate (Eq. 3). On the one hand, higher cooling rates – linked to higher exhumation rates – lead to an increase in T_{c}, and hence a deepening of z_{c} (Eqs. 8 and 9). On the other hand, the advective perturbation of the geotherm due to exhumation, which forces an upward deflection of isotherms, leads to a shallowing of z_{c} for any T_{c} (Eq. 5). The degree of advective perturbation of the geotherm is characterized by the nondimensional Péclet number: $\mathit{Pe}=\dot{e}L/\mathit{\kappa}$ (e.g., Braun et al., 2006); the predicted (surface) geothermal gradient thus increases with increasing exhumation rate (Fig. 6). With higher exhumation rates, the effect of upward, advective perturbation of isotherms on z_{c} dominates over the effect of the increasing T_{c} on z_{c}. The scatter in the amount of overestimation, in particular for the lowertemperature AFT and AHe systems, is linked to the effect of including Δh, which is more important for shallower z_{c} (Eq. 3).
But how important are these differences in the method of calculating exhumation rates relative to the uncertainties in any calculated rate? The uncertainties in reported ages are just one component of the total uncertainty that one can consider in an exhumationrate calculation, but the direct propagation of age uncertainty into the uncertainty on an inferred exhumation rate provides a simple means of comparison (Fig. 5b). Because of the nonlinear relationship between age and exhumation rate, the uncertainties in exhumation rates are asymmetric, with ${\dot{e}}_{\text{max}}\dot{e}>\dot{e}{\dot{e}}_{\text{min}}$. The bulk of the relative uncertainties in exhumation rates associated with age uncertainty lie between 10 % and 50 %, and they are not strongly dependent on exhumation rate. Highertemperature systems (ZHe, ZFT and MAr) are generally associated with lower exhumationrate uncertainties (<10 %) because of the smaller age uncertainties associated with these systems. In contrast, AFT data can have uncertainties that exceed 100 %, because low track counts due to low Ucontents and/or young ages yield large age uncertainties. Some large relative uncertainties in the AHe and ZHe systems at lower exhumation rates (<1 km Myr^{−1}) are probably associated with larger intergrain scatter in ages due to compositional and grainsize effects that become more important at lower cooling and exhumation rates (e.g., Whipp et al., 2022, and references therein). Overall, however, the bulk of the uncertainties in exhumation rate are smaller than the differences between the initial and final exhumation rates shown in Fig. 5a for exhumation rates ≳0.5 km Myr^{−1}. This comparison implies that the thermal effects of exhumation significantly affect inferred exhumation rates in tectonically active areas.
The importance of including samplespecific information in exhumationrate calculations is illustrated in Fig. 7a and b. Our comparison of exhumation rates calculated with a constant surface temperature, T_{s}, versus those calculated with T_{s} dependent on elevation shows a relatively small effect, with differences mostly less than 10 %. However, for the lowtemperature thermochronometers AHe and AFT at exhumation rates <1 km Myr^{−1}, differences can reach 20 % (Fig. 7a). The effect of the local relief correction, Δh, for each sample is generally more important. Although the magnitude of Δh tends to be reduced for the lowertemperature systems (because their closure isotherms more closely mimic surface topography), any given Δh has a stronger impact on exhumation rates for lowtemperature systems (with shallower z_{c}) compared to hightemperature systems (Fig. 7b and Eq. 3). Moreover, the effects are asymmetric: negative Δh values lead to a larger correction in exhumation rates compared to positive Δh values. For example, a Δh of +1 km will lead to a ca. 20 % change in calculated exhumation rate for the AFT system, whereas a Δh of −1 km will lead to a 30 % to 50 % change (Fig. 7b). This asymmetry results from the nonlinear effect of exhumation rate on z_{c}: positive Δh values will lead to increased predicted exhumation rates, but these increases will be partly offset by the resulting advective perturbation of the geotherm. In contrast, the decreased predicted exhumation rates for negative Δh values will be less modified by advective effects. The importance of including Δh when calculating exhumation rates is further emphasized when considering that samples are more commonly collected from valley bottoms (with negative Δh values) than ridgetops. Our Himalayan example dataset bears this out: the histogram of Δh values is skewed toward negative values, with a median Δh of −0.53 km (Fig. 7b inset).
We next explore the sensitivity of calculated exhumation rates to crustal parameters, including the model thickness L (Fig. 7c) and the initial, unperturbed geothermal gradient G_{init} (Fig. 7d). These plots show the percent change in predicted exhumation rates when changing these two parameters to either a higher or a lower value relative to the default values of L=30 km and G_{init}=25 ${}^{\circ}\mathrm{C}\phantom{\rule{0.125em}{0ex}}{\mathrm{km}}^{\mathrm{1}}$. Decreasing L from 30 to 20 km leads to higher predicted exhumation rates (by up to ∼40 %), whereas increasing L from 30 to 40 km leads to lower predicted exhumation rates (by up to $\sim \mathrm{20}$ %), with the magnitude of the effect increasing with exhumation rate (Fig. 7c). This behavior can be understood by considering the effect L has on the advective perturbation of the geotherm, through the Péclet number (see above): the Péclet number is linearly dependent on L so that, for constant exhumation rate $\dot{e}$ and diffusivity κ, increasing L will lead to a stronger perturbation of the geotherm and thus a shallower closure depth for any thermochronometer. The sensitivity of the predictions to G_{init} is of similar magnitude when considering changes from 25 to 30 ${}^{\circ}\mathrm{C}\phantom{\rule{0.125em}{0ex}}{\mathrm{km}}^{\mathrm{1}}$ or from 25 to 20 ${}^{\circ}\mathrm{C}\phantom{\rule{0.125em}{0ex}}{\mathrm{km}}^{\mathrm{1}}$ (Fig. 7d), but in this case, the effect is strongest for relatively low exhumation rates and thus relatively unperturbed geothermal gradients.
4.2 Limitations and recommended use
The assumptions underlying the model limit its use, strictly speaking, to settings where both topography and exhumation rates are temporally constant throughout the time period considered and exhumation is mainly vertical. The requirement for constant exhumation rates is due to both the use of a steadystate solution for the advectively perturbed geothermal gradient and the use of the Dodson (1973) approach to estimate closure temperatures. These assumptions break down in cases where thermal and exhumation histories are more complex, in particular when they include phases of burial and heating. However, as our Himalayan example shows, our approach can provide firstorder information on spatiotemporal patterns of exhumation, highlighting regions of accelerating versus decelerating exhumation. This result comes with the caveat that the model will systematically underestimate the change in exhumation rate, as discussed above. The Himalayan example also shows that firstorder results can be obtained in a setting where horizontal advection of rocks is important; however, in this case, the interpretation of possible accelerations or decelerations in exhumation rate should take the regional structure and kinematics into account. For instance, accelerated exhumation may be due to rocks moving over a flattoramp transition in a crustalscale décollement, rather than a temporal change in tectonic or climatic drivers.
Our analysis of the importance of including advective perturbation of the geotherm in the thermochronometric age predictions shows that this effect is not significant for exhumation rates ≲0.5 km Myr^{−1} (Fig. 5a). At these relatively low rates of exhumation and cooling, kinetic effects also become important in controlling thermochronometric ages (e.g., Whipp et al., 2022) and these are not included in our model. The model is thus best suited for the analysis of regional datasets from rapidly and continuously exhuming regions, e.g., tectonically active mountain belts. The relief correction included in the model makes it suitable to handle data that were collected at widely varying elevations.
The model assumptions of a constant basal temperature together with an input model thickness are unlikely to be valid over long timescales, and in many cases can only be estimated roughly. However, the speed with which exhumation rates can be calculated from our model enables users to easily investigate the sensitivity of their results to these estimated values. Moreover, while these thermal parameters change the absolute values of the predicted exhumation rates, they affect all predictions similarly (if not equally). Therefore, their influence on spatial patterns in exhumation rates or the correlation of exhumation rates with other metrics will be limited.
We provide three different versions of the model in the form of MATLAB scripts on Zenodo (van der Beek and Schildgen, 2022): (1) a basic version that takes a single age–Δh pair as input and returns a single exhumation rate; (2) a version for which a range of thermochronologic ages and Δh values are provided and that returns a lookup table of exhumation rates (used in Sect. 3.1 and Figs. 2 and 3); and (3) a version that reads an input file of sample locations, elevation, thermochronologic system, age, and uncertainty and returns a table of exhumation rates with uncertainty, closure depths, and surface steadystate geotherms for each sample (used in Sect. 3.2 and Fig. 4). A correctly formatted input file is also included in the Zenodo repository. We anticipate the latter version to be most useful, and therefore we also provide a Python script with that same functionality (Schildgen and van der Beek, 2022a). Alternatively, Figs. 2 and 3 of this paper can be used to simply look up appropriate exhumation rates for a given age–Δh combination, but note that these figures are plotted for particular values of the input parameters G_{init}, L, and κ.
4.3 Concluding remarks
The model presented here, age2exhume, enables a firstorder, synoptic view of spatial and temporal variations in exhumation rates, calculated in a rapid, selfconsistent manner from different thermochronometers. The main advantage of our approach over the version of age2edot presented by Willett and Brandon (2013) is that our model does not require the final geothermal gradient as input, but only the initial, unperturbed geotherm. This aspect of our model makes it easily applicable to regions with strongly varying exhumation rates, which are expected to have a wide range of modern geothermal gradients. The modern geothermal gradient, when known, adds an additional constraint to the model solution. However, for many regions of the world, particularly for mountain belts, modern geothermal gradients are essentially unknown. In our entire Himalayan study region, for instance, the global heatflow database (https://ihfciugg.org, last access: 10 October 2022) does not contain a single data point. Although there are data both for the Tibetan Plateau to the north and the Ganges foreland basin to the south, these are not useful for assessing the perturbed geotherm within the mountain belt. Our model does provide the predicted steadystate surface geotherm as output, so it can be compared to any potential measurements (Fig. 6).
Our model assumes steadystate exhumation, unchanging topography, and vertical exhumation pathways, so it is only appropriate for obtaining firstorder, synoptic overviews of exhumationrate patterns in regions of relatively rapid, continuous exhumation. Nevertheless, in the case where ages from multiple thermochronometers are available from individual samples or from samples in close proximity to one another, differences in exhumation rates derived from those ages can be used to map out where changes in exhumation rates have likely occurred, and thus highlight regions where more advanced thermal modeling could be used to extract nonsteadystate exhumation histories. The rapidity with which our model calculates regional patterns of exhumation rates also allows testing its sensitivity to the different input parameters.
To calculate Δh, Willett and Brandon (2013) suggest calculating a mean value of elevation for a circle that has a radius equal to π×z_{c}, where z_{c} is the closure depth of the system. This calculation can be done with standard operations in a geographic information system (GIS) or other tools designed to work with continuous raster datasets. The following instructions can be followed to calculate Δh values within ArcMap from ESRI (version 10.8.1) or within QGIS (version 3.26). We have not tested if the instructions are easily applicable to earlier versions of the software. Nevertheless, small modifications to these procedures can likely be found by searching on the names of the functions described below. Importantly, regardless of the software package used to calculate Δh, the spatial extent of the DEM should extend beyond the limits of the sample points, with a buffer zone at least equal to the highest radius that will be considered. For example, the DEM should extend at least ca. 40 km beyond the spatial extent of the sample data to prevent edge effects from affecting Δh calculations for the MAr system.
ESRI ArcMap
In ESRI's ArcMap version 10.8.1, the mean elevation can be calculated using the Focal Statistic function, found within the “Spatial Analyst Tools – Neighborhood” tools in Arc Toolbox. The Focal Statistic function provides an option to average values over a moving circular window with a radius defined by map units or by a number of pixels. For example, for a standard 90 m resolution Shuttle Radar Topography Mission (SRTM) DEM, and for a desired z_{c} of 2000 m (e.g., for the AHe system), the radius of the circle should be 6280 m, which is approximately equivalent to 70 pixels. To efficiently calculate Δh for all samples in a large dataset, it is practical to take advantage of the “Raster Calculator” (Spatial Analyst Tools – Map Algebra) and the “Extract Values to Points” functions (Spatial Analyst Tools – Extraction). The Raster Calculator can be used to subtract the smoothed DEM calculated in the previous step from the modern DEM. This operation will produce a continuous raster dataset of Δh values. The “Extract Values to Points” function samples a raster at the position of each sample data point and adds the extracted value to a new column (“field”) in the attribute table of the shapefile. Although the exact procedure described here may differ for other versions of ArcMap, general functions to calculate focal statistics, perform arithmetic operations on raster datasets, and automatically extract values from rasters at the location of sample points can be found in many versions of the software.
QGIS
In QGIS 3.26, a procedure to find the average elevation over a defined circular search area can be accomplished with the SAGA plugin, which can be installed directly from the “Plugin” menu and then “Manage and install plugins”. After installation, the SAGA tools can be found within the “Processing” menu, then “Toolbox”. Within SAGA, go to the “Raster Filter” options and then select “Simple filter”. The filter option should be set to “Smooth” (to calculate an average value), and the Kernel type set to “Circle”. The radius should be set to the number of pixels that will provide the correct radius length. Like in the example above, for a standard 90 m resolution DEM and a desired z_{c} of 2000 m (for the AHe system), the radius of the circle (π×z_{c}) should be 6280 m, which is approximately equal to 70 pixels. Next, the Raster Calculator within QGIS can be used to calculate Δh values over the extent of the DEM, by subtracting the smoothed DEM calculated in the previous step from the original DEM. Finally, to extract the calculated Δh values for each sample point, within the standard Processing Toolbox under the “Raster analysis” heading is the function “Sample raster values”. With this tool, the point layer containing the sample points should be given as the “Input layer” and the raster of Δh values should be given as the “Raster layer”. The output point file includes all of the attributes of the original point layer but adds a column containing the extracted Δh value for each point. That file, by default, is only saved to memory. To save it permanently, the small squareshaped icon to the right of the layer name can be clicked to bring up a dialog box that allows saving the file to a defined location.
The MATLAB scripts for three versions of the age2exhume code, together with an input file, are included in the Zenodo repository: age2exhume MATLAB scripts (https://doi.org/10.5281/zenodo.7341603, van der Beek and Schildgen, 2022). The Python version of age2exhume, together with an input file, can be downloaded from the Zenodo repository: age2exhume Python script (https://doi.org/10.5281/zenodo.7341690, Schildgen and van der Beek, 2022a).
Data used in the example dataset were compiled from the sources listed in the Supplement. An Excel file containing the full dataset and calculated exhumation rates is included in the Zenodo repository: Thermochronology dataset for Himalaya (https://doi.org/10.5281/zenodo.7053115, Schildgen and van der Beek, 2022b).
The supplement related to this article is available online at: https://doi.org/10.5194/gchron5352023supplement.
Both authors contributed equally to the development of the code, the compilation of the Himalayan dataset, sensitivity analyses, and writing of the manuscript. TS wrote the Python version of the code.
The contact author has declared that neither of the authors has any competing interests.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Taylor F. Schildgen acknowledges support from the ERC Consolidator grant no. 863490 GyroSCoPe; Peter van der Beek acknowledges support from the ERC Advanced grant no. 834271 COOLER. We thank reviewers Matthew Fox and David Whipp for constructive comments that improved the quality of the manuscript and allowed us to refine some aspects of the code. We thank Associate Editor Marissa Tremblay for efficient handling of the manuscript. We also thank David Whipp for his freely accessible online “GeoPython” course, which provided the base training needed for writing the Python version of the code.
This research has been supported by the European Research Council, H2020 European Research Council.
This paper was edited by Marissa Tremblay and reviewed by David Whipp and Matthew Fox.
Ballato, P., Landgraf, A., Schildgen, T. F., Stockli, D. F., Fox, M., Ghassemi, M. R., Kirby, E., and Strecker, M. R.: The growth of a mountain belt forced by baselevel fall: Tectonics and surface processes during the evolution of the Alborz Mountains, Iran, Earth Planet. Sc. Lett., 425, 204–218, https://doi.org/10.1016/j.epsl.2015.05.051, 2015.
Bernet, M., Zattin, M., Garver, J. I., Brandon, M. T., and Vance, J. A.: Steadystate exhumation of the European Alps, Geology, 29, 35–38, https://doi.org/10.1130/00917613(2001)029<0035:sseote>2.0.co;2, 2001.
Bernet, M., van der Beek, P., Pik, R., Huyghe, P., Mugnier, J.L., Labrin, E., and Szulc, A.: Miocene to Recent exhumation of the central Himalaya determined from combined detrital zircon fissiontrack and $\mathrm{U}/\mathrm{Pb}$ analysis of Siwalik sediments, western Nepal, Basin Res., 18, 393–412, https://doi.org/10.1111/j.13652117.2006.00303.x, 2006.
Bernet, M., Brandon, M., Garver, J., Balestrieri, M. L., Ventura, B., and Zattin, M.: Exhuming the Alps through time: clues from detrital zircon fissiontrack thermochronology, Basin Res., 21, 781–798, https://doi.org/10.1111/j.13652117.2009.00400.x, 2009.
Bracciali, L., Parrish, R. R., Najman, Y., Smye, A., Carter, A., and Wijbrans, J. R.: PlioPleistocene exhumation of the eastern Himalayan syntaxis and its domal “popup”, EarthSci. Rev., 160, 350–385, https://doi.org/10.1016/j.earscirev.2016.07.010, 2016.
Brandon, M. T., RodenTice, M. K., and Garver, J. I.: Late Cenozoic exhumation of the Cascadia accretionary wedge in the Olympic Mountains, northwest Washington State, Geol. Soc. Am. Bull., 110, 985–1009, https://doi.org/10.1130/00167606(1998)110<0985:lceotc>2.3.co;2, 1998.
Braun, J.: Quantifying the effect of recent relief changes on age–elevation relationships, Earth Planet. Sc. Lett., 200, 331–343, https://doi.org/10.1016/s0012821x(02)006386, 2002.
Braun, J., van der Beek, P., and Batt, G. E.: Quantitative Thermochronology: Numerical methods for the interpretation of thermochronological data, Cambridge University Press, 271 pp., https://doi.org/10.1017/CBO9780511616433, 2006.
Braun, J., van der Beek, P., Valla, P., Robert, X., Herman, F., Glotzbach, C., Pedersen, V., Perry, C., SimonLabric, T., and Prigent, C.: Quantifying rates of landscape evolution and tectonic processes by thermochronology and numerical modeling of crustal heat transport using PECUBE, Tectonophysics, 524–525, 1–28, https://doi.org/10.1016/j.tecto.2011.12.035, 2012.
Butler, R. W. H.: Tectonic evolution of the Himalayan syntaxes: the view from Nanga Parbat, Geol. Soc. London Spec. Publ., 483, 215–254, https://doi.org/10.1144/sp483.5, 2019.
Chirouze, F., Huyghe, P., van der Beek, P., Chauvel, C., Chakraborty, T., DupontNivet, G., and Bernet, M.: Tectonics, exhumation, and drainage evolution of the eastern Himalaya since 13 Ma from detrital geochemistry and thermochronology, Kameng River Section, Arunachal Pradesh, Geol. Soc. Am. Bull., 125, 523–538, https://doi.org/10.1130/b30697.1, 2013.
Clark, M. K. and Bilham, R.: Miocene rise of the Shillong Plateau and the beginning of the end for the Eastern Himalaya, Earth Planet. Sc. Lett., 269, 336–350, https://doi.org/10.1016/j.epsl.2008.01.045, 2008.
Clubb, F. J., Mudd, S. M., Schildgen, T. F., van der Beek, P. A., Devrani, R., and Sinclair, H. D.: Himalayan valleyfloor widths controlled by tectonics rather than water discharge, Research Square [preprint], https://doi.org/10.21203/rs.3.rs2065309/v1, 11 October 2022.
Coutand, I., Whipp, D. M., Grujic, D., Bernet, M., Fellin, M. G., Bookhagen, B., Landry, K. R., Ghalley, S. K., and Duncan, C.: Geometry and kinematics of the Main Himalayan Thrust and Neogene crustal exhumation in the Bhutanese Himalaya derived from inversion of multithermochronologic data, J. Geophys. Res., 119, 1446–1481, https://doi.org/10.1002/2013jb010891, 2014.
Coutand, I., Barrier, L., Govin, G., Grujic, D., Hoorn, C., DupontNivet, G., and Najman, Y.: Late MiocenePleistocene evolution of IndiaEurasia convergence partitioning between the Bhutan Himalaya and the Shillong Plateau: New evidences from foreland basin deposits along the Dungsam Chu section, eastern Bhutan, Tectonics, 35, 2963–2994, https://doi.org/10.1002/2016tc004258, 2016.
Curry, M. E., van der Beek, P., Huismans, R. S., Wolf, S. G., Fillon, C., and Muñoz, J.A.: Spatiotemporal patterns of Pyrenean exhumation revealed by inverse thermokinematic modeling of a large thermochronologic data set, Geology, 49, 738–742, https://doi.org/10.1130/g48687.1, 2021.
Dal Zilio, L., Hetényi, G., Hubbard, J., and Bollinger, L.: Building the Himalaya from tectonic to earthquake scales, Nat. Rev. Earth Environ., 2, 251–268, https://doi.org/10.1038/s43017021001431, 2021.
Dodson, M. H.: Closure temperature in cooling geochronological and petrological systems, Contrib. Mineral. Petr., 40, 259–274, https://doi.org/10.1007/bf00373790, 1973.
Ehlers, T. A., Chaudhri, T., Kumar, S., Fuller, C. W., Willett, S. D., Ketcham, R. A., Brandon, M. T., Belton, D. X., Kohn, B. P., Gleadow, A. J. W., Dunai, T. J., and Fu, F. Q.: Computational tools for lowtemperature thermochronometer interpretation, Rev. Mineral. Geochem., 58, 589–622, https://doi.org/10.2138/rmg.2005.58.22, 2005.
Farley, K. A.: Helium diffusion from apatite: General behavior as illustrated by Durango fluorapatite, J. Geophys. Res., 105, 2903–2914, https://doi.org/10.1029/1999jb900348, 2000.
Fox, M., Herman, F., Willett, S. D., and May, D. A.: A linear inversion method to infer exhumation rates in space and time from thermochronometric data, Earth Surf. Dynam., 2, 47–65, https://doi.org/10.5194/esurf2472014, 2014.
Gallagher, K.: Transdimensional inverse thermal history modeling for quantitative thermochronology, J. Geophys. Res., 117, B02408, https://doi.org/10.1029/2011jb008825, 2012.
Gallagher, K. and Brown, R.: The Mesozoic denudation history of the Atlantic margins of southern Africa and southeast Brazil and the relationship to offshore sedimentation, Geol. Soc. London Spec. Publ., 153, 41–53, https://doi.org/10.1144/gsl.sp.1999.153.01.03, 1999.
Govin, G., van der Beek, P., Najman, Y., Millar, I., Gemignani, L., Huyghe, P., DupontNivet, G., Bernet, M., Mark, C., and Wijbrans, J.: Early onset and late acceleration of rapid exhumation in the Namche Barwa syntaxis, eastern Himalaya, Geology, 48, 1139–1143, https://doi.org/10.1130/g47720.1, 2020.
Hames, W. E. and Bowring, S. A.: An empirical evaluation of the argon diffusion geometry in muscovite, Earth Planet. Sc. Lett., 124, 161–169, https://doi.org/10.1016/0012821x(94)000794, 1994.
Herman, F., Copeland, P., Avouac, J.P., Bollinger, L., Mahéo, G., Le Fort, P., Rai, S., Foster, D., Pêcher, A., Stüwe, K., and Henry, P.: Exhumation, crustal deformation, and thermal structure of the Nepal Himalaya derived from the inversion of thermochronological and thermobarometric data and modeling of the topography, J. Geophys. Res., 115, B06407, https://doi.org/10.1029/2008jb006126, 2010.
Hubbard, M., Mukul, M., Gajurel, A. P., Ghosh, A., Srivastava, V., Giri, B., Seifert, N., and Mendoza, M. M.: Orogenic segmentation and its role in Himalayan mountain building, Front. Earth Sci., 9, 641666, https://doi.org/10.3389/feart.2021.641666, 2021.
Ketcham, R. A.: Forward and inverse modeling of lowtemperature thermochronometry data, Rev. Mineral. Geochem., 58, 275–314, https://doi.org/10.2138/rmg.2005.58.11, 2005.
Kirstein, L. A., Fellin, M. G., Willett, S. D., Carter, A., Chen, Y. G., Garver, J. I., and Lee, D. C.: Pliocene onset of rapid exhumation in Taiwan during arccontinent collision: new insights from detrital thermochronometry, Basin Res., 22, 270–285, https://doi.org/10.1111/j.13652117.2009.00426.x, 2010.
Kohn, B. P., Gleadow, A. J. W., Brown, R. W., Gallagher, K., O'Sullivan, P. B., and Foster, D. A.: Shaping the Australian crust over the last 300 million years: insights from fission track thermotectonic imaging and denudation studies of key terranes, Aust. J. Earth Sci., 49, 697–717, https://doi.org/10.1046/j.14400952.2002.00942.x, 2002.
Lang, K. A., Huntington, K. W., Burmester, R., and Housen, B.: Rapid exhumation of the eastern Himalayan syntaxis since the late Miocene, Geol. Soc. Am. Bull., 128, 1403–1422, https://doi.org/10.1130/b31419.1, 2016.
Lang, K. A., Glotzbach, C., Ring, U., Kamp, P. J. J., and Ehlers, T. A.: Linking orogeny and orography in the Southern Alps of New Zealand: New observations from detrital fissiontrack thermochronology of the Waiho1 borehole, Earth Planet. Sc. Lett., 552, 116586, https://doi.org/10.1016/j.epsl.2020.116586, 2020.
Mancktelow, N. S. and Grasemann, B.: Timedependent effects of heat advection and topography on cooling histories during erosion, Tectonophysics, 270, 167–195, https://doi.org/10.1016/s00401951(96)00279x, 1997.
Reiners, P. W. and Brandon, M. T.: Using thermochronology to understand orogenic erosion, Annu. Rev. Earth Pl. Sc., 34, 419–466, https://doi.org/10.1007/9783540486848, 2006.
Reiners, P. W., Spell, T. L., Nicolescu, S., and Zanetti, K. A.: Zircon (UTh)$/$He thermochronometry: He diffusion and comparisons with ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ dating, Geochim. Cosmochim. Ac., 68, 1857–1887, https://doi.org/10.1016/j.gca.2003.10.021, 2004.
Robert, X., van der Beek, P., Braun, J., Perry, C., Dubille, M., and Mugnier, J. L.: Assessing Quaternary reactivation of the Main Central thrust zone (central Nepal Himalaya): New thermochronologic data and numerical modeling, Geology, 37, 731–734, https://doi.org/10.1130/g25736a.1, 2009.
Sambridge, M.: Geophysical inversion with a neighbourhood algorithm – I. Searching a parameter space, Geophys. J. Int., 138, 479–494, https://doi.org/10.1046/j.1365246x.1999.00876.x, 1999a.
Sambridge, M.: Geophysical inversion with a neighbourhood algorithm – II. Appraising the ensemble, Geophys. J. Int., 138, 727–746, https://doi.org/10.1046/j.1365246x.1999.00900.x, 1999b.
Schildgen, T. F. and van der Beek, P. A.: Age2exhume python script, Zenodo [code], https://doi.org/10.5281/zenodo.7341690, 2022a.
Schildgen, T. F. and van der Beek, P. A.: Thermochronology dataset for Himalaya, Zenodo [data set], https://doi.org/10.5281/zenodo.7053115, 2022b.
Schildgen, T. F., van der Beek, P. A., Sinclair, H. D., and Thiede, R. C.: Spatial correlation bias in lateCenozoic erosion histories derived from thermochronology, Nature, 559, 89–93, https://doi.org/10.1038/s4158601802606, 2018.
Stüwe, K., White, L., and Brown, R.: The influence of eroding topography on steadystate isotherms. Application to fission track analysis, Earth Planet. Sc. Lett., 124, 63–74, https://doi.org/10.1016/0012821x(94)000689, 1994.
van der Beek, P. A. and Schildgen, T. F.: age2exhume Matlab scripts, Zenodo [code], https://doi.org/10.5281/zenodo.7341603, 2022.
van der Beek, P. A., Thiede, R. C., Gahalaut, V. K., and Schildgen, T. F.: Topographic and thermochronologic constraints on the geometry of the Himalayan décollement and its lateral variations, in: Himalaya, Dynamics of a Giant, vol. 1, edited by: Cattin, R. and Epard, J. L., Wiley/ISTE Editions, in press, 2023.
Whipp, D. M., Kellett, D. A., Coutand, I., and Ketcham, R. A.: Short communication: Modeling competing effects of cooling rate, grain size, and radiation damage in lowtemperature thermochronometers, Geochronology, 4, 143–152, https://doi.org/10.5194/gchron41432022, 2022.
Willett, S. D. and Brandon, M. T.: Some analytical methods for converting thermochronometric age to erosion rate, Geochem. Geophy. Geosy., 14, 209–222, https://doi.org/10.1029/2012gc004279, 2013.
Willett, S. D., Herman, F., Fox, M., Stalder, N., Ehlers, T. A., Jiao, R., and Yang, R.: Bias and error in modelling thermochronometric data: resolving a potential increase in PlioPleistocene erosion rate, Earth Surf. Dynam., 9, 1153–1221, https://doi.org/10.5194/esurf911532021, 2021.
Zeitler, P. K., Meltzer, A. S., Brown, L., Kidd, W. S. F., Lim, C., and Enkelmann, E.: Tectonics and topographic evolution of Namche Barwa and the easternmost Lhasa block, Tibet, Geol. Soc. Am. Spec. Pap., 507, 23–58, https://doi.org/10.1130/2014.2507(02), 2014.