Study of fault plane solutions and stress drop using local broadband network data: the 2011 Sikkim Himalaya earthquake of Mw 6.9 and its aftershocks

1 Geosciences & Technology Division, CSIR-North East Institute of Science & Technology (NEIST), India 2. Academy of Scientific and Innovative Research (AcSIR), CSIR-NEIST Campus, India 3. Dept. of Geosciences, University of Malta, Malta 4. National Center for Seismology, Ministry of Earth Sciences, New-Delhi, India 5. Geohazards Research Center, King Abdulaziz University, Jeddah, Saudi Arabia 6. National Research Institute of Astronomy and Geophysics, Cairo, Egypt


Introduction
A powerful destructive earthquake of Mw 6.9 at a depth of ~47 km (GCMT report), struck on September 18, 2011 in north-west Sikkim, at the border between India and Nepal (Figure 1). The earthquake was widely felt in Sikkim, Assam, Meghalaya, northern parts of West Bengal, Bihar and parts of other eastern and northern regions of India. The earthquake caused severe damages to property and also caused 60 casualties. This is the strongest instrumentally recorded earthquake in the Sikkim Himalaya during the last two centuries [Pradhan et al. 2013]. The most significant signatures of the 2011 Sikkim earthquake are landslides and landslide-induced loss of life as well as damage to the economy and infrastructure due to which the maximum intensity is reached to VII (MM scale) ] Poor construction of the houses and the hilly landslide-prone terrain are believed to be the main cause of such extensive damage [e.g., INTACH 2012, Rai et al. 2012]. Al-together there were 292 aftershocks (Mw>1.0) reported by a local broadband close-spaced microearthquake network of the National Geophysical Research Institute (NGRI) [Ravi et al. 2012], but the national network of the India Meteorological Department (IMD) as well as the global network (USGS, NEIC) reported only four aftershocks Mw≥4.0 ( Figure 1). The IMD reported that the main shock was preceded by a foreshock of Mw 4.9 at 26 km depth, on June 3, 2011 ( Figure 1).
The devastating 2011 Mw 6.9 Sikkim earthquake highlighted the urgent need to reassess the seismotectonic of the region. Using past global seismological data, GPS data, and the fault plane solutions of the 2011 Sikkim earthquake, several authors have tried to review the seismotectonics of the Sikkim Himalaya [e.g. Rajendran et al. 2011, Dasgupta et al. 2013, Pradhan et al. 2013, Paul et al. 2015, Baruah et al. 2014and Mukul et al. 2014]. The transverse tectonics on the steeply dipping near vertical mantle reaching Tista fault zone was suggested to generate the main shock. Finite fault model analysis indicates that the surface level Peak Ground Acceleration (PGA) values near to source region might reach more than 0.3g [Raghukanth et al. 2012]. Ravi Kumar et al. [2012] made a detailed study of the aftershocks recorded by the NGRI local network, and concluded that the near vertical mantle reaching Tista fault caused the main shock and aftershocks. A rupture model using simulation of near-field and farfield accelerograms was presented by Joshi et al. [2012]. BARUAH ET AL.
2 Figure 1. Tectonic map of the Sikkim Himalaya [after De and Kayal 2003] showing the epicentres of EHB relocated earthquakes . The 18 September 2011 main shock and four significant aftershocks (M≥4.0) are shown by red stars. The foreshock of the earthquake is shown by filled black star. The open black stars show the two damaging strong earthquakes (M≥5.9) with CMT fault-plane solutions. The solution of the 1988 earthquake is from Banghar (1991)based on first motion plot, and the waveform inversion solutions of two smaller events of 2002 are from Baruah et al. [2012a,b]. MCT: Main Central Thrust, MBT: Main Boundary Thrust, HFT: Himalayan Frontal Thrust, and GL: Goalpara lineament. The blue line indicates the river system. The red beach ball indicates the GCMT solution while blue beach ball indicates the solutions obtained in this study. The fault plane solutions of the four aftershocks (named as 1, 2, 3, 4) are also shown as blue beach ball.
3 STUDY OF FAULT PLANE SOLUTIONS OF 2011 SIKKIM EARTHqUAKE Chopra et al [2013] determined source parameters of the main shock by strong ground motion modelling in terms of PGA.
In this study, we have determined fault plane solutions of the 2011 main shock and its four small aftershocks (Mw 3.9-4.5) by waveform inversion using the local broadband data of the IMD national network. Although the main shock GCMT was immediately reported, but fault plane solutions of the four small aftershocks were not reported earlier. We have critically examined these fault plane solutions of the main shock and the four selected aftershocks using the IMD permanent broadband network data, and we have also estimated the source parameters of these events by spectral analysis. These results, with their tectonic implications are presented here.

Tectonic setting
The formation of the Himalaya resulted due to the continent-continent collision of Indian plate with Eurasian plate at ~50 Ma [e.g. Allegre et al. 1984, Baccheschi andD'Amico 2014 andreference therein]. In the Sikkim Himalaya, the major tectonic features are the Main Boundary Thrust (MBT) and the Main Central Thrust (MCT) parallel to the Himalayan trend, and the NNW-SSE trending Tista and Gangtok lineaments, the WNW-ESE trending Goalpara lineament and the SW-NE trending Kanchanjangha lineament transverse to the Himalayan trend [GSI 2000]. The MCT in the Sikkim Himalaya has taken a curvilinear loop shape indicating thicker sediments (Figure 1). A conceptual tectonic model of the Himalaya was first suggested by Seeber et al. [1981] based on the known geology, geophysics and teleseismic hypocentral data. In this model, the two active thrusts, the MBT and MCT, are contemporaneous in nature. However, Ni and Barazangi [1984] has suggested that the MCT is dormant while the MBT is active. The interface between the subducting Indian slab and the Himalayan sedimentary wedge is named the "Plane of Detachment" in the proposed model [Seeber et al. 1981]. Further north, below the MCT lies the "ramp", referred to as "Basement Thrust Front" (also known as Basal Detachement/Decollement or Main Himalayan Thrust), that accumulates the tectonic stress due to northward movement of the Indian plate; abrupt release is believed to be the main cause of earthquakes on the plane of detachment at a shallower depth (0-20 km).
Based on local temporary network data, Kayal [2001 and2010] argued that the Himalayan shallow thrust faulting earthquakes in the western Himalaya fits fairly well with the most accepted conceptual tectonic model, but the deeper (30-60 km) strike-slip faulting earthquakes in the eastern Himalaya cannot be explained as plane of detachment earthquakes. The focal depth of the recent well-recorded 1988 strong earthquake (Mw 6.9) at the Bihar-Nepal border foothills region was given at 50-60 km by different agencies such as USGS, ISC and ERI. Furthermore, Monsalve et al. [2006] reported deeper (30-60 km) source zones of earthquakes in the foothills region; these results were based on permanent digital network data of Nepal, central Himalaya. After the 1988 strong foothills earthquake in the region, the eastern Himalaya produced the 2011 strong Sikkim earthquake (Mw 6.9) further north at the surface trace of the MCT (Figure 1). The 1988 event as well as the 2011 strong earthquakes in this region, have given us the opportunity to understand the earthquake processes better in the eastern Himalaya region. Unlike the 1988 foothills earthquake, the 2011 Sikkim Himalaya earthquake provided local broadband network data that may shed new light to our understanding of the eastern Himalaya earthquakes at the Himalaya Seismic Belt (HSB) zone.

Hypocenter Locations
In this study, the earthquake data recorded by the broadband seismic network of IMD in the North-east India region are used for analysis. Each station is equipped with tri-axial REFTEK 151B Broadband velocimeters with a sampling rate of 100 Hz. Several locations have been reported for the main shock of the 2011 Sikkim earthquake by different national and international agencies ( Table 1). The reported source depth was the most uncertain parameter, ranging from 10 km to 50 km. We have reanalysed the hypocentral parameters using nine local broadband station data of the IMD. We used the HYPOCENTER programme of Lineart et al. [1986] and a velocity model proposed by Acton et al. [2011] for the region. The estimated location errors are ±2.1 km and ±5.1 along the horizontal and depth respectively, and the root mean square error 1.3 sec. The same programme is also used for locations of the aftershocks.

Spectral Analysis
Source parameters are studied by spectral analysis. For spectral analysis, first we applied a base-line correction and a cosine taper to the vertical component of the P-wave seismograms. We then transformed the seismograms from velocity to displacement. We selected a time window of several seconds beginning just Table 1. Hypocentral and Focal parameters of 18 th September, 2011 Sikkim earthquake and its foreshock and aftershocks.  6.9 6.9 6.9 6.9 6.9 6.9 6.9 1.  before the first arrival, usually 5 sec before the onset of the first arrival, and calculated the amplitude spectrum using Fast Fourier Transformation (FFT). All spectra presented in this paper have been corrected for instrument response and for attenuation using the results of Thirunavukarasu et al. [2016] where q=100. We characterised the spectra at a constant level A 0 for the lower frequencies, and with a fall-off above the corner frequency (fc). Earthquake source parameters are estimated following the Brune [1970] model as modified by Hanks and Wyss [1972].
In agreement with most used theoretical models of seismic sources, the far-field displacement Ω(f) can be described by a -corner frequency model: ( 1) where, Ω 0 is the low-frequency spectral level, f C is the corner frequency, n is the high-frequency spectral falloff, and γ is a constant. If γ=1, Equation (1) is the spectral shape proposed by Brune [1970].
The average value of corner frequency 〈f c 〉 is expressed as follows: ( 2) where f c i is the corner frequency at the i th station and N is the number of stations used in this study. The antilog defines the inverse of logarithm.
Source parameters (seismic moment, moment magnitude, stress drop, and source radius) are then calculated for the main shock and aftershocks. The seismic moment (Mo) is estimated from the low-frequency level (Ω 0 ) using the relation of Keilis-Borok [1959], where, ρ =2700 kg/m 3 is the density, Vp is the average Pwave velocity, R is the source to receiver distance, and R θφ is the a co-efficient accounted for the combined correction for radiation pattern, free surface amplification and site effects. The root mean square average of radiation pattern co-efficient R θφ =0.52 is used for all the events [Boore and Boatwright 1984].
An average seismic moment 〈M 0 〉 from all of the stations is determined from the average of the logarithmic values using the following relation [Archuleta et al. 1982]: (4) where, N is the number of stations used.
The average source radius 〈 R(P)〉 is computed using the modified Brune [1970] formula for P-wave spectrum [Hanks and Wyss 1972]: (5) Similarly, stress drop is estimated using the following equation: We calculated moment magnitude from seismic moment following the definition of Hanks and Kanamori [1979]: (7) where M 0 is measured in dyne.cm.
For each seismic event, we determined A0 and fc and the corresponding source radius r and seismic moment M 0 for each seismic station with signal-to-noise ratio greater than 1. The final M 0 and R are the average of these values. Stress drop, average displacement and moment magnitude are then calculated from these averaged M 0 and R. The resulting source parameters for the investigated earthquakes are summarized in Table 2 3.3 Focal Mechanisms and Stress Drop Focal mechanisms are obtained by the moment tensor inversions. We have used the Fortran-Matlab based ISOA-GUI package developed by Sokos and Zahradnik [2008] for moment tensor inversion. In this process, the complete velocity records are used without selecting any particular phases. The 3-component digital records are first converted into displacement waveforms. The displacement waveforms are then filtered between 0.01 and 9 Hz using four pole ban-pass Butterworth filters. The high frequency components are excluded because it is difficult to model high frequency components as it requires a precise knowledge of detailed sub-surface crustal velocity model. Subsequently resampling of records from a frequency of 100 Hz to 33 Hz is carried having with transfer function of the seismometer. Necessary DC removal and Trend line removal are performed. Green's functions are then computed in the complex spectral domain using the suitable crustal velocity model of Acton et al. [2011] pertinent to the region at a point source by Discrete Wavenumber (DW) method [Bouchan 1981]. The Green's functions are then convolved with appropriate instrument response and The inversion scheme uses a frequency band of 0.02 to 0.11. A fine grid search of the strike, dip and rake is performed for the best depth and moment. The final validation of the best fit solutions is accomplished by comparing the observed and synthetic seismograms. The 1D velocity model that is used may not be perfect and a small time shift is required for maximum correlation between the observed and synthetics [Marzooqi et al. 2008]. In the inversion technique, the shift is not known and it becomes a non-linear model parameter which has to be inverted along with the moment tensor and optimal source depth. To account for the horizontal mislocation, the synthetics are shifted relative to the observed by changing the origin time a few seconds before/after the origin time during the inversion. In this study, the inversion is carried out for full moment tensor. The preferred solution is obtained by a simple grid search over the focal depth with certain steps and also over the origin time between 1-2 sec. before/after the location origin time, since the origin time trades off with focal depth. The solution that has a large percentage of variance and double couple component is selected.
There is no doubt that the usage of large station data with good azimuthal coverage yield a reliable focal mechanism solution. But in wave-form inversion technique, single station data or stations with less azimuthal coverage are sufficient to get an accurate solution. Several authors have undertaken comprehensive studies on the focal mechanism solutions derived from single or double station data e.g.: Rao [2009], Delouis and Legrand [1999], Fojtíkova and Zahradnik [2014], Zahradník et al. [2015] etc. In this study, we have used three components broad-band seismograms of nine stations e.g., TAD, DHU, TRA, GAU, ZIR, JOR, MOK, IMP, and DIB operated by IMD in the eight North-eastern states of India (Figure 2). The epicentral distances of the stations are 21.78 km, 219.38km, 281.58km, 348.52km, 529.42km, 576.74km, 612.92km, 619.17km, and 634.34km respectively. The distances are compatible with the ISOLA software for local data (epicentral distances ≤1000km). A best waveform match between observed and synthetic seismograms of 18th September, 2011 Sikkim earthquake are shown in Figure 3. Seismic moment 1.2e+026 dyne-cm (Mw 6.7) is obtained through averaging the seismic moments from all the components of the recording stations. Finally a preferred  solution is obtained (red beach ball) based on maximum correlation between source position and time-shift (Figure 4). It can be seen that the preferred solution (red beach ball) differs from the nearby solutions in terms of DC percentage and correlation value (Figures 4, 5). As illustrated in Figures 4 and 5, it is ascertained from the correlation value that the preferred solution is well correlated at a depth of 47 km for a time-shift of 0.6 sec having a DC% of 90. This result is in conformity with the results of GCMT data. Similarly, we have determined the focal mechanism solutions of the four aftershocks using the same technique. The final results are illustrated in Table 1.
In order to check the stability of our solutions we also applied the CAP waveform inversion technique Helmberger 1996, D'Amico et al. 2010]. Using this method, it is possible to separate the entire records into body waves and surface waves and model them separately avoiding influence by shallow crustal heterogeneities [D'Amico et al. 2010[D'Amico et al. , 2011. The CAP method minimizes the misfit between the observed and synthetic seismogram using a grid search to obtain the best moment magnitude, source depth and focal mechanism [Zhu and Helmberger 1996]. It has been proven to be a stable and reliable method for computing moment tensor solution for small and moderate events [D'Amico et al. 2010[D'Amico et al. , 2011[D'Amico et al. , 2012[D'Amico et al. , 2013[D'Amico et al. , 2014D'Amico 2014, Orecchio et al. 2016. In this case, uncertainties of focal mechanisms and depths were estimated by the methods and the statistical approaches described by Bevington and Robinson [2003]. In general, we found that the uncertainties associated with the focal mechanism solution parameters are in the order of a few degrees.

STUDY OF FAULT PLANE SOLUTIONS OF 2011 SIKKIM EARTHqUAKE
The spectral parameters are estimated by fitting every P-wave displacement spectrum with theoretical general model for displacement spectra as expressed by Equation (1). The waveform analysis tools in seismology SeismoGRAPHer by Abdelwahed [2012] is used to estimate the spectral parameters. In this technique the Marquard Linearized Least-square method is used to solve the non-linear problem [Press et al. 1989]. The best fits of the observed displacement spectra to the theoretical one at a number of stations are plotted in Figure 6. The av-BARUAH ET AL.
8 Figure 4. Plot of correlation vs time-shift, source depth and focal mechanism for the multiple source inversion for 18 th September, 2011 Darjeeling-Sikkim earthquake. The largest correlation was obtained for source depth 47 km and 0.6 sec time shift. The preferred solution is depicted by a red beach ball.  From our calculation it is seen that for stations with epicentral distance >200km, the value of Mw estimated from spectral analysis is lower than the Mw value estimated from waveform inversion. This discrepancy may be due to the attenuation factor considered in the calculation of spectral parameters. It is noted that  estimated the spectral parameters 2011 Sikkim earthquake by using different accelerograph records.
The reported values of corner frequency, seismic moment, stress drop and source radius were 0.12 Hz, 115 bars, 9.68 km, respectively, which are consistent with our results, except the stress-drop value. In addition to these, we have estimated the spectral parameters of all the four aftershocks which also show low value of stress-drop (Table 3).

Discussion and Conclusions
The most widely accepted tectonic model to understand the Himalayan large earthquakes has been the steady state tectonic model envisaged by Seeber et al. [1981] and Ni and Barazangi [1984]. This model suggested that in the Himalayan Seismic Belt (HSB), earthquakes between the MBT and MCT are confined above the plane of detachment. However, several authors [e.g. Mukhopadhyay and Dasgupta 1988, De and Kayal 2003, Baruah et al. 2014 have reported that in the plane of detachment, shallow thrust faulting earthquakes as envisaged in the HSB tectonic model, are not evident in the Sikkim Himalaya and its foot-hills region. The Moho in the Sikkim Himalaya has been interpreted by de la Torre et al. [2007] to lie at 40-80 km depth. The earthquakes at 20-40 km depth between Main Himalayan Thrust and the Moho contain several strikeslip earthquakes [de la Torre et al. 2007], which suggest the presence of strike-slip faults in the region that extend below MHT or the basal decollement (basal detachment) of the Sikkim Himalayan wedge. Mukul et al. [2014] suggested that both MHT and Lesser Himalayan Duplex (LHD) as active structures in Sikkim Himalaya. They emphasized that the Darjeeling-Sikkim area, may be releasing accumulated strain through frequent moderate and micro earthquakes and the seismicity in general appears to be located between MCT and MBT. Furthermore, they suggest shortening in the area accommodated along active east-west trending faults in LHD in Tista Half Window (THW). This view is supported by the presence of a thrust component in the focal mechanism solution of one of the aftershocks of 2011 Sikkim earthquake. The mainshock with dextral strike-slip movement may be related with the incessant Duplex activity of the weaker well foliated meta-sedimentary rock sequences of the region, combined with the active convergence of Indian and Eurasian plates.
In the present study, a plan and sectional view of a part of the Darjeeling-Sikkim-Tibet Himalayan wedge [modified after Mottram et al. 2011, Mukul et al. 2014 shows the indicative dextral strike-slip movement (  The GCMT solution of the 2011 Sikkim earthquake in the eastern Himalaya shows that it occurred at greater depth (~ 47 km) by strike-slip mechanism (Figure 1). A similar strike slip solution is obtained in this study using local broadband network data with an estimated depth of 47 km (Figure 1). We have also examined fault plane solutions of the four significant aftershocks (Mw 3.9-4.5) using the local network data. The waveform inversion results show that three aftershocks occurred by strike-slip mechanisms and one by thrust faulting with strike slip component. Additionally, fault plane solutions of the past two medium to strong earthquakes (1965 and 1980) in main shock epicenter area also show dominantly strike slip solutions (Figure 1). In all these solutions of the 2011 main shock, aftershocks and the past events (1965 and 1980), the NNW-SSE striking nodal plane is comparable with the Tista fault.
The agreement between the different solutions has also been numerically checked by applying the Kagan [1991] method [D'Amico et al. 2011 and reference therein]. The Kagan angle measures the rotation that should be applied to one earthquake source doublecouple to make it coincident with another one. The agreement is very good between being in the order of 85-90%. The small differences among the solutions can be attributed to a combination of several factors for example the use of different velocity model, or different frequency bands. However, we are confident that the final focal mechanisms obtained in this study are robustly determined.
The geologically mapped Tista lineament (GSI 2000) was reported as Tista fault by Kayal et al. [2011] and Baruah et al. [2014]. It may also be noted that the 1988 strong earthquake (Mw 6.9) in the foot-hills region also shows strike slip mechanism with an estimated depth 55-60 km (Figure 1) [Banghar 1991]. A NW-SE transverse fault, termed the West Patna fault, was inferred to be the causative fault for this strike slip earthquake (GSI 2000).
The source parameters that are estimated from the displacement spectra of the broadband seismograms show that the stress drop of the main shock and aftershock sequence has a low value, ranging from 14-38 bars. Similar observation is reported in geodetic measurements by Pradhan et al. [2013]. They have inferred an average coseismic displacement of ~11 mm at Phodong and ~9 mm at Taplejung station near the main shock epicentre. It is important to note that low stress drop in other parts of the world are interpreted as the result of re-rupturing after short healing times [Kanamori et al. 1993, Berberian et al. 1999]. Here, this argument, however, may not hold true because no large earthquake is reported in the Sikkim Himalaya during the recorded history except two medium magnitude earthquakes, one in 1965 (mb 5.8) and the other in 1980 (mb 6.0) respectively. The eastern Himalayan seismicity differs from that of the western Himalaya. The western Himalayan earthquakes are shallow and occur by thrust mechanisms, whereas the Sikkim Himalayan or eastern Himalayan earthquakes are deeper and occur by strike slip mechanisms. Several authors [e.g. Mukhopadhyay 1984, Engdhal et al. 1998, Kayal 2001, De and Kayal 2003, Monsalve et al. 2006and Baruah et al. 2014 reported that the earthquakes occur in the upper half of the crust but much activity is also observed in the lower crust and at the crust-mantle boundary beneath Sikkim and western Nepal. The aftershocks recorded by the local microearthquake network in the Sikkim Himalaya show a vertical structure down to 50 km depth [Ravi Kumar et al. 2012]. This structure is comparable with the Tista fault, which is a mantle reaching vertical fault , Baruah et al. 2014]. The 2011 main shock occurred on this vertical structure by strike slip mechanism in a lower friction condition than those shallow thrust faulting Himalayan earthquakes on the plane of detachment in the western Himalaya. In the eastern Himalaya, the long NE and NW trending transverse structures/faults from foredeep to the high Himalaya, or even beyond (Figure 1), may predate the birth of the Himalaya, producing intersecting patterns, and accommodating the plate convergence by conjugate shear failure on some of these faults [e.g. Mukhopadhyay 1984, Dasgupta et al. 1987]. Such shear failures may be explained by oblique transverse faulting to the direction of regional compression as it is happening in the eastern Himalaya. The mantle reaching Tista fault that existed before the two plates came into collision, reactivated in contrast to the thrust kind of slip of Himalayan seismic belt; this may cause lower stress drop. Presence of transverse structures in different parts of the Himalayan arc, including the Sikkim Himalaya, indicates that the deformation in the long arc is not at all continuous [e.g. Pandey et al. 1999]. Receiver function estimates to the Moho shows that the entire crust beneath the Sikkim Himalaya is seismogenic and the seismicity is diffused in nature [Acton et al. 2011]. The diffused seismicity may be due to the presence of more than one seismogenic fault in the Main Himalayan Thrust system, either out-of-sequence and or reactivated faults [Mukul et al. 2014]. While the GPS study indicates that the rate of convergence in the Sikkim Himalaya and the western part of Bhutan is about 12mm/year, out of which 5 mm/year of aseismic slip is being accommodating at the north of the Tista lineament [Mukul et al. 2014]. Moreover, the variable bvalue in the range of (0.74-1.18) in the Sikkim Himalaya region suggests heterogeneities in the lithosphere [Mishra et al. 2010], which in turn probably makes it favourable for low-friction conditions. Earthquakes due to deep-rooted transverse tectonics in the study region seem to indicate new trends to be followed for the plate tectonics study of India-Eurasia plate conversion.