the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
| 17 Nov 2025
Geothermal implications of the lithosphere’s thermal structure in northern Pakistan
Abstract. 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 to the Nanga Parbat Massif (NPM), Kohistan arc, and Karakoram terrane resulting from translation of heat conduction due to exhumation of blocks. Steady-state results show strong dependence of geotherms on crustal radiogenic heat production (RHP): in the NPM, upper-crustal enrichment (4–5 μWm−3) yields surface heat flow of 85–120 mW m−2, whereas Kohistan produces lower values (50–85 mW m−2) due to its mafic-dominated crust. Karakoram yields intermediate heat flow (65–103 mW m−2), with RHP concentrated in the batholith and metamorphic complexes. ID transient exhumation models demonstrate that uplift rates of 2–3 mm y−1 in the NPM can further amplify geotherms, producing surface heat flow up to 220–250 mW m−2 and inverting deep geotherms at 20 km when RHP is high. Two-dimensional thermal simulations 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−2 (moderate scenarios) to nearly 180 mW m−2 (high exhumation). Crustal differentiation indices further indicate strong upper-crustal enrichment in the NPM and Karakoram, indicating the redistribution of heat-producing elements during crustal thickening and partial melting. The models demonstrate that the region can sustain anomalously high heat flow through the interplay of RHP, exhumation, and crustal differentiation. For northern Pakistan, this provides a robust geoscientific basis for understanding the origin of widespread hydrothermal activity and underscores the region’s significant geothermal potential, positioning it as a promising target for future exploration and sustainable energy development.
- Preprint
(1815 KB) - Metadata XML
-
Supplement
(491 KB) - BibTeX
- EndNote
Status: open (until 18 Apr 2026)
- RC1: 'Comment on egusphere-2025-5252', Tariq Feroze, 07 Feb 2026 reply
-
RC2: 'Comment on egusphere-2025-5252', Anonymous Referee #2, 10 Mar 2026
reply
Geothermal implications of the lithosphere’s thermal structure in northern Pakistan
I have read the article above mentioned article about geothermal energy potential of Pakistan through different advanced methodology and find out some questions related to methodology and geology of northern areas of Pakistan. I really appreciate the authors to a hot current energy related topic for meeting future energy problems of Pakistan .I also thankful to editor for sending such a knowledge full research article for review.
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.
- How you calculate surface heat flow through calculation of surface temperature or used geothermometers.
- 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.
- 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.
- Correct the sentence despite numerous Cenozoic intrusions, the absence of active volcanism and low 3He concentrations suggest a primarily crustal origin of anomalous heat.
- Why not you use geothermometers for actual calculation of subsurface temperature, also use oil wells data for geothermal for geothermal gradient data.
- Heat flow models showing 1200 Centigrade versus 250 km is it possible for economic usage of geothermal energy.
- 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.
- Please explain heat flow anomalous variation in near collision point and rest of the surrounding areas.
Citation: https://doi.org/10.5194/egusphere-2025-5252-RC2 -
AC1: 'Reply on RC2', Muhammad Anees, 11 Mar 2026
reply
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.
Citation: https://doi.org/10.5194/egusphere-2025-5252-AC1
Viewed
| HTML | XML | Total | Supplement | BibTeX | EndNote | |
|---|---|---|---|---|---|---|
| 314 | 171 | 34 | 519 | 66 | 19 | 23 |
- HTML: 314
- PDF: 171
- XML: 34
- Total: 519
- Supplement: 66
- BibTeX: 19
- EndNote: 23
Viewed (geographical distribution)
| Country | # | Views | % |
|---|
| Total: | 0 |
| HTML: | 0 |
| PDF: | 0 |
| XML: | 0 |
- 1
The paper handles a very updated and recent topic. Accepted for publicaiton.