Solving crustal heat transfer for thermochronology using physics-informed neural networks

We present a deep-learning approach based on the physics-informed neural networks (PINNs) for estimating thermal evolution of the crust during tectonic uplift with a changing landscape. The approach approximates the temperature field of the crust with a deep neural network, which is trained by optimizing the heat advection–diffusion equation, assuming initial and boundary temperature conditions that follow a prescribed topographic history. From the trained neural network of temperature field and the prescribed velocity field, one can predict the temperature history of a given rock particle that can be used to compute the cooling ages of thermochronology. For the inverse problem, the forward model can be combined with a global optimization algorithm that minimizes the misfit between predicted and observed thermochronological data, in order to constrain unknown parameters in the rock uplift history or boundary conditions. We demonstrate the approach with solutions of one- and three-dimensional forward and inverse models of the crustal thermal evolution, which are consistent with results of the finite-element method. As an example, the three-dimensional model simulates the exhumation and post-orogenic topographic decay of the Dabie Shan, eastern China, whose post-orogenic evolution has been constrained by previous thermochronological data and models. This approach takes advantage of the computational power of machine learning algorithms, offering a valuable alternative to existing analytical and numerical methods, with great adaptability to diverse boundary conditions and easy integration with various optimization schemes.

1 Introduction

Thermochronology has been used extensively in Earth sciences for estimating ages and rates of geological and landscape evolution processes. To interpret the data, it must be understood how the thermal structure of the crust has changed over time and space. In the thermochronology community, analytical and numerical methods have advanced the analysis and interpretation of data by providing quantitative estimates of temperature fields, tectonic deformation, and landscape changes under various scenarios (e.g., Stüwe et al., 1994; Braun, 2003). However, these methods also involve strict assumptions and limitations which hinder their application to geological problems in a broader range of geological and landscape settings. The analytical methods are computationally efficient but almost only applicable to problems with constant uplift rates. They also require either the assumption of a steady-state temperature field (Stüwe et al., 1994; Brandon et al., 1998) or additional temperature data to quantify the transience of the field (Willett and Brandon, 2013). Numerical methods have the advantage of simulating the transient state of the thermal evolution (Braun, 2003; Braun et al., 2012), but implementing complex boundary conditions in numerical models is a significant challenge. Solving inverse problems with numerical models has a high computational demand which increases exponentially with the dimensionality of the space of unknown parameters.

Recently, the deep-learning method of physics-informed neural networks (PINNs) was proposed for solving partial differential equations, and its application has been demonstrated in many research fields including Earth sciences (e.g., He and Tartakovsky, 2021; Rasht‐Behesht et al., 2022). This mesh-free method benefits from the increasing capabilities of machine learning algorithms. Here we show that PINNs can be used to solve the heat transfer equation in the crust for given rock uplift and landscape histories. Using example models, we show that PINNs can provide estimates of the thermal models with good agreement with numerical solutions while allowing more flexibility in the model configuration. Moreover, PINNs also support the simultaneous optimization of the forward and inverse problems. Therefore, these advantages of PINNs indicate the potential of PINNs for solving problems with more intricate boundary conditions. Our examples are implemented using the TensorFlow 2 library (Abadi et al., 2016) and global optimizers from the SciPy library (Virtanen et al., 2020), but we expect that similar results can be achieved using other machine learning libraries and optimization schemes.

2 From thermochronological ages to exhumation rates: existing methods

In this section we provide a brief overview of the widely used methods for interpretation of thermochronological data, highlighting in particular how they deal with the subsurface thermal field and its influence on the observed data.

2.1 Age–elevation relationship and thermal history modelling

Some conventional methods interpret thermochronological data without quantifying the thermal structure in the crust. The age–elevation relationship (AER) is an approach in low-temperature thermochronology for estimating the rock exhumation rates (e.g., Fitzgerald et al., 1995; Fitzgerald and Malusà, 2019), using the apparent ages collected at different elevations on a vertical profile. Although this method is based on the concept of “closure temperature” (Dodson, 1973), it only assumes that the timer, i.e., a thermochronometer, starts at the same depth in the crust and does not require a calculation of the values of the closure temperature or its depth. Therefore, as rocks of different elevations, from higher to lower, should have passed the closure depth (i.e., depth of the closure temperature isotherm) consecutively, a positive correlation between the cooling ages and elevations of the samples is expected. For a steady-state cooling history (Willett and Brandon, 2002), this should result in a linear relationship between apparent ages and sample elevations, with a slope equal to the exhumation rate for the time span indicated by the thermochronological ages. However, as a change in the exhumation rate can perturb the isotherms, the AER approach is unreliable when the cooling rate is not constant or the thermal field is transient. In addition, application of the method is often limited by the difficulty of sampling on (near-)vertical profiles.

HeFTy (Ketcham, 2005) and QTQt (Gallagher, 2012) are widely used programs for estimating thermal history models using thermochronological data. The forward modelling functions in the programs can predict thermochronological data for a given thermal history path which is defined by a finite number of time–temperature points, based on simulating the noble gas production–diffusion or fission-track production–annealing kinetics. The forward model is often used together with a Monte Carlo sampling method to estimate the optimal thermal history that suits observed thermochronological data. When the time–temperature paths of multiple samples on a vertical profile are made available by modelling, they can be used to infer changes in the geothermal gradient over time: assuming that the elevation offset between samples on the profile has remained constant, the optimal temperature offset history can be used to provide constraints on the evolution of the paleo-geothermal gradient. Both QTQt and the recent version of HeFTy have the function to simultaneously sample the time–temperature offset space for multiple samples on a vertical profile; therefore the result will also yield an optimal history for the geothermal gradient (e.g., Jiao et al., 2014; Jepson et al., 2022). This approach requires no assumption of a steady-state or monotonic cooling history. However, neither of the two programs solves the heat transfer equations nor explicitly models the physical mechanisms that cause these changes. Moreover, at any time in the parameter space, the modelling procedure uses linear interpolation to determine the temperatures of samples according to their positions on the profile, which is not consistent with the crustal thermal profile in the case of rapid exhumation.

2.2 Calculating crustal temperatures

To ensure that the estimated exhumation rate is consistent with the physics of heat transfer, the thermal structure of the crust can be solved as a forward model using analytical methods or numerical models in one, two, or three dimensions. For example, the finite-difference code Tc1D (Whipp, 2022) computes the one-dimensional (1D) temperature model of the crust for various initial and boundary conditions, and it can be applied to predict thermochronological ages under different tectonic and geomorphic processes.

To elucidate the exhumation of the Olympic Mountains on the Cascadia margin, Brandon et al. (1998) designed an approach to compute the apatite fission-track age for 1D steady-state exhumation. For a given exhumation rate, the approach uses an analytical approximation (Stüwe et al., 1994) of the heat advection–diffusion equation to solve the thermal profile of the crust, which is then used to find the exhumation rate and the closure temperature (dependent on cooling rate; Dodson, 1973) through a numerical iteration. By using different kinetic parameters for track annealing or noble gas diffusion in the closure temperature equation, this approach has been updated to estimate exhumation rates from various thermochronometers (Reiners et al., 2003; Van Der Beek and Schildgen, 2023). To avoid the assumption of a steady-state scenario, Willett and Brandon (2013) proposed an alternative approach to solve the thermal profile of the Earth over a half-space domain, which requires constraints on the final geothermal gradient to quantify the transience of the thermal profile.

To solve the 1D heat transfer problem, the methods mentioned in the previous paragraph all assume a constant exhumation rate for the time period from the date indicated by a thermochronometer's cooling age to the present day. Fox et al. (2014) developed a linear inversion method to infer the variation in exhumation rates in time and space. To estimate the temporal variation in rates, the method also uses a 1D thermal model and the approximation of Dodson (1973) but considers the closure depth as a summation of rates over a finite number of time intervals. This discretization leads to more unknown parameters (i.e., exhumation rates) than data (i.e., cooling ages) and is thus an underdetermined problem. Therefore, the solution of the problem through linear inversion relies on some independent knowledge of the unknowns (e.g., a priori mean exhumation rate and the variance on this mean) to construct the covariance matrix of the exhumation rates. In the covariance matrix, Fox et al. (2014) also introduce a spatial correlation function in order to smooth the exhumation rates in space.

The 1D thermal models ignore the heat transport in horizontal directions. To capture the perturbation of isotherms by surface topography, instead of a sample's true elevation, the mean elevation filtered by the wavelength of topography (Stüwe et al., 1994; Mancktelow and Grasemann, 1997) has been used in the 1D models for solving the thermal profile of the crust (Brandon et al., 1998; Willett and Brandon, 2013; Fox et al., 2014; Van Der Beek and Schildgen, 2023). Two-dimensional (2D) numerical modelling has also been used to solve the temperature fields under undulating topography during exhumation. Mancktelow and Grasemann (1997) used a 2D finite-difference method to compute the transient isotherms to assess the influence of topography and rock uplift on estimated exhumation rates in active orogens. To constrain the exhumation pattern and history of the central Wasatch Mountains using thermochronological data, Ehlers et al. (2003) coupled a 2D velocity field and a thermal model to compute the temperature field evolution in the mountains; the thermal model is solved using a finite-difference method. To simulate the temperature field around hydrothermal systems, Luijendijk (2019) published a model to solve the subsurface heat flow in the 2D sections using the finite-element method, which has been used to predict apatite fission-track and (U-Th) $/$ He data.

In the low-temperature thermochronology community, simulation of the three-dimensional (3D) crustal thermal field evolution has been mostly performed using Pecube (Braun, 2003; Braun et al., 2012), a finite-element code written in Fortran. This code solves the 3D heat transfer equation in the crust under a prescribed velocity field and evolving topography, and it uses the solution to compute temperature histories of rocks and ages of various thermochronometers. The program has restricted formats for defining boundary conditions and integrating data and therefore is not straightforward to impose complex boundary conditions that are not already implemented. For example, the code uses a simplified kink-band model to compute the rock velocity field driven by the displacement on a single fault or an array of faults with the same strike; therefore it cannot simulate other kinematic models without significant modifications. Pecube can also operate in the inverse mode, which uses the neighbourhood algorithm (Sambridge, 1999a, b), an optimization method coded in C, to search for the best-fit values of parameters that specify the input tectonic or topographic scenarios. However, the search relies on solutions of a large number of forward problems and can thus be very demanding for computational power for high-dimensional problems with a large finite-element grid or many time steps.

3 Physics-informed neural networks

Neural networks perform complex and nonlinear data operations by mimicking the structure and function of biological neurons. In deep learning, multiple neural layers are combined as deep neural networks (DNNs) to decipher high-level information from the raw data. In scientific machine learning, deep learning is incorporated with fundamental physical laws, aiming to develop reliable, scalable, and physics-consistent machine learning models to facilitate new discoveries from scientific data. Physics-informed neural networks (PINNs) were recently proposed to solve forward and inverse problems involving partial differential equations (PDEs) (Raissi et al., 2019). In the context of PINNs, a fully connected neural network is generally employed to approximate the solution of the PDEs, e.g., the temperature, T(z,t), by taking the spatial and temporal coordinates (z and t) as the inputs. As demonstrated in Fig. 1, the neural network is composed of multiple hidden layers with trainable parameters and nonlinear activation functions. The parameters in the network can be learned using a loss function based on the governing equations (i.e., PDEs) and the boundary and initial conditions of the PDEs. For inverse problems, we add a misfit term representing the difference between the network output and the observation. Two separate optimizers are used for the loss function and data misfit. In the following sections we will demonstrate the setup of the neural networks case by case using numerical experiments.

Figure 1 Workflow for forward and inverse modelling of the 1D thermal profile for thermochronology. f_i and f_b are functions for the initial and boundary conditions of the model, respectively.

A key procedure in PINNs is to compute the derivatives of the PDEs. To address this issue, PINNs use automatic differentiation (AD) to represent all differential operators that exist in the PDEs. The AD calculates the derivatives of the outputs with respect to the network inputs based on the chain rule, which is different from numerical computations and avoids discretization and truncation errors. After training with the loss function (and data misfit), the trained network can be used to approximate the solution of PDEs at any location and any time.

PINNs have been applied in various disciplines such as fluid mechanics (Cai et al., 2021a; Raissi et al., 2020; Boster et al., 2023) and have demonstrated great potential in solving governing equations of physical laws in complex domains, especially in solving inverse problems. They have also attracted attention in the community of Earth sciences. For example, Waheed et al. (2021) used PINNs to solve the eikonal equation in 2D for predicting travel times of seismic waves in isotropic and anisotropic media. He and Tartakovsky (2021) used PINNs to solve the advection–dispersion and Darcy flow equations in 1D and 2D fields with spatially varying hydraulic conductivity. Rasht‐Behesht et al. (2022) used PINNs to solve the acoustic wave propagation and full waveform inversions under different and complex boundary conditions.

4 One-dimensional thermal profile: a synthetic model

4.1 Forward solution

4.1.1 Method and problem setup

In our first example, the 1D forward model describes the thermal profile evolution of a 30 km thick crust from 0 to 100 Myr, with the surface and basal temperatures fixed at 0 and 600 °C, respectively. The model exhumation rate remained constant at 0.05 km Myr⁻¹ from 0 to about 60 Myr, after which it increased to 0.6 km Myr⁻¹. The rapid exhumation continued for about 20 million years, and around 80 Myr the exhumation slowed down to 0.1 km Myr⁻¹ and remained at this rate for the rest of the model time (Fig. 2a).

Figure 2 Setup of the 1D thermal model of the crust. (a) Rock uplift function of the model. (b) Collocation (grey dots) and boundary condition points (coloured cross) of the PINN model.

The 1D heat advection–diffusion equation to solve is given by

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial T}{\partial t}}+u{\displaystyle \frac{\partial T}{\partial z}}-\mathit{\kappa }{\displaystyle \frac{{\partial }^{\mathrm{2}}T}{\partial z^{\mathrm{2}}}}=\mathrm{0},\end{array}$$

where T is the temperature ( °C), t is the model time (Myr), u is the rock uplift rate (km Myr⁻¹), κ is the thermal diffusivity here set at 25 km² Myr⁻¹, and z is the vertical position (km) relative to the model base (z= 0). In this model we assume no topographic change over time, so u is equal to the exhumation rate. We solve the model in the space defined by

$$\begin{array}{}\text{(2)}& t\in [\mathrm{0},\mathrm{100}]\mathrm{and}z\in [\mathrm{0},\mathrm{30}],\end{array}$$

for which Dirichlet boundary conditions are imposed as

$$\begin{array}{}\text{(3)}& T(t,\mathrm{0})=\mathrm{600}\mathrm{and}T(t,\mathrm{30})=\mathrm{0},\end{array}$$

with a linear gradient prescribed for the initial condition, i.e.,

$$\begin{array}{}\text{(4)}& T(\mathrm{0},z)=\mathrm{600}-\mathrm{20}z.\end{array}$$

We assume that u is time-dependent and impose its time dependence in terms of logistic functions as

$$\begin{array}{}\text{(5)}& u\left(t\right)=u_{\mathrm{0}}+{\displaystyle \frac{u_{\mathrm{1}}-u_{\mathrm{0}}}{\mathrm{1}+e^{-(t-t_{\mathrm{1}})/\mathrm{\Delta }t}}}+{\displaystyle \frac{u_{\mathrm{2}}-u_{\mathrm{1}}}{\mathrm{1}+e^{-(t-t_{\mathrm{2}})/\mathrm{\Delta }t}}},\end{array}$$

in which u₀=0.05 km Myr⁻¹, u₁=0.6 km Myr⁻¹, u₂=0.1 km Myr⁻¹, t₁=60 Myr, t₂=80 Myr, and Δt=1 Myr. Δt is an inertia factor that imposes a time length for the transition between two rock uplift rates.

To solve the forward model, we approximate T as a deep neural network (DNN) and such that T(t,z) can be learned by minimizing the loss function

$$\begin{array}{}\text{(6)}& \mathcal{L}={\mathcal{L}}_f+\mathit{\alpha }{\mathcal{L}}_{\mathrm{b}},\end{array}$$

where ℒ_f penalizes the residual of the governing equations (i.e., PDEs) and ℒ_b imposes the boundary conditions and initial conditions of the PDEs; α is a weighting parameter. In our experiments presented in the paper, assigning α to 1 ensures effective optimization of the PINN models. Here the initial condition is regarded as a special boundary condition. ℒ_f and ℒ_b are computed by

$$\begin{array}{}\text{(7)}& {\displaystyle }& {\displaystyle }{\mathcal{L}}_f={\displaystyle \frac{\mathrm{1}}{N_f}}\sum _{i=\mathrm{1}}^{N_f}\left|f\right(t_f^i,z_f^i){|}^{\mathrm{2}},\text{(8)}& {\displaystyle }& {\displaystyle }{\mathcal{L}}_{\mathrm{b}}={\displaystyle \frac{\mathrm{1}}{N_{\mathrm{b}}}}\sum _{i=\mathrm{1}}^{N_{\mathrm{b}}}|T_{\mathrm{b}}^i-T(t_{\mathrm{b}}^i,z_{\mathrm{b}}^i){|}^{\mathrm{2}},\end{array}$$

in which f(t,z) is defined for the left-hand side of the heat transfer equation (Eq. 1), $\mathit{\{}{t}_f^i,z_f^i{\mathit{\}}}_i^{N_f}$ depicts the N_f collocation points for f(t,z), and $\mathit{\{}{t}_{\mathrm{b}}^i,z_{\mathrm{b}}^i{\mathit{\}}}_i^{N_{\mathrm{b}}}$ depicts the N_b initial and boundary training data on T(t,z).

In our demonstration, we use the TensorFlow 2 library (Abadi et al., 2016) to build the DNN, which includes three hidden layers with 20 neurons in each layer and an additional transform layer to impose the boundary conditions as hard constraints (Lagaris et al., 1998; Lu et al., 2021). The hyperbolic tangent function (tanh) is used as the activation function, and the Adam optimizer (Kingma and Ba, 2015) is used for minimizing the loss function (Eq. 6). The model is set up with 2000 randomly sampled collocation points. Another 500 points are used to impose the initial and boundary conditions, including 100, 200, and 200 points for the initial thermal structure, the surface temperature, and the basal temperature, respectively (Fig. 2b). For the solution presented here, the DNN is optimized through 50 000 iterations. The initial learning rate is set at 0.001, which decays at a factor of 0.9 for every 1000 steps; i.e., at the ith step the learning rate is $\mathrm{0.001}\times {\mathrm{0.9}}^{i/\mathrm{1000}}$.

As a comparison, we also solve the same problem using the finite-element method (FEM), which provides estimates of the heat transfer equations by discretizing the time–space domain using 101 points on the crustal profile and a time step of 1 Myr. To compare the final solutions of the PINN and FEM, we compute the square and infinity norms of the misfits between the predictions of the two methods as

$$\begin{array}{}\text{(9)}& {\displaystyle }& {\displaystyle }{l}^{\mathrm{2}}=\sqrt{\sum _{i=\mathrm{1}}^n(T_{\mathrm{pinn}}^i-T_{\mathrm{fem}}^i{)}^{\mathrm{2}}}\mathrm{and}\text{(10)}& {\displaystyle }& {\displaystyle }{l}^{\mathrm{\infty }}=\underset{i=\mathrm{1}}{\overset{n}{max}}\left(\right|T_{\mathrm{pinn}}^i-T_{\mathrm{fem}}^i\left|\right),\end{array}$$

in which $T_{\mathrm{pinn}}^i$ and $T_{\mathrm{fem}}^i$ are the temperatures predicted by the PINN and FEM at locations of each element used in FEM and n is the total number of elements. Based on the temperature histories estimated, we also predict the thermochronological ages of the present-day surface rocks using an empirical model of the fission-track annealing (Ketcham et al., 2007) and the noble gas diffusion kinetics (summarized by Reiners and Brandon, 2006) solved by a finite-difference method (Braun et al., 2006). Then the percent errors in ages between the predictions are calculated as

$$\begin{array}{}\text{(11)}& \mathit{\delta }=\left|{\displaystyle \frac{a_{\mathrm{pinn}}-a_{\mathrm{fem}}}{a_{\mathrm{fem}}}}\right|\times \mathrm{100}\mathit{\%},\end{array}$$

in which a_pinn and a_fem are ages calculated by the thermal histories estimated using the PINN and FEM, respectively.

4.1.2 Results

In the first 3000 iterations, the model loss ℒ has decreased rapidly from > 10 000 to < 10 (Fig. 3). After that, the ℒ_b experiences some fluctuations, but ℒ continues to decrease smoothly and becomes more stable after 30 000 iterations. Figure 4 shows the thermal profiles at different model times predicted by the PINN at various stages of the training, in comparison to the FEM solutions. Between the final PINN and the FEM solutions at model times of 40, 70, and 100 Myr, the l² calculated for 101 elements is 3.8, 51.0, and 30.0 °C, respectively, and the l^∞ is 0.6, 11.4, and 3.5 °C, respectively, indicating reasonable consistency between the two methods (Fig. 4). For the rock at the present-day surface, the PINN also predicts a time–temperature path very consistent with the FEM result (Fig. 5a), with a temperature difference < 20  °C for any point on the path and < 10  °C for the most part of the history (Fig. 5b). For the same rock, the predicted difference in cooling ages between the two modelling methods is < 7 % for apatite (U-Th)/He and apatite fission-track and < 3 % for other thermochronometers (Fig. 5d).

Figure 3 Loss curves of the 1D forward model. ℒ_f and ℒ_b are losses from the heat transfer equation and the boundary conditions, respectively, whereas ℒ is the total loss.

Figure 4 Predicted thermal profiles of the crust at different model times. Coloured lines depict the PINN solutions after 2000, 10 000, 20 000, and 50 000 training iterations of training, respectively. The dashed line indicates the solution using the FEM solution.

Figure 5 Comparison of the thermal histories and ages computed from the PINN and FEM solutions of the forward model. (a) Predicted time–temperature paths and cooling ages for rock sample ended at the surface. Coloured lines represent the PINN predictions after 2000, 10 000, 20 000, and 50 000 iterations, respectively. The dashed line indicates the FEM solution. (b) Temperature discrepancy between the PINN and the FEM solutions. Line colour depicts the iteration numbers performed by the PINN, as shown in the legend of panel (a). (c) Thermochronological ages computed from the thermal models solved by the PINN and FEM. (d) Percent error between the ages predicted by the PINN and FEM. AHe represents apatite (U-Th) $/$ He, AFT represents apatite fission track, ZHe represents zircon (U-Th) $/$ He, KsAr represents K-feldspar ⁴⁰Ar $/$ ³⁹Ar, BiAr represents biotite ⁴⁰Ar $/$ ³⁹Ar, and MuAr represents muscovite ⁴⁰Ar $/$ ³⁹Ar.

4.2 Inverse problem

4.2.1 Method and model setup

We transform the forward model to an inverse problem by considering u₀, u₁, u₂, t₁, and t₂ in the rock uplift function (Eq. 5) as unknown parameters. All other parameters defining the thermal model remain the same as in the forward problem (Sect. 4.1). Implementation of the DNN follows a similar setup to that of a forward model, except that the five unknown parameters are updated and optimized during the iterations by reducing the misfit between predicted and observed ages. The age misfit function is defined by

$$\begin{array}{}\text{(12)}& \mathit{\varphi }={\displaystyle \frac{\mathrm{1}}{N_{\mathrm{a}}}}\sum _{i=\mathrm{1}}^{N_{\mathrm{a}}}{\left({\displaystyle \frac{a_{\mathrm{p}}^i-a_{\mathrm{o}}^i}{{\mathit{\sigma }}_{\mathrm{o}}^i}}\right)}^{\mathrm{2}},\end{array}$$

where $a_{\mathrm{p}}^i$ is the predicted age, $a_{\mathrm{o}}^i$ is the observed age, ${\mathit{\sigma }}_{\mathrm{o}}^i$ is the uncertainty on the observation, and N_a is the total number of observed ages. For the problem demonstrated here, synthetic age data predicted from the temperature history computed using the FEM (Fig. 5) are used as data, and a 10 % uncertainty (standard deviation) is assumed for each data point. For calculation of the synthetic observed data, 0.05 km Myr⁻¹, 0.6 km Myr⁻¹, 0.1 km Myr⁻¹, 60 Myr, and 80 Myr are used as “true” values of u₀, u₁, u₂, t₁, and t₂ in the rock uplift function (Eq. 5), respectively.

The Adam algorithm is used to fit $\mathit{\{}{t}_f^i,z_f^i{\mathit{\}}}_i^{N_f}$ for loss functions derived from the physical law and boundary conditions. To circumvent local minima, the search for optimal values of the unknown parameters in the rock uplift function is conducted using a modified Lipschitzian approach, DIRECT (Jones et al., 1993; Jones and Martins, 2021), which is a global, derivative-free optimization algorithm implemented in the SciPy library (Virtanen et al., 2020). DIRECT operates by progressively partitioning the search space into smaller hyperrectangles and evaluating the objective function at the centre of each. In the subsequent iteration, hyperrectangles with optimal solutions are selected and subdivided for further evaluation. This process is repeated to locate the global minimum in the search space.

The configuration of the DIRECT algorithm requires adjusting a single parameter, ϵ, which defines the minimal acceptable improvement in the function value from the current best solution to the next potentially optimal hyperrectangle for division. A larger ϵ guides the search to explore a broader domain, whereas a smaller ϵ directs a more local exploitation. To optimize the memory usage in our application, we use a strategy in which, after a certain number of iterations, if a parameter value of the optimal solution has no significant change (less than 5 % within the search domain), the DIRECT search recommences within a narrowed search space, reducing the domain for this parameter by 10 %.

More specifically, we define the initial search space for the parameters u₀, u₁, u₂, t₁, and t₂ within ranges of 0–2 km Myr⁻¹, 0–2 km Myr⁻¹, 0–2 km Myr⁻¹, 0–100 Myr, and 0–100 Myr, respectively. We start the Adam optimization of the DNN model by assigning the midpoint values of the search space to the unknown parameters. After n_adam iterations, we use the refined DNN model to explore the space of unknown parameters using the DIRECT approach over n_direct iterations. This process is repeated for multiple cycles to optimize the DNN model and seek the optimum parameter values. In the presented example, the initial n_adam is set to 100 and increased by 50 % for each DIRECT cycle, up to a ceiling of 10 000. The learning rate for Adam follows the same schedule used for the forward model; i.e., at the ith iteration the learning rate is $\mathrm{0.001}\times {\mathrm{0.9}}^{i/\mathrm{1000}}$. For each cycle, ϵ is maintained at 0.1, and the maximum number of DIRECT iterations, n_direct, is limited to 30 with the maximum number of function evaluated fixed at 1000. Cumulatively, the optimization encompasses > 55 000 Adam iterations and > 290 DIRECT iterations, based on > 15 000 age misfit function evaluations.

We also address the inverse problem using a Markov chain Monte Carlo (MCMC) approach. In the MCMC sampling process, the space of unknown parameters is searched iteratively, and at each iteration the parameter values with higher likelihood are accepted to the Markov chain. After sampling, the chain is used to estimate the probability distribution of the parameters. To evaluate the parameter values, the forward thermal model is simulated with the FEM based on 101 elements on the profile and a time step of 1 Myr. The log-likelihood of a forward model is calculated by comparing the predicted thermochronological ages with the observation, as

$$\begin{array}{}\text{(13)}& \log \left(L\right)=-\sum _{i=\mathrm{1}}^{N_{\mathrm{a}}}\left({\displaystyle \frac{\ln \left(\mathrm{2}\mathit{\pi }\right)}{\mathrm{2}}}+\ln \left({\mathit{\sigma }}_{\mathrm{o}}^i\right)+\mathrm{0.5}{\left({\displaystyle \frac{a_{\mathrm{p}}^i-a_{\mathrm{o}}^i}{{\mathit{\sigma }}_{\mathrm{o}}^i}}\right)}^{\mathrm{2}}\right).\end{array}$$

In our analysis, the space of five unknown parameters is explored by an affine invariant ensemble sampler (Goodman and Weare, 2010) within predefined domains. The search is conducted by 20 walkers, with each starting from a random location in the parameter space and moving for 50 000 steps. Ultimately, the generated chain comprises 1 000 000 combinations of parameter values from the search space.

4.2.2 Results

The optimization performance is shown as the evolution of the loss from the physical law and boundary conditions (Fig. 6a), as well as misfits between the predicted and observed thermochronological data (Fig. 6b). The loss function curve, despite minor fluctuations between DIRECT cycles (Fig. 6a), shows a reduction pattern similar to that of the forward model (Fig. 3). Figure 7 shows the sampled values used to evaluate the age misfit function. During the DIRECT search, the minimum age misfits decrease progressively until around 200 DIRECT iterations (Fig. 6b), after which the age misfits and parameter values of the optimal solution stabilize (Fig. 7). Finally, the searches of all five parameters are narrowed to ranges centred around the “true” values, except a relatively larger uncertainty on the estimated rate of the rapid uplift, u₁.

Figure 6 Optimization process of the 1D inverse model. (a) Loss function curve of the optimization. (b) Age misfits evaluated by the DIRECT search.

Figure 7 Parameter values sampled by the DIRECT method coupled with the PINN model for the 1D inverse model. Dots are colour-coded according to evaluated function value (age misfit). Dashed horizontal lines indicate the “true” parameter values.

The MCMC sampling coupled with the FEM model has also identified optimal values for the unknown parameters. Among the 20 walkers deployed, 3 failed to reach equilibrium at the end of the sampling, whereas the remaining walkers achieved stationary states after around 17 000 iterations (Fig. 8). Similar to the result of the DIRECT search, a large uncertainty exists for the estimate of u₁ (Fig. 8b).

Figure 9 presents predictions of the optimized PINN model and the “best-fit” FEM model, the model with the maximum likelihood from the MCMC chain. For both models, the predicted cooling paths of the rock at the current surface align closely with the “true” thermal history (Fig. 9a). The temperature discrepancy at any given time between the PINN solution and the “true” thermal history is < 17 °C and typically < 8 °C (Fig. 9b), whereas that between the “best-fit” FEM model and the “true” thermal history is < 25 °C. Moreover, the optimized PINN and “best-fit” FEM models predict similar thermochronological ages for the surface sample, both consistent with the data computed from the synthetic “true” history (Fig. 9c), with the maximum difference at 5.3 % and 3.7 %, respectively (Fig. 9d).

5 Three-dimensional temperature field: post-orogenic decay of the Dabie Shan

In this section, we use an example from the Dabie Shan (Fig. 10) to demonstrate the application of PINNs for solving the 3D temperature field under changing topography. Compared to the 1D model, the 3D model requires configuration of more complex boundary conditions and is computationally more intensive. The Dabie Shan in eastern China is a Mesozoic mountain range that underwent rapid exhumation during the Triassic–Jurassic orogeny along the boundary between the North and South China blocks (Nie et al., 1994; Hacker et al., 1995; Ratschbacher et al., 2006). Based on apatite and zircon (U-Th) $/$ He and apatite fission-track ages, Reiners et al. (2003) suggested that since the Late Cretaceous, the post-orogenic exhumation of the Dabie Shan has been slow and could be constant at a rate of 0.01–0.06 km Myr⁻¹. By investigating the same dataset using numerical methods, Braun and Robert (2005) estimated the post-orogenic exhumation and topographic evolution of the Dabie Shan and inferred that, since the Late Cretaceous, the mountain range has been exhumed at a mean rate of 0.01–0.04 km Myr⁻¹ and that its topographic relief has reduced by a factor of 2.5–4.5 during the last 60–80 million years. Therefore, thermochronological data provide good constraints on the post-orogenic exhumation of the Dabie Shan, which can be simulated by a simple model with constant exhumation rates, making the mountain range an ideal natural laboratory for testing new thermochronological modelling methods (e.g., Fox et al., 2014).

Figure 8 Parameter values sampled by the MCMC method coupled with the FEM model for the 1D inverse model. Dots are colour-coded according to the log-likelihood of the forward models. Dashed horizontal lines indicate the “true” parameter values.

Figure 9 Thermal histories and ages predicted by the optimized PINN model and the “best-fit” FEM model from the MCMC sampling chain. (a) Predicted time–temperature paths and cooling ages for a surface rock sample at the end of the model history. (b) Temperature discrepancy between the predicted time–temperature paths and the “true” model. (c) Thermochronological ages calculated from the PINN solution and the “best-fit” FEM model, compared to the synthetic data from the “true” thermal history. (d) Percent error between the predicted ages and synthetic data.

Figure 10 (a) Location of the Dabie Shan, eastern China. (b) Topography of the Dabie Shan and locations of thermochronological data.

5.1 Forward solution

5.1.1 Method and model setup

Assuming that the rock motion is in the vertical direction, the 3D heat transfer equation is given by

$$\begin{array}{}\text{(14)}& {\displaystyle \frac{\partial T}{\partial t}}+u{\displaystyle \frac{\partial T}{\partial z}}-\mathit{\kappa }({\displaystyle \frac{{\partial }^{\mathrm{2}}T}{\partial x^{\mathrm{2}}}}+{\displaystyle \frac{{\partial }^{\mathrm{2}}T}{\partial y^{\mathrm{2}}}}+{\displaystyle \frac{{\partial }^{\mathrm{2}}T}{\partial z^{\mathrm{2}}}})=\mathrm{0},\end{array}$$

where x, y, and z are spatial coordinates and other parameters are the same as in Eq. 1. We focus on the last phase of the Mesozoic orogeny and the post-orogenic exhumation history of the mountain and therefore solve its thermal history since 150 Ma. The model space is defined as

$$\begin{array}{}\text{(15)}& \begin{array}{rl}& t\in [\mathrm{0},\mathrm{150}],x\in [\mathrm{0},\mathrm{159}],y\in [\mathrm{0},\mathrm{199}],\\ & \mathrm{and}z\in [\mathrm{0},\mathrm{36.5}],\end{array}\end{array}$$

in which the range of the z domain is configured to include a crustal thickness (h_c=30 km) below sea level and the potential maximum elevation of the mountain range during the model run (h_s). The temperature at the base of the model is fixed at 600 °C. The Earth surface temperature is set to 15 °C at sea level and calibrated according to the elevation using an atmospheric lapse rate at −5 °C km⁻¹. Therefore, the boundary conditions are imposed as

$$\begin{array}{}\text{(16)}& \begin{array}{rl}& T(t,x,y,\mathrm{0})=\mathrm{600}\mathrm{and}\\ & T(t,x,y,z_{\mathrm{s}}(t,x,y\left)\right)=\mathrm{15}-\mathrm{5}{h}_{\mathrm{s}}(t,x,y),\end{array}\end{array}$$

in which z_s and h_s are the vertical coordinate and surface elevation, respectively, with a relationship

$$\begin{array}{}\text{(17)}& z_{\mathrm{s}}=h_{\mathrm{s}}+\mathrm{30}.\end{array}$$

To simulate the surface relief change, we follow Braun and Robert (2005) and define a scenario in which between 0 to 70 Myr (equivalent to geological time 150 and 80 Ma) the surface relief was 4 times as high as that of today, and after that the topography has gradually decreased in a linear fashion towards the current level. The topography is imposed as

$$\begin{array}{}\text{(18)}& h_{\mathrm{s}}(t,x,y)=w\left(t\right)h_{\mathrm{s}}(\mathrm{0},x,y),\end{array}$$

in which w is a time-dependent amplification factor for topography (Fig. 11b), prescribed by

$$\begin{array}{}\text{(19)}& w\left(t\right)=\left\{\begin{array}{l}{w}_{\mathrm{0}},t\in [\mathrm{0},t_{\mathrm{d}}]\\ \mathrm{1}+(w_{\mathrm{0}}-\mathrm{1}){\displaystyle \frac{\mathrm{150}-t}{\mathrm{150}-t_{\mathrm{d}}}},t\in (t_{\mathrm{d}},\mathrm{150}],\end{array}\right.\end{array}$$

where w₀, set at 4, is the amplification of topography at the beginning of model (t=0) and t_d, set at 70 Myr, is the model time when the post-orogenic decay of the topography starts. The present-day topography ($h_{\mathrm{s}}(\mathrm{0},x,y)$) is extracted from the global elevation model (GEBCO 2014) and resampled to a 2 km resolution.

Figure 11 Inputs of the 3D thermal evolution model of the Dabie Shan. (a) Rock uplift model. (b) Topographic amplification model.

The initial thermal field is imposed using a linear interpolation of temperatures between the model base and earth surface, i.e.,

$$\begin{array}{}\text{(20)}& \begin{array}{rl}T(\mathrm{0},x,y,z)& =T(\mathrm{0},x,y,\mathrm{0})\\ & -{\displaystyle \frac{T(\mathrm{0},x,y,\mathrm{0})-T(\mathrm{0},x,y,z_{\mathrm{s}}(\mathrm{0},x,y\left)\right)}{z_{\mathrm{s}}(\mathrm{0},x,y)}}z.\end{array}\end{array}$$

Based on conclusions of Reiners et al. (2003) and Braun and Robert (2005), we assume a spatially uniform rock uplift model and define it using a piecewise function (Fig. 11a),

$$\begin{array}{}\text{(21)}& u\left(t\right)=\left\{\begin{array}{l}{u}_{\mathrm{0}},t\in [\mathrm{0},t_{\mathrm{1}}]\\ u_{\mathrm{1}},t\in (t_{\mathrm{1}},\mathrm{150}],\end{array}\right.\end{array}$$

where u₀=0.6 km Myr⁻¹, u₁=0.05 km Myr⁻¹, and t₁=40 Myr. Therefore, the rock uplift remained constant at 0.6 km Myr⁻¹ between 0 and 40 Myr, and then the rate decreased to 0.05 km Myr⁻¹ and remained constant until the end of model.

Similarly to the solution of the 1D model, we approximate T as a DNN and train it using the same loss function (Eq. 6). The loss functions from 3D heat transfer law and boundary (and initial) conditions become

$$\begin{array}{}\text{(22)}& {\displaystyle }& {\displaystyle }{\mathcal{L}}_f={\displaystyle \frac{\mathrm{1}}{N_f}}\sum _{i=\mathrm{1}}^{N_f}\left|f\right(t_f^i,x_f^i,y_f^i,z_f^i){|}^{\mathrm{2}},\text{(23)}& {\displaystyle }& {\displaystyle }{\mathcal{L}}_{\mathrm{b}}={\displaystyle \frac{\mathrm{1}}{N_{\mathrm{b}}}}\sum _{i=\mathrm{1}}^{N_{\mathrm{b}}}|T_{\mathrm{b}}^i-T(t_{\mathrm{b}}^i,x_{\mathrm{b}}^i,y_{\mathrm{b}}^i,z_{\mathrm{b}}^i){|}^{\mathrm{2}},\end{array}$$

in which $f(t,x,y,z)$ stands for the left-hand side of the 3D heat transfer equation (Eq. 14).

The DNN used for the 3D model comprises 16 hidden layers, each with 20 neurons, and an additional layer to transform the outputs. The tanh activation function and Adam optimizer are employed, similarly to those used for the 1D model. An exponential decay function, $\mathrm{0.001}\times {\mathrm{0.9}}^i/\mathrm{10000}$, is used to define the learning rate. We selected 2000 random points to define the initial thermal field (Fig. 12a), and we used 2000 and 10 000 points to impose the temperatures at the model base and on the surface topography, respectively (Fig. 12b). We configured the model with 100 000 collocation points (Fig. 12c), half of which are allocated randomly within the orogenic phase ($t\in [\mathrm{0},\mathrm{40}]$) and the remainder during the post-orogenic phase ($t\in (\mathrm{40},\mathrm{150}]$).

Figure 12 The initial and boundary conditions and collocation points for the 3D thermal model of the Dabie Shan. Control points are projected onto the $\mathit{\{}x,y,z\mathit{\}}$ space. (a) The initial thermal field. (b) Boundary temperatures at the surface and base of the model. (c) Collocation points.

To evaluate the accuracy of the solution of the PINN, we also computed the thermal field evolution of the Dabie Shan using the FEM code, Pecube (Braun, 2003). The FEM model is set up on a $\mathrm{120}\times \mathrm{120}\times \mathrm{31}$ ($x\times y\times z$) grid, with the thermal field solved at a time step of 1 Myr. The rock uplift history and boundary conditions are identical to those in the PINN model.

5.1.2 Results

Figure 13 shows the loss curves through the optimization process, which present a progressive decrease until stabilization after approximately 300 000 iterations. Figure 14 presents the temperature profiles on three transects across the model at 40 and 150 Myr computed by the PINN, which reveal significant perturbations of the crustal thermal structures caused by the rock uplift and surface topography, resembling the patterns of the FEM solution. At 40 Myr, at the majority of locations (90 %), the temperature differences predicted by the PINN and FEM are < 10 °C. On these transects, the l² and l^∞ of the temperature discrepancies between the solutions of the two methods are 240–426 and 26–35 °C, respectively; these values are calculated from the solutions at the 3131 element locations used in the FEM (Fig. 15). At 150 Myr when the rock uplift rate has decreased and the topographic relief is lower, the discrepancies between the PINN and FEM solutions are significantly reduced, with the l² and l^∞ decreasing to 70–112 and 8–10 °C, respectively.

Figure 13 Loss curves for the 3D forward model.

Figure 14 Predicted isotherms on three transects across the Dabie Shan at model times of 40 and 150 Myr. Solutions using the PINN and FEM are compared. Locations of the transects are shown in Fig. 16.

Figure 15 Discrepancies in temperature solutions between the PINN and FEM on three transects across the Dabie Shan at model times of 40 and 150 Myr. Locations of the transects are shown in Fig. 16.

Figure 16 Maps of the Dabie Shan showing the predicted age misfits between the PINN and FEM solutions. A–A^′, B–B^′, and C–C^′ indicate locations of the transects in Fig. 14.

Figure 16 shows the misfits in the AHe and AFT ages for the present-day surface rocks at 42 000 locations, which are calculated from their thermal histories predicted by the PINN and FEM. For both thermochronometers, the predicted age misfits at all locations are < 20 %. At most locations (90 %), the misfits between AHe ages predicted by the two methods are < 6 %, and those between AFT ages are < 10 %. Compared to AHe ages, the higher misfits in predicted AFT ages reflect larger discrepancies in the temperature solutions between the PINN and FEM for the periods of high rock uplift rate and significant topographic relief (Fig. 15). At the sites with observed data, the misfits in predicted AHe, AFT, and ZHe ages between the PINN and FEM solutions are < 13 %, < 15 %, and < 2 %, respectively, which generally fall within the error margins or clustering ranges of the measured data (Fig. 17).

Figure 17 Predicted and observed thermochronological ages versus sample elevations.

5.2 Inverse problem

5.2.1 Model setup

To test the effectiveness of PINNs in constraining the exhumation models using real data, we formulate an inverse problem, aiming to estimate the post-orogenic uplift and topographic evolution of the Dabie Shan using published low-temperature thermochronological ages. In our analysis, we seek to constrain the orogenic rock uplift rate (u₀), the ending time of the rapid uplift (t₁), and the rock uplift rate during the post-orogenic period (u₁), as described in Eq. (21). In addition, we estimate the amplification factor for the topography during the orogenic phase (w₀) and the onset time for topographic decay (t_d) as described in Eq. (19). We confine the search for these five parameters, u₀, u₁, w₀, t₁, and t_d, within the ranges of 0–1 km Myr⁻¹, 0–1 km Myr⁻¹, 0–6, 0–50 Myr, and 50–100 Myr, respectively. To guide the optimization process, thermochronological age data reported by Reiners et al. (2003) (Fig. 10b) are used, considering only the ages < 150 Ma.

Similarly to the 1D model inversion, we use both the Adam and DIRECT algorithms to optimize the DNN model and explore the search space. For the Adam iterations, we set the learning rate using the same exponential decay function applied in the forward 3D model. In the DIRECT search, we use the “restart-and-zoom-in” strategy akin to that employed in the forward 3D model but increase the number of cycles to 40. Each DIRECT cycle includes up to 30 iterations and a maximum of 1000 function evaluations. Altogether, the inversion process comprises > 500 000 Adam iterations and > 1000 DIRECT iterations, evaluating > 40 000 age misfit functions.

As a reference, we also explore the multi-dimensional parameter space using the neighbourhood algorithm (NA) (Sambridge, 1999a), in which the forward model is computed using Pecube (Braun, 2003). It is worth noting that, compared to the modelling employed by Braun and Robert (2005) in the same region, we do not consider the isostatic rebound in response to the exhumation; thus the result may differ slightly. To reduce the computational demands compared to the model used in the forward problem (Sect. 5.1), we reduce the FEM grid size to $\mathrm{80}\times \mathrm{80}\times \mathrm{31}$ by increasing the horizontal spacing between elements. The NA divides the parameter space into Voronoi cells and evaluates the model misfit according to the parameter values at the cell centres. The subsequent iteration generates new samples within a subset of cells containing the best-fit models from the previous iteration. This process is repeated to find the optimal parameter values that fit the observed data. Here the NA search for parameters u₀, u₁, w₀, t₁, and t_d are within the same domains as predefined for the inversion of the PINN model. The search comprises 150 iterations with each running forward models in 60 cells, out of which the 75 % best cells are resampled for the next iteration. In total, the NA inversion includes 9060 forward Pecube runs.

5.2.2 Results

The training process of the PINN is shown in Fig. 18 as the evolution of the thermal model loss and the misfits in predicted ages. The parameter values evaluated by the DIRECT search are shown in Fig. 19. Notably, the loss function curve displays more fluctuation compared to the forward model, but it maintains a general decreasing trend. This fluctuation reduces as the learning rate decreases (Fig. 18a). Upon the stabilization of the loss function after approximately 200 000 Adam and 500 DIRECT iterations, the predicted minimum age misfits have a declining pattern for around 250 DIRECT iterations and show minimal changes thereafter (Fig. 18b). During the parameter search, convergence mostly occurs within the initial 100 and then from 500 to 1000 DIRECT iterations (Fig. 19). Finally, estimates for the parameters u₀, u₁, w₀ , t₁, and t_d are confined to the ranges of 0.7–0.95 km Myr⁻¹, 0–0.18 km Myr⁻¹, 2–5.5, 21–45 Myr, and 74–91 Myr, respectively.

Figure 18 Optimization process of the 3D inverse model of the Dabie Shan. (a) Loss function curve of the optimization. (b) Age misfits evaluated by the DIRECT search.

Figure 19 Evaluated parameter values for the 3D inverse model of the Dabie Shan using the DIRECT algorithm (Jones et al., 1993) and the PINN model.

The NA search based on forward Pecube models exhibits patterns comparable to those observed in the optimization of the PINN model. Specifically, the search for u₀ also converges rapidly to a range of 0–0.1 km Myr⁻¹, while the other four sampled parameters show significant convergence only after approximately 80 iterations. In the NA search the parameters u₀, u₁, w₀, t₁, and t_d are ultimately confined within the ranges of 0.8–1 km Myr⁻¹, 0.01–0.13 km Myr⁻¹, 4.1–5.1, 33–44 Myr, and 88–97 Myr, respectively (Fig. 20). Except for the parameter w₀ being confined to a narrower range by the NA, the search results by the two inverse models overlap over significant ranges.

Figure 20 Evaluated parameter values for the 3D inverse model of the Dabie Shan using the neighbourhood algorithm (Sambridge, 1999a) coupled with the Pecube model.

6 Discussion

6.1 Comparison between the PINN and numerical method

We have demonstrated that, in the study of thermochronology, PINNs are comparably as effective as conventional numerical methods. In both the 1D and 3D forward models we presented, the mean temperature difference across all points at any given time between the PINN and FEM solutions is < 0.5 °C, with the maximum difference at specific points ranging between 10 and 30 °C, which appear to occur during periods of rapid rock uplift and exhumation (Figs. 5b and 15). This level of discrepancy is minor compared to the typical uncertainties in thermal histories derived from thermochronological data. For 90 % of the surface samples in our examples, temperature discrepancies between the PINN and FEM solutions resulted in < 10 % misfits in the predicted ages. Larger misfits up to 20 % were found in cooling ages of samples that record rapid rock exhumation in areas with high topographic relief (Figs. 10b and 17b). These misfits in age estimates are at comparable magnitude to inherent uncertainties in measured age data, and they can be less significant than the data dispersion caused by the variability in the kinetics of noble gas diffusion or fission-track annealing.

For simple forward models, conventional numerical methods do not require optimization, thus having an advantage in the computational efficiency over PINNs. For instance, solving the 1D forward model in Sect. 4.1 using the FEM required a few milliseconds on a modern CPU (Intel Xeon Gold 5420+), whereas training the PINN model with 80 neurons for 50 000 iterations on a GPU (NVIDIA RTX 6000 Ada) took several hundred seconds, which is 5 orders of magnitude longer than the time required for the FEM solution. This contrast diminishes as the model's dimensionality, size, and complexity increase. To illustrate this, simulating the 3D forward model of the Dabie Shan in Sect. 5.1 on a $\mathrm{120}\times \mathrm{120}\times \mathrm{31}$ grid over 150 time steps took approximately 2400 s using Pecube, whereas training the PINN with 320 neurons over 400 000 iterations required around 32 000 s, merely 1 order of magnitude longer than the FEM solution.

For inverse problems, the computational time escalates as the parameter space expands, which necessitates the evaluation of a larger number of numerical models. In Sect. 4.2, the MCMC sampling which required the evaluation of 1 000 000 FEM models consumed approximately 8000 s. This duration represents an increase by 6 orders of magnitude from a single forward solution. In contrast, the inverse analysis using DIRECT and the PINN model, which evaluated over 15 000 age misfit functions, was completed within 7000 s, marking an increase in time by only 1 order of magnitude from the PINN forward model solution. For the 3D Dabie Shan model (Sect. 5.2), the NA inversion, using 30 CPU cores and coupled with a Pecube model on an $\mathrm{80}\times \mathrm{80}\times \mathrm{31}$ grid, required approximately 300 000 s of clock time to explore the search space and evaluate 9060 misfit functions. The cumulative computation time across all processors amounted to ∼ 9 000 000 s, which is 4 orders of magnitude greater than the time used for solving a single forward model. Conversely, the DIRECT search and the PINN model optimization evaluated 38 130 age misfit functions within 560 000 s, which is increased by 1 order of magnitude from the 3D forward model.

Although the computational performance metrics presented herein are not derived from rigorous benchmarks, they provide a preliminary comparison of the efficiency of PINNs and traditional numerical methods. It is evident that, for inverse problems, PINNs have an advantage due to their seamless integration of parameter search and forward model optimization within one single training process. This integration prevents PINNs from an exponential increase in computational time as the dimensionality of the inverse problem increases, illustrating the method's potential for inverse analysis of expansive models with complex boundary conditions, which can be unattainable with traditional numerical methods.

6.2 Directions for future development

In our applications, we have modelled the rock uplift as a vector in the vertical direction independent of spatial locations. However, thermo-tectonic research often requires us to consider rock uplift and exhumation rates that vary both temporally and spatially. Previous work (e.g., Cai et al., 2021b) has demonstrated the feasibility of incorporating a comprehensive 3D velocity field in the DNN. However, the effectiveness of PINNs in scenarios with abrupt changes in the velocity field requires further exploration or improvement (e.g., Rasht‐Behesht et al., 2022). This capability is vital for studying tectonically active regions, where fault activities can significantly disrupt exhumation and geothermal patterns, and thus will be a focus of our future development of the method.

Our current approach to the inverse problem yields an optimized solution of the heat transfer model. Although the search process results in a collection of sampled parameter values, they are not appraised to quantify the uncertainties associated with the optimized model. To address this, a Bayesian neural network can be employed in the PINN framework. This way, the parameters in the neural network would follow Gaussian distributions in contrast to the deterministic parameters of a traditional fully connected network (Yang et al., 2021). This modification would allow the Bayesian PINNs to accommodate noisy data and provide uncertainty estimates for the inferred parameters in the PDEs.

We have outlined the basic configurations for optimizing PINNs and their integration with the DIRECT search algorithm for inverse analysis. However, we expect that model domains of greater complexity, which require a denser array of collocation points, may need more fine-tuning of the training strategies. Enhancing the algorithm and its configuration for the efficient optimization of PINNs is an active field of research (Jagtap et al., 2020; Shukla et al., 2021). For the inverse analysis, the flexibility in problem formulation with PINNs facilitates the incorporation of various search and optimization algorithms. Our ongoing efforts include integrating and testing these algorithms, to ensure stable training outcomes from PINNs.

7 Conclusions

We have presented applications of physics-informed neural networks (PINNs) in solving crustal heat transfer problems for thermochronology. By harnessing the computational power of deep learning, PINNs offer a flexible and efficient alternative to traditional analytical and numerical techniques. The method allows a straightforward integration of initial and boundary conditions, as well as observed data, and can be used together with gradient-based and derivative-free optimization methods to solve both forward and inverse problems. By applying PINNs to a 1D synthetic model and a 3D natural model of the Dabie Shan under varying tectonic and topographic conditions, we have demonstrated the potential of this approach for estimating the thermal and exhumation histories of the crust during tectonic and landscape evolution. Future work will focus on improving the optimization strategies for inverse models, quantifying the model uncertainties, and expanding the application of PINNs to more intricate geomorphic, geologic, and geodynamic scenarios.