Nonlinear Processes in Geophysics Predicting heat flow in the 2001 Bhuj earthquake ( M w = 7 . 7 ) region of Kachchh ( Western India ) , using an inverse recurrence method

Abstract. Terrestrial heat flow is considered an important parameter in studying the regional geotectonic and geodynamic evolutionary history of any region. However, its distribution is still very uneven. There is hardly any information available for many geodynamically important areas. In the present study, we provide a methodology to predict the surface heat flow in areas, where detailed seismic information such as depth to the lithosphere-asthenosphere boundary (LAB) and crustal structure is known. The tool was first tested in several geotectonic blocks around the world and then used to predict the surface heat flow for the 2001 Bhuj earthquake region of Kachchh, India, which has been seismically active since historical times and where aftershock activity is still continuing nine years after the 2001 main event. Surface heat flow for this region is estimated to be about 61.3 mW m−2. Beneath this region, heat flow input from the mantle as well as the temperatures at the Moho are quite high at around 44 mW m−2 and 630 °C, respectively, possibly due to thermal restructuring of the underlying crust and mantle lithosphere. In absence of conventional data, the proposed tool may be used to estimate a first order heat flow in continental regions for geotectonic studies, as it is also unaffected by the subsurface climatic perturbations that percolate even up to 2000 m depth.


Introduction
Terrestrial heat flow is recognized as an important parameter for understanding the thermal state of the earth's crust and and the upper mantal (Cermak et al. 1991;Nyblade and Pollack, 1993;Pollack et al., 1993;Rudnick et al., 1998;Pandey and Agrawal, 1999;Jaupart and Mareschal, 1999).
Correspondence to: N. Vedanti (nimisha@ngri.res.in)It is closely related to the deeper temperature regime and thus has a direct relevance to regional geotectonic studies and underlying evolutionary processes (Pandey, 1991;Pollack et al., 1993).Despite a large number of heat flow measurements in the past, their distribution is still quite uneven, especially in continental regions which are dominated by North America and Europe (Sclater et al., 1980;Pollack et al., 1993;Neumann et al., 2000;Davies and Davies, 2010).Well over half the earth's land area is still unsampled.Its measurement essentially involves knowledge of the temperature gradient, preferably in a borehole and thermal conductivity of the in-situ rocks pierced by the same borehole (Gupta and Rao, 1970;Jessop, 1990;Pandey, 1981).A large number of boreholes drilled every year cannot be used either because of regulations requiring them to be sealed or because of financial considerations, under which casing is pulled out and the borehole collapses often before the thermal equilibrium is re-established, which is a pre-requisite and requires a long time.In spite of some sporadic attempts, high cost of drilling inhibits heat flow measurements in desired areas.
The most important areas lacking heat flow information are primarily those which are of low economic value.Such areas are often of great geodynamic importance, for example areas associated with active seismicity like Bhuj area of the Kachchh (Gujarat), which is a unique site of intraplate seismic activity.This zone is located in western part of India (23.3 • -23.7 • N; 69.9 • -70.2 • E) (Fig. 1).It has been seismically active since historical times, having experienced two large magnitude earthquakes, one in 1819 (Kachchh earthquake, M w = 7.7) and the other in 2001 (Bhuj earthquake, M w = 7.7) (Rajendran et al., 2008).The latter event caused widespread damage killing over 20 000 people, with aftershock activity still continuing nine years after the main shock (Mandal et al., 2007).This seismic region is especially known for lower crustal seismic activity (Mandal and Pandey, 2010).In the past, this region has been subjected to several tectonic, thermal and magmatic phases, which are manifested through a complex and heterogeneous crustal velocity structure (Biswas, 1987;Gambose, et al., 1995;Kayal et al., 2002;Mandal and Pujol, 2006;Karmalkar et al., 2008;Mandal and Pandey, 2011).In spite of immense importance of this region to seismotectonic studies, no heat flow data is yet available, thereby inhibiting modeling of underlying rheological and thermal regime.
A need was therefore felt to find a suitable method through which a reasonable first order heat flow estimate can be made in such areas.In the past, it has been done using indirect methods like empirical relationship between heat flow and (i) lithospheric thickness (Negi et al., 1987;Chapman and Pollack, 1977), (ii) Pn velocity and seismic models of the crust and upper most mantle (Kubick,1988, Shapiro andRitzwoller, 2004), and (iii) regional geology (Chapman and Pollack, 1975;Sclater et al., 1980;Pollack et al., 1993;Davies and Davies, 2010).In the present paper, we first develop a methodology in the form of an inverse recurrence method to predict a first order surface heat flow for any particular region.We test this methodology in several geotectonic provinces around the globe and then using this approach, we finally predict the expected heat flow for the Kachchh earthquake region of Gujarat (western India).This method requires the knowledge of depth to the lithosphereasthenosphere boundary (LAB), seismic structure of the underlying crust, and the radioactivity of exposed crystalline rocks.
Due to advent of broadband seismic and receiver function techniques, the depth to the LAB is now known with sufficient accuracy in several continental regions.It is considered a first-ordered structural discontinuity that differentiates overlying cool and more rigid lithosphere from Fig. 2. Distribution of mantle solidus with depth (Modified from Gass et al., 1978).Lithosphere-asthenosphere boundary (LAB) has been defined as the intersection point of the mantle solidus and the estimated continental geotherm.the underlying hotter asthenosphere, which flows more easily.Artemieva and Mooni (2001) have adopted the thermal thickness of the lithosphere to be defined as the intersection of geotherm with a mantle adiabat of 1300 • C. It is now understood that the solid state creep and convective heat transfer starts dominating over conduction at around 1100 • C only (Pollack and Chapman, 1977).Therefore, we adopt the base of the lithosphere as being the mantle solidus, which is the intersection of the continental geo-isotherm with the peridotite incipient melting point curve (Gass et al., 1978) (Fig. 2).This is characterized by a temperature of about 1100 • C at surface to about 1360 • C at the depth of about 200 km.Global geotherms usually intersect the mantle solidus at a depth coincident with the top of the seismic low velocity zone, i.e. the base of the lithosphere (Chapman and Pollack, 1977).Similarly, seismic structure of the crust is also now fairly well known in many areas through deep seismic profiling and receiver function analysis.These studies provide detailed velocity discontinuities, which give accurate information on the nature, composition, and thickness of different crustal and uppermost mantle layers.The only other parameter, radioactivity of the upper crustal rocks, can be easily deciphered by measuring U-Th-K concentrations of the crystalline basement rocks in the lab with the help of geochemical instruments like Inductively Coupled Plasma Mass Spectrometer (ICP-MS) and X-Ray Fluorescence (XRF) spectrometer.The measured accuracy and precision of results in cases of ICP-MS are well within 20 % RSD of the certified value (Balaram and Rao, 2003); and in XRF, it is between 1 and 2 % (Krishna and Govil, 2007).Measured U-Th and K concentrations can then be converted into heat generation using following equation (Birch, 1954): HG(µW m −3 ) = ρ(0.035K + 0.097 U + 0.026 Th) where HG = heat generation, ρ = density of the rock (g cm −3 ), K concentration in (%) and U, Th concentrations in ppm.
Heat generation is quite high in upper crustal rocks consisting of granites and gneisses, but comparatively it is far lower in the middle and lower crust and almost negligible (around 0.01 µW m −3 ) in the upper most mantle.

Methodology
We know that heat flow of any layer (say n) can be given as (Lachenbruch, 1970) Where q n+1 is heat flow of layer (n+1) t n is the thickness of layer (n) A n is radiogenic heat generation of layer (n) If we consider a lithospheric model, then the lithosphereasthenosphere boundary (LAB) will be the last interface in the model, thus in this case q n+1 = q LAB Now for layer (n) just above the LAB, we can write that And then for the next upper layer (n-1) Or Thus, the heat flow of any layer can be written in terms of q LAB .Writing these expressions recursively for next upper layers, until the top layer, we get the surface heat flow as Where "n" is no. of layers and q s is the surface heat flow.This implies that if we know the q LAB , radiogenic heat generation and thickness of each layer we can compute the surface heat flow and heat flow of every intermediate layer.Now to compute q LAB we start with the standard equation of temperature computation using constant heat generation model which is given as (Jaeger, 1965) Where T n+1 is temperature of layer (n+1) T n is temperature of layer (n) q n is heat flow of layer (n) t n is the thickness of layer (n) k n is the thermal conductivity of layer (n) A n is radiogenic heat generation of layer (n) For a given "n" layered lithospheric model, we know the LAB, thus, in Eq. ( 4) we know the temperature at LAB i.e.T LaB which is nothing but the temperature of the last layer.
T n+1 = T LAB Thus we can write as It is mentioned earlier that the heat flow of any layer can be written in terms of q LAB (Eq.2).Thus in Eq. ( 5) we have two unknowns T n and q LAB We rewrite Eq. ( 5) as Or substituting the value of T n in terms of T n−1 using Eq.(4), we get Or Or we can rewrite Eq. (8) as Or Thus q LAB for any layered model can be computed by using Eq. ( 10) Substituting Eq. ( 10) in Eq. ( 3) we get equation to compute the surface heat flow q S for any region where LAB is known Substituting Eq. ( 10) in Eq. ( 2) heat flow of any layer can be computed.

Application to known heat flow locations
Broadband seismic information, detailed crust-mantle structure, and the depth to the LAB is now known in many geological provinces of the world.To test the proposed methodology, we have chosen eleven separate locations worldwide, where all the required inputs (viz., measured surface heat flow, crustal seismic structure, and depth to the LAB) are known.Then, considering the effectiveness of the proposed methodology, we moved on to predict the surface heat flow for the Bhuj earthquake region of Kachchh (Gujarat, India) (Fig. 1) as a case study.

New Madrid Seismic Zone (USA)
This region forms the largest and most seismically active region of the eastern USA, having witnessed three major earthquakes (M > 7.0) during 1811-1812.It is associated with an ancient intraplate rift zone, commonly known as Reelfoot rift.On the basis of regional network data and seismic reflection surveys, Mooney et al. (1983) proposed a five-layer crustal model for the region, which consists of a 14 km thick altered lower crust.Details of crustal layers, together with the relevant geothermal parameters (Mooney et al., 1983;Liu and Zoback, 1997;Bedle and Lee, 2006), are given in Table 1.In this zone, measured heat flow is about 60 mW m −2 and depth to the LAB is about 72 km, as obtained below the Reelfoot rift through tomographic study (Bedle and Lee, 2006).

New Hampshire (USA)
This region has been extensively studied as far as surface heat flow and crustal radioactivity is concerned (Jaupart et al., 1982).Based on the migrated Ps waveform images, a sharp and very well defined LAB has been mapped at the depth of 95-100 km below the HRV seismic station (42.51 • N; 71.56 • W).Moho discontinuity is reported to be at 30 km (Rychert et al., 2005).At this particular station, seismic phases have been very clearly observed.Further, HRV is situated at the centre of a quadrant bounded by 42  (Jaupart et al., 1982).Heat production of the upper crustal rocks is about 2.16 µW m −3 .Crustal and upper mantle velocity structure below HRV is also known from the joint inversion of receiver function and surface-wave dispersion of Rayleigh group and phase velocity (Fnais, 2004).

Basin and Range Province (USA)
Basin and Range Province forms one of the most studied geotectonic units of the western United States.It is characterized by anomalous thermal and crustal structure.Here, the measured heat flow is quite high (average ∼ 80 mW m −2 ) (Blackwell, 1971).The LAB lies at an extremely shallow depth of about 60 km (Li et al., 2007).Crustal structure of this region is well defined (viz., Catchings and Mooney, 1991).Like LAB, Moho in this region is also elevated to about only 30 km.Estimated heat production for the upper crust is 2.23 µ W m −3 (Blackwell, 1971).

Colorado Plateau (USA)
This plateau forms a large crustal block in the SW United States that has been raised to nearly 2 km above msl due to underlying geological forces.Seismic structure of this region is well studied (Keller et al., 1976;Parsons et al., 1996;Gilbert and Sheehan, 2006;Roller, 1965).Moho below this region is about 40 km deep and depth to the LAB is 105 km (Biswas and Knopoff, 1974;Chapman and Pollack, 1977).Measured heat generation in the upper crust is 3.07 µW m −3 (Decker, 1969).Surface heat flow in this region is 59-64 mW m −2 (Reiter and Clarkson, 1983).

Kaapvaal craton (South Africa)
This craton, situated in South Africa, is a very large craton covering an area of about 1.2 million sq.km.It is one of the oldest cratons having been formed and stabilized between 3.7 and 2.6 Ga.LBTB seismic station (25.02 • S, 25.60 • E) lies in the northern part of the cratonic shield beneath which Moho as well as LAB is well defined.Moho lies at the depth of about 42 km, while the LAB is extremely deep at 293 km depth (Kumar et al., 2007).Crustal discontinuities (Table 1) are adopted from Midzi and Otlemolter (2001), Durrheim and Green (1992), and Kumar et al. (2007).Reported heat flow at LBTB (LOBATSE) is 30 mW m −2 , with an uncertainty of ±15 % (Ballard et al.,1987).

Yilgarn craton (western Australia)
It is also a quite large Archean (3.2-2.8Ga) craton, constituting the bulk of western Australia.In this craton, 3-D seismic structure of the underlying lithosphere (Goleby et al., 2006) indicates the LAB to be at the depth of 220 km.Crustal structure for this region is also known, with Moho depth ranging between 35 and 38 km (Goleby et al., 2006).The upper crustal radioactive layer, which extends to only about 4.5 km depth, has an average heat production of 3.02 µW m −3 (Sass and Lachenbruch, 1979).Measured heat flow for this region is 37 ± 6 mW m −2 (Jaupart and Mareschal, 1999).

Baltic Shield (Fennoscandia)
This shield is located mainly in Fennoscandia (Norway, Sweden, and Finland).It forms an integral part of the east European craton containing Archean-Proterozoic rocks.Crustal and upper mantle structure beneath this shield is known (Lund and Slunga, 1981).Estimated depth of the LAB is 200 km in the northern part of this shield (Olsson et al., 2007), where measured heat flow and upper crustal heat generation is 45-55 mW m −2 and 3.5 µW m −3 , respectively (Pasquale et al., 1991).

Bohemian Massif (Central Europe)
Bohemian Massif is one of the largest stable outcrop of Variscan rocks in central Europe, located on the territory of Czech Republic, Germany, Poland, and Austria.This region has also been well studied seismically and geothermally.
Crustal discontinuities are shown in Table 1.Moho lies at the depth of about 39 km (Hrubcova et al., 2005), while the LAB is at 111 km depth (Geissler et al., 2010).Measured average surface heatflow and heat generation are, respectively, at 86 mW m −2 and 4.64 µW m −3 (Cermak et al., 1991).More than half of the observed surface heat flow over this region is contributed by upper parts of the crust.

Killari (Latur) earthquake region (India)
This earthquake region, situated in Latur district of Maharashtra, has been associated with one of the deadliest intraplate earthquakes (M w = 6.2) that killed about 10 000 people in 1993.This event led to a large number of geoscientific investigations, including a specially drilled 617 m deep borehole (KLR-1) to understand the seismotectonics of this region.A broadband seismic station was subsequently set up in the earthquake affected region in order to delineate crustal shear velocity structure.Teleseismic receiver function and surface wave phase velocity data are used to constrain the shear wave velocity, which was found to be extremely high for the entire crustal column (Rai et al., 2003).
From about 2 km depth downwards, the Vs hoovers between 3.8 and 4.0 km s −1 .Moho was detected at 37 ± 2 km depth.Measured heat flow in KLR-1 borehole is about 43 mW m −2 (Roy and Rao, 1999).We have also measured the heat generation in the vertical column of the crystalline basement encountered in KLR-1 borehole, between the depths 350 and 620 m, which represent mid crustal amphibolite-granulite facies rocks (Pandey et al., 2009).Average heat generation was found to be about 0.78 µW m −3 .At this location, depth to the LAB is about 110 km (Priestley et al., 2006).

Tummalapalli (Cuddapah basin, India)
This late Archean-early Proterozoic intracratonic region is situated in the eastern part of the Dharwar craton (southern India).At this location, measured heat flow is 51 mW m −2 , while crustal heat production of the gneissic basement is 2.20 µW m −3 (Roy and Rao, 2000).Adapted crustal structure is based on detailed deep seismic sounding study (Mall et al., 2008).Depth to the LAB beneath this locality is 80 km (Kumar et al., 2007).

Adopted thermal parameters
For all the localities discussed above, adopted radioactivity and thermal conductivity of the mid to lower crust and the mantle lithosphere is given in Table 2. Heat generation and thermal conductivity for sediments as well as of underplated magmatic crust are taken from Liu and Zoback (1997).
Mid crustal heat generation is an important parameter in elucidation of thermal regime inside the earth, which cannot be measured directly.However, in order to understand the genesis of 1993 Killari earthquake (M w = 6.3), a 617 m deep borehole (KLR-1) was drilled in the epicentral region that penetrated 280 m of late Archean mid crustal amphibolite -granulite facies basement beneath 338 m of Deccan volcanic rocks (Pandey et al., 2009).An attempt was made to measure the radio elemental concentration of the basement samples in the laboratory by ICP-MS and XRF instrument.Based on these measurements, we got a mean heat generation of 0.78 µW m −3 for the basement, which we adopt uniformly.Similarly, heat generation for the granulitic lower crust is taken from a detailed study by Ray et al. (2003), who found a range of 0.14 to 0.20 µW m −3 , with an average of 0.16 µW m −3 .This value is quite close to the reported estimate of 0.18 µW m −3 for the lower crust by Rudnick and Fountain (1995).Similarly, the adopted radioactivity and thermal conductivity for ultramafic mantle (0.01 µW m −3 and 3.0 W/m • C, respectively) corresponds to periodotitic composition.Further, as mentioned earlier, the LAB temperature is taken as mantle solidus temperature (i.e, the intersection of the geotherm with the peridotitic incipient melting point curve (Gass et al., 1978;Fig. 2).
Parameters given in Tables 1 and 2 are taken as input to the developed methodology to compute surface heat flow.In Table 3, we have included measured as well as predicted heat flow by the proposed inverse recurrence method for all the locations discussed in the text.It is obvious that the heat flow estimated from the proposed method conforms quite well with the measured heat flow (Table 3).Therefore, we use this methodology to estimate the surface heat flow in the Bhuj earthquake region of Kachchh (western India), which will be useful to our understanding of intraplate lower crustal seismogenic process.

Predicting heat flow for Kachchh earthquake region
(western India): a case study

Seismotectonics of the Kachchh region
The Kachchh region is part of a network of rifts that developed along the western margin of India following the breakup of Gondwanaland (Biswas, 1987; Fig. 1).The Quaternary/Tertiary sediments, Deccan volcanic rocks, and Jurassic sandstones resting on Archean basement characterise the geological sequence of the Kachchh region (Gupta et al., 2001).Major structural features of the Kachchh region include several E-W trending faults/folds as shown in Fig. 1.The rift zone is bounded by a north dipping Nagar Parkar fault in the north and a south dipping Kathiawar fault in the south.Other major faults in the region are the E-W trending Allah Bund fault, Island belt fault, Kachchh mainland fault, and Katrol Hill fault.In addition, several NE and NW trending small faults/lineaments are also observed (Biswas, 1987).Latest paleoseismological data substantiate that this region has been experiencing large earthquakes of M > 7.0 since 325 BC (Rajendran et al., 2008).
Recently, the joint inversion of P-receiver function and surface wave group velocity dispersion data depicts a variation in the thickness of crust (35-43 km) and lithosphere (62-77 km) below the epicentral zone of 2001 Bhuj earthquake (Mandal and Pandey, 2011).

Seismic network
The

Seismic structure
Detailed crustal structure of this region is now available, based on a specially designed seismic network that recorded about 5000 aftershocks of M w ≥ 2.0 between 2001 and 2008.
Out of this, 13 862 P-and 13 766 S-wave high quality arrival times were picked up from the seismograms of 2303 well located aftershocks (Mandal and Chadha, 2008).To this, the data collected by the digital networks of CERI (USA) and Hirosaki University (Japan) were also added, thereby covering as many as 58 seismic stations (Fig. 1), which resulted in an excellent sampling of the crustal volume below the aftershock zone (Mandal and Chadha, 2008).This data was used to determine three-dimensional P-and S-wave velocity structure by constructing Vp and Vs tomograms for the different depth bands (0-2, 2-6, 6-8, 8-18, 18-20, 20-26, 26-34, and 34-42 km) (Mandal and Pandey, 2011).An average crustal model constructed from such tomograms over an area of 23.3 • to 23.7 • N and 69.9 • to 70.2 • E is shown in Fig. 3, which indicates an extremely complex crust mantle structure (Mandal and Pandey, 2010) below the earthquake affected region, somewhat similar to New Madrid Seismic Zone (Liu and Zoback, 1997).Thermal parameters and crustal structure used in heat flow estimation are included in Tables 1 and 2.

Lithosphere -Asthenosphere boundary (LAB)
Receiver function from the sites situated in the rift zone, where the 2001 Bhuj earthquake occurred, show a strong negative phase at 8-10s, which corresponds to LAB (Fig. 4).This has been confirmed by the joint inversion of the stacked radial receiver functions and surface wave group velocity dispersion data, which has been modeled with a velocity decrease at a sub-crustal depth of about 63-77 km in an area of 126 × 80 km following Julia et al. (2000).Such an example for one of the seismic stations 'VJP' is shown in Fig. 4 (Mandal and Pandey, 2011).In general, the depth to the LAB beneath the fifteen broadband seismic stations in the Kachchh region are found to be quite shallow, with an average of about 70 km (Mandal and Pandey, 2011).Corresponding to this depth, mantle solidus temperature would be ∼1150 • C (Fig. 2).Moho at this location is 34 km deep (Figs. 3 and 4).

Heat flow estimate
Using the present approach, we arrive at the surface heat flow of 61.3 mW m −2 for the central part of the Kachchh seismic zone.No heat flow measurement is yet available for this region.However, we do have some indirect information.For example, the Kachchh region is known to contain mantle xenoliths in the alkali rocks (Karmalkar et al., 2000).
Based on mantle derived paleo-geotherms on these xenoliths, a minimum heat flow between 60 and 70 mW m −2 has been projected for this region by Mukherjee and Biswas (1988), following the method of calculation from Pollack and Chapman (1977).Similarly, empirical relationships given by Negi et al. (1987) and Chapman and Pollack (1977) would also result in a surface heat flow of about 60 to 75 mW m −2 against a lithospheric thickness of 70 km, which is quite accurately known (Fig. 4).(Julia et al., 2000).Correlation between the receiver functions obtained from observed data and joint inversion are shown by solid grey and dotted black lines, where, "a" and "p" represent the Gaussian width factor and horizontal slowness, respectively, which are used to compute receiver functions; (b) correlation between group velocity dispersion curves obtained from observed data (solid grey line with solid square symbols) and joint inversion (dotted black line); and (c) final model obtained from joint inversion.M represents the Moho and L represents the LAB.
Our calculation would also indicate an extremely high input of heat flow from the mantle (about 44 mW m −2 ) and the presence of high temperatures (∼630 • C) at the Moho depth of about 34 km, below the earthquake affected zone.This may apparently be related to thermal restructuring of the crust and mantle lithosphere and influx of large amount of heat, aqueous gasesous fluids, and ultramafic melts from the mantle during the course of lithospheric stretching and consequent subcrustal erosion due to rise of isotherms.Measured heat flow averages around 77 mW m −2 in the adjoining areas like Cambay graben (Gupta, 1981).

Discussion of results
Using the proposed inverse recurrence method, we arrive at a surface heat flow of 61.3 mW m −2 for the central part of the Kachchh seismogeneic region.Present study thus demonstrates (Table 3) that it is possible to estimate a first order surface heat flow for any continental region using this method, provided the detailed crustal structure and the depth to the LAB are known with sufficient accuracy through seismic studies.The only other crucial parameter needed is the radioactive heat generation of the upper crustal rocks, which can be easily deciphered by measuring U, Th, and K by XRF and ICP-MS instruments, accessible in most of the geochemical labs as explained earlier.Other parameters can be www.nonlin-processes-geophys.net/18/611/2011/ Nonlin.Processes Geophys., 18, 611-625, 2011 adopted from Table 2. Except for the upper crustal heat generation, all other parameters usually vary in a narrow range from place to place.Depending on the in-situ velocity distribution, the crust can be subdivided into different geologic segments (as defined in Table 2), wherein the geothermal parameters can be assigned accordingly.Heat flow estimated using the proposed method appears to have an added advantage over the conventional measurements, as majority of them come from shallow (<1000 m) boreholes (viz.Jessop et al., 1976).Recent study of Gosnold et al. (2005), which is based on the analysis of 1500 deep boreholes located throughout the Fennoscandian shield, east European platform, Danish basin, Germany, Czech Republic, and Poland, indicate an increase in heat flow with depth.This phenomena has also been seen in North American data (Gosnold et al., 2005).Similar inferences have been obtained in other studies too (Pandey, 1981;Safanda and Rajver, 2001;Kukonen and Joeleht, 2003;Majorowic et al., 2008).It would mean that the heat flow may have been underestimated for many areas.This may be especially true for mid to higher latitudes where warming of the ground surface since the Pleistocene has resulted in a transient reduction in surface temperature gradient which often persists to a depth of up to about 2 km.Such signal is also clearly visible in KTB borehole (Germany) (Clauser et al., 1997), which is one of the deepest drilled ever.In this borehole, measured present day heat flow is systematically too low, down to ∼1500 m depth (Clauser et al., 1997).It is common knowledge that in shallow boreholes, more recent climatic perturbations upto the Holocene period are important (Beck, 1977).Thus, heat flow estimated using the proposed method, in principle, should be devoid of subsurface climatic perturbations.
Further, observed surface heat flow has often been correlated with the thickness of the lithosphere (Chapman and Pollack, 1977;Crough and Thompson, 1976;Negi et. al, 1987).However, such plots have shown lots of scatter, primarily due to the fact that on an average almost 40 % of the observed surface heat flow is contributed by the crust itself (Pollack and Chapman, 1977), especially the outer half which is highly enriched in radioactive elemental concentrations compared to lower half and upper most mantle.From many segments of the earth, the outer part of the earth's crust has been considerably eroded, thus the obtained heat flow at the surface does not truly represent deeper thermal regime.Instead, mantle heat flow estimate can be a better guide to the physical and tectonothermal conditions underneath.It can be deduced by subtracting the heat flow arising in the earth's crust from the observed heat flow at the surface.Figure 5 shows the plot between mantle heat flow and the depth to the LAB, which shows a quite good fit between the two parameters.
Further, we all know that the crustal and lithospheric mantle layers are lithologically and compositionally heterogeneous (Rudnik and Fountain, 1995).That makes it difficult to assign the input values for temperature-depth calculations.This includes characterization of radioactive heat generation and thermal conductivity data for different crustal layers.Calculated temperatures do vary with the choice of input data values (Jokinen and Kukonen, 1999).For the present study, we have chosen only such inputs which are usually adopted in such calculations by the global researchers.Still, to have an idea of possible errors in estimated heatflow, we have varied the radioactive heat generation of the layers 2, 3, and 4 by 5 to 20 % in the case of Bhuj earthquake region (Table 1).We found a variation of ±2 mW m −2 in cases of 20 % changes in adopted heat generation values.Thus, the likely uncertainities in adopted values may not have much effect on the end results.
We would like to emphasize here that the purpose of this paper is to provide a simple, usable method for estimating reasonable first order heat flow at the surface for regional geotectonic studies and not to show the superiority of this method over the conventional heat flow measurements.One should always prefer conventional heat flow value if it is derived from a sufficiently deep, cored, and thermally stabilized borehole; however, this may not always be the case.
aftershock activity of the 2001 Bhuj earthquake has been monitored by National Geophysical Research Institute, Hyderabad, during 2001-2002 with a local digital network consisting of eight 24-bit recorders with an external hard disk (2 GB) and GPS timing system.Out of these, six were equipped with short period seismometers (frequency range 1-40 Hz) and two were equipped with broadband seismometers (natural period 30-100 s).Some station locations were shifted during 2001-2002, making a total of 12 station sites.Further, we used data from seismic network of CERI, Memphis (eight stations equipped with Kinemetrics K2 digital recorders and velocity sensor, from 13 February 2001 to 27 February 2001), and Hirosaki University, Japan (nine short-period seismograph stations, from 28 February 2001 to 7 March 2001).This allows us to have a total of 28 stations data (at a time maximum 18 stations).During 2002-2005, NGRI monitored the aftershock activity with a network of five broadband seismograph stations.Again in 2006, NGRI deployed another eight seismographs equipped with broadband sensors (natural period 120 s).Thus, NGRI has been monitoring this activity with a network of 12 broadband seismographs since 2006.The distance between station and epicenters varies from 14 to 90 km.Recording was done in a continuous mode at 100 sps.The above-discussed deployment of stations enables us to use a total of 58 stations data (at a time maximum 18 stations) for our 3-D velocity tomography, which resulted in an excellent sampling of the crustal volume below the aftershock zone.

Fig. 4 .
Fig. 4. (a)Estimation of crustal and lithospheric thickness at VJP seismic station using joint inversion of radial receiver functions and surface wave group velocity dispersion data(Julia et al., 2000).Correlation between the receiver functions obtained from observed data and joint inversion are shown by solid grey and dotted black lines, where, "a" and "p" represent the Gaussian width factor and horizontal slowness, respectively, which are used to compute receiver functions; (b) correlation between group velocity dispersion curves obtained from observed data (solid grey line with solid square symbols) and joint inversion (dotted black line); and (c) final model obtained from joint inversion.M represents the Moho and L represents the LAB.

Fig. 5 .
Fig. 5. Relationship between mantle heat flow and the depth to the lithosphere boundary (LAB),for different geotectonic units studied in the present study.

Table 1 .
Thickness of lithospheric layers and relevant thermal parameters of different geotectonic units used in the prediction of heat flow using inverse recurrence method.

Table 2 .
Adapted geothermal parameters for the crust and underlying mantle lithosphere in estimating surface heat flow.

Table 3 .
Measured and predicted heat flow for different geotectonic locations.Data source is given in the text.