P-Wave Teleseismic Tomography: Evidence of Imprints of Deccan Mantle Plume below the Kachchh Rift Zone, Gujarat, India

Abstract

The Indian plate had experienced the Deccan volcanism at 65 Ma when it moved over the Re-union hotspot, which has altered lithospheric structure below the Kachchh rift zone (KRZ). To quantify the influence of Deccan volcanism on the crust-mantle, the present chapter focuses on the delineation of the upper mantle structure below the KRZ, through the modeling of crust corrected P-residuals and P-wave teleseismic tomography. The crust corrected normalized P-residuals suggest dominant negative residuals associated with the central KRZ, indicating crustal and lithospheric thinning below the KRZ. A low velocity down to a depth of 170 km below the central KRZ is detected through the teleseismic tomography using these P-residuals. However, these residuals also show positive values for the surrounding un-rifted zones. Note that a low shear velocity zone extending from 100–120 km to 170–220 km depth beneath the central KRZ has already been revealed by the modeling of P-RFs. This reduction in seismic velocity in the upper mantle could be explained by the presence of trapped carbonatite/partial melts related to the Deccan volcanism. The influx of volatile CO2 emanating from the carbonatite melts in the asthenosphere might be generating lower crustal earthquakes occurring in the KRZ.

1. Introduction

At ~860 Ma, Indian plate was broken from Antarctica and then it moved along with African plate until 184 Ma [1, 2]. Subsequently, it moved toward North until it got separated from Madagascar at 88 Ma [3]. At 65 Ma, it moved over the Reunion hotspot, which led to outpouring of Deccan volcanism resulting in the Deccan volcanic province occupying an area of 0.5 million km² with 2000 m thick basaltic lava in the western and central parts of the Indian subcontinent within a short duration of time [4].

Occurrence of the initial outpouring of Deccan volcanism in Pakistan (at ~72–73 Ma) has been evidenced by the available ⁴⁰Ar-³⁹Ar age data of the Reunion Island like alkali lavas [3]. The presence of relatively younger basaltic intrusions (~68.5 Ma) in the northern India [5] has supported the fact that India moved toward south. Also, younger basalts have been found further south [6]. Finally, at 65 Ma, the 95% basaltic outpouring of the Deccan volcanism took place in the western and central India [7]. This Deccan plume model has got further support from the presence of a low-velocity anomaly in the upper mantle extending up to a depth of 600 km below the north-western India as modeled by the regional mantle P-wave seismic tomography [8]. The presence of high ³He/⁴He ratio in Rajasthan has provided further support for the Deccan plume model [5]. Adding credence to this plume theory, local earthquake tomography, surface wave dispersion and receiver function modeling studies have imaged crustal mafic underplating, Moho upwarping and asthenospheric thinning underlying the 2001 Bhuj earthquake region [9, 10]. Recently, the shear-wave splitting study suggested that the upper mantle anisotropy in the KRZ is contributed both by lithospheric frozen anisotropy and asthenospheric flow induced anisotropy, which could be inherited from the plume-lithosphere interaction during the Deccan/Reunion plume episode (~65 Ma) [11]. These observations suggest that imprints of the Deccan/Reunion mantle plume are still present in the crust and upper mantle below the north-western region.

The Kachchh region has been experiencing earthquakes since historical times [12]. The region has already experienced seven M6 earthquakes including two M_w 7.7 earthquakes in 1989 and 2001. The latter event has claimed a death toll of 20,000 people. The aftershock activity of this 2001 earthquake is continuing until today, with regular occurrences of M_w 3 events and occasional occurrences of M_w 4 events. The aftershock activity of the 2001 Bhuj earthquake is still continuing that includes 15 M_w ≥ 5, about 300 M_w ≥ 4 and about 6000 M_w ≥ 3 events. We feel that the enigmatic seismicity associated with the Kachchh rift zone is linked with its above-discussed unique Geodynamic history. Aiming at understanding the influence of crustal-mantle structure in the genesis of uninterrupted occurrences of earthquakes since the occurrence of the 26 January M_w 7.7 Bhuj earthquake, in this chapter, the crust corrected P-wave residuals are estimated at 14 broadband stations and a 3-dimensional P-wave teleseismic tomography is performed using the estimated crust corrected residuals. Finally, modeling results are interpreted concerning the geodynamical processes responsible for generating intraplate earthquakes occurring in the Kachchh region.

2. Seismic network and data

For the present study, broadband digital waveforms of 241 teleseismic events from 14 three-component broadband stations in Kachchh, Gujarat are used (Figure 1a and b). The sampling rate of recording is 50 sps. Station spacing in the above seismic network is 30–100 km, however, for our study we use data from 20 broadband stations consisting of 14 NGRI and 6 ISR stations as shown by a square in Figure 2. The station spacing for the network consisting of these 20 stations is 30–60 km. We selected 241 teleseismic events with magnitude ranging from 5.9 to 8.2, with epicentral distances of 29–90° (Figure 1b) with reference to the center of the network.

Figure 1.

(a) Elevation (in m) map showing station distribution in Kachchh, Gujarat. Filled red triangles mark the broadband seismograph stations (SIV—Sivlok; VJP—Vajepar; KNM—Kanmer; BCH—Bhachau; BEL—Bela; GDD—Gadhada; JUM—Jumkunaria; MTP—Motapaya; NGR—Nagor; BHU—Bhuj; NPR—Narayanpar; TPM—Tapar Mundra; MND—Mandvi; VGH—Vaghura; TPR—Tapar Anjar; and NDD—New Dudhai). The inset shows the key map for the area, where the study area is shown by a red open square. K and KL mark Kachchh, Gujarat and Killari, Maharashtra, respectively. An arrow shows the location of Cambay. The black filled portion marks the areal extent of Deccan volcanic province (DVP) in India. (b) Epicentral plot of 241 teleseismic events of M_w 6.0–8.4, whose broadband data are used for our P-wave teleseismic tomography study. A red triangle and green diamond symbols mark the center of our network (lat. 70°, long. 23°) and epicenters of selected teleseismic events. The size of the diamond symbols vary depending on their sizes.

Figure 2.

(a) Crust corrected relative P-residuals in Gujarat. (b) Directional mean of P-residuals by subtracting an average residual from the eight stations in Kachchh. Red open circles mark negative residuals while blue open circles show positive residuals. CRZ marks the Cambay rift zone, which is shown by dotted lines. Black dotted elliptical area mark the central Kachchh rift zone characterized by negative residuals. And, black square area shows our study area considered for our P-wave teleseismic tomography.

First, each trace of an event is used to pick the arrival time of the first P-wave maximum amplitude (either peak or trough). Then, these picks are correlated within the network. Following this procedure, the first P-wave onset times are picked from selected highest-quality traces recorded at different stations for one event. Here, the uncertainty of the picking is used to decide a quality factor for each measurement. Note that most of the measurements are found to be having an uncertainty of ±0.05 s. Finally, these quality factors are used to estimate the average data error, which is found to be ±0.06 s.

3. Calculation of crust-corrected P-relative residual

For picking P, the measured waveforms from our network are bandpass filtered using the inbuilt WWSSN filter of the seismic handler. This narrow frequency band allows us to study the uppermost mantle without caring about finite frequency effects [13]. In order to calculate relative residuals, only first strong and clear phase (peak or trough) are picked from different stations, with the help of high magnification on the computer screen. Since the sampling rate used for data recording is 50 sps, thus, phases could be picked manually with a precision of at least 0.08 s. Thus, the most of the estimated relative residuals are found to have an uncertainty of less than 0.1 s.

The analysis of broadband waveforms of teleseismic events are used here to estimate travel time residuals, as well as deviations in slowness and back azimuth about the standard earth model iasp91 [14]. A travel time residual can be defined as

$${\mathrm{r}}_{\mathrm{ij}}={{\mathrm{t}}_{\mathrm{ij}}}^{\mathrm{m}}-{{\mathrm{t}}_{\mathrm{ij}}}^{\mathrm{t}}$$E1

where t_ij^m is the measured arrival time of a teleseismic phase j from an event at ith station and t_ij^t is the theoretical arrival time derived from an earth model (e.g., iasp91). Now, a systematic bias to r_ij could be resulted from following sources:

1. the uncertainties in origin time and hypocenter location of the teleseismic events

2. the source side heterogeneity.

To isolate information underneath the receiver network, the average of all measurements (n) of a phase j is subtracted from the estimated relative residual using the following equation,

$${{\mathrm{r}}_{\mathrm{ij}}}^{\mathrm{r}}={\mathrm{r}}_{\mathrm{ij}}-\left(\mathrm{1/n}\right)\underset{i=1}{\overset{n}{\Sigma }}\mathit{rij}$$E2

Thus, the residual part, which is common in all measurements, is eliminated. Therefore, the remaining part of traveltime residual, which is caused by heterogeneity below the receivers, defines the relative residual. This normalization (Eq. (2)) produces an equal amount of positive and negative residual times per phase with a zero mean. Now, to use these relative residuals in tomographic inversions, a weighting scheme following Evans and Achauer [15] has been used here. The weights w_ij are assigned while picking the arrival times based on arrival time uncertainties (d) and SNR of the picks. We use four different weights (1, 2, 3, and 4): 1 = 1.0 for d ≤ 0.05 s, 2 = 0.2 for 0.05 s < d ≤ 0.1 s, 3 = 0.1 for 0.1 s < d ≤ 0.2 s and 4 for d ≤ 0.2 s. The picks with quality 4 are excluded from the further analysis.

3.1 Estimation of relative travel-time residual

After estimating relative travel time residuals using Eq. (2), these residuals are corrected from already known heterogeneous crustal anomalies and the travel time anomalies from the top sediments, which generally cause unwanted smearing effects of the residuals along the steep teleseismic ray paths [15, 16]. An available detailed crustal velocity model is used for this purpose and then these residuals are corrected from crustal travel time anomalies (r_ij^crust) relative to the background earth model using the relation:

$${{\mathrm{t}}_{\mathrm{ij},\mathrm{c}}}^{\mathrm{m}}={{\mathrm{t}}_{\mathrm{ij}}}^{\mathrm{m}}-{{\mathrm{r}}_{\mathrm{ij}}}^{\text{crust}}$$E3

3.2 Crustal correction

It is well-established fact that significant travel time residuals can be caused by the heterogeneous continental crust. Further, large delays of the order of 1 s can be resulted from the low-velocity sediments in the rift basins (like Kachchh) or/either the Moho topography, which laterally replaces mantle (Vp ~ 8 km/s) and crustal (Vp ~ 6.5 km/s) material. Thus, the relative residuals should be corrected from such known crustal residuals to study the deep mantle structure. Therefore, these crustal anomalies should be determined a priori and subtracted [16, 17]. Following Martin et al.’s [16] approach, here the crustal corrections are applied using the iasp91 reference earth model consisting of a 20 km thick upper crust with a Vp of 5.8 km/s and a 15 km thick lower crust with a Vp of 6.5 km/s. To apply crustal correction, top sedimentary layer is assumed to be composed of Cenozoic sediments (Vpsed ~ 2.41–4.60 km/s) and Jurassic/Mesozoic sediments (Vpsed ~ 4.61–5.58 km/s). The estimated corrections of the travel time residuals for such 1-D known crustal velocity models in the region vary from 0.29 to 0.79 s for P-waves.

A total of 241 teleseismic events (M_w 6.0–8.0) have been analyzed to estimate crust corrected P-residuals using Eq. (3), they are found to be stable for stations in the Kachchh network (Figure 2). From Figure 2, we notice that there is no clear anisotropic signal in the P-residuals in most of Kachchh region except the Motapaya (MTP) station, which is the northernmost station in Kachchh. Figure 2 shows a plot of estimated crust corrected relative P-residuals at 14 stations as a function of slowness (p) (0–10 s/o along the radial axis) and back azimuth (BAZ 0–360°) to the event. In general, early P-arrivals are noticed for events in the east of the network while delayed P-arrivals are observed for events in the west. Most of the stations, which are lying in the Deccan basalt covered regions in the Kachchh and Cambay rift zones, show negative residuals. This could be attributed to the presence of significant crustal and lithospheric upwarping below these rift zones [18, 19]. At MTP, travel time residuals show images of the spatial variation of wave propagation anomalies of P-waves, which can be attributed to an anisotropy structure of the mantle domains (Figure 2a). This anisotropic nature of P-waves is not seen for other stations due to fewer numbers of observations.

The directional term of relative crust-corrected P-residuals is calculated by dividing the 0–360 to 18 quadrants, and then we calculate the average residuals of stations falling in a particular quadrant. And finally, we estimate the arithmetic mean of the average residuals of stations falling in all 18 quadrants. Subsequently, this mean of residuals is being subtracted from the relative crust-corrected P-residuals for different stations for estimating the directional term of relative crust-corrected P-residuals. Following the above procedure, we also computed directional term of relative crust-corrected P-residuals, which are plotted in Figure 2b. The directional mean of P-residuals, which are obtained by subtracting average residuals of all stations in Kachchh, is shown in Figure 2b. While Figure 2b represents the directional mean of P-residuals, which are calculated by subtracting average residual of 8 stations in Kachchh. Most interestingly, two distinct zones of negative P-residuals associated with Kachchh and Cambay rift zones are noticed from Figure 2b. These zones of negative P-residual probably mark the regions of marked lithospheric thinning. This model gets further support from available estimates of crustal and lithospheric thicknesses, suggesting a marked crustal and lithospheric thinning below the Kachchh and Cambay rift zones [18, 19].

4. P-wave teleseismic tomography

The teleseismic tomography code developed by Weiland et al. [20] and later modified by various investigators [21, 22, 23, 24] has been used here to estimate the 3-D P-velocity structure down to 250 km underlying the Kachchh rift zone (KRZ). The data collected by the digital networks of NGRI (Hyderabad, India) and ISR (Gujarat, India) are combined, resulting in a dataset from a total of 59 stations (Figure 2b). Here, the P-wave teleseismic tomography is performed using relative crust-corrected P-residuals as estimated above. First velocity perturbations in this method are calculated in the orthogonal net of nodes approximating the volume under consideration. And, the trilinear interpolation is used to calculate velocity at adjacent nodes [25]. The initial velocity model of the upper mantle and theoretical travel times are set according to a reference Earth model IASP91 [14]. We run the tomographic code with three-dimensional ray tracing implemented by the Simplex method [25]. The kernel matrix is inverted by Singular Value Decomposition (SVD). The basic equation of the inversion of the tomographic code (i.e., TELINV) can be written as:

$${\mathbf{m}}_{\mathbf{est}}={\left({\mathbf{G}}^{\mathrm{T}}{\mathbf{W}}_{\mathrm{D}}\mathbf{G}+{\mathrm{\epsilon }}^2\mathbf{I}\right)}^{-1}{\mathbf{G}}^{\mathrm{T}}{\mathrm{W}}_{\mathrm{D}}\mathbf{d}$$E4

where m_est and G are estimated model parameters (i.e., velocity perturbations) and matrix of partial derivatives with respect to the model parameters, respectively. W_D is weighting matrix of data, where weights are set according to quality factors assigned to individual arrival time picks. ε² and I are a damping factor and a unit diagonal matrix, respectively. And, d is data vector (i.e., relative residuals) [26]. Several iterative cycles are performed for satisfying assumptions behind the linearization of the inversion. New ray-paths inside of the area studied are traced at each cycle using the improved velocity model retrieved in the previous step. Here, four iteration cycles are applied to reach such data variance which does not decrease with further iterations noticeably and stay above the twice the average data error.

4.1 Model parameterization

The primary study area extends about 320 km in E-W and 240 km in N-S, which covers the E-W trending Kachchh rift zone (Figures 1a and 2b). Here, the volume below the central study area is approximated by the 3D (i.e., x, y and z directions) grid of nodes (17 × 18 × 15). For teleseismic tomography, a much large area covering 10,000 km × 10,000 km × 530 km, is assumed to stabilize the velocity perturbation outside the main study region. Seventeen X-nodes are distributed at distances of −5000 −240 −200 −160 −120 −80 −40 0 40 80 120 160 200 520 560 600 5000 km while eighteen Y-nodes are located at distances of −5000 −360 −320 −280 −240 −200 −160 −120 −80 −40 0 40 80 120 160 200 240 5000 km. And, fifteen nodes in Z directions are at depths of −5 10 70 50 90 130 170 210 250 290 330 370 410 450 490 km. The center of the array is assumed to be at latitude 23.3° and longitude 70.3°. Here, the teleseismic tomographic inversion is performed using a total of 1788 P-residuals, which are estimated using vertical component of broadband seismograms of 241 good teleseismic events (with epicentral distances between 30° and 90°) recorded at 59 stations. Here, the travel time residuals range from −2.276 s to 2.096 s. This study results in an average data error of 0.062 s. For the present study, a travel-time residual of 2.5 s is used while a ε² damping factor of 100 is used for the inversion. To minimize effects of potential inaccuracies of the 3D crustal model applied in the inversion, the first inverted layer of nodes is assumed at 50 km depth. Several different vertical parameterizations with irregular and smaller spacing are being tested, but the variance reduction, as well as the diagonal terms of the resolution matrix, decreased rapidly. In total, the inversion is performed for 980 model parameters.

The crustal structure cannot be resolved through teleseismic tomography due to the sub-vertical directions of the incoming rays. However, the heterogeneous, complex crustal structure can affect the inverted travel-times significantly [27, 28, 29]. Further, the size of crustal heterogeneities is often exceeded by the spacing between stations. Therefore, crustal corrections are inevitable to apply on the travel-time residuals, before inverting the residuals for the upper mantle velocity structure. Here, following the standard procedure [17, 29] crustal corrections as discussed above are applied to the dataset. Modeled dVp (%) tomograms at different depths are shown in Figure 3a–d.

Figure 3.

P-wave velocity perturbations at horizontal slices imaged at (a) 50, (b) 90, (c) 130, and (d) 170 km depths, by the P-wave tomography using the crust corrected relative P-residuals. Major faults (solid white lines): ABF—Allah bund fault; IBF—Island belt fault; KMF—Kachchh mainland fault; KHF—Katrol hill fault; NPF—Nagar Parkar fault; BF—Banni fault. And, NWF (north Wagad fault), the causative fault for 2001 Bhuj earthquake and Gedi fault, are shown by dotted white lines. Locations of mobile broadband stations are marked by gray open triangles. A black dotted elliptical area mark the central KRZ.

4.2. Resolution analyses

Resolution analyses of tomographic models are essential for distinguishing real structures from artifacts caused by methods used and for identifying well-resolved model parameters. The reliability of tomographic images have been analyzed by performing a checkerboard test, which is performed with the ray geometry of the Kachchh network, with 3D ray-tracing, and with the damping factor of 100.

To verify sensitivity in the whole volume studied, a checkerboard test is performed [30]. For this, a net of alternating anomalies of +5% and −5% in nodes is constructed at depths of 50, 90, 130, and 170 (Figure 4a–d), leaving the layers in between them unperturbed as well as in the remaining parts of the model [30]. The input anomalies are recovered well at depths ranging from 50 to 170 km (Figure 4e–h). The low-velocity anomalies in the central part are not resolved successfully at ~210 km depth, which could be attributed to the well-known vertical smearing dominated in the inversions of teleseismic data.

Figure 4.

Input 3-D synthetic P-wave velocity model with alternative grids of high (±5%) and low (−5%) velocity for the checkerboard test at (a) 50, (b) 90, (c) 130, and (d) 170 km depths. Results of checkerboard test showing retrieval of synthetic velocity perturbations at (e) 50, (f) 90, (g) 130, and (h) 170 km depths.

5. Results and discussions

Large negative travel time residuals associated with the central Kachchh rift zone are detected (marked by a black dotted elliptical area shown in Figure 2b), suggesting relatively thin lithosphere underlying the central KRZ while surrounding unrifted parts show positive travel time residuals suggesting relatively thick lithosphere as also observed in the estimates of lithospheric thicknesses [18]. Large negative travel time residuals are also found in the Cambay rift zone (CRZ) (Figure 2b) suggesting relatively thin lithosphere underlying the CRZ as also modeled by an S-RF imaging and surface wave dispersion studies [19, 31]. Note that the estimated directional means of relative P-residuals show a good spatial correlation with the major tectonic regional features. The most negative values (though with large standard deviations) follow the track of the Cambay rift or Deccan plume and basaltic region of Kachchh. This observation shows the robustness of our modeling results.

The horizontal slices of dVp tomograms at different depths (50, 90, 130, and 170 km) suggest that a prominent low-velocity zone (~−2 to −4%) is present below the central Kachchh rift zone extending from 50 to 170 km depths (Figure 3a–d). At 50 km depth, a high-velocity structure is also noticed in the region north-east of the 2001 mainshock location (Figure 3a), which vanishes at deeper depths. But, the low-velocity zone in the KRZ becomes a bit subdued (~−1 to −3%) between 150 and 170 km depths (Figure 3c and d). However, this low-velocity zone also exists at 210 km depth, which might be an artifact due to the poor sampling of rays for the deeper zone. Tomograms also clearly bring out the fact that in the upper mantle (at 70–170 km depth) the P-wave velocity is found to be higher in the western part in comparison to that on the eastern part of the Kachchh region. In the shallower depths (at 50 km), tomogram suggests that the central KRZ, which is characterized by the crustal as well as lithospheric thinning and negative traveltime residuals (Figure 2b; [18]), is found to be associated with the low dVp anomaly. This low-velocity zone associated with the central KRZ also exists between 90 and 170 km depths (Figure 3b–d), which could be attributed to the presence of carbonatite/partial melts related to the Deccan volcanism of 65 Ma [19, 31]. Note that the presence of carbonatite melts in the shallow upper mantle depths below the KRZ has also been proposed by an isotopic ratio study of xenoliths [7] and the modeling of SKS splitting [11]. Thus, the volatile CO₂ (emanating from the crystallization of carbonatite melts in the asthenosphere), can reach lower crustal depths, through deep-seated faults, which might be facilitating the deeper circulation of metamorphic fluids/volatile CO₂, thereby, the generation of lower crustal earthquakes occurring in the Kachchh rift zone [18].

The most robust outcome of the tomography, which does not vary with regularization or parameterization applied, is that low-velocity perturbations prevail in the Kachchh rift zone (as shown by dotted elliptical area in Figure 2b) down to ~170 km depth (see Figure 3a–d). Similarly, earlier studies [8, 10] indicated that the upper mantle beneath the entire KRZ is characterized by lower velocities, about the surrounding area down to these depths. These studies do not exploit only teleseismic data but also include data from regional events. Moreover, regional and global teleseismic tomography studies [8, 32] show the upper mantle beneath the Kachchh region as a part of an extensive low-velocity heterogeneity located in western India, extending down to about 700 km. The regional teleseismic tomography from the western India indicates relatively small velocity variations (±1%), both in size and in amplitude. Such details are below resolution level of global or regional studies which invert data of permanent observatories only. Thus, our results can be viewed as being in agreement with these larger-scale investigations.

6. Conclusions

The variation of crust corrected travel time residuals correlate well with our estimates of crustal and lithospheric thicknesses in the study region. The large negative residuals are found to be associated with the Kachchh and Cambay rift zones, which are characterized by a marked crustal and asthenospheric thinning while positive residuals describe the unrifted regions.

Our tomograms reveal a distinct low dVp anomaly at 50–170 km depths underlying the central Kachchh rift zone. However, this anomaly vanishes for deeper regions (>210 km) due to poor ray sampling. At 90–130 km depths, the resulting tomographic model reveals a slower P-wave speed on the eastern part of Kachchh compared to that of the western part. A slow dVp anomaly is mapped at 50–170 km depths below the central KRZ, which is also characterized by the crustal as well as lithospheric thinning and negative traveltime residuals. This upper mantle low dVp anomaly could be attributed to the presence of carbonatite/partial melts related to the Deccan volcanism of 65 Ma. We propose that the volatile CO₂, which is emanating from the crystallization of carbonatite melts in the asthenosphere, is playing a crucial role in generating lower crustal earthquakes occurring in the Kachchh rift zone [18].

The checkerboard test results suggest that our tomographic model is quite well resolved for the central KRZ at 50–170 km depths while tomographic model at more than 210 km depth is poorly resolved due to poor ray sampling. The checkerboard test reveals that 85–90% of the synthetic velocity anomaly could be retrieved at 50–170 km depths below the central KRZ (Figure 4e–h). Thus, the velocity model below the KRZ seems to be well-resolved. However, the retrieval of synthetic velocity for the surrounding regions are not very good (<80%). To increase the resolution of the tomographic images of the upper mantle velocities, a future study could be planned to combine Kachchh dataset with those from other broadband stations in Gujarat and Rajasthan, using more teleseismic events with epicentral distances between 30 and 90^o and back-azimuths between 180 and 270^o.

Acknowledgments

The author is grateful to the Director, Council of Scientific and Industrial Research—National Geophysical Research Institute (CSIR-NGRI), Hyderabad, India, for his support and permission to publish this work. This study was supported by the Indo-Czech (CSIR-ASCR) collaborative project no, IND 2012/19.