The architecture of the continental lithospheric mantle is vital for understanding the wide-ranging processes that govern the creation of Earth’s continents. The present investigation delineates the crustal and lithospheric boundaries underneath northwest India, which provides essential constraints as an imprint of tectonomagmatic events under the distinct tectonic setting. This study integrates three data sets (gravity, geoid, and topography data) and thermal parameters in a three-dimensional inversion algorithm to northwest India. Inversion results are constrained by available passive seismic crustal thickness, which we are using as a priori information. Our modelling results show crust down to a depth of 33 to 39 km in the Saurashtra and Kachchh region, which increases gradually to about 40 km beneath the Jaisalmer region. The deepest Moho at about 42-45 km is found under the Aravalli-Delhi fold belt and Indo-Gangetic plain. A comparable crustal thickness of 41-44 km is observed beneath the Bundelkhand craton. The mean crustal density varies laterally between 2850 and 2890 kg/m-3 in the study area. The lithosphere-asthenosphere boundary (LAB) is 115-145 km deep under Saurashtra and Kachchh regions, which remains comparable under the Barmer-Pokharan region (130-160 km). A moderate LAB depth under the Aravalli-Delhi fold belt (140-180 km) and relatively deeper under the Bundelkhand craton (175 to 185 km) reaches over 200 km below the Indo-Gangetic plain. The thinned lithosphere beneath the Kachchh, Saurashtra, and Barmer-Pokhran regions, crisscrossed by two different rift systems, together with the dominance of the volcanic centres of two major magmatic episodes, and the possible presence of a low-velocity zone in the subcrustal upper mantle region suggest a genetic link between the evolution and the multistage thermal modification of the mantle lithosphere of this region.

The lithosphere is a chemical, thermal, and mechanical boundary layer that protects the continental crust from the Earth’s hotter and more dynamic asthenosphere [1, 2]. It developed in the Archaean by methods similar to contemporary tectonics and has remained stable beneath the larger cratons ever since [3]. When this cold, thick (~200 km), chemically depleted, and a buoyant protecting layer is disrupted or removed by thermomechanical processes, the heat budget of the crust alters, resulting in uplift, magmatism, and creation of large mineral deposits [4]. Understanding the mechanisms of continental formation and evaluating the long-term stability of continents require an understanding of the structure, composition, and thermal evolution of the Archaean continental crust and underlying mantle lithosphere [5]. The Archaean northwestern Indian shield is one of the geologically and tectonically complex regions that has undergone multistage thermal perturbations and hence requires in-depth investigations (Figure 1).

Growth of the northwestern Indian shield, south of the Indo-Gangatic belt, is evidenced by the crustal generation and accretion of Paleoarchaean Bundelkhand craton (BC) and Late Neoproterozoic Marwar block (MB) sandwiching the Palaeo- to Neoproterozoic Aravalli-Delhi fold belt (ADFB) [6]. These three tectonic domains together make the Northern Indian block welded with the Southern Indian block consisting of Archaean Dharwar, Bastar, and Singhbhum cratons along the Central Indian Fold Belt (CIFB) [7, 8]. Northwestern India including the Kachchh and Saurashtra peninsula is thus one of the segments of the old Earth’s crust showing collision and suturing of distinct terranes in the geological past, which has significantly contributed to the subcontinent’s evolution. Post stabilization of the lithosphere, the northwestern edge of the Indian shield bounded between 20-30°N latitude and 68-78°E longitude was under the influence of two large igneous provinces, namely, Neoproterozoic Malani Igneous Suite (MIS) and K-T boundary Deccan volcanism (Figure 1). Indeed, all these tectonic and magmatic events together with reactivated rifted basin systems [9] must have had a significant impact on the lithospheric growth and its modification. Delineation of the present-day lithospheric structure of northwestern India is thus of fundamental importance for understanding its evolution, growth, and possible modification.

Several investigations using different geophysical proxies have contributed considerably towards the understanding of the deep crustal structure of northwest India [1016]. However, the geophysical investigations for depth of the lithosphere-asthenosphere boundary (LAB) remain scanty [1722]. Most of the studies focused on a particular area using a single geophysical proxy that is often debated because of its inherent limitations [1, 2], whereas all these proxies measured using direct or indirect geophysical investigation are interrelated and depend on pressure, temperature, composition, and the physical state of matter [2328].

To fully comprehend the distinctive density distribution, it may be imperative to critically examine the various geopotential proxies that have tell-tale imprints of crust-mantle interactions. The availability of dense satellite gravity anomaly and geoid undulation data, along with the high-resolution digital topography model ETOPO1 [29], opens new possibilities for integrated lithospheric density modelling. These data enable us to evaluate the significant involvement of the crust and subcrustal lithosphere in preserving surface topography and, as a result, simulate mantle density structure. This study uses the 3D inversion technique based on Bayesian approaches to present density models of the deep crust-mantle structure and lithospheric thickness variation. Program inverts the free-air gravity anomaly, geoid, and topography data to produce the image of Moho and LAB depths (with reference to mean sea level) and average crustal density of the region. Parameter damping and smoothing as well as the use of a priori information, which includes average crustal density and Moho depths from seismological investigations, stabilize the inversion process. The derived upper mantle structure would eventually address the different geological and geophysical features developed in northwest India, and how the deep and surface processes are interconnected.

Northwestern India is a collage of five main lithotectonic domains, namely, the Bundelkhand craton (BC), the Marwar block, the Deccan volcanic province (DVP), the Aravalli-Delhi fold belt (ADFB), and the Central Indian Fold Belt encompassing Son-Narmada Fault (SNF). These Precambrian lithological and tectonic domains are covered by Neoproterozoic acid volcanic rocks of MIS and the K/T boundary flood basalt of DVP (Figure 1). The Paleoarchaean Bundelkhand craton consists of three distinct lithological units viz. (1) an extensively deformed older gneisses-greenstone of Archaean (3.59-2.6 Ga), (2) the undeformed multiphase granitoid plutons of the Early Proterozoic (2.5-2.1), and (3) clusters of younger mafic dykes [30]. Bijawar, Sonrai, and Gwalior are three isolated and narrow basins that exhibit Neoarchaean to Paleoproterozoic volcano-sedimentary units overlying the basement rocks of Bundelkhand craton [31]. The Palaeoproterozoic Vindhyan basin overlying the Bundelkhand craton is bounded by the Great Boundary Fault in the west, the Phanerozoic Indo-Gangetic plain to the northeast, and the Late Cretaceous DVP to the southwest. A few ENE-WSW and NE-SW trending younger dykes also exist, crossing all Archaean lithologies [32].

The NE-SW trending ADFB preserves a long geological history spanning from the Palaeoproterozoic Aravalli fold belt in the east to the Meso- to Neoproterozoic Delhi fold belt in the west [33]. These intensely folded, deformed, and metamorphosed Proterozoic rocks were deposited over the Archaean Banded Gneissic Complex or Bhilwara Gneissic Complex (BGC) of 3.3 Ga age [34]. The 1800 Ma Aravalli orogeny, the 1100 Ma Delhi orogeny, and the 850-750 Ma post-Delhi magmatic events are the four major tectonomagmatic and metamorphic events in this region [35].

The Marwar block (also known as the Trans-Aravalli Terrane) is a geographical region west of the ADFB having a Neoproterozoic to recent geological history. Based on U-Th-Pb (total) monazite age determination data, Chatterjee et al. [36] argued that the Phulad shear zone located along the western margin of ADFB is the most plausible location along which the Greater India landmass accreted with the Marwar block at about ~810 Ma. The Precambrian-Cambrian Marwar supergroup and Phanerozoic sedimentary rocks are largely exposed in the southern section, while a substantial portion of the Marwar block is buried beneath the Quaternary sand deposits of Thar Deserts.

Down south the Kachchh and Saurashtra regions forms an important site of Mesozoic and Cenozoic sedimentation widely covered by Late Cretaceous Deccan basaltic flows overlain by a thin cover of Neogene as well as Quaternary sediments. The Kachchh region is a distinctive terrain due to the existence of active faults and associated seismicity, as well as its evolutionary history. The important structural features that have played an important role in Kachchh evolution include a group of E-W trending uplifts surrounded by residual depressions [9], whereas the Saurashtra peninsula is dotted by several volcanic plugs of alkaline complexes. The Proterozoic crystalline basement is the well-known geological features in these two areas [37]. According to Wang et al. [38], the Marwar block encompasses the Kachchh and Saurashtra regions since the Gondwana reconstruction.

The Indo-Gangetic Plain towards the north divides the Himalayan ranges from Peninsular India. This large alluvial plain, formed by the Ganga, Indus, and Brahmaputra rivers, is integrally linked to the development of the Himalayan orogeny. The Himalayan uplift resulted from the convergence of Indian and Eurasian continental plates, producing the neo-tectonically active Indo-Gangetic Plain [39, 40]. Due to the thick mask of alluvial deposits with the Siwalik sediments, the subsurface geology of the Indo-Gangetic plain is mostly inaccessible, and speculations were made concerning the continuity of important subsurface structural and tectonic features of the Indian peninsular shield, such as the ADFB and Vindhyan group of rocks [41].

The MIS is one of the most extensive anorogenic acid volcanic basalts of India that spans over the Marwar block, western Rajasthan from south of Sirohi to the north of Pokaran and from east of Jodhpur to the border of Pakistan. To de Wall et al. [42], MIS is the northeastern continuation of the Neoproterozoic (800-700 Ma) magmatic belt extending from northern Madagascar, Seychelles into north-west India [43, 44]. The MIS magma erupted around 770-750 million years ago [45, 46] and was possibly triggered by a mantle plume [47]. The Deccan volcanic province (DVP) is the youngest and largest igneous province in India that is often blamed for the Cretaceous-Tertiary K/T boundary mass extinction. Multiple strata of tholeiitic flood basalt making up the bulk of the trap complex were produced during India’s northward migration when it passed over the Réunion hotspot. According to detailed radiometric dating of the Deccan lava, roughly 80% of total basaltic lava erupted rapidly within a short period (less than 1 Ma; [48]). The mafic alkalic complexes of Mer-Mundwara, Sarnu-Dandali, and Tavidar volcanic rocks of a Late Cretaceous to Early Tertiary age show the transition from MIS to Deccan volcanism in southern Rajasthan ([49] and references therein). These two major thermotectonic processes must have had complex mantle-crust dynamics in northwestern India.

The structure of the Earth and its density distribution at various levels are intimately manifested in gravity, geoid, and topography. We applied fast generalized linear iterative 3D inversion techniques of Motavalli-Anbaran et al. [24] to determine the lateral changes in crustal thickness, mean crustal density, and LAB. Our working model can be expressed by the structure and the density parameters, which divided the study area into fixed-size rectangular columns in the N-S (Y) and E-W (X) directions. The depth (Z) of each column is divided into four different sections: seawater (if it exists, with defined thickness), crust, lithospheric mantle, and the asthenosphere.

At the outset, the lithospheric temperature is determined using a 1D technique to facilitate the inversion technique to go faster. The constant temperature at the surface (typically 10°C) and the LAB (usually 1300°C) are used to determine the geotherms in-between. The continental crust is divided into two parts of uniform width equivalent to the upper and lower crust. The upper half of the continental crust has a standard heat production and thermal conductivity of 1 W/m3 and 2.5 W/(Km), whereas the lower half and oceanic crust have 0.2 W/m3 and 2.2 W/(Km), respectively; the mantle has no heat production and a thermal conductivity of 3.3 W/(Km) [50, 51].

The density of seawater (1030 kg/m3) and that of the asthenosphere (3200 kg/m3) is kept uniform and constant. The crust is treated as a single layer in which density increases continuously with depth from the Earth’s surface (or seafloor) to the Moho. Moreover, there is no distinction between sediments, upper crust, and lower crust. To minimize the consequences of this shortcoming, the program enables users to model vertical density gradients in two distinct ways: first is to fix the density at the bottom of the crust, and the program reverts with density differences at the surface, resulting in varying vertical density gradients. Another possibility would be to set the density difference between the surface and the Moho, after which the program reverts to the mean crustal density. Here, we set the density at the Moho to 3000 kg/m3, and the surface density is inverted. The density of the lithospheric mantle varies linearly with temperature [52]:
(1)ρT=ρT01αTT0,
where T0 is the reference temperature, which is 1300°C at the lithosphere-asthenosphere boundary, ρT0 is the density at this temperature, and α is the coefficient of thermal expansion.
With the derived density distribution in the lithosphere, we can calculate the gravity, geoid, and topography response of the upper mantle using a forward modelling approach. The program calculates elevation in one dimension and assumes local isostasy as a function of thicknesses [53]:
(2)ε=fρaρlρaZl+εcal,
where Zl is the depth of the LAB, ρl is the mean lithospheric density, ρa is the asthenospheric density (considered as 3200 kg/m3), and εcal is the calibration constant which coincides with the depth of a free asthenosphere (-2380 m; [54]). For positive topography, the factor f equals to one and for negative topography equals to ρa/ρaρw assuming that the topographic depression is occupied by seawater having density ρw=1030kg/m3. Gallardo-Delgado et al. [55] and Fullea et al. [23] equations are used to compute the gravity and geoid effects, respectively.
For the linearized problem, the forward modelling describes the relationship between the data and model parameter as follows:
(3)D=fp,
wherein D is the data vector, which includes gravity, geoid, and topography data, p is the vector of the modelling variable, which comprises average crustal density, and crust and lithospheric thicknesses, and f is a function of gravity, and geoid and topography. The inversion equation of the potential field data and elevation becomes ill-posed because of the solutions’ nonuniqueness and instability. Problem related to the instability of the solution, particularly different uncertainties for different data sets, is resolved through model optimization. During inversion, the following data misfit function is used in the least-squares sense with the weighted norm:
(4)Ed=dfPTCd1dfP,
wherein Cd is a diagonal matrix that contains the variances of each data point while normalising the variability of the various data types. To overcome the nonuniqueness (ill-posedness) problem of inversion, a linearization term is incorporated for the model parameters making the inverted matrix nonsingular. Even though the geological nonuniqueness is not rectified, this approach brings a mathematically unique solution.
(5)Ep=PP0TCp1PP0.
Here, Cp is a parameter covariance matrix taken as a diagonal matrix having parameter variances σ2pi on the diagonal. P0 is a vector comprising the starting model’s parameters. Another significant application of such standardization parameters is using them as a definition of a priori knowledge. Suppose a priori information is available at certain places, like depth of crust-mantle boundary from seismic investigations or subsurface densities through well-log observation. Such parameters could be utilised to build the initial model and confine the solutions at such places. It is also essential to somewhat smooth the models to prevent geologically inappropriate models. The smoothing of a model is defined as the following:
(6)Es=PTCsP,
where Cs is the matrix describing the model parameters’ finite-difference first derivative [24].
Throughout this inversion operation, a cost function must be optimized and defined as follows after integrating the data misfit, model modification, and model smoothing:
(7)C=Ed+λEp+μEs.

Here, λ and μ are variables that influence the weighting of parameter constraints as well as the smoothing procedure during the entire inversion with regard to the data fit.

The preceding matrix formula is calculated using an iterative solution procedure to reduce the above problem to its simplest form [24, 56]:
(8)Pk+1=ATCd1A+λCp1+μCs1ATCd1dfPk+Pk.

The Frechet matrix A contains the data derivatives regarding parameters to each data point, and k is the iteration count.

To provide an accurate depiction of both crustal as well as lithospheric architecture, we extracted data from publicly available worldwide sources. Free-air gravity (Figure 2(a)) and topography (Figure 2(b)) were retrieved using the TOPEX global data sets ([57]; ftp://topex.ucsd.edu/pub), respectively. The free-air gravity anomaly over the study regions shows variation ranging from less than -140 mGal towards the northwestern part to more than +110 mGal over the Pokhran, Barmer, and ADFB. The average topography of about 200-800 m is observed in the continental part of India. Minimum and maximum topography was observed in the southwestern (offshore) and northeastern parts of the study region, respectively. The topography in the west is typical of low-lying areas, and this pattern continues throughout the study area, along the ADFB at the map’s middle region, where the maximum height is almost 1000 meters above sea level. Topography eventually leads to the Indo-Gangetic plains in the east, which are presumably controlled by topographic ridges and valleys underlying the alluvium. The geoid anomaly is extracted from the EGM2008 global model [58], which is complete to spherical harmonic degree and order 2159 and can provide 10 arc minutes (~18.5 km) resolutions to our acquired model. This geoid anomaly is dominated by the long-wavelength (sublithospheric) signatures, and to remove these signatures from the observed geoid anomaly, the EGM2008 geoid anomalies were filtered equivalent to the lower spherical harmonics up to degree and order 12 (Figure 2(c)). The obtained geoid anomaly shows maximum (5 m) and minimum (-11 m) values situated inside the northern part of the study region in Pakistan. In India, the minimum geoid undulation (-9 m) is situated in the north and east of ADFB and the maximum residual geoid undulation (2 m) is at Rajkot. Barmer and southern part of the ADFB are also characterized by high geoid undulations of 2 m and 3 m, respectively.

The three data sets (topography, free-air gravity, and residual geoid) projected onto a Cartesian coordinate grid of 10×10 km, yielding 112×123 grid points. Each block’s data has been merged into vertical columns of 30×30 km at 10×10 km intervals. The model is comprised of columns of 30×30 km size, 37 blocks in east-west and 41 in the north-south direction, providing 1517 columns and 4551 unknowns (for every single column mean crustal density, Moho depth, and LAB depth). Using only topographical data and geoid anomaly, we built a starting model (Figures 3(a)–3(c)) using a 1D technique described by Fullea et al. [59]. We selected a fixed crustal density of 2860 kg/m3 for this 1D model, which is near to the average crustal density expected for a 40 km thick crust [60]. Crustal densities were permitted to vary with depth in the subsequent 3D inversion, and we chose to alter density gradually from 3000 kg/m3 at the crust-mantle boundary, decreasing upward. To remain stable and also to minimize the nonuniqueness in the inversion process, we used two distinct approaches. Wherever passive seismic data on crustal thickness exists, we used those values like a priori information to minimize the uncertainty of Moho depth in the iteration process for the respective columns. A total of 64 Moho points, mostly from earlier passive seismic studies, are used [22, 25, 6166]. The average of the Moho depths is considered when different Moho depths were available for a single column. Besides, we applied the slight lateral smoothing factors of the model parameters (Table 1), such as the data uncertainties as well as the variance of the model parameters. Table 2 lists the thermal properties considered in temperature calculation; the asthenosphere temperature is set to 1300°C, while the surface temperature is set to 10°C. We have run a series of inversions with changing values for critical parameters like data and model parameter variances, smoothing, and a priori Moho depths within their uncertainties. Finally, preferred a model based on data misfit or that best fit the observed gravity and geoid anomaly and topography data from a set of more than 88 models created this way. For the chosen model, the disparities between observed and estimated data for free-air gravity, geoid, and topography are shown in Figures 2(d)–2(f), respectively. The gravity data has a mismatch standard deviation of 7.61 mGal, the geoid data has a misfit standard deviation of 0.26 m, and the topography data has a misfit standard deviation of 67.32 m. A fairly good agreement is achieved except for short-wavelength mismatches, which can be related to differences in density and shallow crustal features at scales lower than 30 km. The applied local isostasy notion for such a short-wavelength condition would contribute to the misfit.

The parameter uncertainty (Figures 4(a)–4(c)) associated with the propagation of data inaccuracies into parameter values is usually relatively low. Since the uncertainty resulting from altering inversion parameters may be substantially more significant than the uncertainty resulting from data mismatch, the typical posterior variance matrix [56] is insufficient to evaluate model parameter uncertainty. As a result, we tested 88 inversions with varying starting models by changing various parameters for the 1D inversion, different constant densities in the lower crust, Moho a priori data in their uncertainty ranges, smoothing factors, and damping parameters. Finally, we computed standard deviations from each parameter’s results, showing that the parameter uncertainties (Figures 4(a)–4(c)) in the continental region are in the order of 0.5 to 1.0 km for the Moho depth, a 10-15 km for the lithosphere-asthenosphere boundary, and 30-45 kg/m3 for mean crustal densities. However, discrepancies in Moho thickness variations up to 2 to 3 km and LAB depth uncertainty ranging between 15 and 20 km have been observed in the oceanic area. Similarly, the average crustal density has maximum errors of 55-80 kg/m3 in oceanic zones. It appears that our data in oceanic areas are not adequately constrained by the a priori crustal parameters. As shown in Figures 2(d)–2(f), the model with the minimum data misfit is considered the final model. The program also allows us to assess the quality of the inversion result. The resolution matrix and the a posteriori covariance matrix have been used in this approach. The resolution matrix demonstrates how well each parameter can be resolved using only existing data. The parameter resolution is shown in Figure 5 in the form of maps. The variances between Moho prior knowledge and the final Moho depth in the relevant blocks are an average of ±2.6 km, which is acceptable given the typical seismic Moho depth uncertainty of 2-3 km. Figures 3(d)–3(f) illustrates our resulting Moho depth, LAB depth, and average crustal density after 19 iterations, with data misfits of 7.61 mGal, 0.26 m, and 67.32 m for gravity, geoid, and topography, respectively (Table 3). Figure 6 illustrates the evolution of these misfits during the iteration phase.

Our best-fitting model-derived Moho depth ranges from ~30 km in the adjoining Arabian Sea to ~50 km in the Indo-Gangetic Plain (Figure 3(d)). The increase in crustal thickness from the coastal region to the continental region is consistent with the previous research in the area [20, 22, 25, 6165]. The thinnest crust is found beneath the oceanic region, with thicknesses ranging from 30 to 33 km. Our estimate of crustal thickness of 33 km along the coastline of Saurashtra is close to the average mean sea level crustal thickness computed for India [67]. In the Saurashtra and Kachchh regions, the crustal thickness ranges from 33 to 38 km, with the northern section of the region being thicker. Comparable crustal depths of 33-37 km in the Saurashtra region [68] and 31-33 km in the Cambay basin were reported by deep seismic data [69] and 3D gravity modelling [70]. The shear velocity structure obtained by inverting regionalized group velocity dispersion data in the Saurashtra and the south of the ADFB regions revealed that the crustal thickness is about 37-39 km with an almost flat Moho [21]. In the Marwar block, the derived crustal thickness gradually increases from 38 km to 43 km in SW to NE direction and agrees with the seismological observations [71, 72]. The crustal thickness image (Figure 3(d)) illustrates Moho variation ranging from 38 to 45 km in the DVP region (east of Surat), which corroborates well with passive seismic studies [22, 61, 73] and geopotential investigations [7476]. Our inversion result suggests the ADFB be a thick crust zone that varies between 35 and 44 km and conforms to seismic refraction/wide-angle reflection investigation derived crustal thickness of 45-47 km [77]. Using the common reflection surface stacking approach, the reconstructed Nagaur-Jhalawar DSS profile revealed crustal thickness variations of 38 to 50 km beneath the ADFB [15]. Inverted depth analysis based on receiver function studies at locations near the ADFB’s northern boundary, i.e., 40 to 44 km [62, 78, 79]. Our study indicates that the Bundelkhand craton is characterized by a uniform crust of 41-44 km thickness. Joint inversion of teleseismic receiver function and surface wave group velocity estimated Moho at a depth of 38 km beneath the Bundelkhand craton [80]. In the Indo-Gangetic plain, the derived average crustal depth ranges from 42 to 45 km, similar to the crustal depth estimated by Mitra et al. [81] and Caldwell et al. [82] using passive seismic investigations.

The derived LAB changes significantly throughout the study region (Figure 3(e)). The thinnest portion of the lithosphere is ~100 km beneath the oceanic region, which gradually increases to around 110-125 km along Saurashtra’s continental margin. In the mainland of Saurashtra, the LAB depth estimates range from 115 to 140 km whereas using regional S-wave tomography, the LAB is identified at an average depth of about 122 km in the same region [21]. As the surface heat flow observed in the region is high, Negi et al. [17] estimated the thin thermal lithosphere Saurashtra region. Our LAB depth in the Kachchh region remains the same at the depth of about 135-140 km while the LAB delineated by Sharma et al. [21] is normally at 118 km which further shallows to 82 km in the Kachchh seismic zone. Mandal and Pandey [83] suggested a weak lithosphere beneath the Kachchh seismic zone due to the interaction of the Deccan/Reunion mantle plume with the Indian lithosphere. The Cambay basin is again underlain by a comparative continental lithosphere ranging from 135 to 140 km. Negi et al. [84] also reported a shallow Curie depth and thermal lithosphere in the Cambay and Kachchh regions. The LAB depth from our study is therefore consistent with prior seismological estimations (e.g., [20, 21, 71, 85, 86]). In the Marwar block, the Barmer and Pokhran gravity highs (Figure 2(a)) are characterized by shallower LAB at 130 and 140 km, respectively, else, it varies between 150 and 170 km depth towards the north. A comparative LAB depth (140-180 km) is also observed beneath the ADFB. The heat flow in the ADFB varies between 37 and 74 mW/m2, and the computed temperature-depth distribution places the thermal LAB at roughly 95 km [84]. A shallower LAB (135-150 km) characterizes the Narmada-Tapti region. The eastern region of DVP is underlain by a moderate lithosphere ranging from 150 to 170 km, thickening gradually under the northern regions attaining a depth of over 180 km beneath the Bundelkhand craton. Negi et al. [84] analysed the relationship between heat flow and lithospheric thickness for the Indian subcontinent, which reveals the LAB at a depth of about 101 km in the Deccan trap region. Further north, we observe about ~180 to 200 km thick LAB in the Indo-Gangetic plain, mostly occupied with thick alluvial sediments, which indicate the Indian lithosphere under thrusting beneath Eurasia [20].

The analysis of gravity-based geological models needs a comprehensive understanding of the density distribution in the continental crust for a better experience of geodynamic processes. In the present study region, inversion of gravity anomalies produces a nonuniform average crustal density distribution (Figure 3(f)). The average crustal density of 2890 kg/m-3 in the northern part of the study area is higher than that in the southern part of 2850 kg/m-3. However, the southern region shows some significant pockets of high density, which is coherent with the regional geology observed at the surface. Several high-density areas can be seen over Saurashtra, which is related to known Deccan volcanic plugs. Similarly, the high density over Kachchh indicates the existence of intrusive, similar to Saurashtra’s volcanic plugs [87]. The intrusive mass sitting at the midcrustal level or magmatic underplating at the base of the crust are plausible causes of the observed high density in the Narmada-Tapti area, with a maximum value of ~2890 kg/m3 above the Navsari [7476]. Another significant density highly correlates well with the ADFB, with surrounding low density on either side associated with the Erinpura granite and the Marwar supergroup of rocks in the northwest and the Berach granite and the Hindoli group of rocks in the southeast [12]. The MIS is defined by a north-south high density of ~2910 kg/m3 along 72°N longitude, concentrated on the Pokharan-Phalsund-Saranu region, and about ~2890 kg/m3 to the west of Barmer [16]. Another relatively small high-density anomaly is observed at the eastern flank of the Siwana Ring Complex. The Proterozoic Vindhyan sediments or the presence of granitic batholiths over the Bundelkhand craton can be attributed to low-density felsic rocks, which are intruded by high-density mafic and ultramafic Bijawar rocks. Some of the widespread density highs in the Indo-Gangetic plain showed no clear trend or association with the exposed geology and were initially believed to be originated by high-density rocks emplaced in the crust as a result of crust-mantle interaction.

The broad spectrum of measured Poisson’s ratio (0.22-0.36) in the Kachchh region is however attributed to crustal heterogeneity with the possible presence of partial melt and/or the crust rich in mafic composition [64]. The intracrustal low velocities encompassing the high-density eruptive centres of MIS are interpreted to be the compositional imprints of the felsic magma associated with the bimodal Malani volcanism [72]. The Poisson’s ratio of 0.26 to 0.29 under the MIS in the Marwar block signifies thermal and/or compositional perturbations in the lithosphere for the presence of mafic material at the crust base [72]. Even higher values of Poisson’s ratio (≈0.30) may be related to the existence of a high-velocity/density lower crust with a mafic or ultramafic composition [16, 88].

Lithospheric roots from the Archaean and Proterozoic periods are cold, thick, and buoyant that cannot be delaminated; instead, they may be thermomechanically disaggregated (lithospheric thinning and/or rifting) or infiltrated with upwelling fertile mantle material to be modified [89]. These growth or modification processes lead to density differentiation vis-à-vis geopotential anomalies that are extremely useful for interpreting the buoyancy of the lithosphere supporting the present-day topographic load. Generally, the thicknesses of the continental lithosphere are ~200 km in the Southern Hemisphere and ~300 km in the Northern Hemisphere, respectively [90]. Our integrated geopotential model shows a normal lithosphere in the eastern part of the DVP (155-165 km), Bundelkhand craton (175-185 km), the northern part of the Marwar block (160-175 km), and the Indo-Gangetic plain (>200 km). In contrast, an atypical lithosphere is observed underneath the Saurashtra and Kachchh regions (115-140 km), Bay of Cambay and Cambay Basin (130-140 km), Narmada-Tapti Region (135-150 km), the southern part of the Marwar block (130-160 km), and ADFB (140-180 km) (Figure 3(f)). Our derived lithospheric density model in these regions is thus thinner than normal Archaean cratons. A 5-10% reduction in S-wave velocity for the uppermost mantle in the DVP, Kachchh area, and ADFB is noticed as compared to the standard IASP91 model and concluded that the younger processes might have had reworked the Archaean and Proterozoic lithospheric keel of the northwest Indian shield [91]. The Archaean subcontinental lithospheric mantle was replaced through rifting, with concurrent upwelling of fertile asthenospheric material along with breaches in the Archaean root, resulting in thinning and spread of the Archaean subcontinental lithospheric mantle [92]. As a result, it is critical to understand the likely thermomagmatic mechanisms responsible for the possible lithospheric modification of the western Indian shield.

The Neoproterozoic (771-751 Ma) MIS event has three episodes of mantle-derived mafic magma, granites, and rhyolites along with some dolerite dykes [93]. The MIS granitoids of Mt. Abu in southwestern ADFB are 509-514 Ma old and appear to be related to Pan-African events [94], after a long hiatus emplacement of alkaline basalt magma at Mundwara Igneous Complex and Sarnu-Dandali Alkaline Complex in southern Rajasthan took place between ca. 69 and ca. 67 Ma. Further down south, a massive eruption of tholeiite Deccan basalt took place at ca. 66-65 Ma covering a wide region of west-central India including the Kachchh and Saurashtra. This widespread magmatism over a prolonged period of thermal overprinting must have had modified the lithosphere of the entire region considerably.

The presence of alkaline complexes between MIS and DVP regions and their close relationship with Deccan magmatism led Kennett and Widiyantoro [18] to speculate that the cylindrical zone of low P-wave velocity observed in the upper mantle (80-150 km depth) might be linked to the fossil mantle source associated with the 65 Ma Deccan volcanic episode which erupted further south of the MIS. The event may likely have overprinted or obliterated the mantle signatures of the Neoproterozoic (750 Ma) tectonomagmatic event. This low-velocity zone may also reflect the influence of a thermal anomaly beneath the rift systems of northwestern India.

Continental growth of the northwestern Indian shield is examined by delineating the present-day structure of the crust and lithosphere using integrated inversion of geopotential (geoid and gravity) and topography data. Assuming the state of local isostatic equilibrium, the 3D inversion results are constrained by existing Moho depth data from passive seismic studies as a priori information. The findings of the integrated 3D inversion modelling lead to the following conclusions:

  • (i)

    Crustal thickness from 33 to 39 km underneath the Saurashtra and Kachchh region increases gradually to about 40 km beneath the Jaisalmer region. The deepest Moho at about 42-45 km is found under the ADFB and Indo-Gangetic plain. A uniform crustal thickness of 41-44 km is observed beneath the Bundelkhand craton

  • (ii)

    The average crustal density (2890 kg/m-3) in the northern part of the study area is higher than that in the southern part (2850 kg/m-3). However, several pockets of high-density colocated with the tectonically modified zones are observed in this region

  • (iii)

    Significant variations in the lithospheric mantle thicknesses are discerned in the resulting lithospheric structure. Our LAB depth is ~115-140 km under Saurashtra and Kachchh regions, which remains the same under the Barmer basin. A moderate LAB depth (140-180 km) under the Aravalli-Delhi fold belt increases to 175 to 185 km beneath the Bundelkhand craton and reaches ~200 km below the Indo-Gangetic plain

  • (iv)

    Our findings go by geodynamic notions of a genetic connection between growth and thermomechanical modification of the mantle lithosphere of this region

This article’s worldwide data sets may be found at Topography: ftp://topex.ucsd.edu/pub/srtm15_plus/ [97]; Free-air: ftp://topex.ucsd.edu/pub [57]; and Geoid: http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm2008/index.html [58].

The authors declare that they have no conflicts of interest.

The authors express their gratitude to the Director, CSIR-National Geophysical Research Institute, for his unwavering support and permission to publish this paper. The Generic Mapping Tools software is used to create three figures [98]. CSIR-NGRI contribution NGRI/Lib/2021/Pub-157 is under the INDEX program (PSC-0204).

Exclusive Licensee GeoScienceWorld. Distributed under a Creative Commons Attribution License (CC BY 4.0).