Technical note: Incorporating topographic deflection effects into thermal history modelling

This contribution describes a set of equations and relations to calculate accurate cooling paths through the 2D temperature field of an exhuming region with periodic topography. A 1D model adequately captures the time-varying component of the system regardless of the complexity of the cooling path, making the computation efficient. A series of 2D finite element models demonstrate how temperatures below the periodic mid-slope, or mean topography, can be mapped to those below ridges and valleys, and how these transitions vary with topographic period and amplitude and the ratio of the near-surface geotherm to the atmospheric lapse rate. These new calculations are implemented into HeFTy to support simultaneous modelling of samples collected along topographic profiles, particularly for terranes with long-lived topography that exhumed through an inflected temperature field. Users should bear in mind the limitations of this simplified approach, notably that it presumes a single dominant topographic wavelength, and that it does not encompass other advective processes such as faulting.

1 Introduction

High degrees of exhumation are characteristically accompanied by development of variable topography. During mountain building different areas of the crust rise or fall unevenly due to fault movements and block rotations. Erosion counteracts these incompletely, as erosion rates are a function of relief, and topography approaches a steady state condition while deformation continues (Adams, 1980; Willett and Brandon, 2002). Undulating topography in turn deflects isotherms below, increasing the space between geotherms beneath local highs and compressing them beneath local lows (Fig. 1).

Figure 1 Isotherms for periodic topography with wavelength 24 km, amplitude 1500 m, and steady-state geotherm of 20 °C km⁻¹, and atmospheric lapse rate of 5 °C km⁻¹. (a) Steady-state condition. (b) After 20 million years of exhumation at 0.6 mm yr⁻¹. Calculated using FETKin (Almendral et al., 2015).

It is thus natural that the effects of topography on cooling rates, and how they can be used to infer exhumation rates, has long been of interest to thermochronologists. Early treatments of the effect of sinusoidal topography on steady-state thermal structure were provided by Birch (1950), Carslaw and Jaeger (1959), and Turcotte and Schubert (1982). Stüwe et al. (1994) expanded this work to provide equations for steady-state isotherms given a constant erosion rate. Mancktelow and Grasemann (1997) further broadened the set of available analytical solutions to approximate the time-dependent development of isotherms under constant sinusoidal topography after the onset of erosion at a constant rate. To study more general cases Mancktelow and Grasemann (1997) utilized a 2D finite difference model, and most subsequent work has gone into the development of 2D and 3D numerical solutions (e.g., Almendral et al., 2015; Braun, 2003), although further development of analytical solutions has continued as well (Fox et al., 2014).

Thermal history modelling consists of posing a large number of variable time-temperature (t-T) paths, calculating the response of thermochronometric systems, and comparing the results to measured data to investigate which histories the data are compatible with (Gallagher, 2012; Ketcham, 2005, 2024). The histories posed are typically complex, with cooling rates that vary rather than remaining constant.

A recent update to the HeFTy thermal history modelling software provided for the capability to perform inversions in time-depth space (t-z), rather than time-temperature, by utilizing a 1D finite difference model to convert t-z paths to t-T, under the assumption that changes in depth correspond to erosion or deposition at the surface (Ketcham, 2024). It was specifically intended to help with simultaneous modelling of samples along topographic profiles. It allows for topographic development, and relative vertical rotation (tilting) of sample positions with respect to each other during exhumation (see Fig. 7 in Ketcham, 2024). The 1D model is computationally fast enough that models can be calculated interactively using the graphic interface.

A shortcoming of this implementation is that it neglected the effects of topography on thermal structure, thus treating exhumation through ridges and valleys as equivalent. These effects could be roughly approximated by including topographic development, starting a model as a peneplain and then increasing topographic amplitude during exhumation (e.g., Hestnes et al., 2024; Mackaman-Lofland et al., 2024). However, this shortcut is non-ideal for studying areas with long term or steady state topography.

This contribution presents a method to approximate the thermal deflection effects of periodic topography for complex exhumation and topographic development histories, using only the 1D model calculation. It builds upon, and is informed by, previous work described by Turcotte and Schubert (1982) and Mancktelow and Grasemann (1997). It is partly empirical, utilizing a series of 2D finite element models to map out how variations in various parameters (wavelength, amplitude, geotherm, exhumation rate, topographic development) affect thermal structure to develop a set of relations encompassing the range of conditions likely to be encountered in HeFTy-type thermal history modelling.

2 Methods

Thermal calculations for this work are performed using the FETKin 2D finite element model (Almendral et al., 2015). Grid spacing was 1 km, except for models with topographic wavelengths of 8 km, for which accuracy was improved by reducing to 0.5 km. All models used a constant basal flux boundary condition (50 mW m⁻², resulting in a steady-state gradient of 20 °C km⁻¹), began at steady state, and then run for 12 km of exhumation, all with thermal conductivity K=2.5 W m⁻¹ K⁻¹, thermal diffusivity $\mathit{\kappa }=\mathrm{9.26}\times {\mathrm{10}}^{-\mathrm{7}}$ m² s⁻¹, and atmospheric lapse rate β=0.005 °C m⁻¹.

3 Model Development and Description

3.1 Steady state geotherm with constant periodic topography

For steady-state thermal conditions beneath periodic topography, Turcotte and Schubert (1982, Eq. 4-66, excluding heat production) provide:

$$\begin{array}{}\text{(1)}& T\left(z\right)=T_{\mathrm{0}}+{\displaystyle \frac{q_{\mathrm{m}}z}{K}}+\left({\displaystyle \frac{q_{\mathrm{m}}}{K}}-\mathit{\beta }\right)H_{\mathrm{0}}\cos \left({\displaystyle \frac{\mathrm{2}\mathit{\pi }x}{\mathit{\lambda }}}\right)e^{-\mathrm{2}\mathit{\pi }z/\mathit{\lambda }}\end{array}$$

where z is depth with mean surface = 0 and positive upwards¹, T₀ is surface temperature at mean elevation, q_m is mantle heat flow (W m⁻²), H₀ is the amplitude of topography above and below the mean (m), λ is topographic wavelength (m), and x is position along the surface (m). Their solution is based on assuming the separability of the different components of the heat flow equation: the first term is the mean surface temperature, the second represents heat flow from depth (the steady-state geotherm $g=q_{\mathrm{m}}/K$), and the third reflects the transition in thermal gradient from the atmospheric lapse rate at the surface to the steady-state geotherm as a function of wavelength, dropping by factor of $\mathrm{1}/e$ for every $\mathit{\lambda }/\mathrm{2}\mathit{\pi }$ increase in depth.

This equation is built upon a simpler one that posits a flat Earth surface with periodically varying temperature. To incorporate topography, they make the simplifying assumption that it is shallow (Turcotte and Schubert, 1982, p. 153), and use only the first term of the Taylor series expansion to extrapolate from the mean surface. Although adequate for illustrating crust-scale patterns, this assumption does not reproduce the temperature at and immediately below the crests and valleys of periodic topography (important areas for thermal history modeling), because the derivation assumes the geotherm is linear at z=0 while at the same time has the topographic effect exponentially diminish with depth (and amplify with elevation) away from z=0.

The solution is improved by introducing a new variable z_x representing depth below the local surface:

$$\begin{array}{}\text{(2)}& T\left(z_x\right)=T_{\mathrm{0}}+{\displaystyle \frac{q_{\mathrm{m}}z}{K}}+\left({\displaystyle \frac{q_{\mathrm{m}}}{K}}-\mathit{\beta }\right)H_{\mathrm{0}}\cos \left({\displaystyle \frac{\mathrm{2}\mathit{\pi }x}{\mathit{\lambda }}}\right)e^{-\mathrm{2}\mathit{\pi }{z}_x/\mathit{\lambda }}\end{array}$$

where

$$\begin{array}{}\text{(3)}& z_x=z+H_{\mathrm{0}}\cos \left({\displaystyle \frac{\mathrm{2}\mathit{\pi }x}{\mathit{\lambda }}}\right)\end{array}$$

Or in other words, for the purpose of topographic deflection, the Taylor series originates at the local surface rather than at mean topography.

Another consequence of utilizing separability is that it assumes that the geothermal gradient at the topographic midpoint (z=0) is the same as the geothermal gradient without any topography at all. The finite element modelling in this study shows that this is not quite the case, as, in effect, valleys cool the crust more than hills help retain heat. This effect is usually small, but grows with the severity of topography (increasing H₀ and/or decreasing λ).

Figure 2 (a) Temperature difference in the steady state geotherm between finite element model and Eq. (2) solution below the crest, mid-slope, and valley of topography with wavelength 24 km, amplitude 1500 m and heat flow corresponding to a steady-state geotherm of 20 °C km⁻¹. (b) Mid-slope anomalies for models with varying wavelength, amplitude, and geotherm. (c) Residuals for maximum temperature offset estimate. (d) Residuals for depth scaling estimate.

These two effects can be seen by comparing Eq. (2) to a finite element model. Figure 2a shows the difference in the steady state geotherm between model and analytical solution below the crest, mid-slope, and valley of topography with wavelength 24 km and amplitude 1500 m. The z_x shift in variable allows the temperatures to closely replicate in the near-surface, but the periodic topography cools the lower crust by ∼ 4 °C in the FE model compared to the analytical solution. All 1D solutions for areas of varying topography share this bias, but a simple correction can remove most of it.

Figure 2b shows the offset of the mid-slope geotherm for a number of additional scenarios varying in topographic wavelength and amplitude, and geothermal gradient. All are fit well by the relation $T_{\mathrm{off}}\left(z\right)=T_{\mathrm{off},\mathrm{\infty }}(\mathrm{1}-\exp (-z/D\left)\right)$, where T_off,∞ is the temperature offset at infinite depth, and D scales its onset with depth.

Inspection indicates that T_off,∞ can be estimated as a function of amplitude, wavelength, and difference between geotherm and lapse rate as T_off,∞(est) = −2.65 $H_{\mathrm{0}}^{\mathrm{2}}$ ($g-\mathit{\beta })/\mathit{\lambda }$ (Fig. 2c), and D is a function of wavelength and amplitude, and can be estimated empirically as D(est) = 0.073 λ+0.22 H₀ (Fig. 2d). The former correction is imprecise at larger offsets, but even in the worst, severe case (8 km wavelength, 1.5 km amplitude) the correction is only off by 1.1 °C.

The depth correction for 1D solutions of temperature at the topographic mid-slope (subscript “m”) is thus:

$$\begin{array}{}\text{(4)}& T_{\mathrm{m}}\left(z\right)=T_{\mathrm{1}\mathrm{D}}\left(z\right)-{\displaystyle \frac{\mathrm{2.65}{H}_{\mathrm{0}}^{\mathrm{2}}\left(g-\mathit{\beta }\right)}{\mathit{\lambda }}}\left[\mathrm{1}-e^{-z/\left(\mathrm{0.073}\mathit{\lambda }+\mathrm{0.22}{H}_o\right)}\right]\end{array}$$

The resulting correction is 0–10 °C across conditions appropriate for modelling in HeFTy. Although fairly small in absolute temperature terms, this correction also improves the calculation of the geothermal gradient at the mid-slope position, which turns out to be useful in the next section.

3.2 Steady state geotherm with erosion and constant topography

In general, achieving a steady state temperature regime from a standing start (i.e., going from no erosion to a constant rate) takes tens of millions of years, at least for erosion rates capable of creating and sustaining topography. It can be taken as a given that a complete steady state is seldom achieved. But, once conditions are set up, it is useful as a marker for where isotherms are converging toward.

Stüwe et al. (1994) and Mancktelow and Grasemann (1997) derived analytical solutions for steady-state thermal structure beneath constant periodic topography and constant erosion. The Mancktelow and Grasemann (1997) formulation for erosion occurring at rate u, omitting radiogenic heating and keeping the z_x shift, contains the same three elements as Eq. (2), in altered form:

$$\begin{array}{}\text{(5)}& T=T_{\mathrm{0}}+\mathit{\gamma }\left[\mathrm{1}-e^{-uz/\mathit{\kappa }}\right]+\left({\displaystyle \frac{u}{\mathit{\kappa }}}\mathit{\gamma }-\mathit{\beta }\right)H_{\mathrm{0}}\cos \left({\displaystyle \frac{\mathrm{2}\mathit{\pi }x}{\mathit{\lambda }}}\right)e^{m_{\mathrm{2}}{z}_x}\end{array}$$

where κ is thermal diffusivity, and

$$\begin{array}{}\text{(6)}& \mathit{\gamma }={\displaystyle \frac{T_{\mathrm{S}}-T_L}{\mathrm{1}-e^{-Lu/\mathit{\kappa }}}}\end{array}$$

enforces a constant basal temperature boundary condition where T_S and T_L are the temperatures at the mean surface and at depth L, respectively. Either a constant temperature or constant flux basal boundary can be used for this problem, with the difference being that a constant flux condition will allow the basal temperature to rise toward a new limit (Braun, 2025). Also:

$$\begin{array}{}\text{(7)}& m_{\mathrm{2}}=-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left[{\displaystyle \frac{u}{\mathit{\kappa }}}+\sqrt{{\left({\displaystyle \frac{u}{\mathit{\kappa }}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{\mathrm{4}\mathit{\pi }}{\mathit{\lambda }}}\right)}^{\mathrm{2}}}\right]\end{array}$$

The m₂ root term shows that the transition from surface to geothermal gradient occurs over a shallower interval than the no-erosion case, as for any positive u it has a higher magnitude than the corresponding exponential term in Eq. (2). At the same time, the magnitude of the isotherm displacement rises, as the steady state gradient given by uγ/κ is higher than the no-erosion condition.

Figure 3 Offsets of temperature with depth below local surface at periodic topography ridges (negative values) and valleys (positive values) compared to mid-slope (mean) elevation, at 1 Myr intervals for model with parameters as in Fig. 1b. (a) Oblique view. (b) Temperature offset versus depth. (c) Temperature offset versus time.

3.3 Geotherm evolution at crests and valleys with erosion and constant or growing topography

The different offsets from the analytical solution for the geotherm between the ridge and valley of topography in Fig. 2a indicate that different corrections are necessary for each location. Figure 3 shows different views of the temperature difference at the same z_x between the mid-slope versus the ridge (left; negative values) and valley (right; positive), at 1 Myr intervals for a 2D finite element model with 24 km wavelength, 1500 m deflection, and 0.6 mm yr⁻¹ erosion. The constant displacement between the curves at 0 km reflects the unchanging topography and lapse rate. The bulging of the curves with increasing time since onset of erosion reflects the bunching of isotherms toward the surface due to exhumation, which diminishes with depth but grows slowly through time as the bulging propagates downward.

The diagrams illustrate that the change in temperature with depth and time due to topographic deflection has three components: maximum magnitude at great depth, magnitude variation with depth, and magnitude evolution with time. These are separable. The depth of the transition from the surface to maximum magnitude offset does not change through time and is dominantly controlled by the topographic wavelength (Fig. 1b). The offset magnitude at depth increases through time, and shows no sign of converging over these time scales.

The progressive bunching of isotherms also occurs at the mid-slope position, suggesting that the mid-slope geotherm can be used as a normalizing factor. Dividing the offsets by the current local geotherm at the same z-t position beneath the mid-slope location ($\mathrm{d}T/\mathrm{d}{z}_{\mathrm{m}}$ (z)) leads to Fig. 4.

Figure 4 Normalized offsets with depth below local surface at periodic topography ridges (negative values) and valleys (positive values) compared to mid-slope (mean) elevation, at 1 Myr intervals for model with parameters as in Fig. 1b. (a) Oblique view. (b) Normalized offset versus depth. (c) Normalized offset versus time.

It is evident that the time evolution of the offsets is almost entirely characterized by the time evolution of the geotherm at the mid-slope (mean elevation) point. There is a slight evolution in the nearest-surface points, but this is due to the evolution of the geotherm relative to the constant lapse rate. Figure 4 also clarifies the asymmetry between the ridge and valley offsets, indicating that they require separate treatment.

Figure 5 As Fig. 4, for model with developing topography. (a) Oblique view. (b) Normalized offset versus depth. (c) Normalized offset versus time.

The same test can be run with growing topography, starting from no topography and evolving to amplitude 1500 m after 20 Myr (Fig. 5). Normalizing using the mid-slope geotherm again results in a regular pattern for the anomalies at the ridge and valley positions, with depth relationships closely resembling, and converging toward, their counterparts in Fig. 4. This result suggests that, regardless of the history of topographic development, the current state of topography and the current mid-slope geotherm adequately describe the time-varying component of the system.

The asymmetric offsets in Figs. 4b and 5b can be fitted well with an equation of the form:

$$\begin{array}{}\text{(8)}& T_{\mathrm{asym}}\left(z\right)=\left(c_{\mathrm{0}}+c_{\mathrm{1}}{e}^{-\frac{z}{c_{\mathrm{2}}}}\right){\left[{\displaystyle \frac{\mathrm{d}T}{\mathrm{d}z}}\right]}_{\mathrm{m}}\left(z\right)\end{array}$$

in which c₀ defines the maximum normalized offset at depth, c₁ the normalized offset at the surface, and c₂ is the depth scale of the transition. It is readily apparent that c₀ corresponds to H₀, positive for valleys and negative for ridges, as it reflects the difference in temperature at depths where isotherms no longer deflect. Similarly, c₁ corresponds to the surface temperature offset, which is determined by H₀ and β. The set of relations for ridges (subscript r) and valleys (subscript v) are thus:

$$\begin{array}{}\text{(9)}& \begin{array}{rl}& c_{\mathrm{0},v}=H_{\mathrm{0}}\\ & c_{\mathrm{0},r}=-H_{\mathrm{0}}\\ & c_{\mathrm{1},v}=-H_{\mathrm{0}}\left(\mathrm{1}-\mathit{\beta }/{\left[{\displaystyle \frac{\mathrm{d}T}{\mathrm{d}z}}\right]}_{\mathrm{m}}\left(\mathrm{0}\right)\right)\\ & c_{\mathrm{1},v}=H_{\mathrm{0}}\left(\mathrm{1}-\mathit{\beta }/{\left[{\displaystyle \frac{\mathrm{d}T}{\mathrm{d}z}}\right]}_{\mathrm{m}}\left(\mathrm{0}\right)\right)\end{array}\end{array}$$

with H₀ representing current amplitude at a given time (i.e. during topographic development) in units of km, and geotherm and lapse rate in °C km⁻¹. Equations (8) and (9) combined provide:

$$\begin{array}{}\text{(10)}& T_{\mathrm{asym}}\left(z\right)=\pm H_{\mathrm{0}}\left\{\mathrm{1}-\left(\mathrm{1}-\mathit{\beta }/{\left[{\displaystyle \frac{\mathrm{d}T}{\mathrm{d}z}}\right]}_{\mathrm{m}}\left(\mathrm{0}\right)\right)e^{-\frac{z}{c_{\mathrm{2}}}}\right\}{\left[{\displaystyle \frac{\mathrm{d}T}{\mathrm{d}z}}\right]}_{\mathrm{m}}\left(z\right)\end{array}$$

where the sign is positive for valleys and negative for ridges, and the c₂ variable changes between these cases.

To explore the depth scaling (c₂), a series of 20 FETKin models varying in wavelength, amplitude, background steady-state geotherm, erosion rate, and developing versus constant topography was used. Model settings are listed in Table 1, as well as the corresponding c₂ values for ridges and valleys fitted to the model outcomes.

Table 1 FETKin models used to derive empirical relations between physical parameters and isotherm deflection.

^a c: constant topography, g: growing topography. ^b Sensitivity test; model not used for deriving empirical relations for c₂

Figure 6 Residual plots for estimates of depth scaling parameter for ridges and valleys.

The depth scaling is dominantly a function of the topographic wavelength, and secondarily of topographic amplitude. A good empirical fit is given by:

$$\begin{array}{}\text{(11)}& \begin{array}{rl}& c_{\mathrm{2},\mathrm{v}}=\mathrm{0.166}\mathit{\lambda }-\mathrm{0.473}{H}_{\mathrm{0}}\\ & c_{\mathrm{2},\mathrm{r}}=\mathrm{0.156}\mathit{\lambda }+\mathrm{0.805}{H}_{\mathrm{0}}\end{array}\end{array}$$

As shown in Fig. 6, these formulas do an excellent job of reproducing c₂ for a wide range of conditions. The principal apparently systematic departure is for growing versus static topography, with the former increasing the depth scaling by ∼0.5 km compared to similar constant-topography models.

Table 1 also lists the maximum and median absolute temperature differences over the top 10 km of the crust between finite element models and estimates using the FE-calculated mid-slope geotherm and Eqs. (2), (10), and (11). All errors at all times are less than 3 °C, with an average maximum error among models of 1 °C, and the average median error is 0.6 °C. The growing-topography scenarios also show small errors, suggesting that adjusting for such scenarios is a second-order consideration.

Although there is always a danger that empirical relations can degrade if applied outside of the model parameter space used to develop them, the scenarios used for generating Eq. (11) (wavelengths of 8–36 km, total relief of 1–3 km, exhumation rates up to 0.6 mm yr⁻¹ or km Myr⁻¹) reasonably bracket the use cases appropriate for modelling in HeFTy. The conditions modelled also effectively account for changes in conductivity and diffusivity, as the former offsets the geotherm (Eq. 1) and the latter offsets the exhumation rate (Eqs. 5, 6, 7).

To test the effects of venturing beyond these limits, an extreme case was run that generates 3 km of relief with exhumation rates of 6 mm yr⁻¹ over 2 Myr (Table 1, model 21). While the maximum error exceeds 14 °C in valleys, the high dynamic geothermal gradients (>80 °C km⁻¹) means that this corresponds to a vertical offset of less than 200 m or a timing error of ∼30 kyr, which may be within tolerances. Still, in such severe cases investigators should consider using alternative tools that provide full 2D or 3D solutions.

3.4 Varying exhumation rates and interpolating across topography

The topographic correction model also achieves low errors when exhumation rates vary. For example, a model akin to that shown in Fig. 3 in which erosion proceeds for 10 Myr at 0.6 mm yr⁻¹ and then slows to 0.2 mm yr⁻¹ over the subsequent 30 Myr yields a maximum error of 1.2 °C and a median error of 0.8 °C, only marginally above the similar constant-rate model in Table 1 (model 17). This indicates that the mid-slope 1D model adequately captures the evolution of thermal structure even for complex t-z histories.

Interpolating deflection between ridges and valleys is done via the cosine term in Eq. (1). For ideal sinusoidal topography elevation E varies around mean elevation E_m as:

$$\begin{array}{}\text{(12)}& E\left(x\right)=E_{\mathrm{m}}+H_{\mathrm{0}}\cos \left({\displaystyle \frac{\mathrm{2}\mathit{\pi }x}{\mathit{\lambda }}}\right)\end{array}$$

allowing the cosine term to be calculated as:

$$\begin{array}{}\text{(13)}& {\displaystyle \frac{E\left(x\right)-E_{\mathrm{m}}}{H_{\mathrm{0}}}}=\cos \left({\displaystyle \frac{\mathrm{2}\mathit{\pi }x}{\mathit{\lambda }}}\right)\end{array}$$

The complete equation for temperature asymmetry then becomes:

$$\begin{array}{}\text{(14)}& \begin{array}{rl}{T}_{\mathrm{asym}}\left(z,x\right)& =\left(E_{\mathrm{m}}-E\left(x\right)\right)\\ & \cdot \left\{\mathrm{1}-\left(\mathrm{1}-\mathit{\beta }/{\left[{\displaystyle \frac{\mathrm{d}T}{\mathrm{d}z}}\right]}_{\mathrm{m}}\left(\mathrm{0}\right)\right)e^{-\frac{z}{c_{\mathrm{2}}}}\right\}{\left[{\displaystyle \frac{\mathrm{d}T}{\mathrm{d}z}}\right]}_{\mathrm{m}}\left(z\right)\end{array}\end{array}$$

Figure 7 Effects of topographic correction. Each model shows the set of t-T histories assuming the control sample (sample a in all cases) is exhumed from 7 km depth to the surface over a 40 Myr span. For all HeFTy models, MSL temperature is 20 °C, basal heat flux is 60 mW m⁻², thermal conductivity 3 W m⁻¹ K⁻¹, heat capacity 800 J kg⁻¹ K⁻¹, density 2700 kg m⁻³, maximum depth 120 km, grid spacing 2 km, and β 6.5 °C km⁻¹. For the periodic correction, λ=54 km, H₀=915 m, and E_m=1415 m. In all cases, the black curve is for sample position a in (i), and the teal curves are for the others, with sample e always the most distant. Models shown in (a), (c), (e), and (g) do not use the periodic topography correction, which can be compared with models (b), (d), (f), and (h) which do.

4 Implementation in HeFTy

The foregoing calculations have been implemented into HeFTy (version 2.2.0 and higher; version 2.3.1 provided in repository, along with finite element results (Ketcham, 2025)) to provide a better basis for thermal history modelling across elevation transects. In the new version, the position of the 1D thermal model has been changed from the “control” sample location, against which the positions of the other samples are determined in models involving topographic development or tilting (Ketcham, 2024), to the mid-slope. This reduces the impact of choice of control sample. Whereas part of the need for carefully choosing the control sample stemmed from avoiding generating exhumation paths that result in some samples being sent above the ground surface prior to the end of the model, the program now flags these cases during forward modelling, and avoids them automatically during inverse modelling.

4.1 Guidelines for data entry

As with other multi-sample modelling in HeFTy, it is best to fit individual models for each sample prior to combining them, as an initial check for their respective viability and mutual consistency. Regional topographic parameters of wavelength, amplitude, and mean elevation are specified independently, and all samples should lie within the topographic range specified (mean elevation ± amplitude). The wavelength and amplitude should correspond to the full extent of the region being modelled, even if the sample set does not include the peak and valley positions. Although position along the transect is entered via geographic coordinates, it is only used for display purposes; HeFTy uses sample elevation to define its location along the topographic period for the thermal correction according to Eqs. (13) and (14). Users should bear in mind how this data entry scheme may be affected by shorter-wavelength complexities, as, for example, a small-scale reversal in topography will effectively switch sample ordering.

This correction and its implementation are designed with the principal anticipated use case being a single traverse spanning a single major elevation change. However, topography can have multiple periodic components, and users should be mindful that dates from lower-closure-temperature thermochronometers will be sensitive to shorter topographic wavelengths (Stüwe et al., 1994). One can estimate a critical wavelength below which topographic deflection effects will be irrelevant to a thermochronometer as the ratio of the closure temperature (T_c) to the geothermal gradient (Braun, 2002), and this may be used as guidance for additional exploration of lower-T_c systems by modelling subsets of samples across smaller-scale variations. A Fourier deconvolution (e.g., Glotzbach et al., 2015) may also provide helpful insight into the expected effects of 3D topography and superposed wave components in a given study area.

4.2 Demonstration of topographic correction

Figure 7 shows a series of models results following from those used by Ketcham (2024, Fig. 8) without and then with the periodic topographic correction. In Fig. 7a, simple exhumation under constant (static) topography results in a series of parallel paths when there is no correction, as each path differs only in its surface boundary condition, defined by the lapse rate. Applying the topographic correction in Fig. 7b allows the specimens to transition from depth, where they are separated by ∼25 °C, to the surface where their separation is 9 °C. Note that at the starting depth of 7 km the steady-state isotherms are not fully flat (see Fig. 1a), and thus the temperature offset among samples in the initial condition (steady state) does not quite correspond to the geothermal gradient. If topography begins flat at 40 Ma and develops during exhumation, the difference between the non-corrected and corrected path sets (Fig. 7c and d, respectively) is very subtle, with a maximum temperature difference of only about 2 °C.

Results are similar when tilting is included, in this case with a 10° clockwise rotation from past to present such that the topographic high corresponds with the uplifted side of the profile. Under constant topography the periodic correction again makes a difference of up to 10 °C, only this time the dispersion of paths decreases at depth because the topographic effect is inverted (Fig. 7e, f). If topography starts flat, then differences are more subtle, basically because there is no correction to make at the model start, and models all converge to the same ending point. Combined, results in Fig. 7 demonstrate that topographic development has served as a reasonable proxy for periodic topography effects in previous studies using topographic profile modelling in HeFTy (Hestnes et al., 2024; Mackaman-Lofland et al., 2024), but that the models cannot distinguish whether topographic development was really required to fit the data.

With regard to tilting, it should be noted that the 1D model does not incorporate the additional deflection of isotherms up or down due to block rotation; HeFTy only projects samples higher or lower in the crust. The error stemming from this omission is minimized by having the mid-slope be the rotational pivot rather than the control sample, and is expected to usually be second-order.

5 Discussion

As shown in the example, the absolute temperature changes enabled by the periodic correction developed here will tend to be modest, in many cases less than ±10 °C (or net 20 °C between ridge and valley), although the absolute value will depend on the severity of topography. The primary utility in this correction, aside from increasing the overall accuracy of the solution, is to better allow topographic effects to be “felt” at depth during exhumation, in turn increasing the ability for the data to discern when topographic evolution began. In the context of multi-sample thermal history modelling, even relatively mild relative excursions among samples that are geographically linked can significantly affect their ability to have their thermochronometric data be matched in tandem.

It should also be stressed that this functionality is meant to test relatively simple scenarios and provide a first-order approximation of the geological processes being considered; it should not be seen as depicting or reproducing reality. In particular, in studies where the exhumation history is more complex, such as involving spatially variable and/or significant 3D topographic development, or is directly impacted by additional advective processes such as folding, faulting, or fluid movement, more sophisticated modelling tools will be preferable for posing questions with the appropriate level of fidelity and detail (e.g., Almendral et al., 2015; Braun, 2003; Bernard et al., 2022).