Role of Intraplate Strike-Slip Earthquakes in Accommodating Convergence Across the Eastern Himalayan Plate Boundary System

The Eastern Himalayan plate boundary comprises one-third of the Himalayan orogen, produced by the ongoing Indo-Eurasia continent–continent collision [1-3]. To its east, the Indian plate obliquely subducts beneath the Burma microplate and results in the formation of the Indo-Burman fold-thrust belt (IBFB). These two plate boundaries meet at the Eastern Himalayan Syntaxis. The region immediately south and west of these plate boundaries is known as North-East India (NE India), comprising the Himalayan foreland basin (Brahmaputra Valley), the Shillong Plateau and Mikir Hills (uplifted blocks of Precambrian basement rocks), and the foreland basin of the IBFB (the Bengal Basin) [4, 5]. This unique geological setup of NE India is replete with intraplate earthquakes attesting to the distributed deformation within the continental interiors [6, 7] (Figure 1(a)).

The active intraplate structures in NE India are the near vertical dextral strike-slip Kopili Fault Zone (KFZ), which lies between the Shillong Plateau and the Mikir Hills; the south-dipping Oldham Fault, which is buried beneath the northern edge of the Shillong Plateau; and the north-dipping Dawki Fault, which bounds the Shillong Plateau to the south. It has been conjectured that the dextral strike-slip Dhubri-Chungthang Fault (DCF) Zone beneath the Sikkim Himalaya extends southward into the Himalayan foreland basin and bounds the Shillong Plateau along its western edge. Recent moderate intraplate crustal earthquakes have originated in the KFZ [7]. The Oldham Fault was modeled to have hosted the 1897 (M_w 8) Shillong earthquake [8]. The Dawki Fault forms a $∼$1 km high south-facing cliff along the southern margin of the Shillong Plateau and down-faults the Indian continental crust by $∼$18 km, juxtaposing the Bengal Basin sediments against the southern margin of the Plateau [9]. The 2011 Sikkim earthquake (M_w 6.9) occurred on a dextral strike-slip fault at a depth of $∼$53 km beneath northern Sikkim and was followed by moderate-to-small aftershocks which extended in the southeast direction [10]. The GANSSER seismological network revealed a linear belt of micro-seismicity extending further SE beneath the Himalayan foreland basin up to the western edge of the Shillong Plateau [11]. All these earthquakes were associated with the NW-SE-oriented DCF [10]. However, this fault does not have any geomorphological expression in NE India.

The Bengal Basin, south of the Shillong Plateau, is covered by 18–20 km thick sedimentary layers overlying a thinned (transitional) crust of $∼$20 km thickness [9]. The thickness of the sedimentary layers decreases toward the western and northwestern margin of the Bengal Basin with an increase in the crustal thickness to $∼$35 km beneath the Western Bengal Basin [12]. This region of transition is marked by a steep gradient in gravity and a gravity high on the western (continental) side [13]. This demarcates the passive continental margin of eastern India in Eocene before the deposition of the sediments derived from the Himalaya and is known as the Eocene Hinge Zone (EHZ) [4]. The EHZ trends NE-SW and the Eocene Sylhet limestone is the marker formation of the continental margin shelf-slope transition [14]. Thinned transitional crust SE of the EHZ is characterized by paleo-passive-margin rift faults [14, 15]. A recent study of two-dimensional seismic profiles and well-log data from the northern Bengal Basin have been used to constrain the timing of uplift of the Shillong Plateau by modeling the subsidence history of the basin [16]. Seismic profiles across the basin revealed the change from a passive margin having southward thickening sedimentary formations to a flexural basin with northward-thickening strata, as a consequence of uplift and loading by the Shillong Plateau. Moderate strike-slip earthquakes are observed within the Bengal Basin transitional crust between depths of 20 and 40 km [7, 9]. These have been conjectured to originate in these basin margin rift faults, beneath the thick pile of sediments, reactivated by the ongoing NNE-SSW compression across NE India [9].

Recent GPS geodetic measurements have suggested rotational strain in isolated rigid lithospheric blocks across NE India as a mechanism to accommodate the ongoing convergence. These include the clockwise rotation of the Shillong block and the Assam block (separated by the KFZ) with respect to the Indian plate, about a pole situated a few hundred kilometers west of the Shillong Plateau. The observed rotation rates for Shillong block and the Assam block have been modeled to be $∼$1.15° my⁻¹ and 1.13° my⁻¹, respectively [3]. This has led to $∼$3 mm y⁻¹ of dextral motion along the KFZ and eastward increasing convergence rate across the Dawki Fault from $∼$3 mm y⁻¹ to $∼$6 mm y⁻¹ (Figure 1(b)). Previous studies of earthquake source mechanism from NE India [6, 7, 17, 18] have revealed a large number of moderate-to-strong (5.0 < M_w < 6.8) intraplate oblique slip and strike-slip earthquakes. These earthquakes are not associated with either of the two plate boundaries and had mostly originated within the Indian lower crust. Recent deep events in the Brahmaputra Valley and beneath the Indo-Burman Ranges have similar strike-slip mechanism, obtained from the Global Centroid Moment Tensor (G-CMT) catalog [19, 20]. These intraplate earthquakes provide an opportunity to study the role of strike-slip faults in accommodating this distributed intraplate deformation. In this study, we model the source properties of recent moderate-to-strong (5.0 ≤ M_w ≤ 6.7) intraplate strike-slip earthquakes in NE India (Figure 1(a)) and use source directivity for the largest event (M_w 6.7) to identify the fault plane from the auxiliary plane and then through back-projection of high-frequency energy we establish the spatio-temporal history of the rupture. These mechanisms are combined with results from previous studies of earthquake source [6, 7, 10, 18], GPS geodetic velocity vectors [3], and strain-rate field to build a kinematic model of active tectonics in NE India.

We model the source mechanism of six recent earthquakes in NE India using broadband teleseismic waveform data recorded by the Global Digital Seismograph Network stations. The stations lie in the epicentral distance (∆) range of 30° to 80° (Figure 2(a)), and the data have been downloaded from the Incorporated Research Institutions for Seismology Data Management Center. The chosen distance range ensures that we avoid the effect of triplication from nearer distances and interference from the core-mantle refracted shear-wave (SKS) phases for further distances. The three-component (vertical, north-south, and east-west) waveforms have been deconvolved with their respective instrument responses and then reconvolved with a band-pass filter having corner frequencies of 0.066 Hz (15 seconds) and 0.01 Hz (100 seconds), to reproduce the long-period response of the World-Wide Standard Seismograph Network instruments. The horizontal components (north-south and east-west) are rotated to the great-circle-arc path to produce the radial and transverse components, respectively. We use the P-waveform on the vertical component and the SH-waveform on the transverse component for our analysis. The moment tensor inversion algorithm of McCaffrey [21, 22] has been used to invert the waveform data to obtain a double-couple fault plane solution (strike, dip, and rake of the nodal planes), hypocentral depth, and the seismic moment of the earthquake. An inversion window of 60 seconds for P- and SH-waveforms, containing the depth phases have been used for the inversion. The inversion requires starting source parameters, which are provided from the G-CMT. The velocity structure of the source-zone is a crustal layer with thickness and average V_s (one-dimensional shear wave velocity) taken from Mitra [9], underlain by a mantle half-space. Synthetic waveforms are computed for a point source, embedded in the velocity structure, by convolving a source-time function (STF) with the computed Green’s functions for direct arrivals (P or S), near-source reflections (pP, sP, or sS), and multiples. The STF is defined by a series of overlapping isosceles triangle with half-width duration of 1.5 seconds. The integral over time of this STF gives the seismic moment. The polarity of the first motion of direct P- and SH-phases constrains the orientation of the nodal planes. This requires a good azimuthal distribution of the stations (Figure 2(a)). The program minimizes the weighted squares of residuals between the observed and synthetic P- and SH-waveforms, using an iterative approach. As the waveform amplitude and shape are nonlinearly related to the source parameters, these parameters are iteratively perturbed toward minimizing the misfit. Relative weighting is done by providing the SH-waveforms half the weight of the P-waveforms due to its larger amplitude. The observed waveforms are azimuthally weighted by assigning a higher weight to seismograms from isolated azimuths compared to a cluster of stations in a given azimuth. Realignment of the inversion window is done by adjusting the position of the direct P-phase first arrivals in the synthetics with respect to the data. The best-fitting source parameters, comprising strike, dip, and rake of the nodal planes, hypocentral depth, and seismic moment, are obtained by minimizing the misfit of amplitude and shape between all observed and synthetic waveforms (Figure 2(b)). Detailed description of this technique for the computation of source mechanism using long-period teleseismic waveform can be found elsewhere in the literature [23, 24].

Uncertainties for the modeled parameters have been obtained by computing the misfit for fixed values of each parameter around its best-fitting solution, while keeping the other parameters free. Upon fixing the value of one parameter, within a range bracketing the best-fitting solution, the inversion is re-run for other source parameters and its effect on the focal mechanism and the fit to the waveform is noted. The uncertainty is considered within the interval where there is no change in the misfits of other source parameters. This shows the trade-off between the fixed parameter and other free parameters of the focal mechanism, while allowing an estimate of the uncertainty associated with the fixed parameter (Figure 2(c)). The misfit values are approximated as a Gaussian distribution and the ±1 standard deviation (±1σ) bound is estimated from the Gaussian fit. This uncertainty analysis has been performed for hypocentral depth, strike, dip, and rake of the nodal planes (Table 1).

Source parameters of six recent earthquakes, which occurred within the Indian continental crust in NE India, have been modeled using the above procedure. Three events are in the Brahmaputra Valley (events b, d, and e), one beneath the Himalayan foothills (event c), one in the northern edge of the KFZ (event a), and one beneath the Indo-Burman ranges (event f) (Figure 3(a) and Table 1). We present the results of our analysis of these earthquakes along with well-constrained source mechanisms of previous intraplate earthquakes, in each of these regions. The three earthquakes in the Brahmaputra Valley are of moderate-magnitude, with two in the western valley having predominantly strike-slip motion (events d and e) and the third in the central valley having oblique-thrust motion (event b). The strike-slip earthquakes originate at depths of 14 km and 35 km, respectively, and have steep nodal planes (dip > 70°) oriented in the NW-SE and NE-SW directions. The oblique-thrust fault earthquake (event b) lies north of the Shillong Plateau at a depth of 38 km, and its nodal planes strike nearly E-W. In the region of the two strike-slip earthquakes (events d and e), Chen and Molnar [18] modeled a thrust fault earthquake at a depth of 25 km (event k). These two strike-slip earthquakes (events d and e) are likely to be associated with the NW-SE-oriented (DCF which has been conjectured to extend into the Himalayan foreland in this region [11]. The earthquake beneath the Himalayan foothills (event c) originated at a depth of 8 km, possibly within the underthrusting Indian crust and has predominantly strike-slip motion. Further east, at the northern edge of the KFZ, the M_w 6.0 earthquake occurred at a depth of 36 km (event a). This is one of the largest events to have occurred in the Brahmaputra Valley in the recent past and has an oblique-thrust fault mechanism. The nodal plane oriented along the KFZ strikes NW (325°) and dips 55°. In the past, a number of similar moderate-sized oblique-thrust fault earthquakes had occurred along the northern margin of the Mikir Hills (events g, h, i, and j) and one beneath the Eastern Higher Himalaya (event l). All these events originated in the Indian lower crust at depths of 33–42 km. The NW-SE-oriented nodal plane for all these earthquakes, including the two beneath the Eastern Himalaya, is subparallel to the KFZ. The largest earthquake to have occurred in NE India in the twenty-first century is the M_w 6.7 Manipur earthquake (event f) beneath the Indo-Burman Ranges (IBR). This event originated at a depth of 52 km and has strike-slip mechanism. The nodal planes are oriented NNW-SSE and WSW-ENE with high-angle dip. A number of similar strike-slip earthquakes ranging in depth from 36 to 57 km have occurred beneath the IBR, originating within the subducted Indian plate (events v, w, x, y, and z) [7]. To understand the relationship of these earthquakes to the KFZ, it is important to isolate the fault plane from the nodal planes. The M_w 6.7 Manipur earthquake is sufficiently large to perform such analysis of rupture directivity and spatio-temporal evolution (Sections 2.1 and 2.2). East of these events, within the Bengal Basin, a number of strike-slip to oblique-thrust fault earthquakes had been modeled in previous studies (events m, n, o, p, q, r, s, t, and u). It had been conjectured that these events occurred on reactivated paleo-rift faults within the transitional crust [9].

Modeling the far-field displacement of P- and SH-waveforms constrain the geometry of the nodal planes and the orientation of the slip vector, but a fundamental ambiguity exists in distinguishing the fault plane from the auxiliary plane. Faulting in strong-to-great earthquakes displays directivity effects in rupture propagation, where the frequency content of the far-field displacement exhibits a strong azimuthal dependence, in the absence of strong frequency-dependent attenuation [25, 26]. The stations in the direction of the rupture propagation exhibit higher frequency ground motion represented by a shorter STF pulse duration and larger displacement amplitude (as the area under the pulse is constant and scales with the scalar seismic moment M₀). Whereas stations in the linearly opposite direction display a lower frequency content (longer STF pulse duration) and smaller displacement amplitude. This enhancement or depletion in frequency content in the first P-wave pulse can be ascertained by modeling the source spectra [27] and using the corner frequency variation as a function of azimuth [10]. This analysis works best on near-vertical strike-slip earthquakes, where the strikes of the nodal planes are nearly orthogonal and directivity effects will have a distinctive azimuthal pattern. We adopt this technique to distinguish the fault plane of the Manipur Earthquake (event f in Figure 3(a)). We use seismograms from stations within 30° epicentral distance. The seismograms are corrected for instrument response and transformed to displacement traces. A quality check of the data is performed based on the signal-to-noise (SNR) ratio analysis. A signal window, containing the direct P-wave arrival, and a pre-signal noise of the same time duration is manually picked on the vertical component. The amplitude spectra for both the time series are computed by fast-Fourier transformation. The SNR is calculated by taking the ratio of the RMS spectral amplitude of the signal and the noise. P-wave windows with SNR > 2 are selected for the analysis. Data from 63 stations pass the quality check and contribute to a good azimuthal coverage (Figure 4). The displacement spectrum for the P-wave is computed and fit using the Brune source model [27] to estimate the corner frequency (f₀) at each station. For further details of the analysis, the reader is referred to [10]. The distribution of f₀ shows a strong azimuthal dependence (Figure 4(a)). Two clusters of low and high f₀ are observed almost 180°–200° apart from each other. The first cluster segment from azimuth 107° to 180° contains the station Qiongzhong (station code- QIZ) which exhibits the lowest value of f₀ (Figure 4(a)). The second cluster in the azimuthal range of 300°–340° has station Shaartuz (station code-SHAA) with the highest f₀ value. The azimuth of QIZ is 107° and SHAA is 306° and the two are 199° apart. Therefore, the NW-SE-oriented nodal plane is the fault plane of the Manipur earthquake. The plot of P-wave spectra, normalized to the maximum amplitude, for the stations QIZ and SHAA confirms that SHAA has a relatively higher amplitude than QIZ (Figure 4(b)). One notable feature about the f₀ distribution is its bimodal high f₀ along both these directions, which occurs in bilateral ruptures (Figure 4(a)). This aspect of the rupture is further explored in the next section.

In order to further corroborate the identification of the fault plane and to study the nature of the rupture, we compute the spatio-temporal evolution of rupture for the Manipur earthquake. The time-reversal property of seismic waves [28-35] enables us to back project the energy onto a source plane at the location and depth of the hypocenter. Waveform data from broadband seismic arrays at different distances and azimuths are used to ascertain the source location and rupture evolution [36]. The P-wave group arrivals are corrected for three-dimensional effects, its energy is stacked as a function of time and back-projected onto a potential source grid [37]. We use vertical-component data from 78 stations of the Australian Network and 326 stations from the European Network (Figure 5(a)) for the analysis. These arrays of stations are within the epicentral distance range of 30°–90°. The vertical component seismograms are decimated to 20 samples per second (sps) and filtered in the frequency range of 0.2–5 Hz. The chosen frequency range allows us to sample the initiation of the rupture from the low-frequency energy. Additionally, the high-frequency information highlights the detail of the spatio-temporal variation occurring during the rupture propagation. The back-projection analysis has been performed following the methodology of Kiser and Ishii [36]. We set up the potential source grid points in two-dimensions keeping the hypocentral depth fixed at 52 km (obtained from source mechanism modeling result) in order to better resolve the laterally varying behavior. Thus, the potential source is represented by a plane parameterized with 0.1° × 0.1° square grids. Each grid point on the parameterized hypocentral plane is assumed to be a potential source and the P-phase theoretical travel-time is calculated using the ak135 velocity model [38]. The teleseismic P-waveforms from each station are windowed starting from 60 seconds prior to the calculated direct P arrival to 120 seconds after the phase arrival. For both the arrays, the time shift in P-wave arrival time generated by lateral variation in the velocity structure is estimated via cross-correlation. The windowed P-waveform for every station is cross-correlated with a station lying at the center of the array (hereafter referred to as the reference station). Next, the waveforms are weighted by their respective cross-correlation coefficients and the inverse of the density of stations within an array. The weighted stacked energy is back-projected onto the source grid. We repeat this procedure for all potential source grids, for each array. Results from our analysis suggest that rupture occurred over a duration of ~ 37 seconds. Though the STF consists of two distinct pulses at ~ 5 seconds and ~ 32 seconds, the majority of energy is released in the first ~10 seconds (Figure 5(c)). The sharp rise of the STF marks the initiation of the rupture. The shape of the low-amplitude second peak near ~32 seconds indicates a small release of energy for a duration of ~10 seconds. The energy distribution pattern on the hypocentral plane suggests that as the rupture initiates at the focus it first propagates to the SE and then to the NW, revealing a bilateral rupture (Figure 5(b)).

We compute the horizontal strain-rate field and rotation-rate field for this region following the technique given by Howell [39] and using the GPS-velocity field of Vernant [3] (in the India-fixed reference frame) (Figure 1). In order to do this, we separate the horizontal velocity gradient tensor into its symmetric and antisymmetric components (the strain-rate and rotation-rate tensors), where:

In the above equation, u_i is displacement, $δui⋅δxj$ is the velocity-gradient tensor, and ε̇_ij and α̇_ij are the strain-rate and rotation-rate tensors, respectively, and the rotation-rate tensor α̇_ij is half the vorticity tensor ω̇_ij. The GPS-velocity field which is shown by the red arrows (Figure 6) was surfaced using splines and then smoothed with a Gaussian filter, using a range of tension factors between 0 and 1 during the surfacing process. We tested diameters between 30 and 360 km for the Gaussian filter. It is seen that changing the tension factor only changes the smoothness of the velocity field and not its overall nature; therefore, it does not really affect the computed strain-rate and rotation-rate fields. Changing the Gaussian width has important effects because the Gaussian smoothing affects the magnitudes of calculated strain rates. If the width is smaller than 100 km, the strain-rate field becomes dominated by differences in velocity between adjacent sites, which are often smaller than the errors at those sites; these differences in velocity can predict small compressional strains in regions of extension if the width is too small. We used a tension factor of 0.5 and a Gaussian width of 150 km in our calculations.

We discuss the observed focal mechanisms vis-a-vis the strain field (Figure 6). The magnitudes of the principal strains of the horizontal strain-rate field show that the strain rate decreases from north to south from the Himalayan mountain-front to the Brahmaputra Valley and finally the Bengal Basin. The orientations of the principal compressive strain axes changes from N-S in the Himalayan mountains to nearly E-W along the Bengal Basin. Seismicity in the regions north of the Shillong Plateau is consistent with the observed strain-rate field as the earthquakes there require the pressure axes to be oriented grossly in a N-S direction. The earthquakes in the Bengal Basin also indicate N-S-oriented pressure-axes but the strain-rate field shows the orientations of the principal compressive strain axes to be grossly E-W. This possibly implies that N20°E-directed compression from the India-Eurasia collision drives the deformation within the seismogenic continental crust (north of the Dawki Fault) and the transitional crust (beneath the Bengal Basin sediments). However, the surface deformation in the Bengal Basin is dominated by E-W compressive stresses arising from the shortening across the Indo-Burman Convergence Zone and is possibly decoupled from the deformation of the basement.

One of the goals of this article is to understand how the strike-slip and oblique-thrust faulting across NE India accommodate the convergence in this intraplate zone of diffused deformation. This zone lies immediately south of the Eastern Himalaya, which is presently accommodating the N20°E convergence between India and Eurasia; and west of the Indo-Burman convergence zone. None of the earthquakes studied here are associated with either of the two plate boundaries and occur within the mid-to-lower crust of the Indian Plate [9]. This reveals that the entire Indian crust in NE India is seismogenic. This is similar to the rest of the Indian subcontinental crust, which is conjectured to be composed dominantly of dry granulites with temperature colder than 600 °C [40].

To discuss the kinematics of faulting, it is necessary to identify the fault plane from among the two nodal planes. We have done this for the larger magnitude event by using rupture directivity and back-projection. We use this result along with the nearby active KFZ and the DCF zone to identify the fault plane of the earthquakes in the region to the north of the Shillong Plateau. Only in the case of the pure dip-slip faulting (event k) (Figure 3), the direction of motion between the two plates is insensitive to the choice of fault plane. From GPS data, it has been shown that the Kopili Fault has segmented NE India into the Assam and Shillong blocks, which both undergo clockwise rotation in the India fixed frame, with convergence of 3–6 mm yr⁻¹ across the Dawki Fault [3]. The Assam block, in the east, rotates faster than the Shillong block, to its west, and results in dextral strike-slip motion along the NW-SE trending KFZ. The dextral DCF zone which extends into the Himalayan foreland in the Brahmaputra Valley has been conjectured to be the western margin of the Shillong Block [10, 11]. We choose the NW-SE-oriented planes of the strike-slip and oblique-thrust earthquakes (events a, b, c, d, e, k, l, g, h, i, and j) (Figure 3) to be the fault planes. These earthquakes originated in the Brahmaputra Valley and the Mikir Hills, in the vicinity of the DCF and the KFZ, and have similar mechanisms and depth distribution. Hence, we hypothesize that the region between KFZ and DCF is traversed by a set of subparallel faults within the Indian Plate which have dextral strike-slip motion [7]. The EHZ, which runs through the Bengal Basin in a NE-SW direction, is known to be the passive continental margin of the Indian continent during Eocene time [9]. Such continental margins are associated with paleo-rift normal faults which dip in the basin-ward direction. These paleo-rifts represent preexisting weak planes present in the crust of the Bengal Basin which is overlain by a thick pile of sediments (⁠$∼$20 km). These rifts possibly got reactivated in response to the stress generated by the northward convergence of the Indian plate and manifest as strike-slip and oblique-thrust earthquakes (events m, n, o, p, q, r, s, t, and u) (Figure 3). From this understanding, we choose the NE-SW striking nodal planes to be the fault planes as they are almost subparallel to the hinge zone and accordingly, they exhibit sinistral motion.

The fact that the Indian basement has preexisting transverse structures which often extend into the Himalaya and Tibet had been recognized as early as 1976 by Valdiya [41]. He conjectured that such faults in the Himalaya are broadly parallel to the great faults and ridges in the basement of the Ganga Basin. Dasgupta [42-44] proposed that the Himalaya is traversed by conjugate shear fractures in many regions and movement across these fractures is capable of generating large earthquakes, for instance, the M_w 6.7 1988 Udaypur earthquake and later on, the M_w 6.9 2011 Sikkim earthquake [10, 11]. The 1988 event was presumed to have occurred due to reactivation of the East Patna Fault underneath the foothills of the Nepal Himalaya [44]. Some of these transverse structures are known to have structurally segmented the Himalayan arc. For example, the Gish Transverse fault (GTF), which is a NNE-SSW trending sinistral strike-slip fault in the Sikkim Himalaya, is an extension of the Kishanganj fault at the eastern boundary of the Munger-Sahasra ridge of the Ganga Basin. The GTF further continues into the Yadong-Gulu rift to the north [45-48]. Similarly, Duvall [49] used data from industry seismic lines to map NE-striking strike-slip faults in the foreland of eastern Nepal. These faults, several of which bound the Munger-Saharsa ridge, penetrate from the basement up through the Quaternary Upper Siwalik lithounits. Recent studies [50-52] have shown that these basement faults correlate spatially with lateral variations in pressure–temperature–time evolution of the Greater Himalaya, lateral ramps in the Himalayan fold-thrust belt, and the location of graben in southern Tibet. Two key conclusions are shared by all of these studies. First, these faults and terrane boundaries are typically at a high angle to the regional strike of the Himalaya and are confined to the underthrusting Indian crust. Second, based on the orientation of these transverse structures with respect to the present-day stress field, these are reactivated as oblique-slip to strike-slip faults, beneath the Himalayan foreland basin. This is also evident from our present study in NE India.

In Figure 3(b), we have compiled the horizontal projections of the slip vectors of most of the earthquakes based on our choice of the fault plane. In the region North of the Shillong Plateau, the slip vectors dominantly point in the S to SE direction and in most instances show a close match with the GPS velocity vectors. In the region South of the Shillong Plateau, the slip vectors point in the SW direction and are also subparallel to the GPS vectors at most of the sites (Figure 6). We therefore envision a tectonic model where there are two sets of transverse faults, one oriented NW-SE in the region north of the Shillong Plateau, and the other oriented NE-SW in the Bengal Basin. These two sets of faults accommodate the shortening within this intraplate diffused deformation zone [7]. Both sets of faults and their slip vectors are oblique to the overall strike of the deforming zone. In this situation, these two sets of faults function as a “shortening array” as termed by Allen [53], where deforming blocks rotate about a vertical axis. If the deformation is entirely taken up by pure strike-slip faulting on both sets of en echelon faults, the shortening across this zone would most likely be accompanied by lateral extrusion of material parallel to the deforming zone, in order to preserve the surface area (Figure 3(a) and (b) of Jackson and Molnar [54]). However, a number of these earthquakes have oblique-slip (thrust) faulting. This implicitly indicates that there is also a component of crustal thickening (thrusting) involved. This means that oblique convergence across this zone must indeed be accommodated in its entirety by very small rotations of rigid blocks across the deforming zone but in a manner similar to the “pinned” block model of McKenzie and Jackson [55, 56]. In this model (Figure 3(c) and (d) of Jackson and Molnar [54]), where the edges of the blocks are pinned to the margins of the zone of deformation, the loss of surface area is accounted for by crustal thickening or thrusting in between the blocks. The overall convergence in this case is partitioned into two broad components: “W × a,” the component which is parallel to the deforming zone and “2 × T × a,” the component perpendicular to it. Here “a” is the width of the zone of deformation, “W” is the rotation rate of the blocks, and “2 × T” is the rate of crustal thickening. Therefore, in response to the N20°E direction of the India-Asia convergence, the blocks in the region north of the plateau rotate counterclockwise and exhibit dextral motion and those in the Bengal Basin rotate clockwise and show sinistral motion. This represents the simplest fault geometry which is capable of absorbing the distributed deformation (Figure 7). Similar two-dimensional block models of strike-slip faulting accommodating shortening across other broad deforming regions of the Alpine Belt (in Turkey, Iran, and Greece) are known [56-58] and hence adopted for this region.

The tectonics of NE India, though primarily guided by the collision between India and Eurasia, is complicated by the deformation surrounding the Shillong Plateau and the eastward subduction of the Indian Plate beneath the IBR. In this section, we provide a plausible framework to explore the forces. The earthquakes in the regions, north and south of the Shillong Plateau, have strike-slip, oblique-thrust, and thrust motions. From the fault plane solutions, it can be inferred that the pressure axes of all these events in the Assam, Brahmaputra Valley as well as those in the Bengal Basin are oriented predominantly in the N-S direction. This implies that the earthquake faulting is driven by the compressive forces from N-S convergence of India and Eurasia. Of particular interest here is the thrust fault earthquake in the Brahmaputra Valley having a depth of $∼$25 km. This is shallower than the thrust fault earthquakes beneath the Western Ganga Basin [59-61], where the Indian continental crust bends as a single elastic layer in response to the loading from the Himalayan Mountains. The upper crust is under extension (with normal fault earthquakes) and the lower crust is under compression (with reverse faulting earthquakes). The neutral surface lies approximately at a depth of 25 ± 5 km. This model is well supported by observations of a small number of well-recorded M_w > 5 normal-faulting earthquakes having nodal planes parallel to mountain fronts between the surface and 20 km depth and reverse faulting earthquakes with depths of 30–50 km in the forelands of Himalaya and the Andes [62-64]. A similar pattern has been seen in the earthquakes in the outer trench rise and outer trench slope of subduction zones [65]. However, the shallow thrust earthquake beneath NE India indicates that the neutral fiber is at a shallower level. This can be achieved by in-plane compressive stresses [65]. Alternatively, lower in-plane compressive forces, for instance, in the foreland of the Algerian Atlas Mountains serve to shift the neutral fiber deeper and results in normal faulting earthquakes down to a depth of 31 km [66]. We infer that the shallowing of the neutral fiber in NE India could be due to increased stresses at the interface of the plates between India and Eurasia. Furthermore, it is also possible that the loading by the Shillong Plateau may play a role in altering the depth of the neutral surface.

It is worth noting that with N-S-oriented pressure axes (i.e. N-S-directed compression), we would expect dip-slip faults with fault planes striking roughly E-W. However, with the exception of the one shallow thrust earthquake (event k, $∼$25 km depth) that we have mentioned earlier, all the other events are strike-slip or oblique-thrust in character. It is therefore inferred that these fault planes are possibly inherent weak structures in the underthrusting Indian Plate as suggested earlier. In the region north of the Shillong Plateau, these structures are oriented roughly NW-SE following the orientation of the Kopili Fault and the Dhubri-Chungthang Fault. Beneath the Bengal Basin, south of the Shillong Plateau, the earthquakes occur on NE-SW-oriented reactivated paleo-rifts. Hence, the N-S compressive stresses between India and Eurasia cause spatially distributed brittle deformation of the seismogenic layer. South of the Bhutan Himalaya, the Himalayan foreland basin has a thickness of less than 1 km [9]. Therefore, the crystalline basement of the underthrusting Indian crust is shallower than regions to its west, and the width of the foreland basin is no more than 30 km (Figure 1). This is consistent with the fact that from west-to-east there is a decrease in the role of flexure and the flexural wavelength as was observed by Jordan and Watts [67]. From such a line of reasoning, we infer that stress state in the lithosphere of this region is dominated by the net compressional forces relating to the India-Asia collision over flexural effects. This is also corroborated by the strain-rate field to the north of the Shillong Plateau, with the roughly N-S trending principal compressive strain axes (Figure 6). Additionally, we cannot rule out the possibility that the shallowing of the crystalline basement in NE India can also be due to the thrust/crustal thickening/uplift component associated with the movement along KFZ, DCF, and the other faults.

The presence of the uplifted Shillong Plateau is this region is also a testament to the high in-plane compressional forces across this region. Possibly, the only analogue to this kind of an uplifted basement structure are the Laramide structures of the Rocky Mountains in Wyoming, US, though those structures are much smaller in scale compared to the Shillong Plateau. Formation of the Laramide structures involved considerable crustal shortening and horizontal compression which caused a thrust uplift in which the crust behaved as a rigid plate [68, 69]. The uplift of the Shillong Plateau was driven by the localized reactivation of a paleo-rift margin at its southern edge which was caused by the continued northward movement of the Indian Plate. The southern boundary of the plateau (essentially the Dawki Fault) juxtaposes normal continental crust to its north against attenuated (or transitional) continental and oceanic crust in the south [2, 9].

The state of stress south of the Shillong Plateau at the surface and at depth is more complex as the eastward subduction underneath the IBR probably plays an active role in the deformation. The IBR segment of the plate boundary is associated with the E-W convergence and oblique subduction of the Indian Plate beneath the Burma sliver plate at about $∼$46 mm yr⁻¹ [70]. It is estimated that $∼$13–17 mm yr⁻¹ of plate convergence takes place on an active, shallow dipping and locked megathrust below the IBR [70]. The IBR are the surface expression of a very wide accretionary prism or wedge comprising Paleogene-Neogene sediments of the Ganges delta and the Bengal Basin. Maurin and Rangin [71] have proposed that the Indo-Burmese Wedge exhibits diffuse strain partitioning with dextral shearing in the inner wedge and E-W shortening in the outer wedge. This outer wedge is characterized by roughly N-S trending folds and thrusts, anticlinal ridges, and synclinal valleys [72]. The events f, v, w, and x are within the outer wedge and the events y and z are within the inner wedge. Event f has clearly been demonstrated to have occurred through dextral strike-slip faulting. The other events are possibly of similar character. Event f may also be related to the KFZ which would mean that the KFZ extends far beyond its mapped limit, throughout the Indian crust, as far as the IBR [7]. Because of the diffused nature of the strain partitioning, the region in the outer wedge should undergo both E-W shortening and N-S strike-slip deformation, according to Maurin and Rangin [71]. The dextral motion across these regions is concentrated on several strike-slip to oblique-slip faults. In terms of forces, we can infer that the slab pull from the downgoing Indian Plate is indeed a dominant driving force which was suggested by Steckler [70] as the focal mechanisms of all these earthquakes indicate E-W trending tension axes (downdip tension). The slab-pull force deforms the slab and drives its motion in the trench-normal direction. This along with the strike-parallel components of the inter-plate coupling resistance and the mantle drag forces sets up a lateral shear within the downgoing slab as it descends. This lateral shear is responsible for the strike-slip earthquakes at shallow to midcrustal depths within the slab [73, 74].

The Bengal Basin is characterized by a thin “transitional” or “attenuated” crust, which has many rift faults and this crystalline crust is overlain by sediments whose thickness increases progressively southward (from 12 to 14 km in the north to 18–20 km in the south) [9]. As mentioned in the previous sections, the earthquakes in this region occur within a depth range of 20–40 km, lying within the crystalline basement layer and occur on NE-SW trending fault planes which are conjectured to be the paleo-rifts (parallel or subparallel to the EHZ). These are reactivated as a result of the N-S-directed compressive forces from the Indo-Eurasia collision. Hence, it is clear that the earthquake faulting throughout the region both north and south of the Shillong Plateau is driven by the N-S compression. However, in the Bengal Basin, the surface deformation (as seen from the strain-rate field) (Figure 6) is dominantly E-W shortening and it is probably decoupled from the earthquake faulting at depth. The E-W-oriented compression in the outer wedge is seen to extend throughout the Bengal Basin and it is most likely developed due to the resistance forces that manifest at subduction zones, which resist the relative movement between the two plates [75]. The Tripura fold-thrust belt in the outer wedge is primarily the result of this east–west directed compressive stresses which caused a combination of buckling and fault-propagation folding. This gave rise to the eastward dipping thrusts that separate the successive accretionary wedges [76]. It is interesting to note that this rather sharp change in the deformation field occurs almost exactly across the Shillong Plateau. In summary, the N20°E motion of the Indian Plate with respect to Eurasia is partitioned into the following broad components—a N-S compression which drives the earthquake faulting north and south of the Shillong Plateau, which also manifests as dextral shear in the IBR, and an E-W-oriented compression south of the Shillong Plateau which drives the surface deformation in the outer accretionary wedge in the form of the fold-thrust belt. Finally, there is also the slab pull which contributes to the deformation beneath the IBR causing the Indian Plate to rupture along preexisting weak structures due to tensional stresses.

We study the source mechanism of moderate-to-large intraplate earthquakes in North-East India to understand their role in accommodating the ongoing shortening across this wide zone of distributed deformation. All earthquakes originate in the mid-to-lower crust of the Indian Plate, which reveals cold brittle seismogenic Indian crust. These events continue beneath the Eastern Himalaya and the Indo-Burman convergence zone due to the ongoing underthrusting and subduction of the Indian plate, northward and eastward, respectively. For the largest event we use rupture directivity to identify the fault plane from the auxiliary plane, and through back-projection of high-frequency energy we investigate the spatio-temporal evolution of the rupture. These results along with the GPS derived strain-rate field are used to unravel the seismotectonics. The N20°E motion of India with respect to Eurasia occurs across this wide zone of NE India and is partitioned into two components. The first is the N-S-directed compression, which results in earthquake faulting on transverse structures north and south of the Shillong Plateau. Crustal blocks north of the plateau rotate counterclockwise and exhibit dextral motion along NW-SE-oriented faults; and those in the Bengal Basin rotate clockwise and show sinistral motion along NE-SW trending faults. Also high in-plane compressive stresses in the underthrusting Indian Plate are observed in this zone. The second is the E-W-oriented compression which is responsible for the surface deformation in the outer wedge of the IBR, and the Bengal Basin (in the form of the Tripura fold-thrust belt). Slab-pull force from the subducting Indian Plate contributes to the deformation in the IBR by causing the plate to rupture under tensional stresses.

All global seismic data were downloaded through the EarthScope Consortium Web Services (https://service.iris.edu/). Data preprocessing and part of the analysis were performed using Seismic Analysis Code, version 101.6a [77]. All plots were made using the Generic Mapping Tools version 5.0 [78].

The authors declare that there is no conflict of interest regarding the publication of this paper.

This research did not receive any specific grant from funding agencies.

D.C. and R.B. acknowledge Department of Science and Technology (DST), India - Innovation in Science Pursuit for Inspired Research (INSPIRE) for PhD Fellowship. S.M. acknowledges IISER Kolkata Academic Research Fund (ARF). We acknowledge the assistance from Sankha Subhra Mahanti with the moment tensor inversion. We acknowledge Alex Copley for his help with the strain-rate computation and his insightful comments. We thank James Jackson for his valuable inputs. We sincerely acknowledge the critical reviews by two anonymous reviewers and the section editor, which significantly improved the quality of the manuscript.

The best-fitting solutions for the focal mechanisms modelled (events b, c, d, e and f) in this study can be found in the file, online supplementary S1.