“ CRUSTAL STRUCTURE OF ANDAMAN ISLANDS USING JOINT INVERSION OF RECEIVER FUNCTIONS AND SURFACE WAVE DISPERSION MEASUREMENTS „

We investigate the crustal structure of Andaman Islands, central part of Burma-Sunda-Java subduction complex, through joint inversion of receiver functions and surface wave dispersion measurements. For this study, we used teleseismic earthquakes recorded over 13 temporary broadband seismographs operated in two different phases: pre and post 2004 Sumatra earthquake. Beneath the Andaman forearc region, the crustal thickness varies between ~24 and 28 km. The uppermost crust consists of ~4 to 7 km thick soft accretionary sediments. Average Vp/Vs ratio, between ~1.79 and 1.83, suggests accretionary hydrated oceanic crust with different level of saturation. Combining derived crustal structure with global seismicity and CMT fault plane solutions indicate a complex convergence along the arc, with transitional (continental-oceanic) type in the north to oceanic type in the south Andaman.


INTRODUCTION
The Andaman Islands are part of the Andaman-Nicobar Ridge (ANR) which hosts Andaman and Nicobar group of Islands. These Islands mark the eastern margin of the Indian plate and form an important transitional tectonic link between the eastern Himalayan syntaxis in the north and Sunda arc in south. These represent the central part of the ~5000 km long Burma-Sunda-Java subduction complex, displaying major tectono-stratigraphic elements striking approximately parallel to the trend of subduction trench. The tectonic framework of the region has been reviewed by various researchers [e.g., Curray, 2005;Lay et al., 2005;Kamesh Raju et al., 2012 and reference within].
Global plate tectonic reconstructions suggest the complex convergence of Indian plate along Southeast Asian margin, which has resulted in clockwise rotation of the subduction zone and increase in the obliquity [Replumaz et al., 2010]. Subduction in this region is pre-sumed to have started in lower Cretaceous [Scotese et al., 1988], which occurs all along the Sunda arc and extends from the eastern Himalayan syntaxis to Banda arc [Curray, 2005]. The age and thickness of the subducted oceanic crust and convergence rate increase from Andaman towards Java along the arc . This change is observed in the increasing dip and depth of penetration of the wadati-Benioff zone, causing change in subducting slab geometry. Oblique, but predominantly thrust motion occurs in the Andaman trench with a convergence rate of about 1.4 cm/yr . The convergence varies from continental type (Indian continental plate vs Burmese plate) in the Burmese arc to oceanic type (Indian oceanic plate vs Burmese plate) in the Andaman arc [Subrahmanyam et al., 2008]. It is believed that early subduction was comparatively fast which subsequently progressed in multiple episodes as short-and long-lived subduction zones converging at different rates [Richards et al., 2007]. This convergence (at variable rate, increased obliquity, dip and rotation) has resulted in development of a plate sliver that is referred as Andaman or Burmese mi-croplate, which is a sheared off plate parallel to the subduction zone from Myanmar to Sumatra. Prominent morphological features in the region includes several sea mounts, volcanic arc, Barren Island (BaI), Narcondam Island (NaI), Andaman Backarc Spreading Center (ABSC), Alcock and Sewell seamount complexes, Backarc basin ( Figure 1).
Geologically, Andaman Islands are exposed tectonostratigraphic units of an accretionary prism in an outer arc setting and turbidities of a forearc setting. Rocks in this region belong to the upper Cretaceous to Tertiary, with Oligocene flysch covering the western part and the Paleocene to Eocene sediments of Mithakari covering the eastern part. Some Ophiolites belonging to the Mesozoic lower Cenozoic also occur in small patches in different places [Pal et al., 2003].
Seismotectonics, nature of faulting and stress distribution pattern for the Andaman arc have been studied earlier [e.g. Mukhopadhyay, 1984;Guzman-Speziale and Ni, 1996;Radhakrishna and Sanu, 2002;Dasgupta et al., 2003;Khan, 2005]. Predominant seismic activity in Andaman region is attributed to subduction of In-MISHRA ET AL.  ham et al., 2005;Engdahl et al., 2007]. The depth of the earthquakes in the region ranges between ~150 and 300 km; which seems to increase from ~150 km in Andaman region (at around 12° N) to almost ~300 km in the south (at around 4° N) [Sorensen et al., 2007)]. Tomographic images suggest subhorizontal tear in the subducting slab below Burma [Replumaz et al. 2010;Pesicek et al., 201;Mishra et al, 2020]. Curray [2005] and Mishra et al. [2011] suggest the presence of strong heterogeneity in crust and upper mantle, as the lithosphere in the Andaman sea region underwent complex tectonic deformations during the Neogene-Quaternary period. An earthquake of Mw 9.3 (popularly known as Sumatra-Andaman Earthquake) hit the region on 26 December 2004. Several studies were carried out after this mega thrust earthquake highlighting shift in seismogenic coupling zone and ~1400 km rupture zone extending up to Andaman Islands [Kennett and Cummins, 2005;Ammon et al., 2005;Ishi et al., 2005;Grevemeyer and Tiwari, 2006;Shapiro et al., 2008]. This mega event also reactivated Barren volcanic Island which indicates the tight link between tectonics and magmatism of this complex system [Franke et al., 2008]. Ever since, this area has witnessed more than 17,000 aftershocks [Mishra et al., 2011].
Despite being one of the most active terrain and its tectonic linkage to Sumatra subduction system, so far, no systematic study related to crustal structure has been carried out particularly beneath Andaman Islands. Geographical limitation, logistic difficulties, dense forest, aboriginal reserve areas could be possible reasons. To study the crust -upper mantle structure and seismicity of the Andaman region, we operated 13 broadband seismic stations (Figure 1b). In the present work, we derive the crustal structure of Andaman Islands using joint inversion of receiver functions and surface wave dispersion measurements. This information was combined with global seismicity and CMT fault plane solutions to provide some insight into subduction model beneath Andaman Islands.

DATA AND METHODOLOGY
The data used in the present study was recorded by 13 broadband seismic stations operated in two different phases. During November 2003 to February 2004 (phase 1), we deployed 5 stations (orange triangles in Figure   1b), and during January to May 2005 (phase 2) we deployed 8 stations (red triangles). In phase 2, just after the 2004 Sumatra earthquake, we re-occupied 3 previous locations as in phase 1 (pink triangle). The seismological stations configurations included Guralp CMG-3T sensors with a flat velocity response between 0.008 -50 Hz and REFTEK 130-01 data loggers. Data were continuously recorded at 50 samples per second and the corresponding Global Positioning System (GPS) time was logged.
Receiver function (RF), a well-known and established technique, was used to study the crustal structure of Andaman Islands. Receiver functions are radial and transverse waveforms created by deconvolving the vertical component from the radial and transverse components of the seismogram to isolate the receiver site effects from the other information contained in a teleseismic P-wave [Langston, 1979;Ammon, 1991;Liggoria and Ammon, 1999].
To compute receiver functions, teleseismic waveforms of the earthquakes with magnitudes above 5.5 and epicentral distances between 30 o to 95 o were selected. This epicentral range avoids multiple arrivals in the direct P wave field occurring at distances less than ~30° due to triplications caused by the rapid velocity increase in the upper mantle transition zone, and complications at distances greater than ~95° resulting from the core-mantle boundary. Figure 2 shows the location of earthquakes used in this study and recording station network. Due to geographical location of the study region, most of the earthquakes are from NE and SE directions.
We computed receiver function using iterative time domain deconvolution approach [Liggoria and Ammon, 1999]. Different Gaussian filter widths were tested, however, as our objective was to find the first order discontinuities, we preferred RFs with 1.6 Gaussian width, corresponding to low pass filter with a corner frequency at ~0.8 Hz. To further control the quality of waveform, we used only the RF with variance reduction cut-off above 80%. To equalize the effect of variable distances, RFs were moveout corrected to a reference slowness of 6.4 s/deg corresponding to an epicentral distance of 67 o [Yuan et al., 1997] using IASP91 model. Out of total 1329 RFs calculated, only 269 good quality RFs were selected for further analysis. Figure 3 shows moveout corrected radial receiver functions calculated at individual seismic stations. IMD, RGT and HVL stations which operated during both the phases of experiment have maximum RFs, but due to technical problems HVL station could not record much earthquakes and produced lesser RFs.
Based on geographical setting, the whole study region is divided in three segments (North, Middle and South Andaman; Figure 1b); and the RF analysis is presented for these three segments.

NORTH ANDAMAN: [KGT]
RFs obtained at KGT (Kalighat) show a strong coherent phase (positive polarity) arriving at ~3.8 s. Other prominent observed phases (conversion and multiples) are at ~12-14 s and show strong back azimuthal dependence ( Figure 3).

MIDDLE ANDAMAN: [TGP, RGT, BTG]
Both TGP (Tugapur) and RGT (Rangat) seismic stations in this part of Andaman show similarity in their waveform. Receiver functions at these stations show prominent phases (positive polarity) at ~2 s, ~5 s and ~13-14 s ( Figure 3). Seismic station at BTG (Baratang), close to a mud volcano, exhibits complex RFs ( Figure 3). This mud volcano erupted a huge volumes of slurry ma-terial after the main Sumatra-Andaman earthquake and the eruption has been associated with ascending fluids and gases at shallow subsurface level. At BTG, we found that the receiver functions from northern backazimuth (310 o to 50 o ) events produces a visible shift in direct P. This shift in P could possibly be due to these low velocity layer at the subsurface below BTG. Other prominent visible phases are at ~4-5 s and ~10 s ( Figure 3).

SOUTH ANDAMAN: [TIR, IMD, HTN, BTG, HVL, NIL]
The western station, TIR (Tirur) shows relatively simpler receiver functions with a strong coherent phase arriving at ~3.5 s, other phases (conversion and multiples) at ~5-6 s and ~15 s (not clear on all the RFs) (Figure 3). At central stations IMD (Portblair) and HTN (Hopetown) in this part of Andaman, we observed the Ps conversions at ~3 s, and this phase get shifted (~4.5 s) for NW direction (>270°) events ( Figure 3). Other conversion at these stations include a strong phase at ~14 s. CHI (Chidiyatapu), the southernmost station, shows strong conversions at ~2 s and ~12 s (Figure 3). At the eastern stations HVL (Havelock) and NIL (Neil) in the south Andaman, we observed coherent phase at ~3-3.5 s with its multiple (PpPms) around ~12 s ( Figure 3). One more strong conversion is observed at ~7 s with its multiple at ~20 s ( Figure 3). In general, the multiples are masked by various other phase(s), which makes it difficult to be distinguished at most of the stations. This is an inherent problem arising in complex subduction regions e.g., Cocos [Chang and Baag, 2007], Chilean [Dzierma et al., 2012], Korean [Kim et al., 2010], Sumatra [Macpherson et al., 2012]. Hence, we have restricted our analysis mainly based on clear conversions.
To study the nature of Moho conversions (Ps phase) and other intra-crustal phases, we have projected few RFs along the two profiles AA' and BB' (Figure 4a). These RFs were calculated for the earthquakes from SE direction (Figure 4b), recorded on majority of the stations and traversing similar ray path. These receiver functions were corrected for the distance moveout of the Moho converted Ps phase referenced to 67°. The average of the summed receiver function is presented (left and upper panel of Figure Figure 4d). As we move away from trench to easternmost station (HVL), we observed a strong conversion arriving at around ~3 s. Due to limited data at this station we are unable to track its continuity. However, in consultation with neighbouring stations, we believe 3 s to be possible Ps from top crust and a strong conversion at 7 s could be from subducting Indian slab.

JOINT INVERSION
To obtain the seismic crustal structure beneath Andaman Islands, we have inverted the obtained radial receiver functions (RFs). The inversion of receiver functions is essentially a non-linear approach. Being a complex terrain, any additional weights to constrain inversion solution can provide reliability to results while converging towards the final solution. RFs are sensitive to the shear wave velocity contrast and has no control on the absolute shear velocity, whereas surface wave dispersion measurements constrain average shear velocity that reach deeper structure with increasing period. Since both RFs and surface wave dispersion are primarily sensitive to shear wave velocity, these can be inverted jointly to obtain a reliable estimate of shear velocity structure. Merging information from both receiver functions and surface wave dispersion into a single inversion algorithm limits the non-uniqueness inherent in receiver function and provides better constraints on shear wave velocity measurements.
Due to insufficient backazimuthal coverage of the earthquake, geographical location and duration of experiment, we could not analyze effect of dip, anisotropy, scattering etc. and restricted ourselves for the determination of simple 1-D velocity model beneath every station. Only selected receiver functions, sampling the similar ray path from a particular backazimuth (red stars in Figure 4b) were used for the joint inversion. Further, considering the complexity of terrain, we have inverted all individual receiver function at each station to obtain an average velocity structure. We used relatively small time window of 12 s for modeling to avoid any contamination due to multi pathing and other multiple phases.
Receiver functions and surface wave dispersion measurements were jointly inverted using iterative linearized damped least-square scheme [Julia et al., 2000;Herrmann and Ammon, 2004] which incorporates a priori smoothness constraints for velocities in adjacent layers. The 15 -45 s period fundamental mode Rayleigh wave group velocities curves were extracted for each individual station in the Andaman region from surface wave tomography results of Acton et al. [2010]. Group velocity measurements represent an average picture of the region and cannot reflect the rapid change in structure of the region. Hence, a careful consideration was given for the sensitivities of each dataset towards underlying crustal structure and the size of the sampling region. Greater weight is given to fit the receiver function data during the joint inversion to obtain local structure beneath each seismic station. The joint inver-sion was performed for a range of weights and final models were selected based on the best fit to the receiver function, whilst maintaining an adequate fit to the dispersion data. In particular, model group velocities at short periods are allowed greater deviation from the observed group velocities to reflect the rapidly varying structures in the shallow part of crust.
In order to avoid biasing the inversion model, the starting model for the inversion at each of the sites was the same and consisted of the AK135 [Kennett et al., 1995] velocity model with the crust and sub-Moho mantle replaced with 4.48 km/s layers upto 60 km depth. The thickness of starting model was parameterized as homogeneous and isotropic layers of 1 km thickness until 34 km (2 km layer thereafter) to account for any velocity variation in the model for any possible variation in the subsurface velocity. A constant Vp/Vs value of 1.73 was assumed during each inversion as a starting point. Since receiver functions and surface wave dispersion measurements were inverted primarily for retrieving crustal structure, the model damping was chosen to be large for deeper mantle. Christensen and Mooney [1995] and Christensen [1996] has shown that shear wave velocity (Vs) in the lower crust cannot exceed 4.3 km/s, and Vs above this indicates the presence of lithology of mantle composition. Therefore, the depth where Vs >4.3 km/s was used as a marker to compute the Moho depth.
In order to test the robustness of our inversion results, the starting (input) model was perturbed by 5% of input velocity model and final output results were compared. As shown in Figure 5, perturbing the starting model has negligible effect on the final velocity models by the inversion process. In both the cases of velocity perturbation (increase or decrease), we could hardly find any significant variation in the output velocity models down to ~50 km depth; which shows the stability of the final inversion results.
The inversion results obtained from joint inversion of receiver function and surface wave dispersions at stations in different sectors of Andaman are presented in Figure 6. The velocity structure at seismic station (KGT, north Andaman) is well constrained with RFs and surface wave dispersion measurements (Figure 6, north Andaman). In Middle Andaman, TGP shows good matching between observed and modeled surface wave measurements, however RFs are not matched properly. In contrast, at RGT, RFs (observed and modeled) are well matched while, surface waves measurements have diversions at longer periods, however they are in ±1s error limits (Figure 6, middle Andaman). At seismic stations in south Andaman, RFs and surface wave disper-sion measurements (observed and modeled) are well matched and better constrained by the velocity model ( Figure 6). At TIR (in east part of south Andaman) both RFs and surface waves dispersion curves are well modeled ( Figure 6). At HVL (west part of south Andaman, we had lesser number of receiver functions. While surface wave measurements are well matched, RFs are not matched properly during modeling.

H-Vp/Vs STACKING METHOD
To quantify the crustal thickness and Vp/Vs ratio in the vicinity of each station, we modelled the amplitude and travel times of P-to-S conversions at the Moho (Ps) and its crustal multiples (PpPms and PpSms + PsPms) in the radial receiver function using the grid search algorithm [Zhu and Kanamori, 2000]. This algorithm has been successfully applied for obtaining Moho depths and Vp/Vs ratios for subduction zones e.g., Cocos [Chang and Baag, 2007], Korean [Kim et al., 2010] and Chilean [Dzierma et al., 2012]. This algorithm exploits the fact that arrival times and amplitude of specific Moho converted phases and multiples appearing on radial receiver functions are determined by known functions of Moho depth (H), Vp/Vs ratio. For a near true combination of H and k value, the quantity S(H, Vp/Vs) is defined as the weighted sum of the receiver function amplitudes at the calculated times of predicted arrivals of Ps, PpPms and PpSms+PsPms phases would be expected to be maximum.
where r j (t) is the amplitude of receiver function for the jth event, t 1 , t 2 , t 3 are predicted Ps, PpPms and PpSms+PsPms arrival times corresponding to Moho depth H and Vp/Vs. Since travel times used for crustal receiver function analysis are much sensitive to Vs than to Vp, we assume an average Vp (6.31 km/s, Pesicek et al. 2010) for the entire crust and perform the grid search for a large number of crustal models with varying thicknesses H (10 -80 km) and varying Vp/Vs (1.6-2.0). Being a complex region, it was difficult to find welldefined global maximum and a strong trade-off existed between crustal thickness and Vp/Vs ratio (Figure 6d, for each seismic station), therefore Moho depth obtained by joint inversion was used to obtain Vp/Vs ratio beneath each station. Figure 6d, shows the results of H -Vp/Vs for each station along with inversion results. From the present study, Vp/Vs ratio ranges from 1.79 to 1.83 for the entire study region. Due to limitation of the 9 datasets, it is very difficult to conclude the nature of curst, hydration or its level of saturation with fluids. However, possibilities of the presence of fluids generated by the metamorphism at deeper levels of curst in subduction zones [Peacock, 1990] cannot be ignored.

RESULTS AND DISCUSSION
To follow the geometry of the subducting slab, crustal thickness estimates obtained from joint inversion is interpreted together with hypocentral locations, fault plane solutions and piercing points (using ray path of event-station pair). Hypocentral locations are obtained from high precision relocated earthquake from EHB catalogue (Engdahl et al. 1998(Engdahl et al. , 2007, while fault plane solutions are from Harvard CMT (Centroid Moment Tensor) solutions (ISC 2010). Figure 7 shows the depth cross-sections along the profile EE' (North Andaman), DD' (Middle Andaman), CC' (South Andaman) (as shown in Figure 4a) obtained from seismicity, superimposed with our crustal depth estimates and focal mechanism solutions. In the following section, we have discussed these results by diving Andaman Islands into three segments (North, Middle and South Andaman Islands).  KGT (Kalighat), the northern seismic station in the Andaman Islands, shows a crustal thickness of ~28 km from joint inversion result. Top ~6-7 km is overlaid by thick low shear velocity (Vs ~2.5 km/s) Andaman flysch sediments (Figure 6d, North Andaman). Satellite gravity anomaly modeling and qualitative interpretation of RF results shows 30 km thick double oceanic crust, comprising 9 km Indian crust followed by 21 km Burmese crust (Rao et al., 2011). With gross crustal thickness (~28 km), our modeling results is in close agreement with Rao et al. [2011]. Depth cross-section plots along profile EE' (as shown in Figure 4a) was generated by combining the current modeling results with ray paths, seismicity and fault plane solutions ( Figure  7a). As shown in figure 7a (North Andaman), the majority of the earthquakes occur in the depth range of ~10-30 km with main concentration around ~20-30 km. This shallow seismicity has mainly thrust dominated strike-slip motion. This could be due to the seismogenic coupling zone between subducting Indian slab and overriding Burmese crust, as shown by Grevemeyer and Tiwari [2006] and Shulgin et al. [2013] in the northern Sumatra. Focal mechanism solution reveals East (to North-East) gently dipping seismicity trend, which coincides with the direction of subducting Indian slab. The subducting slab has relatively low dip up to ~70 km from the trench and gets steeper further away from the trench. The shallow seismicity cluster is suggestive of intense tectonic deformation the subducting slab undergoes before descending into mantle. In absence of observed seismicity after ~70 km depth, except one earthquake at ~130 km depth, we can only guess the slab penetration in this part of Andaman.

MIDDLE ANDAMAN
In middle Andaman, both the seismic stations TGP (Tugapur) and RGT (Rangat) exhibit almost similar seismic structure (Figure 6d, Middle Andaman). At both the stations, two distinct interfaces (~15 and ~24 km depths) are observed. Based on RF analysis and modeling results, ~24 km is considered as the Moho depth ( Figure  6d, Figure 7b). Below the Moho, shear velocity drops further (Vs ~ 4.0 km/s) and reaches back ~4.5 km/s at a depth of ~55 km (Figure 6d Figure 4a), by combining modeling results with receiver function ray paths, known seismicity and fault plane solutions. The ~55 km interface could be signature from subduction Indian slab, being tracked down to a depth of ~120 km (Figure 7b). The low velocity between ~24 and 55 km could be interpreted as hydrated material sandwiched between the Burmese crust and subducting Indian slab. The earthquakes at a depth of about 100 km in the subducting Indian slab [Dasgupta et al., 2003] can be related to a phase transition in the slab causing fluid release and/or partial melting of the oceanic crust.

SOUTH ANDAMAN
Inversion result from TIR (Tirur), the western and closest station to the trench, shows the crustal thickness of ~28 km. HTN (Hopetown) and IMD (Portblair), in the central part of south Andaman segment, show crustal thickness of ~24 km. HVL (Havlok), the easternmost and farthest station from trench, shows a crustal thickness of ~26 km. Depth cross-section in this segment ( Figure 7c, along profile CC' in Figure 4a) shows Indian slab is penetrating down to ~160 km. The slab convergence is more of an oceanic type as suggested by Subrahmanyam et al. (2008). Focal mechanism of earthquakes at depth >100 km is of normal type with low dipping eastward dominating strike, which can be related to plate bending deformation due to strong tensile (slab-pull) forces acting on subducting slab [Astiz et al., 1988;Conrad et al., 1999].
Crustal thickness obtained from joint inversion beneath the Andaman forearc region varies between ~24 and ~28 km. Joint inversion of receiver functions and surface wave dispersion measurements beneath Sumatra region also shows similar crustal thickness (Macpherson et al. 2012). Subsurface low shear velocity (Vs ~1.3 -2.5 km/s) layer, with an average thickness between ~4 and 7 km, is observed at almost all the seismic stations. This low velocity layer can be interpreted as thick soft Andaman flysch sediments. Magnetotelluric study in the region [Gokarn et al., 2006] also report similar sedimentary thickness. Average Vp/Vs ratio from ~1.79 to 1.83 represents accretionary hydrated oceanic crust with different level of saturation due to water released from downgoing slab as observed in other subduction zones [Peacock, 1990;Chang and Baag, 2007;Audet et al., 2009;Kim et al., 2010]. Earthquakes along the Andaman arc are mostly attributed to the ongoing subduction of Indian slab beneath the Burmese plate. The hypocenter distribution indicates a variable dip of the subducting Indian slab, similar to previous studies [Mukhopadhyay, 1988;Radhakrishna et al., 2008]. Seismicity, dip angle and focal mechanism solutions indicate a complex oblique convergence along the arc with transitional (continental -oceanic) type in the northern segment to oceanic type in the southern Andaman. Further south, in the 2005 Sumatra earthquake rupture area, the deep structural image using seismic tomography modeling of wide-angle ocean bottom data also shows differences in the structure of the subduction system through crustal-scale changes (including crustal thickness change) and possible effects on the variations in the seismogenic behavior [Shulgin et al., 2013].

CONCLUSION
We investigate the crustal configuration of the Andaman Island using joint inversion of receiver functions and Rayleigh wave dispersion measurements. The crustal thickness varies between ~24 and 28 km in this region. In the entire region, a subsurface low velocity (Vs ~1.3 -2.5 km/s) layer with varying thickness (~4 and 7 km) represents Andaman flysch sediments. Middle Andaman is more heterogeneous and has one prominent discontinuity at ~15 km. By combining derived crustal seismic structure with seismicity and CMT fault plane solutions; it is observed that the seismicity in the region is mostly related to the subduction Indian slab beneath the Burmese plate. The Indian slab is subducting at variable dip and down to variable depth. Along the Andaman arc, the convergence is complex; with transitional (continental-oceanic) type in the north Andaman to oceanic type in the south Andaman. In the present study, we provide the information of the crustal seismic structure and subduction geometry in the Andaman region, however considering the complexity of the Andaman region, we find our dataset is limited and suggest a large scale multidisciplinary approach to examine and understand the 3-D structure, deep subduction and intrinsic processes beneath the Andaman and Nicobar Islands.
The presented results are the part of Ph.D. thesis of Santosh Mishra. Special thanks to Prof. Keith Priestley (Cambdridge, U.K.) for providing surface wave dispersion measurements and Prof. Bob Engdahl (Univ of Colorado, Boulder, U.S.A) for EHB dataset. The research was supported by Department of Science and Technology, New Delhi. We thank the Editor Prof. Saskia Goes and anonymous reviewers for their suggestions and comments for improving the quality of our manuscript.