Geothermal exploration in Gujarat: case study from Dholera

Abstract

Utilization of Geothermal resources is an alternate, sustainable energy, which needs to be harnessed in India . Exploration of the resource is the first step to understand potential leads and prospects in the subsurface. An attempt has been made to understand the subsurface picture using geoscientific data, for example, Gravity survey, Landsat imagery, magnetotelluric (MT) survey and water chemistry. Dholera, Gujarat is a potential site for geothermal exploration and exploitation, which is the chosen study area. Dholera also falls under proposed Delhi Mumbai Industrial Corridor Development Corporation Limited (DMICDC), making it a suitable candidate for setting up geothermal power plant in future. The energy potential of prospects identified in Dholera through the above-mentioned studies may help industries in creating geothermal power plant and utilization of energy in space heating & district heating/cooling. Dholera is located 30 km away in southwest direction from Dhandhuka village of Ahmedabad district and 60 km away in north direction from city of Bhavnagar (Fig. 1c). Hot springs exist over gravity high, which is the surface manifestation of deep and shallow water sources. Present research focuses on vegetation index for the study area and land surface temperature; alteration of water chemistry in association with subsurface rocks and understanding of subsurface models using Gravity and MT surveys. 2D geoelectric maps and gravity contours have helped in understanding the aquifer geometry and dimensions. 2D resistivity sections depict the aquifer connectivity across Dholera at depths of 3 and 4 km. MT survey suggests that the volcanic rocks are exposed at places. However, it also depicts locations where sediments exist and are sandwiched between basaltic lavas. MT model is also constrained using residual gravity data and corroborates well with the MT data. Present study will help in harnessing the geothermal energy entrapped within the geobodies and also to understand connectivity between the geobodies. These studies were carried out one after the other, to delineate geothermal reservoir and increase the probability of success. The best location for further exploration and exploitation by drilling has been identified by integration of all these methods.

Sircar et al. Geotherm Energy (2015) 3:22 DOI 10.1186/s40517-015-0041-5

Geothermal exploration inGujarat: case study fromDholera

Anirbid Sircar1*, Manan Shah1, Shreya Sahajpal1, Dwijen Vaidya2, Shubhra Dhale2 and Anjali Chaudhary2

*Correspondence: [email protected]

1 School of Petroleum Technology, Pandit Deendayal Petroleum University, Gandhinagar, Gujarat 382007, IndiaFull list of author information is available at the end of the article

Background

Dholera is located 30 km away in southwest direction from Dhandhuka village of Ahmedabad district and 60km away in north direction from city of Bhavnagar. Dholera thermal springs are located along the margin of Saurashtra Peninsula falling in the vicinity of Western Marginal fault of Cambay Basin as shown in Fig.1c (Sharma 2013).

2015 Sircar et al. This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/

Web End =http:// http://creativecommons.org/licenses/by/4.0/

Web End =creativecommons.org/licenses/by/4.0/ ), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Sircar et al. Geotherm Energy (2015) 3:22

Page 2 of 25

Terrain in Dholera is covered by recent to alluvium and mud ats. The area is also occupied by quaternary soils deposited in subsiding area by the side of Cambay Basin to a thickness of about 100m over Tertiary sediments resting on Deccan Traps at a depth of about 500600m. The general slope of the ground is from west to east in the investigation area. The study area falls under western margin of the Cambay Basin. In Dholera we have presence (Sharma 2013) of old mud ats, ood plains and salt at areas. The soil in Dholera region mainly consists of alternate layers of gravels, ne to coarse grained sand and clay. Chemically the soil is loamy, mixed montmorillonitic, calcareous and mostly saline. Dholera springs are also located over the high gravity area (400400m) indicating a shallow mantle and presence of granitic basement rocks (Sharma 2013). A total of four springs were demarcated in a radius of 4km Dholera, Uthan, Swaminarayan temple and Bhadiyad. Dholera springs have the highest geothermal ow rate in Gujarat (Vaidya etal. 2015).

This work attempts to study the geothermal reservoir existing in shallow and deep crustal levels. The reason for believing that such reservoir exists is because of the fact that the region has registered many shallow level hot springs. Hot springs are the surface manifestations of these reservoirs. Present study is supported by remote sensing (RS) analysis, physical and chemical characterization of thermal springs, audio magnetotelluric (AMT) survey, magnetotelluric survey (MT) and gravity surveys. In Dholera no major geothermal studies have been carried out. RS analysis followed by collection of water samples to determine major chemical parameters and trace elemental analysis. Physical parameters like color, pH and conductivity of water samples were measured insitu and in the laboratory.

Sircar et al. Geotherm Energy (2015) 3:22

Interpretation of geoscientic data helped us to identify resistivity distribution below the subsurface and also the subsurface structures and anomalies. It also identies the area covered by low resistivity below subsurface.

Tectonic frame work ofmainland Gujarat

Gujarat state carries distinction of diverse and unique combination of geological, geo-morphological and physiographic attributes as shown in Fig.1b. The interplay of geologic, tectonic and climate conditions have resulted in three physiographic domains in the region, i.e., Mainland Gujarat, Saurashtra Peninsula and Kachchh Peninsula. Biswas (1987) has coined names for the three depositional basins, namely Kachchh, Cambay and Narmada (Fig.1b), representing the physiographic divisions of land formed due to the events of rifting during TriassicEarly Jurassic, Early Cretaceous and Late Cretaceous times.

Three dominant tectonic trends have controlled Saurashtra Peninsula; viz, Delhi trend (NESW), SonNarmadaTapti (SONATA) trend (ENEWSW) and Dharwar trend (NNWSSE) (Fig.1a). These are regional trends and have conned the Saurashtra Peninsula into a huge horst comprising mainly of the Deccan ows with contemporary intrusion of alkaline nature as well as basalt dierentiates (Biswas 1987).

Saurashtra Peninsula has strong tectonic control and is bounded on all sides by well-established faults, namely the Kim fault in the shelf of Saurashtra coast and Marginal faults in directional continuity with Narmada rift system trending ESEWNW. The western coastal margin is dened by the NWSE trending marginal fault; northern margin is controlled by the North Kathiawar fault. Over and above these there are many prominent lineaments trending NWSE, NS and EW hinting to subsurface fault systems (Biswas 1988). The important thermal springs Tulsishyam and Savarkundla are along EW and NESW lineaments, respectively. Lalpur is in the proximity of the North Kathiawar fault. The high thermal gradients in these areas can be attributed to shallow mantle and radiogenic heat from granitic basement rocks (Merh 1995).

General stratigraphy ofarea

There are two geological columns (stratigraphy) shown in Fig.2, namely Saurashtra Peninsula and Cambay Basin. Dholera region is located in Saurashtra Peninsula and adjacent to Cambay basin. To understand the overall subsurface chrono-stratigraphic events of the study area it is required to study both stratigraphic columns.

Figure2 shows the generalized stratigraphy of Saurashtra Peninsula in the right column and North Cambay Basin in left column (Kulkarni 1985; Merh 1995), which exhibits chrono-stratigraphic markers. The formation of Cambay Basin began with extensive outpour of Deccan Trap of Late Cretaceous age. During Paleocene, the basin continued to a shallow depression and the end of deposition of the Olpad formation is marked by a prominent unconformity. During Early Eocene, the deposition of thick, dark gray shales known as Cambay shale took place which has been divided into Older and Younger Cambay shale with an unconformity in between. During Late Eocene to Early Oligocene, the deposition of Tarapur shale took place over large area of the North Cambay Basin. During Miocene, the depocenters continued to decrease the deposition of sediments such as Alluvium, silts, gravels, sands and clay.

Page 3 of 25

Sircar et al. Geotherm Energy (2015) 3:22

Page 4 of 25

The Saurashtra Basin is mainly a clastic basin and development of limestone is noticed over Saurashtra arch (Fig. 1c). Limestone sequences of Eocene to Miocene age attain considerable thickness. The Proterozoic rocks are not exposed here. The Cretaceous Jurassic marine and uvio-marine sedimentaries overlie the basement granites and other Proterozoic rocks. There was a period of non-deposition represented by unconformity between these sedimentary rocks and Trappean ows laid during the Upper Cretaceous and Lower Eocene. During Palaeocene, presence of laterites over these sedimentaries and volcanic rocks conrms the period of non-deposition. The Upper Tertiary rocks were deposited in the peripheral basins of the Peninsula followed by the Quaternary deposits of varied depositional environments. At some places in Dholera region, upper formations have been eroded and Deccan traps can be encountered at very shallow depths.

Methods

Remote sensing analysis

Remote sensing, also known as earth observation, enables us to obtain information about objects or areas on Earths surface without being in direct contact with the object or area (Calvin etal. 2005). RS is used in previous exploratory phases in the study area,

Sircar et al. Geotherm Energy (2015) 3:22

as well as in geothermal exploration. This would pave the way for detailed geochemical and geophysical surveys for possible delineation and exploitation of geothermal resources (Shah etal. 2015). RS techniques are emerging as useful reconnaissance tools for mapping geology, detecting anomalous structure, surface temperature anomalies and identifying geothermal indicators including hydrothermal alteration of minerals (clays, sulfates), sinter and tu in prospective geothermal areas (Calvin etal. 2005).

The goal of RS is to locate surface indicators, which may indicate suitable traps at depth, lineament mapping, thrust zone mapping, micro-seepage maps of hydrocarbons, and stain maps of hydrocarbons. It also gives information about the location of dierent hot springs.

Detailed RS analysis was carried out to check land surface temperature and identify hot springs in Gujarat (Fig.3). Split window algorithm has been applied during night, moderate-resolution imaging spectro-radiometer (MODIS) imagery for last 1 year to understand regional-scale scenario of land surface temperature anomaly, which in turn helps us in locating geothermal anomalous zones (Calvin etal. 2005). Nighttime land surface temperature (LST) data were used so that sun reection does not aect the observation and we get actual emissivity data from the earth surface (Huang and Pollack 1997). During daytime earths emissivity is very high all over the surface like deserts, cities, barren land, etc., but at night time real values can be obtained as the earth surface temperature comes into equilibrium. The result is the heat of earth, which is measured. It is then cross-checked with LST values with reference to Normalize Dierence Vegetation Index (NDVI) map as shown in Fig.3 (Srivastava and Gupta 2013).

Researchers categorized NDVI data (Srivastava and Gupta 2013) into classes ranging from no vegetation to high vegetation. NDVI is a numerical indicator, which indicates the presence of green vegetation in the target area using visible and near-infrared bands of the electromagnetic spectrum. Low to moderate vegetation with high LST regions was understood to have chance of geothermal potential. Hence, Dholera having low to

Page 5 of 25

Sircar et al. Geotherm Energy (2015) 3:22

moderate vegetation and high LST is considered for further investigation based on RS prospects.

Geochemical analysis

Study was undertaken to measure physical and chemical characterization of the thermal spring water from Dholera hot springs. The main objective of the study was to detect and explain the possible changes in chemical composition and properties of thermal uids over time. Geochemical analysis is useful to evaluate the reservoir temperatures and equilibrium conditions using geothermometers (Pasvanoglu 1998).

For geochemical analysis, samples from thermal springs and ground water were taken from Dholera. These samples were analyzed for color, pH, hardness, salinity, conductivity and other physical properties (Table1). Detailed chemical analysis was carried out to check the presence of cations (Ca2+, Mg2+, Na+, K+) and anions (CO3, HCO3, Cl,

SO42) presented in Table 1. Table 1 also depicts the elemental chemistry of thermal springs.

Dholera water samples have shown moderate to high salinity values. The salinity value ranges from 4500 to 5000mg/l. Water chemistry inherits high salinity owing to the water interaction with granitic basement having high Na+, K+, Cl (Sharma 2013). Conductivity values are moderate in nature. Total dissolved solids (TDS) values are in agreement with the salinity values. Lower carbonate hardness is understandable, since the dominant anion species in these samples are chloride and sulfate. Water category is sodium, potassium, chloride and sulfate type. Low levels of carbonate and bicarbonate components indicate that the water has risen up from deeper levels. Boron is detected in water which indicates that it has mixed components from the formation of older Mesozoic sediments. Fluoride is present in all samples in variable concentration levels. High uoride levels indicate that these waters are in interaction with subsurface mica and apatite-bearing rocks. Silicon (Si) present in sample varies from 12 to 5mg/l (Sharma 2013). Si contribution is expected to happen at deeper levels, where the thermal gradient is high. Si gets dissolute from alkali feldspar, quartz and amorphous silica.

It can be concluded that these waters have meteoric origin. Due to the regional tec-tonic structure and up-welled mantle condition, high temperature gradients have heated

Table 1 Chemical andphysical parameters withtrace elemental chemistry ofthe thermal water samples fromDholera (Sharma 2013)

Sr. no. Chemical parameters Physical parameters Trace elemental chemistry

Parameter Value Parameter Value Parameter Value

1 CO3 (mg/l) Odor B mg/l 4

2 Cl (mg/l) 3798 pH 7.9 F mg/l 0.7

3 HCO3 (mg/l) 72 Conductivity Ms/cm 8370 SiO2 mg/l 11.48

4 SO4 (mg/l) 10.47 Salinity mg/l 6860 Pb mg/l 1

5 Ca (mg/l) 124.5 TDS mg/l 5608 Ni mg/l 0.2

6 Mg (mg/l) 31.6 CO3 Hardness mg/l 144

7 Na (mg/l) 2000

8 K (mg/l) 18

Page 6 of 25

Sircar et al. Geotherm Energy (2015) 3:22

this water during their downward movement along deep channel ways of tectonic origin. The water has risen back to the surface after reaching the boiling point.

Gravity survey

Gravity method is a passive technique that involves the measurement of the acceleration due to the Earths gravitational eld (Rivas 2009). The variations in gravity are due to density lateral change of the subsurface rocks (Gupta and Roy 2007).

The objective of the gravity survey was to identify anomalies having density dierences between anomaly and background (Mickus etal. 1991). The data collected were corrected for surcial elevations and tilt of the gravimeter. Gravity values collected were base referenced and was plotted on a 2D map. Once necessary corrections (free air, latitude, Bouguer, etc.) were applied to the data, Bouguer gravity is separated into regional and residual. The regional was assumed to be a 2D spline which represents the tectonic framework (Cooper et al. 2010). Dierence between Bouguer anomaly and regional gravity provided the residual Bouguer gravity. Residual gravity interpretation gives an understanding of the subsurface anomaly. Data were collected along the MT proles or on its osets. Gravity values were contoured and shown in Fig.4a. The transverse fault bisects the observed gravity low in the Bouguer gravity map (Fig.4a). This shows that there is a less dense body beneath the surface, perhaps an indication of highly conductive body below the subsurface (Mariita 2007). AA is a gravity prole prepared using corrected Bouguer gravity data (Fig.4b). The regional trend is tted to the data (2D spline) to prepare the residual gravity data (Fig.4c). The residual gravity data were inverted and Fig.4d shows the match between calculated and inverted gravity data. Figure suggests that there is a good t between calculated residual gravity data and inverted gravity data. Coefficient of correlation between residual and tted gravity data is found to be 0.89.

Upward continuation (UC) ltering was carried out on the gravity data (Fig.5). UC serves to smooth out near-surface eects (shorter wavelength anomalies) after calculating the gravity eld at an elevation higher than that at which gravity is measured. On computation of Bouguer residual gravity data, six parallel lines were selected to get the geologic model. Residual Bouguer gravity was inverted using the algorithm of Talwani etal. (1959) and a contrast density of 0.2g/cm3 between basalts and shale. The geologic model is shown in Fig.5 with control points. Basaltic lava has a density of 2.4g/cm3

whereas the anomalous body has a lower density of 2.2g/cm3. Geologic model corroborates well with the MT model described later.

Magnetotelluric method (MT)

MT method is a passive surface measurement of earths natural electric eld and magnetic eld in orthogonal directions. It can be seen that relationship between the horizontal orthogonal magnetic and electric eld depends on the subsurface resistivity structure. Based on sample analysis, AMT and MT survey were performed at 66 sounding stations along with six MT proles shown in Fig.6. MT generally refers to recording of 100.001kHz or as low as 0.0001Hz (Bahr and Simpson 2005; Vozo 1991). Gravity measurement was conducted at several places, either on MT stations or at osets. It is, therefore, used to determine the conductivity of earth ranging from a few tens of meters to several hundreds of kilometers (Sahajpal etal. 2015).

Page 7 of 25

Sircar et al. Geotherm Energy (2015) 3:22

Page 8 of 25

MT measurement allows detection of resistivity anomalies associated with geothermal structures, faults and the presence of a cap rock, and allows for estimation of geothermal reservoir temperatures at various depths. Geological materials are generally poor electrical conductors and have a low conductivity. Thus, change in conductivity is useful to infer the subsurface geology and estimate the subsurface material composition (Constable etal. 1987).

Sircar et al. Geotherm Energy (2015) 3:22

Page 9 of 25

Source elds ofMT signals

The MT signals are generated from two sources:

1. First source is featured by lower frequencies, generally less than 1Hz. Origin of these signals is from the interaction of solar wind with the earths magnetic eld. Solar wind emits highly charged ions which travel into space and disturb earths magnetic eld and produce low-frequency electromagnetic energy that penetrates in the solid earth (Cagniard 1953).

Sircar et al. Geotherm Energy (2015) 3:22

2. Second source is characterized by higher frequencies, generally greater than 1Hz. It is created by thunderstorm activities. Energy created by these thunderstorms travels around the earth in a wave form between the earths surface and the ionosphere, with part of energy penetrating into the earth (Keller and Frischknecht 1966).

MT acquisition

Field AMT and MT measurements were performed in Dholera at 66 MT/AMT sounding stations along six proles (Fig.6). Orientation of ve proles was WSWENE and one normal to ve proles. Frequency of the MT/AMT data is in the range of 0.001 10,000Hz. Simultaneously synchronized measurements on reference station located in Kamalpura (Fig.1c) were carried out.

Remote reference in MT method is used to eliminate or reduce the eects of articial electromagnetic noise measurements. MT/AMT eld working along these six proles were made using system2000.net of Phoenix Geophysics, Canada, based on electric receivers: V8-6R, RXU-3E, MTU-5A and magnetic sensors: AMTC-30, MTC-50 and MTC-50H. Measuring array comprises two perpendiculars Ex, Ey electric dipoles of 100m length and three magnetic sensors recording magnetic eld for horizontal components Hx, Hy and vertical Hz. V8 and MTU receivers and sensors were calibrated insitu before survey.

The acquired electric and magnetic vectors were converted to apparent resistivity curves with respect to frequency. MT/AMT raw data (time series les) acquired by geophysical instruments were recorded on compact ash card and stored on PC hard disk, DVD and nally processed (PBG Report 2014).

The variations of electromagnetic eld were recorded in a broadband. The data were acquired in four frequency bands named as ts2, ts3, ts4 and ts5. Sampling rate was as following:

24kHz for ts2 band

2400Hz for ts3 band

150Hz for ts4 band

15Hz for ts5 band

Bands ts2, ts3 and ts4 are used in AMT (from 10kHz to 1Hz recorded with AMTC-30 magnetic coils), bands ts3, ts4 and ts5 correspond to MT (from 320 to 0.0005 Hz recorded with MTC-50, MTC-50H magnetic coils).

Data processing

MT/AMT data have been recorded during eld measurements and processing was carried out with the use of professional software SSMT2000 and MT-Editor. The SSMT2000 program takes the raw time series les, calibration les and site parameter les as input. In an intermediate step, it produces Fourier coefficients, which are then reprocessed with data from reference sites, using robust routines. The outputs are MT Plot les containing multiple cross-powers for each of the frequencies analyzed. MT-Editor program takes as input the MT Plot les created by SSMT2000 and displays resistivity and phase

Page 10 of 25

Sircar et al. Geotherm Energy (2015) 3:22

curves as well as the individual cross-powers that are used to calculate each point on the curves.

Dierent plots were prepared such as apparent resistivity, phase, impedance, tipper magnitude, coherency between channels, and strike direction across the full frequency range that was acquired.

Results anddiscussion ofmagnetotelluric data

MT/AMT data interpretation

1D Occam algorithm is a typical 1D inversion method, in which classical layered model is obtained. For better visualization, such layered models are usually presented after being interpolated in cross-sectional manner. Occam algorithm gives smooth models with sharp boundaries between geoelectrical layers and leads to a simple model containing the essential properties of all possible models tting the eld data (Constable etal. 1987). A large number of geoelectrical models could match the observed data, some of which may be highly complex. Resistivities of layers vary between eld and calculated curves.

2D smooth model inversion routine nds regularized solutions (Tikhonov Regularization) to the 2D inverse problem for MT data using the method of non-linear conjugate gradients. Forward model simulations are computed using nite dierence equations generated by network analogs to Maxwells equations. The program inverts the data for user-dened 2D mesh of resistivity blocks, extending laterally and downwards beyond the central detailed zone, and incorporating topography (Rodi and Mackie 2001).

Qualitative interpretation

Before quantitative interpretation of MT curves, a quality control, dimensionality analysis and editing is usually made. There are several parameters to analyze dimensionality such as skew, tipper magnitude and polar diagrams. The apparent resistivity of three MT soundings D2_12, D2_14 and D2_16 is shown in Figs.7, 8 and 9 respectively with phase curve, tipper magnitude, azimuth and skew.

Skew is a measure of ratio of amplitudes of diagonal impedance elements to the o-diagonal impedance elements. Lower values, usually less than 0.3 are typical for 1D and 2D data. Tipper normally indicates lateral variations and can be used as an indicator of presence of 2D and 3D structure.

Based on the analysis of distribution of skew of impedance tensor, it can be said that for almost whole frequency band, the survey area is characterized by the geological structure equivalent to 1D or 2D geoelectrical model (skew values for the whole area are less than 0.3).

Frequencies greater than 0.001Hz tends to be shallower, thus from the depth corresponding to a frequency of 0.001Hz up to the surface the subsoil would be 3D. Polar diagrams hardly help to dierentiate between 2D and 3D terrain.

Analysis of the shape of polar diagram provides information about the level of 3D distortion and/or noise that may occur within data. For 1D geoelectric structure, the principal impedance polar diagrams are circular in shape. For 2D or 3D structure, they elongate in a direction either parallel or perpendicular to strike, depending on the position of the node with respect to the heterogeneity.

Page 11 of 25

Sircar et al. Geotherm Energy (2015) 3:22

Page 12 of 25

By analyzing polar diagrams, from frequencies 10kHz to about 0.107Hz it can be considered as 1D models where impedance Zxy is almost circular and Zxx is nearly invisible (Fig.10). For frequency 0.011Hz (Fig.11), it is generally considered as 2D model. Polar diagrams (Zxy) are elongated in the NNWSSE direction (azimuth about 23 degrees)

and Zxx is small. Shapes of the polar diagrams indicate that there are two possible axes of homogeneity. It is approximately parallel to Deccan Trap and Western Marginal Fault. For frequency 0.0011Hz (Fig.12), the polar diagram shows 3D model.

Quantitative interpretation1D Occam interpretation ofDholera prole1D Occam inversion was conducted considering that earth is composed of uniform stratied layers (Constable etal. 1987). Resistivity of layers is varied until a good t is

Sircar et al. Geotherm Energy (2015) 3:22

Page 13 of 25

obtained between the calculated and observed data. 2D inversion is carried out using non-linear conjugate gradients (Rodi and Mackie 2001). Forward model simulations are computed using nite dierence equations generated using Maxwells equation or computed using Divided Dierence method using Taylor series ( PBG Report 2014).

Data interpretation was based on 1D inversion Occam algorithm using WinGLink software. The results of inversion algorithm by Occam have been interpolated and are shown as resistivity section. Distribution of resistivity is presented for two types of models: shallow up to 1000m b.s.l. and deep up to 10km b.s.l.

Figure13 is the 1D deep resistivity cross section along the second prole (D2) for deep imaging. It can be described as three-layered geoelectrical model. The rst geoelectrical layer has resistivity of about 17ohmm, which is clearly visible in the shallow part of the model depicted in Fig.14. The thickness of this layer is about 200m in WSW and

Sircar et al. Geotherm Energy (2015) 3:22

Page 14 of 25

becomes 1000m in ENE. Resistivity of the second layer is more than 150ohmm. Minimum thickness of this layer occurs in WSW and increases towards ENE. The third layer has very high-resistivity values, which is interpreted as the basement.

2D geoelectric deep models

Resistivity distribution for geoelectric deep model is up to 10 km below sea level as shown in Figs.15, 16, 17, 18 and 19. Resistivity models can be described as two layered. Bottom of the rst geoelectrical layer has resistivity of about 16ohmm which is indicated by thick dark line in Figs.15, 16, 17, 18 and 19. The thickness of this layer is about 200m in the western parts of the research area and rises up to about 1000m in the NE

Sircar et al. Geotherm Energy (2015) 3:22

Page 15 of 25

part of area. Below this layer there are some high-resistive anomalies. These anomalies are characterized by reduced resistivity up to 400ohmm. This layer can be linked to Quaternary and Tertiary sediments, mixed with salt water. Generally, sediment thickness decreases in the direction SE where Deccan volcanic rocks appear on the surface, which are clearly visible in the cross sections shown in Figs.15, 16, 17, 18 and 19. Drilling a borehole will help in validating this inference. It will help to determine the morphology on the basis of isolines resistivity from 1D or 2D inversion.

Figure15 shows a resistivity cross section based on 2D NLCG (non-linear conjugate gradient) (PBG Report 2014) inversion along prole D1 (Deep Model). From station D2_14 to D2_22 in ENE direction and from station D2_8 to D2_02 in WSW direction at depth of 3000m onwards high-resistivity (4001000ohmm) zones can be seen, which may be an indication of high-density basement of Precambrian basalts. These high-resistivity zones can be seen along all ve parallel proles (D1D5) from Figs.15, 16, 17, 18 and 19, at a depth of up to 3000m, low-resistivity zones (up to 75ohmm) are visible. It might be an indication of low-density sedimentary formations. A signicant

Sircar et al. Geotherm Energy (2015) 3:22

Page 16 of 25

geophysical anomaly having resistivity of 1000ohmm is clearly visible starting from a depth of 4000m. In the second prole, with increasing depth there is a decrease in resistivity. At 7000m depth ultralow-resistivity zone (16ohmm) is depicted. Here, it shows a strong possibility of magmatic reservoir, which is shown as a hatched zone in Fig.15. The temperature of the reservoir can be estimated based on the geothermal gradient of the respective region, which is 3.5C/100m in this case. However, actual temperature can be only known by drilling a geothermal parametric well. This acts as a source of hot water, which eventually migrates to surface in the form of hot springs.

Figure16 shows a resistivity cross section based on 2D NLCG inversion along prole D2 (Deep Model) (PBG Report, 2014). Along prole D2, below stations D1_16 and D1_18, at a depth of 5000m and onwards, a signicant geophysical anomaly is visible, which is an ultralow-resistivity zone. Again, similar to prole D1 it may be an indication of a hot uid or magmatic body. Furthermore, in this case the anomaly is surrounded by high-resistivity zones, i.e., high-density basalts (traps). This might restrict the lateral ow of hot water.

Figure17 indicates resistivity cross section based on 2D NLCG inversion along prole three (Deep Model). It shows similar characteristics as prole D2 (central prole passing through hot springs). Here, also low-resistive geophysical anomaly is sandwiched

Sircar et al. Geotherm Energy (2015) 3:22

Page 17 of 25

Sircar et al. Geotherm Energy (2015) 3:22

Page 18 of 25

between high-resistive basements of over 1000ohmm. It may also be an indication of a geothermal system similar to central prole.

Figure18 shows a resistivity cross section based on 2D NLCG inversion along prole D4. Along prole D4, a geophysical anomaly is found at 6000m depth between stations D4_6 and D4_8 which is sandwiched between high-resistivity blocks. The zones up to 6000m depth between D4_6 and D4_8 stations have moderate resistivity which may be an indicator of sedimentary formations.

Figure19 shows a resistivity cross section based on 2D NLCG inversion along prole D5. No major geophysical anomaly is found. The resistivity distribution is also random, which rules out any possibility of signicant reservoir presence at any depth.

Sircar et al. Geotherm Energy (2015) 3:22

Page 19 of 25

2D geoelectric shallow models

2D geoelectric shallow models (up to 1500 m depth) show the resistivity distribution (Figs. 20, 21 and 22). Low-resistivity zones are clearly visible. Shallow models can be divided into two parts, namely the upper part and the lower part. The Upper part up to 400m is considered as low-resistivity zone. Resistivity value at the upper layer is below 15ohmm in each prole shown in the Figs.20, 21 and 22. This also indicates that there is a possible groundwater zone. The lower part of each prole shows resistivity value below 60ohmm. It indicates that there is a conductive anomaly beneath the surface. Here the only dierence is in the thickness of upper layer at each prole. Therefore, the thickness of the upper layer is variable in the area of research (Srinivas etal. 2012).

Sircar et al. Geotherm Energy (2015) 3:22

Page 20 of 25

Minimum thickness is in the NW part of the area and is of approximately 150m at prole D1. The thickness increases towards NE to about 750m. Thickness of this layer also increases towards south.

Thickness of low-resistivity layer at prole D3 increases from 500m in the WSW to 1150m in the ENE. This very low-resistivity layer can be linked to ocean sediments inltrated by salt water. The thickness of low-resistivity layer at prole D2 increases from 400m in the WSW to 1000m in the ENE. Surface manifestation above the high is a hot water spring. A deep well may conrm the interpretation and provide the geothermal gradient of the subsurface process condition.

Sircar et al. Geotherm Energy (2015) 3:22

Page 21 of 25

Sircar et al. Geotherm Energy (2015) 3:22

Resistivity sections

Figures23 and 24 depict resistivity distributions at depths of 3 and 4km below sea level, respectively. Low-resistivity zones are shown as shaded area in the resistivity sections. The resistivity sections have been prepared based on 2D data. As can be seen in the gures resistivity ranges from 20 to 87ohmm in the zone of interest. The resistivity closures can be seen in both resistivity sections, which is a true validation of the model. The area covered by these closures is about 30km2. These low-resistivity areas might be an indication of a potential geothermal reservoir.

Conclusion

At Dholera geothermal site, extensive exploration activities have been carried out to delineate the geothermal prospect. The Landsat imagery study was carried out to identify the target zones based on Low Vegetation Index and high land surface temperature. Geochemical analysis of hot water from Dholera spring was carried out to determine the physical and chemical properties of water which is helpful to identify nature of the geothermal reservoir.

As a part of geophysical exploration, 2D MT and AMT surveys were carried out along six proles at Dholera. Based on the values of Skew and Tipper, it is concluded that the structure is 1D and 2D in nature. The polar diagram shows that the structure is 3D in nature at some places. Both shallow and deep geoelectric maps depict that the reservoir

Page 22 of 25

Sircar et al. Geotherm Energy (2015) 3:22

Page 23 of 25

is shale/sandstone body sandwiched between high-resistivity basalts. Based on 2D data, cross-sections were prepared at shallow and deep levels.

In conjunction to MT survey, gravity data were collected along the MT proles with oset at some stations. Corrections were applied and forward modeling of residual Bouguer gravity was done. The gravity-derived subsurface picture shows low-density zones sandwiched between high-density zones. In the representative MT cross sections presented in this paper, low-resistivity geophysical anomalies exist between high-resistivity basaltic zones. Integrating both Gravity and MT interpretation, the model supports that beneath the surface manifestations (hot springs), low-density and low-resistivity geophysical anomaly is present which might be an indication of the existence at depth of a geothermal reservoir.

Resistivity closures around the hot springs are seen in both shallow and deep cross sections which is a true validation of the model postulated. Resistivity sections prepared at 3 and 4km depth support the idea. These results suggest that Dholera is a promising site from geothermal point of view. The results will help in harnessing the geothermal energy for domestic and commercial uses. Drilling of wells and data recorded with the help of temperature log will help in estimating geothermal gradient inside the subsurface.

Abbreviations

DMIDC: Delhi Mumbai Industrial Corridor Development Corporation Limited; AMT: audio magnetotelluric; MT: magnetotelluric; RS: remote sensing; MODIS: Moderate-resolution imaging spectro-radiometer; LST: land surface temperature; NDVI: Normalize Dierence Vegetation Index; TDS: total dissolve solids; UC: upward continuation; NLCG: non-linear conjugate gradient.

Sircar et al. Geotherm Energy (2015) 3:22

Authors contributions

AS is the head of Centre of Excellence for Geothermal Energy (CEGE), which is mandated for creation of 1 MW electrical energy using geothermal energy. MS and SS are faculties at School of Petroleum Technology, PDPU and are also working as research scientists at CEGE. DV, SD and AC are employees of CEGE and are involved in R&D activities carried out in CEGE. AC also helped in drafting the gures which are used in this manuscript. All authors read and approved the nal manuscript.

Author details

1 School of Petroleum Technology, Pandit Deendayal Petroleum University, Gandhinagar, Gujarat 382007, India. 2 Centre of Excellence for Geothermal Energy, Pandit Deendayal Petroleum University, Gandhinagar, Gujarat 382007, India.

Acknowledgements

Authors acknowledge the technical support provided by M. S University, Vadodara, Gujarat, India for assisting in conducting the Geochemical Analysis at Dholera. The authors are thankful to Mr. Gaurav Negi and acknowledge the technical support provided by him. Authors are also thankful to PBG Geophysical Exploration Ltd., Poland for conducting the MT Survey at Dholera site and providing technical support. Authors acknowledge the support provided by School of Petroleum Technology, Pandit Deendayal Petroleum University, Gandhinagar, Gujarat, India. Authors are thankful to the school for giving permission to publish this research.

The authors are planning to put up an Organic Rankine Cycle (ORC) in the identied site in the next 2 years for electricity generation.

Competing interests

The authors declare that they have no competing interests.

Received: 2 July 2015 Accepted: 5 November 2015

References

Bahr K, Simpson F. Practical magnetotellurics. Cambridge: Cambridge University Press; 2005.

Biswas SK. Regional tectonic framework, structure and evolution of the western marginal basins of India. Tectonophysics.

1987;135:30727.

Biswas SK. Structure of the western continental margin of India and related activity in Deccan ood basalt. Memoir Geol

Soc India. 1988;10:37190.

Cagniard L. Basic theory of the Magnetotelluric method of geophysical prospecting. Geophysics. 1953;18:60535.

Calvin WM, Coolbaugh M, Kratt C, Vaughan RG. Application of remote sensing technology to geothermal exploration. In:

Rhoden HN, Steininger RC, Vikre PG, editors. Geological society of Nevada symposium 2005: window to the world, Reno, Nevada. 2005. p. 108389.

Constable SC, Parker RL, Constable CG. A practical algorithm for generating smooth models from electromagnetic sounding data. Geophysics. 1987;52:289300.

Cooper SM, Tianyou L, Mbue IN. The Empirical Mode Decomposition (EDM), a new tool for potential eld separation. J

Am Sci. 2010;6(7):1837.

Gupta H, Roy S. Geothermal energyan alternative resource for the 21st Century. UK: Elsevier; 2007. p. 1271. Huang S, Pollack HN. Late quaternary temperature changes seen in worldwide continental heat ow measurement.

Geophys Res Lett. 1997;24(15):194750.

Keller GV, Frischknecht FC. Electrical methods in geophysical prospecting. New York: Pergamon; 1966. p. 517.

Kulkarni VN. Geology of Gujarat senior Geologist Engineering Research Institute. Gujarat State: P.W.D; 1985.

Mariita N. The gravity method. Presented at Short Course II on surface exploration for geothermal resources, Kenya 2007. p. 19.

Merh S. Geology of Gujarat. Banglore: Geological Society of India;1995. p. 224.

Mickus KL, Aiken CLV, Kennedy WD. Regional-residual gravity anomaly separation using the minimum-curvature technique. Geophysics. 1991;56:27983.

Pasvanoglu S. Geochemical study of the geyser geothermal eld in Haukadalur, S-Iceland. Iceland: Geothermal Training

Programme; 1998. p. 281318.

PBG Report. 2D Magnetotelluric survey for four locations for geothermal exploration namely (No. of MT soundings) Unai (66), Gandhar (66) and Dholera (66) (Excluding Tulsishyam, Tuwa and Chabsar). Poland: Exploration Geophysics Service; 2014. p. 124. (CEGE Report).

Rivas J. Gravity and magnetic methods. Short Course on Surface Exploration of Geothermal Resources;2009. p. 113. Rodi W, Mackie RL. Nonlinear conjugate gradients algorithm for 2-D magnetotelluric inversion. Geophysics.

2001;66:17487.

Sahajpal S, Sircar A, Singh A, Vaidya D, Shah M, Dhale S. Geothermal exploration in Gujarat: case study from Unai. Int J

Latest Technol Eng Manag Appl Sci. 2015;4(5):3847.

Shah M, Sircar A, Vaidya D, Sahajpal S, Chaudhary A, Dhale S. Overview of geothermal surface exploration methods. Int J

Adv Res Innov Ideas Educ. 2015;1(4):5564.

Sharma N. Physico-chemical characterization of the thermal spring waters occurring in Gujarat, India;2013. p. 116. (CEGE

Report).

Srinivas Y, Raj SA, Hudson OD, Muthuraj D, Chandrasekar N. Geoelectrical inversion and evaluation of lithology based on optimized adaptive neuro fuzzy inference system (ANFIS). Stud Geophys Geod. 2012;57:115.

Page 24 of 25

Sircar et al. Geotherm Energy (2015) 3:22

Srivastava PK, Gupta DK. Study of geothermal prospects in Gujarat through remote sensing approach;2013. p. 113

(CEGE Report).

Talwani M, Worzel JL, Landisman M. Rapid gravity computations for two-dimensional bodies with applications to the

Mendocino submarine fracture zone. J Geophys Res. 1959;64:4959.

Vaidya D, Shah M, Sircar A, Sahajpal S, Dhale S. Geothermal energy: exploration eorts in India. Int J Latest Res Sci Technol. 2015;4(4):123.s.

Vozo K. The magnetotelluric method. Sydney: Centre for Geophysical Exploration Research Macquarie University; 1991. p. 641709.

Page 25 of 25