GIS-, geochemistry-, and geophysics-based identification of geothermal energy resources in sedimentary basins: a case study of Zindapir Anticlinorium Zone, Eastern Sulaiman Fold and Thrust Belt, Pakistan

Globally, fold and thrust belts are considered promising zones for geothermal energy exploration due to their complex structural geometries. These zones provide favorable conditions that contribute to the trapping and circulation of geothermal resources. However, their inherent structural complexity poses significant challenges in delineating geothermal reservoirs. This study analyzes the geothermal reservoir potential of the Zindapir Anticlinorium Zone, which is part of the Eastern Sulaiman Fold and Thrust Belt in Pakistan. The objectives of this research were achieved by utilizing an integrated approach that involved remote sensing and geochemical, gravity, magnetic, seismic, and petrophysical analyses. The remote sensing data indicated that regional and local tectonic activity, along with the uplift of Eocene and Paleocene rocks, play a significant role in the development of geothermal reservoirs. Additionally, higher land-surface temperature anomalies, active seismicity, and high-gradient streams passing through hot springs are key indicators of geothermal resource potential. The higher ionic concentration of Na⁺ in the geochemical analysis indicates a deep circulation of water that affects thermal activities. Temperature estimates suggest that the Gulki hot spring lies in a zone characterized by immature thermal water, whereas the hot springs of Jaffarabad and Zindapir originate from subsurface zones with partially equilibrated water. Interpretations of gravity and magnetic data suggest that the anticlines of the Zindapir Anticlinorium Zone have less sedimentary thickness and favorable heat flow conditions. Moreover, seismic and petrophysical analyses indicate that the selected reservoir formations (Pab Sandstone, Lower Goru Formation, and Chiltan Limestone) have the potential to act as geothermal reservoirs. The comprehensive evaluation index (CEI), which incorporates temperature, porosity, and permeability parameters, suggests that the Dhodak, Affiband, Rodho, and Zindapir anticlines are ideal locations for the extraction of economical geothermal reservoirs in the Zindapir Anticlinorium Zone. These anticlines have less sedimentary thickness and compressional tectonic deformation, which allow the up transfer of heat from deeper sources. It is suggested that the workflow of this integrated approach can be implemented in any sedimentary basin to identify potential economic locations of geothermal reservoirs.

1 Introduction

Worldwide interest in the exploration and development of unconventional renewable clean energy sources is increasing due to environmental degradation and the depletion of conventional oil and gas reserves (Abbasi and Abbasi, 2012; Kurek et al., 2021; van der Meer et al., 2014; Rezaie and Aghajani, 2013). Geothermal energy is a naturally occurring resource that can be extracted without producing chemical pollutants or waste in its emissions (Lund and Toth, 2021). Across the world, operational geothermal power plants produce 16,318 MW from nearly 198 fields (Bertani, 2012). The United States, Indonesia, the Philippines, and Turkey are the main such countries, and they produce 96,552 GWh of electricity from geothermal energy. In addition, New Zealand, Mexico, Kenya, Italy, Iceland, and Japan have also developed their geothermal capacity (Gutiérrez-Negrín, 2024). These operational geothermal fields mostly use resources from hydrothermal and conventional reservoirs via nearly 3,700 wells (Gutiérrez-Negrín, 2024). However, to meet the urgent objective of limiting the rise in average global temperature below the 1.5 °C threshold, the expansion of geothermal sources should be considered (Gutiérrez-Negrín, 2024).

Geothermal sources are commonly classified based on hydrothermal activities, geopressured zones, and hot dry rock mechanisms developed in the subsurface. The boundaries of tectonic plates are considered favorable locations for the exploration of economic geothermal reservoirs (Saemundsson et al., 2009). Most high-enthalpy geothermal resources are found along seismic belts linked to crustal weakness zones, including plate boundaries and volcanic hotspots (Birkle and Bundschuh, 2007; Soltani et al., 2021). Fractures and faults are mostly developed along structurally weak zones due to seismic activities, creating pathways that facilitate the movement, retention, and accumulation of geothermal energy (Wang et al., 2021; Wang and Kuang, 2022; Wu and Zhou, 2023; Morata et al., 2023).

The Geographic Information System (GIS) and remote sensing (RS) are highly efficient and cost-effective approaches for assessing the suitability of geothermal exploration. For identifying subsurface geothermal reservoirs, surface features such as anomalous land surface temperature (LST), hot springs, lithological distributions, seismicity, and tectonic activities can be systematically integrated using the Analytical Hierarchy Process (AHP), a widely used method within GIS for combining multiple geospatial parameters. The AHP is a valuable decision-making tool that effectively integrates both qualitative and quantitative data in a structured and hierarchical framework (Li et al., 2023). In regions that lack heat flow data, geochemical methods such as helium isotope ratios and silica concentrations can be used to indirectly estimate subsurface temperatures and the origin of hot springs (Swanberg and Morgan, 1978; Sano and Wakita, 1985).

Gravity, magnetic, and seismic techniques are geophysical methods that can be utilized to detect subsurface geothermal potential by revealing the spatial distribution of various geological features (Gessner et al., 2016; Xu et al., 2019; Xu et al., 2021). Gravity can be used to simulate geothermal systems by analyzing crustal variations which help interpret subsurface structures that are favorable for geothermal energy accumulation (Hinze et al., 2013). Magnetic methods are helpful for quantifying the Curie point depth (CPD), geothermal gradient (GG), and heat flow (Q) (Bhattacharyya and Leu, 1977; Blakely, 1988; Bouligand et al., 2009; Espinosa-Cardeña and Campos-Enriquez, 2008; Ross et al., 2004; Shuey et al., 1977; Spector and Grant, 1970; Tanaka et al., 1999). Seismic techniques are among the most effective geophysical methods for generating detailed subsurface geological models (Janjua et al., 2016; Iqbal et al., 2019; Eleslambouly et al., 2024). Additionally, temperature, porosity, and permeability values calculated from well logs can be overlaid on these models to visualize petrophysical variations in potential geothermal reservoirs (Nyakilla et al., 2025; Rashid et al., 2022). Geophysical techniques, including gamma-ray logging, nuclear magnetic resonance (NMR), and fluid temperature profiling, play an important role in understanding subsurface conditions. These methods are frequently integrated with approaches such as drilling, borehole geophysical surveys, direct subsurface sampling, and petrophysical property measurements. Borehole logging and inter-bore tomography are highly effective for delineating fracture networks and identifying fluid migration pathways (Bense et al., 2013). Additionally, borehole investigations often employ temperature, electrical conductivity, and flow velocity logs to identify hydraulically active fractures within bedrock aquifers under natural conditions (Agbotui et al., 2025).

Pakistan is a developing country that faces energy shortage, so efforts are required to develop new economic energy resources to meet its growing demand (Mir et al., 2023). The tectonic framework of Pakistan, which includes the traverse of the global seismic belt and widespread development of alteration zones, suggests that it has a considerable number of exploitable geothermal energy sources, as Quaternary volcanism, fumaroles, and hot springs are present there (Gondal et al., 2017; Zaigham et al., 2009; Ahmad et al., 2007). These features are classified into geopressured systems, seismo-tectonic systems, and Neogene–Quaternary volcanism systems (Wahid et al., 2022; Zaigham et al., 2009). The literature indicates that geothermal resources in Pakistan have only been identified to a limited extent, and their classification remains general and lacks detailed subsurface characterization. The narrow Sulaiman and Kirthar fold and thrust belts, which extends along the Indus River, is indicated as a potential zone for geothermal exploration (Ahmad and Rashid, 2010). In the Sulaiman Fold and Thrust Belt (SFTB), active tectonics, abnormal thermal gradients (3–4 °C/100 m), and several hot springs classify the area as a geopressured system. This system generally belongs to a medium enthalpy framework, where temperatures range from approximately 90 °C to 150 °C. These resources are primarily utilized for direct heating purposes and can also generate electricity through binary cycle technologies such as the organic Rankine cycle. In contrast, high-enthalpy systems exceed 150 °C and are typically found in volcanic regions or active tectonic zones. Such systems are mainly used for power generation using flash steam or dry steam plants (Wishart, 2020). Primarily, the Zindapir, Jaffarabad, and Gulki hot springs are present in the SFTB along the Zindapir Anticlinorium Zone. High levels of seismic activity with magnitudes ≥Mw 6 at shallow depths are also occasionally reported in the region (Figure 1c) (Reynolds et al., 2015). However, no thorough studies have been conducted to evaluate the economic viability of geothermal energy in the SFTB. The objective of this research is to evaluate the potential of geothermal reservoirs along the Zindapir Anticlinorium Zone in the SFBT. To achieve this objective: (i) an AHP prediction model is generated for delineating suitable geothermal potential areas; (ii) the hydrogeochemical properties of hot spring and stream water are examined to understand the impact of water–rock interactions and the origin of hot springs; (iii) crustal variations and heat flow scenarios are investigated to evaluate geothermal potential zones in the subsurface; (iv) the lateral continuity and petrophysical modeling of hydrocarbon reservoir formations are analyzed to identify the potential of geothermal resources in these formations. The general aim of this research is to integrate spatial, geochemical, and geophysical data to determine geothermal prospect zones for exploration. These methods help reduce exploration risk and support the sustainable utilization of geothermal energy which can benefit socio-economic development in the region.

Figure 1

Figure 1. (a) Geological map of the Central Indus Basin. (b) Map of the Zindapir Anticlinorium Zone showing the locations of 2D seismic lines, oil wells, gravity data points, magnetic data points, hot springs, and geological and geographical features (Raza et al., 1989; Karaca et al., 2021; Nazeer et al., 2013). (c) Map showing earthquake distribution with different depth ranges in the study area.

2 Study area and geological setting

This research work was carried out along the Zindapir Anticlinorium Zone situated in the Sulaiman Fold and Thrust Belt (SFTB), which is part of the Central Indus Basin in northwestern Pakistan (Figure 1a). The Central Indus Basin consists of the Punjab Platform to the east, the Sulaiman foredeep in the center, and the SFTB to the west (Kadri, 1995). The area features various tectonic features due to the collision of the Indo-Pakistan and Eurasian Plates during the post-Cretaceous orogenic event (Javed et al., 2021; Mateen et al., 2022; Powell, 1979; Jadoon et al., 2024). The SFTB is a major geological feature located near the western edge of the Himalayas in Pakistan. It has a lobate geometry and spans approximately 400 km in both width and length (Karaca et al., 2021). It mostly consists of open and upright folds that run parallel to subparallel. The Kingri, Domanda, Sulaiman, Chaudawan, and Mughal Kot faults are high-angle thrust and strike-slip faults trending in a N–S direction (Kazmi and January 1997). Normal faults occur locally because of flexural slip and gravitational collapse along steep fold limbs. They are typically secondary features that accommodate localized extension within the predominantly compressional tectonic regime (Kazmi and January 1997). The length of the SFTB mountain front is approximately 1,000 km, bounded by a monoclinal structure and folds such as Sui in the south and Zindapir in the east. It is divided into two zones: the Sulaiman inner and outer folded zones (Reynolds et al., 2015). The Zindapir Anticlinorium Zone is present in the Sulaiman inner folded zone, also known as the “eastern margin” of the SFTB. It is a gentle eastward-leaning fold with a convex axis across the nearly vertical, west-dipping axial surface. The Zindapir Anticlinorium Zone is bounded by the Sulaiman foredeep and the Sulaiman outer folded zone to the east and west, respectively, whereas the Pezu uplift is situated to the north and the Mari Bugti Zone is to the south (Qadri et al., 2022; Solangi et al., 2014). The Zindapir Anticlinorium Zone spans approximately 5,600 km², comprising a north–south-oriented fold system that includes the Dhodak, Rodho, Affiband, and Zindapir anticlines (Nazeer et al., 2013). Eocene rocks are exposed in the cores of the Dhodak, Rodho, and Affiband anticlines, while the Zindapir anticline consists of Paleocene and Upper Cretaceous strata in its core (Nazeer et al., 2013; Karaca et al., 2021). Pliocene, Miocene, and Neogene molasse deposits are founded along the limbs of these anticlines (Figure 1b). It is a promising area for hydrocarbon exploration due to its multiple gas condensate and oil reservoirs (Ahmed et al., 2023; Wajid et al., 2021). Previously, various multidisciplinary studies have been carried out in the SFTB which were primarily focused on its tectonic and sedimentary evolution for hydrocarbon exploration (Jadoon et al., 2024; Nazeer et al., 2018; Zaidi et al., 2021). Numerous studies have also been conducted on geological aspects, deformation mechanics, earthquake analysis, seismic reflection interpretation, and gravity modeling in the region. According to these, the region is influenced by thin-skinned tectonics with preferential convergence along a weak decollement (Sarwar, 1979; Quittmeyer et al., 1984; Banks and Warburton, 1986; Jadoon et al., 1993; Davis and Lillie, 1994).

The stratigraphy of the SFTB is divided into three main groups: molasse deposits from the late Oligocene to Holocene, shallow to deep marine rock units from the Eocene to Permian, and Khojak flysch deposits from the Late Eocene to early Oligocene (Saif-Ur-Rehman et al., 2019). The thickness of these units varies 5,000–10,000 m, with an average of 7,000 m in Mesozoic and early Tertiary formations (Ali et al., 2022; Raza et al., 1989). Paleozoic sedimentary rocks overlie in the basement in the SFTB (Figure 2). However, it is assumed that the Infra-Cambrian to Permian rock sequence is also present, based on wells drilled in the Punjab Platform (Raza et al., 1989; Humayon et al., 1991; Jadoon et al., 1994). The Jurassic and Cretaceous successions were deposited in a pre-collisional condition in the region. The Jurassic strata primarily comprise limestone and shale, while the dominant Cretaceous formations include shelf carbonates deposits of the Parh Formation and near-shore clastic deposits of the Mughal Kot Formation. Additionally, channelized sand deposits are present in the upper Cretaceous Pab sandstone (Khan and Scarselli, 2021). The Zindapir Anticlinorium has proven reservoirs in the Chiltan limestone (Jurassic period), Lower Goru Formation (Early Cretaceous), and Pab sandstone (Late Cretaceous). Additionally, hydrocarbon reservoirs are also found in the Ranikot and Dunghan Formations from the Paleocene epoch (Nazeer et al., 2013; Asim et al., 2014).

Figure 2

Figure 2. Generalized stratigraphic column of the Sulaiman Fold and Thrust Belt with source, seal, and proven reservoir formations (Asim et al., 2014; Nazeer et al., 2013; Saif-Ur-Rehman et al., 2019).

3 Materials and methods

The methodology of this research was based on integrated approaches that included GIS and remote sensing, geochemical analysis, and geophysical studies to identify the potential of geothermal resources in the study area. The methodology to achieve the research objectives is outlined in Figure 3.

Figure 3

Figure 3. Flowchart showing the methodology adopted for identifying geothermal potential zones in the SFTB.

3.1 GIS and remote sensing

In this analysis, different datasets were incorporated using the ArcGIS 10.8 platform. This included infrared remote sensing data, earthquake data, geological data (lithology and faults), and geomorphological data (elevation and stream order), which were considered major evidence that might influence geothermal activities. For calculating land surface temperature (LST), Landsat 8 imagery with 30 m resolution was utilized. This imagery consists of 11 bands in which Band 10 was used for thermal anomaly detection, and Bands 4 and 5 were used for calculating the Normalized Difference Vegetation Index (NDVI). The raster calculator in ArcMap software was utilized to calculate the top-of-atmosphere spectral radiance from the thermal infrared digital numbers (Band 10) of the imagery using the rescaling factor (Equation 1). The radiance at the top of the atmosphere was then used for calculating satellite brightness temperature (Equation 2). Furthermore, the classification of the NDVI was calculated using Equation 3. In this study, the error caused by NDVI in surface temperature was found to be less than 0.1 K after radiometric calibration. Hence, the temperature ratio and specific surface emissivity of land surfaces are strongly linked to the amount of vegetation covering the ground. Estimation of vegetation coverage was determined using Equation 4. Land surface emissivity was calculated using Equation 5. To calculate LST, the results of Equation 2 and Equation 5 were used in Equation 6.

${\mathrm{T}\mathrm{O}\mathrm{A}}_{\mathrm{r}\mathrm{a}\mathrm{d}}={\mathrm{R}}_{\mathrm{T}\mathrm{I}\mathrm{R}\mathrm{S}10}\times \mathrm{D}\mathrm{N}+{\mathrm{B}}_{\mathrm{r}\mathrm{a}\mathrm{d}}-{\mathrm{C}}_{\mathrm{T}\mathrm{I}\mathrm{R}\mathrm{S}10},(1)$

${\mathrm{B}\mathrm{T}}_{\left(\mathrm{T}\mathrm{O}\mathrm{A}\right)}=\raisebox{1ex}{$\mathrm{K}2$}\!\left/ \!\raisebox{-1ex}{$\ln $}\right.\left(\raisebox{1ex}{$\mathrm{K}1$}\!\left/ \!\raisebox{-1ex}{${\mathrm{T}\mathrm{O}\mathrm{A}}_{\mathrm{r}\mathrm{a}\mathrm{d}}$}\right.+1\right)-273.15,(2)$

$\text{NDVI}=\left(\text{Band\hspace{0.17em}}5-\text{Band\hspace{0.17em}}4\right)/\left(\text{Band}+\text{Band\hspace{0.17em}}4\right),(3)$

${\mathrm{P}}_{\mathrm{V}}={\left(\text{NDVI}-{\text{NDVI}}_{\min }/{\text{NDVI}}_{\max }-{\text{NDVI}}_{\min }\right)}^2,(4)$

$\mathrm{\epsilon }=0.004\times {\mathrm{P}}_{\mathrm{V}}+0.986,(5)$

$\text{LST}={\text{BT}}_{\left(\text{TOA}\right)}/\left(1+\mathrm{\lambda }\times {\text{BT}}_{\left(\text{TOA}\right)}/\mathrm{c}2\times \ln \left(\mathrm{\epsilon }\right)\right).(6)$

Here, ${\text{TOA}}_{\text{rad}}$ is the radiance of top of atmosphere. ${\mathrm{R}}_{\text{TIRS}10}$ is the multiplicative radiance of the thermal infrared sensor. DN is the digital number consisting of constant values for Band 10. ${\mathrm{B}}_{\text{rad}}$ is the radiance add. ${\mathrm{C}}_{\text{TIRS}10}$ is the correction value for the band, which is −2.1 k. ${\text{BT}}_{\left(\text{TOA}\right)}$ is the satellite brightness temperature of the top of the atmosphere in degrees Celsius. k1 and k2 are the constant values of Band 10. ${\mathrm{P}}_{\mathrm{V}}$ is the proportion of vegetation cover. ${\text{NDVI}}_{\min }$ and ${\text{NDVI}}_{\max }$ are the minimum and maximum values of the $\text{NDVI}$ layer. For Equation 4, the calculation of vegetation cover assumes that $\text{NDVI}>{\text{NDVI}}_{\max }$ and pixel value equals 1; if $\text{NDVI}<{\text{NDVI}}_{\min }$, the pixel value equals 0. $\mathrm{\epsilon }$ is the land surface emissivity. $\text{LST}$ is land surface temperature. $\mathrm{\lambda }$ is the wavelength of emitted radiance. $\mathrm{c}2$ is the constant value (14,388 µmk).

Earthquake data were acquired from the Pakistan Meteorological Department (PMD) and the United States Geological Survey (USGS) for identifying seismic activity in the region. Data spanning 1980–2024 were collated, including details of earthquake epicenter coordinates, magnitude, and depth. For this study, the magnitude (M) of shallow-source earthquakes (<35 km) was used, assuming that they have a greater impact on surface features. The Gutenberg–Richter relationship was utilized to estimate the b-values of 43 seismic events using Equation 7. For this purpose, the maximum likelihood estimation (MLE) method was used. First, the lowest magnitude earthquake was acquired from the dataset, which indicates that all events are reliably recorded above this. This helps to calculate the mean magnitude through filtering and the dataset of earthquakes that were equal to or greater than the lowest magnitudes.

$\mathrm{Log}\left({\mathrm{\lambda }}_{\mathrm{m}}\right)=\mathrm{a}-\mathrm{b}\times \mathrm{M},(7)$

Where ${\mathrm{\lambda }}_{\mathrm{m}}$ is the number of seismic events and M is their magnitude, “a” is the distribution density of seismic events in magnitude, and b is the relationship between the relative frequencies and magnitudes of small versus large earthquakes. Inverse distance weighting (IDW) spatial interpolation was used to generate the b-value map in the study area.

Geological maps (1:500,000 scale) consisting of lithological distributions and faults were obtained from the Geological Survey of Pakistan (GSP), which were further digitized and converted to raster format. The lithology in the SFTB region was divided into seven specific categories, in which varying weight values were assigned to each category depending on the specific rock units consisting of hot springs and their probable individual influence on geothermal resources. A digital elevation model (DEM) of 12.5 m resolution was extracted from the Advanced Land Observing Satellite Phased Array type L-band Synthetic Aperture Radar (ALOS PalSAR) to extract the elevation and stream order of the study area (Wahid et al., 2016a; Gentana et al., 2019). The surface runoff and overland flow model was generated in ArcGIS to establish the flow direction. For this purpose, filling depressions, flow directions, flow accumulation, and stream order were analyzed. Overall, water system distribution analysis was conducted to relate their connections with hot springs in the research area. Furthermore, multi-ring buffer analysis was performed on faults and streams by establishing different buffer distances confined to the boundary of the study area.

Lastly, all thematic layers were integrated spatially as evidence weight factors to generate a geothermal energy potential map using the Analytical Hierarchy Process (AHP). For this purpose, the values were adjusted based on different parameters in the input raster to align with a consistent assessment level. Hence, land surface temperature, seismic activity, lithological distributions, distance to faults, distance to streams, and elevation were considered as the input rasters for this study. Each input raster pixel value was multiplied by its corresponding importance weight. Furthermore, the input weights of the rasters were combined to generate the result using Equation 8.

$\text{Output}={\mathrm{w}}_1\times {\text{Input}}_1+{\mathrm{w}}_2\times {\text{Input}}_2+\dots {\mathrm{w}}_{\mathrm{n}}\times {\text{Input}}_{\mathrm{n}}.(8)$

The final comprehensive result is the output, with ${\mathrm{w}}_1$, ${\mathrm{w}}_2$, $\dots$ ${\mathrm{w}}_{\mathrm{n}}$ representing the corresponding weights, and ${\text{Input}}_1$, ${\text{Input}}_2$, $\dots$ ${\text{Input}}_{\mathrm{n}}$ denoting the input data. The accuracy of the prediction model was validated with hot springs points, providing a scientific basis for decision-making in geothermal resource development.

3.2 Geochemical analysis

This analysis was conducted on samples of hot springs and stream water in the study area. Water samples were acquired in polythene bottles from the hot springs at Zindapir, Gulki, and Jaffarabad and from streams running along them. The temperature, conductivity, salinity, and resistivity of the hot springs were also measured using the INESA REX DDBJ-350F portable conductivity meter. These samples were sent to the Pakistan Institute of Nuclear Science and Technology (PINSTECH), Islamabad, for laboratory analysis. Inductively coupled plasma-mass spectrometry (ICP-MS) analysis was conducted on representative samples to identify the ionic concentration. The analysis involved filling a 100 mL high-density polyethylene bottle (HDPE) and acidifying the sample with 1% v/v HNO₃. The water sample bottles had been thoroughly rinsed three times before filling. Origin software was then used to generate Schoeller, Piper plots, and Na–K–Mg ternary diagrams for identifying ionic concentration, hydrothermal regime, and water–rock equilibrium conditions, respectively.

3.3 Geophysical methods

3.3.1 Gravity method

The gravity survey in the field was conducted using the Scintrex CG-5 Autograv system, which has a resolution of 0.001 mGal. A total of 239 gravity measurements were taken at intervals of 2.5–3 km to interpret the crustal variations in the study area. Due to the limitation of absolute gravity measurement in the CG-5 gravity meter, the measured gravity values were connected to the gravity base station located at Rawalpindi (Pakistan), with known absolute gravity used for this research. This is a reference base station that is part of a worldwide network known as the International Gravity Standardization Net-1971 (IGSN-71). A local base station was established at the University of Azad Jammu and Kashmir based on the gravity difference between the local and reference base stations. This base station was used as an auxiliary base station in the study area to help in connecting the acquired gravity values with the international gravity network. Furthermore, standard corrections were applied to measure the complete Bouguer anomaly values at each station. These include drift, terrain, tidal effects, latitude, and elevation corrections which affect the gravity field readings. The gravity data were reduced after applying the corrections using Microsoft Excel. In this study, the minimum curvature technique was adopted, and a 1,500 m data interval was used to grid the data in Excel to generate Bouguer anomaly maps. Next, the regional and residual parts were separated from the Bouguer anomalies, which were indicated by longer and shorter effective wavelengths, respectively, using the upward continuation method. This filtering method was applied using Geosoft Oasis Montaj software. After isolation, the regional anomalies were subtracted from the Bouguer anomalies to obtain the residual anomalies using the grid math tool.

3.3.2 Magnetic method

A regional total magnetic intensity (TMI) map was generated using satellite-based data acquired from the World Magnetic Model (https://www.ncei.noaa.gov/) developed by the National Centers for Environmental Information (NCEI) and the British Geological Survey (BGS). The data were taken from the same stations from which the gravity data were acquired. The TMI grid was prepared in Geosoft Oasis Montaj software at a scale of 1:30,000 with a contour interval of 20 nT using the minimum curvature gridding technique. To correlate size and thermal characteristics, the research area was divided into ten blocks of 22 km × 22 km each to estimate the Curie point depth (CPD). To achieve this, 2D spectral analysis of each block was performed individually (Spector and Grant, 1970) to estimate the depth of the magnetized rectangular prisms (${\mathrm{Z}}_{\mathrm{t}}$) from the slope of the log power spectrum. The depth of the centroid of the magnetic source bodies (${\mathrm{Z}}_{\mathrm{o}}$) was also estimated. The bottom depth of the magnetic bodies (${\mathrm{Z}}_{\mathrm{b}}$) was calculated using spectral analysis. Furthermore, windows of adequate width were employed to identify the most profound magnetic source. Blocks with an overlapping length of 11 km were used, and their midpoint was utilized as the representation of ${\mathrm{Z}}_{\mathrm{b}}$. The plot of the radial average power density spectra was generated to determine ${\mathrm{Z}}_{\mathrm{t}}=\left(\ln \left[{\mathrm{P}\left(\mathrm{k}\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\right]\right)$ and ${\mathrm{Z}}_{\mathrm{o}}=\left.\left(\ln {\left[\mathrm{P}\left(\mathrm{k}\right)\right.}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}/\mathrm{k}\right]\right)$. The top of a magnetic source (${\mathrm{Z}}_{\mathrm{t}}$) was determined by fitting a straight line to the middle-high wavenumber part of the power density spectral curve, and the depth to the centroid (${\mathrm{Z}}_{\mathrm{o}}$) of a magnetic source was derived from the lower wave numbers of the power density spectra. Equation 9 was used to obtain an estimate for the basal depth:

${\mathrm{Z}}_{\mathrm{b}}=2{\mathrm{Z}}_{\mathrm{o}}-{\mathrm{Z}}_{\mathrm{t}},(9)$

Where (${\mathrm{Z}}_{\mathrm{b}}$) is the calculated depth of the magnetic source, assumed to be the CPD at which ferromagnetic minerals are transformed into paramagnetic minerals at a temperature of approximately 580 °C (853 K). To estimate the geothermal gradient (GG), Equation 10 was used:

$\text{GG}=580℃/\text{CPD}.(10)$

It is assumed that the top of the magnetic source is determined with the heat flow, as magnetic minerals lose their magnetization at the CPD. Furthermore, the Curie temperature of magnetite is approximately 580 °C. Lastly, heat flow (Q) was estimated using Equation 11:

$\mathrm{Q}=\mathrm{\lambda }\left(580℃/\text{CPD}\right),(11)$

where the thermal conductivity coefficient is denoted by λ. The average thermal conductivity of 2.5 W m^–1 C⁻¹ was used as the standard for igneous rocks.

3.3.3 Seismic and petrophysical investigation

For seismic and petrophysical investigation, 2D seismic lines, wireline logs (gamma ray, bulk density, neutron porosity, temperature, and resistivity), and well tops were utilized to interpret the lateral variations of potential geothermal reservoirs and their petrophysical characteristics, respectively. The 2D seismic lines and well data were thus acquired from the Directorate General of Petroleum Concessions (DGPC), Pakistan. The 2D seismic lines and well data were interpreted using Petrel (Version 2017; Schlumberger) and Techlog (Version 2015.3, Schlumberger) software. Different horizons and faults were interpreted on the 2D seismic section in the time domain. The reflectors on the seismic section were marked based on the coherence and continuity of the seismic amplitude. The interpreted time horizons were converted into the depth domain using the velocity functions. Equation 12 was used to convert time horizons to depth horizons.

$\text{Depth}={\mathrm{V}}_{\text{average}}\times \mathrm{t}/2.(12)$

Petrophysical analysis was performed on different wells. In this study, the results of Zindapir-01 are shown due to commercial confidentiality. However, this representative well was drilled to a significant depth of the proven hydrocarbon reservoirs in the study area. The wireline logs and well-top data were utilized to evaluate the temperature distributions, lithological identification, porosity, and permeability calculations in potential geothermal reservoirs. Equations 13–18 were used to estimate these parameters. The Horner plot was utilized to estimate the corrected temperature, considered the most dependable indicator of equilibrium temperature:

${\text{Temperature}}_{\text{Bottom}}={\text{Temperature}}_{\text{Surface}}+\left(\text{Depth}/\text{Geothermal\hspace{0.17em}gradient}\right),(13)$

$\text{Volume\hspace{0.17em}of\hspace{0.17em}shale\hspace{0.17em}}\left(\text{Vsh}\right)={\text{GR}}_{\log }-{\text{GR}}_{\min }/{\text{GR}}_{\max }-{\text{GR}}_{\min },(14)$

$\text{Density\hspace{0.17em}Porosity}={\text{density}}_{\text{matrix}}-{\text{density}}_{\log }/{\text{density}}_{\text{matrix}}-{\text{density}}_{\text{fluid}},(15)$

$\text{Neutron\hspace{0.17em}Porosity}={\text{Neutron}}_{\log },(16)$

$\text{Average\hspace{0.17em}Porosity}=\left(\text{\hspace{0.17em}Density\hspace{0.17em}Porosity}+\text{Neutron\hspace{0.17em}Porosity}\right)/2,(17)$

$\text{Permeability}=100\times \left({\left(\text{Porosity}\right)}^2\times \left(1-{\text{Sw}}_{\text{irr}}\right)/{\text{Sw}}_{\text{irr}}\right),(18)$

Where ${\text{GR}}_{\log }$ is the gamma ray log value and ${\text{GR}}_{\max }$ and ${\text{GR}}_{\min }$ represent the maximum and minimum values in the gamma ray log used for lithological variations integrated with well tops. Irreducible water saturation is denoted by ${\text{Sw}}_{\text{irr}}$.

Lastly, a 3D grid was generated by incorporating the petrophysical analysis on depth models of selected reservoir formations. For this purpose, the results of temperature distribution, porosity, and permeability calculations from different wells were upscaled (Iqbal et al., 2017; Salim et al., 2015). The average well log values in each cell were used based on arithmetic averaging for continuous logs and the “most of” method for discrete logs. Subsequently, preprocessing was utilized to generate the lateral variations of these parameters using the Kriging method—considered the best unbiased geostatistical estimation technique (Wahid et al., 2016b; Falivene et al., 2007). Furthermore, the CEI was calculated to evaluate favorable locations for geothermal energy extraction in the selected reservoirs using Equation 19:

$\text{CEI}={\left({\text{Temperature}}_{\text{Reservoir}}\times {\text{Porosity}}_{\text{Reservoir}}\times {\text{Permeability}}_{\text{Reservoir}}\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$3$}\right.}.(19)$

4 Results and discussion

4.1 Geothermal distribution assessment

Land surface temperature, seismic history, lithological variations, fault distributions, elevation, and stream data were utilized for geothermal distribution assessment. The LST map illustrates that land surface temperature range 25–50 °C in the study area (Figure 4a). The Zindapir Anticlinorium Zone has higher surface temperature anomalies, whereas temperature decreases toward the east. Moreover, the hot springs are situated in high temperature regions, suggesting a direct relationship between geothermal anomalies and LST. The Gutenberg–Richter b-value map indicates that the study area is seismically active, with frequent higher and lower magnitude earthquakes (Figure 4b). The map shows that the hot springs are situated in areas with medium to high b-value ranges. Higher b-values indicate a higher proportion of smaller magnitude earthquakes due to lower stress. Lower b-values indicate areas of higher stress which result in a significant number of larger magnitude earthquakes. The study area is an active seismic region where small magnitude earthquakes occur frequently. It is assumed that these seismic activities influence crustal fractures and water movement, creating conditions favorable for geothermal energy development. Moreover, the movement and retention of energy are also impacted by lithological distributions and faults. The research area consists of sedimentary rocks, in which Paleocene and Eocene rocks are exposed in the core and younger sedimentary rocks are deposited along the limbs of the Zindapir Anticlinorium Zone (Figure 4c). The hot springs of Zindapir and Jaffarabad are situated in Paleocene rocks, whereas the Gulki hot spring is present in Eocene rocks. This zone also consists of dense local faults bounded by the regional Sulaiman Range Fault and Domanda Fault to the east and west, respectively (Figure 1b). The buffer zones of 300 m, 500 m, and 800 m were applied to these faults to assess their impact on geothermal activity (Figure 4d). The map indicates that the hot springs are situated within the 300 m buffer zones. Based on this, it is inferred that there is a direct relationship between faults and geothermal activities in the study area.

Figure 4

Figure 4. Major evidence factors: (a) land surface temperature (LST); (b) Gutenberg–Richter b-value map; (c) lithostratigraphic units exposed at the surface; (d) distance to faults; (e) digital elevation model; (f) distance to streams utilized for geothermal distribution assessment. (g) The predicted model shows different ranges from very high to very low geothermal potential zones in the study area.

The grading of water systems was established based on morphology and stream flow. It consists of dense streams, which were divided into primary and secondary streams. The primary streams are higher grade streams (stream order 4–5), and the secondary streams are lower grade streams (stream order 1–3). The secondary streams mostly flow from higher elevations (west) and recharge the primary streams toward the east. The primary streams also pass along the hot springs in the study area (Figure 1b). Different buffer zones of 300 m, 500 m, and 800 m were assigned to stream channels. The analysis indicates that all the hot springs are in lower elevated areas and present within the 300 m stream buffer zone (Figures 4e,f). Based on this, it is assumed that these streams might contribute to recharging the hot springs.

Finally, the geothermal potential map was generated by integrating these input layers (thematic maps) using AHP analysis. This is based on the weightage assigned to these layers, which might impact the occurrence of geothermal activities. The research area was divided into five classes, categorized from very high to very low geothermal potential. The very high to high geothermal potential classes are indicated by red and orange, respectively, whereas green, blue, and dark blue colors similarly indicate intermediate, low, and very low geothermal potential classes (Figure 4g). The geothermal potential map indicates that the Zindapir Anticlinorium Zone has very high to high geothermal potential classes compared to adjoining areas. Eocene and Paleocene clastic rocks and carbonates are thrust in this area, indicating that, in this zone, the porosity in clastic rocks and dissolution of carbonate rocks due to fluid movement enhance the surface area and permeability, which contributes to heat exchange and might result in hot springs. The surface temperature and the Gutenberg–Richter b-value map also suggest that the Zindapir Anticlinorium Zone and its surroundings have the highest surface temperature ranges and that this is a seismically active region. Regional and local faults and primary water systems are significantly developed throughout this zone. It is assumed that regional faults provide the heat from the subsurface, which further reaches the surface through pathways provided by local faults. The elevation and stream order maps indicate that water flows from higher to lower elevations, whereas high-grade streams pass through the hot spring areas. This might contribute to recharging the hot spring water. The eastern portion of the Zindapir Anticlinorium Zone has areas of very low to low geothermal potential values. It is observed that all hot springs projected within the model are situated in the very high geothermal potential class, supporting the model’s precision and reliability.

4.2 Groundwater origin and temperature evaluation

Water samples were acquired from the hot springs of Gulki, Jaffarabad, and Zindapir, along with their associated streams, to identify their origin and reservoir temperature; these hot springs have average temperature ranges of 58 °C, 40 °C, and 50.3 °C, respectively. Based on the geochemical analysis, the pH values of the water range from 7.45 to 8.17, indicating hyperthermal slightly alkaline water (Table 1). These variations in value might be due to changes in temperature and the flow paths of the water. The Gulki hot spring (GHS) is situated at a higher elevation (299 m) with a TDS range of 2,370 mg/L, whereas its associated stream water (GSW) has TDS ranges of 1,030 mg/L. The Zindapir hot spring (ZHS) is situated at an elevation of 270 m with a TDS range of 4,065 mg/L; for its associated stream water (ZSW), the TDS is 1,319 mg/L. Moreover, the TDS range of 3,655 to 3,500 mg/L is present in the Jaffarabad hot springs (JHS-1 and JHS-2). It is inferred that the water recharge of the Zindapir and Jaffarabad areas is different from the Gulki area. The higher TDS values in ZHS compared to JHS-1 and -2 might be due to higher temperature ranges. In hot springs, electric conductivity (EC) ranges between 1,750 and 6,570 μS/cm, with an average of 4,912.5 μS/cm. The higher pH, TDS, and EC values in the hot springs indicate the saline nature of the water. This higher salinity is due to the dissolution of carbonate rocks because of water circulation through fractures and faults. On the other hand, sodium (Na) and bicarbonate (HCO₃) concentrations in ZHS, JHS-1, and JHS-2 are higher (>1,000 mg/L), reaching 1,745 mg/L and 2,366 mg/L, respectively, whereas GHS has low concentrations of Na and HCO₃. The high Na concentrations in the thermal waters are due to hydrolysis reactions (chemical weathering) of feldspar (alkali feldspar) in sediments, derived from the Paleozoic basement and basaltic magma. These higher concentrations indicate an increase in groundwater residence time and enhanced rock–water interaction. The high bicarbonate concentrations are due to reactions between circulating meteoric waters and limestone, resulting in CO₂ rich waters. However, the high magnesium content in GHS indicates greater cold-water mixing than other springs. The lower Ca concentration in ZHS could be attributed to the induced precipitation of carbonate by intense evaporation. Lastly, the higher SiO₂ content is explained by the weathering of siliceous rocks in the hot springs, causing the transfer of silica from rocks to water in the study area. In ZHS, the results of SO₄ concentration with higher Cl concentrations suggest the mixing of fossil seawater with meteoric water.

Table 1

Table 1. Representative geochemical analysis results of hot springs and their associated stream water samples.

The Schoeller plot indicates that the concentrations of potassium (K), magnesium (Mg), and calcium (Ca) are lower than Na, HCO₃, and Cl in the study area (Figure 5a). To identify the hydrothermal regime, type of water, and rock–water relations, a Piper plot was generated. In Figure 5b, a diamond plot is centrally located and merged with two tri-linear plots on either side. In this plot, water with high calcium sulfate (gypsum ground and mineral seepage) lies in the top quadrant, calcium bicarbonate water (shallow crisp groundwater) in the middle quadrant, sodium chloride water (marine and deep old groundwater) in the right quadrant, and sodium bicarbonate (profound groundwater influenced by particle exchange) at the base quadrant. The diamond diagram illustrates that most samples are a sodium chloride type, whereas ZHS lies in the sodium bicarbonate type. Overall, the samples indicate a composition mainly of Na + K, Ca type, with Na + K, SO₄ as the dominant ions. Based on the Piper plot, it is inferred that the recharge of hot springs is mostly from clastic origins compared to ZHS, which has an origin from carbonate rocks in the subsurface. The correlation plot (Na⁺ + K⁺ vs. Ca⁺² + Mg⁺²) indicates that total ionic salinity (TIS) ranges 20–78 meq/kg in the springs of the study area (Figure 5c). However, the Na⁺ concentrations are associated with the runoff time of hot spring water, which suggests its deep circulation with an impact on thermal activities. Therefore, ZHS, JHS-1, and JHS-2 have deeper circulations than GHS. Furthermore, the high concentration of Na⁺ is due to dissolution during the water–rock process. As a result, the feldspar in the formations reacts with the geothermal water and generates Na⁺ concentrations in the hot springs of the study area. To determine water–rock equilibrium conditions and subsurface temperatures in the studied hot spring water samples, the Na–K–Mg ternary diagram (Giggenbach, 1988) was utilized; this diagram is useful in regions of active tectonics and continental geothermal areas. It specified that the GHS temperature ranges 180–200 °C and lies in the zone of immature water, whereas ZHS, JHS-1, and JHS-2 have temperature ranges of 90–100 °C in the subsurface, indicating a partially equilibrated water zone (Figure 5d).

Figure 5

Figure 5. (a) Schoeller plot showing ion concentrations. (b) Piper plot indicating the origin of water, type of water, and rock–water interaction. (c) Na⁺ +K⁺ vs. Ca⁺² + Mg⁺² correlation plot, indicating total ionic salinity. (d) Na–K–Mg ternary plot, specifying groundwater origin and temperature estimation.

4.3 Crustal variations and heat flow scenarios

The gravity method was used to understand the crustal variations and thickness of sedimentary rocks in the subsurface. These variations result from tectonic movements that release energy, potentially contributing to geothermal energy generation. To address this, a Bouguer anomaly map was generated using field gravity measurements. This map represents the combined effects of regional and residual anomalies based on density contrast (Figure 6a). Overall, the gravity gradient varies from −127.4 to −225.1 mGal, in which contour closures indicate different anticlines. Two contour closures with comparatively high values (ranges from −135 to −147 mGal) in the north indicate the Dhodak and Rodho anticlines; in the south, the contour closures suggest the Affiband and Zindapir anticlines. The decrease in gravity values in the east and west indicates the Sulaiman Depression and the limb of the Bharti syncline, respectively. The highest Bouguer anomalies, ranging from −127 to −142 mGal in the south and southeast, indicate crustal thinning in the study area which may be due to structural deformation or strike-slip movement in the south. Furthermore, regional anomalies were separated using the upward continuation technique. This separation was based on contour trends generated at different heights from the Bouguer anomaly grid. The general contour trend on the regional Bouguer anomaly is toward the SSW–NNE direction. The gravity variations range from −189.8 to −128.3 mGal, as shown in the regional Bouguer anomaly map (Figure 6b). The lowest gravity anomalies in the northwest indicate higher thicknesses of sedimentary rocks. The higher gravity anomalies in the southeast might be associated with basement uplifting.

Figure 6

Figure 6. (a) Complete Bouguer anomaly map indicating the gravity gradient due to crustal variation. (b) Regional Bouguer anomaly map. (c) Residual Bouguer anomaly map.

The residual Bouguer anomaly map represents gravity variations from −38.8 to 24.5 mGal (Figure 6c). The major faults (Sulaiman Range Fault and Domanda Fault) could be responsible for changes in the gravity gradient due to crustal variations and may delineate geological boundaries. The fault trends are generally aligned with the density contours along the Zindapir Anticlinorium Zone, indicating that these are thin-skinned faults. Based on this, it is assumed that thick-skinned faults associated with the basement are not present in the study area. Furthermore, the core of the Zindapir Anticlinorium Zone has the highest density values, indicating anticlinal pop-up structures. The limbs of this zone toward the east and west indicate lower density values, which could correlate with thick sedimentary alluvium deposits toward the limbs. This suggests that synclinal structures have comparatively higher sedimentary thickness than the Zindapir Anticlinorium Zones in the study area. Overall, this zone is associated with faults and uplifting that bring denser material toward the surface.

The Bouguer anomaly data were utilized to estimate sedimentary thicknesses ranging 3.42–12.48 km (Table 2). The maximum and minimum sedimentary thicknesses along the Dhodak anticline reach up to 7.09 km and 4.32 km along its western limb and core, respectively. The Rodho anticline is situated southeast of the Dhodak anticline. The sedimentary thicknesses along this anticline range between 12.48 km and 5.34 km toward its western limb and core, respectively. The sedimentary thickness along the Affiband anticline ranges between 9.55 km and 4.69 km toward its western limb and core, respectively. The Zindapir anticline has the shallowest sedimentary thickness ranges between 4.42 km and 3.42 km toward its limbs and core, respectively. ZHS, JHS-1, and JHS-2 are also located along this anticline. The results indicate that tectonic activities along the Zindapir Anticlinorium Zone caused the uplifting of strata, with the southern part more influenced by tectonism.

Table 2

Table 2. Estimated sedimentary thicknesses and Moho depth along the Dhodak, Rodho, Affiband, and Zindapir anticlines.

Magnetic data were also utilized to perform heat flow measurements through spectral analysis. Magnetic anomalies are based on variations of magnetic minerals, resulting in local geothermal gradients in the subsurface. Therefore, a total magnetic intensity map was generated in the study area. This indicated that magnetic anomalies range from 48,956.4 to 49,182 nT, with a difference of 225.6 nT in magnetic susceptibility (Figure 7). The general contour trend of the TMI map is in the NW–SE direction. In the northern portion, magnetic susceptibility is high, ranging from 49,143.7 to 49,182 nT due to an increase in basement depth. In the southern portion, magnetic susceptibility is low, ranging from 48,956.4 to 49,079.2 nT. This indicates that the basement is shallower than the northern zone, as the results reveal that magnetic minerals are near the surface.

Figure 7

Figure 7. Total magnetic intensity (TMI) map specifies the regional magnetic anomalies caused by susceptible minerals.

The magnetic data were transformed into the frequency domain to perform 2D radially averaged power spectral analysis to determine the top and centroid depths of the magnetic sources. These are average interface depths which were divided into suitable geographical window-size blocks with 50% overlapping. For magnetic source identification, the radially averaged power spectrum and scaled power spectrum were calculated for each block. This spectral analysis is based on susceptibility variations of the crust in the subsurface and helps compute the geothermal sources in the study area. The power spectrum curves for blocks 1 (Figures 8a,b) and 10 (Figures 8c,d) indicate the estimation of top (Zt) and centroid depth (Zo). Based on the calculations, Zt ranges 12.62–15.79 km and Zo 18.76 km–21.88 km (Table 3).

Figure 8

Figure 8. (a) Calculated top depth (Zt) to magnetic source of block 1. (b) Centroid depth (Zo) of magnetic source calculated for block 1. (c) Top depth (Zt) of magnetic source calculated for block 10. (d) Calculated centroid depth (Zo) of the magnetic source of block 10, based on two-dimensional radially averaged power spectral analysis.

Table 3

Table 3. Estimated Curie point depth (CPD), geothermal gradient (GG), and heat flow (Q) values.

The Curie point depth (CPD) estimated from Zo and Zt ranges 24.90–27.97 km. The calculated geothermal gradient and heat flow values range between 23.28 and 20.73 °C/km and 66.37 to 59.08 mW/m², respectively. The calculated CPD, GG, and Q values were plotted at the center of each block, which were re-gridded to generate contoured image maps. The CPD map shows that the Zindapir area, in which two hot springs are present, has a lower CPD than other parts of the study area (Figure 9a). The geothermal gradient and heat flow are also higher in this area (Figures 9b,c). Overall, the CPD, geothermal gradient, and heat flow maps suggest that the study area has potential geothermal reservoirs in the subsurface.

Figure 9

Figure 9. (a) Curie point depth (CPD), (b) geothermal gradient (GG), and (c) heat flow (Q) variation maps in the study area.

4.4 Geothermal reservoir characterization

Seismic interpretation and petrophysical analysis were conducted to evaluate the lateral continuity and geothermal reservoir potential in different formations. The seismic interpretation of the 2D strike line (O-785-SK-04) (Figure 10a) suggests that Late Eocene to Pliocene formations are exposed mainly along the Dhodak anticline toward the north, and the Ghazij Formation is exposed on the surface at the limb of the Zindapir anticline toward the south (Figure 10b). The complete sequence of Triassic to Paleocene formations is present in the subsurface, in which the Zindapir anticline has been mainly uplifted. The parallel reflectors suggest that the formations were deposited in a stable environment. The seismic amplitude anomalies indicate a gentle dip toward the north and folded structural geometry toward the south. The reverse/thrust faults indicate that the area is deformed and uplifted due to compressional tectonics. Seismic interpretation also suggests an increasing thickness of sedimentary formations from east to west. The En-echelon folds caused by shear stresses are responsible for the development of anticlines along the Zindapir Anticlinorium Zone. These stresses are the result of compressional forces from left-lateral transgressive movement in the study area. Overall, the study area is situated at the collision zone of the Indian and Eurasian plates, with the Himalayan fault and thrust system to the north. This collision is mainly responsible for the tectonic activities that resulted in the development of Zindapir Anticlinorium Zone.

Figure 10

Figure 10. (a) Uninterpreted 2D strike line (O-785-SK-04) and (b) the interpreted seismic line consisting of different geological formations deposited in the study area.

In this research, proven hydrocarbon reservoirs are considered potential geothermal zones in the study area. To address this, time and depth contour maps of the Pab Sandstone, Lower Goru Formation, and Chiltan Limestone were generated. Overall, the horizons deepen toward the northwest. In the north, the contour closure indicates an anticlinal feature where the Dhodak gas field is present. The horizons shallow toward the south, indicating the uplifting of sedimentary sequences due to the impact of compressional tectonic stresses in the study area (Figures 11a–c). In the northwest, the depth contour map of the Pab Sandstone and Lower Goru Formation shows maximum depths of 4,000–5,000 m and 4,800–5,500 m, respectively. Similarly, the Chiltan Limestone has a maximum depth range of 5,200–6,000 m in the study area. In the south, the depths of these formations range 400–900 m, 1,430–2,490 m, and 1,900–2,600 m, respectively (Figures 11d–f).

Figure 11

Figure 11. Time contour maps of (a) Pab Sandstone, (b) Lower Goru Formation, and (c) Chiltan Limestone. The depth contour maps of (d) Pab Sandstone, (e) Lower Goru Formation, and (f) Chiltan Limestone.

The depths of these horizons were validated with available well tops in the study area. Moreover, the crustal variations generated from gravity surveys also indicate similar trends in the seismic horizons. To identify geothermal potential in the selected reservoir formations, the temperature gradient, porosity, and permeability were estimated (Figures 12a–d). The temperature gradient shows a linear increase in temperature with respect to depth. In the Zindapir-01 well, the temperature ranges 30–155 °C from the surface to the total drilled depth of 4,500 m. The average porosity is 10.9% in the Pab sandstone, 6.51% in the Lower Goru Formation, and 2.09% in the Chiltan limestone (Table 4). Porosity increases in the lower part of the Chiltan limestone, indicating that it might be secondary porosity due to the impact of tectonic activity in the region (Figure 12d).

Figure 12

Figure 12. (a) Log responses to total depth of Zindapir-01 well. Temperature gradient, lithological variation, porosity, and permeability of (b) Pab Sandstone, (c) Lower Goru Formation, and (d) Chiltan Limestone.

Table 4

Table 4. Petrophysical results of the Zindapir-01 well, including temperature, porosity, and permeability ranges.

The permeability calculations show favorable permeability ranges in the selected reservoirs (Table 4). Overall, based on the temperature, porosity and permeability calculations, these formations can act as economic geothermal reservoirs in the study area. Temperature variations, porosity, and permeability distributions were generated using petrophysical analysis for the property modeling maps. These distributions were overlain on the horizons of potential geothermal reservoirs to indicate their lateral distribution, which help identify economic prospective zones in the study area. In the Pab sandstone, the temperature distribution model indicates a range of 40–120 °C in the study area (Figure 13a). In this formation, porosity and permeability range 6–16% and 0.3–936 mD, respectively (Figures 13d,g). The temperature distribution model of the Lower Goru Formation indicates a range of 120–220 °C (Figure 13b), with porosity ranging 2–11% and permeability 1.1–9,521 mD there (Figures 13e,h). In the Chiltan limestone the temperature distribution model ranges 141–243 °C, porosity 9%–12%, and permeability 0.4–0.9 mD (Figures 13c,f,i).

Figure 13

Figure 13. Property modeling maps of temperature variations, porosity, and permeability distributions. Temperature distribution model of (a) Pab Sandstone, (b) Lower Goru Formation, and (c) Chiltan Limestone. The porosity distribution model of (d) Pab Sandstone, (e) Lower Goru Formation, and (f) Chiltan Limestone. The permeability distribution model of (g) Pab Sandstone, (h) Lower Goru Formation, and (i) Chiltan Limestone.

Comprehensive evaluation index (CEI) analysis of the Pab Sandstone indicates that the anticlines of Dhodak and Rodho have higher CEI values from 0.5 to 0.75 (Figure 14a). These ranges indicate that the Pab Sandstone is a medium quality reservoir in the study area. The Lower Goru Formation has higher CEI values (0.65–0.9) along the Dhodak, Rodho, and Affiband anticlines (Figure 14b). This indicates that these locations have medium- to high-quality geothermal reservoir potential. The Chiltan limestone has the highest CEI ranges (>0.8) along the Dhodak, Rodho, Affiband, and Zindapir anticlines (Figure 14c). These ranges indicate that this formation has the potential to act as high-quality geothermal reservoir in the study area.

Figure 14

Figure 14. Comprehensive evaluation index (CEI) of selected reservoir formations. (a) Pab Sandstone, (b) Lower Goru Formation, and (c) Chiltan Limestone, indicating potential geothermal reservoir locations along the anticlines in the study area.

5 Conclusion

The following conclusions may be drawn by utilizing the integrated approaches for the identification of geothermal resource potential along the Zindapir Anticlinorium Zone.

1. Based on GIS and remote sensing analysis, the geothermal distribution map suggests that the Zindapir Anticlinorium Zone has high to very high geothermal potential, which is validated by the hot springs located in this region.

2. The geochemical analysis indicates that hot springs and associated stream waters are highly saline in nature, where the hot springs are recharged from a carbonate and clastic origin with deep circulation that caused the thermal activities. Crustal fractures help in heat flow through conduction and probably convection with some elements. The vapor-dominated phase heats the aquifers above and develops the steam-dominated zone in the research area.

3. The gravity technique suggests crustal thinning, thin-skinned faults, and the shallowest sedimentary thickness along the Zindapir Anticlinorium Zone, which indicates the uplifting of high-dense material due to compressional tectonics. The total magnetic intensity map indicates that magnetic minerals are near the surface in the south, which might be due to the shallow basement depth. Heat flow (Q) values have a direct relationship with geothermal gradients. Overall, lesser sedimentary thickness and heat flow ranges suggest that the Zindapir Anticlinorium Zone has potential geothermal resources in the subsurface.

4. The seismic interpretation suggests that Triassic to Paleocene formations are present in the subsurface. The contour map indicates that the formations deepen toward the northwest and shallow toward the south because of compressional stress. The petrophysical interpretations and property modeling of selected reservoirs indicate their potential to act as geothermal reservoirs in the study area. Lastly, based on CEI analysis using seismic, petrophysical, and property modeling, it can be concluded that the anticlines are more favorable locations for the extraction of economic geothermal reservoirs in the Zindapir Anticlinorium Zone.

Overall, the study area is classified as upper medium enthalpy, which makes it suitable for direct-use applications and for electricity generation through binary cycle systems. The current study consists of an ideal workflow that integrates surface to subsurface analyses. This comprehensive approach offers a valuable framework for identifying potential geothermal resources. It can significantly enhance academic research by providing a replicable methodology. Furthermore, the study serves as a strong foundation for initiating industry-led geothermal exploration projects.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material; further inquiries can be directed to the corresponding authors.

Author contributions

AR: Software, Investigation, Writing – review and editing, Formal analysis, Writing – original draft, Data curation, Methodology, Validation. AW: Supervision, Writing – original draft, Formal Analysis, Investigation, Writing – review and editing, Resources, Data curation, Software, Methodology, Conceptualization, Visualization, Validation. HT: Validation, Data curation, Funding acquisition, Methodology, Writing – review and editing, Project administration, Writing – original draft, Investigation. FH: Methodology, Formal Analysis, Validation, Data curation, Software, Writing – review and editing, Investigation. SA: Formal Analysis, Data curation, Methodology, Validation, Investigation, Writing – review and editing, Software. MT: Writing – review and editing, Software, Validation, Investigation, Data curation. GK: Validation, Investigation, Writing – review and editing, Data curation, Formal analysis, Visualization, Methodology. AA: Investigation, Writing – review and editing, Validation, Data curation, Methodology, Software. EB: Data curation, Validation, Writing – review and editing, Investigation, Methodology. MS: Project administration, Validation, Methodology, Data curation, Writing – review and editing, Investigation.

Funding

The authors declare that no financial support was received for the research and/or publication of this article.

Acknowledgements

The authors are thankful to the Institute of Geology, University of Azad Jammu and Kashmir, Muzaffarabad, Pakistan, for providing the facilities to complete this research. They also extend their gratitude to the Directorate General of Petroleum Concessions (DGPC) and Geological Survey of Pakistan (GSP), Ministry of Energy, Petroleum Division, Pakistan, for providing the data for this research. Schlumberger is appreciated for providing an academic license for Petrel software. The authors are also grateful to the Pakistan Institute of Nuclear Science and Technology (PINSTECH), Islamabad, for their support in sample analysis.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The authors declare that no Generative AI was used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Abbasi, T., and Abbasi, S. A. (2012). Is the use of renewable energy sources an answer to the problems of global warming and pollution? Crit. Rev. Environ. Sci. Technol. 42 (2), 99–154. doi:10.1080/10643389.2010.498754

CrossRef Full Text | Google Scholar

Agbotui, P. Y., Firouzbehi, F., and Medici, G. (2025). Review of effective porosity in sandstone aquifers: insights for representation of contaminant transport. Sustainability 17 (14), 6469. doi:10.3390/su17146469

CrossRef Full Text | Google Scholar

Ahmad, I., and Rashid, A. (2010). Study of geothermal energy resources of Pakistan for electric power generation. Energy Sources, Part A Recovery, Util. Environ. Eff. 32 (9), 826–838. doi:10.1080/15567030802606210

CrossRef Full Text | Google Scholar

Ahmad, M., Sheikh, M. R., Akram, W., Tasneem, M., Iqbal, N., and Latif, Z. (2007). Investigation of geothermal fields in Himalayan range in Pakistan using isotope and chemical techniques. Pak. Inst. Nucl. Sci. Technol. Available online at: https://docslib.org/doc/5326011/pinstech-202-investigation-of-geothermal-fields-in-himalayan-range-in-pakistan-using-isotope-and-chemical-techniques.

Google Scholar

Ahmed, S., Nazeer, A., Habib, W., Solnagi, S., and Hameed, S. (2023). Assessment of the probability of success and risking in hydrocarbon exploration: a case study from Zindapir anticline, Sulaiman Foldbelt, Indus Basin, Pakistan. J. Phys. Conf. Ser. 2594, 012002. doi:10.1088/1742-6596/2594/1/012002

CrossRef Full Text | Google Scholar

Ali, S. H., Shoukat, N., Bashir, Y., Qadri, S. M. T., Wahid, A., and Iqbal, M. A. (2022). Lithofacies and sedimentology of Baghanwala Formation (Early-Middle Cambrian), Eastern salt range, Pakistan. Pak. J. Sci. & Industrial Res. Ser. A Phys. Sci. 65 (2), 159–168. doi:10.52763/PJSIR.PHYS.SCI.65.2.2022.159.168

CrossRef Full Text | Google Scholar

Asim, S., Qureshi, S. N., Mirza, Q., Saleem, U., Ali, S., Haroon, M., et al. (2014). Structural and stratigraphical interpretation of seismic profiles along Safed koh trend (eastern part of sulaiman fold belt), Pakistan. Univers J. Eng. Sci. 2 (4), 77–95. doi:10.13189/ujes.2014.020403

CrossRef Full Text | Google Scholar

Banks, C., and Warburton, J. (1986). ‘Passive-roof’duplex geometry in the frontal structures of the Kirthar and Sulaiman mountain belts, Pakistan. J. Structural Geol. 8 (3-4), 229–237. doi:10.1016/0191-8141(86)90045-3

CrossRef Full Text | Google Scholar

Bense, V. F., Gleeson, T., Loveless, S. E., Bour, O., and Scibek, J. (2013). Fault zone hydrogeology. Earth-Science Rev. 127, 171–192. doi:10.1016/j.earscirev.2013.09.008

CrossRef Full Text | Google Scholar

Bertani, R. (2012). Geothermal power generation in the world 2005–2010 update report. Geothermics 41, 1–29. doi:10.1016/j.geothermics.2011.10.001

CrossRef Full Text | Google Scholar

Bhattacharyya, B., and Leu, L.-K. (1977). Spectral analysis of gravity and magnetic anomalies due to rectangular prismatic bodies. Geophysics 42 (1), 41–50. doi:10.1190/1.1440712

CrossRef Full Text | Google Scholar

Birkle, P., and Bundschuh, J. (2007). High and low enthalpy geothermal resources and potentials. Cent. Am. Geol. Resources Hazards 2, 705–776. doi:10.1201/9780203947043.ch26

CrossRef Full Text | Google Scholar

Blakely, R. J. (1988). Curie temperature isotherm analysis and tectonic implications of aeromagnetic data from Nevada. J. Geophys. Res. Solid Earth 93 (B10), 11817–11832. doi:10.1029/jb093ib10p11817

CrossRef Full Text | Google Scholar

Bouligand, C., Glen, J. M., and Blakely, R. J. (2009). Mapping Curie temperature depth in the western United States with a fractal model for crustal magnetization. J. Geophys. Res. Solid Earth 114 (B11). doi:10.1029/2009jb006494

CrossRef Full Text | Google Scholar

Davis, D. M., and Lillie, R. J. (1994). Changing mechanical response during continental collision: active examples from the foreland thrust belts of Pakistan. J. Struct. Geol. 16 (1), 21–34. doi:10.1016/0191-8141(94)90015-9

CrossRef Full Text | Google Scholar

Eleslambouly, A., Aldhanhani, O., Fathy, A., Zeynalli, M., and Alsuwaidi, M. (2024). “Geothermal plays feasibility and potential in a transform Basin: a case Study from Southern Vienna Basin,” in Abu Dhabi international petroleum exhibition and conference.

Google Scholar

Espinosa-Cardeña, J., and Campos-Enriquez, J. (2008). Curie point depth from spectral analysis of aeromagnetic data from Cerro Prieto geothermal area, Baja California, Mexico. J. Volcanol. Geotherm. Res. 176 (4), 601–609. doi:10.1016/j.jvolgeores.2008.04.014

CrossRef Full Text | Google Scholar

Falivene, O., Cabrera, L., and Sáez, A. (2007). Optimum and robust 3D facies interpolation strategies in a heterogeneous coal zone (Tertiary as Pontes basin, NW Spain). Int. J. Coal Geol. 71 (2-3), 185–208. doi:10.1016/j.coal.2006.08.008

CrossRef Full Text | Google Scholar

Gentana, D., Sulaksana, N., Sukiyah, E., and Yuningsih, E. T. (2019). Morphotectonics of Mount Rendingan area related to the appearances of geothermal surface manifestations. Indonesian J. Geoscience 6 (3), 291–309. doi:10.17014/ijog.6.3.291-309

CrossRef Full Text | Google Scholar

Gessner, K., Gallardo, L. A., Wedin, F., and Sener, K. (2016). Crustal structure of the northern Menderes Massif, western Turkey, imaged by joint gravity and magnetic inversion. Int. J. Earth Sci. 105, 2133–2148. doi:10.1007/s00531-016-1324-1

CrossRef Full Text | Google Scholar

Giggenbach, W. F. (1988). Geothermal solute equilibria. derivation of Na-K-Mg-Ca geoindicators. Geochimica Cosmochimica Acta 52 (12), 2749–2765. doi:10.1016/0016-7037(88)90143-3

CrossRef Full Text | Google Scholar

Gondal, I. A., Masood, S. A., and Amjad, M. (2017). Review of geothermal energy development efforts in Pakistan and way forward. Renew. Sustain. Energy Rev. 71, 687–696. doi:10.1016/j.rser.2016.12.097

CrossRef Full Text | Google Scholar

Gutiérrez-Negrín, L. C. (2024). Evolution of worldwide geothermal power 2020–2023. Geotherm. Energy 12 (1), 14. doi:10.1186/s40517-024-00290-w

CrossRef Full Text | Google Scholar

Hinze, W. J., Von Frese, R. R., and Saad, A. H. (2013). Gravity and magnetic exploration: principles, practices, and applications. Cambridge University Press.

Google Scholar

Humayon, M., Lillie, R. J., and Lawrence, R. D. (1991). Structural interpretation of the eastern Sulaiman foldbelt and foredeep, Pakistan. Tectonics 10 (2), 299–324. doi:10.1029/90tc02133

CrossRef Full Text | Google Scholar

Iqbal, M. A., Salim, A. M. A., Baioumy, H. M., Gaafar, G. R., and Wahid, A. (2017). Porosity, permeability and pore throat size distribution of Nyalau formation, Bintulu area, Sarawak basin, Malaysia. J. Appl. Environ. Biol. Sci. 7 (10), 151–158. Available online at: https://www.researchgate.net/publication/320179621_Porosity_Permeability_and_Pore_Throat_Size_Distribution_of_Nyalau_Formation_Bintulu_Area_Sarawak_Basin_Malaysia.

Google Scholar

Iqbal, M. A., Salim, A. M. A., Baioumy, H., Gaafar, G. R., and Wahid, A. (2019). Identification and characterization of low resistivity low contrast zones in a clastic outcrop from Sarawak, Malaysia. J. Appl. Geophys. 160, 207–217. doi:10.1016/j.jappgeo.2018.11.013

CrossRef Full Text | Google Scholar

Jadoon, I., Lawrence, R., and Lillie, R. (1993). “Evolution of foreland structures: an example from the Sulaiman thrust lobe of Pakistan, southwest of the Himalayas,”, 74. London, Special Publications, 589–602. doi:10.1144/gsl.sp.1993.074.01.39

PubMed Abstract | CrossRef Full Text | Google Scholar

Jadoon, I. A., Lawrence, R. D., and Lillie, R. J. (1994). Seismic data, geometry, evolution, and shortening in the active Sulaiman fold-and-thrust belt of Pakistan, southwest of the Himalayas. AAPG Bulletin 78 (5), 758–774. doi:10.1306/a25fe3ab-171b-11d7-8645000102c1865d

CrossRef Full Text | Google Scholar

Jadoon, S. U. R. K., Ding, L., Jadoon, I. A., Razzaq, S. S., Javed, M., Afridi, S., et al. (2024). Digital elevation model-based lineament analysis of the Zindapir anticline, Sulaiman fold-and-thrust belt, Pakistan. Geol. J. 59 (7), 1924–1935. doi:10.1002/gj.4977

CrossRef Full Text | Google Scholar

Janjua, O. A., Wahid, A., Salim, A. M. A., and Rahman, M. N. (2016). Prospect evaluation of shallow I-35 reservoir of NE Malay Basin offshore, Terengganu, Malaysia. AIP Conf. Proc. 1705, 020019. doi:10.1063/1.4940267

CrossRef Full Text | Google Scholar

Javed, A., Wahid, A., Mughal, M. S., Khan, M. S., Qammar, R. S., Ali, S. H., et al. (2021). Geological and petrographic investigations of the Miocene Molasse deposits in sub-Himalayas, district Sudhnati, Pakistan. Arabian J. Geosciences 14, 1–24. doi:10.1007/s12517-021-07529-x

CrossRef Full Text | Google Scholar

Kadri, I. B. (1995). Petroleum geology of Pakistan. Pak. Pet. Ltd.

Google Scholar

Karaca, S. O., Abir, I. A., Khan, S. D., Ozsayın, E., and Qureshi, K. A. (2021). Neotectonics of the Western Suleiman Fold Belt, Pakistan: evidence for bookshelf faulting. Remote Sens. 13 (18), 3593. doi:10.3390/rs13183593

CrossRef Full Text | Google Scholar

Kazmi, A. H., and Jan, M. Q. (1997). Geology and tectonics of Pakistan.

Google Scholar

Khan, N., and Scarselli, N. (2021). Seismostratigraphic architecture of the Sulaiman Fold-Thrust belt front (Pakistan): constraints for resource potential of the Cretaceous-Paleogene strata in the east Gondwana fragment. J. Asian Earth Sci. 205, 104598. doi:10.1016/j.jseaes.2020.104598

CrossRef Full Text | Google Scholar

Kurek, K. A., Heijman, W., van Ophem, J., Gędek, S., and Strojny, J. (2021). The contribution of the geothermal resources to local employment: case study from Poland. Energy Rep. 7, 1190–1202. doi:10.1016/j.egyr.2021.01.092

CrossRef Full Text | Google Scholar

Li, H., Zou, G., Zhang, K., Rong, Y., Zhu, J., Liu, Y., et al. (2023). “Selection of hybrid energy storage based on interval analytic hierarchy process,” in 2023 8th Asia conference on power and electrical engineering (ACPEE).

Google Scholar

Lund, J. W., and Toth, A. N. (2021). Direct utilization of geothermal energy 2020 worldwide review. Geothermics 90, 101915. doi:10.1016/j.geothermics.2020.101915

CrossRef Full Text | Google Scholar

Mateen, A., Wahid, A., Janjuhah, H. T., Mughal, M. S., Ali, S. H., Siddiqui, N. A., et al. (2022). Petrographic and geochemical analysis of indus sediments: implications for placer gold deposits, Peshawar Basin, NW Himalaya, Pakistan. Minerals 12 (8), 1059. doi:10.3390/min12081059

CrossRef Full Text | Google Scholar

Mir, A., Khan, M. R., Wahid, A., Iqbal, M. A., Rezaee, R., Ali, S. H., et al. (2023). Petroleum system modeling of a fold and thrust belt: a case study from the Bannu Basin, Pakistan. Energies 16 (12), 4710. doi:10.3390/en16124710

CrossRef Full Text | Google Scholar

Morata, D., Gallardo, R., Maza, S., Arancibia, G., López-Contreras, C., Mura, V., et al. (2023). Hydrothermal Alteration in the Nevados de Chillán Geothermal System, Southern Andes: multidisciplinary Analysis of a Fractured Reservoir. Minerals 13 (6), 722. doi:10.3390/min13060722

CrossRef Full Text | Google Scholar

Nazeer, A., Solangi, S. H., Brohi, I. A., Usmani, P., Napar, L. D., Jahangir, M., et al. (2013). Hydrocarbon potential of Zinda Pir anticline, eastern Sulaiman fold belt, middle Indus Basin, Pakistan. Pak. J. Hydrocarbon Res. 23, 73–84.

Google Scholar

Nazeer, A., Shah, S. H., Murtaza, G., and Solangi, S. H. (2018). Possible origin of inert gases in hydrocarbon reservoir pools of the Zindapir anticlinorium and its surroundings in the middle Indus Basin, Pakistan. Geod. Geodyn. 9 (6), 456–473. doi:10.1016/j.geog.2018.09.003

CrossRef Full Text | Google Scholar

Nyakilla, E. E., Guanhua, S., Hongliang, H., Charles, G., Nafouanti, M. B., Ricky, E. X., et al. (2025). Evaluation of Reservoir porosity and permeability from well log data based on an ensemble approach: a comprehensive study incorporating experimental, simulation, and fieldwork data. Nat. Resour. Res. 34 (1), 383–408. doi:10.1007/s11053-024-10402-9

CrossRef Full Text | Google Scholar

Powell, C. M. (1979). A speculative tectonic history of Pakistan and surroundings. Geodyn. Pak.

Google Scholar

Qadri, H. A., Wahid, A., Siddiqui, N. A., Ali, S. H., El Aal, A. A., Rashid, A. Q. B. A., et al. (2022). Prospect analysis of Paleocene coalbed methane: a case study of hangu formation, trans-indus ranges, Pakistan. Geofluids 2022 (1). doi:10.1155/2022/8313048

CrossRef Full Text | Google Scholar

Quittmeyer, R. C., Kafka, A. L., and Armbruster, J. G. (1984). Focal mechanisms and depths of earthquakes in central Pakistan: a tectonic interpretation. J. Geophys. Res. Solid Earth 89 (B4), 2459–2470. doi:10.1029/jb089ib04p02459

CrossRef Full Text | Google Scholar

Rashid, A., Siddiqui, N. A., Ahmed, N., Jamil, M., EL-Ghali, M. A., Ali, S. H., et al. (2022). Field attributes and organic geochemical analysis of shales from early to middle Permian dohol formation, peninsular Malaysia: implications for hydrocarbon generation potential. J. King Saud University-Science 34 (8), 102287. doi:10.1016/j.jksus.2022.102287

CrossRef Full Text | Google Scholar

Raza, H. A., Ahmed, R., Ali, S. M., and Ahmad, J. (1989). Petroleum prospects: sulaiman sub-basin, Pakistan. Pak. J. Hydrocarbon Res. 1 (2), 21–56.

Google Scholar

Reynolds, K., Copley, A., and Hussain, E. (2015). Evolution and dynamics of a fold-thrust belt: the sulaiman Range of Pakistan. Geophys. J. Int. 201 (2), 683–710. doi:10.1093/gji/ggv005

CrossRef Full Text | Google Scholar

Rezaie, M., and Aghajani, H. (2013). A new combinational terminology for geothermal systems. Int. J. Geosciences 4 (1), 43–48. doi:10.4236/ijg.2013.41005

CrossRef Full Text | Google Scholar

Ross, H., Blakely, R., and Zoback, M. (2004). Testing the utilization of aeromagnetic data for the determination of Curie-isotherm depth. AGU Fall Meet. Abstr.

Google Scholar

Saemundsson, K., Axelsson, G., and Steingrímsson, B. (2009). Geothermal systems in global perspective. Short course Exploration Geothermal Resources, UNU GTP 11.

Google Scholar

Saif-Ur-Rehman, K. J., Ding, L., Jadoon, I. A., Baral, U., Qasim, M., and Idrees, M. (2019). Interpretation of the Eastern Sulaiman fold-and-thrust belt, Pakistan: a passive roof duplex. J. Struct. Geol. 126, 231–244.

CrossRef Full Text | Google Scholar

Salim, A. M. A., Janjua, O. A., Wahid, A., and Zaheer, A. (2015). “Linear regression approach for porosity and permeability calculations from well logs: a case Study in NW Bonaparte Basin, Australia,” in 4th Annual International Conference on Geological & Earth sciences (GEOS 2015). Singapore: Abstract Book.

Google Scholar

Sano, Y., and Wakita, H. (1985). Geographical distribution of 3He/4He ratios in Japan: implications for arc tectonics and incipient magmatism. J. Geophys. Res. Solid Earth 90 (B10), 8729–8741. doi:10.1029/jb090ib10p08729

CrossRef Full Text | Google Scholar

Sarwar, G. (1979). Arcs, oroclines, syntaxis, the curvatures of mountain belts in Pakistan. Quetta, Pakistan: Geological Survey of Pakistan, 341–349.

Google Scholar

Shuey, R., Schellinger, D., Tripp, A., and Alley, L. (1977). Curie depth determination from aeromagnetic spectra. Geophys. J. Int. 50 (1), 75–101. doi:10.1111/j.1365-246x.1977.tb01325.x

CrossRef Full Text | Google Scholar

Solangi, S., Abbasi, S., Ali, A., Asim, S., Lashari, R., Brohi, I., et al. (2014). Study of subsurface structural trend and stratigraphic Architecture using seismic Data-A case Study from Zindapir inner folded Zone, sulaiman sub-Basin, Pakistan. Sindh Univ. Reserv. J. 46, 373–380.

Google Scholar

Soltani, M., Kashkooli, F. M., Souri, M., Rafiei, B., Jabarifar, M., Gharali, K., et al. (2021). Environmental, economic, and social impacts of geothermal energy systems. Renew. Sustain. Energy Rev. 140, 110750.

CrossRef Full Text | Google Scholar

Spector, A., and Grant, F. (1970). Statistical models for interpreting aeromagnetic data. Geophysics 35 (2), 293–302. doi:10.1190/1.1440092

CrossRef Full Text | Google Scholar

Swanberg, C. A., and Morgan, P. (1978). The linear relation between temperatures based on the silica content of groundwater and regional heat flow: a new heat flow map of the United States. Pure Appl. Geophys. 117 (1), 227–241. doi:10.1007/978-3-0348-6525-8_20

CrossRef Full Text | Google Scholar

Tanaka, A., Okubo, Y., and Matsubayashi, O. (1999). Curie point depth based on spectrum analysis of the magnetic anomaly data in East and Southeast Asia. Tectonophysics 306 (3-4), 461–470. doi:10.1016/s0040-1951(99)00072-4

CrossRef Full Text | Google Scholar

van der Meer, F., Hecker, C., van Ruitenbeek, F., van der Werff, H., de Wijkerslooth, C., and Wechsler, C. (2014). Geologic remote sensing for geothermal exploration: a review. Int. Journal Applied Earth Observation Geoinformation 33, 255–269. doi:10.1016/j.jag.2014.05.007

CrossRef Full Text | Google Scholar

Wahid, A., Salim, A. M. A., Gaafar, G. R., and Yusoff, A. W. I. W. (2016a). Tectonic development of the Northwest Bonaparte Basin, Australia by using digital elevation Model (DEM). IOP Conf. Ser. Earth Environ. Sci. 30, 012005. doi:10.1088/1755-1315/30/1/012005

CrossRef Full Text | Google Scholar

Wahid, A., Salim, A. M. A., Gaafar, G. R., and Yusoff, W. I. W. (2016b). The geostatistical approach for structural and stratigraphic framework analysis of offshore NW Bonaparte basin, Australia. AIP Conf. Proc. 1705, 020030. doi:10.1063/1.4940278

CrossRef Full Text | Google Scholar

Wahid, A., Rauf, A., Ali, S. H., Khan, J., Iqbal, M. A., and Haseeb, A. (2022). Impact of complex tectonics on the development of geo-pressured zones: a case study from petroliferous sub-himalayan basin, Pakistan. Geopersia 12 (1), 89–106. doi:10.22059/GEOPE.2021.324799.648618

CrossRef Full Text | Google Scholar

Wajid, A. A., Anees, M., Alam, S. u., Gorchani, J. K., Shahzad, K., Israr, A., et al. (2021). Lineament mapping for a part of the central sulaiman fold–thrust belt (SFTB), Pakistan. Arabian J. Geosciences 14, 1–19. doi:10.1007/s12517-021-07784-y

CrossRef Full Text | Google Scholar

Wang, G., and Kuang, J. (2022). Genetic analysis of geothermal resources in deep-seated fault area in Tonghe county, Northeast China and implications of geothermal exploration. Sustainability 14 (9), 5431. doi:10.3390/su14095431

CrossRef Full Text | Google Scholar

Wang, Y., Ma, F., Xie, H., Wang, G., and Wang, Z. (2021). Fracture characteristics and heat accumulation of Jixianian carbonate reservoirs in the Rongcheng geothermal field, Xiong'an new area. Acta Geol. Sinica-English Ed. 95 (6), 1902–1914. doi:10.1111/1755-6724.14878

CrossRef Full Text | Google Scholar

Wishart, D. N. (2020). “Binary cycle plant design for water-dominated, low enthalpy geothermal system,” in Proceedings world geothermal congress, 1.

Google Scholar

Wu, Y., and Zhou, X. (2023). Structural control effects on hot springs’ hydrochemistry in the northern Red River Fault zone: implications for geothermal systems in fault zones. J. Hydrology 623, 129836. doi:10.1016/j.jhydrol.2023.129836

CrossRef Full Text | Google Scholar

Xu, Z., Sun, F., Xin, W., Sun, N., Li, F., Niu, J., et al. (2019). Formation and evolution of paleoproterozoic orogenic belt in southern Jilin, Jiao–Liao–Ji Belt, North China craton: constraints from geophysics. Precambrian Res. 333, 105433. doi:10.1016/j.precamres.2019.105433

CrossRef Full Text | Google Scholar

Xu, Z.-H., Sun, Z.-J., Xin, W., and Zhong, L. (2021). Geothermal resource potential assessment of Erdaobaihe, Changbaishan volcanic field: constraints from geophysics. Open Geosci. 13 (1), 1053–1063. doi:10.1515/geo-2020-0282

CrossRef Full Text | Google Scholar

Zaidi, S. I. A., Ahsan, N., and Rehman, S. U. (2021). Geochemical attributes of Pliocene Dhok pathan formation, Eastern Sulaiman Fold and Thrust Belt, Pakistan: implications for provenance, paleoenvironmental, and tectonic settings. Arabian J. Geosciences 14, 1–30. doi:10.1007/s12517-021-08657-0

CrossRef Full Text | Google Scholar

Zaigham, N. A., Nayyar, Z. A., and Hisamuddin, N. (2009). Review of geothermal energy resources in Pakistan. Renew. Sustain. Energy Rev. 13 (1), 223–232.

CrossRef Full Text | Google Scholar