Crustal Configuration and Seismic Stability of the Eastern Indian Shield and Adjoining Regions: Insights for Incidents of Great Earthquakes in the Nepal-Bihar-Sikkim Himalaya

The Eastern Indian Shield (EIS) is comprised of the intracratonic (coal-bearing) Damodar Gondwana basin, rift-controlled extensional Lower Gangetic basin (LGB), and the downward flexed Indo-Gangetic basin (IGB). The present study involves the computations and mapping of the basement configuration, sediment thickness, Moho depth, and the residual isostatic gravity anomaly, based on 2-D gravity modeling. The sediment thickness in the area ranges between 0.0 and 6.5 km, and the Conrad discontinuity occurs at ∼17.0–20 km depth. The depth of the Moho varies between 36.0 and 41.5 km, with the maximum value beneath the Upper Gangetic basin (UGB), and the minimum of ∼36 km (uplifted Moho) in the southeastern part beneath the LGB. The maximum residual isostatic anomaly of +44 mGal in the southern part indicates the Singhbhum shear zone, LGB, and Rajmahal trap to be under-compensated, whereas the northern part recording the minimum residual isostatic anomaly of –87.0 mGal is over-compensated. Although the region experienced a few moderate-magnitude earthquakes in the past, small-magnitude earthquakes are sparsely distributed. The basement reactivation was possibly associated with a few events of magnitudes more than 4.0. Toward the south, in the Bay of Bengal (BOB), seismic activities of moderate size and shallow origin are confined between the aseismic 85 and 90°E ridges. The regions on the extreme north and south [along the Himalaya and the equatorial Indian Ocean (EIO)] are experienced moderate-to-great earthquakes over different times in the historical past, but the intervening EIS and the BOB have seismic stability. We propose that the two aseismic ridges are guiding the lithospheric stress fields, which are being further focused by the basement of the EIS, the BOB, and the N-S extended regional fault systems into the bending zone of the penetrating Indian lithosphere beneath the Himalaya. The minimum obliquity of the Indian plate and the transecting fault systems in the Foothills of the Himalaya channelize and enhance the stress field into the bending zone. The enhanced stress generates great earthquakes in the Nepal-Bihar-Sikkim Himalaya, and on being reflected back through the apparently stable EIS and BOB, the stress field creates deformation and great earthquakes in the EIO.

The present study was undertaken at the core of the EIS, which occupies portions of four Indian states-Bihar, Jharkhand, Orissa, and West Bengal-to investigate the link between the regional geodynamic processes and earthquake occurrences. Besides the tectonic features mentioned above, the eastern part of the EIS is associated with the Rajmahal Trap and several faults striking N-S (e.g., Pingla fault, Sainthia Bahamani fault, and GKGF faults). The core (central part) of the EIS has ∼ E-W striking regional shear zones such as the North Purulia, South Purulia, and Singhbhum (Figure 2). Thus, the area is tectonically very important for its location, and is essentially a geological corridor to the Eastern Himalaya and the Indo-Myanmar Ranges. The present study therefore was carried out to delineate the basement configuration and the Conrad and Mohorovičić (Moho) discontinuities based on analysis of the gravity anomalies. The detailed gravity modeling and mapping and the historical earthquake data available for the region were used for exploring the level of seismic stability at different locations of the study region.
The Moho depth for several small blocks of the EIS was determined by several studies (Choudhury, 1974;Verma et al., 1976Verma et al., , 1980Bhattacharya and Shalivahan, 2002;Singh et al., 2004Singh et al., , 2015Kayal et al., 2011;Sharma et al., 2015). All these investigations were constrained by an insufficient coverage of dataset and smaller sizes of the areas. The outcomes of these investigations provide valuable information regarding the crustal structure (discussed below).
The residual gravity analysis over the Raniganj Coalfield area (∼6 × 10 3 km 2 ) shows a gravity low of around -32 mGa1 associated with the Gondwana sediments, with a maximum estimated sediment thickness of ∼2 km (Verma et al., 1976). Gravity field analysis of North Singhbhum (∼36 × 10 3 km 2 area) could delineate the Bouguer anomalies of +4 to −62 mGal, with the estimated thickness of the basin, based on a 2D gravity model, of about 6-7 km (Verma et al., 1980). An analysis of free air, Bouguer, and isostatic anomalies over a part of the shield near the Gaya region shows a 6.6 km thick sediment layer underneath the Lower Gangetic basin (LGB). Compiling datasets from a few stations, a crustal thickness of 42.5 km underneath the northern part of the shield, 31.5 km underneath the Indo-Gangetic basin (IGB), and 81 km underneath the Higher Himalaya of Nepal-Bihar region were also estimated by Qureshy (1969). While the Bouguer gravity analysis, carried out by Choudhury (1974), for the Indo-Gangetic Plain and the northern Himalayan region, shows a crustal thickness of ∼35-37 and ∼70-72 km, respectively. A magnetotelluric study reports 40 km deeper Moho beneath the Singhbhum granite and adjoining regions (Bhattacharya and Shalivahan, 2002). Receiver function analysis of P-wave gives a depth of ∼41-42 km for the Moho and ∼20 km for the Conrad discontinuity below the IIT (ISM), Dhanbad seismic station (Kayal et al., 2011;Sharma et al., 2015). A 2D gravity modeling study shows that the Moho depth varies from ∼38 km below the Rajmahal Traps to ∼40 km near the Raniganj region (Singh et al., 2004). Lithospheric 3D mapping of the EIS shows that the Moho ranges from ∼35 to 40 km near Mahanadi, Damodar, and LGB . For a regional scale, Bouguer gravity anomalies can reflect the changes in mass disseminations between the lower crust and the upper mantle, and the isostatic stability (Tealeb and Riad, 1986). Bouguer anomalies and surface relief are also closely connected with the crustal thickness (Coron, 1969;Woollard, 1969;Pick et al., 1973;Riad et al., 1983;Riad and El-Etr, 1985). In the present study, the Bouguer gravity anomaly data and the geophysical parameters of various studies (Choudhury, 1974;Verma and Ghosh, 1977;Verma, 1985;Singh et al., 2004) were used to deduce the regional structural framework of the shallow basement, overlying sedimentary cover, and the Conrad and Moho discontinuities. The present study also attempts to visualize the crustal balance in terms of isostasy. Isostatic study of continents accounts for the mechanical balance of lithosphere at the time of loading (Watts et al., 1980;Djomani et al., 1992) and links to the deep-density structure, providing constraints on the geological model (Simpson et al., 1986).

TECTONIC SET-UP
The study area comprises the northern portion of the Singhbhum craton, the Chhotanagpur gneissic complex, and the southern portion of the Himalayas (Figure 1). The area is bounded by the Eastern Ghats mobile belt on the south and the Himalayan foothills on the north. The Gangetic alluvium and the FIGURE 1 | Regional map (reconstructed after Khan et al., 2017b;Valdiya and Sanwal, 2017;Altenbernd et al., 2020) showing the location of the study area Quaternary sediments of LGB border its eastern part, whereas the western part is bordered by the Gondwana sediments of the Sone-Mahanadi basin. The Upper Jurassic to lower Cretaceous basalts form the Rajmahal Traps, at the contact zone of the Chhotanagpur Gneiss Complex and the LGB. The Ganges River flows all along the north of the study area and was evolved during the bending of the Indian plate beneath the Himalaya (Figures 1, 2). The rivers Subarnarekha and Koel pass through the southern part of the area, whereas the coal-bearing, faulted Damodar Valley transects the central part of the region from west to east.
The study area has variable relief, with the highest elevation in the Ranchi Plateau, situated at ∼610-1097 m from mean sea level (MSL) (Pandey et al., 2013). Several dissected hills with a number of hill stations constitute the adjacent Hazaribagh Plateau drained by the Damodar valley. The Chhotanagpur Plateau and the Hazaribagh Plateau have outcrops of gneisses, schists, and granite of Precambrian rocks (Ghose, 1983). Some other plateaus, like FIGURE 2 | Seismoteconics and topographic relief Map of EIS (EIS) and adjoining regions (after Dasgupta et al., 2000;Singh et al., 2020). The study area is marked by a red rectangular box. Black dashed lines represent faults, and black solid lines denote lineaments. Blue solid lines show the various permanent rivers. Black solid squares represent the locations of a few important areas. A combination of black and white solid arrows shows the stress regimes of that part within shallow depth (0-20 km). RC, radial compression; PSS, pure strike-slip; CSS, compressive strike-slip. Red arrows represent the Indian plate velocity motion direction with respect to the Eurasian plate. All stars show the seismicity since 1973 to May 2020 (from USGS Catalog). Magenta stars show the location of damaging earthquakes with magnitudes of 6.5 < M ≤ 8.0 (1: M 7.8 Nepal, 2: 2015M 7.3 Nepal, 3: 2015M 6.7 Nepal, 4 2015M 6.6 Nepal 5: 2011M 6.9 Sikkim, 6: 1988M 6.9 Nepal-India Border, 7: 1934M 8.0 Nepal-India Border, and 8: 1833. Maroon stars represent the location of historical earthquakes that occurred in EIS and nearby regions within 5 < M ≤ 6. 3 (9: 1969M 5.7 Bankura, 10: 1964M 5.5 Midnapur, 11: 1963M 5.2 Thethanagar, SE Bihar, 12: 1868M 5.5 Manbhum, 13: 1868M 5.0 Hazaribagh, 14: 1866M 6.3 Bengal, and 15: 1861. Blue stars show the epicenter location of earthquakes between 6 and 6.5 magnitude, red stars represent earthquake location within 5 ≤ M < 6, while the black stars denote epicenter locations of earthquakes within 4 ≤ M < 5. Kaimur Plateau and Rajmahal Hills, can be differentiated from each other by steep and narrow slopes known as scarps. Mutual cyclic interactions of two converging microplates, i.e., Singhbhum and Chhotanagpur (Sarkar, 1982), during the Proterozoic times facilitated subduction, rifting, mantle plume activities, and mafic magmatism in the region (Pandey and Agarwal, 1999;Bose, 2008). The cyclic interactions happened in three consecutive stages with intervening quiescence (Sarkar, 1982), which further simplified the isostatic adjustments with coeval upliftment with intrusions of basic dyke swarms, erosion, and paralic sedimentation (Bose, 1999(Bose, , 2009. During the Precambrian, the northward movement of the Singhbhum block toward the continental margin of the Chhotanagpur block produced a consuming plate boundary (Figure 1). The compressional stress produced by the plate merging caused N-S crustal shortening, resulting in the upliftment on the Chhotanagpur granite-gneissic terrain and beginning the basement reactivation and fresh rupture development (Sarkar and Saha, 1977).
The Chhotanagpur peneplain has been uplifted three times during the Cenozoic due to the flexural responses caused by the subduction of the Indian lithosphere toward the southern edge of the Himalayas, and gave rise to some familiar waterfalls, like Jonha and Hundru, on the scarps. During the Eocene to Oligocene period the first upliftment took place, forming the Pat region; the second one happened during the Miocene period which developed the Ranchi and Hazaribagh Plateau; and the third one took place during the Pliocene and Pleistocene period, which caused upliftment of the outer Chhotanagpur Plateau (Naqvi and Rogers, 1987;Weaver, 1990).
The Chhotanagpur plateau is comprised of high-grade gneisses and migmatites containing enclaves of metasediments, and mafic and ultramafic rocks. The basement of the Chhotanagpur plateau is represented by an uncategorized gneissic complex. Systems of Proterozoic fold belt occur as secluded patches within the Chhotanagpur Granite Gneissic terrain. Later, during late Palaeozoic and Mesozoic, Gondwana sediments were deposited in linear grabens in the Damodar valley. Irregular sedimentation above crystalline metamorphic rocks, harmonized subsidence, and sprinkled periodic accumulation of coal measures were greatly conditioned by climatic variation in the period of sedimentation in the Gondwana basin in the eastern part of the EIS (Sarkar, 1982). The maximum thickness of sediment filling in the basin is estimated to be 2.9 km in the central part (Verma et al., 1980). Intermittent periods of faulting and epeirogenic uplift of the margins of the master basin dismembered the master basin into several smaller sub-basins (Mahadevan, 2002). Various hot springs with numerous features are broadly dispersed over this area, and the area is apparently associated with several fault zones (Mahadevan, 2002). The Chhotanagpur Granite-Gneiss Complex (CGGC) mainly consists of granitic gneisses, migmatized porphyritic granite, and many metasedimentary enclaves. The Central CGGC gneisses and the western CGGC granites show almost similar ages of 1.7 Ga, whereas the migmatites and granulites, and granitic gneisses of northeastern CGGC (Dumka area), show ages of 1.5-1.6 Ga (Chatterjee et al., 2008).

Gravity Modeling and Validation
Gravity modeling is an important tool for delineation of crustal structures. Gravity data show an ordered link between crustal structure, density, and the surface altitude. In this line, a Bouguer anomaly map was reconstructed for the area covering a part of the Gangetic basin and EIS (i.e., between 22.5-25.5 • N latitude and 84.5-88 • E longitude) compiling the data from the second generation gravity map series of India (GMSI, 2006). An inaccuracy of about ± 0.5 mGal for the errors in elevation of 2-3 m was achieved in these data. The area of the new Bouguer anomaly map (Figure 3) has been divided into 99 square-grids, with a dimension of 1 • × 1 • . A shifting window with 75% overlapping area is used for inherent continuity of the data points and better resolution for delineating the basement and the Conrad discontinuity of the study area.
The Moho depth was computed by the spectral analysis (Figure 4) (Odegard and Berg, 1965;Bhattacharyya and Leu, 1975) over 12 blocks, by dividing the Bouguer anomaly map into smaller areas of dimension of 2 • × 2 • in an overlapping adjacent grid. We have also computed the Moho depth based on the Airy-Heiskanen model at the same grid for better resolution and insufficient data points, as the GMSI (2006) gravity map is created from a grid of dimension 5.5 × 5.5 km 2 . The 2D inverse gravity modeling was performed by adding layers incorporating its prior information as explained in Figure 5 and Table 1. Error minimization between the observed and calculated gravity anomalies was followed at various steps of the iteration processes, to optimize the gravity model.
The inverse gravity modeling using Geosoft was done along 13 profiles, considering various Earth models (Choudhury, 1974;Verma et al., 1976;Verma and Ghosh, 1977;Verma, 1985;Singh et al., 2004). Power spectrum analysis was carried out to compute depths of the sharp density discontinuities appearing at the basement, Conrad, and Moho. The efficacy of this method has been demonstrated in the literatures (Spector and Grant, 1970;Green, 1972;Curtis and Jain, 1975;Tselentis et al., 1988;Djomani et al., 1992;Chavez et al., 1999;Ates and Kearey, 2000;Lefort and Agarwal, 2002;Rivero et al., 2002;Carbó et al., 2003;Gomez-Ortiz et al., 2005;Bansal et al., 2006). We have computed the crustal discontinuities based on the spectral analysis (Figure 4) of the Bouguer gravity field. The Moho depth ( Table 2) has been computed from the slope of the energy spectrum at the lower end of the wave-number band, while the basement depth (Table 3) has been computed at a slightly higher wave-number band. Interface depths with different density contrasts are computed using the following equation: where h is the average interface depth and E 1 and E 2 are the power in the spectrum at wave-number (cycles/km) k 1 FIGURE 3 | Bouguer gravity anomaly map of the study area reconstructed from GMSI (2006). White vertical solid lines are the profiles, taken for the 2D gravity modeling. H-1 and H-2 represent the area associated with high Bouguer gravity anomalies, and L-1 represents the area with low Bouguer gravity anomaly. and k 2 (cf. Spector and Grant, 1970;Karner and Watts, 1983;Maus and Dimri, 1995;Indriana, 2008;Bansal and Dimri, 2010;. Most of the geophysical methods, particularly the spectral analysis of potential field data, isostasy analysis, inversion of potential field data, receiver function, and seismic tomography, used for modeling and interpretation purposes, have some limitations as well as advantages. 2-D spectral analysis usually delivers spatially averaged depth for an undulated plane, i.e., an interface where density changes sharply, and depends on the gridding size of the data (Tselentis et al., 1988). In frequency domain, gravity anomalies generated from a multilayered source can be separated from each other by their own segments with the help of gravity anomaly spectrum and interpreted for the mean depth of the interface . However, the spectral peak for each buried mass at different depths is not always sharp, and the overlapping of such peaks occasionally introduces ambiguity in the result. Refinement of higher wave-number components from noise is sometimes difficult due to its smaller amplitude, and the information loss can be minimized to < 10% by taking care of the window size of four to five times of the source depth (Regan and Hinze, 1976).
The Airy-Heiskanen isostatic technique results in reliable image gravity data on a continent-wide scale for the purpose of regional geologic interpretations. Isostatic residual gravity further enhances the short-wavelength anomalies produced by bodies at or near the surface and emphasizes the regional fabrics and trends in the gravity field. Isostatic anomaly gives information about the lateral mass distribution within the crust and mantle, having some geological interest (Simpson et al., 1986). However, crustal density variations are not considered for one fixed density for the entire crust. Gravity data alone cannot reveal the density inhomogeneity in the locally compensated crust. A solely isostatic model cannot provide the root geometry for a large continental region. There may be errors of less than 2 mGal in the calculation of isostatic residual anomaly because of the uncertainty present in gravity measurement, terrain correction, and station elevation (Chapin, 1996). The anomalies near the flanks are affected by the anomalies present at its neighborhood and become difficult to distinguish.
Gravity inversion method is well known for the derivation of appropriate anomalies due to Moho deflection (Prasanna et al., 2013), and its useful technique to detect Moho depth especially in a region where the availability of seismic data is insufficient. Moho interface can be delineated with the help of an iterative FIGURE 4 | Spectral analysis plots illustrating the estimation of interface depths at different grids. Wave number is plotted on the X-axis, and the corresponding power spectra of gravity are plotted on the Y -axis. Values of deepest density interface show Moho depth, and the shallowest interface shows sediment thickness.
process (Shin et al., 2007). However, the filtering of wavelength effects of gravity from crustal inhomogeneities and deep-seated sources is difficult and complex. Inhomogeneities in the intracrustal density or noise can create a stability problem in inversion and provide a wrong estimation of depth.
Seismic receiver function is a simple and unambiguous method for the detection of 1D crustal structure beneath the seismic station. Receiver function consists of direct, reflected, and converted phases of seismic wave trains (Burdick and Langston, 1977;Ammon et al., 1990;Das et al., 2019). Converted phases and shallow surface multiples carry average crustal information about the surrounding of the receiver station. However, this method can only be applied to a nearly horizontal layer and works best for low incidence angles of seismic rays. Regionalization of Moho depth (2D crustal structure) estimation needs highly dense seismic networks. The seismic tomography can give the information of velocity structures, seismotectonics, and characterization of the source area for tectonically complex regions (Zhao et al., 1992). It gives a high-resolution image of the earth's interior. It can image the structures with different physical properties. Seismic tomography can also be done using surface waves and full waveforms (Geller, 1991;Iyer and Hirahara, 1993). However, the resolution of seismic tomography depends on the wavelengths of seismic waves. Since a longer wavelength is used, deeper sources can be imaged, but small features cannot be resolved. Tomography results are also non-unique and need a good coverage of seismic stations with a huge amount of earthquake data, and hence are suitable for seismically active regions.
The EIS region is not very seismically active, and is located to the south of the seismically active Himalayan orogen. Thus, a good coverage of seismic network will be essential for delineation of the regional crustal structure of the study area. Ambient noise tomography is also another option for the detection of crustal structures in a seismically stable region like EIS region; however, it needs a dense network of seismic stations. Thus, the present results of gravity modeling were partially validated by the receiver function studies (Mandal and Biswas, 2016), and complied by the gravity study of Singh et al. (2015). Only two broadband seismic stations, located in the Chhotanagpur plateau, have been operative since 2007: Dhanbad [run by IIT (ISM)] and Bokaro [run by National Centre of Seismology (NCS), New Delhi]. The 1D crustal structure beneath Dhanbad and Bokaro stations has been estimated by several researchers using different initial velocity structures through receiver function analysis. Das et al. (2019) found 43 and 44 km crustal thickness beneath Dhanbad and Bokaro stations. Kosarev et al. (2013) found Moho at a   Table 4).
We use the Airy-Heiskanen model for the calculation of initial Moho depth assuming that the crust is isostatically being compensated. The isostatic imbalance and the delineation of many layers in the crust might be showing differences in Moho depth between the present study and results of receiver function studies. It is also found that the crustal structures obtained from seismic and gravimetric method are quite different depending on the geological and geophysical hypotheses, as well as different datasets (Reguzzoni et al., 2013). For example, well-known receiver function analysis starts their initial model by assuming a constant velocity and Vp/Vs, taking either a single or double layer of the crust for the calculation of 1D crustal structure (Dueker and Sheehan, 1997;Zhu, 2000;Zhu and Kanamori, 2000;Thompson et al., 2010;Kayal et al., 2011;Sharma et al., 2015;Mandal and Biswas, 2016;Song et al., 2017;Das et al., 2019). An incorrect or different initial seismic model may produce different or inaccurate crustal images, especially in the presence of heterogeneous medium with lateral velocity variations (Das et al., 2019). The calculation of crustal structure using the gravimetric method based on local and regional effect of Earth's crust provides important outputs (Woollard, 1969;Tsuboi, 1979;Chapin, 1996).

Isostasy Analysis
According to the theory of isostasy given by the Airy model, the excess mass or root zone around the base of the crust normally supports the topography. Regional-scale low-amplitude and long-wavelength anomalies are often observed in hilly terrain and near the continental shelves (i.e., edges of continents), and are often masked by anomalies due to upper crustal structures (Simpson et al., 1986). Enhancement of gravity anomalies by removing topography-induced regional effects caused by shallow geologic features is done by routine use of polynomial fitting or wavelength-filtering method. The main objectives of most kinds of geologic and tectonic studies are the density distribution in the upper crust and are more precisely revealed by isostatic residual gravity anomaly maps than any other gravity maps. Qureshy  and Warsi (1980) and Kumar et al. (2020) reconstructed the Airy-Heiskanen isostatic anomaly map over the Indian shield and adjoining regions using the world gravity map (Bonvalot et al., 2012). Several other workers reconstructed the isostatic anomaly maps for the Himalayas and IGB (Qureshy, 1969), Ninety-East Ridge (Tiwari, 2003), Deccan Trap (Qureshy, 1981), and for the Indian continent (Mishra et al., 2004). However, this is the first attempt to reconstruct the isostatic anomaly map for the EIS and adjacent regions using land gravity data (GMSI, 2006). Airy isostatic anomalies were computed from the difference between the gravity effect of topography and a crustal root from the Bouguer corrected anomaly, with the Parker (1972) series up to degree 4 to take care of the nonlinear effects (e.g., Lyons et al., 2000). We fixed the Moho depth to 38 km and formed an average Present study Mandal and Biswas, 2016 Present study Mandal and Biswas, 2016 Present study Mandal and Biswas, 2016 Present study of the spectrum analysis result. The density contrasts between surface and MSL and between mantle and crust are presumed to be 1470 and 400 kg/m 3 , respectively.

GRAVITY PROFILING AND CRUSTAL LAYER MAPPING
The depths to the gravity sources in the sub-surface can be estimated from the shapes of the anomalies or the shapes of the power spectra (Figure 4) computed from the potential field data (Spector and Grant, 1970;Curtis and Jain, 1975;Djomani et al., 1992;Chavez et al., 1999;Lefort and Agarwal, 2002;Rivero et al., 2002;Gomez-Ortiz et al., 2005;Bansal et al., 2006). 2D spectral analysis for 1 • × 1 • and 2 • × 2 • to ensemble basement, Conrad depth, and Moho depth was carried out for archiving reliable information of the subsurface (Tables 2-4).
In the present study, 13 gravity profiles (AA to MM , Figure 3) were selected to investigate the configurations of the crustal layers. Most of the profiles are orthogonal to the regional trend of the gravity anomaly. These profiles run from south to north, extending over ∼225 km, and were selected following the general structural trends of the various shear zones, except profiles 3-4 on the east, which are parallel to the trend of the Bouguer gravity anomaly. Profiles 1-4 (AA to DD , Figure 6), located over the low Bouguer gravity anomaly (−35 to −60 mGal), are associated with the Ranchi-Hazaribagh plateau and UGB, while the profiles 9-13 (II to MM ) lie over the highest Bouguer gravity anomaly (i.e., −20 to +33 mGal), associated with the Rajmahal Trap and surrounding areas. The initial model along profiles AA to II was considered on the basis of spectrum analysis, studies of Verma et al. (1976); Verma (1985), Ismaiel et al. (2019), and Krishna et al. (2019), and other geological background data (Mahadevan, 2002). The final model was selected through minimization of error during optimization of the layer parameters and corresponding densities.
The 2D gravity modeled parameters for the study area (84.75-87.75 • E and 23.0-25 • N) covered by all the 13 profiles (AA to MM ) were used for the computation of parameters of different crustal layers. The results of the 2D gravity modeling were integrated into two maps for the Moho discontinuities and sediment thickness mapping (Figures 7, 8). The mapping of the Conrad discontinuity was not carried out because of its apparent uniformity over the area. The isostatic anomaly ( Figure 9) was computed based on the differences between the Bouguer gravity anomalies and the gravity values of the computed root zones estimated from the topography (Heiskanen and Vening Meinesz, 1958).
The Moho depth map (Figure 7) was generated by the U.S. Geological Survey program called AIRYROOT (Simpson et al., 1983). Moho depth was also computed by 2D spectral analysis, sub-dividing the Bouguer anomaly map into 12 zones of dimension 2 • × 2 • (Regan and Hinze, 1976) with a uniform 75% overlap with the adjacent grid ( Table 2). The sediments depth map (Figure 8) was generated by 2D spectral analysis, sub-dividing the Bouguer anomaly map into 99 zones of 1 • × 1 • dimension (Table 3) with a uniform 75% overlap with adjacent grid (Regan and Hinze, 1976) for better resolution and continuity of data points.

Sub-surface Structures Along the Gravity Profiles
The profiles on the west (AA to DD ) selected over a relatively lower Bouguer gravity anomaly (Figures 6A-D), brought out six layers with some pockets of sediments. The first layer from 0 to 3.5 km comprises the Gondwana sediment with a density of 2.40-2.65 gm/cm 3 and quaternary sediment of density 2.20 gm/cm 3 (Choudhury, 1974). The maximum thickness found in the Damodar valley area corroborates the observation proposed by Verma (1985). The thickness of the second layer consisting of a granitic basement complex with a density of 2.67 gm/cm 3 varies from 0 to ∼10 km, showing thinning toward the north. In some places, the sediments are almost absent, exposing the basement on the surface. The third layer of granite gneiss with a density of 2.71 gm/cm 3 occurs from 2 to 10 km depths, thickening toward the north. The fourth layer of density 2.80 gm/cm 3 is apparently a transition layer between the continental and oceanic crust. Ismaiel et al. (2019) and Krishna et al. (2019) also observed a similar layer at ∼8-13 km depth. The boundary between the lower crust (average density: 2.90 gm/cm 3 ) and the transition layer is found at ∼17-20 km depths. The Moho depth along these profiles is found to vary between ∼38 and 40 km, corroborating the observations of Singh et al. (2015) and Mandal and Biswas (2016). Beneath the Moho, the density of the upper mantle is found to be 3.28 gm/cm 3 .
Profiles 5-8 (EE to HH ) lie over the intermediate Bouguer gravity anomaly (i.e., −80 to −15 mGal), and pass through FIGURE 7 | Moho depth map deduced from spectral analysis technique and AIRYROOT program using a contour interval of 0.5 km (with color-coded index on the east). Blue triangle represents the location of broadband seismic stations in the EIS region and below their blue color texts represent the value of Moho depth calculated by receiver function analysis (Mandal and Biswas, 2016). the Bokaro-Dhanbad coal field and the Purulia shear zone. The only difference is the upper crust exposed on the surface at each profile. Here, sediment thickness is 0-2.5 km and the Moho varies from ∼37 to 41 km. The observations corroborate the observations of Verma (1985), Sharma et al. (2015), and Das et al. (2019). Profile 9 (II ) runs through the Raniganj coalfield. The deposition of the Gondwana sediments (i.e., 2.42 gm/cm 3 ) is laterally continued from 75 to 125 km of the profile. Occasionally, the sediment is punctuated by patches of high density (2.80 gm/cm 3 ) material. Here, Moho is found at a depth of 37-40 km, dipping toward the north. This observation is corroborated by the studies of Verma et al. (1976), Singh et al. (2015), and Mandal and Biswas (2016).
Profiles toward the east (JJ to MM ) pass over the relatively higher Bouguer gravity anomaly (Figures 6J-M). A total of six layers with an additional high velocity crystalline basement of density 2.75 gm/cm 3 are identified along these profiles. Alluvium/sediment forming 0-4 km thick layer, with patches of Gondwana sediments, are found overlain on a 2-6 km thick layer of basement complex. A depression in the Bouguer gravity anomaly, noted at 75 km along the profile JJ , is interpreted to be associated with ∼4-5 km basin of Gondwana sediments (density: 2.42 gm/cm 3 ) in the Raniganj coalfield area. The profile KK shows a merger of ∼0.5-2 km thick Gondwana sediments with the alluvium/sediment of LGB (Verma et al., 1976). We find a high velocity layer at ∼6-8 km depth. A metasediment layer with a density of 2.65 gm/cm 3 is sandwiched between the basement complex and the high velocity layer. A transition layer is identified at ∼10-12 km depth. Along profiles LL and MM , a ∼3-6 km thick layer of alluvium/sediment is identified, which corroborates the observations of Qureshy (1969Qureshy ( , 1970 and Singh et al. (2004). We have identified a fault extending from the surface to 8 km depth, at ∼80 km distance toward the south. In profile LL , a high density (2.80 gm/cm 3 ) material is found intruding the basement complex, but the same intruded materials are sandwiched between two sedimentary layers in the profile MM ; the high-density materials are apparently associated with the Rajmahal traps. In these areas (profiles JJ to MM ), Moho varies from ∼36 to 41 km.

Variation of Crustal Structure
The Moho depth (Table 4 and Figure 7) is found to vary from 40 to 41.5 km in the northern part near the IGB, apparently corroborating the observation of Singh et al. (2015), but differs from Qureshy (1969), who reported the crustal thicknesses beneath the IGB to be 37.5 km. The Damodar Gondwana basin shows moderate to low Moho depth of 37-39 km, whereas the Ranchi plateau shows a maximum depth of ∼39-40 km for the Moho; this result is well corroborated by the study of Mandal and Biswas (2016) (Table 4). The Moho depth varies from ∼37 to 37.5 km beneath the Dhanbad and Bokaro region. This observation differs from Singh et al. (2015), who found the Moho at depths between 36 and 40 km in this region. The LGB in the region of Bengal shows a minimum depth of Moho to be 36.5 km, which agrees with the Moho depth of 35-36 km computed by Singh et al. (2015). Mitra et al. (2008) found the Moho boundary at ∼37.5 km depth near the LGB. Kaila et al. (1992) also found the Moho at ∼35 km depth. We further compared the depth of Moho beneath other continents of a similar age around the world ( Table 4) and found compatible results.
The map (Figure 8) illustrates that the thicknesses of sediments vary between 0.0 and 6.5 km. A strip-like region of ∼160 km lateral extent (∼23.5 to ∼24 • N) and a sediment thickness of 0.5-3.0 km is found to be associated with the Damodar Gondwana basin (Verma and Ghosh, 1977). Sediment thickness varies from 0 to 2 km toward the north near the UGB. Qureshy (1969) and Choudhury (1974) also reported similar thicknesses of the Quaternary sediments in these areas. The eastern part of the region exhibits a maximum thickness of sedimentary cover, reaching more than 6.0 km, particularly in the LGB, and corroborates the observations of Kaila et al. (1992Kaila et al. ( , 1996, Mall et al. (1999), andSingh et al. (2004). At several places, on both sides of the Damodar Gondwana basin, the basement is exposed on the surface, usually found in the Achaean terrain (Verma, 1985;Veevers and Tewari, 1995;Bhattacharya and Bhattacharya, 2015;Acharyya, 2018a).
Isostatic anomaly is found to be varied from −87 to 44 mGal (Figure 9), and the minimum value is found on the northern and eastern part, where the elevation from the MSL is minimum, indicating over-compensation along the UGB. The maximum value of 44 mGal, representing the under-compensation, is identified toward the east and south near the LGB, Rajmahal trap, and Singhbhum shear zone. The coal-rich Damodar Gondwana basin (i.e., Bokaro, Dhanbad, and Asansol areas) shows −1 to −30 mGal isostatic anomalies, indicating over-compensation because of subsiding basement of the areas. The Ranchi and Hazaribagh Plateaus show almost compensated parts with values of isostatic anomalies of 0-5 mGal.

DISCUSSION
The reconstructed map of the study area is characterized by extensively low regional gravity in the range of −40 to −10 mGal (Figure 3), with some highs in the central part in the form of oxbow-shape and the well-defined shear zones characterized by high Bouguer gravity of values around −20 mGal at the southern part. A region in the UGB, elongated from west to east between 24 • and 25 • N in the northern part of the study area, is characterized by relatively low Bouguer value of −50 mGal. A very low anomaly of about −80 mGal has been noted toward the northern part of this region, and might be associated with the downward flexing of the Indian basement beneath the Himalayan Foothills. The Rajmahal trap and LGB show some high Bouguer gravity, ranging between −10 and +33 mGal.
Three individual discontinuities (i.e., basement, Conrad, and Moho) show a varying sediment thickness of zero to 6.2 km with significant lateral density variation from 2.2 gm/cm 3 (approximated to be alluvial sediments) to 2.40 gm/cm 3 (for Gondwana sediments). The upper crustal layer with an average density of about 2.67 gm/cm 3 is followed by the crystalline basement of density 2.75 gm/cm 3 and transitional layer of 2.80 gm/cm 3 . The Conrad layer at ∼17 and 20 km is marked by a sharp increase in density of the lower crustal layer (average density is 2.90 gm/cm 3 ). The model indicates that the crustmantle boundary (i.e., Moho) ranges from ∼36 km in the southeastern part of the LGB to ∼42 km toward north, possibly due to the stretching and thinning of the crust toward the south. The mean density is taken as 3.28 gm/cm 3 for the upper mantle layer ( Table 1).
Elevation of the study area varies from zero to ∼700 m with the highest value observed at its southwest part, and the northern FIGURE 9 | Airy Isostatic anomaly map of the study area, generated by removing the gravity effect of topography and its root (Airy) from the Bouguer gravity anomaly. Dark blue solid circles represent the location of the hot springs (after Sankar et al., 1991). Black stars represent the location of local earthquakes with magnitude M > 2.0, and beach balls show the focal mechanisms of some events. and eastern parts associated with the LGB show minimum elevation (Figure 1). Moho depth varies from ∼36 to 41.5 km with minimum value found in the southeast corner of the study area, near the LGB, but the northern part has a maximum depth of Moho at ∼40.5-41.5 km (Figure 7). Moho dips toward the southeast near LGB, whereas toward the north the Moho is dipping more sharply than other parts of the study area. Sediment thickness varies from 0 to 6 km beneath the study area. Maximum area shows zero sediment thickness, i.e., the basement is exposed on the surface. A few pockets of sediments are found to be associated with coal-bearing Gondwana basin. The eastern part of the study area shows a high accumulation of sediment (Figure 8).

Bouguer and Isostatic Anomalies and the Vertical Stability of the Study Area
Isostatic anomaly is the difference between the Bouguer gravity anomaly and the computed anomaly of root zone computed from the topography. Isostasy anomaly in the study area varies from −87.0 mGal to a maximum of +44.0 mGal (Figure 9). The northern part of the study area (i.e., UGB) is associated with a strong negative Bouguer anomaly (∼−55 to −96 mGal) and moderate to high negative isostatic anomaly (∼−44 to −87 mGal) with minimum elevation of topography. This apparently accounts for the over-compensation, and is presumably caused by a deficit of rock-mass in the basin (Abbas and Subramanian, 1984;Lyon-Caen and Molnar, 1985). The UGB was explained to be a down warping foredeep of the Himalayan orogen, and the sediments in this basin, subducted along with the converging Indian plate against the Eurasian plate, evolved into the Siwalik Hills, which further resulted in the flexing zone into the over-compensated area (Lyon-Caen and Molnar, 1985). Maximum positive isostatic anomaly (H-1 and H-2; +5 to +44 mGal) is observed toward the southeast near the LGB and southern part of the study area near the Singhbhum shear zone and South Purulia shear zone. Here, the low to moderate elevation of topography records little negative (∼−16 mGal) to moderate positive values (∼+33 mGal) Bouguer gravity anomalies. The positive Bouguer gravity anomaly may be due to the presence of high-density material in the shear zones, Raniganj Coal-field, and Rajmahal trap, and the positive isostatic anomaly indicates under-compensation in the region. It is clear from the sediment map that the eastern part is loaded with enormous amounts of sediments (Sengupta, 1966;Abbas and Subramanian, 1984;Alam et al., 2003;Roy and Chatterjee, 2015;Ghosh, 2018). The coal-bearing Damodar Gondwana basin accommodated the E-W extending Damodar River (Verma, 1985;Acharyya, 2018b) where the isostatic anomaly varies from ∼−22 to +5 mGal and Bouguer gravity anomaly varies from ∼−47 mGal to −24 mGal. Moderately negative Bouguer anomaly and moderately negative to slightly positive isostatic gravity anomalies account for either over-isostatic compensation or near to zero-compensation for the region. Within the drainage area, it collects huge amounts of sediments denudated from the adjacent Ranchi and Hazaribagh plateaus (i.e., undercompensated), which makes some areas over-compensated. Ranchi and Hazaribagh plateaus have a maximum elevated topography (∼700 m) and show moderate Bouguer gravity anomaly (∼−38 to −31 mGal) with little positive isostatic anomaly (+4 to +8 mGal). This apparently accounts for the little under-compensation of the area. To inquire into the overall stability, both positive and negative values of isostatic anomalies are apparent over the study area (Figure 9). Negative value within UGB corresponds to the over-compensation, reflecting higher value in Moho depth up to ∼41 to 41.5 km, whereas positive value of isostatic anomaly reflects under-compensation and less Moho depth up to 36 km in the LGB. The western part, near the Ranchi and Hazaribagh plateau, apparently falls under positive isostatic anomaly with a normal Moho depth of 39-39.5 km (Table 4). Overall, the Moho is uplifted in the southeast part, whereas it is subsided toward the north of the study area because of the down-going basement beneath the foreland Gangetic basin.

Earthquake Activities and Seismic Stability of the EIS and Adjoining Regions
The study area is tectonically very significant, comprising mainly the EIS and the Singhbhum craton, starting from the northern end of the Eastern Ghats Belt (EGB) and Bengal estuary to the bending zone beneath the IGB. A geophysical study carried out by Sengupta (1966) indicates the disappearance of the Achaean basement of the shield area below the alluvial plain of the LGB. The exposed Achaean basement in the Shillong plateau in the northeast is apparently a continuation of the Peninsular shield through the Garo-Rajmahal gap (Evans, 1964). A series of faults and fractures pass all through the western margin of the LGB (Figure 2) and are likely continued in the Himalayas toward north (Valdiya, 1976). The buried ridges of this area are likely to be isolated in the LGB in the east. A recent study  shows that the lithospheric deformation in the deeper part of the area is apparently active. Further, numerous hot springs in the Damodar Gondwana basin and adjacent areas near the crest of the flexural bulging of the Indian plate are facilitating reactivation of the crust through the triggering of small magnitude earthquakes because of the strain weakening processes (Mahadevan, 2002;Bilham et al., 2003).
The study area experienced a number of moderate magnitude earthquakes in the past (discussed in section "Introduction") (Figure 2). A recent study  has shown that the faulting processes in the EIS are apparently dominated by strike-slip movements. It was further transformed into extension in the northeast part of the area between the junction of the Ganga and the Brahmaputra River basins. The compression tectonics dominated by thrust faulting presumably triggered many great damaging earthquakes in the Nepal-Bihar segment of the Himalayas (Figure 2). Stretching of the basement toward the northeast of the present study area is also quite active, and likely continued with the subduction of the Indian plate toward the east beneath the Myanmar microplate (Khan, 2005;Singh et al., 2020). The present earthquakes are mostly of shallow origin and low magnitude, occasionally reaching to more than 4.0 magnitude, and few occur in the lower crust . It is found that the depths of 20 events are very shallow (<20 km), associated with the upper crust and likely caused by reactivation processes in the basement. The role of high heterogeneities of the crust (Khan et al., 2016Singh et al., 2019) can also be the reason for the occurrences of low to moderate magnitude earthquakes in this part of Eastern India. The recent 2011 M 6.9 Sikkim and 2015 M 7.8, 7.3 Nepal-Bihar earthquakes, and the great aftershocks of magnitudes more than 6.0, encouraged earthscientists to review the seismotectonics of the Nepal-Bihar-Sikkim segment and its southern shield areas.
Toward the southeastern side of the study area, the Indian oceanic plate is subducting below the Eurasian plate along the Andaman-Nicobar margin  and beneath the Myanmar microplate toward east (Khan, 2005). This has made the Bay of Bengal (BOB) seismically active with incidents of earthquakes of magnitudes (m b ) 5. 5, 5.7, and 5.2 in 19645, 5.7, and 5.2 in , 19895, 5.7, and 5.2 in , and 19925, 5.7, and 5.2 in , respectively. The 1964 Midnapore earthquake was dominated by thrust-faulting (Chandra, 1977), and the other two earthquakes were dominated by strike-slip faultings (Rao et al., 2015). An analysis of earthquakes during 1900-2011 reveals that the older structural features and weak zones near the BOB and adjacent regions were reactivated by regional stress field (Murthy et al., 2010;Rao et al., 2015). The NW-SE trending fracture zones in the BOB are found to be correlated with strike-slip faultings. The seismicity in the southern and eastern parts of the BOB is dominated by thrust faulting due to a predominant N-S compression, which had back-propagated stress from the Himalayan orogen (Biswas and Majumdar, 1997).
Himalayan tectonics are likely linked with the zone of intense deformation of the oceanic lithosphere near the equatorial Indian Ocean (EIO) Basin (Geller et al., 1983;Cloetingh and Wortel, 1986). An east-west extended structural break in this area of the northern Indian Ocean is particularly associated with the diffuse segment , and has isolated the Australian plate and the Indian plate in the south to southwest. The intense seismic activities, along with the occurrences of great earthquakes, operative stress fields, and the high heat flow of the area, also support this observation (Weissel et al., 1980;Geller et al., 1983;Bergman and Solomon, 1984;Weins, 1985). The ∼N-S striking aseismic 85 and 90 • E ridges divide the BOB into three main sub-basins (Gopala Rao et al., 1997), apparently guiding the active stress of the lithosphere, linking the Himalayan operative tectonics through the EIS, Singhbhum craton, and the IGB. Several faults and lineaments on both sides of the Munger-Saharsa ridge are intersecting the Foothills of the Himalaya, and these geo-fractures are likely extended all through the Gangetic basin, concentrating and raising the stress in the flexing zone of the subducting Indian crust (Sibson, 1980;Marshak and Paulsen, 1997;Godin and Harris, 2014;Khan et al., 2017a;Khan et al., 2021). The plate obliquity of 0 • to less than 10 • of the Nepal-Bihar-Sikkim area  further enhances and confines the stress field around the flexing zone through interacting fractures and triggering great earthquakes by the stress build-up processes.

CONCLUSION
The Bouguer gravity anomaly map indicates that the study area is characterized by low (−96 mGal) to moderate (+33 mGal) regional anomalies with the maximum value over the Rajmahal trap, the LGB in the east, and over the well-defined shear zone in the south (Figure 3). The minimum Bouguer gravity anomaly values of −50 to −96 mGal are associated with the UGB, elongated along east-west between 24.5 • and 25.5 • N in the northern part of the study area. Such a low anomaly in the area is possibly associated with the downward flexing of the Indian basement beneath the Himalayan Foothills. Ranchi, Hazaribagh, Giridih, Bokaro, Ramgarh, and the western part of the study area are found to be associated with ∼−30 to −40 mGal anomalies. Further, a few pockets recorded low anomalies of ∼−40 to −50 mGal, as seen in Chaibasa in the extreme south. The 2D gravity modeling ( Figure 5) shows a sediment layer of thicknesses of ∼0-6.5 km, with a significant density variation of 2.2-2.65 gm/cm 3 , underlain by a basement of a density of 2.67 gm/cm 3 . The Conrad is marked by an increase in density from ∼2.67 to 2.90 kg/m 3 at depths between 17 and 20 km. The average density of the upper mantle layer is computed to be ∼3.28 gm/cm 3 . A ∼2-3 km thick high density (3.02 gm/cm 3 ) material overlies the Moho toward the east and is apparently connected with the Rajmahal trap (Singh et al., 2004). The maximum depth of the Moho is found to be ∼40.5-41.5 km toward the north (UGB), while Moho near the LGB has the lowest depth of ∼36 km. Ranchi and Hazaribagh plateau are underlain by Moho at ∼39-39.5 km depths. The uplifted Moho and the maximum deposition of alluvium sediments are found beneath the LGB in the southeastern part (Table 4).
A residual isostatic map provides the status of stability of any study area. Based on isostatic anomaly (Figure 9), it is found that the northern part (Gangetic basin), a down-warping foredeep of the Himalayan orogen, is over-compensated. Eastern and southern parts, comprising the LGB and shear zones (e.g., Singhbhum and North and South Purulia, Figure 1), have predominant records of under-compensation. While the coal-bearing Damodar Gondwana basin is partially overcompensated, the Ranchi and Hazaribagh plateau have achieved partial under-compensation. The present study also shows that the exposed Achaean basement of the shield area at various locations disappears below the Gangetic basin and is affected by active stretching toward the northeast of India . Seismic activities of this region are mostly of shallow origin, particularly associated with the basement within the upper crust. Although the region experienced a few moderate magnitude earthquakes in the historical past (Figure 2), occasional incidents of recent earthquakes have magnitudes less than 5.0. Occasional events are also found in the deeper part and might be related to the reactivation processes in the lower crust. The eastern part of the study area is also activated seismically due to the subduction processes of the Indian plate beneath the Myanmar microplate.
The long-term observations of occurrences of moderate to low magnitude earthquakes in the EIS, Singhbhum craton, and the BOB apparently account for seismic stability of the region, extending from the diffused EIO to the Nepal-Bihar-Sikkim Himalayan region. Lithospheric stress in the BOB is guided by the aseismic 85 and 90 • E ridges, partially accumulated along the different fault-guided zones or other tectonic elements and triggering low-magnitude events. Sometimes, the ∼N-S striking regional fault systems (e.g., EPF, WPF, PF, MSRF, MSRMF, etc., Figure 2) also guide the stress field, concentrating and raising it in the bending zones of the downward migrating Indian lithosphere, are triggering occasional great earthquakes in the Himalaya. The back-propagation stress from the Himalayan Mountain possibly facilitates deformation in the EIO basin (Weissel et al., 1980;Geller et al., 1983;Cloetingh and Wortel, 1986;DeMets et al., 1990). The pure compression in the foothills of the Nepal Himalaya and the minimum obliquity of the converging Indian plate corroborate the observation of stress enhancement in the flexing zone of the down-going plate. We propose that the prevailing stress fields in the converging Indian lithosphere, guided by the aseismic 85 and 90 • E ridges, are migrating through the low-to-moderately stable BOB and EIS. The migrating enhanced stress (because of minimum Indian plate obliquity, Figure 9 of Khan et al., 2014) is further guided by the foothill transecting fault systems to be concentrated in the bending zone of the penetrating lithosphere beneath the Himalaya. This increasing accumulated stress occasionally triggers great earthquakes in the Nepal-Bihar-Sikkim sector of the Himalaya.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
All authors listed have made a substantial, direct, and intellectual contribution to the work, and approved it for publication.

ACKNOWLEDGMENTS
RS is thankful to the Director of the Indian Institute of Technology (ISM), Dhanbad, for financial support to complete the present work and Dhananjay Singh, Geophysicist, M & CSD, Geological Survey of India, Kolkata, West Bengal, India, for his useful and critical suggestions during the computation of 2D gravity modeling. Thanks to both of the reviewers for their valuable comments, which improved the manuscript greatly. The authors are grateful to S. P. Mohanty for thorough revision and correction and improvement of the English of the whole manuscript.