Seismogenesis of clustered seismicity beneath the Kangra–Chamba sector of northwest Himalaya: Constraints from 3D local earthquake tomography

Journal of Asian Earth Sciences 62 (2013) 638–646

Contents lists available at SciVerse ScienceDirect

Journal of Asian Earth Sciences journal homepage: www.elsevier.com/locate/jseaes

Seismogenesis of clustered seismicity beneath the Kangra–Chamba sector of northwest Himalaya: Constraints from 3D local earthquake tomography Naresh Kumar a,⇑, B.R. Arora b, Sagarika Mukhopadhyay c, D.K. Yadav a a

Wadia Institute of Himalayan Geology, Dehra Dun 248 001, India Uttarakhand State Council for Science & Technology, 6, Vasant Vihar Phase I, Dehra Dun 248 006, India c Department of Earth Sciences, Indian Institute of Technology, Roorkee 247 667, India b

a r t i c l e

i n f o

Article history: Received 24 March 2012 Received in revised form 3 October 2012 Accepted 4 November 2012 Available online 17 November 2012 Keywords: Velocity structure Earthquake tomography Himalayan seismicity Chamba Nappe Himalaya

a b s t r a c t To investigate subsurface structure and seismogenic layers, 3D velocity inversion was carried out in the source zone of 1905 Kangra earthquake (M8.0) in the northwestern Himalaya. P-wave and S-wave phase data of 159 earthquakes recorded by a network of 21 stations were used for this purpose. Inverted velocity tomograms up to a depth range of 18 km show significant variations of 14% in Vp and Vs and 6% in the Vp/Vs across the major tectonic zones in the region. Synthesis of seismicity pattern, velocity structure, distinctive focal mechanisms coupled with nature of stress distribution allows mapping of three different source regions that control regional seismotectonics. Accumulating strains are partly consumed by sliding of Chamba Nappe to the southwest through reverse-fault movements along Chamba/Panjal/Main Boundary Thrusts. This coupled with normal-fault type displacements along Chenab Normal Fault in the north account for low magnitude widespread seismicity in upper 8–10 km of the crust. At intermediate depths from 8 to 15 km, adjusting to residual compressive stresses, the detachment or lower end of the MBT slips to produce thrust dominated seismicity. Nucleation of secondary stresses in local NE–SW oriented structure interacts in complex manner with regional stresses to generate normal type earthquakes below the plane of detachment and therefore three seismic regimes at different depths produce intense seismicity in a block of 30  30 km2 centered NE to the epicenter of Kangra earthquake. Ó 2012 Elsevier Ltd. All rights reserved.

1. Introduction The occurrences of large numbers of earthquakes along the Himalayan arc are the manifestation of state of strains accumulated due to continued collision of Indian and Eurasian plates. Confinement of earthquakes to narrow latitudinal belt (Ni and Barazangi, 1984), space–time segmentation of earthquakes (Khattri and Tyagi, 1983; Kayal, 2001) and diversities in source mechanisms (Pandey et al., 1995; Bollinger et al., 2004; Arora, 2012) signify complexities of the tectonic setting controlling strain accumulation and release along the entire length of the Himalayan arc. The Kangra–Chamba region in NW Himalaya that formed the center of devastating 1905 Kangra earthquake (M8.0) has persistently been a zone of intense seismicity as revealed by the high frequency of moderate magnitude earthquakes as well as by the amount of energy released (Fig. 1). An earlier work by Kumar et al. (2009) described the design and first results of a campaign experiment operated to obtain better understanding of the seismicity and seismotectonics of the Kangra–Chamba re-

⇑ Corresponding author. Tel.: +91 135 2525474; fax: +91 135 2625212. E-mail address: [email protected] (N. Kumar). 1367-9120/$ - see front matter Ó 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.jseaes.2012.11.012

gion by suitably augmenting seismological network. Optimum 1D velocity model (Table 1) and improved locations of hypocenters were obtained by simultaneous inversion of travel times of 159 local earthquakes recorded by 21 densely distributed stations (Fig. 2). In the back ground of widely scattered seismic events, the distribution of well-determined epicenters show a close cluster confined to narrow pocket block of 30  30 km2, just northeast of the epicenter of great 1905 Kangra earthquake (Fig. 2). The optimum 1-D velocity structure permitted to image the Basement Thrust Fault (BTF) separating the down going India plate from the over-riding wedge of the Himalaya with a low velocity layer (Kumar et al., 2009). Further, stations corrections resulting as a by product of the VELEST inversion algorithm (Kissling, 1988) coupled with fault plane solutions of select earthquakes helped to establish linkage of diverse seismic character with tectonic elements of the region and indicated multiple source mechanism for the clustered seismicity in the Kangra– Chamba region. In the present work, we invert the seismic P-wave and S-wave phase data of Kumar et al. (2009) to deduce the 3-D velocity structure using the local earthquake tomography (LET) technique (Thurber, 1983) to elucidate further the seismotectonics of the Kangra–Chamba region in general and seismic source character of the clustered seismicity, in particular.

639

N. Kumar et al. / Journal of Asian Earth Sciences 62 (2013) 638–646

34

Seismic Events PT

M8 Kangra 1905

7.0 ≤ Μ < 8.0

32

SF

Chamba

Latitude °N

6.0 ≤ Μ < 7.0 5.0 ≤ Μ < 6.0

ST

33

Kin n

Kangra

aur

Ga

31 HF

T

Dhar

Dehradun

30

r hw al

MB T

MC

chul

a

T

29

New Delhi 28 74

75

76

77

78

79

80

81

82

Longitude °E Fig. 1. Map of NW Himalaya showing the cumulative seismic energy released from all magnitude earthquakes during 1551–2005. The spatial distribution is mapped by computing the cumulative energy for circular area of radius 20 km with 10 km overlap. The epicenters of big size earthquakes (M P 5.0), including 1905 great Kangra earthquake, are marked with circles. Well defined zones of intense seismicity are centerd around Kangra–Chamba, Kinnaur, Garhwal and Dharchula regions. Despite very few big earthquakes, the Kangra–Chamba region remains zone of intense seismicity even when examined for shorter data sub sets.

Table 1 Used 1D crustal model of the Kangra–Chamba region by Kumar et al. (2009). Depth (km)

Vp (km/s)

Vs (km/s)

0.0 10.0 15.0 18.0 46.0

5.27 5.55 5.45 6.24 8.25

3.01 3.21 3.05 3.59 4.73

2. Distinctive tectonic setting From geological and tectonic considerations the Kangra–Chamba region exhibits certain distinctive character (Fig. 2). In this segment the typical cross-section of the Himalaya divided into Siwalik Himalaya (SH), Lesser Himalaya (LH), High Himalayan Crystallines (HHC) and Tethyan Himalayan Sequence (THS) separated respectively by Main Boundary Thrust (MBT), Main Central Thrust (MCT) and Southern Tibet Detachment (STD). In the present study area, this generalized section is transgressed by the extended Chamba Nappe (CN) (Thakur, 1998). The weakly metamorphosed sediments of the CN are similar in facies to the THS. These are considered to have resulted from southwestward sliding of the THS from the north over the metamorphic HHC along the south dipping Chenab Normal Fault (CNF). Due to this southward transgression of the nappe, the LH sequence is pinched out to a narrow belt of few kilometers. This belt composed primarily of phyllites, slates and limestone is designated as the Panjal Imbricate Zone (PIZ) (Rautela and Thakur, 1992; Thakur, 1998). The MCT also loses its identity in this sector; instead the Punjal Thrust (PT) over-riding the PIZ from north is collectively designated as the MCT. The CNF separates the CN from the HHC. The southward extent of the CN is a subject of debate. The Chamba Thrust (CT) tracked as the northern border of the Chail formation was diagnosed to represent the southern limit of the CN at least west of Chamba (Fig. 2) whereas Rautela and Thakur (1992) indicated the CN to extend right up to the PT/

MBT. The CN to the NW and to the SE terminate against the HHC, separated by sharp NNE–SSW and NE–SW trending thrusts, across which the Kistwar (KW) and Rampur (RW) windows in the form of large antiform folds are exposed. 3. Method of analysis Well tested and widely used method of LET (Thurber, 1993; Eberhart-Phillips, 1993) is adopted to derive the 3-D velocity model for the Kangra–Chamba region. The method relies on the iterative simultaneous inversion of large numbers of travel times from well recorded 159 local earthquakes. The P-wave and S-wave phase data recorded by a network including 14 three-component digital and 7 vertical single component analog stations were used for this study (Fig. 2). S-wave arrival times are read only from digital records of horizontal components. Fig. 3 presents a map view and depth sections of event-station ray paths. The map view shows that the central block bounded between 31.90°N and 33.20°N, and 75.70°E and 77.00°E is well covered by ray paths, thus, defining the region considered for the velocity inversion. The data of a few events located on the fringe of this central block, where there are only a few ray paths, were supplemented by the phase data from some regional stations immediately bordering the study area (inset in Fig. 2). With this consideration, area of 150  150 km with center at 32.50°N and 76.50°E was finally considered for the inversion. Overall, a data set of 831 P-phases and 734 S-phases with selection criteria of minimum 6 phase readings, station azimuthal gap 6180° and rms error lower than 0.5 s for the individual earthquake were used for inversion. The inversion scheme could resolve velocity structure only for the upper 20 km of the crust as most of the earthquakes are confined in shallow section of upper crust. SIMULP12 inversion algorithm (Thurber, 1993; Eberhart-Phillips, 1993) used yielded estimates of hypocentral parameters together with Vp and Vp/Vs structure. The Vs is obtained at the last step of iteration from the estimated values of Vp and Vp/Vs. The optimum 1D velocity model obtained by Kumar et al. (2009) was proved effective as an initial input model for the 3D inversion attempted here. We first verified the robustness of the various model

640

N. Kumar et al. / Journal of Asian Earth Sciences 62 (2013) 638–646

Fig. 2. Locations of seismic stations and epicenters of local earthquakes used in the present study superimposed on the geology and tectonic map of the Kangra–Chamba region of HW Himalaya (after Thakur, 1998). The star is the epicenter of M8.0 Kangra earthquake of 1905. The configuration of surface grid nodes (⁄) used for 3-D velocity inversions by LET is also marked. The central 100  100 km2 block for which velocity structure is well resolved is shown by dotted box. The lines AB and CD are the profiles referred to in Fig. 7. SH, Siwalik Himalaya; LH, Lesser Himalaya; HHC, High Himalayan Crystallines; KW. Kistwar Window; RW, Rampur Window; STD. South Tibetan Detachment; CNF, Chenab Normal Fault; CT, Chamba Thrust; MCT, Main Central Thrust; MBT, Mian Boundary Thrust; PT, Panjal Thrust. Inset shows the location of the present study area in the framework of NW India and regional seismic stations.

parameters that control the resolution and stability of the inverted 3-D velocity structure from the given data set. Some of the salient features are summarized here, more details can be found in Kumar (2010). Guided by the inter-station spacing, adopted optimum 1-D velocity and initial trial testing, grid spacing (Fig. 2) of 25 km in horizontal direction and 2–5 km in vertical direction were adopted for the tomographic inversions. For the horizontal grid spacing of 25 km, the numbers of rays passing through or near a grid node (hit count or khit) varied between 200 and 1000, independently for P-waves and S-waves. Analogous to ray hit count, the derivative weight sum (DWS) quantifies the relative ray density in the volume of influence around the model node. Weighting the each ray segment as a function of its distance to the model node makes DWS more sensitive measure of the spatial sampling of the model space (Toomey and Foulger, 1989). The DWS for P-wave (S–P) was uniformly greater than 500 (450) for the central part considered for the inversion (Fig. 4a). The spatial distribution of DWS for Vs is mere replica of the S–P and hence not depicted here. Following Toomey and Foulger (1989), as an additional quantitative measure of the estimate the quality of the inversion, the resolution matrix (R) associated with inverted Vp and Vp/Vs were calculated separately. The spatial distribution contour plots of the diagonal element of the resolution matrix associated Vp and Vp/Vs are shown respectively in Fig. 4b and c for different depth levels. The resolution matrix coupled with high density of DWS shows that al-

most the entire central block is well sampled by ray paths and, thus, ensured robust velocity estimation. The solution stability also depends upon the damping parameters that are determined by data distribution and grid spacing (Eberhat-Philips, 1986). Following the well laid procedure in Eberhat-Philips (1986), we empirically obtained the damping parameters corresponding to point of inflexion on the curve defined by plots between data variance (misfit) with respect to model variance, separately for Vp and Vp/Vs respectively (Fig. 5a and b). We adopt the damping values of 40 and 50 for Vp and Vp/Vs respectively. The stability of the inverted velocity and Vp/Vs structure is checked using an efficient and widely used method of Checkerbroad Resolution Test (CRT) by introducing velocity perturbations, both positive and negative, in the modeling space (Zhao et al., 1992). Following the steps outlined in Singh et al. (2012), synthetic arrival time data for the entire 3D working space were generated by adding to the original arrival time, a time delay corresponding to velocity perturbations up to ±10% respectively at alternative grid nodes (Fig. 4d). The synthetic checkerboard arrival time data, both for Vp and Vs, were further superposed with random Gaussian noise that varied from ±0.1 to ±0.2 and were then inverted fresh using the model parameters adopted in the inversion of the original (real) data. The sections of working space where alternate positive and negative pattern are recovered adequately by inversion represent area where velocity structures are well re-

N. Kumar et al. / Journal of Asian Earth Sciences 62 (2013) 638–646

641

physical interpretation in terms of the seismotectonics of this block. 4. Results and discussion 4.1. 3D Velocity structure

Fig. 3. Spatial ray paths, ray paths in depth along longitude and latitude. The earthquake locations (green circle filled with red color), stations (blue triangles and ray paths (black solid lines) in three-dimensions. MSL: Mean Sea Level. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

solved. Fig. 4d shows the results of this CRT for Vp. It is found checkerboard pattern is generally well recovered for the northwest oriented corridor. Determined by tectonic fabric, large part of the observed seismicity is concentrated along this corridor (Figs. 2 and 3a). Near uniform distribution of hypocenters in this corridor rendered better criss-crossings of seismic rays at all depths, as evidenced by high Hit count and DWS, and hence checkerboard patterns are well recovered at all depths especially along this corridor. The high resolution both for Vp, and Vp/Vs is consistent with the high values of the diagonal elements of the resolution matrix at all depths (Fig. 4b and c). However, resolution appears to be reduced at depths greater than 15 km when hit counts and DWS are low due to lesser number of earthquake events at such depths. It is note worthy that even at this depth, the high value of resolution matrix accompanied by high value of Vp and Vp/Vs overlapping with section of clustered seismicity is robust and merit

The 3D variations in Vp, Vs and the ratio Vp/Vs for the well resolved central area by LET are shown in Fig. 6 as horizontal slices corresponding to depths of 0, 2, 5, 10, 15 and 18 km. To exhibit the true nature of even weak lateral velocity variation at any given depth, all velocity parameters are plotted as percentage perturbation relative to the initial 1D layer velocity. Within the resolution limits, these images clearly indicate the presence of strong lateral heterogeneities within the investigated area up to the depth range of upper 20 km. The similar strong heterogeneous velocity character in other parts of the Himalaya has also been inferred beneath the southeastern part of the Garhwal Himalaya by Mukhopadhyay and Kayal (2003) and to the northwestern part in the Hazara arc by Ni et al. (1991). The Vp and Vs tomograms show low velocity zone (LVZ) in the SW part aligned with the MBT and PT. The LVZ consistently imaged from surface to depths of 10 km can be related with the low velocity sedimentary sequence of Siwaliks and its near sub-vertical sub-surface extension. Further given that the LVZ is better developed and outlined in Vs tomograms, part velocity variability may be related to varying degree of pore fluids saturation (Mishra and Zhao, 2003). Away from these LVZs, much of the northern and northernwestern part of the study area is characterized by relatively high velocity zone (HVZ) where both Vp and Vs are higher by 5–6% from surface down to 10 km. In optimum 1D velocity model, the top 10 km thick layer with P-wave velocity of 5.23 km/s, much lower than the typical velocity of upper continental crust, was related to the thick sheet of weakly metamorphosed sediments of CN that extensively cover the study area. Given this correspondence, imaged HVZ constrains the space-depth extent of CN metamorphosed sedimentary sequence. From the surface to the velocity tomograms of 5 km, spatial expanse of the high velocity CN sequence to the south and southeastern part extent is limited by the LVZ. The vertical extent of the low velocity Siwalik, and relatively high velocity CN, do not seem to extend beyond the depth of 10 km as both the LVZ and HVZ show somewhat scattered picture in the velocity tomogram of 10 km. The change in lithocomposition of crust in the depth range of 10 km is supported by the transition in Vp/Vs ratio from a uniformly high value in the upper 10 km to diffused low values below the depth of 10 km. Another distinctive feature of velocity tomograms is a narrow NE–SW aligned LVZ in the eastern part, on the first appearance seen as local extrusion of the primary NW–SE trending LVZ. However, this local LVZ is confined to depth interval of 2–5 km and is underlined by a high velocity zone in the tomograms of 10 km depth and beyond (Fig. 6a and b). The presence of a structural discontinuity striking NE–SW, coincident with the LVZ, was also evident from the contour plots of station corrections obtained as a by-product of Joint Hypocentre Determination performed to obtain optimum ID velocity model by Kumar et al. (2009) using VELEST algorithm. It was seen that zone of positive station delays, coinciding with low-velocity Siwaliks in the SW quadrant, terminates sharply in NE–SW trend. 4.2. Seismicity-velocity structure linkage The efficacy of the LET derived 3D velocity model in providing the improved estimation of earthquake parameters is demonstrated by the fact that RMS of travel time residuals reduced by

642

N. Kumar et al. / Journal of Asian Earth Sciences 62 (2013) 638–646

Fig. 4. Resolution and stability tests for 3D inversion (a) derivative weight sum (DWS) for different depth sections, (b) diagonal elements of the resolution matrix (RDE) for Pwave, (c) RDE for P–S and (d) horizontal slices of synthetic checkerboard model test with 10% P-wave velocity relative to 1-D initial model.

34.38% from 0.32 s to 0.21 s for initial 1-D model to final 3-D model respectively and the accuracy of hypocenter determination is better than 3 km. To establish correlation of resolved earthquake parameters with crustal heterogeneity as revealed by velocity

structures, the distribution of hypocenters in specific depth intervals is superimposed on the velocity tomograms in Fig. 6. It is clear that most of the hypocenters at different depths are confined to the high velocity sections, representing the dry and compact blocks of

N. Kumar et al. / Journal of Asian Earth Sciences 62 (2013) 638–646

643

Fig. 5. Plots showing data variance versus model variance for selecting optimum damping value. The numbers along the curves are the damping values. Damping values of 40 and 50 corresponding to the point of inflexion on the curves are adopted for the inversion for P-wave and S–P waves respectively.

the crust. The low velocity sections are generally free from seismicity and there is cleared tendency for the concentration of hypocenters at the boundary between low–high velocity couplet (Fig. 6a and b). It is a common observation that hypocenters of many large earthquakes in varying tectonic settings are nucleated at the edge of the low velocity bodies, later ascribed to the highly fractured domain saturated with fluids (Zhao et al., 2010; Mishra and Zhao, 2003). To further access insight into the petrological significance of the above noted seismicity-velocity structure linkage, the depth distribution of hypocenters are superposed on the vertical section depicting variation in the true values Vp/Vs ratio along two profiles AB and CD, cutting across the major structural/tectonic grain of the study region (Fig. 7). The Vp/Vs being direct measure of the Poission’s ratio, the maps of Vp/Vs ratio provide tighter constraint in tracking the boundaries related with composition (lithologic) or with physical state of the medium, e.g. saturation state. As seen from Fig. 2 the profile AB cuts the section of the crust where the majority of seismic events are clustered in a narrow space whereas the other profile CD is located in the northwestern part of Chamba region where the seismicity is evenly distributed. Hypocenters of all earthquakes that fall within the 50 km wide corridor centered on the profiles are plotted to bring out most significant trends. On both cross-sections, numbers of well defined hypocenter alignments correlate well with resolved transition from high to low Vp/ Vs ratio. Given their connotation with the surface expressions of tectonic elements of the region, such coincidences of seismic alignments with transitions in Vp/Vs structure can be best seen as subsurface extensions of the active faults/thrusts of the study region. 4.3. Seismotectonic stress patterns To clarify this tectonic–seismic interaction and to decipher the underlying dynamics, we determined seismotectonic stress patterns from the inversion of fault plane solutions (FPSs) of selected earthquakes. To get regional pictures of prevailing stress in the study area, FPS of some 42 well located events (RMS value <0.1 s and M > 2.5) were computed using P-wave first polarity motions recorded at more than ten stations with maximum azimuthal station gap of less than 100°. The events analyzed included 10 earthquakes processed and interpreted by Kumar et al. (2009) and in addition FPS of some past big earthquakes were processed to examine dominant stress pattern associated with large magnitude earthquake than those recorded by the present campaign mode

experiment. Overall picture of stress in the form compression (Paxis) and extension (T-axis), obtained by tensor decomposition of the FPS, arrows is superimposed respectively on the Vp and Vs tomograms in Fig. 6 for specific depth range. The nature of stress and variations in the dip of earthquakes nodal plane on the two profiles AB and CD are projected in Fig. 7 to exhibit characteristics changes across the tectonic fabric of the Himalaya. The orientation of P-axis and T-axis indicates that stress pattern in the region is governed by two distinct tectonics (Fig. 6). Consistent with the dominance of collision tectonics in the Himalaya, compressive stresses prevail around the surface trace of the MBT, PT and CT (Figs. 6 and 7). The axes of maximum compression fall in the range of N20°W–N25°E. The FPS solution of these events indicated that dip of individual events varied from steep dip up to about 35° near surface to sub-horizontal in the depth of 15–18 km. The extensional stresses are dominant along the CNF on the northern fringe of the study region or in a small pocket east of the NE–SW directed structural discontinuity in the southeast corner of the study region. In the later case, the extensional stresses dominate only at depths greater than 15 km and direction of extensional stress as well as nodal plane deviate to align at right angles to the general strike of collision tectonics. In the tectonic cartoon (Fig. 8) depicting the principal elements of the seismotectonic model for Kangra– Chamba region, these above evidences are used to propose under-thrusting of a localized crustal block along the eastern boundary of the CN. It is suggested that interaction of the stresses generated by detach block of this localized under-thrusting with the down-going Indian plate could facilitate displacement on a NE–SW interface to generate earthquake by normal fault movement. In this and in the depth section along the profile AB, hypocenters of these set of earthquakes are marked by cluster ‘N’. 4.4. Elements of seismotectonic model Both on profiles AB and CD in Fig. 7, the lowest velocity structure imaged is characterized by high Vp/Vs and low Vs. This is in agreement with the minimum 1D velocity model, deduced by Kumar et al. (2009), which indicated presence of anomalous low velocity layer in the depth range of 15–18 km. Attributing the low velocity to the presence of fluids and noting that presence of the LVZ was a conspicuous feature of the Garhwal (Arora, 2003) and Nepal Himalaya (Schulte-Pelkun et al., 2005), the top interface of the mid-crustal LVZ was interpreted as brittle–ductile transition. As observed in many seismically active zones, the brittle–ductile

644

N. Kumar et al. / Journal of Asian Earth Sciences 62 (2013) 638–646

Fig. 6. The lateral variations in Vp, Vs and Vp/Vs at different depth sections of the central dotted box, marked in Fig. 3, where the velocity resolution is robust. Velocity parameters are plotted as percentage perturbation relative to the initial 1D layer velocity. The hypocenter distributions in specific depth intervals are superimposed on the velocity tomograms of respective depth. The compression (P-axis) and extensions (T-axis) obtained for intermediate size earthquakes in specific depth range are denoted with the converging and diverging arrows respectively. K and C represent the locations of Kangra and Chamba respectively. White star is the hypocentre location of M8.0 Kangra earthquake shown in depth section of 15 km.

transition defines the sharp lower cut-off depth of crustal seismicity (Sibson, 1986). Except for a narrow section right beneath the clustered zone of epicenters, the top interface of the high Vp/Vs in Fig. 7 traces this brittle–ductile transition as a northeast low angle dipping interface placed at average depth of 15 km near the MBT. The plunges of the axis of maximum stress for moderate to large earthquakes also tend to align with this interface. Such alignment of nodal planes on a single plane had led Ni and Barazangi (1984) to visualize such a linear plane as BTF, separating the down going Indian plate with the southward propagating sedimentary/

metamorphic wedge of the Himalaya. Such inference is supported when hypocentre location of M8.0 Kangra earthquake of 1905 taken from Molnar and Deng (1984) is examined in respect of the inverted tomographic images. It is observed that hypocentre of this earthquake is nucleated at the edge of low velocity zones in both Vp and Vs tomograms (Fig. 6). In the vertical-section, hypocenter at a depth of 15 km is located in the high Vp/Vs ratio (Fig. 7 – upper panel) low Vs (Fig. 6b – 15 km slice) zone. From these collations, it can be surmised that the boundary between high and low velocity bodies are the zones where maximum stress accumulation occurs.

645

N. Kumar et al. / Journal of Asian Earth Sciences 62 (2013) 638–646

Depth (km)

A

MBT PT PMT

5 0

-20 -30

Vp/Vs

P Q2

-10

B

CNF Q1

R R

N -20

0

20

40

60

80

100

Distance (km)

Depth (km)

C

5 0

MBT PT

D

CT

Vp/Vs

-10

S

-20 -30

-20

0

20

40

60

80

100

Distance (km) Fig. 7. Distribution of earthquake hypocenters in relation to velocity structure as indicated by Vp/Vs ratio on the depth sections of velocities along two profiles AB and CD across the major tectonics (location of profiles marked in Fig. 2). Well defined clusters, referred to in the text, are also marked. Double headed arrows indicate nature of stress, compression or extension. The slope of these lines with respect to horizontal denotes dip of the nodal plane provided by FPS. The white lines represent the detachment plane determined by Kumar et al. (2009) and white star represent the hypocentre location of great M8.0 Kangra earthquake. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 8. Tectonic cartoon summarizing principal seismotectonic components of the Kangra–Chamba region. In the collision regime of the Himalaya, clusters labeled as Q1, Q2 are due to thrust motion on the MBT, whereas clusters P and R are manifestation of strain consumption due to the sliding of the Chamba nappe to SW. Cluster N at depth below the plane of Basement Thrust Fault, later defining locus of large earthquake, is consequence related to the interaction of on two orthogonal structures. The confinement of clusters Q1, Q2, P and N in different depth section to a single vertical column explain the major clustering of epicenters, marked in Fig. 2.

For 1999 Chamoli earthquake (M6.3) similar observations were made by Mukhopadhyay and Kayal (2003). Consistent with Fig. 7, such sharp velocity interface can be viewed as seismically active BTF where large/great earthquakes are nucleated and denote extreme events of episodic release of inter-plate strains. In the depth sections above the plane of detachment, most of the hypocenter concentrations are confined to the sharp transition in velocity structures, quiet clearly reflected in Vp/Vs plots, demarcating major litho-tectonic boundaries. For example on the profile AB, the seismicity is concentrated to the LH, which is just few km wide. The three different concentrations are evident (Fig. 7 – upper panel and Fig. 8); two in the upper depth section of 7 km on the sub-parallel northeast dipping PT (cluster P) and MBT (cluster Q1); third on the deeper section of the MBT (cluster Q2). The later cluster Q2 is best developed where the northeast dipping MBT merges

with the detachment zone. In addition, a weak cluster (cluster R) can also be inferred along a plane, which is best seen as sub-surface extension of the southwest dipping CNF (Fig. 7 – upper panel). In comparison to this active nature of the PT and MBT on the profile AB, both these thrusts are seismically dormant along the profile CD. Instead, the hypocenter alignment is seen along the Chamba Thrust (CT) (cluster S). This marked difference in the behavior of seismicity on these two sectors appears to be related to the varying limits of southward extension of the CN. The southward extent of the CN on the western part of the study area is limited by the CT whereas on the eastern part, the CN extends right up to the PT/ MBT. This varying limit of southward extension of the CN is consistent with field evidence of Singh (1994) who noted that the CT bordering the Chail formation to the north is well marked west of Chamba, whereas this contact is not traced east of Chamba. In

646

N. Kumar et al. / Journal of Asian Earth Sciences 62 (2013) 638–646

the eastern part, the CN extends right up to the PT that is geologically supported by the evidence of ductile deformation preserved in the form of mylonitized augen gneiss exposed at the base of the PT (Rautela and Thakur, 1992). The plunge angles of the axes of maximum stress related to the thrust dominated earthquakes on the northeast dipping PT/MBT as well as on the CT, vary from large angles of 25°–35° in the shallow section to low dips of 10°–15° at the depth section of 10–15 km. This coupled with similar decreasing trend in the dip of nodal planes with depth, reported by Kumar et al. (2009), favors the tectonic model that steeply dipping MBT/PT/CT near surface flatten out at depth to merge with the detachment plane zone marking top of the down-going Indian plate. The earthquakes located in line with subsurface extension of the CNF are dominated by normal fault mechanism and expectedly show NE–SW directed extensional stresses. 5. Summary The synthesis of seismicity pattern, velocity structure, distinctive focal mechanisms coupled with the nature of stress distribution allow new insight into the operative dynamics that control the seismotectonics of the Kangra–Chamba region. Most coherent results are summarized in tectonic cartoon (Fig. 8) to define the principal components of seismotectonic model for this segment of the Himalaya. It is found that presence of low velocity layer in the depth range of 15–18 km is a general feature of the Himalayan collision zone. Fluids on the edge of the low–high velocity interface alter the rheological properties facilitating nucleation of strains and facilitating thrust dominated large and great earthquake along a single plane (BTF) between the down-going Indian plate and overriding wedge of the Himalaya. In the nappe dominated tectonics of the Kangra–Chamba region accumulating strains are partly consumed by normal-fault type displacements along the CNF and reverse-fault movements on the CT/PT/MBT, which account for the low magnitude wide spread seismicity in the upper 8–10 km of the crust. At intermediate depths of 8–15 km, adjusting to the residual compressive stresses, the detachment or the lower end of MBT slips episodically in stick slip manner to produce thrust dominated focused seismicity. The local NE–SW oriented structure in the southeast corner of the study area appears to play important role in producing local stress concentration. The nucleation of such secondary stresses at the base of the northeast dipping MBT may lead to extensional stresses producing normal type earthquakes seated below the plane of detachment. Since in small localized block, the three source regions of seismicity at different depths are vertically coincident (clusters Q1, Q2, P and N), the epicenter distribution map shows intense clustering of earthquakes of different mechanisms. Acknowledgements The authors acknowledge with thanks many fruitful discussions with V.C. Thakur, A.K. Dubey, J.R. Kayal, Sinha Roy and R.K. Chadha. The Directors of the host Institutions are thanked for providing the facilities and permission to publish the paper. The Ministry of Earth Science, Government of India is acknowledged for the project (MoES/P.O. (Seismo)/NPEP/15/2009). N.K. thanks ICTP, Trieste, Italy for providing Junior Associateship. Reviewers are thanked for improving the manuscript.

References Arora, B.R., 2003. Seismotectonics of the frontal Himalaya through the electrical conductivity imaging. In: Fujinawa, N., Yoshida, A. (Eds.), Seismotectonics in Convergent Plate Boundary, Terra Pub, Tokyo, pp. 261–272. Arora, B.R., Gahalaut, V., Kumar, N., 2012. Structural control on along-strike variation in the seismicity of the northwest Himalaya. J. Asian Earth Sci. 57, 15–24. Bollinger, L., Avouac, J.P., Cattin, R., Pandey, M.R., 2004. Stress buildup in the Himalaya. J. Geophys. Res. 109, B11405. http://dx.doi.org/10.1029/ 2003JB002911. Eberhat-Philips, D., 1986. Three dimensional velocity structure in northern California Coast Ranges from inversion of local earthquake arrival times. Bull. Seismol. Soc. Am. 76, 1025–12052. Eberhart-Phillips, D., 1993. Local earthquake tomography: earthquake source region. In: Iyer, H.M., Hirahara, K. (Eds.), Seismic Tomography: Theory and Practice. Chapman and Hall, London, pp. 613–643. Khattri, K.N., Tyagi, A.K., 1983. Seismicity patterns in the Himalayan plate boundary and identification of the areas of high seismic potential. Tectonophysics 96, 281–297. Kayal, J.R., 2001. Microearthquake activity in some parts of the Himalaya and the tectonic model. Tectonophysics 339, 331–351. Kissling, E., 1988. Geotomography with local earthquake data. Rev. Geophys. 26, 695–698. Kumar, N., Sharma, J., Arora, B.R., Mukhopadhyay, S., 2009. Seismotectonic model of the Kangra–Chamba sector of NW Himalaya: constraints from joint hypocenter determination and focal mechanism. Bull. Seismol. Soc. Am. 99 (1), doi: 10.1785.0120080220. Kumar, N., 2010. Quantification of seismic regimes of Kangra–Chamba region of NW Himalaya. Ph. D. Thesis. Indian Institute of Technology, Roorkee, India, p. 300. Mishra, O.P., Zhao, D., 2003. Crack density, saturation rate and porosity at the 2001 Bhuj, India, earthquake hypocenter: a fluid driven earthquake? Earth Planet Sci. Lett. 212, 393–405. Molnar, P., Deng, Q., 1984. Faulting associated with large earthquakes and the average rate of deformation in central and eastern Asia. J. Geophys. Res. 89 (7), 6203–6227. Mukhopadhyay, S., Kayal, J.R., 2003. Seismic Tomography structure of the 1999 Chamoli Earthquake source area in the Garhwal Himalaya. Bull. Seismol. Soc. Am 93 (4), 1854–1861. Ni, J., Barazangi, M., 1984. Seismotectonics of the Himalayan collision zone: geometry of the underthrusting Indian plate beneath the Himalaya. J. Geophys. Res. 89, 1147–1163. Ni, J., Ibenbrahim, A., Roecker, S., 1991. Three-dimensional velocity structure and hypocenters of earthquakes beneath the Hazara Arc, Pakistan: geometry of the underthrusting Indian plate. J. Geophys. Res. 96 (B12), 19865–19877. Pandey, M.R., Tandukar, R.P., Avouac, J.P., Lave, J., Massot, J.P., 1995. Interseismic strain accumulation in the Himalayan crustal ramp Inepal. Geophys. Res. Lett. 22, 741–754. Rautela, P., Thakur, V.C., 1992. Structural analysis of the Punjal Thrust Zone, Himachal Himalaya, India. J. Him. Geol. 3 (2), 195–207. Schulte-Pelkun, V., Monssalve, G., Sheehan, A., Pandey, M.R., Sapkota, S., Bilham, R., Wu, F., 2005. Imaging the Indian subcontinent beneath Himalaya. Nature 435, 1222–1225. http://dx.doi.org/10.1038/nature03678. Singh, A.K., Mishra, O.P., Yadav, R.B.S., Kumar, D., 2012. A new insight into crustal heterogeneity beneath the 2001 Bhuj earthquake region of northwest India and its implications for rupture initiations. J. Asian Earth Sci. 48, 31–42. Sibson, R.H., 1986. Brecciation process in fault zone: inferences from earthquake rupturing. Pure Appl. Geophy. 124, 159–175. Singh, K., 1994. Structural framework of the Chamba thrust sheet, Himachal Himalaya, India, Geoscience Journal XV, 85–92. Thakur, V.C., 1998. Structure of the Chamba Nappe and position of the Main Central Thrust in Kashmir Himalaya. J. Asian Earth Sci. 16, 269–282. Thurber, C.H., 1983. Earthquake location and three dimensional crustal velocity structure in the Coyote Lake area, central California. J. Geophys. Res. 88, 8226– 8236. Thurber, C.H., 1993. Local earthquake tomography velocities and Vp/Vs – theory. In: Iyer, H.M, Hirahara, K. (Eds.), Seismic tomography, Champan, Hall, London, 1993, pp. 563–580. Toomey, D.R., Foulger, G.R., 1989. Application of tomographic inversion to local earthquake data from the Hangill-Grensdalur central volcano complex, Iceland. J. Geophys. Res. 94, 17497–17510. Zhao, D., Hasegawa, A., Horiuchi, S., 1992. Tomographic imaging of P and S wave velocity structure beneath northeastern Japan. J. Geophys. Res. 97, 19909– 19928. Zhao, D., Santosh, M., Yamada, A., 2010. Dissecting large earthquakes in Japan: role of arc magma and fluids. Island Arc 19, 4–16.