Anisotropic tomography of eastern Tibet and its uncertainty from 1 hypocentral errors

region, we deem that the low-V p zone reflects hot and wet upwelling flow from the 32 deep asthenosphere, which ascends to the lower crust along the fault zone. At the NE 33 margin of the Tibetan Plateau, significant low-V p anomalies exist in the lower crust 34 and the FVDs are consistent with the motion direction of the Tibetan block revealed 35 by GPS observations. We think that lower crustal flow exists beneath NE Tibet, which 36 controls the plateau expansion toward the northeast. A low-V p anomaly appears at 30 37 km depth beneath the Sichuan Basin. However, as the weighting level increases, the 38 amplitude of this low-V p anomaly decreases by more than 6%, suggesting that this 39 low-V p anomaly has a larger uncertainty than the other features. 40

Since the Cenozoic, as a product of post-collision north-south convergence between India and Eurasia, the Tibetan Plateau has continued to uplift and grow eastward (e.g.Molnar & Tapponnier 1975;Armijo et al. 1986;Yin & Harrison 2000;Wang et al. 2019;Cheng et al. 2022).In eastern Tibet (Fig. 1), GPS velocity measurements have revealed a north-south convergence and significant eastward extrusion with a maximum rate of ~20 mm/yr (Wang & Shen 2020), suggesting that the eastern margin of the plateau is suffering an ongoing strong tectonic deformation at present.However, due to conflicts between various geological, geophysical and geochemical results, the mechanism of expansion and deformation of the Tibetan lithosphere is still under debate.Three end-member models have been proposed, including rigid block extrusion (e.g., Tapponnier et al. 1982Tapponnier et al. , 2001)), vertically coherent lithospheric deformation (e.g., England & Houseman 1986, 1988), and mid-lower crustal channel flow (e.g., Burchfiel et al. 1995Burchfiel et al. , 2008;;Clark & Royden 2000;Royden et al. 1997Royden et al. , 2008)).
In northeast Tibet, many seismic tomography studies have revealed low P and S wave velocity (low-V p , low-V s ) anomalies in the lower crust beneath the Qilian block (e.g., Cheng et al. 2014;Bao et al. 2015;Sun et al. 2019Sun et al. , 2021;;Li et al. 2022), which are interpreted as lower crustal channel flow.However, ambient-noise tomography shows that the low-V s anomalies are small and discontinuous, suggesting that lower crustal flow may be absent below the Qilian block (Zhao et al. 2021).Recent anisotropic tomography revealed significant NE-SW trending fast-velocity directions (FVDs) of V p azimuthal anisotropy in the lower crust beneath the eastern Qilian region (Cheng et al. 2016;Sun & Zhao 2020).In contrast, another result shows weak azimuthal anisotropy, suggesting that compressional deformation is the main cause of the Qilian orogeny (Zhou et al. 2023).Along the Longmenshan fault zone (LMFZ), most studies have consistently found a low-V p anomaly in the uppermost mantle beneath the plateau, but different explanations have been made about the formation mechanism of this anomaly.Some studies suggested that the low-V p anomaly is caused by horizontal lower crustal flow from the interior of the plateau (Liang & Song 2006;Zhang et al. 2010;Li et al. 2012;Wei et al. 2013;Lei et al. 2014), whereas other studies interpreted it to be a mixture of hot upwelling from the deep asthenosphere and the lower crustal flow (Wang et al. 2010a;Wang et al. 2010b;Lei & Zhao 2016).Beneath western Qinling, azimuthal anisotropic tomography revealed a continuous low-V p zone in the lower crust from eastern Tibet to western Qinling, which is inferred as a channel for lower crustal flow extruding eastward (Cheng et al. 2016;Li et al. 2022).However, other seismic studies did not find that feature (e.g.Ye et al. 2018;Sun & Zhao 2020;Fu & Xiao 2020;Li et al. 2021).The main reason of these contradictions is the inherent uncertainty and resolution of the seismic models, which could be influenced by several factors, such as the sparse seismic network, data errors, different model parameterization and inversion methods, etc.
Many previous studies have used the weighted least-squares (WLS) method to solve the tomographic inversion problem, in which the data are weighted unequally by adding a weighting matrix W to the objective function (e.g.Lees & Crosson 1989;Berryman 1989;Meyerholtz et al. 1989;Tarantola 2005;Rawlinson et al. 2014;Zhou et al. 2023).In these works, the estimated picking errors in observed arrival-time data are usually used to construct W to address the influence of measurement accuracy on the inversion result.However, this strategy does not consider the errors in theoretical travel times caused by the hypocentral mislocation and the approximation of forward theory, which may also have a significant influence on the inversion result, especially for local earthquake tomography (e.g.Kissling et al. 1994;Rawlinson et al. 2014;Rawlinson & Spakman 2016).To obtain a robust result, some researchers selected the relocated events by setting a set of fixed thresholds on errors of four hypocentral parameters, i.e. origin time, latitude, longitude, and focal depth (e.g.Zhao & Kanamori 1995;Liu et al. 2013;Zhao et al. 2022).However, there is still a lack of investigation on the uncertainties of a realistic tomographic model caused by various hypocentral error distributions.In eastern Tibet, the crustal velocity structure exhibits strong lateral and vertical heterogeneities (Wei et al. 2013;Laske et al. 2013;Han et al. 2021).Large uncertainties may exist in some relocated hypocentral parameters, especially in a region with significant structural heterogeneity, when a 1-D layered velocity model is adopted as the initial model for a tomographic inversion (Pavlis 1986;Moser et al. 1992;Lomax et al. 2000).If the initial model has some systematic errors, the true hypocentral errors would be further underestimated (Husen & Hardebeck 2010).
In this work, we apply the WLS method to conduct high-resolution V p azimuthal anisotropic tomography in eastern Tibet (Fig. 1) by using a newly collected  2 Method and Data

Azimuthal anisotropy tomography
Assuming hexagonal symmetric weak anisotropy (Backus 1965) and further considering a dipping ray path in a 3-D velocity model, P-wave velocity (V p ) in an anisotropic media could be expressed as (Eberhart-Phillips & Henderson 2004) where θ, λ, h are latitude, longitude and focal depth of an event, respectively.ΔV is isotropic V p perturbation; ΔA and ΔB are perturbations of two parameters of azimuthal anisotropy; p and q represent two sets of 3-D grid for isotropic V p and azimuthal anisotropy, respectively.T ij cal is theoretical travel time from an earthquake i to a station j after earthquake relocation, which is calculated with a 3-D ray-tracing method (Zhao et al. 1992).T ij arrival is the picked arrival time of the first P-wave, and T i ini is the corrected origin time of the event i after relocation.E ij represents higher-order terms of perturbations and errors in the data.The hypocentral parameters, V p perturbations and anisotropic parameters are determined simultaneously during the inversion.The initial hypocentral parameters are derived by using a linear relocation method (Engdahl & Lee 1976;Zhao et al. 1992), which is described in the following subsection.The fast-velocity direction (FVD) ψ and amplitude (α) of V p azimuthal anisotropy are expressed as follows: where V f and V s represent P-wave velocities in the fast and slow directions, respectively; V 0 is the average isotropic velocity.In the tomographic inversion, smoothing regularization is usually applied as a prior information (e.g.Liu et al. 2013;Jia et al. 2022;Jia & Zhao 2023), so the problem can be further expressed as where C is the coefficient matrix of eq. ( 2), and S is the finite-difference Laplacian operator in three dimensions (Lee & Crosson 1989;Jia & Zhao 2023), which is determined by constructing a trade-off curve between the root-mean-square (RMS) travel-time residual and the norm of a 3-D velocity model (Fig. S1).
According to the WLS method (Draper & Smith 1966;Lees & Crosson 1989), the problem described in eq. ( 5) is equivalent to where W is the weighting matrix, which controls the contribution of each datum to the final model and is assumed to be a diagonal matrix in this work.In fact, the W is a simplified treatment of the complete a prior data covariance matrix, which is usually difficult to recover in practice (Bodin et al. 2012;Duputel et al. 2014;Rawlinson et al. 2014).Therefore, the degree of simplification directly affects the accuracy of uncertainty analysis of the inversion results.Previous studies usually used the estimated picking errors in T arrival in eq. ( 2) to construct W (Lees & Crosson 1989;Rawlinson et al. 2014).However, the complete a prior data covariance matrix in seismic tomography is not only related to the data themselves (e.g.arrival times) but also assumptions in the model parameterization and forward theory (Rawlinson et al. 2014;Rawlinson & Spakman 2016;Freissler et al. 2020).

Hypocentral errors and weighting matrix
According to eqs ( 2) and ( 6) and the above descriptions, we can see intuitively that the errors in d are related to T ini and T cal , both of which are influenced by the event relocation results in local earthquake tomography.Hence, in the following, we address how to build W by taking into account hypocentral errors of each earthquake.
According to Engdahl & Lee (1976), each earthquake can be relocated as follows: where Δt is a travel-time residual before earthquake relocation; T ori is the origin time (T ini =T ori +ΔT ori ); θ, λ, h are latitude, longitude and focal depth of an event, respectively.The ordinary least-squares method can be used to solve eq. ( 8), and the standard error vector of the hypocentral parameters for each event can be expressed as (Flinn 1965;Draper & Smith 1966) where B is the coefficient matrix of eq. ( 8), s is the hypocentral parameter vector, and n is the number of arrival-time data used for the relocation.
Furthermore, to evaluate the influence of the initial model on the relocation results, in this work we conduct multiple sets of earthquake relocation by using different initial models.Then a total error vector e is calculated using all the e k derived from different initial models  where k is the number of earthquake relocations.We call e a composite error vector of the hypocentral parameters.Next, we use e to build a weighting coefficient in the range of [0, 1] for each local event, in order to control the relative contribution of different data to the inversion results.Hence, we need to normalize the four hypocentral parameters, for which the min-max normalization method is a simple and effective one (Han et al. 2001;Ciaburro 2018;Ciaburro et al. 2018;Kiran & Vasumathi 2020).As a result, the weighting coefficient w i in W for each event could be expressed as: where e qi is the qth element in e for ith earthquake, and Z is the total number of earthquakes.The parameter β is a predetermined real number greater than or equal to zero, which is called weighting level in this work.With increasing β, the data with larger hypocentral errors have smaller contribution to the results (Fig. 2), which means that velocity anomalies in the final tomographic model mainly controlled by these data will be more severely weakened than others.To investigate the influence of β on the inversion result, we conducted five sets of inversion by setting β= 0, 1, 2, 4, and 6.Obviously, it is just the traditional ordinary least-squares inversion when β= 0.
The inversion results are shown in the followings.Note that the weighting coefficient is assigned to each event, not to each travel-time datum.Thus, all the P-wave arrivals from the same event are weighted equally.
The weighted object function (eq.7) is resolved by using the L-BFGS-B algorithm (Byrd et al. 1992(Byrd et al. , 1995)), which has been applied by several recent tomographic studies and is proven to be a very efficient solver of tomographic Figure 2. Relationship between the weighting coefficients and errors of hypocentral parameters.The weighting coefficients are calculated using eq.11.The red, yellow, green, and blue dots represent four sets of results for all events by using β = 1, 2, 4, and 6, respectively.The thick black lines are logarithmic fitting curves for the 4 sets of results.

Data selection and earthquake relocation
The raw data set used in this work contains Pg, Pn and Sg wave arrival times recorded by a dense seismic network operated by the CEDC during 2008~2022 (Fig. 1).The total number of local earthquakes (M > 1.0) in the study region is 344,532.To obtain a reliable tomographic result, we rigorously checked and analyzed these data.Firstly, we divide the study region into a series of subareas with a size of 0.5° × 0.5°, and the top 20 largest events with more than 10 P-wave arrivals in each subarea are selected.
Then, multiple relocations are conducted by using different initial 1-D models, which are derived from the CRUST1.0model (Laske et al. 2013).V p and V s values in the upper and lower crust are changed synchronously between the minimum and maximum values in the study region provided by CRUST1.0 so as to simulate possible systematic errors in the initial model.Table 1 shows the five combinations we adopted, in which Set-5 is the average model derived from the 3-D CRUST1.0 model and is used for the final tomographic inversion.Then we use the five sets of relocated hypocentral parameters and their uncertainties to calculate the composite hypocentral errors for each event (eq.10).Two velocity discontinuities in the CRUST1.0model, i.e. the Conrad and Moho, are also considered to determine accurate 3-D ray paths.As a result, the total error distributions of all the events are obtained (Fig. 3).The composite hypocentral errors in this region are generally large, in particular, most events have a depth error of 5 to 15 km.
From the relocated events we select earthquakes for the tomographic inversion according to the following criteria: (1) the total number of observed arrival times for each event exceeds 10; (2) the composite errors of the origin time, focal depth, latitude, and longitude are less than 2.0 s, 20 km, 0.04° (~4.45 km) and 0.04° (~3.69 km), respectively (Fig. 3).The weighting coefficients are calculated for each event by using eqs ( 9)-(11) as mentioned above.As a result, 6084 local earthquakes are selected for the tomographic inversion, which generated 175,076 Pg, 76,810 Pn, and 79,182 Sg wave arrival times, leading to a good ray coverage of the study volume (Fig. 4).   3 Results

Grid setting and resolution tests
To investigate the robustness of the tomographic result, we performed a number of  S7).In the vertical direction, the grid meshes are set at depths of 5, 15, 30, 45, 60, 80 km.The initial synthetic model is constructed by setting ±3% V p perturbations and two anisotropic parameters at every two adjacent grid nodes.
Random errors following a Gaussian distribution with a standard deviation of 0.05 s are added to the synthetic travel times before the inversion.Smoothing parameters are also included, which are determined by a trade-off analysis as mentioned above.The CRT results show that the checkerboard model with a lateral grid interval of 0.5°×0.5°can be restored well for both isotropic V p and anisotropy in the whole crust beneath the study region.The best resolved areas are distributed along the boundaries between the Tibetan Plateau and surrounding blocks, which are discussed in the followings.
Note that the hypocentral locations are assumed to be accurate and not relocated in the CRT.Thus, the uncertainties in the final model arising from the hypocentral errors may be underestimated.Fortunately, these uncertainties can be examined through the WLS tomographic results in this work, which are shown and discussed in the following sections.
3.2 V p structure and azimuthal anisotropy In the upper crust (0-15 km depths), our results show a complex pattern of isotropic V p and anisotropy.The average amplitude of the V p anomalies is small and further decreases as the weighting level increases (Figs S8-S11).Some strong anisotropies occur at 5 km depth in the southwestern part of the study region.
However, they become unclear in the weighted tomographic models when the weighting level is large (Figs S10 and S11).Hence, the results in the upper crust are not discussed further.

Variations of weighted tomographic models
Fig. 6 shows four models of V p anisotropic tomography obtained with the WLS method with different weighting levels.All these models show significant heterogeneity of isotropic V p and anisotropy along the eastern margin of the Tibetan Plateau.As the weighting level increases, the amplitude of V p anomalies decreases to varying degrees.While most of them exhibit a decrease of ~3%, some areas show a larger reduction of up to 6% (Fig. S12).A significant feature is a narrow low-V p anomaly surrounded by high-V p anomalies and fault parallel FVDs in the lower crust beneath the western side of the LMFZ.The maximum amplitude of the low-V p anomaly decreases by only 3% and remains above 4% as the weighting level Beneath the Sichuan Basin, a prominent low-V p anomaly appears at 30 km depth in the tomographic model (Fig. 5c).However, as the weighting level increases (Figs 6c1 and d1), the amplitude of this anomaly is reduced by more than 6% (Fig. S12).The relatively large decrease in its amplitude suggests that this low-V p anomaly has a larger uncertainty than the other anomalies.

Reliability of the weighting method
To know whether our weighting method can accurately detect the robust components in the final tomographic model, we also performed an inverse-weighted test by adopting a weighting level of 6, in which an inverse-weighting coefficient is set as w i ' =1-w i for each event i in eq. ( 11).Theoretically speaking, in this inverse-weighted test, the data with large hypocentral errors would have larger contributions to the final tomographic result, that is, more artificial anomalies would be emphasized whereas robust anomalies are weakened.The test results (Figs S11 and S13) show that the significant decrease in the amplitude of the low-V p anomaly and anisotropy beneath the LMFZ could be identified in the inverse-weighted tomography, which has been regarded as a robust feature as mentioned above.Instead, some "new" strong anomalies are more pronounced, such as strong anisotropies at 5-15 km depths beneath southwest of the study region (Fig. S13a) and significant anomalies at 80 km depth (Fig. S13f).However, these "anomalies" are clearly artifacts, due to the poor ray coverage in those areas (Fig. 4).
To distinguish the influences of the smoothing parameter and weighting level on the final tomographic result, we further performed a series of synthetic tests to examine the response of two "typical" V p anomalies in our model (i.e. two low-V p zones under the western side of the LMFZ and the Sichuan Basin) to varying levels of smoothing regularization and weighting.The test results (Figs S14 and S15) show that as the smoothing increases, the maximum amplitude decreases similarly for both low-V p anomalies.However, as the weighting level increases, the amplitude of the low-V p anomaly beneath the Sichuan Basin decreases more significantly than that beneath the LMFZ, particularly for the optimal smoothing parameters (Fig. S15c).All these test results indicate that it is feasible to identify the robustness of anomalies in our WLS tomographic results.As shown in Fig. 7, a significant low-V p anomaly (LV2) is surrounded by high-V p anomalies in the lower crust and uppermost mantle along the western side of the LMFZ.FVDs of azimuthal anisotropy are nearly NE-SW in the lower crust, which are parallel to the strike of the Longmenshan fault.These features are relatively stable for different weighting levels as mentioned above, suggesting that it is a robust feature in our tomographic model.Several previous studies also found a similar pattern in the lower crust beneath this area (Liang & Song 2006;Li et al. 2012;Wei et al. 2013;Lei et al. 2014) and it was interpreted as horizontal lower crustal flow (Burchfiel et al. 1995;Royden et al. 1997;Clark & Royden 2000), which comes from the interior or the northern parts of the plateau and is obstructed by the rigid Sichuan Basin block (Liang & Song 2006).Another explanation from teleseismic tomography is that the low-V p anomaly beneath eastern Tibet (including the LMFZ) may be also related to hot and wet upwelling flow in the upper mantle due to the deep subduction and dehydration of the India slab (Wang et al. 2010;Lei & Zhao 2016).In our model, the anomaly LV2 extends from the lower crust to the uppermost mantle (Figs 7 and 8), being consistent with the teleseismic tomography results.Furthermore, we find that the low-V p anomaly is surrounded by circular high-V p anomalies in the lower crust and does not extend too far into the interior of the plateau (Fig. 8).In addition, strong anisotropy with fault parallel FVDs and weak anisotropy are found in the low-V p and circular high-V p zones, respectively.According to these features and the previous results, we deem that LV2 beneath the western side of the LMFZ is mainly caused by hot upwelling flow from the deep asthenosphere.When the hot materials rise into the lower crust, they migrate along the fault zone, resulting in the fault parallel FVDs.
The strong vertical movement in the deep area could also cause a large uplift but small horizontal shortening in the upper crust of the plateau along the western side of the LMFZ, which has been revealed by the GPS observations (Wang & Shen 2020) and regional geologic and geodetic studies (Shen et al. 2005;Meade 2007).).These features are also visible in the WLS tomographic models with high weighting levels (Fig. 6).Hence, we deem that the lower crustal flow does exist beneath NE Tibet, which controls the northeastward expansion of the plateau.

Low-V anomaly beneath the Sichuan Basin
Beneath the Sichuan Basin, our results show a low-V p anomaly at 30 km depth, which is similar to some previous results (Lei & Zhao 2009;Wang et al. 2010b) with a slightly larger range at 20-30 km depths.This low-V p anomaly in the lower crust has been rarely discussed and regarded as an unreliable feature by some previous studies (e.g.Lei & Zhao 2009).Numerical simulations also showed that a rigid lower crust beneath the Sichuan Basin could accommodate well with the shape and temporal evolution of mountain ranges in SE Tibet (Penney & Copley 2020).Our WLS tomographic results further show that as the weighting level increases, the amplitude of this low-V p anomaly decreases by more than 6%, which is greater than that of most other anomalies (Figs 9,10 and S12).Therefore, we think that this low-V p anomaly in the lower crust beneath the Sichuan Basin has a large uncertainty.Nevertheless, because there is a strong correlation between the low-V p anomaly and the gas-oil field distribution in this region (Wang et al. 2020b), it is highly possible that significant tectonic processes occur in the lower crust, potentially linked to the formation of the gas-oil field.(1) A robust significant low-V p anomaly and fault parallel FVDs (fast-velocity directions) are revealed in the lower crust beneath the Longmenshan fault zone, which are surrounded by circular high-V p anomalies.This feature suggests that hot upwelling flow from the deep asthenosphere occurs beneath the western side of the Longmenshan fault zone, which has controlled uplift of the upper crust in this region.
(2) In NE Tibet, low-V p anomalies and NE-SW FVDs exist in the lower crust (3) A low-V p anomaly is revealed at 30 km depth beneath the Sichuan Basin, which may reflect significant tectonic processes occurring in the lower crust, potentially linked to the formation of the gas-oil field in the region.However, this low-V p anomaly has a larger uncertainty and so should be investigated further by

Figure 3 .O
Figure 3. Histograms of errors of four hypocentral parameters after the earthquake 272 relocation.The errors are estimated from multiple relocation results.The black dashed lines 273show the cut-off uncertainties adopted for the event selection for the tomographic inversion.274 275

Figure 4 .
Figure 4. Distribution of 6084 local earthquakes (colour dots) selected for conducting tomography.Gray lines denote ray paths.The dot colours denote focal depths, whose scale is shown at the lower-left corner in (a).Blue triangles in (a) denote seismic stations.The black dashed line in (b, c) denotes the Moho discontinuity from the CRUST1.0model (Laske et al. 2013) along N-S and E-W profiles passing through the center of the study region.

Fig. 5
Fig.5shows the final tomographic results obtained after three iterations of the WLS tomographic inversion.A significant narrow low-V p anomaly exists at 30-45 km depths beneath the western side of the LMFZ, which has also been found by several previous studies of local tomography(Liang & Song 2006; Zhang et al. 2010;Li et al. 2012;Wei et al. 2013;Lei et al. 2014).Our model further shows that this low-V p anomaly is surrounded by some high-V p anomalies and extends down to the uppermost mantle (Figs5d-f).Anisotropic results reveal significant NE-SW FVDs in this lower crustal low-V p zone (Fig.5d), which are consistent with the strike of the LMFZ.At the NE margin of the Tibetan Plateau, low-V p anomalies and strong azimuthal anisotropy exist in the lower crust beneath the Qilian block, and the FVDs are NE-SW and E-W in the western and eastern parts of the fault zone, respectively (Figs5c and d).Another prominent low-V p anomaly appears at 45 km depth beneath southeast Tibet and the FVDs are nearly NW-SE in the low-V p zone, which is consistent with the direction of absolute plate motion observed by GPS measurements(Wang & Shen 2020).

OFigure 5 .
Figure 5. V p anisotropic tomography determined by the WLS inversion with a weighting level of 2. The layer depth is shown at the lower-right corner of each map.Red and blue colours denote low and high isotropic V p perturbations, respectively, whose scale is shown below (e).The orientation and length of the black bars denote FVDs and amplitudes, respectively, of V p azimuthal anisotropy.The anisotropic amplitude scale is shown below (e).The anisotropic anomalies with amplitudes <0.5% are not shown.The gray lines denote active faults in the study region.

Figure 6 .O
Figure 6.V p anisotropic tomography determined by the WLS inversion with a weighting level of (a) 1, (b) 2, (c) 4, and (d) 6.Other labels are the same as those in Fig. 5.

OFigure 7 .O
Figure 7. V p anisotropic tomography determined by the WLS inversion with a weighting level of 2. The layer depth is shown at the lower-right corner of each map.Red dashed contour lines mark three significant anomalies discussed in the text.Black arrows denote possible directions of the lower crustal flow.Other labels are the same as those in Fig. 5.

Figure 8 .O
Figure 8. Vertical cross-sections of V p tomography along profiles D1-D3 determined by the WLS inversion with a weighting level of 2. The surface topography is shown above each panel.Short green lines denote the approximate location of the Longmenshan fault (modified from Wu, 2014).The black dashed line denotes the Moho discontinuity.The red star denotes the hypocenter of the 2008 Wenchuan earthquake (Ms 8.0).Gray dots on the map denote large earthquakes (M > 5.9).The V p perturbation scale is shown below the map.SCB, the Sichuan Basin.

Figure 9 .OFigure 10 .
Figure 9. Map views of V p tomography at 30 km depth derived from WLS inversions with different weighting levels as shown at the upper-left corner of each map.The black thick contour denotes the approximate boundary of the Sichuan Basin (SCB).Gray thin lines denote active faults.The white line shows the profile D2-D2' in Fig. 10.