Seismic evidence of crustal heterogeneity beneath the northwestern Deccan volcanic province of India from joint inversion of Rayleigh wave dispersion measurements and P receiver functions

Journal of Asian Earth Sciences 128 (2016) 54–63

Contents lists available at ScienceDirect

Journal of Asian Earth Sciences journal homepage: www.elsevier.com/locate/jseaes

Full length Article

Seismic evidence of crustal heterogeneity beneath the northwestern Deccan volcanic province of India from joint inversion of Rayleigh wave dispersion measurements and P receiver functions A.A. Deshpande, G. Mohan ⇑ Department of Earth Sciences, Indian Institute of Technology Bombay, Powai, Mumbai 400076, India

a r t i c l e

i n f o

Article history: Received 6 March 2016 Received in revised form 15 July 2016 Accepted 20 July 2016 Available online 20 July 2016 Keywords: Crustal structure Deccan Rayleigh waves Receiver function Shear velocity Upper mantle

a b s t r a c t The northwestern Deccan volcanic province (NWDVP) of India, encompassing the Saurashtra peninsula and the adjoining Gulf of Cambay, is investigated through joint inversion of surface wave dispersion measurements and teleseismic P receiver functions, to estimate the crustal and shallow upper mantle shear wave velocity (Vs) structure. The Mw
7.7 Bhuj earthquake and the post Bhuj regional events, recorded during the period 2001–2010 at 7 stations along 37 source-receiver paths were used along with 35 teleseismic events. A joint curve fitting inversion technique is applied to obtain a best fit for the fundamental mode Rayleigh wave group velocity dispersion curves for time periods 5–50 s and high quality crustal P wave receiver functions obtained at each station. Significant crustal heterogeneity is observed within the study region with the average crustal Vs ranging from 3.5 km/s to 3.8 km/s with the paths cutting across the Gulf of Cambay exhibiting large reduction in shear wave velocities. Utilizing the average crustal Vs  3.66 km/s estimated for Saurashtra, together with the average crustal P wave velocity (Vp)  6.54 km/s derived independently through deep seismic sounding studies, yields a bulk Vp/Vs ratio of 1.786 or an equivalent crustal Poisson’s ratio of 0.271. A major contribution to the high Poisson’s ratio comes from the 12 to 16 km thick lower crustal layers with shear velocities ranging from 3.8 km/s to 4.19 km/s suggesting widespread magmatic underplating due to emplacement of mafic cumulates in the lower crust. The shallow uppermost mantle shear velocities are in the range 4.2–4.5 km/s averaging 4.36 km/s, which is less than that observed for the Indian shield, indicating the effects of residual thermal anomaly. The variation in the crustal Vs, high Poisson’s ratios and low upper mantle shear velocities reflect the thermal and compositional effects of the Deccan volcanism which are manifested in terms of pervasive presence of mafic dykes, volcano plutonic complexes, lower crustal mafic cumulates and possible partial melt due to elevated geotherms. The significant crustal alterations are possibly a result of weakening of the lithosphere by a mantle thermal source, leading to a thin and weak lithosphere which facilitated the Deccan volcanism through the rift zones of northwestern India. Ó 2016 Published by Elsevier Ltd.

1. Introduction The northwestern Deccan volcanic province (NWDVP) of India (Fig. 1) comprises the Saurashtra horst bounded by the pericratonic rift systems of Kutch, Cambay and Narmada, all of which have varied tectonic history and have evidenced crustal emplacement of magmatic melts. The 65 Ma old Deccan volcanics are predominantly composed of tholeiitic basalts and minor volumes of alkali basalts with carbonatitic lavas and intrusions exposed mostly along rift zones (Sen et al., 2009). Saurashtra is marked by several mafic intrusions in terms of volcanic plugs aligned along ⇑ Corresponding author. E-mail address: [email protected] (G. Mohan). http://dx.doi.org/10.1016/j.jseaes.2016.07.022 1367-9120/Ó 2016 Published by Elsevier Ltd.

the coast (Fig. 1) which are manifestations of the Deccan volcanism, which is hypothesized to be a result of the passage of the Indian plate over the Reunion hotspot (White and Mckenzie, 1989). The composition of the crustal intrusions is varied with several volcanic plugs oriented in the E-W trend of the Narmada rift composed of mafic/ultramafic to felsic rocks while those following the N-S trend of the Cambay rift are composed of alkaline rocks (Chandrasekhar et al., 2002). While the major intrusives (Fig. 1) were demarcated through prominent gravity and magnetic anomalies, several diffuse magnetic anomalies and circular and semi-circular gravity highs could also be linked to buried intrusives which are not delineated and were predating the Deccan event by 3 Ma suggesting a possible incubation (Chandrasekhar et al., 2002). It is expected that such crustal heterogeneities

A.A. Deshpande, G. Mohan / Journal of Asian Earth Sciences 128 (2016) 54–63

55

uppermost mantle beneath the Saurashtra peninsula and the Gulf of Cambay, and, characterize the alterations caused by the Deccan volcanism in northwestern India. 2. Data

Fig. 1. Geological map of northwestern Deccan volcanic province of India comprising of the Saurashtra peninsula and the adjoining rifts. The earthquake source (star) to receiver (inverted triangle) paths are shown for all seven stations. The deep seismic sounding (DSS) profiles are Navibunder (N) – Amreli (A) in Saurashtra and Mehmadabad (M) – Billimora (B) along the Cambay rift. The deep boreholes are located at Lodhika (L) and Dhanduka (D). Inset map of India shows the study region.

translate to variation in the crustal properties which would be evidenced through variations in the seismic velocity structure. Our knowledge about the subsurface P-wave velocity structure in NWDVP is primarily limited to the deep seismic sounding (DSS) studies along an E-W profile across Saurashtra (Kaila et al.,1980; Rao and Tewari, 2005) and a N-S profile along the Cambay rift (Kaila et al., 1981). The P wave receiver function (RF) studies (Praveen Kumar and Mohan, 2014; Rao et al., 2015) in Saurashtra reveal significant variations in the crustal thicknesses and Poisson’s ratios. The upper mantle shear velocity structure estimated through surface wave dispersion study (Prajapati et al., 2011) along the source-receiver paths cutting across northwestern DVP and southcentral DVP is similar to that observed beneath the Indian shield. Traditionally, the surface wave dispersion data and the receiver functions are analyzed and inverted separately to estimate the shear velocity structure with depth. Surface waves are sensitive to the average shear wave velocities while the receiver functions are sensitive to the shear wave velocity contrast and the vertical travel times. Surface wave inversions suffer from a vertical velocity depth trade off which can be minimized by using other seismic data e.g., Receiver functions, which provide a constraint. Combining the surface wave and receiver function data sets can effectively bridge the resolution gap associated with individual data sets. In this study, we attempt to carry out a joint inversion of both the surface wave and receiver function data sets (Ozalaybey et al., 1997; Julia et al., 2000, 2009) to investigate the crustal and shallow

Three component broadband stations were operated in standalone mode in phases in DVP during the period 1999–2010 (Fig. 1). Initially, the stations at MULG and VARE were established during 1999–2003. Later, stations were established in Saurashtra at BHOD, KHER, KAND, SONT and TANA during 2004–2010. All the stations were equipped with 120 s CMG 3-T Broadband seismometers and 24-bit recorders which recorded the earthquake data along with background noise at a sampling rate of 50 sps. The data used for surface wave dispersion measurements are the records of (i) the 26th January 2001 earthquake of Mw
7.7 at Bhuj, (ii) the aftershocks of the Bhuj earthquake, and, (iii) a few regional events in the vicinity of Bhuj, which occurred during the period 2001–2010. The Bhuj earthquake occurred along an ENE-WSW trending south dipping reverse fault at a shallow depth of 16 km (USGS) in the continental intraplate Kutch rift. Studies conducted to explain the source parameters indicate the existence of an approximately 1300 sq. km shallow rupture zone (Bodin and Horton, 2004) giving rise to a reverse thrust fault. The seismically active Kutch region continues to generate mild to moderate earthquakes (Mandal and Horton, 2007). A total of 37 regional events originating around Bhuj, with magnitudes (mb) P 4, were used for the dispersion studies (Fig. 1). Thirty-five high quality teleseismic events with magnitude (mb) P 5.5 in the epicentral distance range of 30°—90° were used for constructing the P receiver functions. 3. Methodology 3.1. Surface wave dispersion Fundamental mode Rayleigh waves are processed using the Multiple Filter Technique (MFT) (Dziewonski et al., 1969) with a routine constructed by Herrmann and Ammon (2002), to measure the group velocities for periods ranging from 5 to 50 s. For a single source-receiver geometry, the arrival times of different wave packets is picked as a function of frequency. Gaussian band-pass filters are then applied to isolate the various modes of dispersive signals. The parameter for this filter is selected based on the length of the path and longer the path, larger is the filter parameter due to the increase in dispersion envelope. The maxima of the envelope of the filtered signal is a marker for the arrival time of the wave packet and is used to determine the velocity of the wave group. This is achieved by computing the quadrature spectrum of the signal and computing the instantaneous amplitude for a given frequency. By repeating this procedure for different frequencies, the group velocity curve is obtained. The fundamental mode is picked manually from a group velocity and spectral amplitude plot. The dispersion measurements were made using the regional events occurring locally around Bhuj for the source–receiver paths to seven stations (Fig. 1). The short paths from Bhuj to the four stations BHOD, KAND, TANA and SONT sample different segments of the Saurashtra peninsula and the long paths from Bhuj to the three stations MULG, KHER and VARE cut across the Gulf of Cambay. The source-receiver paths to each station varied from a minimum of three paths for the pair Bhuj-BHOD to a maximum of 19 paths for the pair Bhuj-MULG. The source-receiver combinations are clustered by station, in order to obtain the average range of group velocities along the paths which lie within a narrow azimuthal range at each station. The dispersion plots for the source receiver

56

A.A. Deshpande, G. Mohan / Journal of Asian Earth Sciences 128 (2016) 54–63

Fig. 2. Fundamental mode Rayleigh wave group velocity dispersion measurements for several source-receiver paths in a narrow azimuthal range at each station. The best fit obtained through surface wave inversion is shown in red color. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

paths to the seven stations are shown in Fig. 2. The observed dispersion curves indicate large variations in the group velocities for the shorter periods up to 25 s revealing the heterogeneous nature of the shallow crustal structure. The best fit curve, indicated by the dispersion curve in red color in Fig. 2, is obtained after picking the fundamental mode dispersion curve from the MFT process and clubbing the multiple values observed at a single station together and inverting these values for a velocity-depth model until a maximum signal power fit for the observed points is reached. This best fit curve is then used as the station representative input curve for the joint inversion of surface wave dispersion curves and receiver functions. 3.2. Receiver functions Receiver function technique is a powerful method to investigate the shear wave velocity structure beneath a station using teleseismic waveforms (Langston, 1979; Owens, 1984; Ammon et al., 1990). P wave receiver functions (RF) can be considered to be source-equalized waveforms comprising of the direct P wave followed by the P to s converted phases (Ps) resulting from the interaction of the incoming P wavefront with near receiver subsurface velocity discontinuities. The RFs also comprise of two reverberations between the free surface and the discontinuity for each velocity interface. In this study, 35 teleseismic waveforms sampled at 50 Hz with high signal to noise ratios P 5, corresponding to earthquakes of mb P 5.5 in the epicentral distance range of 30°–90° corresponding to the slowness range 8.8–4.4 s/deg were used to compute the RFs. The back-azimuthal distribution of the events used is shown in Fig. 3. The waveforms were windowed 5 s before and 25 s after the P wave arrival, detrended, tapered, and decimated to a uniform sampling rate of 20 samples per

second, after low-pass filtering below 5 Hz to avoid aliasing. An additional high pass filter of 0.01 Hz is applied to the traces for improving the stability of the RFs. The instrumental response was deconvolved from the waveforms and based on the backazimuth the waveforms in the ZNE (vertical, north, east) were rotated around the vertical component into radial and tangential components i.e., ZRT system. An iterative time domain deconvolution procedure (Ligorria and Ammon, 1999) was performed to deconvolve the vertical component from the corresponding radial components to yield the receiver functions. The time domain deconvolution procedure applies a Gaussian low-pass filter of width 2.5 (f < 1.2 Hz) to the waveforms prior to deconvolution to limit the bandwidth of the resulting receiver functions. The radial component of the receiver function is plotted as a function of ray parameter for each station (Fig. 4). The arrival time, amplitude and polarity of the Ps converted phases and the associated multiples are a function of the crustal parameters. The variations in the delay time of the P to s conversion on the RFs corresponding to three different ray parameters at all stations are shown in Fig. 4. The RFs are marked by two prominent positives corresponding to the P to s conversions from a shallow and a deep velocity interface. The positive Ps conversion observed within 0.5 s on the RFs at all stations corresponds to the conversion from the base of the Deccan traps. The Pms phase represented by a positive peak corresponding to the P to s conversion from the Moho is observed around the 5 s mark. In order to model the crustal structure from RFs, the procedure adopted is to generate synthetic RFs based on an assumed velocity model with fixed or variable thicknesses, P and S wave velocities, and densities respectively, and obtain a best fit between the observed and synthetic RF. In this study, the RFs are used as a constraint to jointly invert the surface wave dispersion curves and the receiver functions as described below.

A.A. Deshpande, G. Mohan / Journal of Asian Earth Sciences 128 (2016) 54–63

57

Fig. 3. Back azimuthal distribution of teleseismic events used with respect to station TANA.

3.3. Joint inversion of surface wave and receiver function data Jointly inverting multiple sets of data simultaneously generates a velocity model that is more holistic since different geophysical studies are sensitive to different parameters of the earth’s structure. For each of the given two sets of data, surface wave dispersion (for Rayleigh wave group velocities) and receiver function, a synthetic is prepared corresponding to the initial earth structure. The synthetic is then compared to the observed and is rejected if errors exceed the maximum threshold. The goal is to iteratively achieve a minimum error for the functional(S), where

!2  NS Nr  OSj  PSj Ori  Pri 2 ð1  pÞ X p X S¼ þ Nr i¼0 N S j¼0 rr i rSj

ð1Þ

where Ori , observed receiver function at time ti; Pri , predicted receiver function at time ti; rri , standard error of observation at time ti; OSj , j0 th observed surface wave dispersion; PSj , j0 th predicted surface wave dispersion; rSj , standard error of j0 th surface wave observation; Nr , total number of receiver function points; NS , total number of surface wave dispersion points; p, influence factor, 0 P p P 1 (after Julia et al., 2000).

The parameter ’p’ can be altered based on the quality of data available. It alters the influence either set of data will have on the error minimization process. A scenario where p = 1 would generate a solution based entirely on the surface wave dispersion data. Conversely, p = 0 would cause the solution to be exclusively dependent on the receiver function data. Since the quality of both sets of data i.e., surface wave dispersions as well as the receiver functions, have a high SNR and provide information about the velocity and depth structures respectively, the inversion for velocities has a p factor solely dependent on the surface wave dispersion data and the depth inversion routine is entirely dependent on the RF data. The layer weightage is alterable, and to provide a better emphasis on the crustal depth, the weights of the layers binding the Moho are greater (0.8) than those for the upper crustal layers (0.5). The standard errors of data can be estimated if surface wave observations used are averages of group velocity measurements for several events at a station, and, receiver functions from a range of slownesses are stacked. However, the number of events utilized in the study are too few to give statistically reliable estimates of the standard errors. Hence, we preferred to assume the standard errors of 0.05 km/s for the dispersion curves and 0.02 s1 for the receiver functions, which are characteristic error bounds for group velocity measurements in dispersion and receiver function studies respectively (Julia et al., 2000). Normalizing Eq. (1) by the standard

58

A.A. Deshpande, G. Mohan / Journal of Asian Earth Sciences 128 (2016) 54–63

Fig. 4. The P wave receiver functions plotted with respect to delay time (s) for varying ray parameters at each station. The P to s conversion from the Moho is denoted as Pms.

errors corrects for the different physical units of the dispersion (km/s) and receiver function (s1), and avoids dominance of one data set over the other. A linearized, iterative algorithm (Julia et al., 2000) is employed to invert for a one-dimensional layered model by minimizing a weighted least squares norm consisting of surface wave and receiver function data sets, a model roughness norm, and a matrix of differences between inverted and preset model parameters. In this study, instead of generating a single stacked RF waveform, the inversion routine attempts to fit multiple RF waveforms corresponding to different slownesses simultaneously. Similarly, with every iteration, a ’better’ fit curve among the dispersion curves for the given station is generated. The initial model used is a hybrid derived from DSS studies which include both sub-basaltic sediments (Kaila et al., 1980) and a high velocity lower crust (Rao and Tewari, 2005). Higher weightage is applied to the lower crustal and uppermost mantle layers to emphasize the Moho characteristics. The Vp/Vs for each layer is kept fixed based on the Poisson’s ratio of 0.27 obtained from earlier studies in Saurashtra (Praveen Kumar and Mohan, 2014). Inversion was also attempted for fixed Poisson’s ratio of 0.26 and 0.28 to eliminate any bias. The joint inversion yielded the shear wave velocity structure beneath each source-receiver path corresponding to the best fit curves observed in Fig. 5 for surface wave dispersion and RF waveforms respectively. 4. Shear velocity structure The shear velocity models obtained from inversion of solely surface wave dispersion curves, together with the velocity-depth models obtained by jointly inverting the RFs at individual stations

with the corresponding dispersion curves, are shown in Fig. 5. A good fit is observed for the dispersion curves and the RFs at all stations for a four layered crustal shear wave velocity model overlying the upper mantle. Although the velocity-depth models estimated by both methods appear similar, subtle variations are observed in the thicknesses and shear velocities of intracrustal layers respectively (Fig. 5). It is observed that the first layer in the model corresponding to the basaltic and sub-basaltic layers is not well constrained by surface wave inversion due to the absence of reliable short period data with time periods less than 5 s. The thickness of the top layer is estimated to range from 4.26 km to 6.49 km which is overestimated, since the lithologs from the well at Lodhikha in Saurashtra (Fig. 1) clearly indicate that the depth to the basement averages 3.9 km (Sain et al., 2002). In contrast, the joint inversion estimates the thickness of the top layer to be averaging 3.0 ± 1.2 km (Table 1) which places a tighter constraint, primarily due modeling the strong Ps conversion from the base of the top layer observed in the RFs. The average shear velocities of the top layer are estimated to be 2.98 km/s while joint inversion reveals a Vs of 2.86 km/s which is more appropriate as lithologs reveal the presence of sub-basaltic sediments which should contribute to lowering the velocities of the top layer. Surface wave inversion estimates the depth to the Moho along the paths to MULG and TANA to be  39.5 km and the lower crustal thickness to be  18 km, both of which are larger and inconsistent with the known DSS model in Saurashtra (Rao and Tewari, 2005) and known crustal thicknesses (Rao et al., 2015). Similarly, the depth to Moho along the path to BHOD is estimated to be 34 km which is underestimated. Summary of the inversion results yields a larger range in the thicknesses of the lower layer for solely surface wave

A.A. Deshpande, G. Mohan / Journal of Asian Earth Sciences 128 (2016) 54–63

59

Fig. 5. The joint inversion of Rayleigh wave dispersion curves and P wave receiver function data at all stations. The blue and red curves are the observed and synthetic best fits respectively. The velocity depth section indicates the shear wave velocity structure obtained through joint inversion (red) and only surface wave inversion (blue). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

60

A.A. Deshpande, G. Mohan / Journal of Asian Earth Sciences 128 (2016) 54–63

Table 1 The mean and standard deviation of the layer thickness and shear velocity for the velocity-depth models obtained for the short paths to stations in Saurashtra, the long paths to stations across the Gulf of Cambay, and the average for NWDVP. The crustal averages are given in bold. Model layers

Layer 1 Layer 2 Layer 3 Layer 4 Crustal average Upper mantle

Short paths

Long paths

Average

H (km)

Vs (km/s)

H (km)

Vs (km/s)

H (km)

Vs (km/s)

2.70 ± 1.60 6.96 ± 2.53 14.14 ± 2.16 14.31 ± 2.15 38.12 ± 1.22 Half space

2.85 ± 0.25 3.27 ± 0.10 3.67 ± 0.14 3.95 ± 0.23 3.66 ± 0.12 4.43 ± 0.08

3.52 ± 0.36 5.49 ± 1.01 12.15 ± 0.68 14.59 ± 2.12 35.74 ± 2.02 Half space

2.77 ± 0.07 3.15 ± 0.15 3.57 ± 0.06 3.94 ± 0.16 3.57 ± 0.03 4.28 ± 0.07

3.05 ± 1.23 6.33 ± 2.04 13.29 ± 1.90 14.43 ± 1.96 37.10 ± 1.93 Half space

2.86 ± 0.19 3.24 ± 0.13 3.64 ± 0.12 3.94 ± 0.19 3.62 ± 0.10 4.36 ± 0.10

inversion relative to joint inversion where they appear to be better constrained within a smaller range. For instance, the lower crustal thicknesses are approximately in the range 11 km–19 km averaging 15.92 ± 2.81 km compared to 12–16 km, averaging 14.43 ± 1.96 km estimated from joint inversion. The range and standard deviations for the thickness of the lower crust are smaller for the models estimated through joint inversion. The velocitydepth models estimated through inversion of surface waves are subject to depth trade off which can be constrained by using RFs. However, the average crustal velocities and thicknesses obtained by both methods for the entire study area of NWDVP are remarkably similar to each other. The joint inversion provides tighter constraints on the velocity-depth models which are summarized in Table 1 alongwith the average models for the short and long paths which cut across Saurashtra and the Gulf of Cambay respectively. The top layer corresponding to the basaltic layer/basalts + subbasaltic sediments reveals significant diversity across the study region with the Vs ranging from 2.5 km/s to 3.06 km/s and averaging 2.86 ± 0.19 km/s. The shear velocities for the three crustal layers beneath the top basaltic layer range from 3.04 km/s to 3.38 km/ s, 3.47 km/s to 3.79 km/s and 3.7 km/s to 4.19 km/s (Fig. 5) and average 3.24 ± 0.13 km/s, 3.64 ± 0.12 km/s and 3.94 ± 0.19 km/s respectively. The shear velocities of the four crustal layers are clearly distinct from one another along all the source-receiver paths (Fig. 5). The lower crustal Vs are mostly in the range 3.82– 4.19 km/s along all the source receiver paths with the exception of the path to KAND where Vs  3.7 km/s is significantly lower. The paths to BODH, SONT and VARE exhibit the highest lower crustal Vs > 4.1 km/s. The crustal thickness in the region along all paths averages 37 ± 1 km. The average crustal Vs along the sourcereceiver paths in Saurashtra ranges from 3.49 km/s to 3.77 km/s, averaging 3.66 ± 0.1 km/s. The paths cutting across the Gulf of Cambay to MULG, KHER and VARE exhibit lower average crustal Vs ranging from 3.53 km/s to 3.58 km/s relative to those observed in Saurashtra. The shallow upper mantle shear velocity in the study region ranges from 4.2 km/s to 4.5 km/s, averaging 4.36 ± 0.1 km/s. The paths to SONT, TANA, KAND and BODH in Saurashtra exhibit upper mantle Vs ranging from 4.35 km/s to 4.5 km/s. The lowest Vs of 4.2 km/s is observed for the paths cutting across the Gulf of Cambay to KHER, followed by 4.35 km/s for paths to MULG and VARE which are lower than the 4.6 km/s reported by Bhattacharya (1981, 1992) and 4.58 km/s–4.78 km/s reported by Lyon-Cean (1986) for the Indian shield region.

5. Discussion The significant finding revealed by the velocity depth profiles detailed earlier is the wide range of average crustal shear wave velocities along the azimuthally diverse paths in NWDVP, the low crustal shear velocities along paths sampling the Gulf of Cambay, and, the significantly low shear velocities in the uppermost upper mantle. A comparison of the shear velocity – depth structure

estimated for the short paths in Saurashtra and long paths cutting across the Gulf of Cambay respectively (Table 1) indicates consistently lower shear velocities for all the layers for the long paths. The average crustal Vs in NWDVP is observed to be in three ranges i.e., 3.48–3.57 km/s, 3.58–3.67 km/s and 3.68–3.77 km/s respectively (Fig. 6). The Vs varies from a low in central Saurashtra to a high in the eastern and western margins of Saurashtra. The average crustal Vs for the paths across the Gulf of Cambay ranges from 3.53 km/s to 3.58 km/s, exhibiting a reduction in velocities. The estimated seismic heterogeneities complement the density heterogeneities reported for Saurashtra through gravity studies (Mishra et al., 2001; Tewari et al., 2009), which indicate a low in central Saurashtra and Cambay rift region, and Bouguer highs along the margins of Saurashtra. The density heterogeneities are attributed to the variation in the thickness of the sub-basaltic sediments and the presence of volcanic plugs with shallow roots (Tewari et al., 2009). The gravity anomalies primarily represent the shallow crustal density heterogeneities (Tewari et al., 2009) while the increase in average crustal Vs in eastern and western Saurashtra along paths which cut across gravity lows is due to the deeper mafic underplating (Vs > 4.0). The estimated low Vs in central Saurashtra is an average cumulative effect of the shallow subbasaltic sediments coupled with the relatively low Vs in the lower crust. The lithologs from the boreholes at Lodhika and Dhanduka (Fig. 1) indicate the presence of sub-basaltic sediments of varying thicknesses, which contribute to the variation in the shallow upper crustal shear velocities. The rapid and large volume of basalts flooding Saurashtra through the adjoining rift zones has led to several crustal intrusions in terms of dykes and sills as the molten magma rises to the surface through fissures. A variety of compositions would be produced during the ascent of the magma through the continental crust by various differentiation processes such as fractional crystallization, crustal contamination and addition of magmas derived from continental lithosphere (e.g., Sano et al., 2001; Lightfoot et al., 1990; Hooper, 1994). A valuable direct information about subsurface composition is provided by the presence of several volcano-plutonic complexes, mafic dyke swarms containing picrites, sub-alkali basalts and basaltic andesites, and large variation in the rock types including picrite, lamprophyres, diorite and rhyolite (Cucciniello et al., 2015; Melluso et al., 1995; Sheth et al., 2011, 2012). The extensive crustal heterogeneity is evidenced by the variation in the shear velocity structure which also complements and compares well with the variation in the bulk Vp/ Vs ratio  1.8 ± 0.11 across NWDVP (Rao et al., 2015). These Vp/Vs ratios are used as a proxy to constrain the composition of the lower crust under the assumption that a mafic lower crust will enhance the Vp/Vs ratios. The bulk Vp/Vs estimated using the average crustal Vs  3.66 km/s in Saurashtra along with the average P wave velocity (Vp)  6.54 km/s derived independently through deep seismic studies yield a bulk Vp/Vs ratio  1.786 or the equivalent Poisson’s ratio of 0.271. Significant contribution to the bulk Poisson’s ratio comes from the 12 to 16 km thick lower crustal layers with Vs ranging from 3.83 km/s to 4.19 km/s and averaging

A.A. Deshpande, G. Mohan / Journal of Asian Earth Sciences 128 (2016) 54–63

61

Fig. 6. Variation of average shear velocities along paths in NWDVP overlying the Bouguer gravity anomaly map contoured at intervals of 5 mGals (after Tewari et al., 2009).

3.95 ± 0.18 km/s, with the exception of the path to KAND. The high Vs in the lower crust is complimentary to the high average Vp  7.2 km/s obtained from DSS results (Rao and Tewari, 2005) corresponding to the lithologies of an underplated crust. As an average Vp is available from DSS, the Vp/Vs ratios are estimated using the average Vs estimated for the lower crust, which serves as a link to lithological composition. An average Vs  3.95 km/s together with an average Vp of 7.2 km/s for the lower crust results in a Vp/Vs ratio of 1.82 and an equivalent Poisson’s ratio of 0.286. Laboratory measurements of seismic velocities in rock samples (Christensen, 1996) report that for depths between 35 km and 50 km, a Vp/Vs ratio of 1.792 ± 0.130 is observed for felsic granulites, 1.818 ± 0.076 for mafic granulites and 1.807 ± 0.092 for mafic garnet. This suggests that shear velocities up to 4.0 km/s can be accounted for by any of the lithologies and shear velocities ranging between 4.1 and 4.3 km/s can only be accommodated by rocks of mafic composition (Christensen, 1996; Julia et al., 2009). These measurements assume average heat flow gradients. However, the crustal geotherms play a vital role in the reduction of shear velocities and the heat flow values reported for the Cambay rift are high (Roy and Rao, 2000). The high lower crustal shear velocities can be interpreted to be due to emplacement of mafic

cumulates in the lower crust. Correlation of the spatial extent of the mafic cumulates in the deep crust with the pervasive mafic dyke swarms observed on the surface, leads to the conclusion that mafic underplating occurred during the Deccan volcanism about 65 Ma ago. However, a similar joint inversion study of the Indian shield by Julia et al. (2009) indicate high lower crustal velocities in southcentral DVP, as well as in areas which have not been affected by Deccan volcanism. A plausible explanation for such high lower crustal velocities was given by French et al. (2008), who suggested that it is possible that mafic dyke swarms may be buried beneath the thick Deccan basaltic flows and could represent erosional remnants of feeders to continental flood basalts from a Proterozoic event. Possible evidence to such hypothesis also is seen in the borehole at Lodhika wherein basalts of uncertain ages are reported beneath the sub-basaltic Mesozoic sediments (Singh et al., 1997). The average shear velocity – depth model obtained through joint inversion is compared with the shear velocity models (Fig. 7) estimated for the Indian shield (Bhattacharya, 1981), southcentral DVP (Julia et al., 2009) and NWDVP which encompasses parts of the western segment of SCDVP (Prajapati et al., 2011) respectively. As illustrated in Fig. 7, the lower crustal velocities

62

A.A. Deshpande, G. Mohan / Journal of Asian Earth Sciences 128 (2016) 54–63

and Gudfinnsson (2005) suggested that the presence of carbonatite reduces the solidus temperature, resulting in melting and development of low velocity zones. Interestingly, geodetic measurements around Bhuj in the Kutch rift revealed a weak mantle with unusually low viscous strength which is possibly a long lasting result of thermal weakening by the source of the Deccan volcanism (Chandrasekhar et al., 2009). Based on these post-seismic deformation studies, it was proposed that the mantle lithosphere in NW India is weaker than that of a stable continental shield, which is in concurrence with the present results. Therefore, a plausible explanation for the significant crustal alterations in NWDVP is the possibility of weakening of the lithosphere by a mantle thermal source leading to a thin and weak lithosphere which facilitated the Deccan volcanism through the rift zones of northwestern India. 6. Conclusions

Fig. 7. Comparison of average crustal shear wave velocity-depth model estimated from the present study with those estimated earlier for the Indian shield, southcentral DVP (SCDVP), and northwestern DVP (NWDVP) which encompasses parts of the western segment of SCDVP.

are higher, indicating magmatic underplating, and the upper mantle velocities are anomalously low. The shallow upper mantle shear velocities along the source-receiver paths across the Gulf of Cambay and Saurashtra range from 4.2 to 4.5 km/s, averaging 4.36 km/s, which is lower than those reported for the Indian shield and southcentral DVP. The Gulf of Cambay is a zone of confluence of the north-south trending Cambay rift and the east-west trending Narmada rift. The Cambay rift is known to be a zone of high heat flow of about 80 mW which is much higher than the average of 44 mW observed in DVP (Roy and Rao, 2000). The low shear velocities in the uppermost upper mantle are an indication that the residual thermal anomaly plays a dominant role in the reduction of velocities which is in consonance with the P-wave perturbations observed in northwestern India, which extend south into Saurashtra (Kennett and Widiyantoro, 1999). Residual thermal anomalies associated with the 65 Ma old Deccan volcanism are expected to be preserved as per the estimates of Ziegler and Cloetingh (2004) who propose that
35% of a deep seated thermal anomaly would be present even after 60 Ma, while numerical modeling predicts that thermal anomalies can last up to 200 Ma (Campbell, 2007). The findings in the present study also substantiate the indirect inference of low velocities in the upper mantle based on the anomalous 1.0 s delay time of P410s converted phases observed at stations in Saurashtra (Mohan et al., 2012). The overall crustal and uppermost upper mantle shear wave velocity structure reveals that the lithosphere beneath NWDVP is altered due to thermal and compositional changes associated with the Deccan volcanism. Compositional changes may be associated with the variations in the degree of depletion due to magmatic activity and metasomatism (Thybo, 2006). Based on geochemical evidences, Sen et al. (2009) hypothesized that the lithosphere beneath Kutch in NWDVP was subjected to carbonatite metasomatism. Presnall

The average crustal shear wave velocities were estimated to range from 3.5 km/s to 3.8 km/s indicating the presence of significant crustal heterogeneities. Seismic evidence of widespread magmatic underplating is seen in terms of high lower crustal velocities ranging from 3.8 km/s to 4.19 km/s and high lower crustal Poisson’s ratio  0.286. The crustal and uppermost mantle shear velocities for the paths cutting across the Gulf of Cambay are lower than those estimated for the short paths in Saurashtra. The shear wave velocities in the uppermost mantle are anomalously low, averaging 4.36 km/s, which is lower than that estimated for the Indian shield. The low shear velocities are a result of the combined effects of a residual thermal anomaly and compositional changes associated with the Deccan volcanism. The study reveals that the crustal and upper mantle structure beneath NWDVP is significantly altered by the Deccan volcanism. Acknowledgements This work was sponsored by the Deep Continental Studies Program of the Department of Science and Technology, India. The authors thank the anonymous reviewers for their suggestions in improving the manuscript. References Ammon, C.J., Randall, G.E., Zandt, G., 1990. On the non-uniqueness of receiver function inversions. J. Geophys. Res. 95, 15303–15318. Bhattacharya, S.N., 1981. Observation and inversion of surface wave group velocities across central India. Bull. Seismol. Soc. Am. 71, 1489–1501. Bhattacharya, S.N., 1992. Crustal and upper mantle structure of India from surface wave dispersion. Curr. Sci. 62, 94–100. Bodin, P., Horton, S., 2004. Source parameters and tectonic implications of aftershocks of the Mw 7.6 Bhuj earthquake of 26 January 2001. Bull. Seismol. Soc. Am. 94 (3), 818–827. Campbell, I.H., 2007. Testing the plume theory. Chem. Geol. 241, 153–176. http:// dx.doi.org/10.1016/ j.chemgeo.2007.01.024. Chandrasekhar, D.V., Mishra, D.C., Rao, G.V.S.P., Rao, J.M., 2002. Gravity and magnetic signatures of volcanic plugs related to Deccan volcanism in Saurashtra, India and their physical and geochemical properties. Earth Planet. Sci. Lett. 201, 277–292. Chandrasekhar, D.V., Burgmann, R., Reddy, C.D., Sunil, P.S., Schmidt, D.A., 2009. Weak mantle in NW India probed by geodetic measurements following the 2001 Bhuj earthquake. Earth Planet. Sci. Lett. 280, 229–235. Christensen, N.I., 1996. Poisson’s ratio and crustal seismology. J. Geophys. Res. 101, 3139–3156. http://dx.doi.org/10.1029/95JB03446. Cucciniello, C., Demonterova, E.I., Sheth, H., Pande, K., Vijayan, A., 2015. 40Ar/39Ar geochronology and geochemistry of the central Saurashtra mafic dyke swarms: insights into magmatic evolution, magma transport, and dyke-flow relationships in the northwestern Deccan traps. Bull. Volcanol. 77, 45. Dziewonski, A., Bloch, S., Landisman, M., 1969. A technique for the analysis of transient seismic signals. Bull. Seismol. Soc. Am. 59 (1), 427–444. French, J.E., Heaman, L.M., Chacko, T., Srivastava, R.K., 2008. 1891–1883 Ma Southern Bastar-Cuddapah mafic igneous events, India: a newly recognized large igneous province. Precambr. Res. 160, 308–322. http://dx.doi.org/10.1016/ j.precamres.2007.08.005.

A.A. Deshpande, G. Mohan / Journal of Asian Earth Sciences 128 (2016) 54–63 Herrmann, R.B., Ammon, C.J., 2002. Surface waves, receiver functions and the crustal structure in Computer programs in Seismology, version 3.30., 110pp, . Hooper, P.R., 1994. Sources of continental flood basalts: the lithospheric component. In: Subbarao, K.V. (Ed.), Volcanism. Wiley Eastern, New Delhi, pp. 29–53. Julia, J., Ammon, C.J., Herrmann, R.B., Correig, A.M., 2000. Joint inversion of receiver function and surface wave dispersion observations. Geophys. J. Int. 143, 1–19. Julia, J., Jagadeesh, S., Rai, S.S., Owens, T.J., 2009. Deep structure of the Indian shield from joint inversion of P wave receiver functions and Rayleigh wave group velocities implications for precambrian crustal evolution. J. Geophys. Res. 114, B10313. http://dx.doi.org/10.1029/2008JB006261. Kaila, K.L., Krishna, V.G., Tewari, H.C., Sharma, P.L.N., 1980. Crustal structure from deep seismic sounding studies along Navibandar-Amreli profile in Saurashtra, India. Memoirs Geol. Soc. India 3, 218–232. Kaila, K.L., Krishna, V.G., Mall, D.M., 1981. Crustal structure along the MehmadabadBillimora profiles in the Cambay Basin, India, from deep seismic soundings. Tectonophysics 76, 99–130. Kennett, B.L.N., Widiyantoro, S., 1999. A low seismic wavespeed anomaly beneath northwestern India: a seismic signature of the Deccan plume? Earth Planet. Sci. Lett. 165, 145–155. Langston, C.A., 1979. Structure under Mount Rainier, Washington, inferred from teleseismic body waves. J. Geophys. Res. 84, 4749–4762. Ligorria, J.O., Ammon, C.J., 1999. Iterative deconvolution of teleseismic seismograms and receiver function estimation. Bull. Seismol. Soc. Am. 89, 1395–1400. Lightfoot, P.C., Hawkwsworth, C.J., Devey, C.W., Rogers, N.W., Van Calsteren, P.W.C., 1990. Source and differentiation of Deccan Trap lavas: implications of geochemical and mineral chemical variation. J. Petrol. 31, 1165–1200. Lyon-Cean, H., 1986. Comparison of the upper mantle shear wave velocity structure of the Indian shield and the Tibetan plateau and tectonic implications. Geophys. J. Roy. Astron. Soc. 86 (3), 727–749. http://dx.doi.org/10.1111/j.1365246X.1986.tb00657.x. Mandal, P., Horton, S., 2007. Relocation of aftershocks, focal mechanisms and stress inversion: Implications toward the seismo-tectonics of the causative fault zone of Mw7.6 2001 Bhuj earthquake (India). Tectonophysics 429, 67–78. Melluso, L., Beccaluva, L., Brotzu, P., Gregnaniun, A., Gupta, A.K., Morbidelli, I., Traversa, G., 1995. Constraints on the mantle sources of the Deccan Traps from the petrology and geochemistry of the basalts of Gujarat state (western India). J. Petrol. 36, 1393–1432. Mishra, D.C., Singh, Bijendra, Gupta, S.B., Prabhakar Rao, M.R.K., Singh, A.P., et al., 2001. Major lineaments and gravity magnetic trends in Saurashtra, India. Curr. Sci. 80, 1059–1067. Mohan, G., Ravi Kumar, M., Saikia, D., Praveen Kumar, K.A., Tiwari, Pankaj Kr., Surve, G., 2012. Imprints of volcanism in the upper mantle beneath NW Deccan Volcanic Province. Lithosphere. http://dx.doi.org/10.1130/L178.1. Owens, T.J., 1984. Determination of crustal and upper mantle structure from analysis of broadband teleseismic P-waveforms Ph. D. thesis. Department of Geology and Geophysics, The University of Utah. Ozalaybey, S., Savage, M.K., Sheehan, A.F., Louie, J.N., Brune, J.N., 1997. Shear wave velocity structure in the northern Basin and Range Province from the combined

63

analysis of receiver functions and surface waves. Bull. Seismol. Soc. Am. 87 (1), 183–199. Praveen Kumar, K.A., Mohan, G., 2014. Crustal velocity structure beneath Saurashtra, NW India, through waveform modeling: Implications for magmatic underplating. J. Asian Earth Sci. 79, 173–181. Prajapati, S., Suresh, G., Bhattacharya, S.N., 2011. Crustal structure of the northwest Deccan Volcanic Province, India, and the adjoining continental shelf from observed surface-wave dispersion. Bull. Seismol. Soc. Am. 101 (4), 1488–1495. Presnall, D.C., Gudfinnsson, G.H., 2005. Carbonate rich melts in the oceanic low velocity zone and deep mantle. In: Foulger, G.H., Natland, J., Presnall, D.C., Anderson, D.L., (Eds.), Plates, Plumes and Paradigms: Geological Society of America Special Paper, vol. 388, pp. 207–216. Rao, K., Ravi Kumar, M., Rastogi, B.K., 2015. Crust beneath the northwestern Deccan Volcanic Province, India: evidence for uplift and magmatic underplating. J. Geophys. Res. 120, 3385–3405. http://dx.doi.org/10.1002/2014/B011819. Rao, G.S.P., Tewari, H.C., 2005. The seismic structure of the Saurashtra crust in northwest India and its relationship with the Reunion plume. Geophys. J. Int. 160, 318–330. Roy, S., Rao, R.U.M., 2000. Heat flow in the Indian shield. J. Geophys. Res. 105, 25587–25604. Sain, K., Zelt, C.A., Reddy, P.R., 2002. Imaging of subvolcanic mesozoics in the Saurashtra peninsula of India using travel time inversion of wide-angle seismic data. Geophys. J. Int. 150, 820–826. Sano, T., Fujii, T., Deshmukh, S.S., Fukuoka, T., Aramaki, S., 2001. Differentiation processes of Deccan trap basalts: contribution from geochemistry and experimental petrology. J. Petrol. 42 (12), 2175–2195. Sen, G., Bizimis, M., Das, R., Paul, D.K., Ray, A., Biswas, S., 2009. Deccan plume, lithosphere rifting, and volcanism in Kutch, India. Earth Planet. Sci. Lett. http:// dx.doi.org/10.1016/j.epsl.2008.10.002. Sheth, H.C., Choudhary, A.K., Bhattacharya, S., Cucciniello, C., Laishram, R., Gurav, T., 2011. The chogat-Chamardi subvolcanic complex, Saurashtra, northwestern Deccan traps: geology, petrochemistry, and petrogenetic evolution. J. Asian Earth Sci. 41, 307–324. Sheth, H.C., Cucciniello Choudhary, C.A.K., Bhattacharya, S., Laishram, R., Gurav, T., 2012. Geology, petrochemistry and genesis of the bimodal lavas of Osham hill, Saurashtra, northwestern Deccan traps. J. Asian Earth Sci. 43, 176–192. Singh, D., Alat, C.A., Singh, R.N., Gupta, V.P., 1997. Source rock characteristics and hydrocarbon generating potential of Mesozoic sediments in Lodhika area, Saurashtra basin, Gujarat, India. In: Proceedings of 2nd International Petroleum Conference and Exhibition. Petrotech-97, New Delhi, pp. 205–220. Tewari, H.C., Rao, G.S.P., Rajendra Prasad, B., 2009. Uplifted crust in parts of western India. J. Geol. Soc. India 73, 479–488. Thybo, H., 2006. The heterogeneous upper mantle low velocity zone. Tectonophysics 416, 53–79. http://dx.doi.org/10.1016/j.tecto.2005.11.021. USGS. United States Geological Survey, . White, R.S., McKenzie, D.P., 1989. Magmatism at rift zones: the generation of volcanic continental margins and flood basalts. J. Geophys. Res. 94, 7685–7729. Ziegler, P.A., Cloetingh, S., 2004. Dynamic processes controlling evolution of rifted basins. Earth-Sci. Rev. 64, 1–50. http://dx.doi.org/10.1016/S0012-8252(03) 00041-2.