We thank the reviewer for the positive assessment and for recognizing the timeliness and relevance of our work. We are pleased that the manuscript is considered suitable for publication.

We sincerely thank the reviewer for careful reading of the manuscript and for raising constructive and thoughtful questions. We appreciate the reviewer's recognition of the importance of this work. Below, we provide point-by-point responses to each comment.

Response to point 1

“In Pakistan there are three types of geothermal energy environments, Volcanic in Balochistan, Tectonic in Northern area and Geopressuried   in Southern in Indus Basin. You mentioned here radiogenic means source of radioactive particles present in northern area. Why not you classified area on the basis of geothermal gradient variation.”

Our manuscript focuses on the tectonic (collisional) geothermal environment of northern Pakistan, comprising of the Himalaya–Kohistan–Karakoram terranes. In this region, no reliable regional-scale borehole temperature data exist from which to derive geothermal gradients. In contrast, much of the gradient information in Pakistan comes from sedimentary basins (e.g. Indus Basin) that are tectonically and thermally very different from the high-relief orogenic domains we investigate.

Because of this data limitation, we chose to classify the northern area by tectono‑lithologic domains and their radiogenic heat production (RHP), which can be directly constrained from petrophysical and geochemical measurements, and then to predict geothermal gradients and heat flow with 1D and 2D thermal models. 

Response to point 2

“How you calculate surface heat flow through calculation of surface temperature or used geothermometers.”

In our study, surface heat flow was not measured directly from surface temperatures or geothermometers, but was instead calculated as the output of our thermal models. Specifically, in the 1D steady-state conductive models (Section 3), surface heat flow is computed from the temperature gradient at the uppermost model node using Fourier's Law of heat conduction (q = -k × dT/dz), where k is the thermal conductivity of the surface layer and dT/dz is the computed near-surface temperature gradient.In the 1D transient advective-conductive models (Section 4), surface heat flow additionally includes the advective component arising from upward movement of rock during exhumation. In the 2D model (Section 5), surface heat flow is computed from the spatial gradient of the temperature field at the model surface boundary.

We did not employ geothermometers in this study, as geothermometers (e.g., silica, Na-K, or Na-K-Ca geothermometers) estimate reservoir temperatures from hot spring fluid chemistry and are therefore best applied to characterize individual hydrothermal systems. Our study targets the broader lithospheric-scale thermal structure.

Response to point 3

“What is the HDR potential in your three selected areas, because there is no such granite rocks are present such rocks present in Nagar Parker Pakistan.”

Contrary to this view, granitoid rocks are in fact abundant and well-documented in northern Pakistan. In our study area: (1) the Nanga Parbat Massif exposes Proterozoic Indian basement gneisses and leucogranites with radiogenic heat production of 4–5 μW/m³ , which is well within the range of HDR-suitable granitic basement; (2) the Karakoram Batholith is a 600 km long, up to 30 km wide granitoid intrusion (granodiorites, leucogranites, syenites) with RHP of 2.5 μW/m³; and (3) the Kohistan arc, while predominantly mafic, contains the Kohistan Batholith (felsic granitoid component) with lower but non-negligible RHP of ~1 μW/m³. These rocks are described in detail in our previous studies (Anees et al., 2023 & 2024).

Regarding HDR potential specifically: the NPM and Karakoram Batholith represent the most promising HDR targets due to their high radiogenic enrichment, and surface accessibility in deep river valleys such as the Indus and Hunza gorges.

Response to point 4

“Are you discussing about Himalayan Geothermal Belt in introduction, have you gone through the recent development of geothermal energy exploitation activity in Tibet in China and Puga valley India .What is subsurface temperature there and at what depth geothermal anomaly present and also production potential of these fields and correlate it with geothermal areas of Pakistan.”

Thank you for this suggestion. The Himalayan Geothermal Belt (HGB) extends continuously along the whole Himalayas through Tibet and into Nepal, and host numerous hot springs (with varying surface temperatures) as geothermal manifestations. The physical and chemical nature these hot springs vary significantly depending upon local hydrogeological conditions.  The Tibet and the Puga Valley (Ladakh, India) are among the high temperature hydrothermal systems of the HGB, which represent heat convection by deep circulation of meteoric water.

In Tibet, the Yangbajing/Yangbajain geothermal field has a shallow reservoir at 150–165 °C at only 180–280 m depth, while a deep reservoir at 950–2000 m reaches 250–329 °C. In India’s Puga Valley (Ladakh), studies suggest reservoir temperatures of ~200–250 °C at depths around 1–2 km.

Our models for northern Pakistan predict surface heat flow locally >120–180 mW m⁻² and temperatures exceeding 200 °C at depths of approximately 3 km in major valleys of the Nanga Parbat Massif and parts of the Karakoram, even under conservative assumptions. These values fall in the medium‑ to high‑enthalpy range suitable for power generation, especially if HDR/EGS approaches are used.

In the revised Discussion we will add a short subsection explicitly comparing reservoir depths and temperatures in Tibet and Puga with our modeled values, to highlight the relevance of northern Pakistan within the broader Himalayan Geothermal Belt.

Response to point 5

“Correct the sentence despite numerous Cenozoic intrusions, the absence of active volcanism and low 3He concentrations suggest a primarily crustal origin of anomalous heat.”

We appreciate this remark and propose to revise it to:

'Although the region hosts numerous Cenozoic crustal intrusions — the products of collision-induced anatexis — the absence of active arc or rift volcanism and the low ³He/⁴He ratios in hot spring gases collectively indicate that the anomalous heat is of crustal origin, rather than reflecting a direct mantle contribution.' 

Response to point 6

“Why not you use geothermometers for actual calculation of subsurface temperature, also use oil wells data for geothermal gradient data.”

Chemical geothermometers for hot‑spring waters and oil‑well temperature logs are very useful for site‑specific resource assessment, but they are not ideal constraints on the egional background lithospheric thermal structure that we aim to model.

In northern Pakistan, most deep wells with reliable temperature logs are located in the Indus Basin and other southern sedimentary basins, which are tectonically distinct from the high‑relief collision belt we study here; so using those gradients would therefore mix fundamentally different thermal regimes.​

Hot‑spring geothermometers in the Himalaya–Karakoram commonly reflect complex mixing, boiling and re‑equilibration along deep flow paths; they are strongly influenced by local hydrology and permeability and can differ substantially from purely conductive geotherms.​​

For these reasons, we chose to base our lithospheric models on (i) measured radiogenic heat production and thermophysical properties of the main lithologies, and (ii) geophysical constraints on crustal and lithospheric thickness, and then to predict geotherms and surface heat flow self‑consistently.

Response to point 7

“Heat flow models showing 1200 Centigrade versus 250 km is it possible for economic usage of geothermal energy?”

The temperatures of ~1200–1300 °C in our figures refer to the assumed basal temperature at the lithosphere–asthenosphere boundary (LAB), which we place at 150–250 km depth to control the deep boundary condition of the conductive models. These values are not intended to represent exploitable geothermal resources; they simply define the deep thermal state of the lithosphere in line with standard continental geotherm modeling.​

For economic geothermal exploitation, the relevant depths are far shallower: our models predict temperatures exceeding 100°C at depths of approximately 1200–2500 m in NPM valleys (Section 6.3), and temperatures of 200°C or more within 3 km depth under high-RHP or high-exhumation scenarios. These are entirely within the range of current deep drilling technology (Enhanced Geothermal Systems / EGS currently target 3–6 km depths).

Response to point 8

“According to your conclusion heat energy not produce due to magmatisim but due to overburden and topographic effect, but I think main agent of heat in collision part is heat produce due to magmatisim originated due to collision effect and plat boundaries .Overburden heat generated by geopressuried effect of overburden mostly in south Indus Basin. Please explain it.”

We acknowledge that magmatism has played an important role in the thermal evolution of the Himalaya–Karakoram system, and our manuscript already notes the presence of numerous Cenozoic intrusions. Our main point, however, is that present‑day elevated heat flow and geothermal potential in northern Pakistan can be explained without requiring a currently active magmatic body in the upper crust or shallow mantle. Instead, our models and petrological data indicate that crustal thickening and partial melting have redistributed heat‑producing elements into the upper and middle crust, increasing integrated RHP and thereby raising crustal temperatures over tens of millions of years.

Additionally, rapid exhumation advects this radiogenically heated crust upward, further amplifying near‑surface temperature and heat flow. Finally, topography and focused fluid flow then localize this heat in valleys and along fault zones, giving rise to the observed hot‑spring belts.

In contrast, “overburden” or geopressured effects in the southern Indus Basin are related to thick, low‑permeability sedimentary successions and are indeed a different geothermal play type (geopressured aquifers) than the crystalline‑basement‑dominated system we study here.

Response to point 9

“Please explain heat flow anomalous variation in near collision point and rest of the surrounding areas.”

Our models reveal that heat flow anomalies in northern Pakistan are not uniformly distributed but are strongly controlled by three spatially variable factors: (1) radiogenic heat production (RHP), which is highest in the NPM and Karakoram Batholith and lowest in the mafic-dominated Kohistan arc; (2) exhumation rate, which is fastest in the NPM (2–5 mm/yr) and Karakoram, and slowest in Kohistan; and (3) topography, which creates local focusing of geothermal gradients in deep valleys.

The NPM represents the active culmination point of the India–Asia collision, where Indian basement crust is being rapidly extruded upward. The combination of maximum exhumation rates and high crustal RHP produces predicted surface heat flow of 120–180 mW/m² (our models) which is 2–3 times the global continental average of ~65 mW/m². This is consistent with the widespread hydrothermal activity concentrated around the NPM syntaxis.

In the Karakoram, intermediate to elevated heat flow (65–103 mW/m² steady-state; up to 120 mW/m² with exhumation) is predicted, controlled by the Karakoram Batholith's moderate RHP (2.5 μW/m³) and local rapid exhumation. Hot spring clusters along the Karakoram Fault reflect this thermal anomaly.

Dear Dr. Sebastian Oriolo,

We sincerely thank you for your thorough and constructive review of our manuscript. The comments are highly relevant and will lead to meaningful improvements. We address each comment in detail below.

Response to point 1

“You assume radiogenic heat production based on surface geology and use some assumptions to evaluate their spatial continuity (e.g., Section 3.2.2). Don’t you have any xenolith in exposed intrusions that may be useful to evaluate the subsurface geology?”

This is an excellent point that directly addresses one of the key uncertainties in our modelling approach. We acknowledge that the vertical extrapolation of surface-derived RHP values into the subsurface is a significant source of uncertainty, and that xenolith studies could in principle provide direct petrological constraints on mid-to-lower crustal lithology and composition.

To the best of our knowledge, no systematic xenolith studies focused on radiogenic heat production have been conducted in the NPM, Kohistan, or Karakoram intrusions. The region is geologically complex and poorly accessible, and xenolith populations in the exposed intrusions have not been characterised for their thermophysical or geochemical properties in a way that would constrain crustal RHP at depth. The Kohistan arc, however, exposes a near-complete crustal cross-section from the lower crust to the upper crust. We have used RHP values calculated from this exposed stratigraphy (Mukai et al., 1999) and assigned them to depth-equivalent crustal layers, which is arguably as informative as xenolith constraints for this particular terrane.

For the NPM and Karakoram, in the absence of xenolith data, we rely on the global analogues for the continental crust (Hasterok and Chapman, 2011; Jaupart et al., 2016). We will add a sentence to Section 3.2.2 and to Section 6.4 (Modelling Limitations) explicitly acknowledging that xenolith studies from intrusions in the region would provide a valuable independent constraint on deep crustal composition and RHP, and represent an important target for future work.

Response to point 2

“Some parts of the text need further details on criteria to define model parameters (see comments in the PDF). For instance, when you refer to ‘preferred models’, you have to state clearly why are they considered as such.”

We agree that the selection criteria for ‘preferred models’ are not sufficiently explained in the current manuscript. The preferred models were selected based on a combination of the following criteria:

(1) Consistency with the observed surface geology and crustal structure from the published literature, such that RHP values assigned to each layer are consistent with measured surface heat production from our field gamma spectrometry data (Anees et al., 2023) and with published values for the corresponding lithologies.

(2) Moho and mid-crustal temperatures that fall within ranges compatible with regional metamorphic constraints.

(3) Median parameter selection, so that Model 3 in each set represents median values of the explored parameter range, avoiding the end-member scenarios that produce extreme results.

We will revise Sections 3.2.3–3.2.5 to add one to two sentences after the introduction of each ‘preferred model’ explicitly spelling out these criteria, and will add a brief statement in the Methods section clarifying that the preferred models are median, geologically consistent solutions within the explored parameter ranges.

Response to point 3

“Please revise the use of the terms ‘uplift’ and ‘exhumation’, since in some cases they may be mixed up.”

We thank the reviewer for this important terminological point. In the context of our thermal models, the physically relevant parameter is exhumation, and the rates we cite from thermochronological data (e.g., 2–5 mm/yr in the NPM) are exhumation rates derived from cooling age gradients.

We will carry out a full manuscript-wide review and revise the text and figures accordingly to ensure that ‘exhumation’ is used consistently wherever the thermal model process is described, distinguishing it clearly from ‘rock uplift’ where that term is appropriate in its tectonic context.

Response to point 4

“Lines 337-339: These statements are not totally correct. Both Th and U can be concentrated in magmas, as in the case of A-type magmatism (see Oriolo et al. 2026 J Environm Radioactivity and references therein). In fact, your own geochemical data show a general positive correlation between U and Th (so they are not decoupled), except for some gneisses and granites that show high U with low Th. On the other hand, U may be more relevant for radiogenic heat production than Th, and can also be concentrated in zircon, allanite, etc. In migmatites, leucosomes may have the ‘magmatic’ fingerprint, so they may not be necessarily different to granitoids in general.”

We appreciate this detailed and very useful comment and mostly agree that there is a broadly positive correlation between U and Th across the sampled lithologies. Here, we specifically highlight the challenges of estimating RHP at mid-crustal levels in orogenic settings, where processes such as high-temperature (HT) metamorphism and sub-solidus and partial melting prevail. We point to our previous findings from the NPM (Anees et al., 2024), where migmatite gneisses exhumed from mid-crustal levels have high Th/U ratios compared to the leucogranites. Our intent here is to emphasize that, in the current tectonic setting, assigning RHP values at mid-crustal levels remains highly uncertain.

We will revise lines 337–339 to clarify the case for assigning RHP to mid-crustal levels at the NPM, and add statement regarding U–Th co-enrichment, and their relative contributions of U and Th to RHP.

Response to point 5

“The introduction of crustal differentiation first in Section 6.2 of the Discussion should be revised. Some parts should be perhaps included in the Results.”

We agree that the concept of crustal differentiation and the differentiation index (DI) are partly methodological and therefore better introduced earlier. We will move the definition of the differentiation index, the equations used, and the key numerical DI values for each domain (Nanga Parbat, Kohistan, Karakoram) from Section 6.2 into a new short sub-section at the end of the Results (following the 1D and 2D modelling results). The Discussion section will then focus on the interpretation of these values in the context of crustal evolution and geothermal implications.

Response to point 6

“The authors generally consider crustal anatexis as the main mechanism for magma generation and transfer, but keep in mind that a mantle source with subsequent differentiation may also be possible (e.g., Ding et al. 2025 PNAS).”

We fully agree that mantle-derived magmas and subsequent differentiation have also contributed to the magmatic history of the India–Asia collision zone. In the Karakoram, the presence of mid-Cretaceous subduction-related granodiorites and diorites (Hunza plutonic unit) in the Batholith is itself evidence for magmas with at least a partial mantle contribution during the subduction phase. The post-collisional Miocene syenites in the Baltoro and Kande plutonic complexes may also have more complex source signatures. However, for the purposes of our thermal model, the critical question is not the petrogenetic origin of the intrusions but rather their present-day RHP and their spatial distribution in the crust. Whether heat-producing elements were originally derived from crustal anatexis or from differentiation of mantle-derived melts, the current thermal state of the crust is controlled by where those elements reside today.

We will clarify this by stating that the ‘absence of active magmatism’ refers to the lack of currently active volcanic centres and shallow intrusions, not to the historical role of mantle magmatism during the Cenozoic evolution of the belt. This addition will make clear that our crust-focused thermal models do not exclude a mantle contribution to magmatism; rather, they show that present-day elevated heat flow can be maintained even in the absence of ongoing mantle melt supply. We will also cite Ding et al. (2025, PNAS) in this context.

Response to point 7

“Do you have any information to robustly demonstrate the relationship of hot springs with rocks related to RHP? E.g., high radon flow?”

This is a perceptive and important question. Elevated radon in hot spring in spring waters can indicate interaction of the hydrothermal fluid with uranium-bearing radiogenic rocks. In fact, a couple of studies focused on health risk assessment — Ullah et al. (2021) and Muhammad & Haq (2023) — have measured radon concentrations in hot spring waters at Raikot Bridge (on the western flank of the NPM) and in the Hunza–Nagar valley. They found spatially variable concentrations, with the highest value of 304 Bq/L recorded at the Tattapani hot spring at Raikot Bridge, which is consistent with with radiogenic lithologies of NPM. A future detailed study focusing on radon or noble-gas data would allow a robust, site-specific correlation between fluid chemistry and RHP.

We will revise Section 6.3 to more explicitly discuss this published radon data as partial evidence for the fluid–RHP rock interaction.

We are once again grateful to you for the detailed and constructive engagement with our work. We are confident that the proposed revisions will substantially improve the scientific rigour and clarity of the manuscript.

Best regards,

Muhammad Anees

(On behalf of all co-authors)

We sincerely thank the reviewer for the thorough, technically detailed, and constructive review. We have addressed each point below and revised the manuscript accordingly.

MAJOR COMMENTS

1. Geothermal implications versus geothermal resources

"The manuscript  discusses “geothermal potential”, “geothermal resources”, and geothermal exploration targets. However, the presented modelling framework only evaluates the temperature field and associated conductive heat flow. For sure, temperature is a necessary prerequisite for geothermal utilisation. But in practice, other aspects matter even more:

Geothermal development in practice depends a lot on additional parameters including permeability, reservoir (rock) quality, fluid availability, connectivity of fracture networks, stress state, etc. - usually captured in (resource-wide) estimates of extractable heat. These aspects are not the central aspect in the present study. Consequently, several statements in the Discussion and Conclusions appear stronger than I would see justified by the presented study.

The authors should consistently distinguish between: thermal favourability, geothermal resource estimation, and viable geothermal exploration targets.

At present, the manuscript demonstrates elevated temperatures and heat flow scenarios as an important pre-requisite, Suggestion: reframe in this context.”

We acknowledge that our models evaluate only the thermal field and do not address permeability, reservoir quality, fluid availability, fracture connectivity, or stress state, all of which are essential for geothermal resource estimation. The present study therefore demonstrates thermal favourability as a first-order prerequisite rather than a geothermal resource assessment.

To bring more clarity, we have revised the terminology throughout the Abstract, Discussion (Section 6.3), and Conclusions.

Specifically:

1. statements such as “prime targets for geothermal exploration” and “significant geothermal potential” have been reworded in terms of thermal favourability, e.g., “regions of elevated thermal favourability that warrant further site-specific investigation”;
2. we have added an explicit statement in Section 6.3 clarifying that the transition from thermal favourability to resource estimates requires additional characterisation of permeability, stress state, and hydrogeology, none of which are addressed by our models; and
3. the Conclusions have been softened accordingly, presenting the modelled temperature fields as a screening-level basis for future exploration, rather than as demonstrated resources.

2. Quantitative confidence versus available constraints

"The authors present thermal scenarios, but the degree of quantitative confidence sometimes exceeds what is justified by the available data. This is particularly important because the study area lacks direct heat-flow measurements and deep temperature constraints.

The thermal models are primarily constrained by: surface radiogenic heat production measurements, geological interpretations, thermochronological exhumation estimates, and literature-derived crustal structures.

Validation is largely qualitative and relies on geological plausibility, cooling ages, and the occurrence of thermal springs.

The manuscript would benefit from a more explicit discussion of the limited calibration basis, as also mentioned by the previous review comments (so, I’ll leave more details)."

We agree. The models are constrained by surface RHP measurements, geological interpretations, thermochronological exhumation rates, and literature-derived crustal structures owing to the absence of direct borehole heat-flow and temperature data. Validation is therefore qualitative (geological plausibility, cooling age compatibility, hot spring occurrence).

In the revised manuscript, we have added a dedicated paragraph to Section 6.4 (Modelling Limitations) explicitly detailing:

• the data constraints underpinning the models;
• the absence of direct heat-flow calibration; and
• what this means for interpreting the predicted temperatures and heat-flow values as bracketed scenarios rather than deterministic predictions.

We have also revised the abstract and conclusions to reflect this more clearly, replacing confident quantitative statements with appropriate uncertainty ranges and conditional language.

3. Constraints on subsurface RHP distribution

"As also already highlighted by previous reviewers, the largest uncertainty in the modelling framework concerns the vertical distribution of radiogenic heat production.

The authors correctly acknowledge that: surface-derived RHP values cannot be directly extrapolated to depth, no suitable xenolith constraints exist, and global analogues are therefore used.

This is a reasonable approach, but it does not resolve the underlying uncertainty. The predicted temperature fields remain fundamentally controlled by assumptions regarding crustal stratification and the partitioning of RHP with depth.

This limitation should be emphasized more strongly throughout the manuscript."

We agree, and we appreciate that the reviewer recognises our approach to be reasonable given the data situation. As discussed in our response to Dr. Oriolo’s community comment, no xenolith constraints on deep crustal RHP exist for the region, and only the Kohistan arc offers a (partial) substitute through its exposed crustal cross-section. For the NPM and Karakoram, the vertical RHP partitioning is necessarily assumption-driven.

We have strengthened this point in three locations in the revised manuscript:

1. at the end of Section 3.2.2, where the layered RHP scenarios are introduced, we have stated explicitly that the vertical partitioning of RHP is the dominant model uncertainty;
2. in Section 6.1, the discussion of depth-dependent RHP has been reframed to emphasise the multi-scenario approach and the resulting uncertainty;
3. in Section 6.4 and the Conclusions, we have highlighted that the range spanned by the seven scenarios per terrane should be read as an expression of this structural uncertainty rather than as statistical confidence bounds.

4. Extremely high predicted temperatures and heat-flow values

"Several model scenarios predict: surface heat flow approaching 180–250 mW m⁻², temperatures exceeding 600 °C at 10 km depth, and even inverted geotherms.

These values appear high for a non-magmatic collisional setting. Please discuss more explicitly: whether these values are physically realistic, whether they represent transient end-member scenarios, and whether they should be interpreted as localized anomalies rather than regional background conditions.

The distinction between exploratory end-member simulations and plausible present-day thermal states is currently not sufficiently clear. Also interesting could be a (non-dimensional) Pe-number analysis to compare and relate effects of diffusion vs. Advection (through exhumation) more fundamentally (also as a basis for a more solid comparison to other regions)."

This is a very important point and we agree that the distinction between exploratory end-members and plausible present-day states must be made much clearer. The highest values, such as temperatures exceeding 600 °C at 10 km, surface heat flow of 220–250 mW m⁻² and inverted geotherms, arise only in the transient end-member scenarios combining sustained exhumation of 3 mm y⁻¹ over 10 Ma with high radiogenic heat production. These scenarios are exploratory upper bounds, not expected regional conditions. We have clarified this in Sections 4.2 and 6.3:

1. these scenarios are exploratory upper bounds, not expected regional conditions;
2. spatially, such values could only apply to the localized, rapidly exhuming core of the NPM syntaxis, consistent with the spatial confinement seen in the 2D results (Fig. 7); and
3. the moderate scenarios (≤120 mW m⁻²) are our preferred representation of the regional thermal state.

However, we would like to point out that the extreme near-surface thermal gradients have been inferred at Nanga Parbat from evidence of shallow boiling hydrothermal systems and Pleistocene anatexis (Chamberlain et al., 1995; Crowley et al., 2009), but we agree these likely represent localized anomalies rather than background conditions.

We thank the reviewer for the excellent suggestion of a Peclet-number analysis. We have added it to Section 4 (with a brief cross-reference in Section 5). For a characteristic length scale L = 25 km (model crustal column), thermal diffusivity κ ≈ 10⁻⁶ m² s⁻¹, and exhumation velocities of 1–3 mm y⁻¹, the Peclet number Pe = vL/κ ranges from approximately 0.8 to 2.4. This quantifies why exhumation at 1 mm y⁻¹ produces only moderate perturbation of the conductive state (Pe < 1, diffusion-dominated), whereas 2–3 mm y⁻¹ places the system in the advection-influenced regime (Pe > 1), explaining the strong amplification of geotherms and the emergence of inverted temperature profiles in those runs. 

5. Two-dimensional model limitations

"The limitations of the 2D modelling framework deserve a more details discussion, as it: extends only to 10 km depth, inherits its lower boundary conditions from the 1D models, neglects hydrothermal circulation, and represents exhumation using simplified translational motion.

As a result, the 2D simulations are not an independent validation of the 1D results but largely inherit their assumptions. This dependency should be discussed more explicitly.

Furthermore, hydrothermal circulation is not included despite the fact that thermal springs constitute one of the main observational motivations for the study. Consequently, the connection between the simulated conductive temperature fields and observed thermal springs remains indirect. Please clarify."

The 2D model derives its basal (10 km bsl) temperature boundary condition from the median 1D transient models and therefore cannot be regarded as an independent validation of the 1D results. Its added value lies instead in resolving the lateral interaction between adjacent blocks with contrasting RHP and exhumation, and the topographic modulation of shallow isotherms, neither of which the 1D models capture. We have revised the opening of Section 5 and Section 6.4 to state this dependency explicitly and to reposition the 2D model as a lateral/topographic extension of the 1D framework rather than an independent check.

The 10 km depth is intentional and appropriate for the geothermally relevant upper crust, but we agree this limit should be explicitly stated as a modelling choice rather than a general limitation.

Simulating hydrothermal circulation in the 2D model requires additional local hydrogeological and structural constraints (fault permeability, fracture network connectivity, groundwater recharge, etc.) that are not available and also not within the scope of the present regional model. We already acknowledge this in Section 6.4 but have strengthened the discussion to clarify that the conductive temperature field provides a background thermal state. The link between our modelled temperatures and observed springs is therefore indirect as we demonstrate sufficient heat at depth to sustain hydrothermal systems, but not the specific flow pathways.

SECTION-SPECIFIC COMMENTS

Comments on the 1D steady-state model (Section 3)

"The governing equation is correctly nonlinear due to the temperature dependence of thermal conductivity. However, the manuscript would benefit from clarifying why a numerical finite-difference solution is required for what is otherwise a classical 1D conductive geotherm problem for which analytical or semi-analytical solutions exist under common assumptions. (I assume the main reason is the T-dependent thermal conductivity?)

In addition, the implications and physical validity of assuming steady-state conductive conditions in a rapidly evolving collisional orogen should be discussed more critically, especially with respect to the transient solution in the next chapter: why is the steady-state solution needed, at all?

I also encourage the authors to use standard differential notation throughout the mathematical formulations. The current notation is somewhat unconventional and reduces readability.

The choice of Dirichlet boundary conditions at the base of the lithosphere deserves further justification. In lithospheric thermal modelling, basal heat-flow conditions are often preferred because they may be more robust than prescribing a fixed temperature at a poorly constrained depth. The sensitivity of the model to this assumption should be discussed. For a broader discussion on the influence of boundary conditions and model trustworthiness, see Degen et al. (2025, Solid Earth, 16, 477–502)."

(1)  Finite-difference approach:

The motivation for the finite-difference solution is not solely the nonlinear temperature dependence of conductivity. The geological model contains abrupt lithological contrasts in conductivity and heat production. We therefore discretise the conservative operator

∂/∂z (k ∂T/∂z)

directly using coefficients evaluated at half stations. This avoids first applying the product-rule expansion kT_zz+k_z T_z, which requires a differentiable conductivity field and is not classically valid at abrupt material interfaces. The conservative formulation requires no numerical derivative of k and treats smooth, layered, and piecewise-discontinuous coefficient distributions within the same numerical framework. The same formulation also accommodates the nonlinear update of k(T) and alternative basal boundary conditions. We have added a brief statement to Section 3.1 making this explicitly clear.

(2)  Justification for Steady-state models: 

The steady-state models serve three purposes that we now state explicitly at the start of Section 3:

• they provide the reference conductive state against which the transient exhumation effect can be isolated and quantified. So, without the steady-state baseline, the increment factors in Table 3 would have no meaning;
• they supply physically consistent initial conditions for the transient runs (Section 4.2), which is preferable to arbitrary initial geotherms given the evidence that crustal geotherms had largely re-equilibrated over ~35 Ma prior to the onset of rapid exhumation (Fairley, 2016); and
• they permit an efficient sensitivity exploration of structural parameters (LAB depth, HPL thickness, layer RHP) that would be computationally and conceptually redundant to perform in transient mode.

We have included a clarification that in the rapidly exhuming NPM the steady-state assumption is clearly violated at present, which is precisely why the transient and 2D models are introduced.

(3)  Differential notation: 

We agree the current notation is unconventional. We have replaced it with standard partial differential notation (e.g., ∂T/∂t, ∂/∂z(k∂T/∂z)) throughout Sections 3 and 4.

(4)  Basal boundary condition: 

Our rationale for the Dirichlet condition is that the LAB is, by definition, an isotherm-like rheological boundary commonly associated with the 1300 °C mantle adiabat (McKenzie et al., 2005), so prescribing its temperature is physically natural, whereas the basal (mantle) heat flow beneath the study area is itself unmeasured and would have to be assumed with at least as much uncertainty. In cratonic settings, Moho or basal heat-flow conditions are indeed preferable because they can be anchored to measured surface heat flow (e.g., Kumar et al., 2007, 2009); in our case no such anchor exists, which is why we opted for the LAB isotherm.

Importantly, the sensitivity of the results to this assumption is partly addressed by the LAB-depth test in Section 3.2.1 (150–250 km), which shows that a 100 km variation in the position of the fixed-temperature boundary has a negligible effect on upper-crustal temperatures and surface heat flow. This implies that the upper-crustal results are robust against reasonable variations of the deep boundary condition, since the crustal thermal field is dominated by internal heat production rather than basal heat input. We have made this connection explicit in Section 3.2.1, by adding a sentence acknowledging the alternative Neumann (basal heat flow) formulation and its trade-offs, and cited Degen et al. (2025) for the broader discussion of boundary-condition influence on model reliability.

Comments on the transient model (Section 4)

"The formulation of the transient model requires substantial clarification.

The chapter title refers to an “advective-conductive” thermal model. However, the governing equation presented in Section 4.1 only contains conductive heat diffusion and internal heat production. No advection term appears in the governing equation. This is problematic because exhumation is the central process investigated in the transient simulations and Table 3 clearly suggests that advective heat transport is included through exhumation rates.

The governing equation should therefore explicitly include the material advection term associated with exhumation and clearly define the adopted sign convention.

Furthermore: the units of thermal conductivity appear to be incorrect, the transition from Eq. (6) to Eq. (7) seems to omit the ρc term, and the numerical scheme appears to be semi-implicit (or Picard-linearized) rather than fully implicit if the coefficient matrix is evaluated using the previous time step."

(1) Formulation of the transient model:

We agree that Equation (5) was incomplete as presented. The physical problem corresponds to

ρc(∂T/∂t+v_z  ∂T/∂z)=∂/∂z [k(T,z)  ∂T/∂z]+A,

where z is positive downward and upward exhumation therefore has v_z<0.

In the numerical implementation, advection and conduction were treated by operator splitting. The prescribed exhumation displacement was accumulated in time. Whenever it reached one spatial grid interval, the temperature field and the associated material properties (such as thermal conductivity, radiogenic heat production, density, and specific heat capacity) were translated upward by one grid cell. Material leaving the upper boundary was removed and material with prescribed properties was introduced at the lower boundary. The translated state was then advanced through the conductive-radiogenic timestep. Thus, advection was implemented by fixed-grid material remapping rather than by direct finite-difference approximation of v_z ∂T/∂z. We have revised Section 4.1 to state this explicitly.

(2) Units of thermal conductivity:

The units of thermal conductivity in Section 4.1 are given as W m⁻³ K⁻¹, which is a typographical error; the correct units are W m⁻¹ K⁻¹, as stated correctly in Section 3.1. This has been corrected.

(3) Missing ρc term:

We agree that the transition to the matrix notation was insufficiently explicit. In the implementation, both the conductive operator and the radiogenic source term are divided by volumetric heat capacity ρc. We have retained ρc explicitly in the revised derivation and defined the matrix operator and source term dimensionally.

(4) Semi-implicit vs. fully implicit:

The temperature field is advanced using backward Euler, while the temperature-dependent conductivity is evaluated from the preceding temperature state. The method is therefore more precisely described as a linearly implicit, lagged-coefficient scheme rather than a fully nonlinear implicit scheme. We have corrected the terminology in Section 4.1.

Comments on the 2D thermal model (Section 5)

"The statement that a “translational term” is added to the heat conduction equation requires clarification. If this term represents material advection due to exhumation, then the governing equation is no longer purely conductive but rather a heat transport equation.

Because the interplay between thermal diffusion, exhumation-driven advection, and radiogenic heat production is central to the manuscript, I suggest including a Peclet-number analysis. Such an analysis would provide a useful quantitative framework for evaluating the relative importance of conductive and advective heat transport in the investigated scenarios.

In addition, Figure 7 would benefit from an additional panel showing the imposed exhumation/advection field and the geometrical subdivision into the different tectonic blocks."

We thank the reviewer for this important comment. We have revised Section 5 to clarify that exhumation was implemented using the Solid with Translational Motion formulation in COMSOL Multiphysics®. This formulation is solved in the spatial frame and introduces heat transport by moving solid material. Accordingly, the governing equation corresponds to an advective–conductive heat-transport equation with internal radiogenic heat production. No fluid-flow advection was considered.

Following the reviewer's suggestion, we also added a Peclet-number analysis for 2D models to quantify the relative importance of heat transport associated with exhumation and thermal diffusion. For the 2D domain, the relevant length scale is the model depth extent (L ≈ 14 km, i.e., 10 km bsl plus mean topographic elevation) rather than the 25 km crustal column used in the 1D transient analysis; the resulting Peclet numbers are therefore correspondingly lower. Calculated Peclet numbers for Model 2 and 3 are < 1, while Model 1 has Pe > 1 for areas with exhumation of 3 mm y⁻¹. This indicates that conductive heat transfer dominates under low exhumation rates, whereas heat transport associated with exhumation exceeds conduction at high-exhumation scenarios, particularly beneath the Nanga Parbat Massif.

We have included an additional panel in Figure 6 as per the reviewer's suggestion to include the geometrical subdivision into the different tectonic blocks together with the imposed exhumation (translational velocity) field. Since Figure 6 already presents the model setup (geological model and mesh), the imposed velocity field belongs there, keeping Figure 7 purely results. 

Discussion and conclusions

"The manuscript successfully demonstrates the importance of: radiogenic heat production, crustal differentiation, exhumation, and topographic effects for shaping the thermal structure of northern Pakistan.

The result that surface RHP values cannot simply be extrapolated to the middle and lower crust is particularly important and can be highlighted.

The topographic effects identified in the simulations are physically plausible and consistent with previous studies, although they do not constitute a fundamentally new insight. However, their application to the specific geological setting of northern Pakistan is  valuable.

The geothermal implications remain somewhat vague because only the temperature field is analysed. Statements regarding petrothermal resources, geothermal targets, and exploration opportunities should be revised to reflect the fact that reservoir properties, permeability, and hydrogeological conditions have not been investigated."

We thank the reviewer for this assessment. In line with our response to Major Comment 1 above, all statements on petrothermal resources, geothermal targets, and exploration opportunities in Sections 6.3 and 7 have been revised to the level of thermal favourability, with an explicit acknowledgment that reservoir properties, permeability, and hydrogeological conditions were not investigated and are prerequisites for any resource classification. Following the reviewer’s suggestion, we have also given greater prominence, in both the Abstract and the Conclusions, to the finding that surface RHP values cannot be extrapolated to the mid- and lower crust, as it is the most robust and transferable outcome of the study. 

Minor comments

- "Please clarify the basis for the “preferred models” and ensure that the rationale is explicitly stated in the manuscript. I agree with the concerns raised previously on this point."

As detailed in our response to the community comment by Dr. Oriolo (Comment 2), the preferred models (Model 3 per terrane) are median-parameter scenarios selected for

• consistency of layer RHP with our measured surface heat production (Anees et al., 2023) and published lithology-specific values,
• Moho and mid-crustal temperatures compatible with regional metamorphic and geophysical constraints, and
• avoidance of end-member extremes.

These criteria have now been stated explicitly in Sections 3.2.3–3.2.5 and summarised in the Methods.

- "Section 3.2.3 (NPM): thermal conductivity values appear to be missing for one unit."

We acknowledge this typographical omission at line 167, where the sentence reads “thermal conductivity of Wm⁻¹K⁻¹” with the numerical value missing. The lower-crustal granulitic rocks were assigned a thermal conductivity of 2.6 W m⁻¹ K⁻¹, as listed in Table 1. The missing value has been inserted in the text.

- "Section 3.2.3: please explain the origin of the value 5.33 μW m⁻³ used for radiogenic heat production."

The value of 5.33 μW m⁻³ represents the enriched end-member of upper-crustal heat production in the NPM, derived from our in-situ gamma spectrometry dataset (Anees et al., 2023): it corresponds to the mean RHP of the most enriched lithological group exposed in the massif (the high heat-producing gneisses and leucogranites), as opposed to the 4 μW m⁻³ used elsewhere, which is the area-weighted average of all NPM lithologies. Models 5–7 use this enriched value to test the effect of an upper crust dominated by the most radiogenic units. We have added a sentence to Section 3.2.3 stating this derivation explicitly, with the reference to Anees et al. (2023).

- "Please ensure consistent usage of the terms “uplift” and “exhumation” throughout the manuscript."

As committed in our response to the community comment, we have conducted a full manuscript-wide terminology review to ensure consistent and correct usage throughout.

- "Consider adopting more standard differential notation in the equations. "

As noted above, we have revised all mathematical formulations in Sections 3 and 4 to use standard partial differential notation.