the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
An efficient approach for inverting rock exhumation from thermochronologic age–elevation relationship
Yuntao Tian
Lili Pan
Guihong Zhang
Xinbo Yao
This study implements the leastsquares inversion method for solving the exhumation history from the thermochronologic age–elevation relationship (AER) based on the linear equation among exhumation rate, age and total exhumation from the closure depth to the Earth surface. Modeling experiments suggest significant and systematic influence of initial geothermal model, the a priori exhumation rate and the time interval length on the a posteriori exhumation history. Lessons learned from the experiments include that (i) the modern geothermal gradient can be used for constraining the initial geothermal model, (ii) a relatively high a priori exhumation rate would lead to systematically lower a posteriori exhumation and vice versa, (iii) the variance of the a priori exhumation rate controls the variation in the inverted exhumation history, and (iv) the choice of time interval length should be optimized for resolving the potential temporal changes in exhumation. To mitigate the dependence of inverted erosion history on these initial parameters, we implemented a new stepwise inverse modeling method for optimizing the model parameters by comparing the observed and predicted thermochronologic data and modern geothermal gradients. Finally, method demonstration was performed using four synthetic datasets and three natural examples of different exhumation rates and histories. It is shown that the inverted rock exhumation histories from the synthetic datasets match the whole picture of the “truth”, although the temporal changes in the magnitude of exhumation are underestimated. Modeling of the datasets from natural samples produces geologically reasonable exhumation histories. The code and data used in this work are available on Zenodo (https://doi.org/10.5281/zenodo.10839275).
 Article
(7099 KB)  Fulltext XML

Supplement
(1863 KB)  BibTeX
 EndNote
Quantifying rock exhumation from the Earth interior to the surface is important information for better understanding many geological problems, ranging from orogenic growth (e.g., Zeitler et al., 2001; Whipp et al., 2007) and decay (e.g., House et al., 2001; Hu et al., 2006) to resource and hydrocarbon evaluation and exploration (e.g., Armstrong, 2005; Mcinnes et al., 2005), as well as the underpinning endogenic and exogenic processes and their interactions (e.g., Burbank et al., 2003; Fox et al., 2015; Tian et al., 2015). Various experimental and modeling methods have been invented for estimating the rock exhumation at different crustal levels (e.g., Braun, 2003; Reiners and Brandon, 2006; Anderson et al., 2008; Braun et al., 2012; Fox et al., 2014).
One type of the methods for estimating the rock exhumation in the middle and upper crust relies on thermochronologic cooling ages acquired from by noble gas and fissiontrack dating of a series of accessory minerals, such as Ar–Ar, fissiontrack and $\left(\mathrm{U}\text{}\mathrm{Th}\right)/\mathrm{He}$ analyses (Ault et al., 2019, and references therein). Based on the closure temperature theory (Dodson, 1973), assuming monotonic cooling, a thermochronologic age records the time duration that a rock cooled through the corresponding closure temperature, which is a function of the kinematics describing fissiontrack annealing and noble gas diffusion and rock cooling rate (Dodson, 1973). If the depth of the closure temperature isotherm can be estimated from the crustal temperature field, a timeaveraged exhumation rate can be obtained from the cooling age.
Based on the thermochronologic methods and thermoexhumation modeling, many analytical and numerical tools have been implemented for inverting the exhumation and/or the associated cooling history from thermochronologic data. These tools have different functions, such as inverting temperature history (Laslett et al., 1987; Ketcham, 2005; Gallagher, 2012), determining timeaveraged exhumation rates (Brandon et al., 1998; Ehlers, 2005; Willett and Brandon, 2013; Glotzbach et al., 2015; Van Der Beek and Schildgen, 2023), finding spatiotemporal changes in exhumation (Sutherland et al., 2009; Herman et al., 2013; Fox et al., 2014; Willett et al., 2020), and defining the evolution of exhumation in two or three dimensions given a tectonic framework (Batt and Brandon, 2002; Braun, 2003; Van Der Beek et al., 2010; Valla et al., 2011; Braun et al., 2012).
A convincing estimate of the exhumation history for a region requires both a proper sampling strategy for thermochronologic data and a robust modeling approach for exhumation inversion, especially when the rock exhumation and its spatiotemporal changes are tectonically controlled (Ehlers and Farley, 2003; Schildgen et al., 2018). A routine and efficient sampling strategy acquires thermochronologic ages from an elevation transect over a significant relief and a relatively confined spatial distance. Plotting the age versus elevation, i.e., the age–elevation relationship (AER), and analyzing the slope changes in the plot can provide firstorder understanding of the exhumation history (Fitzgerald et al., 1986). Because both the subsurface geothermal field and closure temperature of thermochronometers are functions of the thermal advection and cooling during rock exhumation (e.g., Dodson, 1973; Brandon et al., 1998), as well as the longwavelength topography (Braun, 2002; Ehlers and Farley, 2003; Glotzbach et al., 2015), estimating reliable exhumation rates requires us to account for temporal variations in the thermal field caused by changes in the thermal and kinematic boundary conditions.
Fox et al. (2014) reported a linear inversion modeling method that solves exhumation history from AER, given a combination of a priori exhumation rates and assumed geothermal parameters. However, as shown in that study, the inverted exhumation history depends highly on these a priori values and geothermal assumptions. Building on that study, here we provide a detailed test for the method and report an improved modeling method that makes use of both the AER and the modern geothermal gradient for inverting exhumation history.
Our inversion of exhumation from thermochronologic data followed the linear inversion approach of Fox et al. (2014). Rock exhumation from the closure depth of a thermochronometer, z_{c}, to the Earth's surface can be described as an integral of the exhumation ($\dot{e}$) from the cooling age (τ) to the present (Brandon et al., 1998; Fox et al., 2014). For a set of correlated bedrock samples with a shared history of exhumation rates ($\dot{\mathit{e}}$), their thermochronologic ages (A) and the corresponding closure depths (z_{c}) can be expressed by the following equation:
where A is a model matrix, with n rows (the total number of samples) and m columns (the total number of time intervals). Each row of the matrix is a discretization of a sample age, which is composed of a number of time lengths (Δt) followed by an age residual (R_{i}) and a number of zeros. The $\dot{\mathit{e}}$ is an mlength vector of exhumation rates, and the z_{c} is an nlength vector of closure depths.
This linear equation can be solved using the leastsquares regression approach assuming the Gaussian uncertainties and a priori mean exhumation rate (${\dot{\mathit{e}}}_{\text{pr}}$) and associated variance (σ_{pr}) (Tarantola, 2005; Fox et al., 2014). Such an approach requires an m×msized parameter covariance matrix, C, and an n×nsized data covariance matrix, C_{ε}, which includes the uncertainties in the closure depths. These two matrices can be constructed as Eqs. (2) and (3), respectively.
where ${\dot{e}}_{\text{pr}}$ and σ_{pr} are the a priori exhumation and the associated variance and the ε_{i} is analytical uncertainty in the age data. The construction of the data covariance matrix assumes the age data are uncorrelated. Worth noting is that previous studies used different constructions of the data covariance, changing from using the analytical age uncertainties (Fox et al., 2014, 2015) to constant values (Jiao et al., 2017; Stalder et al., 2020).
Given the above model parameters, the Eq. (1) has a maximum likelihood solution for the exhumation rate vector:
where ${\dot{\mathit{e}}}_{\text{pr}}$ is a nlength vector of ${\dot{\mathit{e}}}_{\text{pr}}$ and z_{c} is the nlength vector of closure depths calculated using a combination of exhumation and geothermal model parameters (see Sect. 3). The ${\dot{\mathit{e}}}_{\text{po}}$ is the posteriori maximum likelihood estimate of the exhumation rate, with a covariance matrix, C_{po}, which provides an estimate of the uncertainties in the model parameters (Eq. 5).
The method also provides a model resolution matrix, R, which gives a measure of how well the model estimates correspond to the true values:
Inversion of the exhumation using the Eq. (1) requires accurate estimates of the closure depths of the thermochronologic ages (z_{c}), i.e., the depth of the closure temperatures (Fig. 1). The latter can be determined by modeling the temperature of the crust using a 1D thermalkinematic model, which accounts for heat conduction, advection and production (Turcotte and Schubert, 2002):
where A_{b} is the heat production (in °C Myr^{−1}). This function can be numerically solved using a Crank–Nicolson time integration with a set of initial and boundary conditions, such as an initial geothermal gradient (G0) at the start time of the model and surface temperature (T_{s}) (Turcotte and Schubert, 2002; Fox et al., 2014).
The closure temperature (T_{c}) of a thermochronometer is a function of cooling rate ($\dot{T}$) at the closure time and kinetic parameters of helium and argon diffusion and fissiontrack annealing in mineral phases (Dodson, 1973):
where Ω and E_{a} are the diffusion frequency factor normalized by the mineral size and geometry and activation energy, respectively. Parameter R is the gas law constant. See reviews by Reiners and Brandon (2006) for the Ω and E_{a} parameter values for different thermochronometers.
The cooling rate ($\dot{T}$) can be computed from the derivative of transient geotherms, T_{m}(t,z) that can be computed using Eq. (7) (Fox et al., 2014):
where $\dot{e}$ is unknown exhumation that can be computed through the Eq. (1).
Combining the Eqs. (7)–(9), the closure depth of a thermochronological system (z_{c,m}) can be numerically computed. This depth also needs a topographic correction because of the topographic perturbation, p, on the isotherms (Braun, 2002; Ehlers and Farley, 2003; Fox et al., 2014; Glotzbach et al., 2015). Such a perturbation can be determined by the following equation (Mancktelow and Grasemann, 1997; Fox et al., 2014):
where γ_{a} is the atmospheric lapse rate and γ_{0} and ${\mathit{\gamma}}_{{z}_{\mathrm{m}}}$ are the thermal gradients at the model surface and at the depth z_{m}. The h(λ) is a cosine function expression of the model surface topography, which can be determined using the discrete fast Fourier transform at the frequency domain. Here we use the SRTM30 data for computing the topography of regions of interest.
Finally, the closure depth of the z_{c} is corrected by the topographic perturbation (e.g., Brandon et al., 1998):
where z_{c,m} is the closure depth calculated using the 1D geothermal model, p and h are the topographic perturbation and elevation difference with respect to the mean elevation at the sample site (Fig. 1), and the i denotes the ith age.
As shown by the Eqs. (7), (8) and (9), the closure depth is a nonlinear function of rock cooling and exhumation. Therefore, the problem of interest is nonlinear, which can be addressed by iterative numerical modeling methods. In this work, the solution of exhumation is approximated by coupling and iterating the linear inversion and closure depth modeling. As shown in Tarantola (2005) and Fox et al. (2014), the algorithm converges in a few iterations and produces stable outputs.
Quantitative model assessment relies on a misfit value, i.e., the difference between observed and predicted ages weighted by the observed analytical uncertainty:
where τ_{obs,i} and τ_{prd,i} are the observed and predicted ith age calculated from the exhumation history and ε_{i} is the uncertainty in the observed ith age. Following Fox et al. (2014), both the a priori and a posteriori misfits, Φ_{τ,pr} and Φ_{τ,po}, are determined for the models. The difference between these two misfit values provides a measure of the model improvements. A smaller posteriori misfit value indicates an improved model result and vice versa.
To evaluate the geothermal parameters, we also determined the misfit value of the predicted to the observed modern geothermal gradient value using the following equation:
where γ_{prd} and γ_{obs} are the predicted and observed geothermal gradients and ε_{γ} is the uncertainty in the observed value. Because the depth–temperature curves are slightly nonlinear, the predicted geothermal gradient (γ_{prd}) is calculated as a mean value for the upper 1 km of the model. Similar to the assessment of age data, we also determined the a priori and a posteriori misfits, Φ_{γ,pr} and Φ_{γ,po} values for assessing the geothermal parameters.
Following Willett and Brandon (2013) and Fox et al. (2014), here we use the published AFT data acquired from Denali Massif (Fitzgerald et al., 1995) for method demonstration (Fig. 2a). A break in slope is shown by the AER at ∼7–6 Ma, indicating a coeval change in slope, i.e., the apparent exhumation rate (Fitzgerald et al., 1995), increasing from 0.17 ± 0.04 to 1.2 ± 0.6 km Myr^{−1} (Fig. 2b). AER regression of young dates from the lower part of the transect (between 4.3–2.0 km) also predicts a closure depth that is the intercept at −3.3 ± 3.4 km (Fig. 2b). However, using the present geothermal gradient (38.9 °C km^{−1}) (Fox et al., 2014), a nominal closure temperature of AFT method (110 °C) (Reiners and Brandon, 2006) and a −12 °C surface temperature (Fox et al., 2014), the closure depth is predicted as ∼3.1 km beneath the mean elevation (∼4 km), which is equivalent to an elevation of ∼0.9 km. This closure depth is significantly higher than the intercept (−3.3 ± 3.4 km). Such a difference indicates the AER slope of the lower part overestimates the exhumation rates since ∼7–6 Ma.
Following the protocol outlined in Fox et al. (2014), the reference inverse model uses the following parameters: a start time at 25 Ma, a time interval (Δt) of 2.5 Myr, a 4020 m mean elevation, a −12 °C surface temperature, an a priori exhumation rate of 0.5 ± 0.15 km Myr^{−1}, a 24 °C km^{−1} initial geothermal gradient, a 38.9 °C km^{−1} present geothermal gradient, a model block with a thickness of 80 km and a 30 km^{2} Myr^{−1} thermal diffusivity.
The exhumation history output of the reference model is shown in Fig. 3. The inversion results reveal an more than 2fold increase in exhumation rate to a value of ∼0.6 km Myr^{−1} at 7.5 Ma (Fig. 3b), consistent with the development of the break in slope in the AER. The model also shows a gradual decrease in exhumation rate from a priori exhumation rate (0.5 km Myr^{−1}) to 0.3 km Myr^{−1} from 25 to 7.5 Ma. The invariant exhumation during the starting stage resulted from the fact that all ages are younger than 17.5 Ma, and thus the data have no resolution for the time span. These results are similar to those of Fox et al. (2014). The a posteriori misfit for the age is 1.88, significantly smaller than that of the a priori model (4.51), suggesting the improvement by the inverse modeling (Fig. 3b). Such a model also provides reasonable fit to the modern temperature field, as shown by the small misfit (0.39) in the geothermal gradient (Fig. 3b).
The resolution of the inverted exhumation history can be assessed by the resolution matrix R (Eq. 6). Imaging of the matrix shows the model provides no resolution for the time period before 17.5 Ma (Fig. 3c), consistent with the fact that the oldest input age is younger than 16.1 ± 0.9 Ma. For the time span between 15 and 5 Ma, the model resolution is high, as shown by the diagonal elements of the matrix, with the highest resolution at 7.5–5 Ma span, including eight age date points (Fig. 3c). The two most recent phases of exhumation (5–0 Ma) are less resolved, as shown by the nearly equal resolution values for the two phases, i.e., the latest four pixels of the matrix (Fig. 3c). This is because no input ages fall into this time span, when the modeled exhumation results are timeaveraged values. The slight decrease in the last stage reflects changes in geothermal gradient.
For assessing the correlation among model parameters, the calculated covariance matrix is scaled by the diagonal covariance matrix (Fox et al., 2014):
The correlation matrix for the reference model is shown in Fig. 3d. The diagonal correlation values are 1 and offdiagonal ones are dominantly negative, indicating anticorrelated uncertainties (Fig. 3d), which suggests exhumation parameters were not resolved independently by the modeling. In fact, it is expected to have the anticorrelation because, given two steps of rock exhumation, decreasing the exhumation during one step would increase that of the other step.
Here we use the Denali data set for demonstrating the influences of (i) the initial geothermal parameters, (ii) the a priori mean and (iii) variance values of the exhumation rates, and (iv) the time interval length on the inverted exhumation history. Also discussed in this section are the solutions for optimizing the model setup for these parameters.
6.1 Dependence on initial thermal model
Different initial model geothermal parameters would lead isotherms to shift either downward to greater depths or upwards to the Earth surface and either compression or expansion among isotherms. Therefore, the initial thermal models have systematic influence on the closure depths and consequently the a posteriori exhumation.
This is demonstrated by modeling experiments presented in Fig. 4. Using a relatively low initial geothermal gradient produces relatively high a posteriori exhumation rates (comparing the models shown in Fig. 4a–f) and vice versa. Such an influence is significant even for the time and elevation intervals with multiple age constraints (10–5.0 Ma). For example, using relatively low geothermal gradients of <22 °C km^{−1} would yield significantly higher average exhumation rates of >0.75 km Myr^{−1} for the last two stages (<5 Ma) (Fig. 4a–c) than those (<0.6 km Myr^{−1}) using higher initial geothermal gradients of ≥26 °C km^{−1} (Fig. 4d–f). Further, it is also shown that models using higher and lower prior geothermal gradients of <20 °C km^{−1} (Fig. 4a and b) and >30 °C km^{−1} (Fig. 4e and f) yield worse misfits (${\mathrm{\Phi}}_{\mathit{\gamma},\text{po}}>\mathrm{1}$) for the observed presentday geothermal gradient than those (${\mathrm{\Phi}}_{\mathit{\gamma},\text{po}}<\mathrm{1}$) using medium initial gradients (22–26 °C km^{−1}) (Figs. 3, 4c and d).
These results highlight the importance of taking geothermal parameters into account in inverting the exhumation history and model evaluation. We proposed running a set of models using different a priori geothermal parameters, especially the initial geothermal gradient, to search for the proper initial geothermal setup that provides reasonable fits to both the ages and the modern geothermal gradient (see Sect. 7 for details).
6.2 Dependence on the a priori exhumation rate
Both the mean and variance of the a priori exhumation rate have important influences on the model solution for the maximum likelihood estimation method. Our modeling experiments show that the mean value of the a priori exhumation has systematic influences on the inverted exhumation. Similar to the reference model, exhumation of the preceding three stages (25–17.5 Ma) without age constraints is the same as the a priori input. For the following stages, a relatively high mean value of the a priori exhumation results in relatively low a posteriori exhumation rates (comparing different models presented in Fig. 5). For example, models using the mean a priori exhumation of ≤0.4 km Myr^{−1} yields an a posteriori exhumation of 0.5–0.9 km Myr^{−1} for the stages <7.5 Ma (Fig. 5a–c), whereas those using a higher a priori value (≥0.6 km Myr^{−1}) result in an a posteriori exhumation of 0.45–0.6 km Myr^{−1} for the same stages (Fig. 5d–f). This is because a relatively high a priori value, which would be used for calculating thermal models, would lead to a quicker increase in geothermal gradient and thus relatively shallow closure depths and relatively low exhumation rates.
The variance of the a priori exhumation rate has an important influence on both the exhumation rates and the posterior variance. Models with lower a priori variances yield fewer variations in the a posteriori exhumation history and vice versa (comparing models in Fig. 6). Further, models using the input variance of the a priori exhumation of 0.2–0.3 km Myr^{−1} (40 %–60 % of the mean value), the variation in the inverted exhumation history becomes stable (Figs. 3, 6c and d). Given the uncertainty in the input age data, which is often 10 %–20 % at two sigma levels, larger variance of the inverted exhumation would be unreasonable (Fig. 6e and f), especially when multiple age data points are available at different elevations.
We proposed running a set of models using different a priori mean values of erosion rates to search for the one that provides appropriate fits to both the ages and the modern geothermal gradient. As to the a priori variance, we propose using a value 30 %–70 % of the a priori erosion rate. Future applications of the method may need to test a set of the variance inputs so as to get a stable exhumation output. Larger a priori variance would lead to larger uncertainties for the exhumation rates, which is unreasonable and not meaningful for geological studies.
6.3 Dependence on time interval length
Constraining the onset time of major changes in exhumation rates is one of the important tasks for inverting the exhumation history from thermochronologic data. Using a large time interval length cannot accurately capture the potential transition time of exhumation rates. As shown in Fig. 7b–d, models using time lengths of ≤3.5 Myr show an abrupt increase in exhumation at 7–6 Ma, consistent with that shown in AER plot. However, the models using a large time interval length (≥4.5 Ma) overestimate the onset time of the enhanced exhumation (Fig. 7e and f). Further, a relatively short time length would smooth temporal changes in exhumation rates, leading to an underestimating of the variations. For example, as shown in the Fig. 7a, the model using a relatively short time length (0.5 Ma) yields an exhumation variation between 0.35–0.60 km Myr^{−1}, significantly lower than those using relatively larger time interval lengths (Fig. 7b–f). In addition, a shorter time length also significantly increases the computational time and resources, especially when processing a large number of vertical transects.
Given the interests in major exhumation changes, we propose that the time interval length (Δt) should be optimized for constraining the transitional time and the associated exhumation changes. Therefore, the time interval length should be set as the absolute uncertainty at two sigma levels at the break point (τ_{b}) (Eq. 15). If the break point is unclear in AER, we suggest using the absolute uncertainty at two and three sigma levels at the median age value ($\stackrel{\mathrm{\u0303}}{\mathit{\tau}}$) (Eq. 15) to focus on the time intervals where ages cluster.
where δ is the relative age uncertainty at two sigma levels, varying between 10 %–20 % among different studies. Following this method, the Denali case should use a time length of ∼1.5 Ma (7 Ma × 20 %), slightly lower than that used in the reference model (Fig. 3).
Following the modeling protocol outlined above, a stepwise modeling guideline is developed for addressing the model dependencies on the initial geothermal parameter, the a priori exhumation rates and time interval length. As illustrated in the Fig. 8, the approach includes the following three steps.
 i.
Estimating a timeaveraged erosion rate. Dividing each nominal closure depth, which can be estimated from the nominal closure temperatures and the modern geothermal gradient, by the corresponding age results in a timeaveraged erosion rate. Following this, a mean value can be determined by averaging the rates. Such a mean value and assumed variance (30 %–50 % in this work) will be used as the a priori erosion rate.
 ii.
Optimizing the fit to the modern geothermal gradient. This step runs a set of inversion models (20 in this work) using different geothermal gradients, ranging from 50 % to 120 % of the modern value, together with the a priori erosion rate estimated in the first step, for determining the initial geothermal gradient that yields the maximum fit to the modern value, i.e., the minimum Φ_{γ} (Eq. 13).
 iii.
Optimizing the fit to both the age data and the geothermal gradient. Given the model dependence on the geothermal parameters (see Sect. 6.1), a comprehensive evaluation of the models should assess not only the age misfit (Φ_{τ}) but also that of the geothermal gradient (Φ_{γ}). In the third step, a set of inversion models (20 in this work) are run using different a priori erosion rates, changing from 10 % to 200 % of the mean value estimated in the first step, together with the estimated geothermal gradient by the second step, to search for the model that provides the best fit to both the age data and the modern geothermal gradient. This study uses the following compound misfit function to evaluate the models:
$$\begin{array}{}\text{(16)}& \mathrm{\Phi}={\mathrm{\Phi}}_{\mathit{\tau}}+{\mathrm{\Phi}}_{\mathit{\gamma}}/\sqrt{N}\phantom{\rule{0.125em}{0ex}},\end{array}$$where Φ_{τ} and Φ_{γ} are misfit values for the age and geothermal gradient calculated using the Eqs. (12) and (13) and N is the number of age inputs. Dividing Φ_{γ} by the square root of N in this equation, as also done for calculating the Φ_{τ} (Eq. 12), means that the modern geothermal gradient is given the same weight as an age input for evaluating the model.
We firstly test our stepwise inversion scheme by using synthetic datasets generated by thermokinematic models modified from Braun et al. (2012) (their Fig. 9). The synthetic age dataset is produced by Pecube using the following parameters: a steadystate topography with a 20 km wavelength and a 2 km relief, a model block thickness of 30 km with a basal temperature of 600 °C, a thermal diffusivity of 25 km^{2} Myr^{−1}, a sea level temperature of 10 °C, and a lapse rate of 5 °C km^{−1}. Worth noting is that these parameters are the same as Braun et al. (2012). For model details, see Braun et al. (2012). For the model setup, see Fig. S1 in the Supplement.
Synthetic AFT and AHe ages (Table S1 in the Supplement) were calculated for both surface and borehole samples for four different exhumation histories. The synthetic models I and II are characterized by a sudden decrease in exhumation rate from 1 to 0.1 km Myr^{−1} (model I, the same as the that shown in Fig. 9 of Braun et al., 2012) and 0.3 km Myr^{−1} (model II) at 5 Ma, respectively. Models III and IV include a sudden increase in exhumation rate from 0.3 km Myr^{−1} (model III) and 0.1 km Myr^{−1} (model IV) to 1 km Myr^{−1} at 5 Ma. All models start from 40 Ma. Except for the synthetic age data (plotted in the first row of Fig. 9), these four models generate modern geothermal gradients of 26.5, 28.6, 35.5 and 34 °C km^{−1} for the uppermost 2 km crust, respectively.
The inversion of rock exhumation history used a start time of 20 Ma and a time interval length of 1.0 Myr for all synthetic datasets, which were assigned with a 6 % uncertainty. As shown by the modeling output visualized in Fig. 9a, our inversion of the rock exhumation from the synthetic dataset I finds an optimal initial geothermal gradient of 22 °C km^{−1} and a priori rate of 0.85 ± 0.25 km Myr^{−1} and yields a decrease in exhumation rates from ∼0.9 km Myr^{−1} (before 6 Ma) to 0.3–0.1 km Myr^{−1} (4–0 Ma) via a gradual decrease during 6–4 Ma. The data have no resolution for the exhumation history before 10 Ma. Compared to the synthetic model (abrupt decrease from 1 to 0.1 km Myr^{−1} at 5 Ma), the rates before 5 Ma are underestimated by 0.1 km Myr^{−1}, whereas the values after 5 Ma are overestimated by 0.1–0.3 km Myr^{−1}.
The inversion for synthetic dataset II results in an optimal initial geothermal gradient of 21.7 °C km^{−1}, an a priori rate of 0.81 ± 0.24 km Myr^{−1} and an increase in exhumation rates from ∼0.85 (before 5 Ma) km Myr^{−1} to 0.4–0.5 km Myr^{−1} (4–0 Ma) via a gradual decrease during 5–4 Ma (Fig. 9b). Compared to the synthetic model (abrupt decrease from 1 to 0.3 km Myr^{−1} at 5 Ma), the rates before 5 Ma are underestimated, whereas the values before 5 Ma are overestimated by ∼0.1–0.2 km Myr^{−1}.
The inversion for synthetic dataset III yields an optimal initial geothermal gradient of 24.3 °C km^{−1}, an a priori rate of 0.55 ± 0.17 km Myr^{−1} and a decrease in exhumation rates from ∼0.45–0.3 km Myr^{−1} (before 5 Ma) to 1.0 km Myr^{−1} (3–0 Ma) via a gradual increase during 5–3 Ma (Fig. 9c). Compared to the synthetic model (abrupt decrease from 0.3 to 1.0 km Myr^{−1} at 5 Ma), the rates during 5–3 Ma are underestimated, whereas the rates before 5 Ma are overestimated by 0–0.15 km Myr^{−1}.
The inversion for the synthetic dataset IV produces an optimal initial geothermal gradient of 24.5 °C km^{−1}, an a priori rate of 0.25 ± 0.08 km Myr^{−1} and an increase in exhumation rates from ∼0.1–0.2 km Myr^{−1} (before 5 Ma) to 1.0 km Myr^{−1} (3–0 Ma) via a gradual decrease during 5–3 Ma (Fig. 9d). Compared to the synthetic model (abrupt decrease from 1 to 0.3 km Myr^{−1} at 5 Ma), the rates before 5 Ma are slightly overestimated, whereas the values during 5–3 Ma are underestimated.
To summarize, the inverted rock exhumation histories for the four synthetic datasets match the whole picture of the synthetic “truth”, but the variations in exhumation are underestimated, and the sharp changes at 5 Ma are smoothed. It is worth noting that inversions using only surface samples produce similar inversion results (Fig. S2 in the Supplement).
Below we use three examples to demonstrate our new method. The Denali data are used again for demonstrating the efficiency of our method in finding the proper initial geothermal gradient and the a priori exhumation rate. Following this, we further test our method using the Himalayan Dhanladar range and KTB borehole (the Continental Deep Drilling Project in Germany) thermochronologic data for representing regions of fast and slow erosion, respectively.
9.1 The Denali transect
Using the stepwise inversion modeling guideline, the Denali transect yields an exhumation history generally similar to that of the reference model (Fig. 10a). Differences in the a priori parameters include that the new inversion finds and uses an initial geothermal gradient of 25.2 °C km^{−1} (slightly higher than that of the reference model), an a priori erosion rate of 0.46 ± 0.23 km Myr^{−1} (slightly lower than that of the reference model) and a time interval length of 1.5 Ma. The combination of these a priori parameters results in a major increase in erosion rate to 0.55–0.6 km Myr^{−1} at 6 Ma, which is 1.5 Myr later than that of the reference model (7.5 Ma). The subtle differences from the reference model mainly result from the time interval length used in these models. Comparing the misfit values, the new model produces slightly better fits than the reference model, with the a posteriori misfit values of 1.81 and 0.11 for the observed age and geothermal data.
9.2 Himalayan Dharladar range transect
AFT and ZHe data from the Dharladar range in the northwestern Himalayas, reported in the publications by Deeken et al. (2011) and Thiede et al. (2017), are used as an example for regions of young cooling ages and fast exhumation. The samples were collected in an elevation range between 1.5 and 4.5 km, covering a topographic relief of 3 km within a spatial distance of ∼15 km on the hanging wall of the main central thrust of the Himalayan fold and thrust belt (Deeken et al., 2011; Thiede et al., 2017). AER slope regression of ZHe and AFT ages performed in Deeken et al. (2011) produced apparent erosion rates of ∼2.8 and ∼0.2 km Myr^{−1} for the time intervals 6.4–14.5 and 1.7–3.7 Ma, respectively, implying a potential increase in erosion rates at ∼3.7–6.4 Ma. Using geothermal gradients of 25–45 °C km^{−1}, timeaveraged erosion rates were estimated as 0.8–2.0 km Myr^{−1} since 3.7 Ma (Deeken et al., 2011).
The modeling of the Dharladar range data uses a modern geothermal gradient constraint of 45 ± 8 °C km^{−1} (Deeken et al., 2011). The relatively large uncertainty is assigned for the geothermal gradient because of the absence of direct geothermal measurements in the study area. Our exhumation inversion for the AER data using the stepwise modeling guideline yields relatively slow rates of 0.1–0.6 km Myr^{−1} and fast rates of 1.2–1.6 km Myr^{−1} before and after ∼3 Ma, respectively (Fig. 10b). The abrupt increase in exhumation rates at ∼3 Ma is generally consistent with the estimates from the slope regression results of Deeken et al. (2011). However, the inverted exhumation rates since 3 Ma are significantly lower than the estimation from the AER slope (∼2.8 km Myr^{−1}), which is likely due to the overestimation of exhumation of the AER slope due to topographic perturbation of isotherms. Such a perturbation is a function of exhumation rates: the higher the exhumation, the larger the perturbation (Glotzbach et al., 2015). The modeling yields a history of the geothermal gradient that gradually increases to a modern value of ∼46 °C km^{−1}, close to the input value (45 ± 8 °C km^{−1}).
9.3 KTB borehole
The KTB borehole yields a large amount of thermochronologic and geochronologic age data (Warnock and Zeitler, 1998; Stockli and Farley, 2004). Previous studies suggest the borehole is truncated by multiple faults, which offset the age–depth relationship (Wagner et al., 1997). Here we use the data at depths shallower than 1 km, where data are abundant and have linear relationship with depths.
The KTB apatite, zircon and titanite $\left(\mathrm{U}\text{}\mathrm{Th}\right)/\mathrm{He}$ (AHe, ZHe and THe) and AFT age data vary greatly between 85–50 Ma. These clustered ages have been interpreted as indicating a late Cretaceous phase of exhumation, followed by slow exhumation (Wagner et al., 1997; Stockli and Farley, 2004), as also shown by previous thermal history reconstructions based on K feldspar ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ data (Warnock and Zeitler, 1998).
Our modeling, using the AER data and a modern geothermal gradient of 27.5 ± 2.8 °C km^{−1} (Clauser et al., 1997), shows that elevated exhumation rates (0.1–0.13 km Myr^{−1}) between 80–50 Ma, followed by slower exhumation rates of ∼0.04 km Myr^{−1} (Fig. 10c), are similar to previous estimates (Wagner et al., 1997; Warnock and Zeitler, 1998; Stockli and Farley, 2004). Associated with changes in exhumation, the geothermal gradient gradually decreases from the peak values at 70–60 Ma to a value of ∼28 °C km^{−1} at the present day.
The a priori information has important effects on the inversion results using the leastsquares inversion method. Our study demonstrates the importance of geothermal gradient and the a priori exhumation rate in estimating the exhumation history from the thermochronology data. To take into account the geothermal data for the exhumation history inversion, we outlined a stepwise inversion method that first searches for the appropriate initial geothermal gradient, which is then used in the modeling searching for the a priori exhumation rate. Our modeling guideline produces an exhumation history and a geothermal gradient that provide reasonable fits for both the observed AER and modern geothermal data, as tested by datasets of both synthetic models and natural samples.
The code and data used in this work are available on Zenodo (https://doi.org/10.5281/zenodo.10839275, Tian, 2024).
The supplement related to this article is available online at: https://doi.org/10.5194/esurf124772024supplement.
YT: conceptualization, methodology, software, data curation, visualization, investigation, writing – original draft preparation. LP: visualization, writing – reviewing and editing. GZ and XY: writing – reviewing and editing.
The contact author has declared that none of the authors has any competing interests.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.
This study is funded by the National Natural Science Foundation of China (grant nos. 42172229, 42488201 and 41772211). Discussions with Jie Hu and Donglan Zeng are greatly appreciated. Comments and suggestions from Gilby Jepson and Christoph Glotzbach clarified many points of this work.
This research has been supported by the National Natural Science Foundation of China (grant nos. 42172229, 42488201 and 41772211).
This paper was edited by Richard Gloaguen and reviewed by Christoph Glotzbach and Gilby Jepson.
Anderson, J. L., Barth, A. P., Wooden, J. L., and Mazdab, F.: Thermometers and Thermobarometers in Granitic Systems, Rev. Mineral. Geochem., 69, 121–142, https://doi.org/10.2138/rmg.2008.69.4, 2008.
Armstrong, P. A.: Thermochronometers in Sedimentary Basins, Rev. Mineral. Geochem., 58, 499–525, https://doi.org/10.2138/rmg.2005.58.19, 2005.
Ault, A. K., Gautheron, C., and King, G. E.: Innovations in $\left(\mathrm{U}\text{}\mathrm{Th}\right)/\mathrm{He}$, Fission Track, and Trapped Charge Thermochronometry with Applications to Earthquakes, Weathering, SurfaceMantle Connections, and the Growth and Decay of Mountains, Tectonics, 38, 3705–3739, https://doi.org/10.1029/2018TC005312, 2019.
Batt, G. E. and Brandon, M. T.: Lateral thinking: 2D interpretation of thermochronology in convergent orogenic settings, Tectonophysics, 349, 185–201, 2002.
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, Bull. Geol. Soc. Am., 110, 985–1009, 1998.
Braun, J.: Quantifying the effect of recent relief changes on ageelevation relationships, Earth Planet. Sc. Lett., 200, 331–343, 2002.
Braun, J.: Pecube: a new finiteelement code to solve the 3D heat transport equation including the effects of a timevarying, finite amplitude surface topography, Comput. Geosci., 29, 787–794, 2003.
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.
Burbank, D. W., Blythe, A. E., Putkonen, J., PrattSitaula, B., Gabet, E., Oskin, M., Barros, A., and Ojha, T. P.: Decoupling of erosion and precipitation in the Himalayas, Nature, 426, 652–655, 2003.
Clauser, C., Giese, P., Huenges, E., Kohl, T., Lehmann, H., Rybach, L., Šafanda, J., Wilhelm, H., Windloff, K., and Zoth, G.: The thermal regime of the crystalline continental crust: Implications from the KTB, J. Geophys. Res.Sol. Ea., 102, 18417–18441, https://doi.org/10.1029/96JB03443, 1997.
Deeken, A., Thiede, R. C., Sobel, E. R., Hourigan, J. K., and Strecker, M. R.: Exhumational variability within the Himalaya of northwest India, Earth Planet. Sc. Lett., 305, 103–114, https://doi.org/10.1016/j.epsl.2011.02.045, 2011.
Dodson, M. H.: Closure temperature in cooling geochronological and petrological systems, Contrib. Mineral. Petr., 40, 259–274, 1973.
Ehlers, T. A.: Crustal Thermal Processes and the Interpretation of Thermochronometer Data, Rev. Mineral. Geochem., 58, 315–350, https://doi.org/10.2138/rmg.2005.58.12, 2005.
Ehlers, T. A. and Farley, K. A.: Apatite $\left(\mathrm{U}\text{}\mathrm{Th}\right)/\mathrm{He}$ thermochronometry: methods and applications to problems in tectonic and surface processes, Earth Planet. Sc. Lett., 206, 1–14, 2003.
Fitzgerald, P. G., Sandiford, M., Barrett, P. J., and Gleadow, A. J. W.: Asymmetric extension associated with uplift and subsidence in the Transantarctic Mountains and Ross Embayment, Earth Planet. Sc. Lett., 81, 67–78, https://doi.org/10.1016/0012821X(86)901019, 1986.
Fitzgerald, P. G., Sorkhabi, R. B., Redfield, T. F., and Stump, E.: Uplift and denudation of the central Alaska Range: A case study in the use of apatite fission track thermochronology to determine absolute uplift parameters, J. Geophys. Res.Sol. Ea., 100, 20175–20191, https://doi.org/10.1029/95JB02150, 1995.
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.
Fox, M., Herman, F., Kissling, E., and Willett, S. D.: Rapid exhumation in the Western Alps driven by slab detachment and glacial erosion, Geology, 43, 379–382, 2015.
Gallagher, K.: Transdimensional inverse thermal history modelling for quantitative thermochronology, J. Geophys. Res., 117, B02408, https://doi.org/10.1029/2011JB008825, 2012.
Glotzbach, C., Braun, J., and van der Beek, P.: A Fourier approach for estimating and correcting the topographic perturbation of lowtemperature thermochronological data, Tectonophysics, 649, 115–129, https://doi.org/10.1016/j.tecto.2015.03.005, 2015.
Herman, F., Seward, D., Valla, P. G., Carter, A., Kohn, B., Willett, S. D., and Ehlers, T. A.: Worldwide acceleration of mountain erosion under a cooling climate, Nature, 504, 423–426, 2013.
House, M. A., Wernicke, B. P., and Farley, K. A.: Paleogeomorphology of the Sierra Nevada, California, from $\left(\mathrm{U}\text{}\mathrm{Th}\right)/\mathrm{He}$ ages in apatite, Am. J. Sci., 301, 77–102, 2001.
Hu, S. B., Raza, A., Min, K., Kohn, B. P., Reiners, P. W., Ketcham, R. A., Wang, J. Y., and Gleadow, A. J. W.: Late Mesozoic and Cenozoic thermotectonic evolution along a transect from the north China craton through the Qinling orogen into the Yangtze craton, central China, Tectonics, 25, TC6009, https://doi.org/10.1029/2006TC001985, 2006.
Jiao, R., Herman, F., and Seward, D.: Late Cenozoic exhumation model of New Zealand: Impacts from tectonics and climate, EarthSci. Rev., 166, 286–298, https://doi.org/10.1016/j.earscirev.2017.01.003, 2017.
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.
Laslett, G., Green, P. F., Duddy, I., and Gleadow, A.: Thermal annealing of fission tracks in apatite 2. A quantitative analysis, Chem. Geol., 65, 1–13, 1987.
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.
McInnes, B. I. A., Evans, N. J., Fu, F. Q., and Garwin, S.: Application of Thermochronology to Hydrothermal Ore Deposits, Rev. Mineral. Geochem., 58, 467–498, https://doi.org/10.2138/rmg.2005.58.18, 2005.
Reiners, P. W. and Brandon, M. T.: Using thermochronology to understand orogenic erosion, Annu. Rev. Earth Pl. Sc., 34, 419–466, 2006.
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.
Stalder, N. F., Herman, F., Fellin, M. G., Coutand, I., Aguilar, G., Reiners, P. W., and Fox, M.: The relationships between tectonics, climate and exhumation in the Central Andes (18–36° S): Evidence from lowtemperature thermochronology, EarthSci. Rev., 210, 103276, https://doi.org/10.1016/j.earscirev.2020.103276, 2020.
Stockli, D. F. and Farley, K. A.: Empirical constraints on the titanite $\left(\mathrm{U}\text{}\mathrm{Th}\right)/\mathrm{He}$ partial retention zone from the KTB drill hole, Chem. Geol., 207, 223–236, https://doi.org/10.1016/j.chemgeo.2004.03.002, 2004.
Sutherland, R., Gurnis, M., Kamp, P. J. J., and House, M. A.: Regional exhumation history of brittle crust during subduction initiation, Fiordland, southwest New Zealand, and implications for thermochronologic sampling and analysis strategies, Geosphere, 5, 409–425, https://doi.org/10.1130/GES00225.1, 2009.
Tarantola, A.: Inverse Problem Theory and Methods for Model Parameter Estimation, SIAM, Philadelphia, https://doi.org/10.1137/1.9780898717921, 2005.
Thiede, R., Robert, X., Stübner, K., Dey, S., and Faruhn, J.: Sustained outofsequence shortening along a tectonically active segment of the Main Boundary thrust: The Dhauladhar Range in the northwestern Himalaya, Lithosphere, 9, 715–725, https://doi.org/10.1130/L630.1, 2017.
Tian, Y., Kohn, B. P., Hu, S., and Gleadow, A. J. W.: Synchronous fluvial response to surface uplift in the eastern Tibetan Plateau: Implications for crustal dynamics, Geophys. Res. Lett., 42, 29–35, https://doi.org/10.1002/2014GL062383, 2015.
Turcotte, D. and Schubert, G.: Geodynamics, Cambridge Univiversity Press, https://doi.org/10.1017/S0016756802217239, 2002.
Valla, P. G., van der Beek, P. A., and Braun, J.: Rethinking lowtemperature thermochronology data sampling strategies for quantification of denudation and relief histories: A case study in the French western Alps, Earth Planet. Sc. Lett., 307, 309–322, https://doi.org/10.1016/j.epsl.2011.05.003, 2011.
van der Beek, P. and Schildgen, T. F.: Short communication: age2exhume – a MATLAB/Python script to calculate steadystate vertical exhumation rates from thermochronometric ages and application to the Himalaya, Geochronology, 5, 35–49, https://doi.org/10.5194/gchron5352023, 2023.
van der Beek, P. A., Valla, P. G., Herman, F., Braun, J., Persano, C., Dobson, K. J., and Labrin, E.: Inversion of thermochronological ageelevation profiles to extract independent estimates of denudation and relief history – II: Application to the French Western Alps, Earth Planet. Sc. Lett., 296, 9–22, https://doi.org/10.1016/j.epsl.2010.04.032, 2010.
Vermeesch, P.: IsoplotR: A free and open toolbox for geochronology, Geosci. Front., 9, 1479–1493, https://doi.org/10.1016/j.gsf.2018.04.001, 2018.
Wagner, G. A., Coyle, D. A., Duyster, J., HenjesKunst, F., Peterek, A., Schröder, B., Stöckhert, B., Wemmer, K., Zulauf, G., Ahrendt, H., Bischoff, R., Hejl, E., Jacobs, J., Menzel, D., Lal, N., Van den Haute, P., Vercoutere, C., and Welzel, B.: PostVariscan thermal and tectonic evolution of the KTB site and its surroundings, J. Geophys. Res.Sol. Ea., 102, 18221–18232, https://doi.org/10.1029/96JB02565, 1997.
Warnock, A. C. and Zeitler, P. K.: ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ thermochronometry of Kfeldspar from the KTB borehole, Germany, Earth Planet. Sc. Lett., 158, 67–79, https://doi.org/10.1016/S0012821X(98)000442, 1998.
Whipp Jr., D. M., Ehlers, T. A., Blythe, A. E., Huntington, K. W., Hodges, K. V., and Burbank, D. W.: PlioQuaternary exhumation history of the central Nepalese Himalaya: 2. Thermokinematic and thermochronometer age prediction model, Tectonics, 26, TC3003, https://doi.org/10.1029/2006tc001991, 2007.
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.
Tian, Y.: An Matlab application for inverting rock exhumation from thermochronologic age data, Zenodo [data set] and [code], https://doi.org/10.5281/zenodo.10839275, 2024.
Zeitler, P., Meltzer, A., Koons, P., Craw, D., Hallet, B., Chamberlain, C., Kidd, W., Park, S., Seeber, L., Bishop, M., and Shroder, J. F.: Erosion, Himalayan Geodynamics, and the Geomorphology of Metamorphism, GSA Today, 11, 4–9, https://doi.org/10.1130/10525173(2001)011<0004:EHGATG>2.0.CO;2, 2001.
 Abstract
 Introduction
 Linear inversion method
 Closure depth and topographic correction
 Model evaluation
 The reference inverse model
 Dependence on model parameters and proposed solutions
 A new modeling guideline
 Synthetic models for testing the new modeling guideline
 Natural examples for testing the new modeling guideline
 Conclusion
 Code and data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Supplement
 Abstract
 Introduction
 Linear inversion method
 Closure depth and topographic correction
 Model evaluation
 The reference inverse model
 Dependence on model parameters and proposed solutions
 A new modeling guideline
 Synthetic models for testing the new modeling guideline
 Natural examples for testing the new modeling guideline
 Conclusion
 Code and data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Supplement