Deep low velocity layer in the sublithospheric mantle beneath India

SUMMARY

Globally, there is now a growing evidence for a low velocity layer in the deeper parts of the upper mantle, above the 410 km discontinuity (hereafter called LVL-410). The origin of this layer is primarily attributed to interaction of slabs or plumes with a hydrous mantle transition zone (MTZ) that results in dehydration melting induced by water transport upward out of the MTZ. However, the ubiquitous nature of this layer and its causative remain contentious. In this study, we use high quality receiver functions (RFs) sampling diverse tectonic units of the Indian subcontinent to identify Ps conversions from the LVL-410. Bootstrap and differential slowness stacking of RFs migrated to depth using a 3-D velocity model reveal unequivocal presence of a deep low velocity layer at depths varying from 290 to 400 km. This layer appears more pervasive and deeper beneath the Himalaya, where detached subducted slabs in the MTZ have been previously reported. Interestingly, the layer is shallower in plume affected regions like the Deccan Volcanic Province and Southern Granulite Terrane. Even though a common explanation does not appear currently feasible, our observations reaffirm deep low velocity layers in the bottom part of the upper mantle and add to the list of regions that show strong presence of such layers above the 410 km discontinuity.

1 INTRODUCTION

Global scale analysis of seismic phases sampling the mantle transition zone (MTZ, Fig. 1) reveals presence of a low velocity layer (LVL) at staggered depths above the 410 km (d410) discontinuity, without a tectonic affinity (Tauzin et al. 2010; Wei & Shearer 2017; Han et al. 2021). Receiver functions (RFs, Vinnik & Farra 2002, 2007; Vinnik et al. 2003, 2010; Fee & Dueker 2004; Jasbinsek & Dueker 2007; Wittlinger & Farra 2007; Tauzin et al. 2010; Schaeffer & Bostock 2010; Jasbinsek et al. 2010; Schmandt et al. 2011; Tauzin et al. 2013; Huckfeldt et al. 2013; Pérez-Campos & Clayton 2014; Morais et al. 2015; Thompson et al. 2015; Bonatto et al. 2015; Liu et al. 2018; Taylor et al. 2018; Zhang & Dueker 2020; Frazer & Park 2021), P/S wave triplications (Song et al. 2004; Gao et al. 2006; Obayashi et al. 2006; Li et al. 2017, 2020; Guohui et al. 2019; Li et al. 2022a), SS precursors (Wei & Shearer 2017) and ScS reverberations (Courtier & Revenaugh 2007; Bagley et al. 2009) are the methods that have been commonly employed to detect this layer. Fig. 2 shows the geographical locations where the LVL-410 is detected in diverse regions of the globe and Table 1 provides a synthesis of the nature of the finding, causative mechanisms and the technique used to detect it.

A popular explanation for these LVLs is partial melting of ascending mantle material caused by slabs subducting into a hydrous MTZ—called the transition zone water filter model (Bercovici & Karato 2003). Experimental results simulating an upwelling peridotite mantle across the d410 discontinuity reveal that a negative Vs anomaly of 4 per cent can be explained by a 0.7 per cent melt fraction in peridotite at the base of the upper mantle, corresponding to a total water content of 0.22 |$\pm$| 0.02 wt per cent in the transition zone (Freitas & Manthilake 2019). On the other hand, a negative discontinuity at |$\sim$|350 km depth, delineated based on S-RFs sampling 10 regions across the globe, was primarily associated with Precambrian platforms that interacted with either Mesozoic or Cenozoic mantle plumes (Vinnik et al. 2007). Precluding the transition-zone-water-filter model, Vinnik & Farra (2007) suggest dehydration of water-bearing silicates as a plausible mechanism (Vinnik et al. 2007). Hydrous silicate melts can be gravitationally stable above d410, however, at their upper limits, they become more buoyant and saturate the olivine-bearing upper mantle (Andrault & Bolfan-Casanova 2022). Positive buoyancy of hydrous melts is revealed experimentally when the density and viscosity of melts are parametrized in a |${\rm MgO}{\text -}{\rm SiO}_2{\text -}{\rm H_2O}$| system as a function of temperature, pressure and composition. Neutral buoyancy can be achieved in certain cases with the addition of certain components like Fe (Drewitt et al. 2022). A global compilation reveals dehydration induced melting as the primary cause of the LVLs above d410 (Table 1). Most of the explanations invoke temperature/compositional changes in the presence of deep diving slabs or mantle plumes and their interactions with the surrounding mantle. Such LVLs are reported as intermittent features in various studies (Table 1) albeit suggestions that it is a global phenomenon (Tauzin et al.2010).

We used the most widely used technique of RFs to investigate deep LVLs beneath the Indian sub-continent. Although in a previous study, the presence of an LVL-410 at some locations beneath the Indian subcontinent was passingly reported (Kumar et al. 2013a), a thorough search for the same has not been hitherto attempted. In this study, we systematically analyse about 61 000 RFs sampling various tectonic domains of the Indian region to explore an LVL above the d410 discontinuity.

2 PROBING THE LVL-410 DISCONTINUITY

2.1 P-Receiver functions

Broadband data from 448 stations located on diverse tectonic units of the Precambrian Indian shield, Deccan Volcanic Province (DVP), Himalaya, Indo-Gangetic plains, Assam foredeep, Shillong Plateau (SP) and Burmese arc are assembled and processed for computation of RFs (Fig. 3). The seismic stations in India (IN: 233), including those in the Arunachal Himalaya (Singh et al. 2021), have been mostly operated by Council of Scientific and Industrial Research-National Geophysical Research Institute (CSIR-NGRI). However, the National Center for Seismology (NCS) maintains a number of stations as a part of the national network (Bansal et al. 2021). The data from stations in Nepal (YL: 29, Sheehan et al. 2001; XF: 182, Nábelek 2002) and Bhutan (XA: 5, Miller 2002) are accessed from the IRIS Data Management Center (Refer Supplemental Material Table S1 for detailed information). High quality waveforms of teleseismic earthquakes (epicentral distance of 30°–100°) with a signal to noise ratio (SNR) |$\ge$| 2.5 are used to compute RFs, using the advanced extended-time multitaper frequency domain technique (EMTRF, Helffrich 2006). The technique is an extension of the multitaper RF calculation method of Park & Levin (2000). The EMTRF method uses a series of short, overlapping, multiple tapers (10s window, 75 per cent overlap and three Slepian tapers) to window the time series across its length and sums the individual Fourier transformed signals to produce an RF estimate. In this manner, the relative amplitude of the Ps conversion and phase information for each subwindow are preserved. A frequency domain low pass cos2 taper with a cut off at 1 Hz (or 0.5 Hz) is applied to all the RFs to avoid Gibbs effect. A total of 61 248 RFs are eventually used for further analysis.

2.2 Bootstrap stacking of depth migrated receiver functions

In order to reliably detect weak conversions from second order discontinuities like the LVL-410, the RFs are migrated and binned in circular regions. First, effects of station elevation are removed and the RF amplitudes are projected along the ray path, using a 3-D velocity model derived from LITHO1.0 (Pasyanos et al. 2014). The steps followed are similar to those described in an earlier paper (Saikia et al. 2020). After migration, the RFs within circular regions having a radius of 75 km, centred at regular horizontal intervals of |$\sim$|55 km (0.5°), in a depth interval of 1 km (between 250 and 800 km) have been binned to enhance conversions from second order discontinuities like the LVL-410. Any bin having less than 25 RFs is discarded from further analysis. The conversion depths and associated uncertainties in each circular region are obtained by applying the bootstrap resampling technique of Efron & Tibshirani (1986). For each bin, at a particular depth, the RF amplitudes are randomly selected based on the number sampling. The whole procedure is repeated 200 times, using entirely different sequences of uniformly distributed pseudo-random numbers. Finally, using the mean and standard deviation of the amplitudes obtained from 200 iterations, a summation trace is constructed from 250 to 800 km, along with the corresponding uncertainties. The depths to the LVL-410 (if detected) and d410 discontinuities are picked from the summation trace. In addition, we performed differential slowness stacking of the RFs in bins where the LVLs are observed, in order to ensure that the observed amplitudes correspond to Ps conversions (Tauzin et al. 2010). For bins where d410 is not clear (e.g. Fig. 4c, 80.0°, 13.0°), the differential slowness stack confirms the presence of both features.

2.3 Observations

High quality moveout corrected P-RFs, arranged and stacked in narrow slowness bins in the range of 4.4 to 8.8 s deg⁻¹ are shown in Fig. 4 for grids spanning the Indian shield, Himalaya, Burmese arc and adjoining regions. In view of the large number of bins, they are grouped based on their tectonic affinity and/or geographical location and shown in different subplots. Fig. 4 reveals a clear signature from the d410 and a low velocity layer prominently seen for grids in the south Indian shield (Fig. 4a) and Himalaya (Fig. 4c). The bootstrap summation stacks corresponding to grids sampling the junction of the western ghats, Deccan Volcanic Province (DVP) and Southern Granulite Terrane (SGT) reveal strong conversions from an LVL at depths of |$\sim$|335 km (Fig. 5). However, such conversions from the LVL-410 are not detectable in a few neighbouring bins. The bins within and in the vicinity of the Eastern Ghats Mobile Belt also show weak arrivals close to d410 (Fig. 6). Interestingly, grids sampling the Shillong Plateau (SP), Mikir Hills and adjacent regions in northeast India, show strong conversions from the LVL-410, with significant lateral variations, straddling a depth range of 320–370 kms (Fig. 7). In all these cases, the veracity of LVL observations is verified using differential slowness stacks to ensure that they are Ps conversions.

The LVLs above d410 are detected in 845 out of the 1293 (nearly 65 per cent) bins spanning India. However, the depth to the LVL-410 varies largely, from 292 to 402 km (Fig. 8b, Table S2). We also observe a variability in the ratio of the Ps conversion amplitude from the LVL and the d410 discontinuity (Fig. 8a). Out of the 845 bins, 385 bins have amplitude ratios greater than one, implying a strong LVL-410. These LVLs are dominant beneath a majority of regions like the Cambay Rift Basin (CRB), parts of western Himalaya and foredeep, Burmese arc, regions south of the Cuddapah basin and within the SGT. In nearly 42 bins, the amplitude ratio is less than 0.5.

3 DISCUSSION

3.1 Robustness and variations of LVL-410

In India, most of the earlier studies pertaining to the sublithospheric mantle were focussed on investigating lateral variations in depths to the upper mantle discontinuities and amplitude variations of the Ps conversions from them (Yuan et al. 1997; Kosarev et al. 1999, 2013; Saul et al. 2000; Ramesh et al. 2005; Kumar & Mohan 2005; Oreshin et al. 2008; Shen et al. 2008; Kiselev et al. 2008; Singh & Kumar 2009; Oreshin et al. 2011; Singh et al. 2012; Kumar et al. 2013a; He et al. 2014). Previously, two global studies that included data from a solitary GEOSCOPE station (HYB) in India reported contrary findings. While Vinnik & Farra (2007) found no evidence for an LVL-410, Tauzin et al. (2010) found the same at depths close to 350 km. Subsequently, a continental scale study of the transition zone discontinuities documented intermittent presence of an LVL-410 beneath India in a few regions like northeast India, Himalayan foredeep and Eastern Ghat Mobile Belt, based on stacking of RFs at multiple stations sampling these (Kumar et al. 2013a). However, this study could not conclusively detect this feature beneath the plume affected southern and northern Deccan Volcanic Province. A recent study combining data from a new experiment in Sikkim Himalaya with the existing ones in the eastern Himalaya, reported strong negative amplitude Ps conversions in RFs sampling the northernmost Himalaya and southern Tibet, corresponding to a depth of 350 km (Kumar et al. 2023).

Our results show clear presence of an LVL-410 beneath various geological provinces of the Indian shield (Fig. 8), with Ps conversions from it being weak in some (Fig. 6) and strong in most other regions (Figs 5 and 7). Even those close to the d410 conversion appear realistic, since identical sidelobes are not seen on either side of d410 conversions. At most of the locations, the LVLs are distinguishable with strong amplitudes, well separated from other phases. The differential slowness stacks further made it possible to verify that the LVL conversions are robust.

The amplitude ratios of PLVL-410s/P410s are quite large in certain cases (>1, Fig. 8a), suggesting strong heterogeneities above the MTZ with a large decrease in Vs, compared to standard velocity models. Synthetic modelling suggests that nearly 9 per cent reduction in shear velocity is required to explain an amplitude ratio of 0.9 with no change in the velocity gradient. The same reduction (9 per cent) will result in smaller amplitude ratios when a gradient is introduced. For a 50 per cent gradient, the amplitude ratios are close to 0.7 for 9 per cent reduction in Vs. A large reduction in Vs (8.9 per cent) with a sharp velocity gradient (⁠|$\lt $|6.4 km) was required to explain LVLs beneath the northern Rocky Mountains (Jasbinsek & Dueker 2007). LVLs atop the MTZ beneath the central Asian orogenic belt were explained with a P-wave velocity drop of |$\sim$|5.7 per cent (Li et al. 2022b). This orogenic belt lies to the immediate north of the Qaidam basin and forms the largest accretionary orogen on Earth (Li et al. 2022b). In extreme cases, a shear wave velocity reduction of nearly 10 per cent was required to best fit the synthetics with the observed data in the Pacific Ocean (Wei & Shearer 2017). It is evident that sharp reductions are reported globally in different tectonic environments akin to our observations for the Indian subcontinent.

Beneath India, the depth of LVLs atop d410 is found to be quite variable, the shallowest being around 290 km (Fig 8b). In most cases, it seems to be within 50 km above the d410 (Fig. 8b). Such variations are commonly reported globally (Tauzin et al. 2010). Theoretical modelling suggests that all possible melt compositions in the |${\rm MgO}{\text -}{\rm SiO}_2{\text -}{\rm H_2O}$| ternary at different temperatures (3000 K and 1800 K) are less dense compared to the density values in the standard earth model (PREM, Dziewonski & Anderson (1981). The hydrous melts produced atop d410 due to transport of bulk |${\rm H}_{2}{\rm O}$| from the transition zone absorbed by olivine or olivine plus hydrous melt, are positively buoyant (Drewitt et al. 2022). The melts may pound atop d410 between a neutral buoyancy depth and d410, where melt may become denser than the solid (Andrault & Bolfan-Casanova 2022). Neutral melt-mantle buoyancy defines the upper depth limits of observed LVLs, after which the melt is positively buoyant and saturates the mantle towards the surface (Andrault & Bolfan-Casanova 2022). The neutral melt-mantle buoyancy layer may vary laterally and could possibly define the variation in depths of these LVLs (Fig. 8b). Drewitt et al. (2022) argued that the upper mantle cannot have a significantly different |${\rm H}_{2}{\rm O}$| content than the MTZ over larger geological time scales (mostly |$\gt $|1000 Ma). However, they do not deny a spatially heterogeneous and temporally episodic differential distribution of |${\rm H}_{2}{\rm O}$| in the upper mantle. The subduction driven models suggest variations in the |${\rm H}_{2}{\rm O}$| content which may affect generation of melts (Drewitt et al. 2022), probably at smaller time scales. Beneath the study region, such theoretical models find significance since they explain the observed variations in the potential drop in shear velocities across the MTZ. The Indian mantle is continuously fed by Tethyan slabs as the Indian Plate moves northwards (Van Der Voo et al. 1999).

3.2 Plumes and slabs affecting the MTZ

Various models proposed to explain the presence of deep LVLs require dehydration induced melting above the MTZ due to the entrainment of hydrated MTZ rock by plumes or return flow associated with a subducting lithosphere (Tauzin et al. 2010; Liu et al. 2016, 2018). The Indian sub-continent is an amalgamation of different geological domains, marked by tectonic episodes that affected it since Precambrian. We find evidence of an LVL-410 beneath the Indian subcontinent, which suggests temperature/compositional changes within the upper mantle. Two dominant possibilities which may alter the normal compositions of upper mantle are presence of Tethyan slabs and plumes like the Réunion and Kerguelen.

The rapid drift of the Indian Plate began nearly 130 Ma, after its break up from the Gondwana supercontinent (Cande & Stegman 2011). Closure of the Tethys ocean which separated the Laurasia and Gondwana supercontinents began with the Gondwana fragments moving northwards. The closure of Neo-Tethys ocean culminated with the collision of India with Asia and other continental blocks, extending from the Mediterranean to the Indonesian archipelago (Hafkenscheid et al. 2006). Tomographic images revealed large scale anomalies at depths deeper than 1000 km beneath India, that are postulated to be remnants of ancient oceanic Tethyan lithospheric slabs (Van Der Voo et al. 1999). The position of these anomalies in the deep mantle (Van Der Voo et al. 1999; Helffrich 2006), clearly suggest regular sinking of the Tethyan oceanic slabs into the Indian upper mantle during the northward voyage of India since |$\sim$|130 Ma. Subducting oceanic slabs can transport water deep into the mantle, with the peak water content in the MTZ reaching values of |$\sim$|0.1 to 1 wt per cent (Karato et al. 2020). When water rich mantle materials are transported out of the MTZ due to a convecting mantle, partial melting occurs. The LVL-410 beneath India could primarily result from dehydration induced melting due to an MTZ enriched by water transported by the Tethyan slabs. However, an important issue to ponder is the retention of water in the MTZ for prolonged periods since Cretaceous. Kuritani et al. (2011) argued in favour of a stable hydrated MTZ over longer time scales (⁠|$\gt $|1 Ga) beneath northeast China. Geochemical characteristics of basalts seem to be affected by at least two hydration events, ancient and recent, related to stagnation of ancient and recent Pacific slabs, respectively (Kuritani et al. 2011). Additionally, Karato et al. (2020) present mineral-physics mechanisms aiding greater hydration in mantle layers above and below the MTZ. Mineral phase transitions of partial melts at d410 and d660 result in a negative buoyancy above and positive buoyancy below the MTZ. Further, water partitioning into the partial melts causes melts to return to the MTZ lending to retention of greater hydration in a steady-state mantle convection. The mechanisms provide a scientific rationale for a wet MTZ in maintaining a global patchwork of LVL-410.

Strong amplitudes of Ps conversions from the LVL-410 are observed near the plume affected CRB, westernmost parts of the DVP and SGT (Fig. 8). Passage of the Indian Plate over the Réunion hotspot resulted in a rapid (⁠|$\lt $|1 Ma) eruption of basalts that constitute the DVP (Courtillot et al. 1986). Seismic tomographic images revealed slow P-wave anomalies (−1.5 per cent) in the mantle that were interpreted as imprints of plume activity (Kennett & Widiyantoro 1999; Singh et al. 2014). Also, the Ps conversion times from the MTZ discontinuities beneath the SGT and DVP are significantly delayed, by about 1.5 s (Kumar et al. 2013b). Frazer & Park (2021) have reported that an ascending dry mantle plume may entrain hydrous MTZ rocks from their surroundings and induce partial melting above the d410, in the context of a deep-rooted Yellowstone plume. Vinnik & Farra (2007) reported LVLs as deep as |$\sim$|350 km beneath the Precambrian platforms and correlated their association with Mesozoic or Cenozoic mantle plumes and dehydration melting of water bearing silicates. In India, strong signatures of an LVL-410 possibly reflect partial melting induced by the Réunion plume, also suggesting a greater degree of coupling between the mantle and overlying lithosphere which retains the effects of a plume since Cretaceous.

4 CONCLUSIONS

A synthesis of results from global studies affirms presence of deep low velocity layers spanning major continents and oceans, without tectonic affinity (Table 1, Fig. 2). Interpretations from various studies primarily highlight the role of water in the MTZ in the genesis of these LVLs. Transportation of water into the deep mantle is facilitated by subducting slabs and retained in the MTZ due to its large storage capacity (⁠|$\sim$|1 wt per cent, Fig. 9). Subsequently, the hydrated MTZ rocks release water once they are transported to regions of relatively lower water bearing capacity. Observations from this study populate the global data base of these LVLs above the 410 km discontinuity and emphasize the role of Tethyan slabs in enriching the water content of MTZ below India. In regions like the DVP and SGT, the water could have been entrained into the mantle by ascending plumes. Absence of LVLs in certain locales is possibly due to very thin melt layers above the d410, which may have been transparent to the RFs sampling them. Strong lateral variations in the depth and sharpness of the LVL-410 can be explained by a neutral melt-mantle buoyancy that defines the upper depth limits and the amount of water in the MTZ, as revealed by theoretical modelling.

AUTHOR CONTRIBUTIONS

M. Ravi Kumar (Investigated the work and wrote most sections of the manuscript, Conceptualization, Data curation, Project administration, Supervision, Validation, Writing–original draft, Writing–review & editing); Arun Singh (Synthesized the results and wrote parts of the manuscript, Data curation, Formal analysis, Investigation, Validation, Visualization, Writing–original draft) and Dipankar Saikia (Wrote the codes and performed the analysis, Data curation, Formal analysis, Investigation, Methodology, Software, Writing–original draft).

ACKNOWLEDGEMENTS

This work was performed under the Main Lab project of CSIR-NGRI. The IRIS data management centre is gratefully acknowledged for making the seismic data available. Figures are made using the GMT software (Wessel & Smith 1998). We gratefully acknowledge the thorough reviews of Prof Jeffrey Park, an anonymous reviewer and those of the Associate Editor and the Editor Prof Yao.

DATA AVAILABILITY

The RFs used in this study can be accessed from https://doi.org/10.5281/zenodo.12772424. Waveform data of IRIS stations can be downloaded from IRIS http://ds.iris.edu/ds/nodes/dmc/.

REFERENCES