Correlation between crustal anisotropy and seismogenic stress field beneath Shillong–Mikir Plateau and its vicinity in North East India

Abstract A systematic study towards understanding the correlation between polarization direction of crustal anisotropy with seismogenic stress field at different locations of the Shillong-Mikir Plateau and its vicinity in North East India is attempted. We used data from a 17-station broadband seismic network. In our earlier work , crustal anisotropic parameters were determined using ANISOMAT + for the 17 seismic stations. In this study, we have estimated stress field around the stations using focal mechanism solutions (FMS). Some 215 FMS are obtained by waveform inversion. These solutions are used for stress tensor inversion to estimate stress field around each location. It is observed that polarization direction of crustal anisotropy is consistent with that of the maximum horizontal stress (Ơmax) as well as the minimum horizontal stress (Ơmin). In addition to this, two orthogonal fast polarizations in some locations are also noted. The bivariate nature of correlations helps us to understand that the major mechanisms of seismic crustal anisotropy are not only due to the regional stress, but active faults and other geological conditions play a significant role in contemporary orientation of seismic crustal anisotropy and seismogenic stress field.


Introduction
The phenomenon of crustal anisotropy may be caused by the preferred alignment of micro-cracks or joints in the sedimentary bedding planes in the formation, or in the highly foliated metamorphic rocks (Margheriti et al. 1997). Sometime the cracks in the rocks are preferentially aligned in the maximum compressive stress direction and the anisotropy in the crust resulting from aligned cracks are used to estimate the crustal stress. Further, anisotropy can be observed near local active faults and lineaments with preferentially oriented cracks, which may indicate the presence of stress field (Pastori et al. 2012). In a region dominated by strike-slip mechanism, ơ 1 is horizontal and is represented by ơ H (maximum horizontal stress) and the minimum stress is represented as Ơ h . Generally, effective elastic properties of rocks and the crustal stress field are related to the alignment of microcracks and the seismic shear waves in the direction of ơ H are expected to travel faster than those in the direction of r h  Kayal 2008). Inset: Map of India showing the study region. (Nur and Simmons 1969;Nur 1971;O'Connell and Budiansky 1974;Hudson 1981). The shear wave splitting is commonly characterized by two parameters: Polarization direction (/) and time delay (dt) between the fast and slow shear waves. Shear wave splitting helps to understand characteristics of anisotropy in the crust below seismic stations. A detailed description of the method, analysis and results is given in our earlier work (Sharma et al. 2017).
The Shillong -Mikir plateau and its surrounding area in North East India (hereafter called, NE India) is considered as one of the most seismically active zones in the world (Kayal 2008). The region comprises of a number of faults with varying kinematics, lineaments, shear zones and intricate composite of rock masses of diverse lithology. The seismicity in this region is caused by the collision tectonics at the Himalayan arc to the north and by subduction tectonics at the Indo-Burmese arc to the east. Several studies during last 50 years in this region reported seismogenesis and the stress field (e.g. Verma and Kumar 1987;Kayal and De 1991;Nandy 2001;Mitra et al. 2004;Kayal et al. 2006;Bhattacharya et al. 2008;Angelier and Baruah 2009;Kayal et al. 2012;Baruah et al. 2013;Baruah et al. 2016).
In a recent work on crustal seismic anisotropy, Sharma et al. (2017) reported local crustal anisotropic characteristics in Shillong-Mikir plateau but shedding no light on the stress tensor behavior. In the present study, we try to investigate the localized seismic stress pattern around 17 seismic stations by Stress Tensor Inversion (STI) technique. Focal mechanism solutions (FMS) of some 118 well located earthquakes are obtained by waveform inversion and 97 FMS are taken from published results. Stress tensors of 215 FMS are studied, and a comprehensive interpretation with the shear wave velocity anisotropy is made here to shed light on its tectonic implications.

Tectonic settings
The NE region of India falls in zone V in seismic zonation map of India (BIS 2006) and evaluated as the highest seismic risk/vulnerable zone of the country. The intraplate zone of the study area is seismically much active due to Himalayan collision tectonics to the north and Indo-Burma subduction tectonics to the east and presence of several active faults in the intra-plate zone of the region (Figure 1). The tectonics of the study region can be characterized by the splays of north dipping major thrust faults formed during the Indian and Eurasian plate's convergence (Pandey et al. 2017). The region has produced two great earthquakes, the 1897 Shillong earthquake (Ms 8.7, revised Mw 8.1) in the intraplate zone and the 1950 Assam earthquake (Ms 8.7,revised Mw 8.4) in the eastern Himalaya syntaxis zone and several large earthquakes (Mw> 7.0) in the surrounding areas (Figure 1). In a recent work, Bilham and England (2001) proposed that the 1897 Shillong plateau earthquake occurred by popup tectonics at the south dipping Oldham fault by thrust faulting (Figure 1). The Shillong plateau, a part of the Indian shield, separated and moved to east along the Dauki fault (Evans 1964). The E-W trending Dauki fault separates the continental crust of the Indian shield (Shillong plateau) from the Cretaceous ocean floor hosting Bengal basin (Chen and Molnar 1990). The Dapsi thrust, northwest extension of the Dauki fault separates the granite gneiss and the Tertiary sediments within the Shillong plateau (Kayal et al. 2012). The NW-SE trending long Kopili fault in the Assam valley, on the other hand, separates the Shillong plateau and Mikir massif, a fragmented part of the plateau. The Kopili fault, previously named Kopili lineament in the literature (Nandy 2001), is identified to be intensively seismogenic where earthquakes occur from a shallow depth down to 40-50 km (Bhattacharya et al. 2008;Kayal et al. 2012). The long transverse structure of the Kopili fault may be existing since a time predating the birth of Himalayas, which accommodates the oblique plate convergence by conjugate shear failure, and the earthquakes occur by strike slipmechanisms (e.g. Mukhopadhyay 1984;Dasguptaet al. 1987;Kayal 2001;Kayal et al. 2012). The 21 September 2009 earthquake (Mw 6.3) and its aftershocks in the Bhutan Himalaya occurred on the N-S segment of the curvilinear MCT. Kayal et al. (2010) Kayal et al. (2012) and Bhattacharya et al. (2010) argued that the intensively active Kopili fault, which had produced two large earthquakes in the past, the 1869 Cachar earthquake Mw 7.4 and the 1943 Assam earthquake Mw 7.1, respectively, is vulnerable for an impending large earthquake in the near future.

Data and analysis
A 17-station broadband seismic network in the study area was under operation by the CSIR-NEIST, Jorhat; CSIR-NGRI, Hyderabad and IMD, New Delhi during 2001-2014. The stations were equipped with Nanometrics Trillium 120P and CMG 3T velocity meters, sampled at 100 Hz; the records were available in real time at the CSIR-NEIST in Jorhat. Some 215 well located earthquakes, (M D 1.7-4.8; M D : Duration magnitude) within 100 km to the nearest station, are selected for the present study. We used the HYPOCENTER programme developed by Lienert et al. (1986) and the velocity model of Bhattacharya et al. (2008) for locating the earthquakes. The estimated epicentre and depth errors are within ±2.1 km and ±5.1 km, respectively. The estimated root mean square (rms) errors, difference between the computed and observed travel times, are smaller than 0.4 s.

Anisotropy analysis
We have discussed the detailed analysis of crustal shear wave splitting in our earlier work (Sharma et al. 2017). Out of 215 well located events, 163 events fulfilled the criteria for the analysis. The results included 77 null measurements where delay time is equal to zero. Null measurements are observed when the original polarization of the shear wave is parallel to the fast or slow directions of anisotropic media (Schutt et al. 1998;Pastori et al. 2009). We used ANISOMATþ (Piccinini et al. 2013), a Matlab-based software developed by Bowman and Ando (1987). An appropriate window length of the horizontal seismograms containing the S-wave arrival was chosen and thereafter, utilizing the cross-correlation (CC) method, the anisotropic parameters were derived that best fits the model in which the S-wave splits into S-fast and Sslow components. In this model, the two components of S-wave have orthogonal polarizations and different velocities but exhibit the same shape (Piccinini et al. 2013). With this assumption, the N-S and E-W components of the seismograms are rotated in the horizontal plane by 1 o step through 180 o . The CC of two components was performed at each step, and the correlation coefficient with its respective time lag was computed. The software used only a small part of the waveform, including the S-wave arrival so that it does not contain converted phase (Piccinini et al. 2013). To improve the signal-to-noise ratio, all waveforms were band-pass filtered between 1and 20 Hz using a band-pass Butterworth fourth-order two-pass filter. We applied this filter based on the dominant frequency for local earthquake shear wave, which ranges from 2 to 15 Hz. A time window length of about 0.3 s was selected, fixing the start of the windows 0.1 s before the S arrival pick. To ensure consistency of the results we considered only those measurements that showed a CC coefficient larger than or equal to 0.7 (Sharma et al. 2017). Figure 2 illustrates the average fast polarization of shear-wave and time delay of slow shear waves.

Focal mechanism solutions
Out of 215 FMS (Figure 3), 118 FMS are determined in this study by waveform inversion technique. Minimum three-station data are used for the FMS. In waveform inversion, single station or two station data may produce fairly precise solution (e.g. Delouis and Legrand 1999;Fojt ıkova and Zahradnik 2014;Zahradn ık et al. 2015;. Various studies suggested that stress field is independent of magnitude of earthquakes (Angelier 1984;Michael 1984Michael , 1987Angelier and Baruah 2009). Fault plane solutions of lower magnitude events are successfully used for stress tensor inversion by many authors (e.g. Fan and Wallace 1991;Dreger and Helmberger 1993;Mancilla et al. 2002;Dreger 2003;Musumeci et al. 2005;Zhu et al. 2006;D'Amico et al. 2010D'Amico et al. , 2011. The FMS of the events are determined using the ISOLA code of Zahradnik et al. (2001) and Zahradn ık and Plesinger (2005). The ISOLA code is based on multiple point-source representation and iterative deconvolution method (Kikuchi and Kanamori 1991). In this method the complete velocity records are used without selecting any particular phase. The 3-component velocity records are first converted into displacement waveforms and then the displacement waveforms are low pass filtered below the corner frequency to remove any offset. The components with high frequency are excluded because it is difficult to model and it requires a precise knowledge of detailed subsurface crustal velocity structure. Necessary DC (double-couple component) removal and trend line removal are also performed. Green's functions are then computed in the complex spectral domain using suitable crustal velocity model pertinent to the study area at a point source by Discrete Wave (DW) number method (Bouchon1981) using the AXITRA program (Coutant 1989). The Green's functions are then convolved with appropriate moment tensor solutions based on a point source model and simple triangular source time function. The inversion is carried out for a frequency band (0.01-0.3 Hz) which is free of noise or has higher signal-to-noise ratio and occur below the corner frequency. A fine grid search of the strike, dip and rake is performed for the best depth and for the best moment solution. The final best fitting solutions are accomplished by comparing the observed seismograms and the synthetic seismograms (Baruah et al. 2016). A suitable velocity model given by Bhattacharya et al. (2008) is used here for particular stations as an input to the waveform inversion program. The figures (I, II, III and IV) in the supplementary material represent two examples showing efficacy of the waveform inversion technique used in this study.

Stress tensor inversion
Stress tensor inversion is carried out using results of all the 215 fault plane solutions. We used Michael (1984Michael ( , 1987 method for the stress tensor inversion analysis. The Michael method's algorithm uses the statistical method of bootstrap resampling and allows determining orientation of three principal stresses (ơ 1 ¼ maximum, ơ 2 ¼ intermediate and ơ 3 ¼ minimum). These parameters are determined by finding the best fit stress tensor to the observed solutions. The assumptions made for the input data are: (i) stress is uniform in the area during the observed time interval, (ii) the earthquakes are shear-dislocations on pre-existing fault(s), and (iii) slip occurs in the direction of the resolved shear stress on the fault plane (Baruah et al. 2016). The results of stress tensor inversion are shown in Figures 4 and 5 and other details are given in Table 1.

Crustal anisotropy
The fast polarization directions obtained by shears wave splitting analysis at the 17 stations, selecting some 163 local earthquakes, exhibit crustal anisotropy beneath the study area (Figure 2). It is observed that fast polarization directions vary from station to station. Local active faults play a significant role in addition to regional stress. Out of 17 stations, 11 stations represent fast shear wave polarization in near NE-SW direction. Among the remaining six stations, the BKD, SHL and TURA stations show two orthogonal fast directions; these stations are situated near two orthogonal active faults (Figure 2). The station BOKO in Assam valley, near the E-W trending Oldham/Brahmaputra fault, shows E-W fast polarization. At the KZI and NLI stations, the fast polarization directions seem to be influenced by the local stress field as well as by the local active faults or geologic conditions.

Focal mechanism solutions
FMS of the 215 earthquakes are illustrated in Figure 3. Out of 215 solutions, 118 FMS are determined in this study as mentioned above and rest 97 solutions are adopted from the published results (refer Table I Figure I of the supplementary material. The event was analyzed with a starting depth of 20 km. Seismic waveform data from four stations namely NGL, JPA, RUPA and TZR are used. In the inversion scheme, it is assumed that the source is located on a horizontal plane with a starting depth of 20 km (24.9 km being the estimated depth) and a spacing of 2 km for 9 trial source positions. The used frequency band for the inversion is fixed at 0.03-0.05 Hz with a cosine tapering to obtain detailed information about the source rupture process. The best waveform match between the observed and synthetic is obtained with estimated Green's function for the NGL, JPA, RUPA and TZR stations (refer Figure I in the supplementary material). The seismic moment at the four respective stations are estimated with an average value 4.39 Â 10 21 dyne-cm, which is basically used to scale the amplitudes. A preferred solution is obtained (red beachball) based on maximum correlation (refer Figure  .432 E, depth 21.8 km, duration magnitude M D 3.6). The event was analyzed taking a starting depth of 18 km. Seismic waveform data from three stations, namely RUPA, JPA and TZR, are used for the inversion. As mentioned above, the source position is assumed on a horizontal plane with a spacing of 2 km for 9 trial source positions. The used frequency band is also fixed at 0.03-0.05 Hz with cosine tapering as mentioned above. The best waveform match between observed and synthetic is obtained with estimated Green's function for the RUPA, JPA and TZR stations (refer Figure

Stress tensor inversion
Out of 17 seismic stations, the stress tensor inversion results are obtained for around 12 stations, using Michael's method, since the data were not sufficient at SIL, NLI, SNG, UDL and BKD stations. The results of the 12 stations are given in Table 1 and are illustrated in Figure 4 and Figure 5. Events mostly with normal faulting occurred around the AGIA station showing extensional stress ơ 2 and ơ 3 in NE-SW and NW-SE directions, respectively. Events with normal faulting to vertical dip-slip results were observed around the BOKO station, where ơ 2 and ơ 3 are directed along NE-SW and NW-SE, respectively. The events with strike-slip faulting solutions show compressional stress ơ 1 in NW-SE and intermediate extensional stress ơ 2 in NE-SW direction around the GAU station. The events with normal faulting solutions are observed around JPA station, where ơ 2 and ơ 3 are trending along NW-SE and NE-SW directions, respectively. The events with strike-slip solutions are observed around KZI station, where compressional stress ơ 1 is directed along NE-SW and extensional ơ 2 along NW-SE. For both the LNK and MND stations, strike-slip solutions are observed with ơ 1 along $ N-S direction and ơ 2 along $ E-W direction. Around the NGL station, strike-slip solutions are observed with compressional stress ơ 1 along NE-SW direction and extensional ơ 2 along NW-SE direction. Around RPA station, the events with strike -slip solutions show ơ 1 and ơ 2 along NE-SW and NW-SE, respectively. The events with strike-slip solutions around SHL station show ơ 1 and ơ 2 along $ NW-SE and $ NE-SW directions, respectively. Around SNG and TZR stations, strike-slip solutions are prevalent with ơ 1 along NE-SW direction and ơ 3 along NW-SE direction. The events with strike-slip solutions around TURA show ơ 1 along $ N-S direction and ơ 2 along $ E-W direction.
It is observed that the regional NE-SW compressional stress ơ 1 and NW-SE extensional stress ơ 3 are not consistent at several stations in the Shillong plateau. The local major tectonic features like the N-S trending Barapani Shear zone, E-W trending Dapsi Thrust and the Dauki Fault, N-S trending Dhubri fault and NW-SE trending Oldham Fault may have influenced the stress tensor inversion results around several adjacent stations. In Shillong plateau, the compressional stress around GAU and NGL stations are along NW-SE and NE-SW directions, respectively and around other stations, namely MND, TURA and SHL, it is within $ NW-SE to $ N-S directions. In Mikir plateau, on the other hand, the stations around KZI and TZR show compressional stress along NE-SW direction and extensional stress along NW-SE direction.

Discussion
In our earlier work (Sharma et al. 2017) we studied crustal anisotropy in the Shillong-Mikir plateau area, NE India using the 17-station broadband seismic network that was operated during 2001-14. Some 163 local earthquakes were selected for determining shear wave splitting by an automated cross-correlation (CC) method. The results are reproduced here in Figure 2. It is observed that fast polarization directions vary from station to station. Singh et al. (2006) reported NE-SW orientation of mean fast polarization in the mantle beneath the Shillong plateau using teleseismic SKS/SKKS phases. The observed variations in crustal anisotropy in this study may be explained by the preferred direction of fluid-filled cracks along the direction of maximum horizontal compressional stress in the region and local active faults with other geological conditions (Sharma et al. 2017). A small rotation of the average fast polarization from northeast trend to east-northeast is noted at the stations from SNG to TZR and SIL (Figure 2). This variation may be explained by clockwise rotation of northern part of the Kopili fault as reported by GPS measurements (Vernant et al. 2014). The E-W trending Brahmaputra fault may provide an additional impact on this rotation. (Sharma et al. 2017). A comparison of fast polarization directions and delay times with other geophysical results is also shown in Figure 2. It is noted that the NE-SW fast polarization at most of the stations is parallel to the plate motion direction obtained by the GPS (Sahu et al. 2006;Socquet et al. 2006;Gahalaut and Gahalaut 2007;Jade et al. 2007;Banerjee et al. 2008;Barman et al. 2017) and also parallel to the maximum horizontal compressive stress reported by stress tensor inversion (Baruah et al. 2013;Baruah et al. 2016). Further, the results of World Stress Map (WSM) project including borehole breakout results (Heidbach et al. 2016) and the GPE derived deviatoric stress (Baruah et al. 2016) are in fair agreement with the crustal seismic anisotropy pattern (Figure 2).
In continuation to the above study, we have analyzed some 215 fault plane solutions obtained by waveform inversion in the study area, and then studied the stress tensor to understand the crustal anisotropy variation in terms of the stress field. The stress tensor inversion results at around 15 stations are shown in Figures 4 and 5. The trend of compressive stress with minimum extensional stress was found around TZR. The regions around GAU, KZI, LNK, MND, NGL, RPA, SHL and TURA show maximum compressive stress with intermediate stress. Dominant extensional stress is observed around the AGIA, BOKO and JPA stations. The observation on crustal anisotropy at TZR, KZI, and NGL stations reveal average fast directions parallel to maximum compressional stress ơ max ( Figure 5). The maximum compressive stress direction around TZR and NGL stations is oriented at 25 0 -30 0 NE and the direction of crustal anisotropy is almost parallel to it. Around the KZI station, maximum compressive stress direction is towards the NW-SE direction; the crustal anisotropy direction is also parallel to it. Thus, is it is clear that around these five stations, the seismogenic stress field govern the crustal anisotropy direction. On the other hand, the fast polarization directions at the JPA and NLI orient parallel to the direction of minimum extensive stress ơ min . At the stations AGIA, GAU, RPA and SHL, the fast polarization directions are parallel with the intermediate extensional stress, i.e. ơ 2. But in BOKO and TURA stations, the fast polarization directions are nearly parallel with the intermediate extensional stress ơ 2. The anisotropy direction of these two stations seems to be influenced by the presence of NW-SE trending Oldham fault near the BOKO station and E-W trending Dapsi thrust near the TURA station. At LNK and MND stations, polarization direction differs by almost 30 0 -60 0 with the maximum stress, ơ max . Again, the anisotropy directions at these stations seem to be influenced by the presence of NW-SE trending Oldham fault near to the MND station and NW-SE trending Kopili fault zone near the LNK station. It may also be due to the presence of two possible anisotropic layers with two distinct rheology or may be due to the presence of one anisotropic layer with micro-cracks oriented in two different directions. It has long been assumed that tectonic stress at local or regional scale causes shear wave splitting (Savage et al. 2016). However, some studies have reported that in certain area, the anisotropy is better explained by structures such as fault fabric (e.g. Kaneshima 1990) and mineral orientation (e.g. Do Nascimento et al. 2004). In the near southern part of Brahmaputra, the fast polarization direction obtained through crustal anisotropy coincides with the minimum horizontal stress axis. Brahmaputra fault, with E-W trend, located in region may have some influence in the stress and anisotropic behavior.

Conclusions
An integrated study of shear wave splitting and stress tensor inversion of local earthquakes shows that the directions of fast polarization and stress tensor vary from station to station within the study area. At some stations, such as TZR and KZI in Mikir Plateau and NGL in Shillong Plateau, the fast polarization directions are parallel to the direction of maximum horizontal stress ơ max ; while the stress pattern is oriented in NE-SW direction around TZR, NGL stations and NW-SE around KZI station. On the other hand, the fast polarization directions at stations JPA in Shillong Plateau is parallel to the minimum horizontal stress, ơ min trending towards NE-SW around the JPA. The direction of fast polarization at stations AGIA, GAU, RPA and SHL are parallel with the intermediate extensional stress direction which trends along the NE-SW direction. The stations BOKO, MND and TURA located in Shillong plateau and LNK in Mikir plateau near the Kopili fault are not conformable with the direction of stress tensor, where local active structure/faults or geologic conditions may play a major role. The pop-up tectonics of the Shillong plateau as well as the transverse tectonics along the Kopili fault in the Assam valley has additional influence on the spatial variation of stress regime in the intraplate zone of our study area. Our observations suggest that crustal anisotropy in Shillong-Mikir plateau and in its vicinity depends not only on the regional stress field but also on the type of faulting, stress induced aligned cracks or on rheology of the rocks. The spatial variation of crustal anisotropy and stress tensor observations may play a key role for developing robust earthquake forecast models, however, with inculcation of other geophysical parameters (Dey et al. 2021). Future study may focus on temporal variation of the stress tensor in the anomalous zones where stress tensors are not conformable with the fast polarization.