Imaging the lithosphere beneath NE Tibet : teleseismic P and S body wave tomography incorporating surface wave starting models

Ceri Nunn,1 Steven W. Roecker,2 Frederik J. Tilmann,3,4 Keith F. Priestley,1 Ross Heyburn,5 Eric A. Sandvol,6 James F. Ni,7 Yongshun John Chen,8 Wenjin Zhao9 and the INDEPTH IV Team 1Department of Earth Sciences, University of Cambridge, Cambridge, UK. E-mail: cn277@cam.ac.uk 2Department of Earth and Environmental Sciences, Rensselaer Polytechnic Institute, Troy, NY 12180, USA 3Helmholtz-Centre Potsdam, GFZ German Research Centre for Geosciences, Potsdam, Germany 4Department of Earth Sciences, Freie Universität Berlin, Berlin, Germany 5AWE Blacknest, Reading, UK 6University of Missouri, Columbia, MO 65211, USA 7New Mexico State University, Las Cruces, NM 88003, USA 8Peking University, Beijing, China 9Chinese Academy of Geological Sciences, Beijing, China


I N T RO D U C T I O N
The high (∼5 km) elevation of the Tibetan Plateau is a consequence of the collision between the Indian and Eurasian plates during the closure of the Tethys Ocean.Based on a change in plate motion, most researchers suggest collision between India and Eurasia probably began 50-55 Ma (e.g.Molnar & Tapponnier 1975;Patriat & Achache 1984;Copley et al. 2010), although Yin & Harrison (2000) suggest it may have begun earlier and Aitchison et al. (2007) prefer a later collision time of ∼34 Ma.The International Deep Profiling of Tibet and the Himalaya, Phase IV (INDEPTH IV) (Fig. 1) was set up to image the northeastern margin of the Tibetan collision zone and explore the mechanisms for the continued growth of the plateau.In this paper, we describe the results of an arrival time tomography study of the upper mantle beneath northeast Tibet using recordings of P and S waves from earthquakes at teleseismic distances from the INDEPTH IV array.
The evolution of the Tibetan Plateau (Fig. 1) involved numerous collisions, with terranes accreted to Eurasia during the closure of the Tethys Ocean.In eastern Tibet, the Lhasa and Qiangtang terranes are separated by the Bangong-Nujiang Suture (BNS), which began as a rift in the Early Ordovician or Carboniferous, and later became a convergent margin.The BNS closed during the late Jurassic (Dewey et al. 1988).The Qiangtang and the Songpan-Ganzi terranes are separated by the Jinsha Suture (JS), which was formed by the closure of the Songpan-Ganzi Ocean (Yin & Harrison 2000) in the late Triassic or earliest Jurassic (Dewey et al. 1988).The northern edge of the Songpan-Ganzi is marked by the major left-lateral strikeslip Kunlun Fault (KF).The KF has major branches on the eastern side of Tibet: the South Kunlun thrust fault (SKF) and the North Kunlun strike-slip fault (NKF).The crust is approximately 75 km thick beneath the Lhasa Terrane and becomes thinner (63-74 km) north of the BNS (Kind et al. 2002;Yue et al. 2012).The thickness of the lithosphere is less well established.One way to estimate the thickness is to look for the lithosphere-asthenosphere boundary (LAB) using negative converters in receiver function studies.Mechie & Kind (2013) note that three studies show the LAB to be at 150-200 km depth underneath northern Lhasa and the Qiangtang (Kumar et al. 2006;Zhao et al. 2011;Yue et al. 2012).
The Kunlun Shan lies to the north of the KF, where it is marked by a steep drop in topography, and bounds the Qaidam Basin in a series of buried thrusts called the North Kunlun Thrust System (NKTS).The 3 km high Qaidam Basin is a desert, with 8 km of sediment, and a total crustal thickness of 50 km (Karplus et al. 2011;Yue et al. 2012).The Qaidam is bound to the west by the Qimen Tagh Thrust Belt, to the northwest by the left-lateral strike-slip Altyn Tagh Fault and to the northeast by the mainly north-dipping high-angle North Qaidam Thrust System (NQTS).The Qilian Shan chain lies to the north of the basin, flanked by the south-dipping Qilian Shan Thrust Belt (Yin et al. 2008).
The role of the Indian lithosphere in the evolution of Tibet remains a subject of debate, but evidence from surface wave tomography (Shapiro & Ritzwoller 2002;Priestley et al. 2006;Lebedev & van der Hilst 2008), body wave tomography (Li et al. 2008;Liang et al. 2012;Zhang et al. 2012), joint body wave and surface wave tomography (Obrebski et al. 2012), receiver functions (Kind et al. 2002;Nábělek et al. 2009;Zhao et al. 2010;Yue et al. 2012), reflection seismology (Hauck et al. 1998) and gravity modelling (Jiménez-Munt et al. 2008) has led to the general acceptance that the Indian lithosphere has subducted underneath at least the southern part of the plateau.Considerable controversy remains concerning the northern edge of the plateau.In particular, Tapponnier et al. (2001) have posited the existence of southward-directed subduction of Asian lithosphere below northern Tibet, and Kind et al. (2002) point to a prominent south-dipping interface in P-receiver functions to corroborate this hypothesis.Subduction of Asian lithosphere is also supported by the identification of two interfaces that generate negative polarity S-wave mode conversions at ∼100 and ∼170 km depth below northern Tibet that were interpreted by Zhao et al. (2011) as the lower boundaries of stacked mantle lithospheres of Tibetan and Asian provenance.These interfaces appear connected to a shallow LAB at 120 km below the Qaidam Basin.At the same time, the receiver function study of Yue et al. (2012) found no candidate interfaces associated with subduction of Asian lithosphere, and the Rayleigh wave tomography study of Ceylan et al. (2012) found no evidence of continuous higher velocity zones that would be consistent with subducting lithosphere from the north.
High P n (i.e.uppermost mantle) velocities are found beneath the southern plateau and major basins, including the Tarim and Qaidam (Hearn et al. 2004;Liang & Song 2006).A belt of low P n velocities stretches across the Songpan-Ganzi and Qiangtang terranes (McNamara et al. 1997), suggesting that the uppermost mantle beneath the Songpan-Ganzi and Qiangtang terranes is warm and may be weak while that beneath the basins and southern Tibet is strong and cold.
In this study, we use P and S body wave teleseismic data, along with a starting model based on surface wave analysis.Teleseismic regional tomography is sensitive to lateral velocity contrasts, and provides resolution to depths comparable to the aperture of the array.However, because it uses relative arrival times, it is insensitive to absolute wave velocities.It also is prone to smearing images in areas lacking crossing rays, particularly in shallow regions such as the crust (e.g.Aki et al. 1977;Priestley & Tilmann 2009).In the absence of any prior information about the distribution of wave velocities in the crust, investigators generally use either station corrections derived from average residuals (e.g.Frederiksen et al. 1998;Graeber et al. 2002) or force as much of the misfit as possible into the crust (e.g.Abers & Roecker 1991) to isolate signal due exclusively to structure in the mantle.However, this approach can underestimate variations of mantle velocities (Waldhauser et al. 2002;Priestley & Tilmann 2009).Results from previous investigations in northeastern Tibet (e.g.Yue et al. 2012) suggest large variations in Moho depth, as well as considerable lateral heterogeneity in crustal wave velocities.
In contrast to body wave tomography, surface wave tomography generally exhibits good depth resolution in the upper 250-300 km depth, but tends to smear short-wavelength lateral anomalies.We use an approach previously used by Rawlinson & Fishwick (2012), Downloaded from https://academic.oup.com/gji/article-abstract/196/3/1724/583135 by guest on 01 May 2019 which allows the surface waves to account for the shallow heterogeneity and to provide absolute velocities and the long-wavelength background structure (down to ∼400 km, albeit at lower resolution below ∼200 km depth).
We use a surface wave model (based on Acton et al. 2010;Priestley et al. 2006) to construct a starting model for our body wave inversion.This model incorporates the major differences in crustal wave velocities in this region, and represents the long-wavelength background structure in the shallow mantle.We compare our results to those starting with a simple 1-D model, and one where we explicitly account for expected differences in Moho depth.Our approach allows us to make a correction for shallow heterogeneity, but is also able to sample lateral heterogeneity within the mantle.
Relative arrival times for both P and S waves were determined using the multichannel cross-correlation (MCCC) method developed by VanDecar & Crosson (1990).P arrivals were correlated on the vertical component after bandpass filtering between 0.4 and 1.8 Hz; S arrivals were correlated on the transverse component after bandpass filtering between 0.03 and 0.2 Hz.The analysis windows for P-wave correlation were assigned manually for each event, with a mean length of 2.8 s.S-wave correlation was less straightforward as cycle-skipping was common.The analysis window for S waves had a mean length of 14.9 s.Each S-wave correlation was checked manually for cycle-skipping or for later phases influencing the correlation, and seismograms exhibiting poor correlation were removed before recorrelating.At least three picks were made for each of the included events.We include example correlations with the Supporting Information.
To provide an estimate of measurement error, for the purpose of weighting the inversion, uncertainties (σ d ) in the relative P-wave arrival times were estimated by first assigning events into one of three quality classes.We then searched for pairs of observations where events within 20 km of each other were observed at the same station.There were over 100 event-pairs with a minimum of 11 stations and a maximum of 65 stations common to both events (resulting in an average of 989 paired observations per event-pair).Using the statistics calculated from the pairs, we assigned an uncertainty to each of the quality classes.2107 picks were assigned an uncertainty of 0.20 s, 11 126 picks assigned 0.28 s and 590 picks assigned 0.47 s.This compares with initial traveltime residuals of 0.51-0.53s for our P models.
There were insufficient data to allow a similarly quantitative estimate of the S-wave arrival time uncertainties.To at least provide a working estimate for purposes of weighting, we calculated the uncertainty for six event-pairs within 20 km of each other.For all events, we then estimated the uncertainty manually by considering how well the waveforms correlated across the network for each event.Finally, the manually estimated uncertainties were scaled using a scaling factor of 0.7, based on a comparison between the statistics for the six event-pairs and the manual estimates for the same events.The mean assigned uncertainty for the 3957 S arrivals used in this study is 0.7 s, the lowest uncertainty 0.3 s and the maximum uncertainty 1.7 s.This compares with an initial traveltime residual rms of 1.06-1.13s for our S models.While these estimates do not appear unreasonable, we emphasize that they are themselves quite uncertain, and, as discussed below, are likely to be too conservative.
We examined 937 P-wave teleseismic events with an epicentral distance ranging from 30 • to 90 • and moment magnitudes (M W ) between 5.0 and 6.9.498 events were excluded due to low waveform quality, and a further 152 to avoid events from the east and southeast dominating the signal.Our assembled data set contains 13 823 P-wave arrivals recorded from 287 events (Fig. 2, red circles).S-wave arrival times were derived from recordings of 78 events with epicentral distances ranging from 25 • to 80 • and moment magnitudes (M W ) between 5.2-6.9 (Fig. 2, blue diamonds).To balance the azimuthal distribution of data for the S tomography, we restricted selection of events from the east and southeast to those with M W > 5.9.After reviewing waveform quality, we retained 3957 S arrivals from an original set of 4846 (Fig. 2).The ray coverage for both P and S is shown in Fig. 3.

M E T H O D
We analyse the P and S arrival times using a method described in detail in Roecker et al. (2004Roecker et al. ( , 2006) ) and Li et al. (2009).Traveltimes are calculated by adding the time from the base of the regional model to a given station to that from the hypocentre to the base of the model.Times within the regional model are computed using the finite difference eikonal equation solver of Hole & Zelt (1995), adapted for teleseisms and spherical coordinates by Li et al. (2009).Times from the hypocentre to the base of the model are taken from traveltime tables for a standard earth model (IASP91: Kennett & Engdahl 1991).By Fermat's principle, the true traveltime is the minimum of these sums.Ray paths are then determined by following the steepest traveltime gradient, and are updated with each iteration.
Two model grids are defined: a fine spaced grid for the finite difference computation of the traveltimes and ray paths, and a coarser grid for the inversion.We choose a fine grid spacing that is 5 km in depth and 0.0625 • (∼7 km) laterally, and a coarse grid which is 25 km in depth and 0.25 • (∼27 km) laterally.The depth range is from 10 km above sea level to 600 km depth.
In order to mitigate the effects of variation in structure outside the model, we remove the mean and work with only relative traveltimes.Because P-and S-wave traveltimes can be influenced in different ways by structure outside the model, we solve for variations in P and S wave velocities independently.The inverse problems are linearized and the partial derivatives for each observation are calculated along the ray paths.The resulting series of linear equations are solved iteratively using the Least-Squares (LSQR) algorithm of Paige & Saunders (1982).Each row of the system of linear equa-tions is weighted to reflect the uncertainties in the observations.Perturbations are regulated in two ways: first, the least-squares solution is damped to prevent overstepping at any single iteration, and, second, short-wavelength perturbations are smoothed using a moving average window prior to modifying the existing model (e.g.Li et al. 2009).The values one chooses for these regulators (size of damper and length of moving window) are somewhat arbitrary and we experimented with a range of values until we found a combination that led to stable convergence while fitting the observations reasonably well without requiring a too complicated model.Based on our trade-off curves (Fig. 4) we chose a damper of 200 s 2 km −2 and smooth over five nodes in each cardinal direction (the target node and two nodes either side), corresponding to 1.25 • laterally and 125 km in depth.

Starting models
The results presented in this paper use recent surface wave models as a starting model for our teleseismic body wave inversion.This is to mitigate the effects of a highly heterogeneous crust and uppermost mantle.In this section, we describe the preparation of this model (SW), and compare it with two other models: a simple 1-D model (AK) and a model where we explicitly account for Moho depth variations (MO).

Preparation of different starting models
The starting models are compared in Figs 5 and 6.Starting model AK is a simple 1-D model that uses global reference model ak135 (Kennett et al. 1995).Starting model MO was generated using the Moho depth model of Yue et al. (2012), which is based on stacking direct and multiple P s converted phases using the slant-stack method of Zhu & Kanamori (2000).Below the Moho, we used the average at the corresponding depth from model SW (see below).Above the Moho, we assumed a P velocity of 6.4 km s −1 and an S velocity of 3.6 km s −1 .We adjust the nodes either side of the Moho, such that the integrated traveltime approximates the traveltime across the Moho discontinuity.
The S-wave version of starting model SW (Fig. 6a) is based on shear wave velocities obtained from prior surface wave analysis.At depths less than 100 km we use the results of Acton et al. (2010), based on the inversion of Rayleigh group velocities at periods of 10-70 s, for the crust and uppermost mantle.For depths larger than 150 km up to the base of our model grid at 600 km depth we use an updated model based on a fundamental and higher mode Rayleigh wave inversion between 50 and 160 s described in Priestley et al. (2006), but which utilizes 11 times as many waveforms.Below 325 km, this model shows only small (<±2 per cent) variation from PREM (Dziewonski & Anderson 1981).Between 100 and 150 km depth, we grade the Acton and Priestley models into each other.The expected depth resolution is ∼20 km above 100 km depth and ∼30-50 km at greater depths.The expected lateral resolution is ∼100-200 km above 100 km depth and ∼200-400 km at greater depths.
The average P-wave velocity for each depth within starting model SW (Fig. 5a) was estimated from the surface wave models by assuming the same V p /V s ratio at a particular depth as ak135.Variations in mantle velocities normally reflect changes in temperature rather than composition; we expect that compositional changes will result in velocity anomalies <1 per cent.In addition, temperature has a greater effect on S velocities: a 100  lowers V s by 0.7-4.5 per cent and V p by 0.5-2 per cent (Goes et al. 2000).Hence, the variation in mantle P velocities should be less than those in S, and we set the P-velocity anomalies to be only half the magnitude of those in the S-wave model.

Comparing results using different starting models
The final results from the three different models show broadly the same features (Fig. 5a  discussion of the models, we add labels to the larger features for ease of reference.The body waves cannot retrieve the shallowest structure, because there is little crossing ray coverage immediately below the stations, and so the largest differences between the models are within the shallow structure.In model AK, anomaly A is strongly smeared vertically (see Fig. 5a for labels).By comparison, in model SW, the depth range is more constrained, and consequently the amplitude of the anomaly is higher.Below 300 km, the models are very similar because of improved crossing ray coverage at these depths, and also because there is less structure in starting model SW.
The map views of the results (Fig. 5b) also show broadly similar features, especially C, D and E. However, the boundary of A (shown by white dashed lines) is 50-100 km further north in model SW.Recovery is affected by the station distribution: there are holes where the coverage is sparser.
Features A-G and I are also recovered in the S-wave models (Fig. 6a).Neither Model AK nor MO retrieve H, although it remains in model SW with reduced amplitude in comparison to starting model SW.Feature A does not reach as far north in model AK and MO as in SW (Figs 6a and b), probably because it is obscured by vertical smearing of feature B.

Reduction of the residuals
All three sets of models show a large reduction in the residuals (Fig. 7).For P, there are small differences for initial fit, with model AK the best, and model SW the worst (Fig. 7a).For S, there are slightly larger differences in the initial fit, with model SW the best, and model MO the worst (Fig. 7b).Surprisingly, in both cases, the models eventually reach nearly the same fit to the residuals.
The different starting models can generate differences in the relative arrival times of up to 1.4 s for S, though typical differences are on the order of one or two-tenths of a second, substantially less than the estimated picking errors.Changing the starting models (which affect only the shallow structure in this case) results in almost no change to the fit of the residuals.What differences there are in calculated traveltimes are absorbed in relatively minor changes in model perturbations at all depths (Figs 5 and 6, third row).We find that our images are robust, with only minor differences between the models below 300 km, despite major differences in the shallow structure that reflect the differences in starting models.
We note that for P, the residuals converge to a value close to our estimate of the picking error (Fig. 7).As discussed earlier (Section 2), we were not able to determine a reliable quantitative estimate of S uncertainty from the correlated waveforms, and the final residuals suggest that our original estimate of 0.7 s was too conservative.The final rms of the fit is 0.52 s, which is still nearly twice that for P, and so would seem to be a reasonable estimate of actual S uncertainty.
The residuals do not indicate that one model is statistically better than the others.However, the main differences between models AK, MO and SW are above 200 km, in the region where surface waves have good resolution.As an additional test on our confidence in model SW, we calculated surface wave group traveltimes through starting model SW, and compared them with the traveltimes for all three final models.We found that including the body waves into model SW introduced only small variations to the surface wave traveltimes, and that SW provided a significantly better fit to the surface waves than the other two models (Supporting Information Fig. S9).Therefore, we prefer to use model SW, since it is compatible with two independent data sets (the body and surface waves).We use it to provide constraints for the shallow structure and the longwavelength deeper structure.

Synthetic tests
We conducted a series of tests to evaluate the resolution and robustness of our results.Using the same method as our forward calculation (Section 3), we calculated traveltimes using an eikonal equation solver within the test model and traveltimes from a standard earth model outside.The data sets we use for these synthetic tests mimic the real data as closely as possible.We use the same source-receiver pairs as in the real distribution, and add random Downloaded from https://academic.oup.com/gji/article-abstract/196/3/1724/583135 by guest on 01 May 2019 We show the starting model (top and second rows), the difference between starting and final model (third row) and the final model (fourth row, with major features labelled).The starting models in the top row are plotted as absolute seismic velocities (in km s −1 ).As the body waves have no sensitivity for layer averages, the remainder of the profiles are plotted as percentage anomalies relative to the average at any particular depth for that particular model.This has the effect of balancing positive and negative anomalies at each depth, and we plot all our results in this way.Topography and stations are plotted as in Fig.  Gaussian errors with the same standard deviation as the picking error (0.26 s for P and 0.57 s for S).

Synthetic test: 'complex' model
We were interested in how well a complex (but known) model can be retrieved.Model SW (after 30 iterations) was used as the input (Fig. 8a).We compare the results when the upper part of the model is provided as part of the starting model versus a much simpler starting model.For one starting model, we used Model AK (Sy-AK, Fig. 8b).For the second starting model we used a combined model which replicated the upper 125 km exactly and used ak135 below that depth (Sy-CO, Fig. 8c).
Below 125 km the features are generally well recovered, and the output models are quite similar.As expected, Sy-AK recovers only   the broadest features above 125 km.Note that A, which is nearly horizontal in the input model, becomes close to vertical, especially on Sy-AK (Fig. 8b), and in both cases the strong fast wave velocities at A and G in the input model smear into each other.

Checkerboard tests
In order to understand the resolution we can achieve, as well as exploring how well models can retrieve an image below a complicated crust and uppermost mantle, we carried out a variation on the standard 'checkerboard' test (Fig. 9).Below 125 km we created a model consisting of alternating positive and negative velocity anomalies of 5 per cent on a background of ak135.The crust and uppermost mantle structure of model SW was used for the upper part of the input model.The checker dimensions were 1.25 • width in latitude and longitude, and 75 or 100 km depth, separated by neutral regions of 0.75 • in latitude and longitude and 50 km in depth in order to be able to better judge smearing.
The P checker pattern (including the crust and upper mantle structure) is alternately inverted using an initial model which replicates the upper 125 km, with ak135 below (PCh-CO, Fig. 9a), or a model which only has the simple ak135 structure (PCh-AK, Fig. 9b).
The checker patterns are well retrieved in both cases, although smearing occurs along the ray paths at the edge of the model, where there is no crossing of the rays.The checkers have marginally improved amplitudes for PCh-CO, and marginally less smearing into the crust (Fig. 9a).PCh-AK reflects the major scale shallow structure, but misses the detail.The S checker pattern is a little more smeared towards the edges, and amplitude recovery a little lower, but the overall pattern is similar (SCh-CO, Fig. 9c).

Synthetic test for probable structures
In order to evaluate how different types of features may be recovered or generate artefacts, we created a model based on some simple shapes designed to represent structures that we may be able to image.We are interested in determining how various features might smear, the effect of uneven coverage, and how dealing with relative anomalies might affect the final picture.
The input model contains a number of rectangular regions, labelled A, B, C, G, H and L (Fig. 10, left-hand column), which represent structures that we find in our models (Figs 11 and 12).Note that the model is plotted as relative percentage anomalies, and that only the labelled structures are included in the input models.
The P results (Fig. 10, centre column) show that features A and C are well retrieved at 175 km, as are G and L at 500 km.However, E and K, and to a lesser extent D, F and I, all appear despite being Downloaded from https://academic.oup.com/gji/article-abstract/196/3/1724/583135 by guest on 01 May 2019 absent from the input model.This is due to a combination of smearing and our imaging of relative anomalies across boundaries.Feature H is not particularly well retrieved, and G becomes connected to A, despite being separated by a 200 km gap in the input model.At 100 km, A has smeared upwards.The southwest of the model is slower than the southeast, consistent with the location of B, but the slow velocities of B are only observed north of the area of A.
The recovery using S waves is similar (Figs 10 m-r), with slightly more vertical smearing (e.g.C in Fig. 10o), and poorer recovery for H (Fig. 10o).Note that we are adding noise based on measurement errors which may be too conservative for the S waves (Sections 2 and 3.1.2),and thus the results may also be too conservative.

R E S U LT S
As discussed earlier, we prefer model SW and base our discussion on it.The images of V p and V s we obtain from modelling P and S-wave arrival times (Figs 11 and 12, equivalent figures for MO and AK are shown in Figs S5-S8) reveal considerable heterogeneity beneath northeastern Tibet both laterally and as a function of depth.We find that there is considerable similarity between the main features for P and S, except that the S model shows less short wavelength structure.Also, we find that our P-wave results move further from the starting point.This behaviour may be due to the greater number of P-wave observations, but could also be a result of the P-wave starting model being less accurate than the S-wave version, because we have to make assumptions on V p /V s when translating the velocities from S to P. Note that the body waves do not sample the regions outside the area bounded by the black lines in these figures, and so the anomalies outside these regions come exclusively from the surface wave starting model.As discussed in Section 3.1.2,the shallow structure is controlled mostly by starting model SW, and the structure below ∼300 km is controlled by the body waves.Between ∼100 and 300 km, both are exerting an influence.
We find that the range of percentage perturbations imposed by the body waves is higher in P than in S. For example, at 175 km, the total percentage range is 9.3 and 6.6 per cent for P and S, respectively.This is likely to be due to the larger P data set, and possibly also due to different choices of damping parameters.However, the combined influence of the body waves and the starting models results in final models which have a greater range in S than in P (since we set the P percentage anomalies to be only half the percentage anomalies of S in the starting models).

P-wave results
Feature A is a shallow dipping, fast anomaly that extends from ∼140-240 km depth (Fig. 11a), and across the southern edge of the model at 175 km (Fig. 11b).The amplitude reduces gradually between the two dashed white lines.B is a slow wave velocity region on the western side of the region, extending from the Kunlun Shan, through the Songpan-Ganzi into the Qiangtang Terrane, up to the dashed white line (Fig. 11g).B extends from mid-crustal depths to ∼120 km on sections L1 and L2 (Figs 11c and a) At 175 km, the eastern side of the Qaidam Basin is faster than the west (C and D in Fig. 11b).E is a slow anomaly extending across the Songpan-Ganzi Terrane and the Kunlun Shan, terminating sharply near the NKTS (Figs 11a, b and d).Its southern boundary lies close to the JS, up to the boundary with A. E is between the fast regions A and C, and so potentially may only be slow in comparison to these areas, rather than exhibiting particularly slow wave velocities.In L2 (Fig. 11a), E blends in the west into the deeper slow anomaly F, which lies beneath the Songpan-Ganzi, Kunlun Shan and Qaidam Basin at depths of ∼250-475 km.
G is a deep fast anomaly beneath the Qiangtang and Lhasa terranes west of 92 • E (Fig. 11i).On L2 (Fig. 11a), it is separated from A by the slow anomaly H.In L1, at the very edge of our coverage (Fig. 11c), there is a connection through a fast structure J, although this is probably due to smearing (Section 3.3).I is another deep fast anomaly that lies beneath the Kunlun Mountains and reaches to the northern edge of the model (L2, Fig. 11a).There are also fast structures to the east and west (L1 and L3), although the shape is less well defined where the coverage is poorer.Along L3, K is a slow anomaly from the IYS to the JS.In the south, it lies close to the base of the model, but is ∼100 km shallower further north (Figs 11d and i).L is a deep linear feature beneath L5 concentrated below 400 km depth (Figs 11e, f and i).

S-wave results
The S-wave results (Fig. 12) are generally similar to the P-wave results (Fig. 11).All of the features described in the P-wave model appear in the S-wave model as well, but with some variation in shape and amplitude.We focus on the important differences here.
Feature A extends 50-100 km further north (Fig. 12b).There is an area on the eastern side which extends ∼200 km north of the NKF, although it is arguable whether or not this forms a continuous structure (it is interrupted by slower velocities on L4, Fig. 12e).C is shallower, and the boundaries are less steep (Fig. 12a).In L2 (Fig. 12a), there is some connection between A and G, but it is low Downloaded from https://academic.oup.com/gji/article-abstract/196/3/1724/583135 by guest on 01 May 2019   , L4 and L5 (c, a, d, e, f) and at 125, 175, 300 and 550 km depth (g, b, h, i).The same features are labelled as in the P-wave models (Fig. 11), although the location may be shifted.The amplitude of A in the S-wave model reduces gradually between the two dashed blue lines on the map view at 175 km, and the white lines show the equivalent for the P-wave results (Fig. 11).White lines indicate the extent of B at 125 km and G and L at 550 km.
amplitude, and they continue to appear as separate structures.On L1 and L3 (Figs 12c and d) A and G are disconnected.There is a slow anomaly along the bottom of L2 (Fig. 12a), which is fast in the P waves (although also slower than G and I).Finally, the strongly defined linear shape of L in the P waves is more diffuse in the S waves (Figs 11i and 12i).
A and H are stronger features.It is likely that H remains from the starting model, rather than being clearly retrieved by the body Downloaded from https://academic.oup.com/gji/article-abstract/196/3/1724/583135 by guest on 01 May 2019 waves (Figs 12a, c, d and h).However, A is also well retrieved by all the S-wave models (Fig 6).

The lithosphere underneath Tibet
Beneath the Lhasa and Qiangtang terranes, our images show a shallow dipping, fast feature at depths of ∼120-220 km, below a high-amplitude slow anomaly from ∼30-110 km (A and B , respectively, in Figs 11 and 12; interpretative cartoon Fig. 13).The overall shape of A comes from starting model SW, but the body waves have sharpened the change from fast to slow, and moved the boundary further south.Although P-and S-wave velocities behave differently in cases of temperature and compositional variations, we find that the shape of the structure is similar in our P and S models.However, as upper-mantle seismic velocities are much more sensitive to temperature than to compositional variations (e.g.Goes et al. 2000;Schutt & Lesher 2006), the fast velocities are likely to be due to colder than average temperatures irrespective of any compositional anomaly which might additionally exist.
Several other studies have imaged a continuous high velocity body stretching from the Himalaya northwards (e.g.Shapiro & Ritzwoller 2002;Lebedev & van der Hilst 2008;Li et al. 2008;Priestley et al. 2008;Hung et al. 2010;Zhang et al. 2012) although the apparent northward extent of this anomaly differs by several 100 km between different methodologies, and also along the strike of the Himalayan Arc.Hence, feature A is commonly interpreted as underthrusting Indian lithosphere (e.g.Kind et al. 2002;Li et al. 2008;Chen et al. 2010).
Alternatively, A could represent lithospheres of two different origins.The southern part could originally come from India, and the northern part could represent an underthrust Lhasa Terrane lithospheric slab.This has been proposed by Yue et al. (2012), based on a shallow dipping (6 • -10 • ) interface that generates mode conversions in receiver functions and dips from ∼80 km below the BNS to ∼120 km below the Kunlun Mountains.A belt of Eocene-Oligocene potassic volcanism beneath the Qiangtang suggests that there may have been earlier continental subduction (Roger et al. 2000;Ding et al. 2003), although it is unclear whether it would have been from the north or south.Our results suggest that there is no present subduction, because there is no connection between fast structures at 140 km and the surface.
A third explanation for A is that pure shear thickening of the entire collision zone has forced some lithospheric material deeper, which will make it cooler than normal material at the same depth, until it can return to thermal equilibrium (McKenzie & Priestley 2008).
A popular model for plateau uplift involves delamination or removal of the bottom part of the lithosphere due to Rayleigh-Taylor  (Craig et al. 2012).Qaidam lithosphere is faster than normal lithospheric mantle.Between the Qaidam and the underthrusting Indian Plate is an area of relatively warm material, but it is not hot enough to cause melt.Below 350 km, we see an area of slower mantle, which coincides with a deepening of the 660 discontinuity (Yue et al. 2012), which could potentially be the remains of a disconnected oceanic slab which has passed through the 410.A depression in the 410 km discontinuity coincides with a slower area of mantle in our results, and may indicate small-scale upwelling.
instabilities (e.g.Molnar 1988;Molnar et al. 1993).We image fast velocities underneath the southern plateau, suggesting that the lithosphere is thickened rather than delaminated.
Assuming a collision age of ∼50 Ma, plate reconstructions from marine magnetic anomalies and fracture zone reconstructions show that there may be up to ∼3200 km of India-Asian convergence accommodated on the eastern side of Tibet, much more than the 450-900 km of documented Himalayan shortening (van Hinsbergen et al. 2011).There is thus sufficient Indian mantle lithospheric material available to account for anomaly A, but seismic methods cannot uniquely identify the provenance of the material in the collision zone.
An explanation for the shallow low velocities (B) was put forward by McKenzie & Priestley (2008) who argue that radiogenic heat production in the thickened crust of Tibet is sufficient to heat the mid-crust, which eventually heats the lower crust and the upper mantle.This causes a temperature inversion, with temperatures highest in the mid-crust, and a negative temperature gradient beneath this region.This heating could result in an appreciable reduction in the shear wave velocities in the thickened mid-crust after only 10 Ma (McKenzie & Priestley 2008).It also could heat the lithosphere underneath (either the underthrusting Indian lithosphere or remnant Qiangtang lithosphere) to ∼120 km depth (Craig et al. 2012).
A number of studies have seen large differences in P n and shear wave splitting between northern and southern Tibet.High P n (i.e.uppermost mantle) velocities are found beneath most of the plateau south of the BNS, contrasting with a belt of low P n velocities that stretches across the Songpan-Ganzi and Qiangtang terranes (McNamara et al. 1997;Hearn et al. 2004;Liang & Song 2006).The significant shear wave splitting observed in northern Tibet contrasts with the isotropic or weakly anisotropic fabric of southern Tibet and the Indian Shield (Chen & Özalaybey 1998;Huang et al. 2000;Chen et al. 2010).These changes to the shear wave splitting and P n have been interpreted as marking the edge of the advancing Indian lithosphere.Instead, we suggest that these changes mark the edge of the region where significant heating occurs (the onset of B).Heating of the uppermost mantle would lower P n and S n .Similarly, Craig et al. (2012) suggest that the transition from weak anisotropy or isotropy to significant anisotropy may reflect where the material becomes hot and weak enough to deform significantly and develop an anisotropic fabric.However, it is also possible that B represents partially asthenospheric material or remnant weak Tibetan lithosphere that has been underthrust by Indian mantle lithosphere.
An E-W contrast, similar to our shallow feature B, also appears in the S-wave tomographic model of Liang et al. (2012) and is also consistent with the P n study of Liang & Song (2006), which found low uppermost mantle velocities in the Songpan-Ganzi Terrane and the northern and western Qiangtang terrane, with particularly low wave velocities corresponding to B. The E-W contrast is also consistent with an area of S n blockage seen by Barron & Priestley (2009).The propagation of the S n phase is reduced in areas where there is a negative shear wave velocity gradient.At high frequencies, S n is blocked across the entire plateau.However, in the mid-frequencies, S n propagates with reduced efficiency in areas west of 93 • E in our area (corresponding to B).As described earlier, a negative shear wave velocity gradient can occur when the mid-crust is heated.At low frequency, S n waves are able to propagate efficiently across the whole of Tibet, suggesting that the low velocities below the Qiangtang are limited to the uppermost few tens of kilometres of the mantle.
Further west (at ∼90 • E), between the IYS and BNS, Tilmann et al. (2003) imaged a fast, steeply dipping body that extends from at least 100-350 km depth.We see a fast anomaly below 400 km, on the southwestern side of our region (G in Figs 11a and i), but it is separated from A by a region of slower velocities.In our area, we observe no well-resolved steeply dipping features, although in areas of poorer coverage, we do see a connection between A and G.However, this is likely to be due to smearing.Our synthetic test (Fig. 10) shows that two anomalies (one shallow and one deep) can easily become connected.

Slow velocities beneath the Songpan-Ganzi and Kunlun Shan
The relatively narrow slow wave velocity region E beneath the Songpan-Ganzi Terrane and the Kunlun Shan at 175 km (Fig. 11b) is also seen in other tomographic studies (e.g.Li et al. 2008;Ceylan et al. 2012;Liang et al. 2012;Obrebski et al. 2012;Zhang et al. 2012).Zhang et al. (2012) suggested this area may be a mantle diapir upwelling between Indian and Asian lithospheres, and Ceylan et al. (2012) suggested that it may be shear heating along the KF.However, there is no recent volcanism in the Songpan-Ganzi, east of ∼92 • E (Searle et al. 2011), so this area is presumably not hot enough to produce melt.It is also likely that the lithosphere is 150-200 km thick beneath the Songpan-Ganzi and Kunlun Shan (Kumar et al. 2006;Zhao et al. 2011;Yue et al. 2012).Therefore, we think that E is slow only in a relative sense.It is sandwiched between two particularly fast areas A and C at the same depth.

No evidence of southward subduction of Asian lithosphere
As discussed in Section 1, a number of investigators have proposed the southward subduction of Asian lithosphere beneath northern Tibet (Tapponnier et al. 2001;Kind et al. 2002;Zhao et al. 2011).If this were the case, we would expect to see a connection between the Qaidam lithosphere and structures in the plateau, and potentially see fast structures underthrusting or dipping into the mantle.Like Liang et al. (2012), Obrebski et al. (2012) and Li et al. (2008), we image relatively slower material here, which is difficult to reconcile with continuous southward subduction.Similarly, the Rayleigh wave tomography of Ceylan et al. (2012) and Lebedev & van der Hilst (2008) and the receiver function images of Yue et al. (2012) and Karplus et al. (in preparation) found no evidence of subduction in this region.

Qaidam Basin
At depths of ∼125 km, the Qaidam Basin is underlain by a fast region, particularly on the eastern side.Similarly, Hearn et al. (2004) and Liang & Song (2006) find high P n (i.e.uppermost mantle) velocities beneath the Qaidam.The LAB is observed at depths of 100-150 km with S-receiver functions (Zhao et al. 2011) and Yue et al. (2012) also infer a shallower Qaidam LAB on their western profile than on their eastern profile.
The Qaidam has been suggested by many workers to act in a similar way to the larger Tarim Basin (e.g.Braitenberg et al. 2000).Both basins are thought to be cratonic in origin, and overlain by deep sediments.Interpretation of the seismicity and GPS observations has shown that the Tarim rotates like a rigid block (Braitenberg et al. 2000;Reigber et al. 2001).In the Qaidam, GPS observations show that it is currently taking up more shortening than the Downloaded from https://academic.oup.com/gji/article-abstract/196/3/1724/583135 by guest on 01 May 2019 Songpan-Ganzi or Qiangtang terranes, albeit considerably less than the mountains of the Qilian Shan (Gan et al. 2007).Yin et al. (2008) show that Cenozoic upper crustal shortening decreases eastwards across the basin from >48 per cent in the west to <1 per cent in the east.Our observation that the east is faster (and presumably colder and more rigid) than the west may explain why there has been less shortening in the east.

Mantle structure beneath 300 km
We see a deep, fast anomaly in the southwest (G in Figs 11i, a and  c) and another in the northwest (I in Figs 11i and a), separated by a region which is slower (in a relative sense) at the same depths.G corresponds to a depression in the 660 km discontinuity (Yue et al. 2012).A reduction of the velocity contrast at the 660 had also been inferred from modelling the triplication of seismic waveforms with bounce points further west, in central Tibet (Chen & Tseng 2007).Following earlier suggestions by the works just cited and others, we therefore interpret D as a fragment of lithosphere that has detached and sunk into the mantle.This lithospheric fragment could either be a part of the detached Tethyan oceanic slab, or could be related to subduction or delamination of lithosphere after the onset of the continent-continent collision.We also see deep slow structures K and L. Both may represent small upwelling structures (with K possibly related to the emplacement of G).Alternatively, our synthetic tests suggest K may be an artefact of G since a low amplitude slow structure is recovered next to G (Figs 10h and n).Deep slow wave velocities have been previously observed in southeastern Tibet, where we see K and L (Koulakov 2011).
Finally, E, F and H (Fig. 11a) combined correspond to slow wave velocities imaged by Wittlinger et al. (1996).Also, Yue et al. (2012) see a thinning of the transition zone (probably due to deepening at the 410 discontinuity) in a location approximately corresponding to F. However, in our synthetic test (Fig. 10), the contrast between A and the rest of the region, as well as smearing from shallow anomaly B, was sufficient to retrieve at least some of the amplitude of slow anomaly at F. Wittlinger et al. (1996) suggested that E, F and H may be a mantle plume of hot upwelling material.However, they note no melt from the plume can have reached the surface, since the Neogene volcanism in the Qiantang is further west than their profile (Yin & Harrison 2000).We think it is unlikely that there is a major plume in the area.

C O N C L U S I O N S
We have constructed images of the P-and S-wave velocity structure beneath northeastern Tibet by analysing the arrival times of teleseismic body waves.We used recent surface wave models as a starting point for our inversions.We find similar features in the S-and P-wave results.
There is a shallow dipping structure which we interpret as the Indian Plate, which is at least 100 km thick and underthrusts below depths of ∼140 km, and dies out gradually between 33 • N and 35 • N.Alternatively, it could represent a region which has been thickened by pure shear or the northern part could represent an underthrust Lhasa Terrane lithospheric slab.The deep structure below 300 km shows colder areas of mantle, consistent with receiver function observations of a thickened transition zone, which could potentially be a fragment of oceanic lithosphere.In this area, we find that this fragment is not connected to faster structures above.We find slow and presumably warm uppermost mantle, particularly on the western side of our model.This coincides with an area of reduced propagation of the P n phase (Barron & Priestley 2009), and is likely to be due to warming of the uppermost mantle due to radiogenic heating in the thickened crust (McKenzie & Priestley 2008).We find that the Qaidam Basin is faster than normal mantle wave velocities to depths of around ∼180 km, with the east faster than the west.The fast region coincides with an area which has seen only small amounts of shortening over the Cenozoic (Yin et al. 2008).We find no connection between structures beneath the Qaidam and further south in the plateau.From this we conclude that there is no evidence for southward subduction of Eurasian lithosphere.

A C K N O W L E D G E M E N T S
We thank all members of the scientific teams of the  (Wessel & Smith 1998).We thank James Mechie, Dan McKenzie, Timothy Craig, Alex Copley and Tuna Eken for their comments on the draft manuscript.The manuscript has been improved by comments from anonymous reviewers and the editor.

S U P P O RT I N G I N F O R M AT I O N
Additional Supporting Information may be found in the online version of this article:

Figure 1 .
Figure 1.A map of the study region showing terranes and major faults (after Yin et al. 2008): barbed lines are thrusts, continuous lines are strike-slip and dashed lines are sutures.Red triangles are INDEPTH IV stations, and lines L1-L5 are the location of our profiles.The location of the study area is shown on the right.Abbreviations on the map are: MFT, Main Frontal Thrust; LT, Lhasa Terrane; QT, Qiangtang Terrane; SG, Songpan-Ganzi Terrane; KS, Kunlun Shan; QS, Qilian Shan; IYS, Indus-Yalu Suture; BNS, Bangong-Nujiang Suture; JS, Jinsha Suture; JRF, Jinsha River Fault; KF, Kunlun Fault; SKF, South Kunlun Fault; NKF, North Kunlun Fault; NKTS, North Kunlun Thrust System; QTFS, Qimen Tagh Fault System; QSTB, Qilian Shan Thrust Belt; NQTS, North Qaidam Thrust System; AT, Altyn Tagh; AB, Alxa Block.Note that there is some dispute about the nature of the SKF and the vergence of the NKTS.

Figure 3 .
Figure 3. P-wave (left-hand column) and S-wave (right-hand column) ray coverage along profiles L1, L2, L3, L4 and L5.Rays are shown within 60 km of each profile.The topography and major faults and sutures are projected onto the profile, along with stations within 50 and 150 km of the profile (blue and grey triangles, respectively).Profile locations in Fig. 1.

Figure 4 .
Figure 4.The trade-off curves between the roughness of the imposed model perturbations and rms residual.(a) P-wave inversions and (b) S-wave inversions after four inversions for model SW.Dashed lines connect the same damping parameter, and solid shapes show the smoothing parameter.We chose damping of 200 s 2 km −2 and smoothing across five nodes for both P and S (circled in red).The star shows the initial residual (the initial roughness is zero because we consider model perturbations).

Figure 5 .
Figure 5.A comparison of P-wave velocity models resulting from three different starting models.Panels in (a) show profiles along L2 for Model AK, MO and SW.We show the starting model (top and second rows), the difference between starting and final model (third row) and the final model (fourth row, with major features labelled).The starting models in the top row are plotted as absolute seismic velocities (in km s −1 ).As the body waves have no sensitivity for layer averages, the remainder of the profiles are plotted as percentage anomalies relative to the average at any particular depth for that particular model.This has the effect of balancing positive and negative anomalies at each depth, and we plot all our results in this way.Topography and stations are plotted as in Fig. 3. Panels in (b) show map views for the starting model (top row), difference (middle row) and final model (bottom row) at 175 km depth.For comparison, dashed white lines indicate the northern edge of A for each model.The map views are also plotted relative to the layer average.Stations are shown as grey triangles.The extent of the body wave ray coverage is indicated by the black line on both the map views and the profiles.
Figure 5.A comparison of P-wave velocity models resulting from three different starting models.Panels in (a) show profiles along L2 for Model AK, MO and SW.We show the starting model (top and second rows), the difference between starting and final model (third row) and the final model (fourth row, with major features labelled).The starting models in the top row are plotted as absolute seismic velocities (in km s −1 ).As the body waves have no sensitivity for layer averages, the remainder of the profiles are plotted as percentage anomalies relative to the average at any particular depth for that particular model.This has the effect of balancing positive and negative anomalies at each depth, and we plot all our results in this way.Topography and stations are plotted as in Fig. 3. Panels in (b) show map views for the starting model (top row), difference (middle row) and final model (bottom row) at 175 km depth.For comparison, dashed white lines indicate the northern edge of A for each model.The map views are also plotted relative to the layer average.Stations are shown as grey triangles.The extent of the body wave ray coverage is indicated by the black line on both the map views and the profiles.

Figure 6 .
Figure 6.A comparison of S-wave velocity models resulting from three different starting models.Figure format as in Fig. 5.

Figure 7 .
Figure 7.The reduction in the residual for (a) P-wave models and (b) S-wave models.Model AK in green, MO in blue and SW in red.The dashed grey lineshows the residual, we might expect given our a priori estimation of the measurement errors.Note that we believe that our estimate of the measurement errors may be too conservative for S.

Figure 8 .
Figure 8. Synthetic test to understand how well we can retrieve a known, complicated model.(a) The input is Model SW (after 30 iterations).(b) Sy-AK: the inversion is carried out using Model AK as the initial model.(c) Sy-CO: using the same synthetics, but with Model CO as the initial model.This has simple structure below 125 km (marked by the dashed white line) but replicates the complicated structure above.Shown along profile L2 (Fig. 11).

Figure 9 .
Figure 9.We created a checkerboard model consisting of a complicated upper 125 km above a checkerboard pattern imposed on a 1-D model.Coloured rectangles show 5 per cent positive and negative anomalies in the input model; areas in between have 0 per cent anomaly in the input model.(a) PCh-CO is the P-wave checkerboard retrieved using a starting model which replicates the upper 125 km (above dashed white line).(b) PCh-AK is the P-wave checkerboard retrieved using Model AK as the starting model.(c) SCh-CO is the S-wave checkerboard retrieved using a starting model which replicates the upper 125 km (the equivalent of Fig. 9a for S).The profile is along 95.5 • (indicated on map views).The map views are at 225 and 550 km.The initial perturbation was ±5 per cent, and the location of each checker is outlined.

Figure 10 .
Figure 10.Input to a synthetic test for probable structures (left-hand column), P results (cental coulmn) and S results (right-hand column).Results shown along L2 and L4, and at 175, 500, 100 and 300 km depth.Note that perturbations were only applied to the labelled features in the input model (white labels).As before, the model is plotted as relative percentage anomalies, and so the background colours change in some cases, most obviously in the 150-250 km depth range.Black text labels either indicate smearing or features not originally in the input model.

Figure 11 .
Figure 11.Results: P-wave velocity anomalies along L1, L2, L3, L4 and L5 (c, a, d, e, f) and at 125, 175, 300 and 550 km depth (g, b, h, i).The model uses teleseismic body wave data and an initial starting model provided by surface waves.The extent of the body wave ray coverage is indicated by the black line on both the map views and the profiles, and so anomalies outside this region come purely from the starting model.Major features are labelled.The amplitude of A reduces gradually between the two dashed white lines on the map view at 175 km.White lines indicate the extent of B at 125 km and G and L at 550 km.As before, both profiles and map views are plotted as relative percentage anomalies (see Fig. 5 for fuller explanation) and topography and stations are plotted as in Fig. 3.

Figure 12 .
Figure 12.Results: S-wave velocity anomalies along L1, L2, L3, L4 and L5 (c, a, d, e, f) and at 125, 175, 300 and 550 km depth (g, b, h, i).The same features are labelled as in the P-wave models (Fig.11), although the location may be shifted.The amplitude of A in the S-wave model reduces gradually between the two dashed blue lines on the map view at 175 km, and the white lines show the equivalent for the P-wave results (Fig.11).White lines indicate the extent of B at 125 km and G and L at 550 km.

Figure 13 .
Figure 13.Interpretation of the crust and upper mantle architecture of the central Tibetan Plateau (92 • -93 • E), modified from Yue et al. (2012).The central box indicates where our coverage is best.We observe a fast region underneath the plateau (dark green), which we interpret as underthrusting Indian lithosphere.Alternatively, it could represent a region which has been thickened by pure shear or the northern part could represent an underthrust Lhasa Terrane lithospheric slab.Radiogenic heating (McKenzie & Priestley 2008) of the thickened crust has heated the lower crust and upper mantle beneath the Qiangang and Songpan-Ganzi terranes, and the top of the underthrusting lithosphere(Craig et al. 2012).Qaidam lithosphere is faster than normal lithospheric mantle.Between the Qaidam and the underthrusting Indian Plate is an area of relatively warm material, but it is not hot enough to cause melt.Below 350 km, we see an area of slower mantle, which coincides with a deepening of the 660 discontinuity(Yue et al. 2012), which could potentially be the remains of a disconnected oceanic slab which has passed through the 410.A depression in the 410 km discontinuity coincides with a slower area of mantle in our results, and may indicate small-scale upwelling.

Figure S1 .
Figure S1.Histograms showing initial and final residuals for P and S models AK, MO and SW.FigureS2.An additional 5939 P-wave arrivals from 152 events (green squares) were used to test the fit to the various models (Fig.S1).FigureS3.Example of a waveform for a typical P-wave event.Figure S4.Example of a waveform for a typical S-wave event; traces are aligned according to relative delays determined by crosscorrelation.Figure S5.Model AK (P waves).Figure S6.Model MO (P waves).Figure S7.Model AK (S waves).Figure S8.Model MO (S waves).Figure S9.A comparison between predicted group traveltimes from final model AK, MO and SW with starting model SW, for a group of surface wave paths.Figure S10.We generate noise-free synthetic P-wave traveltimes through starting model MO (a) and SW (c), and recover these traveltimes against starting model AK (http://gji.oxfordjournals.org/lookup/suppl/doi:10.1093/gji/ggt476/-/DC1).

Figure S2 .
Figure S1.Histograms showing initial and final residuals for P and S models AK, MO and SW.FigureS2.An additional 5939 P-wave arrivals from 152 events (green squares) were used to test the fit to the various models (Fig.S1).FigureS3.Example of a waveform for a typical P-wave event.Figure S4.Example of a waveform for a typical S-wave event; traces are aligned according to relative delays determined by crosscorrelation.Figure S5.Model AK (P waves).Figure S6.Model MO (P waves).Figure S7.Model AK (S waves).Figure S8.Model MO (S waves).Figure S9.A comparison between predicted group traveltimes from final model AK, MO and SW with starting model SW, for a group of surface wave paths.Figure S10.We generate noise-free synthetic P-wave traveltimes through starting model MO (a) and SW (c), and recover these traveltimes against starting model AK (http://gji.oxfordjournals.org/lookup/suppl/doi:10.1093/gji/ggt476/-/DC1).

Figure S3 .
Figure S1.Histograms showing initial and final residuals for P and S models AK, MO and SW.FigureS2.An additional 5939 P-wave arrivals from 152 events (green squares) were used to test the fit to the various models (Fig.S1).FigureS3.Example of a waveform for a typical P-wave event.Figure S4.Example of a waveform for a typical S-wave event; traces are aligned according to relative delays determined by crosscorrelation.Figure S5.Model AK (P waves).Figure S6.Model MO (P waves).Figure S7.Model AK (S waves).Figure S8.Model MO (S waves).Figure S9.A comparison between predicted group traveltimes from final model AK, MO and SW with starting model SW, for a group of surface wave paths.Figure S10.We generate noise-free synthetic P-wave traveltimes through starting model MO (a) and SW (c), and recover these traveltimes against starting model AK (http://gji.oxfordjournals.org/lookup/suppl/doi:10.1093/gji/ggt476/-/DC1).

Figure S4 .
Figure S1.Histograms showing initial and final residuals for P and S models AK, MO and SW.FigureS2.An additional 5939 P-wave arrivals from 152 events (green squares) were used to test the fit to the various models (Fig.S1).FigureS3.Example of a waveform for a typical P-wave event.Figure S4.Example of a waveform for a typical S-wave event; traces are aligned according to relative delays determined by crosscorrelation.Figure S5.Model AK (P waves).Figure S6.Model MO (P waves).Figure S7.Model AK (S waves).Figure S8.Model MO (S waves).Figure S9.A comparison between predicted group traveltimes from final model AK, MO and SW with starting model SW, for a group of surface wave paths.Figure S10.We generate noise-free synthetic P-wave traveltimes through starting model MO (a) and SW (c), and recover these traveltimes against starting model AK (http://gji.oxfordjournals.org/lookup/suppl/doi:10.1093/gji/ggt476/-/DC1).

Figure S6 .
Figure S1.Histograms showing initial and final residuals for P and S models AK, MO and SW.FigureS2.An additional 5939 P-wave arrivals from 152 events (green squares) were used to test the fit to the various models (Fig.S1).FigureS3.Example of a waveform for a typical P-wave event.Figure S4.Example of a waveform for a typical S-wave event; traces are aligned according to relative delays determined by crosscorrelation.Figure S5.Model AK (P waves).Figure S6.Model MO (P waves).Figure S7.Model AK (S waves).Figure S8.Model MO (S waves).Figure S9.A comparison between predicted group traveltimes from final model AK, MO and SW with starting model SW, for a group of surface wave paths.Figure S10.We generate noise-free synthetic P-wave traveltimes through starting model MO (a) and SW (c), and recover these traveltimes against starting model AK (http://gji.oxfordjournals.org/lookup/suppl/doi:10.1093/gji/ggt476/-/DC1).

Figure S7 .
Figure S1.Histograms showing initial and final residuals for P and S models AK, MO and SW.FigureS2.An additional 5939 P-wave arrivals from 152 events (green squares) were used to test the fit to the various models (Fig.S1).FigureS3.Example of a waveform for a typical P-wave event.Figure S4.Example of a waveform for a typical S-wave event; traces are aligned according to relative delays determined by crosscorrelation.Figure S5.Model AK (P waves).Figure S6.Model MO (P waves).Figure S7.Model AK (S waves).Figure S8.Model MO (S waves).Figure S9.A comparison between predicted group traveltimes from final model AK, MO and SW with starting model SW, for a group of surface wave paths.Figure S10.We generate noise-free synthetic P-wave traveltimes through starting model MO (a) and SW (c), and recover these traveltimes against starting model AK (http://gji.oxfordjournals.org/lookup/suppl/doi:10.1093/gji/ggt476/-/DC1).

Figure S8 .
Figure S1.Histograms showing initial and final residuals for P and S models AK, MO and SW.FigureS2.An additional 5939 P-wave arrivals from 152 events (green squares) were used to test the fit to the various models (Fig.S1).FigureS3.Example of a waveform for a typical P-wave event.Figure S4.Example of a waveform for a typical S-wave event; traces are aligned according to relative delays determined by crosscorrelation.Figure S5.Model AK (P waves).Figure S6.Model MO (P waves).Figure S7.Model AK (S waves).Figure S8.Model MO (S waves).Figure S9.A comparison between predicted group traveltimes from final model AK, MO and SW with starting model SW, for a group of surface wave paths.Figure S10.We generate noise-free synthetic P-wave traveltimes through starting model MO (a) and SW (c), and recover these traveltimes against starting model AK (http://gji.oxfordjournals.org/lookup/suppl/doi:10.1093/gji/ggt476/-/DC1).
Figure S1.Histograms showing initial and final residuals for P and S models AK, MO and SW.FigureS2.An additional 5939 P-wave arrivals from 152 events (green squares) were used to test the fit to the various models (Fig.S1).FigureS3.Example of a waveform for a typical P-wave event.Figure S4.Example of a waveform for a typical S-wave event; traces are aligned according to relative delays determined by crosscorrelation.Figure S5.Model AK (P waves).Figure S6.Model MO (P waves).Figure S7.Model AK (S waves).Figure S8.Model MO (S waves).Figure S9.A comparison between predicted group traveltimes from final model AK, MO and SW with starting model SW, for a group of surface wave paths.Figure S10.We generate noise-free synthetic P-wave traveltimes through starting model MO (a) and SW (c), and recover these traveltimes against starting model AK (http://gji.oxfordjournals.org/lookup/suppl/doi:10.1093/gji/ggt476/-/DC1).

Figure S10 .
Figure S1.Histograms showing initial and final residuals for P and S models AK, MO and SW.FigureS2.An additional 5939 P-wave arrivals from 152 events (green squares) were used to test the fit to the various models (Fig.S1).FigureS3.Example of a waveform for a typical P-wave event.Figure S4.Example of a waveform for a typical S-wave event; traces are aligned according to relative delays determined by crosscorrelation.Figure S5.Model AK (P waves).Figure S6.Model MO (P waves).Figure S7.Model AK (S waves).Figure S8.Model MO (S waves).Figure S9.A comparison between predicted group traveltimes from final model AK, MO and SW with starting model SW, for a group of surface wave paths.Figure S10.We generate noise-free synthetic P-wave traveltimes through starting model MO (a) and SW (c), and recover these traveltimes against starting model AK (http://gji.oxfordjournals.org/lookup/suppl/doi:10.1093/gji/ggt476/-/DC1).