Crust and upper mantle velocity structure of the northwestern Indian Peninsular Shield from inter-station phase velocities of Rayleigh and Love waves

We measure the inter-station Rayleigh and Love wave phase velocities across the northwestern Indian Peninsular shield (NW-IP) through crosscorrelation and invert these velocities to evaluate the underneath crust and upper mantle velocity structure down to 400 km. We consider a cluster of three stations in the northern tip of the Peninsula and another cluster of eight stations in the south. We measure phase velocities along 28 paths for Rayleigh waves and 17 paths for Love waves joining two stations with one from each cluster and using broadband records of earthquakes which lie nearly on the great circle joining the pair of stations. The phase velocities are in the period range of 10 to 275 s for Rayleigh waves and of 10 to 120 s for Love waves. The isotropic model obtained through inversion of the phase velocities indicates 199.1 km thick lithosphere with 3layered crust of thickness 36.3 km; the top two layers have nearly same velocities and both constitute the upper crust with thickness of 12.6 km. The upper crust is mafic, whereas the lower crust is felsic. In the mantle lid, velocities increase with depth. The velocities of mantle lid beneath NW-IP is lower than those beneath south Indian Peninsula showing the former is hotter than the later perhaps due to large Phanerozoic impact on NW-IP. The significant upper mantle low velocity zone beneath NW-IP indicates high temperature which could be attributed to the past existence of a broad plume head at the west-central part of the Peninsula.


Introduction
The Indian Peninsula is a combination of a few Precambrian cratonic nuclei and is one of the oldest Archean shield regions of the world. The amalgamation of the cratons occurred at the end of the Archean and all the nuclei stabilized at about 2.5-2.6 Ga [Meert et al. 2010]. The juxtaposition along the Central Indian Tectonic zone (CITZ, Figure 1) took place during the earliest Paleo-Proterozoic [Stein et al. 2004]. Radhakrishna [1989] proposed CITZ as an ancient suture between north and south protocontinents. The south block consists of the north Dharwar craton, its northern part is covered by the Deccan trap flows and heart of the Deccan Volcanic Province (DVP), which is result of interaction between the mantle plume and the overriding continental lithosphere at ~65 Ma [White and Mckenzie 1989]. The north block consists of the Aravalli-Bundelkhand craton which constitutes the basement beneath the Lesser Himalayan successions, and the Quaternary alluvium of the Indo-Gangetic plain [Kailasam 1976]. To the northwest of this block, the Rajasthan orogenic belt had evolution in the Archean and continued to the end of Phanerozoic. The north and south cratonic blocks have a dominantly Proterozoic cover with properties of Archean lithosphere [Mahadevan 2013].
Early use of surface waves to find lithospheric velocity structure beneath Indian region has been reviewed by Bhattacharya [1992]. However, such studies for 1D structure were based on analog data and generally of periods < 100 s. With the availability of broadband digital data, we can now prepare 2D maps of surface wave velocities; however, such maps have been prepared with period only up to of 70 s [Acton et al. 2010]. Inter-station phase velocity measurements using waves from distant earthquakes generate longer period surface wave velocity data and such long period data allow us to evaluate not only the lithosphere but also deep structure of the uppermost mantle; further, the measured inter-station velocities are nearly independent of the errors in epicenter and origin time determinations [Mitra et al. 2006, Bhattacharya et al. 2013.
Many cratons over the globe are underlain by lithospheric roots which extend to depth > 250 km in contrast to the oceanic area and Phanerozoic conti-nents [Mahadevan 2013]. These deep roots may perturb the surrounding mantle and the mantle convection [Thompson et al. 2011]. However, the loss of these roots beneath Indian cratons might have been caused by the Deccan plume [Lehmann et al. 2010]. Thus the topology of Lithosphere-Asthenosphere boundary (LAB) i.e. the bottom of the mantle lid, beneath the Indian Peninsular shield will allow us to model its formation and evolution. Further studies of upper mantle velocity structure beneath the Peninsular shield using long period surface waves are limited [Mitra et al. 2006]. The topology of LAB has been estimated by Priestley and Mckenzie [2006] by converting S-wave velocity (V S ) into temperature profiles and then fitting the geotherms beneath Eurasian plate. They showed LAB depth is ~160 km in central India and increasing northward to ~250 km in the northern tip of the Indian Peninsula. In recent studies, lithosphere thickness beneath the Indian Peninsular shield has been estimated with variations of 70 to 140 km thick on one side using receiver function studies [Kumar et al. 2007, Kumar et al. 2013a; and 150 to 200 km on the other side using Rayleigh wave phase velocities for periods up to 200 s by Mitra et al. [2006], gravity studies by Arora et al. [2012], and receiver function analysis by Bodin et al. [2014]. These results underline the need for studying the uppermost mantle structure further.
So far the upper mantle structure of the northwestern part of the Indian Peninsular shield is poorly resolved due to limited broadband earthquake records in the northwestern part of the shield region. With the installation of a few broadband stations recently, we measure inter-station phase velocities of both Love and Rayleigh waves to evaluate the velocity structure of the northwestern Indian Peninsular shield (NW-IP) region down to 400 km. Surface wave dispersion mainly depends on S-wave velocity structure, which is reliably obtained from surface wave studies. The dispersion depends least on density. Left aside density, Rayleigh wave dispersion depends both on P-wave velocity and S-wave velocity structure, while Love wave dispersion depends on S-wave velocity structure only. Thus use of both Love and Rayleigh waves, as has been done here, will give a reliable result on S-wave velocity structure. Further, a discrepancy to fit an isotropic elastic model to both Love and Rayleigh wave dispersion will show a radial anisotropy, if any, in the structure. The derived results are used to study (i) crustal structure (ii) nature of mantle lid and depth of LAB; and (iii) characteristics of upper mantle LVZ.

Data
We consider records of broadband seismic stations operated by the India Meteorological Department (IMD) and the National Geophysical Research Institute (NGRI) to measure inter-station surface wave dispersion over broad period range. We consider two clusters of stations: the south cluster, consisting of 8 stations, is located in Koyna-Warna region in the NW Dharwar craton (in the western Peninsula) and the north cluster, consisting of 3 stations, is in the northern tip of the Peninsula (Figure 1). The information on the stations is given in Table 1.
For inter-station phase velocity measurement, we consider the records of two earthquakes (Table 2) each of which lies nearly in the same great circle joining a pair of stations with one from each cluster. The angle between the two great circles (1) joining the two stations and (2) the longer of the two great circles joining each station and the epicenter, is considered to be less than 3° [Bhattacharya et al. 2013], in order to minimize the influence of the structure between the earthquake and nearest station; this also minimizes the errors in azimuth arising from refraction, but such errors introduce only second order effects in phase velocity [Yao et al. 2005]. We process 28 inter-station wave paths for Rayleigh waves and 17 paths for Love waves (Table 3) after careful examination for good quality records consisting of least scattered waves and only those waveforms having a signal to noise ratio ≥ 2.5. We consider each wave path between one of the stations of north cluster and another of the stations of south cluster (Table  3). These wave paths are close to each other and mainly cross north-Dharwar craton, narrow Central  Indian Tectonic Zone (CITZ) in the middle and then northwestern cratonic block. Epicentral distances of event 1 (Table 2) are around 1770 km from the northern cluster of stations and around 3070 km from southern cluster of stations. While those of event 2 (Table 2) are around 10,560 km from southern cluster of stations and around 11,860 km from northern cluster of stations.

Measurement of inter-station Phase velocities
Measurement of inter-station phase velocities is based on cross-correlation method [Yao et al. 2005]. A brief description of the methodology implemented through MATLAB is as follows: (a) Recorded data are corrected for its instrument response to get true ground motion records.
(b) The N-S and E-W components are rotated to obtain radial and transverse components; we use the vertical component for Rayleigh wave and the transverse component for Love wave.
(c) For each seismogram, we obtain group arrival times (t gi = D/U i ) of a given period T i , where D is epicentral distance and U i is the group velocity at that period T i . We obtain U i through frequency time analysis (FTAN: Dziewonski et al. [1969], Bhattacharya [1981Bhattacharya [ , 1983). Now consider two stations with a pair of vertical seismograms (for Rayleigh wave) and another pair of transverse seismograms (for Love wave). Using a given pair of seismograms, we evaluate inter-station phase velocity at period T i as described in the next three steps: where m > n. We have normally selected m = 4.5 and n = 4.0; we have kept a small difference between m and n to have small tapering parts on both sides of the window. This window effectively removes the noise of higher mode waves and other spurious signals.
(e) Each windowed seismogram is now filtered with a narrow frequency band centered at frequency 1/T i . For filtering, we use wavelet transform with Morlet wavelet [Wu et al. 2009].
(f ) Finally we cross correlate a pair of such filtered seismogram, one from the north cluster and another from the south cluster. The correlated waveform shows maximum amplitude (crest) at a time when the phase difference is zero; this time corresponds to the phase arrival time and is used to compute phase velocity (= inter-station distance/ arrival time) of period T i . Thus, computed phase velocities for various periods have been plotted as phase velocity image. On this image, at higher periods phase velocities marked by amplitude peaks (crest) are easily discernible and phase  (1) velocities of available period range are obtained in a continuous manner. As an example, Figure 2 shows a pair of vertical component of seismograms of NDI and SKP from the event 2; the cross correlated plot of this pair as in the steps above is shown in Figure 3.

Observed data and comparison with data of the Peninsular India
The measured phase velocity extends for periods from 10 to 275 s for Rayleigh waves and from 10 to 120 s for Love waves (Table 4, Figure 4). We could get phase velocities of surface waves of lower period range (10 to 52 s) from event 1 and higher period range (26 to 275 s) from event 2. Phase velocities of Love waves are restricted up to 120 s due to insufficient signal-to-noise ratio at longer periods. Love wave phase velocity with period is not so smooth ( Figure 4); from 70 s onward the phase velocity is nearly flat with an upward jump at 100 s. Such anomalous part may not be due to noise, because we have obtained these data with many wave paths and the standard deviations at this part are small. We feel that this part is likely to be due to some lateral variation along the wave path. Yoshida [2000] has demonstrated that a lateral variation of Moho depth (at about 35 km) gives rise to a flat phase velocity of Love wave around the period 30 s. In the present case the flat phase velocity observed between 74-100 s indicates that such a lateral variation along the wave path is occurring at a greater depth. As per the partial derivatives of phase velocities of Love waves with respect to S-wave velocity , this period range is most affected by a layer at a depth ~200 km. A down going LAB was noted by Priestley and Mckenzie [2006] from south to north across CITZ. Such anomalous part of phase velocity for Rayleigh wave is not discernible because effect of lateral heterogeneity is small on Rayleigh waves than on Love waves [Lavender 1985]. We shall model the phase velocity data with a laterally homogenous structure assuming that the effect of lateral heterogeneity in a small period range may not affect the result significantly.
For the Bastar craton, on the eastern part of the Peninsula, crust and upper mantle structure was obtained earlier through inversion of inter-station phase velocities of Rayleigh and Love waves [Bhattacharya et al. 2009]. The theoretical curves for this model closely agree with data of Love wave, but not that of Rayleigh wave (Figure 4). The theoretical phase velocity curves are obtained here for an isotropic elastic layered spherical Earth following Bhattacharya [1996Bhattacharya [ , 2009 Table 2). The start time of both the waveform is 03h 04m 08.82s. The phase velocities obtained through cross correlation of this pair of seismograms are shown in Figure 3.  (Table 2). White and black areas show positive and negative amplitude part in the cross correlated waveform at each period. A peak (crest) is shown by a symbol, which is properly chosen as phase velocity.

Inversion and results
We perform a nonlinear iterative search of a model using genetic algorithm (GA) [Suresh et al. 2008]. In GA, we obtain a model minimizing the misfit between observed and theoretical phase velocities. We define misfit as  theoretical phase velocity of a given model at period T i (i =1, 2, …, NR) for Rayleigh wave; d j L are corresponding values of period T j (j = 1, 2, …, NL) for Love wave. We have used misfit as the maximum of misfitR and misfitL to give equal importance to Rayleigh and Love waves. In this case, in the GA solution of a model, the values of misfitR and misfitL are nearly equal.
We consider a layered structure down to 400 km with 3 layered crust, below which there are 8 layers. The model parameters V S , V P /V S and thickness in each of the top 10 layers are considered as variables with possible upper and lower limits as search ranges (Table 5).
The GA begins with a random K models. In each successive generation, we create K new models through crossover and mutation as well as through a few (e.g. 2) elite models, which replace the worst models of current generation of models with the best individuals of previous generation. After L number of generations, we retain the best solution (model). With 30 variables (Table 5), we consider K = 600 models during each generation and best solution is obtained after 200 generations. In each GA operation, we get best solution for 30 variables. GA operation is repeated 24 times and mean of the 24 best solutions is considered as the final accepted solution ( Table 5). The theoretical curves for the accepted model are in Figure 4; the misfit { for this model is 0.0288. The model parameters along with standard deviations of the accepted model are given in Table 5. The accepted model shows that LAB is the top of the 7th layer. From each of the 24 solutions, we have obtained the depth of LAB as 199.1 ± 9.6 km.
It is well known that Rayleigh wave dispersion largely depends on V S compared to V P and Love wave dispersion is independent of V P . The dependence of Rayleigh wave phase velocity to V S and V P at different depths of the accepted model is shown by relative sensitivity in Figures 5 and 6. The relative sensitivity of phase velocity of Love waves to V S is shown in Figure  7. The relative sensitivity of Love wave is maximum in the crust and decreases with depth. The crustal V S affects phase velocity of Rayleigh waves of period below 30 s; above this period, Rayleigh wave sensitivity to V S has peaks in the uppermost mantle ( Figure 5). The sensitivity to V P has maximum in the crust at all periods and decrease with depth ( Figure 6). In Figures 5 and 6, we have plotted relative sensitivity to V S and V P only for 5 representative periods. We may also draw sensitivity for all the 25 periods (Table 4) considered for GA inversion and get two figures each with 25 curves instead of 5 as in Figures 5 and 6: one for V S and other for V P . From each figure, for a given depth, we evaluate the maximum sensitivity among all 25 curves. Figure 8 shows the depth-wise variation of maximum sensitivities to V S and also to V P . In Figures 5 and 6 each curve shows the relative sensitivity of a particular period. While the two curves of maximum sensitivities in Figure 8 show the relative sensitivities to V P and V S in the period range 10-275 s. In this figure we also plot depth-   wise variation of ratio of maximum sensitivity = (maximum sensitivity to V P )/ (maximum sensitivity to V S ). At depths below 10 km the ratio is about 0.1. Thus for Rayleigh waves, nearly at all depths the sensitivity to V P is about 10 per cent of the sensitivity to V S in the period range considered here. This shows that the dependence of Rayleigh wave phase velocity to V P is relatively small but not negligible. Thus using the limits of V P /V S , we have inverted V P , whose determination is less reliable than that of V S .

The crust
Along the present wave paths we have considered a 3-layered crust namely shallow, upper and lower crust based on the findings for DVP [Prajapati et al. 2011, Suresh et al. 2014]. The V S as well as V P in the top two layers are close to each other (Table 5); thus considering these top two layers as an upper crust, we find the thickness of this part of the crust as 12.6 km, while the lower crust is 23.7 km. Such thin upper crust is seen in a northsouth profile crossing whole of the Indian Peninsula [Bhattacharya 1974]. Further, Rao et al. [2002] obtained 13.8 km thick upper crust for central India in a grid search method through ray tracing using the data of body waves at regional distances from 1997 Jabalpur earthquake occurred in Central India. The southern part of the present wave paths are close to those used by Prajapati et al. [2011] across NW-DVP, where upper crust is 15.2 km thick; in addition there is a top Deccan Trap layer of thickness 1.6 km ( Figure 9). Here, we have ignored the thin Deccan Trap since the present wave  (Table 5) at various representative periods. The sensitivity shows partial derivatives ∂c/∂Vs in a layer of 1 km centered at the given depth.  (Table  5) at various representative periods. The sensitivity shows partial derivatives ∂c/∂V P in a layer of 1 km centered at the given depth.  (Table 5) at various representative periods. The sensitivity shows partial derivatives ∂c/∂Vs in a layer of 1 km centered at the given depth.
paths are partly across DVP where the Deccan Trap is thinner; further, the periods considered here are 10 s and above where the effect of a top thin layer may not be significant.
In the upper crust of NW-IP, the V P /V S is between 1.83 and 1.85. Such mafic crust was noted in adjoining West Ganga Basin [Mitra et al. 2011]. However, at NW-IP, the lower crust is felsic (V P /V S~1 .68). Felsic lower crust has also been seen earlier in the Indian Peninsular by Bhattacharya ([1981], IP11, V P /V S = 1.68), Singh et al. ([1999], V P /V S =1.60) and Rao et al. ([2002], V P /V S = 1.70). Zhou et al. [2000] explained such felsic crust in central India as possible compositional differences between Archean and Proterozoic crust.
The lower crust has V S (= 3.87 km/s) similar to that of NW-DVP (V S = 3.89 km/s). Such V S at lower crust is higher than the corresponding values in the adjoining regions like West Ganga Basin [Mitra et al. 2011] and Indus block crust [Suresh et al. 2008], where the crust is thicker due to sedimentary deposits ( Figure 9).

The mantle lid
Down to 164 km, the mantle lid of NW-IP is similar to that of the Bastar craton ( Figure 10). Below this depth there is a slight increase of V S down to LAB, which is at 199 km depth. Such increase is consistent with other observations in cratons, where V S reach higher values with depth in the lithosphere [Lebedev et al. 2009]. For south-IP (S-IP), a significant high velocity CRUST-MANTLE OF NW INDIAN PENINSULA Figure 8. Rayleigh wave maximum relative sensitivities to V S and V P at a given depth considering sensitivities of all the periods (Table 4) considered in inversion. These two curves show relative sensitivities in the period range 10 to 275 s. The figure also shows depth-wise variation of ratio of maximum sensitivities = (maximum sensitivity to V P )/ (maximum sensitivity to Vs). At depths below 10 km the ratio is about 0.1.  was noted in mantle lid with LAB at 155 km depth [Mitra et al. 2006]. Thus the mantle beneath S-IP is colder than that beneath NW-IP. The Phanerozoic impact was feeble in S-IP compared to NW-IP making S-IP lighter and colder lithosphere [Mahadevan 2013]. In the recent study by Bodin et al. [2014] clarified that the LAB at depth about 110 km obtained from S receiver functions by Kumar et al. [2007] and Kumar et al. [2013a] is a mid-lithospheric discontinuity and the actual LAB is seen at depth between 150 and 200 km in central India. In the western side of Dharwar craton, gravity data shows LAB is as deep as 165 km [Arora et al. 2012]. Considering maximum negative V S gradient as an indicator of LAB, the velocity structure of the Bastar craton shows LAB depth as 160 km beneath this craton ]. While in NW-IP, we find LAB is at 199.1 ± 9.6 km depth. Thus, the LAB depth is increasing from south to north agreeing with the findings of Priestley and Mckenzie [2006]. In the mantle lid beneath NW-IP, V P /V S~ 1.69. Low values of V P /V S~ 1.70 are noted in the mantle lid beneath central India [Singh et al. 1999] and beneath southern Africa [Wang et al. 2008]. Low velocity ratios were also noted in the stable cratons around the Tibet plateau [Pei et al. 2011]. In the mantle lid, lower aluminium content might have caused increase in clinopyroxene and orthopyroxene and decrease in garnet and thus lowering velocity ratio in the mantle lid [Wang et al. 2008]. Since an isotropic model satisfied both the measured Rayleigh and Love wave velocities, the radial anisotropy in the Peninsular shield is nearly absent in the crust and uppermost mantle.

The LVZ of upper mantle
Beneath NW-IP, the 164 km thick upper mantle LVZ extend from 199 to 363 km depth ( Figure 10); the lowest velocity in LVZ is 4.383 km/s, which is 6 per cent decrease of V S in LVZ. Although the LVZ thickness closely agrees with other cratons but the decrease of V S in LVZ is slightly larger compared to 4-5 per cent beneath other cratons Mckenzie 2006, Lebedev et al. 2009]. Kumar et al. [2013b] reported delays in P-to-S conversions from the 410 km discontinuity below the Indian shield and suggested low shear wave speeds in the lithospheric and sub-lithospheric mantles due to higher temperatures. Approximately 10 per cent reduction of V S in LVZ has been noted for Tanzanian craton [Weeraratne et al. 2003]; such reduction indicated high temperature and the presence of melt which are due to plume head beneath Tanzanian craton. For the Indian Peninsular shield, Rao and Lehmann [2011] postulated past existence of a large plume head of diameter 2000-2500 km centered at the middle of the Peninsula. Partial melt only decreases V S and not V P [Wang et al. 2008]. In LVZ beneath NW-IP, we find V S has decreased along with V P (Table 5). Thus we do not consider the existence of partial melt in LVZ, which is hot due to past existence of a plume head. In LVZ also we find V P /V S~1 .69, which is below normal value. Low velocity ratio may be caused by lower aluminium content as indicated above for the mantle lid.

Conclusions
1. The 3-layered crust is 36.3 km thick. The top two layers have velocities close to each other and these two layers constitute upper crust, which is 12.6 km thick and is mafic. However, the lower crust is felsic.
2. The theoretical curves of the isotropic elastic structure satisfy the observed phase velocities of both Rayleigh and Love waves. This indicates that the radial anisotropy either is absent or weak in the crust and uppermost mantle.
3. The lithosphere is 199.1 km thick and V S increases downward below the crust. However, V S in mantle lid of NW-IP is lower than S-IP. Relatively low V S in mantle lid shows its high geothermal gradient. This low V S shows that the lithosphere in NW-IP is hotter than S-IP perhaps due to large Phanerozoic impact on NW-IP compared to S-IP.
4. Strong LVZ is interpreted as positive thermal gradient due to which thermal mantle convection may not occur [Eaton et al. 2009]. The high temperature in LVZ is likely due to past existence of a broad plume head centered at middle of the Indian Peninsula.
5. The important result of crust and upper mantle structure obtained here will help to study the evolution and stability of lithosphere and to obtain upper mantle thermal convection model beneath cratons of NW-IP.