Seismic anisotropy inferred from direct S-waves derived splitting measurements and its geodynamic implications beneath southeastern Tibetan Plateau

The present study deals with detecting seismic anisotropy parameters beneath southeastern Tibet near Namche Barwa Mountain using splitting of the direct S-waves. We employed the reference station technique to remove the effects of source side anisotropy. Seismic anisotropy parameters, splitting time delay and fast polarization directions were estimated through analyses on a total of 501 splitting measurements obtained from direct-S waves from 25 earthquakes (>5.5 magnitude) that were recorded at 42 stations of Namchebarwa seismic network. We observed a large variation in time delays ranging 5 from 0.64 to 1.68 s but in most cases it is more than 1 s, which suggests for a highly anisotropic lithospheric mantle in the region. A comparison between direct Sand SKS-derived splitting parameters generally shows a close similarity although some discrepancies exist where null or negligible anisotropy is reported earlier using SKS. The seismic stations with hitherto null or negligible anisotropy are now supplemented with new measurements with clear anisotropic signatures. Our analyses indicate a sharp change in lateral variations of fast polarization directions (FPDs) from consistent ENE-SSW or E-W to NW10 SE direction at the southeastern edge of Tibet. Comparison of the FPDs with global positioning system (GPS) measurements, absolute plate motion (APM) directions and surface geological features signify that the observed anisotropy and hence inferred deformation patterns are not only due to asthenospheric dynamics but it is a combination of lithospheric deformation and sublithospheric (asthenospheric) mantle dynamics. Splitting measurement using direct-S waves proves their utility to supplement the anisotropic measurements in the study region and fills the missing links that remain rather illusive due to lack of SKS 15 measurements.

neglected by searching for optimum splitting parameters resulting in the minimum difference between the S-wave signals corrected for receiver-side anisotropy beneath the reference and target stations. Signals used for that comparison are that of the reference station previously corrected for known reference receiver-side anisotropy and of the target station whose receiverside splitting parameters are desired to be estimated. In this technique we utilize seismic anisotropy parameters, which were previously inferred from the SKS measurements by Sol et al. (2007) as the reference knowledge of the receiver-side anisotropy 5 beneath the reference station. The Reference Station Technique has been earlier successfully tested through both synthetic and observed data collected along the northeastern and southwestern parts of the Tibetan Plateau (Eken and Tilmann, 2014;Singh et al., 2016). Present study focuses on the southeast part of the Tibet near Namche Barwa (Figure 1). Study region is located between and around the Indus-Tsangpo Suture Zone (ITSZ) and Bangong-Nuijiang Suture Zone (BNSZ). Our major motivation is to find out the S-wave derived seismic anisotropic parameters that may have a potential link to current tectonic 2 Tectonics of the region 15 The formation of Tibetan plateau and the Himalayan mountain belts are due to collision and post collision processes of the Indian and Eurasian plates since 50 million years ago (Argand, 1924;Garzanti and Van Haver, 1988;Molnar and Tapponnier, 1975;Yin and Harrison, 2000;Royden et al., 2008). Underthrusting of the Indian lithosphere beneath the Eurasian lithosphere is the main responsible reason for the formation of the Himalayan and Karakorum ranges (Nelson et al., 1996;Kumar et al., 2006;Tseng et al., 2009) along with the formation of the Central Tibetan region (Argand, 1924;Nelson et al., 1996;Li et al., 20 2008). Understanding the development of the northern and eastern Tibetan plateau region are a big challenge (Karplus et al., 2011;Royden et al., 2008). McKenzie and Priestley (2008) discuss the development of the northern Tibetan lithosphere as an accreted one. Royden et al. (2008) argued that the Tibetan plateau has been evolved due to the subduction of Indian lithosphere beneath Eurasia, which is also responsible for the thickening of Tibetan crust and afterward extrusion of the Tibetan lithosphere towards east. 25 Various models have been developed regarding the deformation of Tibet (Royden et al., 1997;Molnar and Tapponnier, 1975;Houseman and England, 1986Tapponnier et al., 1982Tapponnier et al., , 2001Shen et al., 2001;Holt et al., 1995Holt et al., , 2000Replumaz and Tapponnier, 2003;Flesch et al., 2001) but no single model can explain the whole scenario. The debate is open for understanding the deformation mechanism, crust and mantle deformation patterns and ongoing geodynamics. Multistage subduction of Tethys oceanic plate and Indian plate below the Eurasian plate resulting in highly heterogeneous and anisotropic 30 lithosphere is also discussed in recent works (Jagoutz et al., 2015;Van Hinsbergen et al., 2012). The Tibetan and Himalayan region is mainly dominated by thrust and strike slip faulting (Figure 1). Suture zones are extended in the E-W direction and take a sharp turn around eastern Himalayan Syntaxis (EHS). Strike slip faulting become more dominant to east of the EHS. Figure 1 signifies that the eastern portion of the subducting Indian plate are found up to the EHS (León Soto et al., 2012) where the structural and topographical features takes a sharp trend in N-S striking from nearly E-W striking.

Data and Method
In this study we used a total of 5285 waveforms with the direct-S waves extracted from 161 teleseismic events with magnitudes 5.5 within an epicentral distance range from 30 to 90 ( Figure 2). These teleseismic earthquakes used in this study are 5 recorded at 47 seismic stations of the XE network, which was operated in between (Sol et al., 2007 (Figure 3). Out of 47 stations we used only those 36 seismic stations where we have knowledge of seismic anisotropy (SKS splitting parameters measured by Sol et al. (2007) to be used when correcting for receiver-side anisotropy beneath the reference stations. Prior to the data analysis, we removed the instrument response from the original seismograms. At the stage of the preprocessing, a band-pass filter between 0.03 and 0.2 Hz was performed to enhance S signals and resampled the seismograms at 20 samples 10 per seconds. Signals with any possible contaminations with other phases such as ScS, SKKS, and SKS are omitted from the analysis. We selected only those waveforms with a good signal to noise ratio (SNR) on transverse and radial components for further analysis. Waveform selection was achieved by performing a visual inspection manually that allowed only 40% of the direct-S waveforms. All the selected good quality waveforms show the characteristic of the splitting with clear energy on the transverse component and elliptically polarized particle motion. We started data analysis by determining station pairs over 15 the entire area. Out of 355649 station pairs from 161 teleseismic events, we were left with 22816 station pairs within the interstation spacing of 300 km. We can use the Reference Station Technique (Eken and Tilmann, 2014) at a station pair only when four horizontal components of the same event are available at reference and target station. To minimize the effects of the coda waves and converted phases, we used a 45s time window starting 15s before the theoretical onset of the direct-S waves on the basis of the IASP91 1-D radial earth velocity model of Kennett (1991). 20 The approach used in the present study avoids the source side anisotropy by minimizing the misfit function between the corrected seismic waveforms at the reference and target stations. At the first stage an inverse splitting operator that consists of a backward angular rotation with two horizontal components and a time shift and the reversal of the backward rotation is employed to correct the reference station for known receiver-side anisotropy (generally inferred from SKS splitting analyses) when estimating of direct S-derived individual splitting parameters. Following the correction of the reference station S-signals 25 at the target stations were corrected for arbitrary splitting parameters ( and t) in a grid search manner. Corrected S-wave signals at reference and target stations are compared to each other for each pair of arbitrary splitting parameters. Such comparison also allows for arbitrary time shifts and amplitude corrections to account for the lateral heterogeneities and differences in site response between these stations by optimizing the time shift (4t) and amplitude factor (a). In the end, we selected final splitting parameters (FPD and TD) beneath the target station by minimizing the misfit function that simply represents the 30 difference between corrected reference and target station traces.
The Reference Station Technique relies on two most important underlying assumptions: i) the distance between receiver and target stations is small compared to epicentral distance by virtue of that the ray path at these two stations can be considered Add citation to Eken and TIlmann, 2014 These are not arbitrary. These are best fit from a grid search You present results for station pairs, but what if you use a different reference at a specific target? equivalent in deeper mantle part and near the source-side region, ii) waveform differences between the receiver and target stations are only due to difference in anisotropic structure after correcting any waveform differences in time and amplitude presumably due to the lateral heterogeneities and differences in site response between these stations. Any potential difference between the thickness of the crust and sedimentary layers will also cause the timing and amplitude of converted phase but Eken and Tilmann (2014) earlier showed numerically its influence on expected splitting parameters would be negligible. The 5 optimum interstation distance less than 300 km was previously validated for regional arrays by several application of the technique (e.g., Eken and Tilmann, 2014;Singh et al., 2016). During the application of the technique we let and t to vary from 0 to 180 with an increment of 1 and 0 to 3 s with an increment of 0.05 s, respectively. We performed the inverse F-test error analysis of Silver and Chan (1991) for uncertainty estimates of obtained splitting parameters. In this process, we checked reliability of the individual splitting parameter by comparing variation in the residual energy distribution away 10 from the minimum with the variation according to preset confidence level. An example of the basic steps of the Reference Station Technique can be found in Figure 4 for target stations ES01 and ES03 respectively. Figure 5 presents the examples of the obtained splitting parameters with null at target stations ES19 and ES33, where we get null measurements by Sol et al. (2007). Null splitting may arise due to the two reasons: i) if the incoming polarization direction below an anisotropic layer becomes parallel to the fast or slow axis; ii) if the region is isotropic in nature, due to complex anisotropy (e.g., Saltzer 15 et al., 2000;Wüstefeld and Bokelmann, 2007). Following Eken and Tilmann (2014) and Singh et al. (2016) we also benefited from the F-test with the null-split rejection criteria to be able to avoid the contamination of the null measurement for the good splitting measurements. In this process, we theoretically calculate the residual energy under the assumption of the null measurements and compare this with the observed residual energy at the minimum. Figure 6 shows two examples of the null splitting measurements at target stations ES29 and ES37 respectively. To ensure the stability of the results we have performed a 20 stepwise quality assessment criterion before calculating the average splitting parameters at each station. To achieve this aim we considered only those waveform pairs that have (i) normalized residual energy (4E) smaller than 0.5; (ii) amplitude correction factor parameter (a) in between 0.4 and 0.6; and (iii) 95% confidence level for null splitting rejection. We rejected the waveform pairs that have an FPD error greater than 25 and delay time error greater than half of the delay time itself. After these quality assessments we were left with only 3231 waveform pairs. At this stage of the processing we again performed 25 another visual inspection manually to enhance the quality our estimates that finally allowed only 501 waveform pairs. These waveform pairs were extracted from only 25 teleseismic events and used to calculate the average splitting parameter at each station. We apply the Von Mises approach (Cochran et al., 2003) to calculate the circular mean at each target stations for and an arithmetic mean is used for t.

30
We present here 501 splitting measurements observed for 42 seismic stations of the XE network between the years 2003 and 2004. Angular average of individual direct S-derived splitting parameters ( s and t s ) at each station is given in Table 1. Stationaveraged splitting parameters usually reflect that the significant anisotropy is present in the region. We observe a large variation in delay times ranging from 0.64 s to 1.68 s. By using the reference station technique, we revealed significant anisotropy in Because so much data is thrown out, you really don't improve sampling that much "Significant" requires a statistical test. Perhaps you mean non-neglible? or reference to SKS splitting amplitudes areas where negligible anisotropy was reported earlier depending on only the SKS wave splitting measurements. For example at station ES31 direct-S waves provided relatively large time delay (1.23 s) although SKS splitting analysis performed by Sol et al. (2007) earlier resulted in a much smaller time delay about 0.3 s. In general, we observe the NE-SW to E-W trend in s before the edge margin of the southeastern Tibetan region. A consistent change in variation of s is observed further east when orientations take a sharp change from ENE-SSW or E-W to NW-SE direction ( Figure 7). We find the significant 5 anisotropy ( 0.64 s) at seismic stations ES19, ES20, ES22, ES32, ES33, ES34, ES42 and ES45, where previously hitherto null or negligible anisotropy were obtained with splitting of the SKS waves (Sol et al., 2007). This could stem from multi-layered anisotropic orientations or insufficient amount of SKS-derived splitting measurements. This rule out the consideration of the isotropic nature of the Indian lithosphere and signifies the complex 3D nature of more complex deformation pattern of the EHS.
The scatter plot in Figure 8 exhibits a comparison between estimated splitting parameters ( s and t s ) from the analysis of 10 direct-S wave and previous SKS splitting measurements (Sol et al., 2007). The overall trend of the obtained splitting parameters ( s ) is consistent with the previous SKS measurement except at some places while generally observe the larger delay time of the S wave as compared to SKS wave.
We combined our splitting study with the existing geodetic measurements (GPS velocity vector, APM velocity vector). We

15
The APM velocities were calculated via a web-based plate motion calculator (https://www.unavco.org/software/geodeticutilities/plate-motion-calculator/plate-motion-calculator.html) that is based on an integrated global plate motion model (GSRM v1.2) originally developed by Kreemer et al. (2003). Figure 9 simply presents a complex map, which reflects the correlative analysis of splitting parameters by using direct S-phases, APM directions, and GPS measurements in our targeted region. The observation suggests that the observed anisotropy is not only due to lithospheric deformation or due to asthenospheric dynamics 20 at the base of the lithosphere, but it is a combined effect of both

Origin of anisotropy in the Southeastern Tibetan region
Our resultant splitting measurements varying in a range of t s from 0.64 to 1.68 s suggest presence of very complex deformation types in the region. The fast polarization directions are rather consistent and matching well with the surface geology similar to 25 those observed from the SKS phases in Sol et al. (2007). The fast axis closely follow the follow the strike of the major sutures like BNSZ and ITSZ and surface strain fields as observed through GPS and are under the influence of bending at Eastern Himalayan Syntaxis (Figure 7). Vertically coherent deformation of the crust and upper mantle are associated to the seismic anisotropy directions parallel to the surface geologic features such as faults etc. (e.g., Savage, 1999;Flesch et al., 2005) and earlier has been invoked as a possibility to explain the anisotropic character in the eastern and northeastern Tibet (Holt et al.,30 2000; León Soto et al., 2012;Eken et al., 2013;Eken and Tilmann, 2014). Mechanically coupled crust-mantle interface or coherent deformation of whole lithosphere with similar boundary conditions are used to exaplin the observed anisotropy in eastern Tibet (Holt et al., 2000;Sol et al., 2007). In absence of any compelling evidence for crust-mantle coupling, we argued  (Flesch et al., 2005;Sol et al., 2007;Holt et al., 2000).
Observed larger time delays (>1s), in this study, reflect a highly anisotropic region with similar deformation patterns at depths. Presence of a more complex anisotropic structure (e.g. double-layer) with different orientation in fast axis at various depths may result in smaller delay times (Saltzer et al., 2000). In the western Himalayan region, Vinnik et al. (2007) observed 5 the different orientation of the fast direction of seismic azimuthal anisotropy from 60 at 80-160 km to 150 at 160-220 km depth region by using the joint inversion of the SKS particle motion and P receiver functions. This provides an argument to explain the null or negligible anisotropy as reported from the same region using SKS phases (Sandvol et al., 1994). Smaller time delays in Nepal Himalaya and Sikkim Himalaya were attributed to the combined effect of shear at the base of lithosphere due to APM related strain of Indian plate and ductile flow along the collision front due to compression, with possibly completely 10 different orientations . Null measurements were reported at few stations used in the present study by Sol et al. (2007) using SKS/SKKS phases, although no specific reasons are provided. The transition between deformation types at the boundaries of Indian and Eurasian lithospheric plates was considered to be the main reason for observed null or negligible anisotropy further west beneath southern Tibetan Plateau Chen and Ozalaybey, 1998;Barruol and Hoffmann, 1999;Zhao et al., 2014). The lack of anisotropy beneath southern Tibet was mainly explained by the isotropic 15 nature of Indian tectonic plate or lack in ability of SKS phases to sample the anisotropy due to sub-vertical mantle shear strain field created due to downwelling Indian lithosphere Sandvol et al., 1997). However the isotropic nature of Indian lithosphere was contradicted in various studies (Singh et al., 2006Kumar and Singh, 2008) and significant anisotropy is reported beneath Tibet in region of null (Gao and Liu, 2009;Singh et al., 2016). Indian lithosphere is considered to be more isotropic in nature to explain the null or negligible anisotropy beneath southern Tibet (Chen and Ozalaybey, 1998;20 Barruol and Hoffmann, 1999;Chen et al., 2010). Taking forward the same argument to explain the null measurements at few specific stations beneath Namchebarwa Mountains requires the northern extent of the Indian lithospheric mantle north of ITSZ.
The recent tomographic studies (Griot et al., 1998;Huang et al., 2003;Zhou and Murphy, 2005;Yao et al., 2008;Priestley et al., 2006;Singh et al., 2014;Pandey et al., 2014) suggest that in the western Tibetan side where the N-S extension is less, the Indian lithosphere is supposed to be up to the Jinsa River Suture Zone (JRSZ) (Zhao et al., 2010), while in the eastern Tibet side 25 the Indian lithosphere is extending up to the ITCZ Zhao et al., 2010). Combined study using seismic anisotropy and Bouguer gravity anomalies placed the Indian mantle front up to 33 N in central Tibet . In the present segment of the Himalaya-Tibet collision zone northern limit of the Indian lithospheric mantle does not seem to extend beyond ITSZ . The lack of anisotropy reported using SKS/SKKS phases (Sol et al., 2007) at few seismic stations may be due to insufficient measurements. By adding considerably large amount of measurements from direct-S waves we observe

Comparison between direct S and SKS derived splitting parameters
The scatter plot in Figure 8 exhibits a comparison between estimated splitting parameters ( s and t s ) from the analysis of direct-S wave and previous SKS splitting measurements (Sol et al., 2007). The obtained splitting parameters ( s ) from the direct-S wave measurements depict that it is consistent with the previous SKS measurements except at some places. Overall consistency between splitting parameters inferred from SKS and direct-S waves is most likely due to the fact that both are 20 exposed to the same sort of large-scale anisotropic structures. Large differences between SKS-and S-derived t s that appear as a move-out in the scatter plot is occur due to the two reason: i) longer travelling paths of the direct-S waves as compared to the SKS exposed for the same type of the of large scale anisotropy and, ii) large number of the events from different azimuths contributing to the direct-S wave measurements as compared to the SKS. The accuracy of the measurements relies on more measurements and direct-S waves splitting measurements supplement the splitting measurements from the region. Comparative 25 studies of the splitting parameters obtained from the direct-S waves are in close agreement to the SKS splitting measurements as observed in the Himalaya-Tibet collision zone (e.g., McNamara et al., 1994;Singh et al., 2016;Huang et al., 2011;Eken et al., 2013) and in the Indian shield (Saikia et al., 2010)

This is not well supported by event distribution
This is a key point, but I misread it at first southeastern or eastern Tibetan region. We observe a sharp transition in the spatial distribution of s due to change in the fast polarization direction from the nearly east west direction to nearly NE to SW or NNE to SSW near the edge margin of the southeastern Tibet. Structural and topographical feature such as major suture zones, mountain belts, etc. tends to rotate around the EHS from nearly E-W or ENE-WSW to N-S or NE-SW direction (Hallet and Molnar, 2001;Booth et al., 2004).
The observed FPDs and GPS velocity vectors follow nearly the similar trend. The absolute plate motion (APM) directions are 5 well consistent with the present ongoing asthenospheric flow (Vinnik et al., 1992(Vinnik et al., , 1995Vinnik and Montagner, 1996). The FPDs overall do not match with the APM velocity direction of Indian plate with respect to Eurasian plate or with the APM velocity of the Eurasian plate and Indian plate with respect to fixed frame or no-net rotation frame. The discrepancy between the FPD and APM may indicate that the obtained splitting and hence the anisotropic behavior of the study area is not only due to asthenospheric dynamics but it is a combined effect of the lithospheric deformation and asthenospheric dynamics. The by Royden et al. (1997Royden et al. ( , 2008. 15 The present day GPS measurement does not reveal the deformation of the whole crust but could be rather associated to the deformation of the shallow crust (Chen et al., 2013). and surficial features. According to Sol et al. (2007)Sol et al. (2007 coupling between the crustal and mantle material for the entire southeastern Tibetan region existed as similarly observed in the northeast Tibetan plateau (e.g., Eken et al., 2013;Eken and Tilmann, 2014). The seismic anisotropy directions that were previously obtained from the inversion of receiver functions, however, do not suggest vertically coherent deformation of the crust (e.g., Ozacar and Zandt, 2004;Sherrington et al., 2004). 25 One possibility for this mismatch could be due to the nature of receiver functions analyses that are prone to be sensitive to rather fossilized fabrics representing complex deformation history of the crust in Tibet although crust is deformed coherently.
On the basis of driving forces, two kinematic models have been proposed to explain the coupling-decoupling of the crust and lithospheric mantle. First one is simple asthenospheric model (Richardson, 1992) proposed to explain the decoupling of the crust and mantle by intruding the weak mechanical layer such as asthenosphere and other in the crust. Whenever a weak

Do you mean NW-SE?
Based on what evidence?
Why are RFs only sensitive to fossil structures? (hint, they aren't!) of the buoyancy forces from crust to the mantle. This model requires a rigid lower part of the crust. However recent the low shear wave velocity anomalies resolved in various tomographic studies recently (Huang et al., 2002(Huang et al., , 2009Hung et al., 2010;Yao et al., 2008;Li et al., 2009) have indicated a weak layer in the deeper region of the crust beneath the SE Tibetan region by contradicts this. It is noteworthy to mention that we avoid making any comment on the possible linkage between deformation and coupling of crust and underlying lithospheric mantle by only using splitting parameters inferred from direct-S waves and 5 geodetic measurements.

Conclusions
Our splitting observations using direct-S waves add new constraints in understanding the deformation pattern and their causes in the southeastern Tibetan region near Namche Barwa. We can list main concluding remarks from the present study as follows: 1. The observed splitting analysis suggests for a highly deformed crust and lithospheric mantle.  5 Code availability 20 The multisplit C++ code used for carrying out splitting measurements of the direct S waveforms is available with a General Public License (GPL) license at http://github.com/ftilmann/multisplit.

Data availability
The waveform data used in this study are downloaded from the Incorporated Research Institutions for Seismology Data Management Center (IRIS-DMC) data archive system Namche Barwa Tibet network code XE (2003XE ( -2004. 25 Acknowledgements. The IRIS data management centre and Project team of Nameche Barwa Seismic network (PASSCAL) are gratefully acknowledged for making the seismic data available. This work has been performed under the ISIRD project (ASI) of IIT Kharagpur. T.
Eken appreciates the Alexander-von-Humboldt (AvH) foundation for the postdoctoral scholarship as well as equipment subsidy that have contributed greatly at different stages of the development of the main software (multisplit) used in the present study. igure 3. Earlier SKS/SKKS measurements in the area (Sol et al., 2007;Wüstefeld et al., 2008). The length of line at each seismic station represents the delay time ( t s ) and their orientation denotes the fast polarization direction ( s ). For clarity, the seismic stations used in the present study are shown by yellow filled rectangles along with the SKS splitting measurements of Sol et al. (2007). Seismic stations where null or negligible anisotropy is reported in earlier studies (see Wüstefeld et al., 2008) are shown by gray filled rectangles. igure 9. Lateral variation of anisotropic, geodetic, absolute plate motion data shown over topographic and tectonic features of the study area.
(a) Map view comparison between the splitting measurement inferred from direct-S waves (this study shown by red bars) and SKS phase (Sol et al., 2007, in blue bars) measurements. (b) The Global Positioning System (GPS) velocity (mm/yr) vectors (black arrows) around SE Tibetan region calculated with respect to the South China reference frame and stable Eurasia. GPS data is compiled from several studies including Chen et al. (2000); Zhang et al. (2004); Shen et al. (2005); Sol et al. (2007). (c) Absolute plate motion (APM) velocities calculated from https://www.unavco.org/software/geodetic-utilities/plate-motion-calculator/platemotion-calculator.html by using GSRM v1.2 (2004) model of Kreemer et al. (2003) Note that arrows in orange, green, and purple represent the APM velocities of the Eurasian plate in no net rotation frame, of the Indian plate with respect to the Eurasian plate, and the motion of the Indian plate in no net rotation frame, respectively.