- Email: [email protected]
Crack density, saturation rate and porosity at the 2001 Bhuj, India, earthquake hypocenter: a fluid-driven earthquake?
Available online at www.sciencedirect.com R
Earth and Planetary Science Letters 212 (2003) 393^405 www.elsevier.com/locate/epsl
Crack density, saturation rate and porosity at the 2001 Bhuj, India, earthquake hypocenter: a £uid-driven earthquake? O.P. Mishra , Dapeng Zhao Geodynamics Research Center, Ehime University, Matsuyama 790-8577, Japan Received 10 December 2002; received in revised form 12 May 2003; accepted 14 May 2003
Abstract Quantitative estimates of the crack density (O), saturation rate (h) and porosity (i) parameters from seismic velocities (Vp , Vs ) and Poisson’s ratio (c) in the 2001 Bhuj earthquake area in western India indicate that the 2001 Bhuj earthquake hypocenter is associated with high-O, high-h and high-i in the depth range of about 23^28 km, extending 15^30 km laterally. These anomalies may be due to a fluid-filled, fractured rock matrix, which might have contributed to trigger the 2001 Bhuj earthquake in the intraplate, stable continental region of the Indian Peninsula. This feature is similar to that of the 1995 Kobe earthquake [D. Zhao et al., Science 274 (1996) 1891^1894]. High-i areas are generally consistent with high-O areas, but high-h areas have a wider distribution, indicating that microcracks exist in more localized areas within the crust, and that permeation of fluids in the hypocenter zone might have occurred extensively through the intergranular and fractured pores due to hidden intersecting fault geometry. Here we suggest the possibility that earthquake occurrence is closely related to in situ material heterogeneities, rather than stress conditions alone. @ 2003 Elsevier Science B.V. All rights reserved. Keywords: crack density; saturation rate; porosity; £uids; crustal earthquakes
1. Introduction The January 26, 2001 Bhuj earthquake (Mw 7.6) in western India has been investigated by several researchers using multidisciplinary approaches, making it one of the best studied Indian earthquakes in the annals of seismology. It was one of the most catastrophic earthquakes in In-
* Corresponding author. Tel.: +81-89-927-8257; Fax: +81-89-927-8167. E-mail addresses: [email protected] (O.P. Mishra), [email protected] (D. Zhao).
dia, which rocked Gujarat province at Bhuj (23.4‡N, 70.3‡E) in western India, with an estimated death toll of over 20 000 [2]. The 2001 Bhuj mainshock was followed by a series of aftershock events (over 4000), magnitude v 1.0, which were recorded during the period from January 30 to April 15, 2001 [3]. Arrival times generated from 368 selected aftershocks have been used (Fig. 1) for tomographic estimations of three-dimensional (3-D) P-wave (Vp ) and S-wave (Vs ) velocity and Poisson’s ratio (c) structures in the hypocentral area of the Bhuj earthquake to understand what may have triggered the Bhuj mainshock [4,5]. To better understand the relationship between cracks/
0012-821X / 03 / $ ^ see front matter @ 2003 Elsevier Science B.V. All rights reserved. doi:10.1016/S0012-821X(03)00285-1
EPSL 6696 2-7-03 Cyaan Magenta Geel Zwart
394
O.P. Mishra, D. Zhao / Earth and Planetary Science Letters 212 (2003) 393^405
£uids and the earthquake nucleation processes, we have attempted to estimate the 3-D distribution of crack density (O), saturation rate (h) and porosity (i) parameters in the source volume of the 2001 Bhuj earthquake, by applying the crack theory [6] to the values of Vp , Vs and c derived from a high resolution, 3-D tomographic inversion [7]. The spatial and temporal distribution of the events in an aftershock sequence contains information about the aftershock-generating process in particular, and by extrapolation, the earthquake-generating process in general, because the distribution of aftershocks is closely linked to the geometry of the fault zone that ruptured to produce the mainshock, which in turn depends on the physical properties of the fault zone as well as ambient conditions, especially the distribution of strength, stress, and temperature [8]. In this study, aftershock events are used as a source for imaging the subsurface structure from a particular seismic network design (Fig. 1). Aftershock events have previously been used by other researchers to image the 3-D structure in the source zone of di¡erent earthquakes elsewhere in the world to understand what may have triggered these earthquakes, and how the rupture developed after initiation [9^ 13]. Complex crustal structure would have certainly in£uenced the initiation, propagation and termination of the earthquake rupture. If variations in material properties in£uence the seismic velocity, then we should be able to image them in a velocity model [13,14]. Fault zone £uids are widely believed as a key to understanding earthquake generation. Speci¢cally, it has been hypothesized that high £uid pressure in the fault zone acts to reduce the frictional strength of the fault zone, and that time variations in £uid pressure control the timing of earthquakes [15,16]. Our study shows that £uid involvement in the fault zone may be a causative agent in triggering large earthquakes as a consequence of unique £uid pressure, rather than lithologic conditions [10,12^15]. Evidence to date for £uid-driven seismicity for the 2001 Bhuj earthquake includes appreciable seismic velocity (Vp , Vs ) and Poisson’s ratio (c) contrasts in the source zone [4,5], which are well corroborated with the presence of conductive crustal structure as inferred from magne-
totelluric measurements [17]. Evidence for the presence of £uids in other seismogenic fault zones elsewhere in the world is also reported [9^12,18]. 3-D estimates of crack density, saturation rate and porosity parameters in this study may provide additional information that will add to our understanding of the rupture initiation process for the 2001 Bhuj earthquake besides a supportive evidence for £uid-driven earthquake.
2. Tectonic setting and geology The Bhuj mainshock occurred in the Kutch Rift Basin. The Kutch region is underlain by a Mesozoic rift system associated with three major basins (Kutch, Combay, Narmada) in the western margin of the peninsular India craton (Fig. 2a). Mesozoic sediments were uplifted, folded, intruded and covered by Deccan Trap basaltic £ows in Late Cretaceous and early Paleocene time [19]. Faults within such rift systems are known to have the potential to generate large earthquakes [20], and the area has witnessed several large historical earthquakes in the past [21]. The source region of the 2001 Bhuj earthquake is ¢lled with sediments ranging in age from Middle Jurassic to present. A recent study [22] suggests a marine origin for the Jurassic^Cretaceous unit of the Kutch Rift Basin, in which Eocene and older sediments are distributed in highlands and lowlands as shown in Fig. 2b.
3. Data and method Fig. 1 shows the epicentral distribution of 368 Bhuj aftershocks, which were selected from a total of over 4000 aftershocks (M v 1.0). These 368 events generated about 4000 phase data that consist of 1960 P- and 1880 S-wave arrivals, which were used in tomographic imaging. The reliability of P- and S-wave arrival time data were carefully checked using digital and analog recorders [4]. Our reading accuracy varies from 0.01 to 0.1 s for P and 0.05 to 0.3 s for S arrivals. The majority of the events are located by more than four seismic stations from a seismic network consisting of
EPSL 6696 2-7-03 Cyaan Magenta Geel Zwart
O.P. Mishra, D. Zhao / Earth and Planetary Science Letters 212 (2003) 393^405
395
Fig. 1. Epicentral distribution of the 368 Bhuj aftershocks used in this study. The blue circles (aftershocks); the red star (Bhuj mainshock). Lines AB and CD show the locations of the cross-sections in Figs. 3a^c and 7a^c and Fig. 3d^f and 7d^f, respectively. The red triangles denote portable seismic stations. The variation of size of blue circles indicates the variation of magnitude.
12 stations (Fig. 1). The average accuracy of hypocentral locations is better than 4 km. Depth sections of the events show that most of the aftershocks occurred at a depth range of 20^40 km (Fig. 1). We have considered 25 km as the accurate source depth of the 2001 Bhuj mainshock as determined by using converted depth phase data [3^5]. We conducted detailed checkerboard resolution tests [7] for both Vp and Vs at di¡erent grid spac-
ing to assess the resolution of our data set. We found that the present data set is able to image O, h and i (as derived from Vp , Vs and c) structures with sizes of 30 km horizontal and 5 km vertical near the Bhuj mainshock hypocenter [5]. The crack theory of O’Connell and Budiansky [6] (hereafter we call it OB74) is considered as a concise, self-consistent, isotropic theory, capable of estimating crack density and saturation rate in stressed rock from seismic velocity data [12].
EPSL 6696 2-7-03 Cyaan Magenta Geel Zwart
396
O.P. Mishra, D. Zhao / Earth and Planetary Science Letters 212 (2003) 393^405
EPSL 6696 2-7-03 Cyaan Magenta Geel Zwart
O.P. Mishra, D. Zhao / Earth and Planetary Science Letters 212 (2003) 393^405
A set of equations relating changes in Vp and Vs to c [23], and also changes in Vp , Vs and c to O and h are used to estimate O and h. The rock parameters are derived from OB74 on the assumption that £at-circular cracks are aligned randomly in the crust and are ¢lled with £uids or air [6,12]. Changes in Vp and Vs to c can be expressed as: ðV p =V s Þ2 ¼ 2ð13c Þ=ð132c Þ:
O ¼ NGa f;
study region necessarily behave like sedimentary rocks, but rather that the bulk properties of the highly fractured, £uid-saturated regions probably behave like a porous medium. The porosity theory of P63 states that i is directly proportional to the product of the seismic velocities (Vp and Vs ), which can be expressed as: i ¼ V p UV s ;
ð6Þ
ln i ¼ ln V p þ ln V s ;
ð7Þ
d i = i ¼ dV p =V p þ dV s =V s ;
ð8Þ
ð1Þ
According to OB74, the crack density parameter O can be de¢ned as: 3
397
ð2Þ
where N is the number of micro-cracks per unit volume and Gaf is the average radius of circular cracks. If it is assumed that N1 cracks per unit volume are dry and that N2 = N3N1 cracks are saturated, then the saturation rate h is de¢ned as N2 /N. OB74 predicts the following relations: V p2 13c 1 þ c 16 13c 2 ¼ W 13 ð13 h ÞO ð3Þ V 2p 1 þ c 13c 9 132c V s2 32 3 ð13 ¼ 13 c Þ 13 h þ O ð4Þ V 2s 45 23c 45 c 3c 23c ; ð5Þ O¼ 16 13c 2 ð13 h Þð1 þ 3c Þð23c Þ32ð132c Þ
where V p , V s , c are P-wave and S-wave velocities and Poisson’s ratio for cracked volume, respectively, while Vp , Vs , c are the corresponding parameters for uncracked rocks. Thus O and h can be estimated from observed Vp , Vs and c using Eqs. 3^5. The method of [24] based on porosity theory of Pickett [25] (hereafter we call it P63) is used to estimate the porosity parameter (i). The P63 model does not assume that the rocks in the
where di/i is the change of porosity per unit porosity, while dVp /Vp and dVs /Vs are perturbed values of Vp and Vs , respectively. Eq. 8 is used to estimate the distribution of porosity in the 2001 Bhuj earthquake region. The perturbed values of Vp and Vs with respect to an average initial velocity model at di¡erent depth layers (0^30 km) were obtained after setting up a 3-D grid net in the study area [4,5]. We have estimated O, h and i at the same depth layers using the corresponding perturbed values of Vp , Vs and c (Fig. 3). It is imperative to note that high-i denotes low porosity, and the reverse is also true, because i bears a direct relationship to porosity as the product of seismic velocities [25]. Some new and improved models have subsequently been proposed to estimate O and h, due to inapplicability of OB74 for cases involving large crack density and di¡erent pore geometry [26,27]. However, our tomographic results show that the change of seismic velocities is less than 5^6% [4,5], which leads to O 6 0.25, and a corresponding average i perturbation in the range of U 10%. Hence the OB74 and P63 models are viable for this study.
6 Fig. 2. (a) Map showing tectonic features of the study region (compiled from Kayal et al. [3] and Talwani and Gangopadhyay [21]). KRB: Kutch Rift Basin; CRB: Cambay Rift Basin; NRB: Narmada Rift Basin; NPF: Nagar Parkar Fault; IBF: Island Belt Fault; KMF: Kutch Mainland Fault; NKF: North Kathiwar Fault; WCF: West Coast Fault; AB: ‘Allah Bund’ (see text). Epicenters (solid circles) and fault plane solutions of the past signi¢cant earthquakes (M s 5.0) are shown. For fault plane solutions, the solid circle indicates P-axis in the dilatational zone, and the open circle T-axis in the compressional zone. The star indicates the mainshock epicenter of the January 26, 2001 Bhuj earthquake. The insert map shows the major plate boundaries around India and the present study area (adapted from Kayal et al. [4]). (b) Map showing geology of the area (adapted from Talwani and Gangopadhyay [21] after Krishna et al. [22]).
EPSL 6696 2-7-03 Cyaan Magenta Geel Zwart
398
O.P. Mishra, D. Zhao / Earth and Planetary Science Letters 212 (2003) 393^405
Fig. 3. The upper panel (a^c) represents vertical cross-sections along the line AB for Vp , Vs and c, respectively. The lower panel (d^f) same as that of (a^c) but along the line CD. The locations of lines AB and CD are shown in Fig. 1. Small white circles denote the Bhuj aftershocks within 9 km width along the lines AB and CD. Star denotes the hypocenter of Bhuj mainshock; its focal depth is 25.0 km. The perturbation (in %) scales for velocity and Poisson’s ratio are shown at the bottom.
4. Analysis and results We assumed in the theoretical formulations of OB74 that uncracked rocks have the maximum seismic velocity at each depth in the study area. The e¡ects of porosity have been analyzed separately in order to understand the in situ status of pores around the 2001 Bhuj mainshock hypocenter. Average crack porosity for stressed rock is taken as a function of crack orientation and aspect ratio [28]. The distribution of Vp , Vs and c in depths (Fig. 3) along lines AB and CD (Fig. 1) are used to estimate O, h and i. The intersecting fault geometry (Fig. 4) shows bilateral and simultaneous ruptures that propagated along two fault
trends (NE and NW) where the largest (ML 5.7) and second largest (ML 5.5) aftershocks occurred along the NE and NW rupture nucleated trends, respectively [3]. The temporal variation of all aftershocks, recorded by a portable seismic network (Fig. 1) is shown in Fig. 5, upper part, while Fig. 5, three lower parts, show the temporal variations of 368 relocated aftershocks along various sections, which are used in our tomographic study [5]. This ¢gure shows a well-oriented clustering of aftershocks around the mainshock hypocenter. Detailed analyses of aftershock events [3] have established that the 2001 Bhuj mainshock might have been triggered by a deep-seated hidden fault at the base of a paleo-rift zone by reverse faulting.
EPSL 6696 2-7-03 Cyaan Magenta Geel Zwart
O.P. Mishra, D. Zhao / Earth and Planetary Science Letters 212 (2003) 393^405
399
Fig. 4. Epicenter locations of the mainshock and the best located aftershocks showing two trends [3]; the arrows indicate the rupture directions. B-L: Bachau Lineament; A-R-L: Anjar Rapar Lineament. The mainshock is shown by the solid star, and the two largest aftershocks are shown by open stars. The isoseismal IX is shown by curved line (after Shanker and Pande [46]). The fault plane solutions for mainshock (M.S) and for some selected aftershocks at selected depth ranges are shown in the equal area projections. The dark and white zones in focal mechanism solutions denote compression and dilatation, respectively (adapted from Kayal et al. [3]).
This supports our interpretation that the 2001 Bhuj mainshock hypocenter is associated with high-O, high-h and high-i, because the intersecting fault geometry might have rendered the paleorift zone as a suitable ‘£uids-trapper’ with increased £uid pressure, increased permeability, decreased seismic strength and subsequently high seismicity in the continental interior [29,30]. Tomographic images of O, h and i at 25 km depth, the source zone for the mainshock hypocenter, are shown in Fig. 6. The cross-sectional images of O, h and i along lines AB and CD (Fig. 1) are presented in Fig. 7. These results reveal strong dependence of these parameters on the seismogenic strengths and material properties of the in situ
rocks within the 2001 Bhuj mainshock hypocenter. The average estimated value of O is 0.1 except for the Bhuj mainshock hypocenter and tectonically active areas of the Kutch Rift Basin, where O approaches 0.2. This suggests that micro-cracks may exist in more localized areas, while £uids permeated through wider areas of the crust. The O assumes lower values from the surface to 15^ 20 km depth, and gradually acquires higher values from 20 km downward and becomes highest at the source depth of the 2001 Bhuj mainshock hypocenter (25 km). This suggests that the in situ rocks near the source zone of the mainshock may be associated with a fractured rock matrix. The Bhuj region is devoid from volcanic or geothermal
EPSL 6696 2-7-03 Cyaan Magenta Geel Zwart
400
O.P. Mishra, D. Zhao / Earth and Planetary Science Letters 212 (2003) 393^405
anomaly, and hence changes in Vp , Vs and c (Fig. 3) may be attributable mainly to cracks and £uids in the crust vis-a'-vis evidence from other geophysical observations (e.g. magnetotelluric sounding), as reported by [17]. Lithological variations may contribute to a small proportion of the changes in the estimated velocity, but a detailed investigation of site-relevant lithology in this area would be necessary in order to con¢rm this. A signi¢cant discrepancy between O and h exists in the present study at certain depth layers. The actual reason for this is not known, although it may be that some cracks whose aspect ratio is lower than 1033 may not be imaged in our 3-D distribution of the O. Experimental studies have shown that cracks with aspect ratio less than 1033 might be missed in imaging or analyses because the longest observed crack length of 565 Wm corresponds to the aspect ratio value of 5U1035 [28]. We therefore surmise that the distribution of micro-cracks at certain depths (Figs. 6 and 7) may have lower aspect ratio values than our assumed value (1033 ) in model OB74, so that we might have missed their imaging in the present study, and as a consequence observe a negative correlation between O and h.
5. Interpretation and discussion The seismic velocities are governed by temperature and lithological changes and thus the application of crack theory to the tomographic results is a matter for detailed investigation of these parameters [12]. Studies revealed that the Bhuj mainshock hypocenter is located in a distinctive zone characterized by high-Vp , low-Vs and high-c (Fig. 3) [4,5]. The actual cause of high-Vp within the source area of the Bhuj earthquake and its gradual variation into the lower crust is still unclear. It could be due to several reasons, such as the presence of an igneous intrusion, high compressibility and high pore pressures in the prevailing source rock. The interpretation of low-Vs with an increase of £uids is based on the fact that Vs is more sensitive to £uids than Vp [1,6,31,32]. The sensitivity of Vs to £uids is a function of pore pressure, crack density and saturation rate. Once
Fig. 5. Temporal distribution of (upper panel) the aftershocks recorded during January 30^April 15, 2001. Three lower panels: The best relocated events (368) plotted depthward, recorded in di¡erent months. Solid circles denote aftershock events and solid star denotes the 2001 Bhuj mainshock hypocenter.
the pores get saturated due to the presence of £uids, the pore pressure may continue to decrease. This causes high-c due to high-Vp /Vs ratio, in which Vs gets reduced strongly, while Vp remains una¡ected. The reverse is also true: that c or Vp /Vs ratio becomes low when rocks become undersaturated due to dilatancy, which will strongly reduce Vp , but will have little e¡ect on Vs [33,34]. High-O, high-h and high-i at the 2001 Bhuj mainshock hypocenter (Fig. 6) support our tomographic interpretation that the mainshock may have been initiated by a £uid-¢lled, fractured
EPSL 6696 2-7-03 Cyaan Magenta Geel Zwart
O.P. Mishra, D. Zhao / Earth and Planetary Science Letters 212 (2003) 393^405
401
Fig. 6. (a^c) Distribution of O, h, and i perturbations (in %) at 25 km depth. The depth range of aftershocks (open circles) are shown at the top of the map. Solid lines show active faults and lineaments with their abbreviations correspond to the same as that of Figs. 2a and 4. The solid star denotes the Bhuj mainshock. The perturbation scales are shown at the bottom of each image.
rock matrix, and hence a £uid-driven earthquake occurs within the stable continental region of the Indian Peninsula. However, the source of the fault zone £uids, and the level of in situ pressures remain unknown. These are important scienti¢c goals for future deep drilling at the source zone [15]. The anomalous zone exists in the depth range of 23^28 km and extends 15^30 km laterally, and is interpreted as a £uid-¢lled, fractured rock matrix that contributed to the initiation of the 2001 Bhuj earthquake. In volcanic areas, variations in seismic velocity may be caused mainly by variations in temperature. For the Bhuj region, we have analyzed all the available evidence (e.g. geothermal anomaly) and conclude that variations in seismic velocity and Poisson’s ratio are mainly due to the presence of £uid-¢lled micro-cracks [4,5]. This may be due
to regional tectonic faulting, rifting, or the e¡ects of magmatic intrusion and subsequent cooling at the time of Deccan volcanic eruption (65 Ma). Extensive fracturization, high permeability, and the lowering of the seismic strength of the rocks at the Bhuj mainshock hypocenter are expected due to intersecting fault geometry (Fig. 4) in the continental interior. The temporal and spatial distribution of aftershocks (Fig. 5) indicates that almost all selected events are concentrated in the rupture nucleated zones, which propagate bilaterally in NE and NW directions, with a prominent NE trend (Figs. 1 and 5). These NE and NW trends correspond to the Anjar Rapar Lineament (A-R-L) and Bachau Lineament (B-L), respectively (Figs. 4 and 6). It can be inferred from our results that the hypocentral locations are a function of time and therefore
EPSL 6696 2-7-03 Cyaan Magenta Geel Zwart
402
O.P. Mishra, D. Zhao / Earth and Planetary Science Letters 212 (2003) 393^405
Fig. 7. The same as Fig. 3 but for O, h and i distributions along cross-sectional lines AB (a^c) and CD (d^f). The corresponding perturbation scales are shown in the middle of the two panels (a^c) and (d^f). The star denotes the 2001 Bhuj mainshock hypocenter.
demonstrate the role of £uids, because seismogenic permeability is a characteristic value of fractured rocks, where seismicity is associated with increased pore pressure [35]. In order to substantiate our ¢ndings, analyses of hydrologic parameters (e.g. hydraulic di¡usivity and permeability)
may be a logical step in the near future. Albeit it is not uncommon for an earthquake to cause fracturing and an accompanying increase in fracture density and permeability, it is thereby perturbing the hydrologic conditions in the hypocentral region as a post-earthquake phenomenon, instead
EPSL 6696 2-7-03 Cyaan Magenta Geel Zwart
O.P. Mishra, D. Zhao / Earth and Planetary Science Letters 212 (2003) 393^405
of during the pre-earthquake period [36,37]. But we believe that our macroscopic estimates of O, h and i may be related to in situ rocks of the 2001 Bhuj hypocenter in which the accumulation of cracks may have been further reinforced by ongoing micro-seismic activity due to active tectonics [38]. A few studies have also been carried out to investigate the in situ status and distribution of cracks and £uids in active fault zones and earthquake source areas [1,24,39,40], and these have found that in situ material heterogeneities have a strong in£uence over earthquake initiation, propagation and termination processes. Rock porosities have also been related to i. Relative £uctuations in i may be interpreted as £uctuations in porosity. However, at a certain depth, it is believed that all pores and cracks are closed and porosity approaches zero, which may be the cause of the loss of correlation between the porosity parameter (i) and porosity [41]. Vertical cross-sections show a wider distribution of O, h and i at varying depths (Fig. 7). At the depth range of 10^15 km, where very few earthquakes occur, O and h show low values (Fig. 7b and e). High-h areas are generally consistent with high-O but their distribution is comparatively wide. For most of the study area, h is larger than 50% (0.5). The i shows a similar pattern of distribution in the earthquake region (Fig. 7c and f), the shallower layers are marked with average to low-i (high porosity), while the deeper layers possess average to high-i (low porosity). The evolution of porosity depends on the saturation state of pore £uids. Super-saturation leads to solute precipitation (cementation), and undersaturation leads to free surface dissolution [42]. Our result shows that the 2001 Bhuj hypocenter region (20^25 km) is associated with high-h, and the rocks of the hypocenter may have acquired excessive strain due to high pore pressure. This suggests that the e¡ect of porosity is reduced due to compaction and cementation at deeper levels (20^25 km), because high-h may have a profound e¡ect on porosity loss at high strain rate due to high e¡ective pore pressure (total pressure minus pore pressure of saturated and dry cracks), primarily through its control on cementation [42,43]. Thus the 2001 Bhuj hypocenter possesses
403
high-i (low porosity), despite its high-O and highh. These anomalies are expected to be regional in scale because of the occurrence of a large earthquake at the base of an ancient paleo-rift [3^5,30]. The aftershocks are less concentrated at shallower depths (0^10 km), and are mostly di¡used around the mainshock epicenter area (Figs. 1, 5 and 7). This probably suggests the presence of brittle layers at these depths, which might have started fracturing following the mainshock. The present study vindicates the tomographic ¢ndings that the Bhuj hypocenter is associated with a £uid-¢lled, fractured rock matrix and the prevailing rocks su¡ered from intense fracturization (high-O) with high-h and high-i. It is intriguing to note that the hypocentral depth (25 km) of the 2001 Bhuj earthquake is greater than most intraplate earthquakes ( 6 15 km) and those earthquakes occurring in regions of granitic lithology. The evidence for basaltic layer and the widespread presence of Deccan basalt is also suggested by travel time tomography [44], with a lack of a clear boundary between granitic and basaltic layers in the study region. However, our tomographic results are derived from the 1-D velocity model of [3], in which 20.5 km is considered to be a distinct boundary between granitic and basaltic layers. This indicates that the 2001 Bhuj mainshock may not have been initiated at the transition between these two lithologies (granitic and basaltic), but in the middle to relatively deep crust. Still we believe that lithological variability may lead to variations in the tomographic images, and corresponding di¡erences in the frictional properties would in£uence the rupture process. However, lateral variations in lithology in the middle to deep crust are unknown and di⁄cult to evaluate at present. This problem may properly be addressed by applying other tools, such as seismic attenuation and shear wave splitting, and siterelevant lithological analyses [13,45]. The velocity patterns revealed by seismic tomography [4] may suggest the in£uence of cracks and £uids, pending more de¢nitive tests to investigate whether other factors such as lithology have an insigni¢cant effect, but they further demonstrate that earthquake occurrence is closely related to in situ material heterogeneities. Our present results (Figs. 6 and
EPSL 6696 2-7-03 Cyaan Magenta Geel Zwart
404
O.P. Mishra, D. Zhao / Earth and Planetary Science Letters 212 (2003) 393^405
7) corroborate the velocity (high-Vp , low-Vs ), Poisson’s ratio (high-c) (Fig. 3) [4,5] and magnetotelluric models (conductive layer at a depth range of 10^17 km) [17]. This suggests that the 2001 Bhuj hypocenter is associated with a £uid¢lled, fractured rock matrix, and that the rocks existing at the mainshock hypocenter may have su¡ered from intense fracturization (high-O) with high-h and high-i.
6. Conclusions The 2001 Bhuj mainshock hypocenter is located in a distinctive zone characterized by high-O, highh, and high-i. This may indicate a £uid-¢lled, fractured rock matrix due to the presence of intersecting fault geometry near the 2001 Bhuj mainshock hypocenter, and can be taken as additional supportive evidence for a £uid-driven earthquake. Our tomographic results depict a well-oriented cluster of aftershock events in the vicinity of the Bhuj mainshock which were triggered by the fractured rock matrix. The presence of a £uid¢lled, fractured rock matrix, with intense fracturization, high permeability and lowered seismic strength in the source zone might have a¡ected the long-term structural and compositional evolution of the fault zone, its strength, and the local stress regime in the source area, accelerating the concentration of stress in the seismogenic layer. Consequently, this may have led to mechanical failure of asperities and barriers and thus might have contributed to the initiation of the 2001 Bhuj earthquake. It is inferred from this study that earthquake occurrence is closely related to in situ material heterogeneities, rather than stress conditions alone. Our results are comparable to those of the 1995 Kobe earthquake [1,12]. The use of magnetotelluric, gravity, and geodetic measurements along with hydrological observations, all coupled to seismic tomography, can yield a better image of the fault zones containing £uids. Future advances in deep drilling, seismic wave attenuation, and shear wave splitting may provide us authentic information on lithological heterogeneity in the crust, which can be used to better explain our tomographic results.
Acknowledgements We thank Scott D. King, Manika Prasad and three anonymous reviewers for their critical and useful comments for improvements of the manuscript, J.R. Kayal and A. Yamada for help in data processing and discussion, R.P. Rapp for proof reading the manuscript and providing us helpful comments. This work is partially supported by a grant (Kiban-B No. 11440134) from the Japan Society for the Promotion of Science to D.Z.[SK]
References [1] D. Zhao, H. Kanamori, H. Negeshi, Tomography of the source area of the 1995 Kobe earthquake: evidence for £uids at the hypocenter?, Science 274 (1996) 1891^ 1894. [2] K. Rajendran, C.P. Rajendran, M. Thakkar, M.P. Tuttle, The 2001 Kutch (Bhuj) earthquake: Coseismic surface features and their signi¢cance, Curr. Sci. 80 (2001) 1397^1405. [3] J.R. Kayal, R. De, S. Ram, B.V. Srirama, S.G. Gaonkar, Aftershocks of the January 26, 2001 Bhuj earthquake and its seismotectonic implications, J. Geol. Soc. India 59 (2002) 395^417. [4] J.R. Kayal, D. Zhao, O.P. Mishra, R. De, O.P. Singh, The 2001 Bhuj earthquake: Tomographic evidence for £uids at the hypocenter and its implications for rupture nucleation, Geophys. Res. Lett. 29 (2002) 10.1029/ 2002GL015177. [5] O.P. Mishra, D. Zhao, J.R. Kayal, The 2001 Bhuj earthquake: Seismic image of the source zone and its implications for the rupture nucleation, Phys. Earth Planet. Inter., submitted. [6] R.J. O’Connell, B. Budiansky, Seismic velocities in dry and saturated cracked solids, J. Geophys. Res. 79 (1974) 5412^5426. [7] D. Zhao, A. Hasegawa, S. Horiuchi, Tomographic imaging of P and S wave velocity structure beneath northeastern Japan, J. Geophys. Res. 97 (1992) 19909^19928. [8] C. Kisslinger, Aftershocks and fault-zone properties, Adv. Geophys. 38 (1996) 1^36. [9] J.M. Lees, Tomographic P-wave velocity images of the Loma Prieta earthquake asperity, J. Geophys. Res. 17 (1992) 1433^1436. [10] P.A. Johnson, T.V. McEvilly, Park¢eld seismicity: Fluiddriven?, J. Geophys. Res. 100 (1995) 12937^12950. [11] D. Eberhart-Phillips, A.J. Michael, Three-dimensional velocity structure, seismicity, and fault structure in the Park¢eld region, central California, J. Geophys. Res. 98 (1993) 15737^15758.
EPSL 6696 2-7-03 Cyaan Magenta Geel Zwart
O.P. Mishra, D. Zhao / Earth and Planetary Science Letters 212 (2003) 393^405 [12] D. Zhao, T. Mizuno, Crack density and saturation rate in the 1995 Kobe earthquake region, Geophys. Res. Lett. 26 (1999) 3213^3216. [13] D. Zhao, H. Negishi, The 1995 Kobe earthquake: seismic image of the source zone and its implications for the rupture nucleation, J. Geophys. Res. 103 (1998) 9967^ 9986. [14] D. Zhao, O.P. Mishra, R. Sanda, In£uence of £uids and magma on earthquakes: seismological evidence, Phys. Earth Planet. Inter. 132 (2002) 249^267. [15] E. Roelo¡s, The Park¢eld, California earthquake experiments: An update in 2000, Curr. Sci. 79 (2000) 1226^ 1236. [16] S.A. Miller, Fluid-mediated in£uence of adjacent thrusting on the seismic cycle at Park¢eld, Nature 382 (1996) 799^802. [17] B.R. Arora, G. Rawat, A.K. Singh, Mid-crustal conductor below the Kutch Rift Basin and its seismogenic relevance to the 2001 Bhuj earthquake, DCS-DST News, Department of Science and Technology, Govt. of India, New Delhi, August 2002. [18] H.K. Gupta, S.V. Sarma, T. Harinarayana, G. Virupakshi, Fluids below the hypocentral region of Latur earthquake, India: Geophysical indicators, Geophys. Res. Lett. 23 (1996) 1569^1572. [19] S.K. Biswas, Regional tectonic framework, structure and evolution of the western marginal basins of India, Tectonophysics 135 (1987) 307^327. [20] A.C. Johnston, L.R. Canter, Earthquakes in stable continental crust, Sci. Am. 262 (1990) 68^75. [21] P. Talwani, A. Gangopadhyay, Tectonic framework of the Kachchh earthquake of 26 January 2001, Seism. Res. Lett. 72 (2001) 336^345. [22] J. Krishna, I.B. Singh, J.D. Howard, S.A. Jafar, Implications of new data on Mesozoic rocks of Kuchchh peninsula, western India, Nature 305 (1983) 790^792. [23] T. Utsu, Seismology, Kyoritsu Press, Tokyo, 1984, 310 pp. [24] J.M. Lees, H. Wu, Poisson’s ratio and porosity at Coso geothermal area, California, J. Volcanol. Geotherm. Res. 95 (2000) 157^173. [25] G.R. Pickett, Acoustic character logs and their applications in formation evaluation, J. Petrol. Tech. 15 (1963) 650^667. [26] Y. Takei, E¡ects of pore geometry on Vp/Vs: from equilibrium geometry to crack, J. Geophys. Res. 107 (2002) 10.1029/2001JB000522. [27] K. Yamamoto, M. Kosuga, T. Hirasawa, A theoretical method for determination of e¡ective elastic constants of isotropic composites, Tohoku Geophys. J. 28 (1981) 47^67. [28] K. Hadley, Comparison of calculated and observed crack densities and seismic velocities in westerly granite, J. Geophys. Res. 81 (1976) 3484^3493. [29] P. Talwani, Fault geometry and earthquakes in continental interiors, Tectonophysics 305 (1999) 371^379.
405
[30] S. Marshak, T. Paulsen, Structural style, regional distribution, and seismic implications of Midcontinental faultand-fold zones, United States, Seismol. Res. Lett. 68 (1997) 511^520. [31] M.N. Toksoz, C.H. Cheng, A. Timur, Velocities of seismic waves in porous rocks, Geophysics 41 (1976) 621^ 645. [32] D. Moos, M.D. Zoback, In situ studies of velocity in the fractured crystalline rocks, J. Geophys. Res. 88 (1983) 2345^2358. [33] C.H. Scholz, L.R. Sykes, Y.P. Aggarwal, Earthquake prediction: A physical basis, Science 181 (1973) 803^809. [34] M. Zimmer, M. Prasad, G. Mavko, Pressure and porosity in£uences on Vp^Vs ratio in unconsolidated sands, Leading Edge 21 (2002). [35] P. Talwani, Seismogenic properties of the crust inferred from recent studies of reservoir-induced seismicity ^ Application to Koyana, Curr. Sci. 79 (2000) 1327^1333. [36] S. Rojstaczer, Permeability changes associated with large earthquakes: An example from Loma Priets, California, Geology 20 (1992) 211^214. [37] S. Rojstaczer, S. Wolf, R. Michael, Permeability enhancement in the shallow crust as a cause of earthquake-induced hydrological changes, Nature 373 (1995) 237^ 239. [38] P. Malign, The seismology of extensional hydrothermal systems, Geotherm. Resour. Counc. Trans. 18 (1994) 17^22. [39] G. Vlahovic, C.A. Powell, Three-dimensional S wave velocity structure and Vp/Vs ratios in the New Madrid Seismic Zone, J. Geophys. Res. 106 (2001) 13501^ 13513. [40] I. Serrano, D. Zhao, J. Morale, 3-D crustal structure of the extensional Granada basin in the convergent boundary between the Eurasia and African plates, Tectonophysics 344 (2002) 61^79. [41] W.F. Brace, E. Silver, K. Hadley, C. Goyze, Cracks and pores: A closer look, Science 178 (1972) 162^164. [42] W. He, A. Hajash, D. Sparks, A model for porosity evolution during creep compaction of sand stones, Earth Planet. Sci. Lett. 197 (2002) 237^244. [43] M. Prasad, Correlating permeability with velocity using £ow zone indicators, SEG Expanded Abstracts, 69th Annual SEG Meeting 587, 1999. [44] B.L.N. Kennet, S. Widiyantoro, A low seismic wavespeed anomaly beneath northwestern India: a seismic signature of the Deccan plume, Earth Planet. Sci. Lett. 165 (1999) 145^155. [45] M. Prasad, M.H. Manghanani, E¡ects of pore and di¡erential pressure on compressional wave velocity and quality factor in Bera and Michigan sandstones, Geophysics 62 (1997) 1163^1176. [46] R. Shanker, P. Pande, Geoseismological studies of Kutch (Bhuj) earthquake of 26 January, 2001, J. Geol. Soc. India 58 (2001) 203^208.
EPSL 6696 2-7-03 Cyaan Magenta Geel Zwart
Earth and Planetary Science Letters 212 (2003) 393^405 www.elsevier.com/locate/epsl
Crack density, saturation rate and porosity at the 2001 Bhuj, India, earthquake hypocenter: a £uid-driven earthquake? O.P. Mishra , Dapeng Zhao Geodynamics Research Center, Ehime University, Matsuyama 790-8577, Japan Received 10 December 2002; received in revised form 12 May 2003; accepted 14 May 2003
Abstract Quantitative estimates of the crack density (O), saturation rate (h) and porosity (i) parameters from seismic velocities (Vp , Vs ) and Poisson’s ratio (c) in the 2001 Bhuj earthquake area in western India indicate that the 2001 Bhuj earthquake hypocenter is associated with high-O, high-h and high-i in the depth range of about 23^28 km, extending 15^30 km laterally. These anomalies may be due to a fluid-filled, fractured rock matrix, which might have contributed to trigger the 2001 Bhuj earthquake in the intraplate, stable continental region of the Indian Peninsula. This feature is similar to that of the 1995 Kobe earthquake [D. Zhao et al., Science 274 (1996) 1891^1894]. High-i areas are generally consistent with high-O areas, but high-h areas have a wider distribution, indicating that microcracks exist in more localized areas within the crust, and that permeation of fluids in the hypocenter zone might have occurred extensively through the intergranular and fractured pores due to hidden intersecting fault geometry. Here we suggest the possibility that earthquake occurrence is closely related to in situ material heterogeneities, rather than stress conditions alone. @ 2003 Elsevier Science B.V. All rights reserved. Keywords: crack density; saturation rate; porosity; £uids; crustal earthquakes
1. Introduction The January 26, 2001 Bhuj earthquake (Mw 7.6) in western India has been investigated by several researchers using multidisciplinary approaches, making it one of the best studied Indian earthquakes in the annals of seismology. It was one of the most catastrophic earthquakes in In-
* Corresponding author. Tel.: +81-89-927-8257; Fax: +81-89-927-8167. E-mail addresses: [email protected] (O.P. Mishra), [email protected] (D. Zhao).
dia, which rocked Gujarat province at Bhuj (23.4‡N, 70.3‡E) in western India, with an estimated death toll of over 20 000 [2]. The 2001 Bhuj mainshock was followed by a series of aftershock events (over 4000), magnitude v 1.0, which were recorded during the period from January 30 to April 15, 2001 [3]. Arrival times generated from 368 selected aftershocks have been used (Fig. 1) for tomographic estimations of three-dimensional (3-D) P-wave (Vp ) and S-wave (Vs ) velocity and Poisson’s ratio (c) structures in the hypocentral area of the Bhuj earthquake to understand what may have triggered the Bhuj mainshock [4,5]. To better understand the relationship between cracks/
0012-821X / 03 / $ ^ see front matter @ 2003 Elsevier Science B.V. All rights reserved. doi:10.1016/S0012-821X(03)00285-1
EPSL 6696 2-7-03 Cyaan Magenta Geel Zwart
394
O.P. Mishra, D. Zhao / Earth and Planetary Science Letters 212 (2003) 393^405
£uids and the earthquake nucleation processes, we have attempted to estimate the 3-D distribution of crack density (O), saturation rate (h) and porosity (i) parameters in the source volume of the 2001 Bhuj earthquake, by applying the crack theory [6] to the values of Vp , Vs and c derived from a high resolution, 3-D tomographic inversion [7]. The spatial and temporal distribution of the events in an aftershock sequence contains information about the aftershock-generating process in particular, and by extrapolation, the earthquake-generating process in general, because the distribution of aftershocks is closely linked to the geometry of the fault zone that ruptured to produce the mainshock, which in turn depends on the physical properties of the fault zone as well as ambient conditions, especially the distribution of strength, stress, and temperature [8]. In this study, aftershock events are used as a source for imaging the subsurface structure from a particular seismic network design (Fig. 1). Aftershock events have previously been used by other researchers to image the 3-D structure in the source zone of di¡erent earthquakes elsewhere in the world to understand what may have triggered these earthquakes, and how the rupture developed after initiation [9^ 13]. Complex crustal structure would have certainly in£uenced the initiation, propagation and termination of the earthquake rupture. If variations in material properties in£uence the seismic velocity, then we should be able to image them in a velocity model [13,14]. Fault zone £uids are widely believed as a key to understanding earthquake generation. Speci¢cally, it has been hypothesized that high £uid pressure in the fault zone acts to reduce the frictional strength of the fault zone, and that time variations in £uid pressure control the timing of earthquakes [15,16]. Our study shows that £uid involvement in the fault zone may be a causative agent in triggering large earthquakes as a consequence of unique £uid pressure, rather than lithologic conditions [10,12^15]. Evidence to date for £uid-driven seismicity for the 2001 Bhuj earthquake includes appreciable seismic velocity (Vp , Vs ) and Poisson’s ratio (c) contrasts in the source zone [4,5], which are well corroborated with the presence of conductive crustal structure as inferred from magne-
totelluric measurements [17]. Evidence for the presence of £uids in other seismogenic fault zones elsewhere in the world is also reported [9^12,18]. 3-D estimates of crack density, saturation rate and porosity parameters in this study may provide additional information that will add to our understanding of the rupture initiation process for the 2001 Bhuj earthquake besides a supportive evidence for £uid-driven earthquake.
2. Tectonic setting and geology The Bhuj mainshock occurred in the Kutch Rift Basin. The Kutch region is underlain by a Mesozoic rift system associated with three major basins (Kutch, Combay, Narmada) in the western margin of the peninsular India craton (Fig. 2a). Mesozoic sediments were uplifted, folded, intruded and covered by Deccan Trap basaltic £ows in Late Cretaceous and early Paleocene time [19]. Faults within such rift systems are known to have the potential to generate large earthquakes [20], and the area has witnessed several large historical earthquakes in the past [21]. The source region of the 2001 Bhuj earthquake is ¢lled with sediments ranging in age from Middle Jurassic to present. A recent study [22] suggests a marine origin for the Jurassic^Cretaceous unit of the Kutch Rift Basin, in which Eocene and older sediments are distributed in highlands and lowlands as shown in Fig. 2b.
3. Data and method Fig. 1 shows the epicentral distribution of 368 Bhuj aftershocks, which were selected from a total of over 4000 aftershocks (M v 1.0). These 368 events generated about 4000 phase data that consist of 1960 P- and 1880 S-wave arrivals, which were used in tomographic imaging. The reliability of P- and S-wave arrival time data were carefully checked using digital and analog recorders [4]. Our reading accuracy varies from 0.01 to 0.1 s for P and 0.05 to 0.3 s for S arrivals. The majority of the events are located by more than four seismic stations from a seismic network consisting of
EPSL 6696 2-7-03 Cyaan Magenta Geel Zwart
O.P. Mishra, D. Zhao / Earth and Planetary Science Letters 212 (2003) 393^405
395
Fig. 1. Epicentral distribution of the 368 Bhuj aftershocks used in this study. The blue circles (aftershocks); the red star (Bhuj mainshock). Lines AB and CD show the locations of the cross-sections in Figs. 3a^c and 7a^c and Fig. 3d^f and 7d^f, respectively. The red triangles denote portable seismic stations. The variation of size of blue circles indicates the variation of magnitude.
12 stations (Fig. 1). The average accuracy of hypocentral locations is better than 4 km. Depth sections of the events show that most of the aftershocks occurred at a depth range of 20^40 km (Fig. 1). We have considered 25 km as the accurate source depth of the 2001 Bhuj mainshock as determined by using converted depth phase data [3^5]. We conducted detailed checkerboard resolution tests [7] for both Vp and Vs at di¡erent grid spac-
ing to assess the resolution of our data set. We found that the present data set is able to image O, h and i (as derived from Vp , Vs and c) structures with sizes of 30 km horizontal and 5 km vertical near the Bhuj mainshock hypocenter [5]. The crack theory of O’Connell and Budiansky [6] (hereafter we call it OB74) is considered as a concise, self-consistent, isotropic theory, capable of estimating crack density and saturation rate in stressed rock from seismic velocity data [12].
EPSL 6696 2-7-03 Cyaan Magenta Geel Zwart
396
O.P. Mishra, D. Zhao / Earth and Planetary Science Letters 212 (2003) 393^405
EPSL 6696 2-7-03 Cyaan Magenta Geel Zwart
O.P. Mishra, D. Zhao / Earth and Planetary Science Letters 212 (2003) 393^405
A set of equations relating changes in Vp and Vs to c [23], and also changes in Vp , Vs and c to O and h are used to estimate O and h. The rock parameters are derived from OB74 on the assumption that £at-circular cracks are aligned randomly in the crust and are ¢lled with £uids or air [6,12]. Changes in Vp and Vs to c can be expressed as: ðV p =V s Þ2 ¼ 2ð13c Þ=ð132c Þ:
O ¼ NGa f;
study region necessarily behave like sedimentary rocks, but rather that the bulk properties of the highly fractured, £uid-saturated regions probably behave like a porous medium. The porosity theory of P63 states that i is directly proportional to the product of the seismic velocities (Vp and Vs ), which can be expressed as: i ¼ V p UV s ;
ð6Þ
ln i ¼ ln V p þ ln V s ;
ð7Þ
d i = i ¼ dV p =V p þ dV s =V s ;
ð8Þ
ð1Þ
According to OB74, the crack density parameter O can be de¢ned as: 3
397
ð2Þ
where N is the number of micro-cracks per unit volume and Gaf is the average radius of circular cracks. If it is assumed that N1 cracks per unit volume are dry and that N2 = N3N1 cracks are saturated, then the saturation rate h is de¢ned as N2 /N. OB74 predicts the following relations: V p2 13c 1 þ c 16 13c 2 ¼ W 13 ð13 h ÞO ð3Þ V 2p 1 þ c 13c 9 132c V s2 32 3 ð13 ¼ 13 c Þ 13 h þ O ð4Þ V 2s 45 23c 45 c 3c 23c ; ð5Þ O¼ 16 13c 2 ð13 h Þð1 þ 3c Þð23c Þ32ð132c Þ
where V p , V s , c are P-wave and S-wave velocities and Poisson’s ratio for cracked volume, respectively, while Vp , Vs , c are the corresponding parameters for uncracked rocks. Thus O and h can be estimated from observed Vp , Vs and c using Eqs. 3^5. The method of [24] based on porosity theory of Pickett [25] (hereafter we call it P63) is used to estimate the porosity parameter (i). The P63 model does not assume that the rocks in the
where di/i is the change of porosity per unit porosity, while dVp /Vp and dVs /Vs are perturbed values of Vp and Vs , respectively. Eq. 8 is used to estimate the distribution of porosity in the 2001 Bhuj earthquake region. The perturbed values of Vp and Vs with respect to an average initial velocity model at di¡erent depth layers (0^30 km) were obtained after setting up a 3-D grid net in the study area [4,5]. We have estimated O, h and i at the same depth layers using the corresponding perturbed values of Vp , Vs and c (Fig. 3). It is imperative to note that high-i denotes low porosity, and the reverse is also true, because i bears a direct relationship to porosity as the product of seismic velocities [25]. Some new and improved models have subsequently been proposed to estimate O and h, due to inapplicability of OB74 for cases involving large crack density and di¡erent pore geometry [26,27]. However, our tomographic results show that the change of seismic velocities is less than 5^6% [4,5], which leads to O 6 0.25, and a corresponding average i perturbation in the range of U 10%. Hence the OB74 and P63 models are viable for this study.
6 Fig. 2. (a) Map showing tectonic features of the study region (compiled from Kayal et al. [3] and Talwani and Gangopadhyay [21]). KRB: Kutch Rift Basin; CRB: Cambay Rift Basin; NRB: Narmada Rift Basin; NPF: Nagar Parkar Fault; IBF: Island Belt Fault; KMF: Kutch Mainland Fault; NKF: North Kathiwar Fault; WCF: West Coast Fault; AB: ‘Allah Bund’ (see text). Epicenters (solid circles) and fault plane solutions of the past signi¢cant earthquakes (M s 5.0) are shown. For fault plane solutions, the solid circle indicates P-axis in the dilatational zone, and the open circle T-axis in the compressional zone. The star indicates the mainshock epicenter of the January 26, 2001 Bhuj earthquake. The insert map shows the major plate boundaries around India and the present study area (adapted from Kayal et al. [4]). (b) Map showing geology of the area (adapted from Talwani and Gangopadhyay [21] after Krishna et al. [22]).
EPSL 6696 2-7-03 Cyaan Magenta Geel Zwart
398
O.P. Mishra, D. Zhao / Earth and Planetary Science Letters 212 (2003) 393^405
Fig. 3. The upper panel (a^c) represents vertical cross-sections along the line AB for Vp , Vs and c, respectively. The lower panel (d^f) same as that of (a^c) but along the line CD. The locations of lines AB and CD are shown in Fig. 1. Small white circles denote the Bhuj aftershocks within 9 km width along the lines AB and CD. Star denotes the hypocenter of Bhuj mainshock; its focal depth is 25.0 km. The perturbation (in %) scales for velocity and Poisson’s ratio are shown at the bottom.
4. Analysis and results We assumed in the theoretical formulations of OB74 that uncracked rocks have the maximum seismic velocity at each depth in the study area. The e¡ects of porosity have been analyzed separately in order to understand the in situ status of pores around the 2001 Bhuj mainshock hypocenter. Average crack porosity for stressed rock is taken as a function of crack orientation and aspect ratio [28]. The distribution of Vp , Vs and c in depths (Fig. 3) along lines AB and CD (Fig. 1) are used to estimate O, h and i. The intersecting fault geometry (Fig. 4) shows bilateral and simultaneous ruptures that propagated along two fault
trends (NE and NW) where the largest (ML 5.7) and second largest (ML 5.5) aftershocks occurred along the NE and NW rupture nucleated trends, respectively [3]. The temporal variation of all aftershocks, recorded by a portable seismic network (Fig. 1) is shown in Fig. 5, upper part, while Fig. 5, three lower parts, show the temporal variations of 368 relocated aftershocks along various sections, which are used in our tomographic study [5]. This ¢gure shows a well-oriented clustering of aftershocks around the mainshock hypocenter. Detailed analyses of aftershock events [3] have established that the 2001 Bhuj mainshock might have been triggered by a deep-seated hidden fault at the base of a paleo-rift zone by reverse faulting.
EPSL 6696 2-7-03 Cyaan Magenta Geel Zwart
O.P. Mishra, D. Zhao / Earth and Planetary Science Letters 212 (2003) 393^405
399
Fig. 4. Epicenter locations of the mainshock and the best located aftershocks showing two trends [3]; the arrows indicate the rupture directions. B-L: Bachau Lineament; A-R-L: Anjar Rapar Lineament. The mainshock is shown by the solid star, and the two largest aftershocks are shown by open stars. The isoseismal IX is shown by curved line (after Shanker and Pande [46]). The fault plane solutions for mainshock (M.S) and for some selected aftershocks at selected depth ranges are shown in the equal area projections. The dark and white zones in focal mechanism solutions denote compression and dilatation, respectively (adapted from Kayal et al. [3]).
This supports our interpretation that the 2001 Bhuj mainshock hypocenter is associated with high-O, high-h and high-i, because the intersecting fault geometry might have rendered the paleorift zone as a suitable ‘£uids-trapper’ with increased £uid pressure, increased permeability, decreased seismic strength and subsequently high seismicity in the continental interior [29,30]. Tomographic images of O, h and i at 25 km depth, the source zone for the mainshock hypocenter, are shown in Fig. 6. The cross-sectional images of O, h and i along lines AB and CD (Fig. 1) are presented in Fig. 7. These results reveal strong dependence of these parameters on the seismogenic strengths and material properties of the in situ
rocks within the 2001 Bhuj mainshock hypocenter. The average estimated value of O is 0.1 except for the Bhuj mainshock hypocenter and tectonically active areas of the Kutch Rift Basin, where O approaches 0.2. This suggests that micro-cracks may exist in more localized areas, while £uids permeated through wider areas of the crust. The O assumes lower values from the surface to 15^ 20 km depth, and gradually acquires higher values from 20 km downward and becomes highest at the source depth of the 2001 Bhuj mainshock hypocenter (25 km). This suggests that the in situ rocks near the source zone of the mainshock may be associated with a fractured rock matrix. The Bhuj region is devoid from volcanic or geothermal
EPSL 6696 2-7-03 Cyaan Magenta Geel Zwart
400
O.P. Mishra, D. Zhao / Earth and Planetary Science Letters 212 (2003) 393^405
anomaly, and hence changes in Vp , Vs and c (Fig. 3) may be attributable mainly to cracks and £uids in the crust vis-a'-vis evidence from other geophysical observations (e.g. magnetotelluric sounding), as reported by [17]. Lithological variations may contribute to a small proportion of the changes in the estimated velocity, but a detailed investigation of site-relevant lithology in this area would be necessary in order to con¢rm this. A signi¢cant discrepancy between O and h exists in the present study at certain depth layers. The actual reason for this is not known, although it may be that some cracks whose aspect ratio is lower than 1033 may not be imaged in our 3-D distribution of the O. Experimental studies have shown that cracks with aspect ratio less than 1033 might be missed in imaging or analyses because the longest observed crack length of 565 Wm corresponds to the aspect ratio value of 5U1035 [28]. We therefore surmise that the distribution of micro-cracks at certain depths (Figs. 6 and 7) may have lower aspect ratio values than our assumed value (1033 ) in model OB74, so that we might have missed their imaging in the present study, and as a consequence observe a negative correlation between O and h.
5. Interpretation and discussion The seismic velocities are governed by temperature and lithological changes and thus the application of crack theory to the tomographic results is a matter for detailed investigation of these parameters [12]. Studies revealed that the Bhuj mainshock hypocenter is located in a distinctive zone characterized by high-Vp , low-Vs and high-c (Fig. 3) [4,5]. The actual cause of high-Vp within the source area of the Bhuj earthquake and its gradual variation into the lower crust is still unclear. It could be due to several reasons, such as the presence of an igneous intrusion, high compressibility and high pore pressures in the prevailing source rock. The interpretation of low-Vs with an increase of £uids is based on the fact that Vs is more sensitive to £uids than Vp [1,6,31,32]. The sensitivity of Vs to £uids is a function of pore pressure, crack density and saturation rate. Once
Fig. 5. Temporal distribution of (upper panel) the aftershocks recorded during January 30^April 15, 2001. Three lower panels: The best relocated events (368) plotted depthward, recorded in di¡erent months. Solid circles denote aftershock events and solid star denotes the 2001 Bhuj mainshock hypocenter.
the pores get saturated due to the presence of £uids, the pore pressure may continue to decrease. This causes high-c due to high-Vp /Vs ratio, in which Vs gets reduced strongly, while Vp remains una¡ected. The reverse is also true: that c or Vp /Vs ratio becomes low when rocks become undersaturated due to dilatancy, which will strongly reduce Vp , but will have little e¡ect on Vs [33,34]. High-O, high-h and high-i at the 2001 Bhuj mainshock hypocenter (Fig. 6) support our tomographic interpretation that the mainshock may have been initiated by a £uid-¢lled, fractured
EPSL 6696 2-7-03 Cyaan Magenta Geel Zwart
O.P. Mishra, D. Zhao / Earth and Planetary Science Letters 212 (2003) 393^405
401
Fig. 6. (a^c) Distribution of O, h, and i perturbations (in %) at 25 km depth. The depth range of aftershocks (open circles) are shown at the top of the map. Solid lines show active faults and lineaments with their abbreviations correspond to the same as that of Figs. 2a and 4. The solid star denotes the Bhuj mainshock. The perturbation scales are shown at the bottom of each image.
rock matrix, and hence a £uid-driven earthquake occurs within the stable continental region of the Indian Peninsula. However, the source of the fault zone £uids, and the level of in situ pressures remain unknown. These are important scienti¢c goals for future deep drilling at the source zone [15]. The anomalous zone exists in the depth range of 23^28 km and extends 15^30 km laterally, and is interpreted as a £uid-¢lled, fractured rock matrix that contributed to the initiation of the 2001 Bhuj earthquake. In volcanic areas, variations in seismic velocity may be caused mainly by variations in temperature. For the Bhuj region, we have analyzed all the available evidence (e.g. geothermal anomaly) and conclude that variations in seismic velocity and Poisson’s ratio are mainly due to the presence of £uid-¢lled micro-cracks [4,5]. This may be due
to regional tectonic faulting, rifting, or the e¡ects of magmatic intrusion and subsequent cooling at the time of Deccan volcanic eruption (65 Ma). Extensive fracturization, high permeability, and the lowering of the seismic strength of the rocks at the Bhuj mainshock hypocenter are expected due to intersecting fault geometry (Fig. 4) in the continental interior. The temporal and spatial distribution of aftershocks (Fig. 5) indicates that almost all selected events are concentrated in the rupture nucleated zones, which propagate bilaterally in NE and NW directions, with a prominent NE trend (Figs. 1 and 5). These NE and NW trends correspond to the Anjar Rapar Lineament (A-R-L) and Bachau Lineament (B-L), respectively (Figs. 4 and 6). It can be inferred from our results that the hypocentral locations are a function of time and therefore
EPSL 6696 2-7-03 Cyaan Magenta Geel Zwart
402
O.P. Mishra, D. Zhao / Earth and Planetary Science Letters 212 (2003) 393^405
Fig. 7. The same as Fig. 3 but for O, h and i distributions along cross-sectional lines AB (a^c) and CD (d^f). The corresponding perturbation scales are shown in the middle of the two panels (a^c) and (d^f). The star denotes the 2001 Bhuj mainshock hypocenter.
demonstrate the role of £uids, because seismogenic permeability is a characteristic value of fractured rocks, where seismicity is associated with increased pore pressure [35]. In order to substantiate our ¢ndings, analyses of hydrologic parameters (e.g. hydraulic di¡usivity and permeability)
may be a logical step in the near future. Albeit it is not uncommon for an earthquake to cause fracturing and an accompanying increase in fracture density and permeability, it is thereby perturbing the hydrologic conditions in the hypocentral region as a post-earthquake phenomenon, instead
EPSL 6696 2-7-03 Cyaan Magenta Geel Zwart
O.P. Mishra, D. Zhao / Earth and Planetary Science Letters 212 (2003) 393^405
of during the pre-earthquake period [36,37]. But we believe that our macroscopic estimates of O, h and i may be related to in situ rocks of the 2001 Bhuj hypocenter in which the accumulation of cracks may have been further reinforced by ongoing micro-seismic activity due to active tectonics [38]. A few studies have also been carried out to investigate the in situ status and distribution of cracks and £uids in active fault zones and earthquake source areas [1,24,39,40], and these have found that in situ material heterogeneities have a strong in£uence over earthquake initiation, propagation and termination processes. Rock porosities have also been related to i. Relative £uctuations in i may be interpreted as £uctuations in porosity. However, at a certain depth, it is believed that all pores and cracks are closed and porosity approaches zero, which may be the cause of the loss of correlation between the porosity parameter (i) and porosity [41]. Vertical cross-sections show a wider distribution of O, h and i at varying depths (Fig. 7). At the depth range of 10^15 km, where very few earthquakes occur, O and h show low values (Fig. 7b and e). High-h areas are generally consistent with high-O but their distribution is comparatively wide. For most of the study area, h is larger than 50% (0.5). The i shows a similar pattern of distribution in the earthquake region (Fig. 7c and f), the shallower layers are marked with average to low-i (high porosity), while the deeper layers possess average to high-i (low porosity). The evolution of porosity depends on the saturation state of pore £uids. Super-saturation leads to solute precipitation (cementation), and undersaturation leads to free surface dissolution [42]. Our result shows that the 2001 Bhuj hypocenter region (20^25 km) is associated with high-h, and the rocks of the hypocenter may have acquired excessive strain due to high pore pressure. This suggests that the e¡ect of porosity is reduced due to compaction and cementation at deeper levels (20^25 km), because high-h may have a profound e¡ect on porosity loss at high strain rate due to high e¡ective pore pressure (total pressure minus pore pressure of saturated and dry cracks), primarily through its control on cementation [42,43]. Thus the 2001 Bhuj hypocenter possesses
403
high-i (low porosity), despite its high-O and highh. These anomalies are expected to be regional in scale because of the occurrence of a large earthquake at the base of an ancient paleo-rift [3^5,30]. The aftershocks are less concentrated at shallower depths (0^10 km), and are mostly di¡used around the mainshock epicenter area (Figs. 1, 5 and 7). This probably suggests the presence of brittle layers at these depths, which might have started fracturing following the mainshock. The present study vindicates the tomographic ¢ndings that the Bhuj hypocenter is associated with a £uid-¢lled, fractured rock matrix and the prevailing rocks su¡ered from intense fracturization (high-O) with high-h and high-i. It is intriguing to note that the hypocentral depth (25 km) of the 2001 Bhuj earthquake is greater than most intraplate earthquakes ( 6 15 km) and those earthquakes occurring in regions of granitic lithology. The evidence for basaltic layer and the widespread presence of Deccan basalt is also suggested by travel time tomography [44], with a lack of a clear boundary between granitic and basaltic layers in the study region. However, our tomographic results are derived from the 1-D velocity model of [3], in which 20.5 km is considered to be a distinct boundary between granitic and basaltic layers. This indicates that the 2001 Bhuj mainshock may not have been initiated at the transition between these two lithologies (granitic and basaltic), but in the middle to relatively deep crust. Still we believe that lithological variability may lead to variations in the tomographic images, and corresponding di¡erences in the frictional properties would in£uence the rupture process. However, lateral variations in lithology in the middle to deep crust are unknown and di⁄cult to evaluate at present. This problem may properly be addressed by applying other tools, such as seismic attenuation and shear wave splitting, and siterelevant lithological analyses [13,45]. The velocity patterns revealed by seismic tomography [4] may suggest the in£uence of cracks and £uids, pending more de¢nitive tests to investigate whether other factors such as lithology have an insigni¢cant effect, but they further demonstrate that earthquake occurrence is closely related to in situ material heterogeneities. Our present results (Figs. 6 and
EPSL 6696 2-7-03 Cyaan Magenta Geel Zwart
404
O.P. Mishra, D. Zhao / Earth and Planetary Science Letters 212 (2003) 393^405
7) corroborate the velocity (high-Vp , low-Vs ), Poisson’s ratio (high-c) (Fig. 3) [4,5] and magnetotelluric models (conductive layer at a depth range of 10^17 km) [17]. This suggests that the 2001 Bhuj hypocenter is associated with a £uid¢lled, fractured rock matrix, and that the rocks existing at the mainshock hypocenter may have su¡ered from intense fracturization (high-O) with high-h and high-i.
6. Conclusions The 2001 Bhuj mainshock hypocenter is located in a distinctive zone characterized by high-O, highh, and high-i. This may indicate a £uid-¢lled, fractured rock matrix due to the presence of intersecting fault geometry near the 2001 Bhuj mainshock hypocenter, and can be taken as additional supportive evidence for a £uid-driven earthquake. Our tomographic results depict a well-oriented cluster of aftershock events in the vicinity of the Bhuj mainshock which were triggered by the fractured rock matrix. The presence of a £uid¢lled, fractured rock matrix, with intense fracturization, high permeability and lowered seismic strength in the source zone might have a¡ected the long-term structural and compositional evolution of the fault zone, its strength, and the local stress regime in the source area, accelerating the concentration of stress in the seismogenic layer. Consequently, this may have led to mechanical failure of asperities and barriers and thus might have contributed to the initiation of the 2001 Bhuj earthquake. It is inferred from this study that earthquake occurrence is closely related to in situ material heterogeneities, rather than stress conditions alone. Our results are comparable to those of the 1995 Kobe earthquake [1,12]. The use of magnetotelluric, gravity, and geodetic measurements along with hydrological observations, all coupled to seismic tomography, can yield a better image of the fault zones containing £uids. Future advances in deep drilling, seismic wave attenuation, and shear wave splitting may provide us authentic information on lithological heterogeneity in the crust, which can be used to better explain our tomographic results.
Acknowledgements We thank Scott D. King, Manika Prasad and three anonymous reviewers for their critical and useful comments for improvements of the manuscript, J.R. Kayal and A. Yamada for help in data processing and discussion, R.P. Rapp for proof reading the manuscript and providing us helpful comments. This work is partially supported by a grant (Kiban-B No. 11440134) from the Japan Society for the Promotion of Science to D.Z.[SK]
References [1] D. Zhao, H. Kanamori, H. Negeshi, Tomography of the source area of the 1995 Kobe earthquake: evidence for £uids at the hypocenter?, Science 274 (1996) 1891^ 1894. [2] K. Rajendran, C.P. Rajendran, M. Thakkar, M.P. Tuttle, The 2001 Kutch (Bhuj) earthquake: Coseismic surface features and their signi¢cance, Curr. Sci. 80 (2001) 1397^1405. [3] J.R. Kayal, R. De, S. Ram, B.V. Srirama, S.G. Gaonkar, Aftershocks of the January 26, 2001 Bhuj earthquake and its seismotectonic implications, J. Geol. Soc. India 59 (2002) 395^417. [4] J.R. Kayal, D. Zhao, O.P. Mishra, R. De, O.P. Singh, The 2001 Bhuj earthquake: Tomographic evidence for £uids at the hypocenter and its implications for rupture nucleation, Geophys. Res. Lett. 29 (2002) 10.1029/ 2002GL015177. [5] O.P. Mishra, D. Zhao, J.R. Kayal, The 2001 Bhuj earthquake: Seismic image of the source zone and its implications for the rupture nucleation, Phys. Earth Planet. Inter., submitted. [6] R.J. O’Connell, B. Budiansky, Seismic velocities in dry and saturated cracked solids, J. Geophys. Res. 79 (1974) 5412^5426. [7] D. Zhao, A. Hasegawa, S. Horiuchi, Tomographic imaging of P and S wave velocity structure beneath northeastern Japan, J. Geophys. Res. 97 (1992) 19909^19928. [8] C. Kisslinger, Aftershocks and fault-zone properties, Adv. Geophys. 38 (1996) 1^36. [9] J.M. Lees, Tomographic P-wave velocity images of the Loma Prieta earthquake asperity, J. Geophys. Res. 17 (1992) 1433^1436. [10] P.A. Johnson, T.V. McEvilly, Park¢eld seismicity: Fluiddriven?, J. Geophys. Res. 100 (1995) 12937^12950. [11] D. Eberhart-Phillips, A.J. Michael, Three-dimensional velocity structure, seismicity, and fault structure in the Park¢eld region, central California, J. Geophys. Res. 98 (1993) 15737^15758.
EPSL 6696 2-7-03 Cyaan Magenta Geel Zwart
O.P. Mishra, D. Zhao / Earth and Planetary Science Letters 212 (2003) 393^405 [12] D. Zhao, T. Mizuno, Crack density and saturation rate in the 1995 Kobe earthquake region, Geophys. Res. Lett. 26 (1999) 3213^3216. [13] D. Zhao, H. Negishi, The 1995 Kobe earthquake: seismic image of the source zone and its implications for the rupture nucleation, J. Geophys. Res. 103 (1998) 9967^ 9986. [14] D. Zhao, O.P. Mishra, R. Sanda, In£uence of £uids and magma on earthquakes: seismological evidence, Phys. Earth Planet. Inter. 132 (2002) 249^267. [15] E. Roelo¡s, The Park¢eld, California earthquake experiments: An update in 2000, Curr. Sci. 79 (2000) 1226^ 1236. [16] S.A. Miller, Fluid-mediated in£uence of adjacent thrusting on the seismic cycle at Park¢eld, Nature 382 (1996) 799^802. [17] B.R. Arora, G. Rawat, A.K. Singh, Mid-crustal conductor below the Kutch Rift Basin and its seismogenic relevance to the 2001 Bhuj earthquake, DCS-DST News, Department of Science and Technology, Govt. of India, New Delhi, August 2002. [18] H.K. Gupta, S.V. Sarma, T. Harinarayana, G. Virupakshi, Fluids below the hypocentral region of Latur earthquake, India: Geophysical indicators, Geophys. Res. Lett. 23 (1996) 1569^1572. [19] S.K. Biswas, Regional tectonic framework, structure and evolution of the western marginal basins of India, Tectonophysics 135 (1987) 307^327. [20] A.C. Johnston, L.R. Canter, Earthquakes in stable continental crust, Sci. Am. 262 (1990) 68^75. [21] P. Talwani, A. Gangopadhyay, Tectonic framework of the Kachchh earthquake of 26 January 2001, Seism. Res. Lett. 72 (2001) 336^345. [22] J. Krishna, I.B. Singh, J.D. Howard, S.A. Jafar, Implications of new data on Mesozoic rocks of Kuchchh peninsula, western India, Nature 305 (1983) 790^792. [23] T. Utsu, Seismology, Kyoritsu Press, Tokyo, 1984, 310 pp. [24] J.M. Lees, H. Wu, Poisson’s ratio and porosity at Coso geothermal area, California, J. Volcanol. Geotherm. Res. 95 (2000) 157^173. [25] G.R. Pickett, Acoustic character logs and their applications in formation evaluation, J. Petrol. Tech. 15 (1963) 650^667. [26] Y. Takei, E¡ects of pore geometry on Vp/Vs: from equilibrium geometry to crack, J. Geophys. Res. 107 (2002) 10.1029/2001JB000522. [27] K. Yamamoto, M. Kosuga, T. Hirasawa, A theoretical method for determination of e¡ective elastic constants of isotropic composites, Tohoku Geophys. J. 28 (1981) 47^67. [28] K. Hadley, Comparison of calculated and observed crack densities and seismic velocities in westerly granite, J. Geophys. Res. 81 (1976) 3484^3493. [29] P. Talwani, Fault geometry and earthquakes in continental interiors, Tectonophysics 305 (1999) 371^379.
405
[30] S. Marshak, T. Paulsen, Structural style, regional distribution, and seismic implications of Midcontinental faultand-fold zones, United States, Seismol. Res. Lett. 68 (1997) 511^520. [31] M.N. Toksoz, C.H. Cheng, A. Timur, Velocities of seismic waves in porous rocks, Geophysics 41 (1976) 621^ 645. [32] D. Moos, M.D. Zoback, In situ studies of velocity in the fractured crystalline rocks, J. Geophys. Res. 88 (1983) 2345^2358. [33] C.H. Scholz, L.R. Sykes, Y.P. Aggarwal, Earthquake prediction: A physical basis, Science 181 (1973) 803^809. [34] M. Zimmer, M. Prasad, G. Mavko, Pressure and porosity in£uences on Vp^Vs ratio in unconsolidated sands, Leading Edge 21 (2002). [35] P. Talwani, Seismogenic properties of the crust inferred from recent studies of reservoir-induced seismicity ^ Application to Koyana, Curr. Sci. 79 (2000) 1327^1333. [36] S. Rojstaczer, Permeability changes associated with large earthquakes: An example from Loma Priets, California, Geology 20 (1992) 211^214. [37] S. Rojstaczer, S. Wolf, R. Michael, Permeability enhancement in the shallow crust as a cause of earthquake-induced hydrological changes, Nature 373 (1995) 237^ 239. [38] P. Malign, The seismology of extensional hydrothermal systems, Geotherm. Resour. Counc. Trans. 18 (1994) 17^22. [39] G. Vlahovic, C.A. Powell, Three-dimensional S wave velocity structure and Vp/Vs ratios in the New Madrid Seismic Zone, J. Geophys. Res. 106 (2001) 13501^ 13513. [40] I. Serrano, D. Zhao, J. Morale, 3-D crustal structure of the extensional Granada basin in the convergent boundary between the Eurasia and African plates, Tectonophysics 344 (2002) 61^79. [41] W.F. Brace, E. Silver, K. Hadley, C. Goyze, Cracks and pores: A closer look, Science 178 (1972) 162^164. [42] W. He, A. Hajash, D. Sparks, A model for porosity evolution during creep compaction of sand stones, Earth Planet. Sci. Lett. 197 (2002) 237^244. [43] M. Prasad, Correlating permeability with velocity using £ow zone indicators, SEG Expanded Abstracts, 69th Annual SEG Meeting 587, 1999. [44] B.L.N. Kennet, S. Widiyantoro, A low seismic wavespeed anomaly beneath northwestern India: a seismic signature of the Deccan plume, Earth Planet. Sci. Lett. 165 (1999) 145^155. [45] M. Prasad, M.H. Manghanani, E¡ects of pore and di¡erential pressure on compressional wave velocity and quality factor in Bera and Michigan sandstones, Geophysics 62 (1997) 1163^1176. [46] R. Shanker, P. Pande, Geoseismological studies of Kutch (Bhuj) earthquake of 26 January, 2001, J. Geol. Soc. India 58 (2001) 203^208.
EPSL 6696 2-7-03 Cyaan Magenta Geel Zwart