Geothermal implications of the lithosphere's thermal structure in northern Pakistan

Conventional geothermal resources are typically associated with volcanically active plate boundaries, yet collisional orogens can also sustain elevated heat flow through radiogenic enrichment, crustal thickening, and rapid exhumation. Northern Pakistan, encompassing the Himalaya, Kohistan, and Karakoram terranes, hosts numerous hot springs aligned with major fault zones despite the absence of active volcanism. The origin of this anomalous heat remains debated, reflecting the lack of surface heat-flow measurements and limited geophysical constraints on the lithosphere. To address this gap, we apply 1D steady-state, 1D transient, and 2D advective-conductive thermal models with radiogenic heat production and exhumation-driven heat transport to the Nanga Parbat Massif, Kohistan arc, and Karakoram terrane. These models, in the absence of direct heat-flow measurements in the region, remain uncalibrated; however, they bracket a range of plausible thermal states. Steady-state models show a strong dependence of geotherms on crustal radiogenic heat production (RHP): in the NPM, upper-crustal enrichment (4–5 µW m⁻³) yields surface heat flow of 85–120 mW m⁻², whereas Kohistan produces lower values (50–85 mW m⁻²) due to its mafic-dominated crust. Karakoram yields intermediate heat flow (65–103 mW m⁻²), with RHP concentrated in the batholith and metamorphic complexes. These models also reveal that surface-derived RHP values cannot be extrapolated to the mid- and lower crust: doing so yields Moho temperatures inconsistent with petrological and geophysical constraints, implying redistribution of heat-producing elements during partial melting and differentiation. 1D transient models suggest that exhumation of 2–3 mm yr⁻¹ in the NPM can further amplify geotherms, producing surface heat flow up to 220–250 mW m⁻² and inverting deep geotherms at 20 km in exploratory end-member scenarios. 2D models capture the combined effects of radiogenic enrichment, exhumation, and rugged topography. Isotherms are compressed beneath valleys and expanded beneath peaks, with the strongest thermal anomalies localized in the NPM and Karakoram. Surface heat-flow patterns reflect these contrasts, ranging from ∼120 mW m⁻² (moderate scenarios) to nearly 180 mW m⁻² (fast exhumation). Péclet numbers ranging from below unity to greater than unity indicate that exhumation progressively increases the contribution of advective heat transport. Crustal differentiation indices calculated from the model-derived heat-flow estimates are consistent with strong upper-crustal enrichment of heat-producing elements in the NPM and Karakoram, indicating their redistribution during crustal thickening and partial melting. The models suggest that the region can sustain anomalously high heat flow through the interplay of RHP, exhumation, and crustal differentiation. For northern Pakistan, these results provide a physical basis for the widespread hydrothermal activity and indicate elevated subsurface temperatures in the NPM and Karakoram, identifying them as regions that warrant further site-specific investigation of reservoir properties and hydrogeological conditions.

1 Introduction

Shallow crustal regions of elevated heat flow are typically associated with active volcanism along plate boundaries, representing conventional geothermal resources (Jolie et al., 2021). Recently, ancient cratonic basements have gained attention as potential unconventional geothermal targets (Moeck, 2014). Crystalline rocks such as gneisses and granites often contain high concentrations of radiogenic elements, enhancing heat production and surface heat flow, which makes them promising candidates for hot-dry rock geothermal systems (Hillis et al., 2004). Temperatures in these systems can reach ∼240 °C at 3.5 km depth in high heat-producing granitic basements (3.8–8.7 µW m⁻³), as observed in the Cooper Basin, Australia (Hillis et al., 2004; Beardsmore, 2005; Meixner et al., 2014). Crystalline basements and granitic intrusions, commonly forming the cores of orogenic belts, are large, accessible, and associated with regional orogenic heat flow, making them viable for geothermal exploration (McCay and Younger, 2017; Gnojek et al., 2018; Wanner et al., 2019).

Northern Pakistan lies within the Indian-Asian collision zone, encompassing the Himalaya, Karakoram, and Tibet. More than 500 hot springs occur here, forming the Himalayan Geothermal Belt (Tong and Zhang, 1981). This ∼3000 km long E-W trending belt extends from the western Himalayas to northwest India and includes parts of the Pamir and Tibet (Fig. 1a). Although the region hosts numerous Cenozoic crustal intrusions, the absence of active arc or rift volcanism and the low ${}^{\mathrm{3}}\mathrm{He}{/}^{\mathrm{4}}\mathrm{He}$ ratios in hot spring gases are consistent with a predominantly crustal contribution to the observed thermal anomaly rather than a direct mantle contribution (Hochstein and Zhongke, 1995). Geothermal activity is concentrated along regional faults and suture zones, with hot springs discharging in deeply incised valleys (Hochstein and Regenauer-Lieb, 1998). Interaction of meteoric water with rapidly exhuming hot rocks during ongoing collision contributes to the hydrothermal activity (Chamberlain et al., 1995).

Figure 1 (a) Distribution of hot springs of the Himalaya Geothermal Belt (after Hochstein and Zhongke, 1995) within the hinterland of the India-Asia collision zone. The black rectangle shows the location of (b). (b) Distribution of geothermal manifestations (hot springs) in the hinterland of the western Himalaya and Karakoram and associated heat-flow map (modified after Shanker, 1988; Thussu, 2002; Craig et al., 2013). Heat flow is calculated from geothermal fields in Indian Kashmir and has been extrapolated towards northwest in Pakistan (hatched pattern) by Craig et al. (2013). Abbreviations: KF – Karakoram Fault; ISZ – Indus Suture Zone; MBT – Main Boundary Thrust; MCT – Main Central Thrust; MFT – Main Frontal Thrust; MKT – Main Karakoram Thrust; MMT – Main Mantle Thrust; STD – South Tibet Detachment.

Comprehensive surface heat-flow data are scarce in the Himalaya-Karakoram, limiting knowledge of subsurface geothermal gradients and the lithosphere's thermal state. Available subsurface data, such as boreholes, seismic, gravity, and magnetic surveys are limited, hindering imaging of lithospheric and crustal architecture. While high heat flow (>100 mW m⁻²) has been reported from geothermal sites in Indian Kashmir (Fig. 1b; Shanker, 1988; Craig et al., 2013), these values may not represent background heat flow due to advective transport in hydrothermal systems. Accurate lithospheric structure and rock properties are essential for thermal modelling, and without them, models only yield broad estimates (Goes et al., 2020).

This study estimates the lithospheric thermal state of the Himalaya, Kohistan, and Karakoram using 1D and 2D thermal models. The models incorporate surface radiogenic heat production, literature-based geophysical parameters, and thermal properties. We assess the influence of variable crustal radiogenic heat production, crustal thickness, and exhumation on geotherms and surface heat flow. Using the 1D models and petrophysical parameters from Anees et al. (2024), 2D models along a cross-section illustrate upper-crustal isotherms and their geothermal implications.

2 Tectonic and lithospheric framework

The western Himalayan-Karakoram orogen, forming most of northern Pakistan, is tectonically subdivided from north to south into the Karakoram block (Asian Plate), Kohistan (-Ladakh) arc, and the Indian Plate (Nanga Parbat Massif) (Gaetani et al., 1996) (Fig. 2a). The Kohistan arc, largely intra-oceanic in origin, was sandwiched between the Asian and Indian plates during the closure of the Neotethys along two suture zones, followed by their collision in the Cenozoic era (Tahirkheli, 1979). These tectonic boundaries are the Main Karakoram Thrust (MKT) or Shyok suture zone, separating the Karakoram and Kohistan, and the Main Mantle Thrust (MMT) or Indus suture zone, marking the suture between the Kohistan arc and India (Gansser, 1980; Kazmi et al., 1984).

Figure 2 (a) Regional tectonic map of the western Himalaya and surrounding regions (modified after Beck et al., 1996; Badshah et al., 2000; Dipietro and Pogue, 2004; Faisal et al., 2018). The Kohistan and Ladakh arcs (in green) are sandwiched between the colliding Indian (blue) and Asian (orange) plates. The study area comprises the Karakoram, Kohistan-Ladakh arc, and Nanga Parbat Massif. (b) Regional lithospheric structure of Himalaya, Karakoram, and Pamir orogenic belts (after Mechie et al., 2012).

The study area encompasses three tectonic domains: the Nanga Parbat Massif (NPM) of the Himalaya, the Kohistan (-Ladakh) batholith, and the Karakoram batholith. The NPM, the northernmost part of the Himalayas in Pakistan, exposes Proterozoic Indian basement up to ∼8 km elevation. It has experienced exceptionally high exhumation of up to 30 km during the Himalayan orogeny (15 km in the last 3 Ma), forming an N-S striking syntaxis thrusting onto the Kohistan arc (Butler et al., 1989; Zeitler et al., 1993).

The Kohistan-Ladakh arc extends ∼700 km in an E-W direction and is positioned between the Karakoram (in the north, Asian Plate) along the MKT and the Indian Plate (in the south) along the MMT. This Cretaceous intra-oceanic island-arc complex has been deformed and partly metamorphosed due to obduction onto the Indian Plate. Exhumation of the NPM eroded overlying Kohistan crust, splitting it into the Kohistan arc in the west and the Ladakh arc in the east (Tahirkheli, 1979; Petterson, 2010). The Kohistan arc preserves a complete crustal section, from mantle to uppermost crustal rocks, whereas the Ladakh arc predominantly exposes upper-crustal volcanic and plutonic rocks (Treloar et al., 1990).

The Karakoram terrane represents the southern margin of the Asian Plate and is subdivided into the Northern Karakoram sedimentary belt, Karakoram Batholith, and Karakoram Metamorphic Complex (Searle and Khan, 1996). The northern Karakoram forms a back-arc sedimentary basin with Paleozoic carbonates and clastic deposits extending northward toward the southern Pamir (Palin et al., 2012). The Karakoram Batholith, an E-W trending, 600 km long, and up to 30 km wide axial intrusion, separates the northern sedimentary sequence from the southern metamorphic complex (Searle et al., 1989; Crawford and Searle, 1992). Lithologies include mid-Cretaceous subduction-related granodiorites and diorites in the west and central parts (Hunza plutonic unit) and Jurassic diorite gneisses intruded by post-collisional leucogranites and Miocene syenites in the east (Baltoro plutonic unit, Kande plutonic complex) (Villa et al., 1996; Searle et al., 2010).

The present lithospheric structure reflects these complex subduction and suturing processes (Fig. 2b). Crustal thickening during the Himalayan orogeny occurred via underthrusting of the Indian Plate beneath Kohistan and Karakoram. Moho depth increases from 30–40 km beneath the Indian Shield (Acton et al., 2010) to 50–60 km beneath the Himalayan front and 70–80 km within the wedge (Rai et al., 2006; Hazarika et al., 2014; Priestley et al., 2019), with the Indian Plate subducting to >90 km beneath Kohistan and Ladakh (Kufner et al., 2016). Lithosphere-asthenosphere boundary estimates remain variable, between 150–250 km (Li and Mashele, 2009; Kumar et al., 2022; Jadoon et al., 2021), and major suture zones appear decoupled from the lithospheric mantle (Kind and Yuan, 2010; Zhao et al., 2010).

Seismic studies divide the Himalayan crystalline crust into four layers based on P-wave velocities (Guerra et al., 1983; Bhukta et al., 2006). The orogenic core consists of upper to mid-crustal metamorphic rocks corresponding to mid-Archean to Neoproterozoic shields, underlain by a granulitic lower crust (Grujic et al., 2011; Warren et al., 2011). The Kohistan crust is northward-tilted, exposing upper-crustal granitoids and volcanic rocks to the north and mid- to lower-crustal amphibolites and granulites to the south. Its present thickness is ∼10–25 km (Malinconico, 1986; Burg et al., 2006; Burg and Bouilhol, 2019).

Subducted Indian lower crust extends as far north as Tibet and the Pamir (Wittlinger et al., 2009; Kumar et al., 2022; Zhao et al., 2010), although its geometry beneath Karakoram and western Tibet remains unresolved. The Karakoram crust comprises a lower crust of unradiogenic felsic and mafic granulites and a radiogenically fertile gneissic middle crust (Searle and Hacker, 2019), with cold, dry granulitic Indian crust likely underlying western Tibet (Priestley et al., 2008). The upper- and mid-Indian crust levels beneath Karakoram are believed to have been largely accreted into the Himalaya via lateral material transfer (Razi et al., 2014).

3 1D steady-state conductive thermal model

This section describes the physical and mathematical background of steady-state heat transport in solid material as well as the scenarios explored in this work. The steady-state models here provide the reference conductive state against which the transient effect of exhumation can be isolated and quantified; thereby providing the baseline for the increment factors in Table 3. Additionally, they permit an efficient sensitivity exploration of structural parameters such as LAB depth, heat-producing layer thickness and layer-wise RHP which would be computationally and conceptually redundant in transient mode. We note that, in the rapidly exhuming Nanga Parbat Massif the steady-state assumption is clearly violated at present, which is precisely why the transient and two-dimensional models are introduced in Sects. 4 and 5.

3.1 Mathematical formulation

A nonlinear, steady-state, one-dimensional conductive heat-transfer equation is solved as follows:

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial }{\partial z}}\left[k(T,z){\displaystyle \frac{\partial T}{\partial z}}\right]+A\left(z\right)=\mathrm{0},\end{array}$$

where T is temperature [K], k(T,z) is a temperature-dependent thermal conductivity [$\mathrm{W}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$], and A(z) is radiogenic heat production [W m⁻³ in the governing equations (model input values are reported in µW m⁻³)], both expressed as functions of depth. Depth z is positive downward, and radiogenic heat production A is defined as a positive volumetric heat source. Accordingly, upward exhumation corresponds to a negative vertical material velocity in the transient models of Sect. 4.

Temperature-dependent thermal conductivity is modelled following Jaupart et al. (2016), with an additional temperature-dependent radiative component at temperatures >1000 K (Jaupart and Mareschal, 2005). The governing equations are:

$$\begin{array}{}\text{(2)}& \begin{array}{rl}{k}_{\text{lat}}(T,z)& =\mathrm{2.26}-{\displaystyle \frac{\mathrm{618.251}}{T}}\\ & +k_{\mathrm{o}}\left(z\right)\left({\displaystyle \frac{\mathrm{355.576}}{T}}-\mathrm{0.30247}\right),\end{array}\end{array}$$

and

$$\begin{array}{}\text{(3)}& \begin{array}{rl}& k_{\text{rad}}\left(T\right)=\mathrm{0.37}\times {\mathrm{10}}^{-\mathrm{9}}{T}^{\mathrm{3}},T\ge \mathrm{1000}\mathrm{K};\\ & k_{\text{rad}}\left(T\right)=\mathrm{0},T<\mathrm{1000}\mathrm{K},\end{array}\end{array}$$

where k_o is the thermal conductivity of the rock under surface conditions, which differs between lithological layers (Table 1), k_lat is the lattice (non-radiative) contribution and k_rad the radiative contribution, which is negligible below 1000 K. The total thermal conductivity is therefore:

$$\begin{array}{}\text{(4)}& k(T,z)=k_{\text{lat}}(T,z)+k_{\text{rad}}\left(T\right).\end{array}$$

Table 1 Parameters of median steady-state conductive thermal models of the Nanga Parbat Massif, Kohistan arc, and Karakoram based on previous studies (Anees et al., 2023, 2024) and published literature (Hasterok and Chapman, 2011; Roy and Mareschal, 2011). Ind-1 and Ind-2, respectively, refer to high and low heat-producing crustal layers. Koh-F and Koh-M, respectively, refer to the upper felsic and lower mafic Kohistan crust overlying the Indian crust. KB-Karakoram batholith, KMC-Karakoram metamorphic complex.

Equation (1) is solved using Dirichlet boundary conditions with fixed temperatures at the surface and the base of the lithosphere. The model allows an arbitrary number of layers, each with distinct k_o and A. To handle abrupt changes in properties, a finite-difference “half-station” scheme is employed (Hindle and Besson, 2023). A tridiagonal matrix algorithm combined with fixed-point iteration solves the nonlinear system:

$$\begin{array}{}\text{(5)}& \left(\mathbf{K}\right[k]T{)}_i={\displaystyle \frac{k_{i+\mathrm{1}/\mathrm{2}}(T_{i+\mathrm{1}}-T_i)-k_{i-\mathrm{1}/\mathrm{2}}(T_i-T_{i-\mathrm{1}})}{\mathrm{\Delta }{z}^{\mathrm{2}}}}.\end{array}$$

The geological model contains abrupt lithological contrasts in conductivity and heat production, so the conservative operator $\partial /\partial z$ ($k\partial T/\partial z$) is discretized directly at half stations. For a uniform grid, this gives Eq. (5), where K[k] denotes the discrete conservative conduction operator and $k_{i+\mathrm{1}/\mathrm{2}}$ and $k_{i-\mathrm{1}/\mathrm{2}}$ are conductivities evaluated at the cell faces.

No numerical derivative of k is required. This matters because first expanding the operator as $k{\partial }^{\mathrm{2}}T/\partial z^{\mathrm{2}}+(\partial k/\partial z)(\partial T/\partial z)$ implicitly assumes that k is differentiable, whereas at an abrupt lithological boundary $\partial k/\partial z$ is not an ordinary function. The half-station formulation instead works with conductive fluxes across cell faces and preserves the conservative structure of the equation. The same operator is used in the transient calculations of Sect. 4.

3.2 Modelling scenarios and results

Several scenarios were tested to explore the sensitivity of crustal temperatures to key parameters such as lithospheric thickness and upper-crustal radiogenic heat production.

3.2.1 Effect of lithosphere–asthenosphere boundary’s (LAB) depth

A basal temperature of 1300 °C is prescribed at the base of the lithosphere. A Neumann (basal heat flow) condition is often preferred in lithospheric modelling and can be anchored to measured surface heat flow in cratonic settings (e.g. Kumar et al., 2007, 2009); however, no such anchor exists for our case, and mantle heat flow beneath the study area would itself have to be assumed. A temperature of approximately 1300 °C is commonly used as a reference temperature for the base of the thermal lithosphere (McKenzie et al., 2005), making a prescribed basal temperature a physically motivated boundary condition.

While the Moho depth is estimated at 60–70 km (Priestley et al., 2019), LAB estimates vary between 150 and 250 km (Li and Mashele, 2009; Jadoon et al., 2021; Kumar et al., 2022). The effect of this 100 km variation on geotherms was tested with three LAB depths: 150, 200, and 250 km (Fig. 3a). A shallower LAB (150 km) slightly increases temperatures at depth, but its effect on upper-crustal temperatures and surface heat flow is negligible. An intermediate LAB depth of 200 km is therefore used in subsequent models.

Figure 3 1D steady-state conductive geotherms showing the effect of thickness of lithospheric and heat-producing layer in the upper crust (a) Geotherms showing variations in temperature (solid lines) and heat flow (dashed lines) with varying lithospheric thickness from 150 to 250 km. (b) Geotherms showing variations in temperature and heat flow with varying thicknesses of heat-producing layer in the upper crust from 5 to 15 km.

3.2.2 Effect of heat-producing layer's (HPL) thickness

Radiogenic heat production generally decreases with depth (Jaupart et al., 2016), but vertical and lateral heterogeneity in the upper crust complicates its prediction within distinct lithological units (Vilà et al., 2010). Additionally, systematic xenolith studies characterizing deep-crustal radiogenic heat production in the exposed intrusions are lacking in the region studied here. Consequently, the vertical partitioning of RHP represents the dominant model uncertainty. Therefore, for upper-crustal radiogenic heat production, scenarios with variable HPL thicknesses of 5, 10, and 15 km were tested using a mean heat production of 4 µW m⁻³ (weighted average of the NPM). The crust is modelled as a 25 km thick upper crust (including HPL) overlying a 45 km thick lower crust, with standard literature-based heat production values of 0.02, 0.4, and 2 µW m⁻³ for mantle, lower crust, and upper crust, respectively (Hasterok and Chapman, 2011). Boundary temperatures are 10 °C at the surface and 1300 °C at the LAB (McKenzie et al., 2005). Temperature-dependent thermal conductivity and radiative transfer are included.

Results (Fig. 3b) show that a 10 km variation in HPL thickness strongly affects upper-crustal geotherms, whereas a 100 km variation in LAB has minimal effect. Calculated temperatures at 10 km depth range from 337 °C (HPL=5 km) to 430 °C (HPL=15 km), increasing to 578–686 °C at 20 km depth. Corresponding surface heat flows are 84, 94, and 103 mW m⁻² for 5, 10, and 15 km thick HPLs, respectively. This highlights the significant role of a thickened heat-producing upper crust in controlling crustal temperatures and surface heat flow.

3.2.3 Nanga Parbat massif (NPM)

For the 60 km thick crust of the NPM, the lower 30 km are assumed to consist of granulitic rocks with radiogenic heat production of 0.4 µW m⁻³ and thermal conductivity of 2.6 $\mathrm{W}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$ (Roy and Mareschal, 2011). The crustal thickness and layer geometry were assigned according to the approximate crustal architecture beneath the modelled terranes, as represented by the grey columns in Fig. 4a. To explore the effect of radiogenic heat distribution in the upper and middle crust, seven models with varying parameters were tested (Fig. S1 and Table S1 in the Supplement).

Figure 4 (a) Crustal-scale N–S cross-section (along XY line in Fig. 5a) of the Nanga Parbat Massif, Kohistan, and Southern Karakoram (modified after Searle et al., 2010). The structure of the crust is based on various cross-sections from the literature (Searle et al., 1989; Bhukta et al., 2006; Burg et al., 2006; Burg and Bouilhol, 2019). Grey columns represent the approximate locations of 1D models. Panels (b)–(d) show 1D steady-state conductive geotherms, and panels (e)–(g) show surface heat flow for different modelling scenarios with varying parameters in the Nanga Parbat Massif (b, e), Kohistan arc (c, f), and Karakoram (d, g).

Models 1–4 examine the effect of decreasing the thickness of the heat-producing upper crust (4 µW m⁻³) while increasing the thickness of the mid-crust (with 2 µW m⁻³). Models 5 and 6 consider a highly enriched upper crust (5.33 µW m⁻³) with thicknesses of 5 and 10 km, respectively. Model 7 assumes a 10 km thick enriched upper crust (5.33 µW m⁻³) combined with a 20 km thick slightly depleted mid-crust. Model 3, representing median parameters and geologically consistent solutions within the explored parameter ranges (rather than best-fitting solution in statistical sense), is selected as the preferred scenario (Table 1).

The value of 5.33 µW m⁻³ represents the enriched end-member of upper-crustal heat production in the NPM. It corresponds to the mean RHP of the high heat-producing gneisses and leucogranites as reported in our in situ gamma spectrometry dataset (Anees et al., 2023). The value 4 µW m⁻³ used elsewhere is the area-weighted average across all NPM lithologies.

The 1D crustal temperature and surface heat-flow profiles for the preferred model are shown in Fig. 4b and e (solid blue line), with variations from the other six models shaded. Temperatures at 10 km depth range between 343 and 506 °C across all models, with 405 °C for the preferred model. Surface heat flow varies between 85 and 120 mW m⁻², with 100 mW m⁻² for the median scenario.

Model 3 is selected on the basis of three criteria: (i) the layer-wise RHP values are consistent with the observed surface geology and crustal structure reported in the literature, and with measured surface heat production from our field gamma spectrometry data (Anees et al., 2023); (ii) the resulting Moho and mid-crustal temperatures fall within ranges compatible with regional metamorphic constraints; and (iii) it represents median values of the explored parameter range, avoiding end-member scenarios that produce extreme results. Similar criteria are also used for selecting preferred models for Kohistan and Karakoram.

3.2.4 Kohistan arc

The Moho beneath the northern Kohistan arc is estimated at ∼70 km depth (Priestley et al., 2019). While the cumulative thickness of the arc rocks is 50–55 km (Jagoutz and Schmidt, 2012; Petterson, 2019), the present-day structural thickness is estimated between 10–25 km (Malinconico, 1986; Petterson, 2019). For modelling purposes, the upper crust is assumed to be 20 km thick (Fig. 4a).

Seven scenarios were tested with variable crustal layer thicknesses and radiogenic heat production (Fig. S2 and Table S2 in the Supplement). In Models 1–5, the upper crust includes felsic (Kohistan batholith) and mafic (amphibolite–granulite) layers, each 10 km thick, with weighted average heat production of 1 µW m⁻³ for felsic rocks and 0.08 µW m⁻³ for mafic rocks (Mukai et al., 1999). The Indian crust, exposed in the NPM to the south, is assumed to underthrust the Kohistan arc and occupy mid-to-lower crustal depths (Searle and Hacker, 2019). The lower crust for all models is similar to the NPM: 30 km thick with heat production of 0.4 µW m⁻³. However, due to uncertainty in parameters of mid-crustal layers (see Fig. 4a), two basement components of Indian origin with variable radiogenic heat production are tested. Thus, the representative 1D thermal column comprises ∼20 km of Kohistan arc crust, a 15 km mid-crustal layer representing underthrust Indian crust, and a 30 km lower crust, giving a total modelled thickness of 65 km.

Model 3, with median parameters within the explored parameter ranges, is selected as the preferred model for the Kohistan arc due to geologically consistent solutions (Table 1). The 1D crustal temperature and heat-flow profiles for the preferred model are shown in Fig. 4c and f (solid green line), with variations from the other six models shaded. Temperatures at 10 km depth vary from 211–413 °C across all models, with 299 °C for the preferred scenario. Surface heat flow ranges from 50–85 mW m⁻², with 65 mW m⁻² for the median model.

3.2.5 Karakoram

To account for present-day uncertainties in the crustal structure of the Karakoram, seven models were tested with varying parameters (Fig. 4a; Fig. S3 and Table S3 in the Supplement). The total crustal thickness is assumed to be 70 km (Priestley et al., 2019), with a lower depleted crust (0.2 µW m⁻³) for all models except Model 2. For the upper crust, weighted average radiogenic heat production values of 2.5 and 1 µW m⁻³ are assigned to the Karakoram Batholith (KB) and the Karakoram Metamorphic Complex (KMC), respectively. Some models also explored the inclusion of fertile Indian crust at mid-crustal levels, which led to exceptionally high Moho temperatures (∼1200 °C), requiring readjustment of radiogenic heat production values. Model 3, representing median parameters and geologically consistent solutions, is selected as the preferred scenario (Table 1).

The 1D crustal temperature and heat-flow profiles for the preferred model are shown in Fig. 4d and g (solid green line), with other models depicted as shaded regions. Temperatures at 10 km depth vary between 274–460 °C across all scenarios, with 385 °C for the median model. Surface heat flow ranges from 65–103 mW m⁻², with 80 mW m⁻² for the median model.

4 1D transient advective-conductive thermal model

This section presents the advective–conductive heat transport formulation, its discretisation, and the exhumation scenarios tested.

4.1 Mathematical formulation

A nonlinear, transient, one-dimensional advective–conductive heat equation with internal radiogenic heat production is solved as:

$$\begin{array}{}\text{(6)}& \mathit{\rho }c\left({\displaystyle \frac{\partial T}{\partial t}}+v_z{\displaystyle \frac{\partial T}{\partial z}}\right)={\displaystyle \frac{\partial }{\partial z}}\left[k(T,z){\displaystyle \frac{\partial T}{\partial z}}\right]+A\left(z\right),\end{array}$$

where z [m] is positive downward, so that upward exhumation corresponds to v_z<0 [m s⁻¹]. T is temperature [K], k(T,z) the temperature-dependent thermal conductivity [$\mathrm{W}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$], A(z) the radiogenic heat production [W m⁻³] with model input values are reported in µW m⁻³, ρ the density [kg m⁻³] and c the specific heat capacity [$\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$].

Advection and conduction are treated by operator splitting on a fixed grid. The prescribed exhumation displacement is accumulated in time; whenever it reaches one grid interval Δz, the temperature field and the associated material properties, such as reference thermal conductivity, radiogenic heat production, density and specific heat capacity, are translated upward by one node. Material leaving the upper boundary is removed, and material carrying the properties of the lowermost layer (Table 2) is introduced at the base. We denote by Tⁿ the temperature field after this remapping step. Advection is therefore represented by fixed-grid material remapping rather than by direct finite-difference approximation of $v_z\partial T/\partial z$. Equation (6) states the physical conservation law, whereas the remapping procedure specifies how its material-transport component is represented numerically.

Table 2 Parameters used for the 1D transient advective-conductive thermal model to test the effects of exhumation on geotherms and heat flow.

The temperature-dependent conductivity is lagged as follows:

$$\begin{array}{}\text{(7a)}& k^n=k(T^n,z),\end{array}$$

then the nodal update:

$$\begin{array}{}\text{(7)}& {\displaystyle }{\mathit{\rho }}_i{c}_{p,i}{\displaystyle \frac{T_i^{n+\mathrm{1}}-{\stackrel{\mathrm{̃}}{T}}_i^n}{\mathrm{\Delta }t}}=\left(\mathbf{K}\right[k^n]{\mathit{T}}^{n+\mathrm{1}}{)}_i+A_i,\text{(8)}& {\displaystyle }(\mathbf{C}-\mathrm{\Delta }t{\mathbf{L}}^n){\mathit{T}}^{n+\mathrm{1}}=\mathbf{C}{\stackrel{\mathrm{̃}}{\mathit{T}}}^n+\mathrm{\Delta }t\mathit{A}\end{array}$$

with

$$\begin{array}{}\text{(8a)}& \mathbf{C}=\text{diag}\left({\mathit{\rho }}_i{c}_{p,i}\right),{\mathbf{L}}^n=\mathbf{K}\left[k^n\right].\end{array}$$

Here ${\stackrel{\mathrm{̃}}{\mathit{T}}}^n$ is the temperature field immediately after the remapping step, $\mathbf{C}=\text{diag}\left({\mathit{\rho }}_i{c}_{p,i}\right)$ is the volumetric heat-capacity matrix and Lⁿ is the tridiagonal conservative conduction operator assembled from the half-station conductivities of Eq. (5), evaluated at the preceding temperature state. Volumetric heat capacity is therefore retained explicitly and varies between layers (Table 2). The temperature field is advanced with backward Euler while the conductivity is lagged, so the scheme is linearly implicit with lagged coefficients rather than a fully nonlinear implicit solve. Dirichlet conditions are imposed at the top and base of the column by modifying the corresponding rows of C and Lⁿ. The system is solved at each timestep with a tridiagonal matrix algorithm.

4.2 Modelling scenarios

To model the transient effect of exhumation on geotherms in a dynamically evolving orogenic belt, the ideal initial condition would represent the pre-collision state. One approach is to consider the present-day cratonic India as an analogue for the pre-collision Indian margin and to simulate subduction of the Indian passive margin and its sedimentary cover beneath Kohistan from 50 Ma onward, thereby increasing crustal thickness. However, uncertainties in the dynamic and kinematic conditions between 50 and 10 Ma before the onset of rapid exhumation make this approach challenging (Zeitler, 1985). It is also likely that, over ∼35 Ma, crustal geotherms re-equilibrated to a near-steady-state condition (Fairley, 2016).

Zircon fission-track data (Zeitler, 1985) show very young cooling ages (<5 Ma) in the NPM and Karakoram, versus >10 Ma in Kohistan, indicating rapid exhumation (∼7 mm yr⁻¹) (Fig. 5a). Alternative interpretations suggest slower exhumation between 50–10 Ma followed by constant rates of 3–4 mm yr⁻¹ over the last 10 Ma (Whittington, 1996). A recent petrochronologic and thermal-modelling study suggested an ultrafast exhumation pulse in the NPM (∼9–13 mm yr⁻¹ beginning at ∼1 Ma), preceded by sustained rapid exhumation of 2–5 mm yr⁻¹ (Guevara et al., 2022). However, this ultrafast exhumation is probably localized to the core of the NPM because river-sediment flux estimates suggest 3–5 mm yr⁻¹ as the regional exhumation rate (Clift et al., 2022).

Figure 5 (a) Zircon fission track cooling ages (after Zeitler, 1985) in the study area. Young cooling ages indicate rapid exhumation in the Nanga Parbat Massif and Karakoram compared to Kohistan. The X-Y line shows the surface profile of the 2D models in Fig. 6. (b) 1D transient advective-conductive geotherms and (c) surface heat-flow profiles showing the effects of variable exhumation rates after 10 Ma for low- (dashed lines) and high-RHP (solid lines) models.

Figure 6 Top: Geological model along the X-Y profile in Fig. 5a, showing the subsurface structure of the major lithological units to 10 km b.s.l. (modified after Searle and Khan, 1996). Centre: Two-dimensional triangular finite-element mesh containing 72 643 elements, created in COMSOL Multiphysics^® for the dynamic thermal modelling. Bottom: Subdivision of the model domain into kinematic blocks with the imposed translation velocity ranges used in Models 1–3.

For this work, we do not consider the ultrafast exhumation scenario and test only the effect of varying regional exhumation rates and RHP on temperature and surface heat flow. Because material is continuously removed at the surface and introduced at the base (Sect. 4.1), the column represents a fixed window through which solid material passes rather than a deforming block: its total thickness remains 25 km, while the fraction occupied by the upper, more radiogenic layer decreases and that occupied by the lower layer increases as exhumation proceeds. Over 10 Ma, this corresponds to 10, 20 and 30 km of throughflow at 1, 2, and 3 mm yr⁻¹, respectively, so that in the two faster scenarios the original upper layer is entirely exhumed before the end of the simulation and the elevated temperatures are sustained by advection rather than by the initial heat-production architecture.

The transient model is intentionally restricted to the uppermost 25 km of the crust and is not a transient simulation of the entire lithosphere. Its purpose is to isolate the sensitivity of the geothermally relevant upper crust to sustained exhumation, at a resolution and computational cost that permit systematic variation of exhumation rate and RHP. The base of the domain is fixed at 800 °C, corresponding to the metamorphic conditions inferred beneath the NPM (Crowley et al., 2009; Guevara et al., 2022), and is consistent with the 578–686 °C obtained at 20 km in the steady-state models of Sect. 3. Holding this boundary fixed assumes that the thermal structure below 25 km is unaffected by exhumation. In a fully coupled calculation, the deeper geotherm would also be advected upward, supplying material warmer than 800 °C to the base of the window; the transient warming reported here is therefore conservative in this respect.

Figure 5b and Table 3 summarize the effects of low and high RHP scenarios with exhumation rates of 1–3 mm yr⁻¹. Geotherms after 10 Ma show significantly higher temperatures than the initial geotherms. For the high RHP models at 3 mm yr⁻¹, the geotherms invert at a depth of 20 km depth. Increasing exhumation rates in low RHP models by 1 mm yr⁻¹ produces geotherms comparable to high RHP models. At 10 km depth, low RHP geotherms increase by a factor of 1.22, while high RHP models increase by 1.14 per 1 mm yr⁻¹ increase in exhumation.

Table 3 Temperatures at 10 km depth (°C) from the 1D transient advective-conductive model for varying exhumation rates over 10 Ma.

Surface heat flow also rises with exhumation. For low RHP models, surface heat flow increases from 75 mW m⁻² (initial) to 100, 160, and 220 mW m⁻² for exhumation rates of 1, 2, and 3 mm yr⁻¹, respectively. High RHP models show an increase from 110 mW m⁻² to 135, 180, and 250 mW m⁻² for the same exhumation rates (Fig. 5c). The highest values in the transient end-member scenarios combining sustained 3 mm yr⁻¹ exhumation over 10 Ma with high RHP represent exploratory upper bounds rather than expected regional conditions, and could apply spatially only to the localized, rapidly exhuming core of the NPM syntaxis, consistent with the spatial confinement seen in the 2D results (Fig. 7).

Figure 7 2D simulation results showing the subsurface temperature distribution and surface heat flow for combinations of exhumation rates and radiogenic heat production scenarios in the Nanga Parbat Massif, Kohistan and Karakoram. Basal temperatures are derived from the 1D transient advective–conductive models.

Because the basal temperature is held fixed while internal heat production continues, material advected upward can exceed the imposed basal value, so that temperatures within the column locally exceed 800 °C and the profile inverts below the resulting maximum. The inverted geotherms obtained in the high-RHP, 3 mm yr⁻¹ case are therefore a joint consequence of the prescribed basal condition and the imposed heat production, and mark the limit of validity of the fixed-boundary, fixed-RHP formulation rather than a predicted physical state.

The Péclet number provides a dimensionless measure of the relative importance of advective and conductive heat transport and can be expressed as $\mathit{\text{Pe}}=vL/\mathit{\kappa }$, where v is the exhumation velocity, L a characteristic length scale, and κ the thermal diffusivity. For the 25 km model domain and $\mathit{\kappa }\approx {\mathrm{10}}^{-\mathrm{6}}$ m² s⁻¹, exhumation rates of 1, 2 and 3 mm yr⁻¹ yield Pe values of approximately 0.8, 1.6 and 2.4, respectively. Thus, advection is comparable to conduction at 1 mm yr⁻¹ and increasingly important at higher exhumation rates, consistent with the progressive amplification of the geotherms.

5 2D thermal model

Two-dimensional thermal modelling was conducted to investigate the lateral distribution of temperature close to the surface. Because the basal temperature boundary at 10 km below sea level (km b.s.l.) is derived from the median 1D transient models, the 2D simulations do not constitute an independent validation of the 1D results but inherit their assumptions. Their added value lies instead in resolving the lateral interaction between adjacent blocks with contrasting RHP and exhumation rates, and the topographic modulation of shallow isotherms, neither of which the 1D models capture; the 2D model is therefore a lateral and topographic extension of the 1D framework. The 10 km depth limit is a deliberate modelling choice reflecting the geothermally relevant upper crust, not a general constraint of the approach.

Two-dimensional thermal modelling was performed with the Solid with Translational Motion formulation in COMSOL Multiphysics^® (COMSOL AB, 2024a, b). Unlike the standard Solid formulation, which is solved in the material frame, this formulation is solved in the spatial frame and adds a translational-motion contribution representing heat transport by the moving solid material. The governing equation is therefore the transient advective–conductive heat equation with internal radiogenic heat production,

$$\begin{array}{}\text{(9)}& \mathit{\rho }{C}_p\left({\displaystyle \frac{\partial T}{\partial t}}+\mathit{v}\cdot \mathrm{\nabla }T\right)=\mathrm{\nabla }\cdot \left[k\right(T\left)\mathrm{\nabla }T\right]+A,\end{array}$$

where ρ [kg m⁻³] is density, C_p [$\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$] specific heat capacity, v [m s⁻¹] the prescribed translational velocity of the solid material, k(T) [$\mathrm{W}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$] the temperature-dependent thermal conductivity and A [W m⁻³] the radiogenic heat production. Petrophysical properties, especially thermal conductivity and RHP, are assigned per lithological domain from the parameters in Table S5 in the Supplement rather than as a function of depth as in Sects. 3 and 4 (Anees et al., 2023, 2024). No fluid-flow heat transport is included; heat advection represents the translation of crustal blocks during exhumation. The equation is solved by the finite-element method. To account for temperature-dependent thermal conductivity, the formulation of Jaupart et al. (2016) was implemented in COMSOL Multiphysics^® through an analytical function. The final triangular finite-element mesh consists of 72 643 elements, covering an area of 2255 km² (Fig. 6).

Boundary conditions include thermally insulated lateral margins and a fixed temperature of 0 °C at the surface. This value corresponds approximately to the elevation-weighted mean annual surface temperature along the profile, whose mean elevation is close to 4200 m. It differs from the 10 °C used in the 1D models (Sects. 3 and 4), where the surface value is a nominal boundary condition whose precise choice has a negligible effect on lithospheric-scale geotherms; in the 2D model the surface temperature enters directly into the shallow isotherm geometry of interest, and a profile-representative value was therefore preferred. At the basal boundary (10 km b.s.l.), laterally variable temperatures were imposed, derived from the median 1D transient models for the NPM, Kohistan, and Karakoram (Table S4 in the Supplement). Three scenarios were evaluated to test the influence of exhumation and RHP.

The modelled geotherms show contrasts in thermal structure depending on the imposed parameters (Fig. 7). In all cases, the 2D models show that topographic relief strongly modifies the geometry of shallow isotherms. This is particularly evident in the NPM and Karakoram, where isotherms are closer to the surface beneath deeply incised valleys and deeper beneath topographic highs.

Model 1 (fast exhumation, low RHP), with basal temperatures exceeding 600 °C at an exhumation rate of 3 mm yr⁻¹, generates the warmest upper crust beneath the NPM. Here, relatively rapid exhumation drives strong advective heat transfer, whose effect exceeds that of the weak radiogenic contribution. The anomaly is spatially confined to the NPM, consistent with focused exhumation.

Model 2 (slow exhumation, high RHP), with basal temperatures ∼520–570 °C beneath the NPM, also produces high upper-crustal temperatures, but the elevated geotherms are maintained by radiogenic heating. The anomaly is broader and less localized than in Model 1, reflecting the distributed radiogenic heat source and subsequent conductive heat transfer.

Model 3 (slow exhumation, moderate RHP), with ∼450–480 °C basal temperatures beneath the NPM, represents an intermediate state. Neither exhumation nor RHP alone generates the extreme anomalies of Models 1 and 2, yielding a more subdued thermal field.

These contrasts are reflected in the modelled surface heat-flow distributions (Fig. 7). At the NPM, values range from ∼120 mW m⁻² in Model 3 to nearly 180 mW m⁻² in Model 1, capturing the strong sensitivity of surface heat flow to exhumation and RHP. By comparison, Kohistan exhibits relatively uniform values (60–90 mW m⁻²) across scenarios, while the Karakoram displays a modest increase to 100–120 mW m⁻². Collectively, the results show that either rapid exhumation or high RHP can independently produce elevated geotherms beneath the NPM, though their effects differ in spatial expression and magnitude.

To quantify the relative importance of exhumation-driven advective heat transport and thermal conduction, the Péclet number, $\mathit{\text{Pe}}=vL/\mathit{\kappa }\left(T\right)$, was evaluated, where v is the exhumation velocity, L [m] is a characteristic vertical length scale, and $\mathit{\kappa }\left(T\right)=k\left(T\right)/\left(\mathit{\rho }Cp\right)$ [m² s⁻¹] is thermal diffusivity. The material properties employed in the model yield volumetric heat capacities (ρCp) of approximately 1.9×10⁶–2.3×10⁶ $\mathrm{J}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{K}}^{-\mathrm{1}}$. Owing to the temperature dependence of thermal conductivity, thermal diffusivity varies spatially, with representative values of approximately 30–40 km² Myr⁻¹. The model extends to 10 km b.s.l. and incorporates topography reaching nearly 8 km elevation in the Nanga Parbat Massif, yielding characteristic vertical length scales between 10 and 18 km. Because the basal boundary is located 10 km b.s.l. and the mean surface elevation along the profile is approximately 4 km, the representative vertical distance between the surface and the basal boundary is ∼14 km. For exhumation rates between 1 and 3 mm yr⁻¹, the corresponding Péclet numbers are approximately 0.3–1.2 for the average-topography case and up to about 0.5–1.8 for the high-relief end-member. These values indicate that conductive heat transfer dominates at slow exhumation rates, whereas heat transport associated with exhumation exceeds conduction in the fast-exhumation scenario, particularly beneath the Nanga Parbat Massif. Because the 2D domain is shallower than the 25 km column used in Sect. 4, the corresponding length scale and Péclet numbers are lower.

6 Crustal differentiation index

The distribution of heat-producing elements in the crust provides an important window into its differentiation history. Perry et al. (2006) proposed a crustal differentiation index (D_I) to indicate vertical stratification in the crust in response to the enrichment of upper-crustal layers. They defined D_I as the ratio of mean surface heat production (A_o) to total mean crustal heat production (A_c);

$$\begin{array}{c}{\displaystyle }{D}_{\mathrm{I}}={\displaystyle \frac{A_{\mathrm{o}}}{A_{\mathrm{c}}}},\\ {\displaystyle }{A}_{\mathrm{c}}={\displaystyle \frac{Q_{\mathrm{s}}-Q_{\mathrm{m}}}{h}},\end{array}$$

where Q_s [mW m⁻²] is the surface heat flow, Q_m [mW m⁻²] is the heat flow from the mantle, and h [m] is the crustal thickness.

In the absence of surface heat-flow data, values derived from the 1D steady-state models are used to calculate D_I. For the NPM, the total mean crustal heat production (1.5 µW m⁻³) is obtained by subtracting the mean mantle heat flow (8 mW m⁻² at Moho in case of preferred steady-state conduction model) from the mean surface heat flow (100 mW m⁻²) and dividing it by total crustal thickness (60 km). The ratio of mean surface heat production (4.6 µW m⁻³) to total mean crustal heat production (1.5 µW m⁻³) yields a differentiation index of approximately 3. Similarly, the D_I for Kohistan and Karakoram are calculated to be 1.1 and 2, respectively.

Because the 1D models are steady-state and conductive, Q_s−Q_m equals the depth-integrated heat production of the assumed crustal column, so A_c is equivalent to the thickness-weighted mean of the RHP values imposed in Table 1. The differentiation index calculated here is therefore the ratio of the measured mean surface heat production to the assumed mean crustal value, and it does not constitute an independent test of the assumed vertical RHP structure. Its value lies in quantifying the degree of depth-dependent depletion required for the modelled geotherms to remain consistent with petrological and geophysical constraints on Moho temperature, expressed in a form directly comparable with published values for other terranes.

Figure 8 shows a plot of differentiation index (D_I) vs. bulk mean crustal heat production (A_c) for the NPM, Kohistan, and Karakoram and compares them with the Indian (Bundelkhand and Dharwar) and North American (Trans-Hudson, Grenville, and Appalachians) cratons (Perry et al., 2006; Podugu et al., 2017). The plot shows a positive correlation between D_I and mean crustal heat production, with D_I>1 suggesting radiogenic radioelement enrichment in the upper crust, possibly due to magmatic differentiation (Perry et al., 2006). The higher D_I in the NPM (3) and Karakoram batholith (2) are comparable to that of the Appalachians (2.5), indicating enrichment of radioelements in the upper crust. In comparison, the Kohistan batholith shows low D_I (1.1), which is on par with the Dharwar craton, Grenville Province and Trans-Hudson Orogen. However, the Kohistan batholith shows higher average crustal heat production compared to these Indian and North American terranes, which can be attributed to the subducting enriched Indian crust under it. It should be noted that, in the absence of surface heat-flow measurements, the estimations of differentiation index and average crustal heat production should be considered preliminary.

Figure 8 Plot of differentiation index (D_I) vs. average crustal heat production (A_c) for the Nanga Parbat Massif, Kohistan and Karakoram batholith in comparison with Indian shields (after Podugu et al., 2017) and North American geological provinces (after Perry et al., 2006). Values for the NPM, Kohistan and Karakoram are derived from modelled surface heat flow, whereas values for the comparison provinces are based on measured heat flow.

7 Discussion

7.1 Thermal modelling of the lithosphere and radiogenic heat production

The tectonic framework of the Himalaya–Kohistan–Karakoram system exerts a profound influence on its present-day thermal state. Subduction and imbrication of the Indian Plate have produced substantial crustal thickening, while subsequent exhumation has advected heat toward the surface. These processes act in tandem: crustal thickening increases the vertically integrated radiogenic heat production of the crust (Jaupart et al., 2016), whereas exhumation rapidly transports hot rocks upward, generating strong near-surface anomalies (Jamieson et al., 2002). Additionally, thickened crust fosters partial melting at depth, with melts emplaced as plutons in the upper crust (Jamieson and Beaumont, 2011). Collectively, these processes contribute to the thick crust, elevated crustal temperatures, shallow Curie depths (∼580 °C isotherm), and anomalously high background heat flow commonly observed in active orogenic belts (Gao et al., 2017).

Developing reliable lithospheric thermal models therefore requires robust constraints on crustal structure and physical properties. In stable cratonic settings, this is feasible because steady-state models can be anchored to deep borehole heat-flow data, which minimize transient near-surface effects (Fuchs et al., 2022). In contrast, active orogens such as the Himalaya lack reliable surface heat-flow datasets and suffer from poorly constrained geophysical profiles at depth. The present study addresses this challenge by adopting a suite of models with variable parameters, thereby bracketing the plausible thermal states of the lithosphere rather than committing to a single deterministic solution.

Within thermal models of the continental crust, RHP is widely regarded as the dominant parameter controlling crustal geotherms (Hasterok and Chapman, 2011). However, its vertical distribution remains ambiguous. It is generally accepted that lower-crustal melting and melt extraction redistribute heat-producing elements upwards in the crust, leading to depletion in the lower crust and enrichment in the upper layers (McKenzie and Priestley, 2016). At mid-crustal levels, this redistribution is neither complete nor uniform. U and Th are commonly co-enriched in evolved granitic melts and are broadly positively correlated across the lithologies in the study area. Both contribute significantly to RHP, and both may be retained variably in resilient accessory phases (e.g., monazite, zircon, allanite) during near- or sub-solidus melting (Alessio et al., 2018; Yakymchuk and Brown, 2019; Weller et al., 2020). Migmatitic gneisses exhumed from mid-crustal levels in the NPM nonetheless show high $\mathrm{Th}/\mathrm{U}$ ratios relative to the associated leucogranites (Anees et al., 2024), illustrating that partitioning behaviour varies among lithologies within a single terrane and confirming the absence of any systematic negative correlation of RHP with increasing metamorphic grade (Hasterok et al., 2018).

The present modelling approach provides a critical evaluation of these hypotheses. Assigning surface-derived RHP values uniformly through the crust produces unrealistically hot geotherms (Figs. 4 and S2), with Moho temperatures far exceeding petrological and geophysical constraints. This strongly suggests that RHP diminishes with depth in active orogens, likely through fractionation, melt extraction, and redistribution processes. At the same time, localized upper-crustal enrichment by pluton emplacement may sustain high surface values. Thus, the results argue against simplistic uniform-depth assumptions and instead favour models that incorporate depth-dependent variability in RHP. The multi-scenario approach constrains the direction of the required depletion but not its precise depth distribution, which remains the principal unresolved uncertainty of the model.

The similarly high differentiation indices of the NPM and the Bundelkhand craton contrast with their markedly different average crustal heat production values, which are considerably higher in the modelled NPM value (1.5 µW m⁻³) than in Bundelkhand (0.41–0.68 µW m⁻³). This difference reflects composition rather than thickness alone: crustal thickening by imbrication and under-thrusting in the NPM stacks fertile Indian upper- and middle-crustal material, raising the mean volumetric heat production of the resulting column. For a given mean volumetric RHP, increasing crustal thickness increases the vertically integrated radiogenic heat production and therefore the crustal contribution to surface heat flow, so a thickened and compositionally fertile column produces elevated temperatures at mid- to lower-crustal depths. These conditions may promote crustal melting and the redistribution of heat-producing elements toward the upper crust, potentially increasing the differentiation index without substantially altering bulk mean crustal heat production if the elements remain within the crustal column. This redistribution may also moderate lower-crustal temperatures (Perry et al., 2006). This is also corroborated by our modelling results, where higher RHP (≥2 µW m⁻³) layers in the mid-crust generated exceptionally hot geotherms (Figs. 4 and S2).

Recent thermomechanical models further emphasize the role of radiogenic enrichment in the upper crust. Fan et al. (2024) demonstrated that sustained ultra-high-temperature metamorphism in the Greater Himalayan Crystalline Complex was primarily driven by elevated radiogenic heat production (∼4 µW m⁻³) in a partially molten upper crust, with mechanical heating playing only a transient role during the early stages of collision. This highlights that crustal differentiation, whereby heat-producing elements are preferentially concentrated in the upper crust through partial melting and magmatic segregation, can maintain high geotherms over tens of millions of years. Such findings are consistent with our results, where high D_I in the NPM and the Karakoram reflect strong upper-crustal enrichment in radiogenic elements, enhancing the thermal gradients.

7.2 Upper crustal heat flow and geothermal implications

1D modelling results indicate that surface heat flow is strongly controlled by the thickness and intensity of heat-producing layers in the upper crust. While transient 1D models emphasize the influence of exhumation, 2D models suggest that its lateral impact is limited compared to radiogenic heat production (RHP). In the NPM, a high-temperature thermal regime results from the combined effects of exhumation and RHP, producing elevated near-surface geothermal gradients. Consequently, the thermal structure beneath rapidly exhuming regions, particularly the Nanga Parbat Massif, cannot be adequately described by purely conductive models. Péclet numbers of approximately 1.6–2.4 in the 1D rapid-exhumation scenarios indicate that advective transport is comparable to or greater than conductive heat transfer and therefore contributes substantially to the upper-crustal thermal structure. Nevertheless, radiogenic heat production remains an important internal heat source controlling the regional thermal field.

The 2D models indicate that topographic relief modulates the near-surface conductive temperature field, while the concentration of hot springs in deeply incised valleys is consistent with additional structural and hydrogeological controls that are not represented in the model. In addition, unlike purely conductive heat transport, topography-driven deep groundwater flow interacting with heat-producing rocks promotes advective heat transfer to the surface. This advective process can transfer heat to the surface more efficiently and rapidly than conduction alone, particularly in evolving orogens (Chamberlain et al., 1995; Wanner et al., 2020). While extreme near-surface gradients have independently been inferred at Nanga Parbat from evidence of shallow boiling hydrothermal systems and Pleistocene anatexis (Chamberlain et al., 1995; Crowley et al., 2009), these likely represent localized anomalies rather than background conditions.

The Himalayan Geothermal Belt hosts several high-temperature hydrothermal systems in which heat is convected by deep circulation of meteoric water. At Yangbajing in Tibet, a shallow reservoir reaches 150–165 °C at only 180–280 m depth, while a deeper reservoir at 950–2000 m attains 250–329 °C (Fan, 2002; Xue et al., 2023). In the Puga Valley of Ladakh, reservoir temperatures of ∼200–250 °C are inferred at depths of 1–2 km (Shanker et al., 2000; Harinarayana et al., 2006). Our models for northern Pakistan indicate modelled surface heat flow locally exceeding 120–180 mW m⁻² and modelled temperatures above 200 °C at approximately 3 km depth (below surface) in major valleys of the NPM and parts of the Karakoram. We emphasise that, unlike the drilled and geophysically constrained values quoted above for Yangbajing and Puga, these are scenario-based model outputs derived without calibration against measured subsurface temperatures. Subject to that distinction, they still indicate a thermal regime comparable to medium- to high-temperature geothermal systems and underline its relevance within the belt. Under the more conservative exhumation and RHP scenarios, modelled isotherms exceeding 100 °C occur at depths of 1200–2500 m (below surface) in the valleys.

A possible indication of fluid interaction with uranium-bearing lithologies comes from radon measurements in regional hot spring waters. Some recent publications have reported spatially variable radon concentrations, with the highest value of 304 Bq L⁻¹ at the Tattapani spring at Raikot Bridge on the western flank of the NPM (Ullah et al., 2022; Muhammad and Haq, 2023). Because ²²²Rn derives from the ²³⁸U decay series and has a half-life of only 3.8 d, such values indicate possible contact with uranium-bearing rock along the shallow portion of the flow path. However, they provide no direct constraint on the contributions of Th or K (which also contribute to RHP) and no information on a deeper heat source. Both studies were undertaken for health risk assessment rather than geothermal characterisation. A dedicated noble-gas survey, in particular of ${}^{\mathrm{4}}\mathrm{He}{/}^{\mathrm{3}}\mathrm{He}$, would be considerably more diagnostic, since radiogenic ⁴He accumulates directly from U and Th decay and would test the crustal origin of the heat inferred from the low ${}^{\mathrm{3}}\mathrm{He}{/}^{\mathrm{4}}\mathrm{He}$ ratios reported for the belt (Hochstein and Zhongke, 1995).

The results indicate favourable thermal conditions, which constitute one of the major components of geothermal resource potential. Moving from thermal favourability to a resource estimate, and ultimately to viable exploration targets, requires additional characterisation of permeability, fracture-network connectivity, in situ stress and hydrogeological conditions, none of which are addressed by the present models. Nonetheless, the deeply incised valleys of the NPM and Karakoram represent areas of elevated thermal favourability that warrant further site-specific investigation.

7.3 Modelling limitations

The thermal models presented here provide first-order estimates of crustal geotherms and surface heat flow, while accounting for uncertainties arising from limited data and heterogeneous parameters. Boundary conditions represent a further source of uncertainty that is only partially assessed. The sensitivity test presented in Sect. 3.2.1 varies the depth of the basal boundary but not the temperature prescribed at it or the type of condition imposed, and it is a local rather than a global test in the sense of Degen et al. (2025): parameters are varied individually about a reference configuration, so correlations between the basal boundary and other model parameters are not resolved. A systematic assessment, varying LAB depth, LAB temperature and boundary type jointly, would establish more firmly which features of the predicted geotherms are independent of these choices. In these models, intra-layer RHP values are assumed to be constant. However, as noted in Sect. 7.1, assigning surface RHP values to mid-crustal layers leads to unrealistically high temperatures. Accordingly, mid-crustal RHP values were adjusted to produce reasonable geotherms consistent with stratified continental crust models (Rudnick and Gao, 2003; Hacker et al., 2015).

The range spanned by the seven scenarios tested for each terrane should be interpreted as an expression of structural uncertainty in the assumed crustal stratification, not as statistical confidence bounds on a fitted solution. The scenarios are constrained by surface radiogenic heat production measurements, geological interpretations, thermochronologically derived exhumation rates and literature-derived crustal structures; however, they are not calibrated against any direct thermal observation, because no surface heat-flow determinations and no deep borehole temperature data exist within the study area. The models are therefore not validated in the formal sense, since validation requires comparison against independent observations of the modelled quantity, and no such observations exist. They are instead assessed for internal consistency and geological plausibility: the requirement that Moho and mid-crustal temperatures remain compatible with petrological and geophysical constraints provides the principal external check, and is the criterion by which the preferred scenarios were selected. The modelled temperatures and heat flow values should accordingly be read as bracketed scenarios rather than deterministic predictions. In this regard, systematic xenolith studies of the exposed intrusions would provide a valuable independent constraint on deep crustal composition and RHP, and would represent an important target for future work.

In active collisional settings, additional dynamic and transient processes can influence crustal temperatures and heat flow, including exhumation, shear heating, partial melting and melt migration, crustal thickening, and paleoclimate effects. Shear heating is significant in fault zones during orogenic evolution (Burg and Gerya, 2005; Wang et al., 2013; Ai et al., 2021), but its crustal-scale contribution is negligible compared to radiogenic heat (Fagereng and Biggs, 2019). Melt migration from mid-crustal anatexis can advect heat to the upper crust, but, in the present case, its effect is spatially limited and neglected. Crustal thickening increases the total radiogenic content but is also excluded due to its long timescale relative to the model duration.

Lateral and vertical lithological variations, along with deformation intensity, strongly influence rock properties and can significantly affect local heat flow (Abdulagatova et al., 2020; Weinert et al., 2020). Thermophysical parameters used in the 2D model, derived from outcrop analogues, do not fully represent deep-crustal conditions affected by temperature, pressure, and fluids. Fluid circulation and rock permeability are not included; cold meteoric water percolation can cause hydrothermal cooling and locally reduce conductive geotherms (Cao et al., 2019). Conversely, in hot spring areas, upward advection of heat by fluids can elevate near-surface temperatures and disturb conductive heat-flow (Bächler et al., 2003; Koltzer et al., 2019). Consequently, regions with high deformation and fracture-related permeability are most susceptible to advective heat transport.

Simulating hydrothermal circulation would require local hydrogeological and structural constraints, such as fault permeability, fracture-network connectivity, and groundwater recharge. These parameters are not available and lie beyond the scope of a regional model. The conductive temperature field presented here therefore provides a background thermal state: it suggests that sufficient heat is available at depth to sustain hydrothermal systems, but not the specific flow pathways by which that heat reaches the surface. The link between the modelled temperatures and the observed springs is accordingly indirect.

8 Conclusion

This study investigated the interplay of radiogenic heat production (RHP) and exhumation in shaping the thermal structure of the Nanga Parbat Massif (NPM), Kohistan, and Karakoram, with a particular focus on their geothermal implications. Given the limited direct heat-flow measurements and the uncertainty in lithospheric structure, we applied multiple steady-state, transient, and two-dimensional thermal models to evaluate a range of plausible scenarios.

The most robust modelling result is that directly extending surface-derived RHP values into the mid- and lower crust produces Moho temperatures inconsistent with the petrological and geophysical constraints considered here. While depth-dependent depletion of heat-producing elements is itself well established, the models quantify the degree of depletion required in a crust thickened and imbricated to this extent, and show that the requirement is insensitive to the particular vertical partitioning adopted. The approach is transferable in method rather than in result: the test of whether uniform downward extrapolation of measured surface RHP produces geologically permissible Moho temperatures can be applied in any orogen where calibrated heat-flow data are unavailable.

The absolute temperatures and heat-flow values reported here, by contrast, depend on the assumed vertical partitioning of RHP, which remains unconstrained by direct observation. The ranges given express structural uncertainty rather than statistical confidence bounds.

Exhumation enhances heat advection by transporting hot material upward, thereby amplifying surface heat flow even where upper-crustal RHP is modest. However, its thermal influence is spatially localized, primarily confined to regions directly above rapidly exhuming blocks such as the NPM. By contrast, variations in RHP exert a broader and more sustained influence on the crustal thermal structure.

Simulation results based on two-dimensional advective-conductive models further suggest the strong imprint of topography on near-surface thermal fields. Isotherms expand beneath topographic highs and are compressed beneath valleys, producing localized zones with high geothermal gradients. Such structural–topographic interactions suggest that deep valleys, particularly in the NPM and Karakoram, represent areas where elevated temperatures may occur at comparatively shallow depths.

Taken together, these findings indicate that the collisional setting of northern Pakistan sustains anomalously high heat flow through the combined effects of radiogenic enrichment, crustal thickening, and rapid exhumation. The thermal regime revealed by our models provides a geodynamic context for the widespread occurrence of hot springs in the region and underscores the associated complexities and uncertainties. Importantly, valleys within the NPM and Karakoram emerge in the models as areas of particularly high thermal favourability. These results provide a screening-level basis for future exploration in a tectonic setting traditionally overlooked in geothermal assessments. Future studies should focus on reservoir, permeability, and hydrogeological conditions to enable resource characterization and assessment.