- Email: [email protected]
Quantification of crustal geotherms along with their error bounds for seismically active regions of India: A Matlab toolbox
ARTICLE IN PRESS Computers & Geosciences 35 (2009) 2095–2099
Contents lists available at ScienceDirect
Computers & Geosciences journal homepage: www.elsevier.com/locate/cageo
Quantification of crustal geotherms along with their error bounds for seismically active regions of India: A Matlab toolbox$ Kirti Srivastava , Likhita Narain, Swaroopa Rani, V.P. Dimri National Geophysical Research Institute, Council of Scientific and Industrial Research, Hyderabad 500007, Andhra Pradesh, India
a r t i c l e in fo
abstract
Article history: Received 11 September 2008 Received in revised form 24 December 2008 Accepted 29 December 2008
Determination of the thermo-mechanical structure of the crust for seismically active regions using available geophysical and geological data is of great importance. The most important feature of the intraplate earthquakes in the Indian region are that the seismicity occurs within the entire crust. In Latur region of India an earthquake occurred in the upper crust. In such situations, quantifying the uncertainties in seismogenic depths becomes very important. The stochastic heat conduction equation has been solved for different sets of boundary conditions, an exponentially decreasing radiogenic heat generation and randomness in thermal conductivity. Closed form analytical expressions for mean and variance in the temperature depth distribution have been used and an automatic formulation has been developed in Matlab for computing and plotting the thermal structure. The Matlab toolbox presented allows us to display the controlling thermal parameters on the screen directly, and plot the subsurface thermal structure along with its error bounds. The software can be used to quantify the thermal structure for any given region and is applied here to the Latur earthquake region of India. & 2009 Elsevier Ltd. All rights reserved.
Keywords: Stochastic Gaussian Random thermal conductivity Heat conduction Matlab
1. Introduction The thermal structure of the Earth’s crust is influenced by its geothermal parameters such as thermal conductivity, radiogenic heat sources and initial and boundary conditions. Two approaches to modeling are commonly used to estimate the subsurface temperature field: (1) the deterministic approach and (2) the stochastic approach. In the deterministic approach the subsurface temperature field is obtained by assuming that the controlling thermal parameters are known with certainty. However, due to the inhomogeneous nature of the Earth’s interior there are some uncertainties in estimates of the geothermal parameters. Uncertainties in these parameters may arise from the measurement inaccuracies or a lack of information about the parameters themselves. In the stochastic approach the uncertainties in parameters are incorporated and an average picture of the thermal field along with its associated error bounds is determined. To assess the properties of the system at a glance we need to obtain the mean value that gives the average picture, along with the variance or standard deviation that indicate variability or errors associated with the system behavior due to errors in the system input.
$
Code available from server at http://www.iamg.org/CGEditor/index.htm
Corresponding author.
E-mail address: [email protected] (K. Srivastava). 0098-3004/$ - see front matter & 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.cageo.2008.12.013
Subsurface temperatures are seen to be very sensitive to perturbations in the input thermal parameters and hence several studies have been carried out to quantify perturbations in the temperature and heat flow using stochastic analytical and random simulation techniques. Uncertainty in the heat flow was quantified by Vasseur et al. (1985) using a least squares inversion technique incorporating uncertainties in the temperature and thermal conductivities. The effect of variations in heat source on the surface heat flow has also been studied (Vasseur and Singh, 1986; Nielson, 1987). In most studies the stochastic heat equation was solved using the small perturbation method. Srivastava and Singh (1998) used the small perturbation method to analytically solve the steady-state heat conduction equation, incorporating Gaussian uncertainties in the heat sources, and obtained the closed form analytical expressions for the mean temperature field and its variance. The thermal structure is obtained numerically using the random simulation method or a Monte Carlo simulation. Using a Monte Carlo simulation, several authors have studied the subsurface thermal field along with its error structure by incorporating randomness in the controlling thermal parameters such as thermal conductivity, radiogenic heat production and boundary conditions (Royer and Danis, 1988; Gallagher et al., 1997; Jokinen and Kukkonen, 1999a, b). This numerical modeling is very useful when studying nonlinear problems, but sometimes a simple 1-D analytical solution to the mean behavior and its associated error bounds is more useful for quantifying the uncertainty.
ARTICLE IN PRESS 2096
K. Srivastava et al. / Computers & Geosciences 35 (2009) 2095–2099
Stochastic differential equations in other fields are now being solved by yet another approach called the decomposition method (Adomian, 1994; Serrano, 1995). Adomian’s method of decomposition has now been generalized as a general analytic procedure to solve deterministic or stochastic and linear or nonlinear equations. It has been shown to be systematic, robust, and sometimes capable of handling large variances in the controlling parameters. In a recent study, using this new approach the stochastic heat conduction equation was solved by Srivastava and Singh (1999), incorporating uncertainties in the thermal conductivity, where the solution for the temperature field was obtained using a series expansion method. The thermal conductivity was considered to be a random parameter with a known Gaussian colored noise correlation structure. Later, Srivastava (2005) extended the study to obtain the analytical expressions for mean heat flow and its variance. Graphical user interface (GUI) packages in Matlab are becoming very popular with geoscientific researchers. Witten (2004) developed a sequence of Matlab m-files and two graphical user interfaces to display raw or processed geophysical data to produce the final graphics. In this paper, we have developed a simple graphical user interface viewer which is a MATLAB-based software consisting of m-files. This software is the first package that computes the thermal structure along with its error structure, with the important advantage that it is not a compiled code. The m-file in the package is integrated through a GUI and the controlling thermal parameters such as crustal thickness, radiogenic heat production, characteristic depth, surface temperature, surface heat flow, mean thermal conductivity, coefficient of variability in thermal conductivity and correlation length scale are all given on the screen. The program model is intended as a user friendly tool to deal with various aspects of analysis and interpretation of subsurface geothermal structure for any given region.
Variance in the temperature 2 T
s ¼ c1 Term1 þ c2 Term2 þ c3 Term3 þ c4 Term4
(5)
where CK ¼ sK/K is the coefficient of variability in thermal conductivity and c1 ¼ C 2K A20 ð1 rDÞ2 =K
2
2 c2 ¼ C 2K A0 rðrD 1ÞðQ s A0 DÞ=K
c3 ¼ c2 2 c4 ¼ C 2K r2 ðQ s A0 DÞ2 =K
Term1 ¼
1 fðrD 1Þð2z2 2zD D2 e2z=D þ D2 Þ 4ðr 1=DÞ2 4½zðr 1=DÞ þ 1 þ ½zðr þ 1=DÞ ezðrþ1=DÞ þ 1 ðr þ 1=DÞ2 þ ½2zD þ D2 e2z=D D2 g 1 þ fðrD þ 1Þð2z2 2zD D2 e2z=D þ D2 Þ 4ðr þ 1=DÞ2 4 þ ½zðr 1=DÞezðrþ1=DÞ þ e2z=D ezðrþ1=DÞ ðr 1=DÞ2 ½2zD þ D2 e2z=D D2 g
Term2 ¼
1
r2 þ
f2rðz2 D 2zD2 2D3 ez=D þ 2D3 Þ ð1 þ rzÞ ðr þ 1=DÞ2 erz ðr 1=DÞ2
½zðr þ 1=DÞ þ ezðrþ1=DÞ 1 ½zðr 1=DÞ þ ezðr1=DÞ 1g
2. Mathematical formulation
1 fðr 1=DÞðz2 D 2zD2 2D3 ez=D þ 2D3 Þ ðr 1=DÞ2 zðr 1=DÞ þ 1 ½rz þ erz 1 2
The heat conduction equation with random thermal conductivity can be expressed as d dT (1) ðK þ K 0 ðzÞÞ ¼ A0 ez=D dz dz
þ ½zD þ D2 ez=D D2 g 1 þ fðr þ 1=DÞðz2 D 2zD2 2D3 ez=D þ 2D3 Þ ðr þ 1=DÞ2
Term3 ¼
where T is the temperature (1C), A(z) the radiogenic heat source (mW/m3), and K(z) the thermal conductivity (W/m1 C) expressed as the sum of a deterministic component K and random component K0 (z). K0 (z) has mean zero and a Gaussian colored noise correlation structure represented by EðK 0 ðzÞÞ ¼ 0 EðK 0 ðz1 ÞK 0 ðz2 ÞÞ ¼ s2K
(2) erjz1 z2 j
(3)
where sK2 is the variance in thermal conductivity, r the correlation decay parameter (or 1/r the correlation length scale) and z1 and z2 the depths. The analytical expressions derived by Srivastava et al. (2006) for the temperature depth distribution along with its error bounds for two sets of boundary conditions, i.e. (i) surface temperature and surface heat flow Qs (mW/m2) and (ii) surface temperature and basal heat flow Q* (mW/m2), have been used in this study. The solution to Eq. (1) for the first set of boundary conditions is Mean temperature 2
Q A0 D T ¼ EðTðzÞÞ ¼ T 0 þ s z þ K K
z 1 ez=D D
r
þ
ezðrþ1=DÞ
r2
½rz þ erz 1
½zD þ D2 ez=D D2 g Term4 ¼
1 2
r
þ
2 3 ðrz þ 1Þ rz ðrz þ erz 1Þ 3 r2
erz
r2
½rz þ erz 1g
For the second set of boundary condition i.e. surface temperature and basal heat flow Srivastava (2005) have given the solution to the mean heat flow and its variance as Mean temperature
Q A0 D zþ T ¼ EðTðzÞÞ ¼ T 0 þ K K
2
z 1 eL=D ez=D D
Variance in the temperature
s2T ¼ c1 Term1 þ c2 Term2 þ c3 Term3 þ c4 Term4
(7)
where the constants are different and are given as (4)
(6)
c1 ¼ C 2K A20 ð1 rDÞ2 =K
2
ARTICLE IN PRESS K. Srivastava et al. / Computers & Geosciences 35 (2009) 2095–2099 2 c2 ¼ C 2K A0 rðrD 1ÞðQ A0 DeL=D Þ=K
c3 ¼ c2 2 c4 ¼ C 2K r2 ðQ A0 DeL=D Þ2 =K
The terms Term1, Term2, Term3 and Term4 are same as given above.
3. Numerical examples and discussion The analytical solutions given by Srivastava et al. (2006) for the mean and variance of the temperature depth distribution for two different sets of boundary conditions have been used to quantify the mean temperature along with its error bounds for the seismically active Latur region of India. The Latur region in the Deccan volcanic province in central India experienced a devastating earthquake of magnitude Mw ¼ 6.2 on 30th September, 1993. The Deccan volcanic province is covered by a thin layer of continental flood basalts. The earthquake occurred in a seismically quiet zone, highlighting the unstable nature of the highly fragmented Indian Shield. The focal depth of the earthquake as suggested by USGS was around 9 km, and the focal mechanism was a reverse fault with NW–SE-orientated nodal planes. This earthquake was studied by several researchers from different geophysical angles. Singh et al. (1995) attributed the intraplate seismicity in the Deccan volcanic province of India to the presence of fluids in the upper and middle crust. Gupta et al. (1996) carried out various geophysical investigations of the crustal structure in
2097
the region. From magnetotelluric studies they found a fluid-filled conducting zone at 6–10 km depth that would have enhanced the accumulation of stresses in the crust, causing a mechanical failure. Agrawal and Pandey (1999) provided an alternative explanation in which the seismic activity stems from a hot underlying asthenosphere that causes continuous build up of localized stresses. Patro et al. (2006) suggested that the seismicity occurred in a relatively thick and highly resistive upper crust. From a thermal perspective, Roy and Rao (1999) carried out heat flow measurements in four bore holes in the vicinity of the rupture zone. One of the bore holes was drilled in the surface rupture zone and penetrated 338 m in the Deccan basalts and 271 m further down into the Archaean granite-gneiss basement. They found that the surface heat flow in the region is about 43 mW/m2, which is characteristic of a low-heat-flow regime. Gupta et al. (1990) characterized the Dharwar craton by low-heatflow values, ranging from 25 to 50 mW/m2, similar to values in other Archaean terrains in the world. In this study, we have used an exponential model for the heat generation as it has been used extensively in the literature to quantify the thermal structure. The mean radiogenic heat production has been taken to be about 2.6 mW/m3, bearing in mind the dominance of granite in the crust. The characteristic depth has been taken to be about 12 km. From deep seismic sounding results Kaila et al. (1979) found an average crustal thickness of about 37 km for the Dharwar craton. Thermal conductivity measurements for the bore hole samples have been found to vary from 2.6 to 3.5 W/m 1C (Roy and Rao, 1999), so a mean thermal conductivity value of about 3.0 W/m 1C has been taken in this study. The random component parameters are the coefficient of variability in thermal conductivity which is
Fig. 1. Plot of mean temperature and mean 7 1 S.D. for Set 1 when boundary conditions employed are surface temperature and surface heat flow.
ARTICLE IN PRESS 2098
K. Srivastava et al. / Computers & Geosciences 35 (2009) 2095–2099
taken to be 0.3 and the correlation length scale which is considered to be about 15 km. The basal heat flow estimated for this region is about 12 mW/m2 and the surface temperature for the region is considered to be 30 1C. In the GUI the controlling input parameters can be given directly on the screen and both the plot and the results are displayed instantaneously. For the Latur region, we used the controlling parameters as given above. Set 1 uses the surface temperature and surface heat flow as the boundary conditions and Set 2 uses the surface temperature and basal heat flow as its boundary conditions. The coefficient of variability in thermal conductivity is in the range 0–1. To incorporate a 30% error in the thermal conductivity we need to take the coefficient of variability to be 0.3. The correlation length scale should be less than the model length. In many situations, one third of the model length for the correlation length scale is a good approximation. The input-controlling parameters are given in the boxes and the output of the temperature depth profile along with the mean 7 one standard deviation is plotted; also the results of depth, mean temperature, standard deviation, mean plus one standard deviation and mean minus one standard deviation are displayed on the screen. Numerical values for Set 1 Crustal thickness (L) Radiogenic heat production (A) Characteristic depth (D) Surface temperature (T0) Surface heat flow (Qs)
37 km 2.6 mW/m3 12 km 30 1C 43 mW/m2
Random thermal conductivity Mean thermal conductivity K Coefficient of variability Ck Correlation length scale 1/r
3.0 W/m 1C 0.3 15 km
The results obtained are plotted in Fig. 1. From the figure, we observe that the mean temperature at the base of the crust (37 km) is about 295 1C and the standard deviation is about 29 1C. The crustal temperature here varies between 323 and 266 1C. Numerical values for Set 2 Crustal thickness (L) Radiogenic heat production (A) Characteristic depth (D) Surface temperature (T0) Basal heat flow (Q*)
37 km 2.6 mW/m3 12 km 30 1C 12 mW/m2
Random thermal conductivity Mean thermal conductivity K Coefficient of variability Ck Correlation length scale 1/r
3.0 W/m 1C 0.3 15 km
Results obtained are plotted in Fig. 2. From the figure, we observe that at the base of the crust the mean temperature is about 276 1C and the standard deviation is about 24 1C. We observe that there is a difference in the temperatures between
Fig. 2. Plot of mean temperature and mean 7 1 S.D. for Set 2 when boundary conditions employed are surface temperature and basal heat flow.
ARTICLE IN PRESS K. Srivastava et al. / Computers & Geosciences 35 (2009) 2095–2099
Sets 1 and 2. This is due to the different boundary conditions employed. Also the error bounds are seen to be highly dependent on the coefficient of variability and the correlation length scale, as seen from the closed form expressions. The mean temperatures obtained fall well within the estimated deterministic temperatures obtained by earlier workers. Previous studies have examined the thermal structure from a deterministic point of view. In this study, only the thermal conductivity has been considered to be a random parameter. Other parameters such as radiogenic heat sources and boundary conditions can also be considered to be random. Obtaining an analytic solution to a stochastic differential equation with more than two random parameters would be difficult. However, using the perturbation methods or Monte carlo simulation one can obtain the solution to the problem. Sometimes the solution to a simple stochastic equation with one random parameter is essential as it gives an average picture, and the standard deviation (which is the variability indicator) indicates the errors associated with the system behavior due to errors in the system input. Hence, this study has helped with quantifying the errors in the temperature depth distribution, and the software package developed can be used for any given region. 4. Conclusions A package has been developed to compute and plot the error bounds on subsurface temperatures due to errors in the thermal conductivity for a 1-D steady-state conductive earth model for two different sets of boundary conditions and an exponential heat source. The analytical expressions for error bounds on the temperature depth distribution as given by Srivastava et al. (2006) have been used and a graphical user interface has been developed in MatLaboratory. This is a very user-friendly package that allows the controlling input thermal parameters such as crustal thickness, radiogenic heat production, characteristic depth, surface temperature, surface heat flow, mean thermal conductivity, coefficient of variability in thermal conductivity and correlation length scale to be specified and these can be shown directly on the screen. A plot of mean temperature along with its error bounds is also displayed directly on the screen. The developed package has been applied to quantify the conductive geotherms along with their error bounds for the Latur earthquake region of India. The mean temperatures at the base of the crust in the region are seen to be between 260 and 325 1C indicating that the region is characteristic of a cold type of crust. The package can be applied to quantify the thermal state of the crust, along with its error bounds, for any given region.
Acknowledgements The authors wish to thank the students of the Department of Mathematics, Osmania University for their help. The comments
2099
and suggestions made by two anonymous reviewers have helped us in improving the paper.
Appendix A. Supporting Information Supplementary data associated with this article can be found in the online version at doi:10.1016/j.cageo.2008.12.013.
References Adomian, G., 1994. Solving Frontier Problems in Physics—The Decomposition Method. Kluwer Academic, Boston, 372pp. Agrawal, P.K., Pandey, O.P., 1999. Relevance of hot underlying aesthenosphere to the occurrence of Latur earthquake and Indian peninsular shield seismicity. Geodynamics 28, 303–316. Gallagher, K., Ramsdale, M., Lonergan, L., Marrow, D., 1997. The role thermal conductivities measurements in modeling the thermal histories in sedimentary basins. Marine and Petroleum Geology 14, 201–214. Gupta, M.L., Sunder, A., Sarma, S.R., 1990. Heat flow and heat generation in Dharwar cratons and implications for the southern Indian shield geotherm and lithospheric thickness. Tectonophysics 194, 107–122. Gupta, H.K., Sarma, S.V.S., Harinarayana, T., Virupakshi, G., 1996. Fluids below the hypocentral region of Latur earthquake India, geophysical indicators. Geophysical. Research Letters 23, 1569–1572. Jokinen, J., Kukkonen, I.T., 1999a. Random modeling of lithospheric thermal regime: Forward simulation applied in uncertainty analysis. Tectonophysics 306, 277–292. Jokinen, J., Kukkonen, I.T., 1999b. Inverse simulation of lithospheric thermal regime using the Monte Carlo method. Tectonophysics 306, 293–310. Kaila, K.L., Roy Choudhury, K., Reddy, P.R., Krishna, V.G., Narain, H., Subbotin, S.I., Sollogub, V.B., Chekunov, A.V., Kharetchko, G.E., Lazarenko, M.A., Ilchenko, T.V., 1979. Crustal structure along Kavali–Udipi profile in the Indian peninsular shield from deep seismic sounding. Journal Geological Society of India 20, 307–333. Nielson, S.B., 1987. Steady state heat flow in a random medium and linear heat flow heat production relationship. Geophysical Research Letters 14, 318–321. Patro, B.P.K., Nagarajan, N., Sarma, S.V.S., 2006. Crustal geoelectric structure and the focal depths of major stable continental region earthquakes in India. Current Science 90 (1), 107–113. Roy, S., Rao, R.U.M., 1999. Geothermal investigations in the 19993 Latur earthquake area, Deccan Volcanic Province, India. Tectonophysics 306, 237–252. Royer, J.J., Danis, M., 1988. Steady state geothermal model of the crust and problems of boundary conditions: Application to a rift system, the southern Rhinegraben. Tectonophysics 156, 239–255. Serrano, S.E., 1995. Forecasting scale dependent dispersion from spills in heterogeneous aquifers. Journal of Hydrology 169, 151–169. Singh, R.P., Sato, T., Nyland, E., 1995. The geodynamic context of the Latur (India) earthquake, 30, September 1993. Physics of the Earth and Planetary Interiors 91, 245–251. Srivastava, K., 2005. Modeling the variability of heat flow due to random thermal conductivity of the crust. Geophysical Journal International 160 (2), 776–782. Srivastava, K., Singh, R.N., 1998. A model for temperature variation in sedimentary basins due to random radiogenic heat sources. Geophysical Journal International 135, 727–730. Srivastava, K., Singh, R.N., 1999. A stochastic model to quantify the steady state crustal geotherms subject to uncertainty in thermal conductivity. Geophysical Journal International 138, 895–899. Vasseur, G., Lucazeau, F., Bayer, R., 1985. The problem of heat flow density determination from inaccurate data. Tectonophysics 121, 23–34. Vasseur, G., Singh, R.N., 1986. Effect of random horizontal variation in radiogenic heat source distribution on its relationship with heat flow. Journal Geophysical Research 91, 10397–10404. Witten, A., 2004. A Matlab based three dimensional viewer. Computers & Geosciences 30, 693–703.
Contents lists available at ScienceDirect
Computers & Geosciences journal homepage: www.elsevier.com/locate/cageo
Quantification of crustal geotherms along with their error bounds for seismically active regions of India: A Matlab toolbox$ Kirti Srivastava , Likhita Narain, Swaroopa Rani, V.P. Dimri National Geophysical Research Institute, Council of Scientific and Industrial Research, Hyderabad 500007, Andhra Pradesh, India
a r t i c l e in fo
abstract
Article history: Received 11 September 2008 Received in revised form 24 December 2008 Accepted 29 December 2008
Determination of the thermo-mechanical structure of the crust for seismically active regions using available geophysical and geological data is of great importance. The most important feature of the intraplate earthquakes in the Indian region are that the seismicity occurs within the entire crust. In Latur region of India an earthquake occurred in the upper crust. In such situations, quantifying the uncertainties in seismogenic depths becomes very important. The stochastic heat conduction equation has been solved for different sets of boundary conditions, an exponentially decreasing radiogenic heat generation and randomness in thermal conductivity. Closed form analytical expressions for mean and variance in the temperature depth distribution have been used and an automatic formulation has been developed in Matlab for computing and plotting the thermal structure. The Matlab toolbox presented allows us to display the controlling thermal parameters on the screen directly, and plot the subsurface thermal structure along with its error bounds. The software can be used to quantify the thermal structure for any given region and is applied here to the Latur earthquake region of India. & 2009 Elsevier Ltd. All rights reserved.
Keywords: Stochastic Gaussian Random thermal conductivity Heat conduction Matlab
1. Introduction The thermal structure of the Earth’s crust is influenced by its geothermal parameters such as thermal conductivity, radiogenic heat sources and initial and boundary conditions. Two approaches to modeling are commonly used to estimate the subsurface temperature field: (1) the deterministic approach and (2) the stochastic approach. In the deterministic approach the subsurface temperature field is obtained by assuming that the controlling thermal parameters are known with certainty. However, due to the inhomogeneous nature of the Earth’s interior there are some uncertainties in estimates of the geothermal parameters. Uncertainties in these parameters may arise from the measurement inaccuracies or a lack of information about the parameters themselves. In the stochastic approach the uncertainties in parameters are incorporated and an average picture of the thermal field along with its associated error bounds is determined. To assess the properties of the system at a glance we need to obtain the mean value that gives the average picture, along with the variance or standard deviation that indicate variability or errors associated with the system behavior due to errors in the system input.
$
Code available from server at http://www.iamg.org/CGEditor/index.htm
Corresponding author.
E-mail address: [email protected] (K. Srivastava). 0098-3004/$ - see front matter & 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.cageo.2008.12.013
Subsurface temperatures are seen to be very sensitive to perturbations in the input thermal parameters and hence several studies have been carried out to quantify perturbations in the temperature and heat flow using stochastic analytical and random simulation techniques. Uncertainty in the heat flow was quantified by Vasseur et al. (1985) using a least squares inversion technique incorporating uncertainties in the temperature and thermal conductivities. The effect of variations in heat source on the surface heat flow has also been studied (Vasseur and Singh, 1986; Nielson, 1987). In most studies the stochastic heat equation was solved using the small perturbation method. Srivastava and Singh (1998) used the small perturbation method to analytically solve the steady-state heat conduction equation, incorporating Gaussian uncertainties in the heat sources, and obtained the closed form analytical expressions for the mean temperature field and its variance. The thermal structure is obtained numerically using the random simulation method or a Monte Carlo simulation. Using a Monte Carlo simulation, several authors have studied the subsurface thermal field along with its error structure by incorporating randomness in the controlling thermal parameters such as thermal conductivity, radiogenic heat production and boundary conditions (Royer and Danis, 1988; Gallagher et al., 1997; Jokinen and Kukkonen, 1999a, b). This numerical modeling is very useful when studying nonlinear problems, but sometimes a simple 1-D analytical solution to the mean behavior and its associated error bounds is more useful for quantifying the uncertainty.
ARTICLE IN PRESS 2096
K. Srivastava et al. / Computers & Geosciences 35 (2009) 2095–2099
Stochastic differential equations in other fields are now being solved by yet another approach called the decomposition method (Adomian, 1994; Serrano, 1995). Adomian’s method of decomposition has now been generalized as a general analytic procedure to solve deterministic or stochastic and linear or nonlinear equations. It has been shown to be systematic, robust, and sometimes capable of handling large variances in the controlling parameters. In a recent study, using this new approach the stochastic heat conduction equation was solved by Srivastava and Singh (1999), incorporating uncertainties in the thermal conductivity, where the solution for the temperature field was obtained using a series expansion method. The thermal conductivity was considered to be a random parameter with a known Gaussian colored noise correlation structure. Later, Srivastava (2005) extended the study to obtain the analytical expressions for mean heat flow and its variance. Graphical user interface (GUI) packages in Matlab are becoming very popular with geoscientific researchers. Witten (2004) developed a sequence of Matlab m-files and two graphical user interfaces to display raw or processed geophysical data to produce the final graphics. In this paper, we have developed a simple graphical user interface viewer which is a MATLAB-based software consisting of m-files. This software is the first package that computes the thermal structure along with its error structure, with the important advantage that it is not a compiled code. The m-file in the package is integrated through a GUI and the controlling thermal parameters such as crustal thickness, radiogenic heat production, characteristic depth, surface temperature, surface heat flow, mean thermal conductivity, coefficient of variability in thermal conductivity and correlation length scale are all given on the screen. The program model is intended as a user friendly tool to deal with various aspects of analysis and interpretation of subsurface geothermal structure for any given region.
Variance in the temperature 2 T
s ¼ c1 Term1 þ c2 Term2 þ c3 Term3 þ c4 Term4
(5)
where CK ¼ sK/K is the coefficient of variability in thermal conductivity and c1 ¼ C 2K A20 ð1 rDÞ2 =K
2
2 c2 ¼ C 2K A0 rðrD 1ÞðQ s A0 DÞ=K
c3 ¼ c2 2 c4 ¼ C 2K r2 ðQ s A0 DÞ2 =K
Term1 ¼
1 fðrD 1Þð2z2 2zD D2 e2z=D þ D2 Þ 4ðr 1=DÞ2 4½zðr 1=DÞ þ 1 þ ½zðr þ 1=DÞ ezðrþ1=DÞ þ 1 ðr þ 1=DÞ2 þ ½2zD þ D2 e2z=D D2 g 1 þ fðrD þ 1Þð2z2 2zD D2 e2z=D þ D2 Þ 4ðr þ 1=DÞ2 4 þ ½zðr 1=DÞezðrþ1=DÞ þ e2z=D ezðrþ1=DÞ ðr 1=DÞ2 ½2zD þ D2 e2z=D D2 g
Term2 ¼
1
r2 þ
f2rðz2 D 2zD2 2D3 ez=D þ 2D3 Þ ð1 þ rzÞ ðr þ 1=DÞ2 erz ðr 1=DÞ2
½zðr þ 1=DÞ þ ezðrþ1=DÞ 1 ½zðr 1=DÞ þ ezðr1=DÞ 1g
2. Mathematical formulation
1 fðr 1=DÞðz2 D 2zD2 2D3 ez=D þ 2D3 Þ ðr 1=DÞ2 zðr 1=DÞ þ 1 ½rz þ erz 1 2
The heat conduction equation with random thermal conductivity can be expressed as d dT (1) ðK þ K 0 ðzÞÞ ¼ A0 ez=D dz dz
þ ½zD þ D2 ez=D D2 g 1 þ fðr þ 1=DÞðz2 D 2zD2 2D3 ez=D þ 2D3 Þ ðr þ 1=DÞ2
Term3 ¼
where T is the temperature (1C), A(z) the radiogenic heat source (mW/m3), and K(z) the thermal conductivity (W/m1 C) expressed as the sum of a deterministic component K and random component K0 (z). K0 (z) has mean zero and a Gaussian colored noise correlation structure represented by EðK 0 ðzÞÞ ¼ 0 EðK 0 ðz1 ÞK 0 ðz2 ÞÞ ¼ s2K
(2) erjz1 z2 j
(3)
where sK2 is the variance in thermal conductivity, r the correlation decay parameter (or 1/r the correlation length scale) and z1 and z2 the depths. The analytical expressions derived by Srivastava et al. (2006) for the temperature depth distribution along with its error bounds for two sets of boundary conditions, i.e. (i) surface temperature and surface heat flow Qs (mW/m2) and (ii) surface temperature and basal heat flow Q* (mW/m2), have been used in this study. The solution to Eq. (1) for the first set of boundary conditions is Mean temperature 2
Q A0 D T ¼ EðTðzÞÞ ¼ T 0 þ s z þ K K
z 1 ez=D D
r
þ
ezðrþ1=DÞ
r2
½rz þ erz 1
½zD þ D2 ez=D D2 g Term4 ¼
1 2
r
þ
2 3 ðrz þ 1Þ rz ðrz þ erz 1Þ 3 r2
erz
r2
½rz þ erz 1g
For the second set of boundary condition i.e. surface temperature and basal heat flow Srivastava (2005) have given the solution to the mean heat flow and its variance as Mean temperature
Q A0 D zþ T ¼ EðTðzÞÞ ¼ T 0 þ K K
2
z 1 eL=D ez=D D
Variance in the temperature
s2T ¼ c1 Term1 þ c2 Term2 þ c3 Term3 þ c4 Term4
(7)
where the constants are different and are given as (4)
(6)
c1 ¼ C 2K A20 ð1 rDÞ2 =K
2
ARTICLE IN PRESS K. Srivastava et al. / Computers & Geosciences 35 (2009) 2095–2099 2 c2 ¼ C 2K A0 rðrD 1ÞðQ A0 DeL=D Þ=K
c3 ¼ c2 2 c4 ¼ C 2K r2 ðQ A0 DeL=D Þ2 =K
The terms Term1, Term2, Term3 and Term4 are same as given above.
3. Numerical examples and discussion The analytical solutions given by Srivastava et al. (2006) for the mean and variance of the temperature depth distribution for two different sets of boundary conditions have been used to quantify the mean temperature along with its error bounds for the seismically active Latur region of India. The Latur region in the Deccan volcanic province in central India experienced a devastating earthquake of magnitude Mw ¼ 6.2 on 30th September, 1993. The Deccan volcanic province is covered by a thin layer of continental flood basalts. The earthquake occurred in a seismically quiet zone, highlighting the unstable nature of the highly fragmented Indian Shield. The focal depth of the earthquake as suggested by USGS was around 9 km, and the focal mechanism was a reverse fault with NW–SE-orientated nodal planes. This earthquake was studied by several researchers from different geophysical angles. Singh et al. (1995) attributed the intraplate seismicity in the Deccan volcanic province of India to the presence of fluids in the upper and middle crust. Gupta et al. (1996) carried out various geophysical investigations of the crustal structure in
2097
the region. From magnetotelluric studies they found a fluid-filled conducting zone at 6–10 km depth that would have enhanced the accumulation of stresses in the crust, causing a mechanical failure. Agrawal and Pandey (1999) provided an alternative explanation in which the seismic activity stems from a hot underlying asthenosphere that causes continuous build up of localized stresses. Patro et al. (2006) suggested that the seismicity occurred in a relatively thick and highly resistive upper crust. From a thermal perspective, Roy and Rao (1999) carried out heat flow measurements in four bore holes in the vicinity of the rupture zone. One of the bore holes was drilled in the surface rupture zone and penetrated 338 m in the Deccan basalts and 271 m further down into the Archaean granite-gneiss basement. They found that the surface heat flow in the region is about 43 mW/m2, which is characteristic of a low-heat-flow regime. Gupta et al. (1990) characterized the Dharwar craton by low-heatflow values, ranging from 25 to 50 mW/m2, similar to values in other Archaean terrains in the world. In this study, we have used an exponential model for the heat generation as it has been used extensively in the literature to quantify the thermal structure. The mean radiogenic heat production has been taken to be about 2.6 mW/m3, bearing in mind the dominance of granite in the crust. The characteristic depth has been taken to be about 12 km. From deep seismic sounding results Kaila et al. (1979) found an average crustal thickness of about 37 km for the Dharwar craton. Thermal conductivity measurements for the bore hole samples have been found to vary from 2.6 to 3.5 W/m 1C (Roy and Rao, 1999), so a mean thermal conductivity value of about 3.0 W/m 1C has been taken in this study. The random component parameters are the coefficient of variability in thermal conductivity which is
Fig. 1. Plot of mean temperature and mean 7 1 S.D. for Set 1 when boundary conditions employed are surface temperature and surface heat flow.
ARTICLE IN PRESS 2098
K. Srivastava et al. / Computers & Geosciences 35 (2009) 2095–2099
taken to be 0.3 and the correlation length scale which is considered to be about 15 km. The basal heat flow estimated for this region is about 12 mW/m2 and the surface temperature for the region is considered to be 30 1C. In the GUI the controlling input parameters can be given directly on the screen and both the plot and the results are displayed instantaneously. For the Latur region, we used the controlling parameters as given above. Set 1 uses the surface temperature and surface heat flow as the boundary conditions and Set 2 uses the surface temperature and basal heat flow as its boundary conditions. The coefficient of variability in thermal conductivity is in the range 0–1. To incorporate a 30% error in the thermal conductivity we need to take the coefficient of variability to be 0.3. The correlation length scale should be less than the model length. In many situations, one third of the model length for the correlation length scale is a good approximation. The input-controlling parameters are given in the boxes and the output of the temperature depth profile along with the mean 7 one standard deviation is plotted; also the results of depth, mean temperature, standard deviation, mean plus one standard deviation and mean minus one standard deviation are displayed on the screen. Numerical values for Set 1 Crustal thickness (L) Radiogenic heat production (A) Characteristic depth (D) Surface temperature (T0) Surface heat flow (Qs)
37 km 2.6 mW/m3 12 km 30 1C 43 mW/m2
Random thermal conductivity Mean thermal conductivity K Coefficient of variability Ck Correlation length scale 1/r
3.0 W/m 1C 0.3 15 km
The results obtained are plotted in Fig. 1. From the figure, we observe that the mean temperature at the base of the crust (37 km) is about 295 1C and the standard deviation is about 29 1C. The crustal temperature here varies between 323 and 266 1C. Numerical values for Set 2 Crustal thickness (L) Radiogenic heat production (A) Characteristic depth (D) Surface temperature (T0) Basal heat flow (Q*)
37 km 2.6 mW/m3 12 km 30 1C 12 mW/m2
Random thermal conductivity Mean thermal conductivity K Coefficient of variability Ck Correlation length scale 1/r
3.0 W/m 1C 0.3 15 km
Results obtained are plotted in Fig. 2. From the figure, we observe that at the base of the crust the mean temperature is about 276 1C and the standard deviation is about 24 1C. We observe that there is a difference in the temperatures between
Fig. 2. Plot of mean temperature and mean 7 1 S.D. for Set 2 when boundary conditions employed are surface temperature and basal heat flow.
ARTICLE IN PRESS K. Srivastava et al. / Computers & Geosciences 35 (2009) 2095–2099
Sets 1 and 2. This is due to the different boundary conditions employed. Also the error bounds are seen to be highly dependent on the coefficient of variability and the correlation length scale, as seen from the closed form expressions. The mean temperatures obtained fall well within the estimated deterministic temperatures obtained by earlier workers. Previous studies have examined the thermal structure from a deterministic point of view. In this study, only the thermal conductivity has been considered to be a random parameter. Other parameters such as radiogenic heat sources and boundary conditions can also be considered to be random. Obtaining an analytic solution to a stochastic differential equation with more than two random parameters would be difficult. However, using the perturbation methods or Monte carlo simulation one can obtain the solution to the problem. Sometimes the solution to a simple stochastic equation with one random parameter is essential as it gives an average picture, and the standard deviation (which is the variability indicator) indicates the errors associated with the system behavior due to errors in the system input. Hence, this study has helped with quantifying the errors in the temperature depth distribution, and the software package developed can be used for any given region. 4. Conclusions A package has been developed to compute and plot the error bounds on subsurface temperatures due to errors in the thermal conductivity for a 1-D steady-state conductive earth model for two different sets of boundary conditions and an exponential heat source. The analytical expressions for error bounds on the temperature depth distribution as given by Srivastava et al. (2006) have been used and a graphical user interface has been developed in MatLaboratory. This is a very user-friendly package that allows the controlling input thermal parameters such as crustal thickness, radiogenic heat production, characteristic depth, surface temperature, surface heat flow, mean thermal conductivity, coefficient of variability in thermal conductivity and correlation length scale to be specified and these can be shown directly on the screen. A plot of mean temperature along with its error bounds is also displayed directly on the screen. The developed package has been applied to quantify the conductive geotherms along with their error bounds for the Latur earthquake region of India. The mean temperatures at the base of the crust in the region are seen to be between 260 and 325 1C indicating that the region is characteristic of a cold type of crust. The package can be applied to quantify the thermal state of the crust, along with its error bounds, for any given region.
Acknowledgements The authors wish to thank the students of the Department of Mathematics, Osmania University for their help. The comments
2099
and suggestions made by two anonymous reviewers have helped us in improving the paper.
Appendix A. Supporting Information Supplementary data associated with this article can be found in the online version at doi:10.1016/j.cageo.2008.12.013.
References Adomian, G., 1994. Solving Frontier Problems in Physics—The Decomposition Method. Kluwer Academic, Boston, 372pp. Agrawal, P.K., Pandey, O.P., 1999. Relevance of hot underlying aesthenosphere to the occurrence of Latur earthquake and Indian peninsular shield seismicity. Geodynamics 28, 303–316. Gallagher, K., Ramsdale, M., Lonergan, L., Marrow, D., 1997. The role thermal conductivities measurements in modeling the thermal histories in sedimentary basins. Marine and Petroleum Geology 14, 201–214. Gupta, M.L., Sunder, A., Sarma, S.R., 1990. Heat flow and heat generation in Dharwar cratons and implications for the southern Indian shield geotherm and lithospheric thickness. Tectonophysics 194, 107–122. Gupta, H.K., Sarma, S.V.S., Harinarayana, T., Virupakshi, G., 1996. Fluids below the hypocentral region of Latur earthquake India, geophysical indicators. Geophysical. Research Letters 23, 1569–1572. Jokinen, J., Kukkonen, I.T., 1999a. Random modeling of lithospheric thermal regime: Forward simulation applied in uncertainty analysis. Tectonophysics 306, 277–292. Jokinen, J., Kukkonen, I.T., 1999b. Inverse simulation of lithospheric thermal regime using the Monte Carlo method. Tectonophysics 306, 293–310. Kaila, K.L., Roy Choudhury, K., Reddy, P.R., Krishna, V.G., Narain, H., Subbotin, S.I., Sollogub, V.B., Chekunov, A.V., Kharetchko, G.E., Lazarenko, M.A., Ilchenko, T.V., 1979. Crustal structure along Kavali–Udipi profile in the Indian peninsular shield from deep seismic sounding. Journal Geological Society of India 20, 307–333. Nielson, S.B., 1987. Steady state heat flow in a random medium and linear heat flow heat production relationship. Geophysical Research Letters 14, 318–321. Patro, B.P.K., Nagarajan, N., Sarma, S.V.S., 2006. Crustal geoelectric structure and the focal depths of major stable continental region earthquakes in India. Current Science 90 (1), 107–113. Roy, S., Rao, R.U.M., 1999. Geothermal investigations in the 19993 Latur earthquake area, Deccan Volcanic Province, India. Tectonophysics 306, 237–252. Royer, J.J., Danis, M., 1988. Steady state geothermal model of the crust and problems of boundary conditions: Application to a rift system, the southern Rhinegraben. Tectonophysics 156, 239–255. Serrano, S.E., 1995. Forecasting scale dependent dispersion from spills in heterogeneous aquifers. Journal of Hydrology 169, 151–169. Singh, R.P., Sato, T., Nyland, E., 1995. The geodynamic context of the Latur (India) earthquake, 30, September 1993. Physics of the Earth and Planetary Interiors 91, 245–251. Srivastava, K., 2005. Modeling the variability of heat flow due to random thermal conductivity of the crust. Geophysical Journal International 160 (2), 776–782. Srivastava, K., Singh, R.N., 1998. A model for temperature variation in sedimentary basins due to random radiogenic heat sources. Geophysical Journal International 135, 727–730. Srivastava, K., Singh, R.N., 1999. A stochastic model to quantify the steady state crustal geotherms subject to uncertainty in thermal conductivity. Geophysical Journal International 138, 895–899. Vasseur, G., Lucazeau, F., Bayer, R., 1985. The problem of heat flow density determination from inaccurate data. Tectonophysics 121, 23–34. Vasseur, G., Singh, R.N., 1986. Effect of random horizontal variation in radiogenic heat source distribution on its relationship with heat flow. Journal Geophysical Research 91, 10397–10404. Witten, A., 2004. A Matlab based three dimensional viewer. Computers & Geosciences 30, 693–703.