# A global horizontal shear velocity model of the upper mantle from multimode Love wave measurements

Tak Ho, Keith Priestley, Eric Debayle

► To cite this version:

Tak Ho, Keith Priestley, Eric Debayle. A global horizontal shear velocity model of the upper mantle from multimode Love wave measurements. Geophysical Journal International, 2016, 207 (1), pp.542-561. 10.1093/gji/ggw292. hal-02331975

HAL Id: hal-02331975

https://univ-lyon1.hal.science/hal-02331975

Submitted on 8 Sep 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# A global horizontal shear velocity model of the upper mantle from multimode Love wave measurements

Tak Ho,<fnref n="1" /> Keith Priestley<fnref n="1" /> and Eric Debayle<fnref n="2" />

Accepted 2016 July 29. Received 2016 May 23; in original form 2015 December 4

## SUMMARY

Surface wave studies in the 1960s provided the first indication that the upper mantle was radially anisotropic. Resolving the anisotropic structure is important because it may yield information on deformation and flow patterns in the upper mantle. The existing radially anisotropic models are in poor agreement. Rayleigh waves have been studied extensively and recent models show general agreement. Less work has focused on Love waves and the models that do exist are less well-constrained than are Rayleigh wave models, suggesting it is the Love wave models that are responsible for the poor agreement in the radially anisotropic structure of the upper mantle. We have adapted the waveform inversion procedure of Debayle & Ricard to extract propagation information for the fundamental mode and up to the fifth overtone from Love waveforms in the 50–250 s period range. We have tomographically inverted these results for a mantle horizontal shear wave-speed model ($\beta_h(z)$) to transition zone depths. We include azimuthal anisotropy (2$\theta$ and 4$\theta$ terms) in the tomography, but in this paper we discuss only the isotropic $\beta_h(z)$ structure. The data set is significantly larger, almost 500 000 Love waveforms, than previously published Love wave data sets and provides ~17 000 000 constraints on the upper-mantle $\beta_h(z)$ structure. Sensitivity and resolution tests show that the horizontal resolution of the model is on the order of 800–1000 km to transition zone depths. The high wave-speed roots beneath the oldest parts of the continents appear to extend deeper for $\beta_h(z)$ than for $\beta_v(z)$ as in previous $\beta_h(z)$ models, but the resolution tests indicate that at least parts of these features could be artefacts. The low wave speeds beneath the mid-ocean ridges fade by ~150 km depth except for the upper mantle beneath the East Pacific Rise which remains slow to ~250 km depth. The resolution tests suggest that the low wave speeds at deeper depths beneath the East Pacific Rise are not solely due to vertical smearing of shallow, low wave speeds. Four prominent, low wave-speed features occur at transition zone depths—one aligned along the East African Rift, one centred south of the Indian peninsula, one located south of New Zealand and one in the south Pacific Ocean coinciding with the location of the South Pacific Superswell. The low wave-speed features south of New Zealand and south of the Indian peninsula correspond spatially with the two largest negative geoid lows on Earth.

**Key words:** Surface waves and free oscillations; Seismic anisotropy; Seismic tomography; Cratons.

## 1 INTRODUCTION

Surface wave studies in the 1960s (e.g. Anderson 1961; Aki & Kaminuma 1963; McEvilly 1964) provided the first indications that the upper mantle was radially anisotropic; that is, the horizontal shear wave velocity ($\beta_h(z)$) is greater than the vertical shear wave velocity ($\beta_v(z)$). Resolving mantle anisotropic structure is important because it may yield information on deformation patterns and flow in the upper mantle (e.g. Montagner 2007; Becker et al. 2008). Much effort has gone into extracting the radially anisotropic structure of

the upper mantle (e.g. Montagner & Tantimoto 1991; Gung et al. 2003; Panning & Romanowicz 2006; Kustowski et al. 2008; Nettles & Dziewoński 2008; Lekić & Romanowicz 2011; French et al. 2013; Auer et al. 2014; Moulik & Ekström 2014) and the large-scale features of existing models do match between the models; however, the finer scale structures show poorer agreement (Chang et al. 2014).

Rayleigh waves have been studied extensively (e.g. Debayle et al. 2005; Lekić & Romanowicz 2011; Ritsema et al. 2011; Debayle & Ricard 2012; French et al. 2013; Priestley & McKenzie

<footnote n="1">¹Department of Earth Sciences, Bullard Laboratories, University of Cambridge, Cambridge CB3 0EZ, United Kingdom. E-mail: th332@cam.ac.uk</footnote>

<footnote n="2">²Laboratoire de Géologie de Lyon: Terre, Planètes et Environnement, Université Claude Bernard (Lyon1) and ENS Lyon, 6 rue Raphaël Dubois, F-69622 Villeurbanne, France</footnote>
2013; Schaeffer 2013; Auer et al. 2014) and most recent models show reasonable agreement. Less work has focused on Love waves and the models that do exist (e.g. Masters et al. 2000; Gung et al. 2003; Zhou et al. 2005) are less well-constrained than are Rayleigh wave models for several reasons. First, seismic data collected from horizontal seismographs are intrinsically more noisy than are data recorded on vertical seismographs and the number of Love wave measurements is significantly smaller than those for Rayleigh waves; hence, the Love wave path coverage is less dense compared to the Rayleigh wave path coverage, resulting in a lower resolution for the Love wave model. Second, for the same frequency, the Love wave fundamental mode, which has been measured in most past studies, has a shallower depth sampling than does the same frequency fundamental mode Rayleigh wave. As a result, Love waves are more sensitive to the crustal structure than are Rayleigh waves; hence, the upper-mantle $\beta_h(z)$ wave-speed models are more susceptible to contamination from inaccuracies in the crustal model than are $\beta_v(z)$ models. In addition, measurements for lower frequency fundamental mode Love waves, compared to Rayleigh waves, are necessary to cover the same depth range in the Earth and these longer-period measurements are not often available.

Because of these difficulties in making Love wave measurements, the resolution of the $\beta_h(z)$ models is lower than that for the $\beta_v(z)$ models and, therefore, the radial anisotropy models are only as good as the $\beta_h(z)$ models. Love and Rayleigh wave dispersion maps for the fundamental mode exist (e.g. Trampert & Woodhouse 1995, 2000; Larson & Ekström 2001; Visser et al. 2007; Ekström 2011; Pasyanos et al. 2014) and Rayleigh wave maps exist for some higher modes (van Heijst & Woodhouse 1997, 1999; Shapiro & Ritzwoller 2002; Beucler et al. 2003; Visser et al. 2007, 2008). Less work has been done on inverting Love wave measurements for $\beta_h(z)$ models (Li & Romanowicz 1996; Mégnin & Romanowicz 2000; Zhou et al. 2005) and most existing $\beta_h(z)$ upper-mantle models are constrained using only fundamental mode measurements.

This paper focuses on the isotropic $\beta_h(z)$ structure of the upper mantle. The depth range between about 200 km where the sensitivity of the fundamental mode becomes weak and the transition zone where body wave sensitivity increases, has been a particularly opaque region for seismologists. Including higher mode information allows stronger constraints on the $\beta_h(z)$ structure in this depth range. The goal of this paper, improving knowledge of the upper-mantle $\beta_h(z)$ structure, is a prerequisite for better defining the radially anisotropic structure of the upper mantle.

In this study, we have adapted the Rayleigh wave inversion routines of Debayle & Ricard (2012) so as to invert Love waveforms. We show that this procedure can recover velocity structure of the upper mantle to transition zone depths. The path-average $\beta_h(z)$ models are then combined into a 3-D model using the regionalization formulation of Debayle & Sambridge (2004). Here, we discuss the Love wave data set, the analysis procedure, the effect of the crustal model, the depth sensitivity of the waveform inversion and resolution of the tomography. Although we have included both the $2\theta$ and $4\theta$ Love wave azimuthal terms (Montagner & Nataf 1986) in the tomographic inversion, we concentrate on the isotropic $\beta_h(z)$ structure of the Earth from Love wave inversions using information from the fundamental mode and up to the fifth overtone for the period range of 50–250 s from ~500 000 Love wave propagation paths.

# 2 BUILDING THE $V_{sh}$ UPPER-MANTLE MODEL

## 2.1 The data

The data set used in building the 3-D upper-mantle $\beta_h(z)$ model consists of seismograms obtained from the IRIS and EIDA Data Management Centers for the years 1976–2015 for the permanent stations of the Federation of Digital Seismographic Networks (FDSN) and a large number of temporary seismographs. The locations of the stations and events providing data for this study are plotted on Fig. 1. Coverage is best in the Northern Hemisphere and gets progressively worse in the Southern Hemisphere (Figs B1 and B2). We do not employ major arc-path data; we use seismograms with propagation path lengths between 10° and 140°.

## 2.2 Examining the horizontal component noise

We select seismograms for inversion by first examining the signal-to-noise ratio in a series of discrete period bands centred at 50, 75, 110, 165 and 250 s. We test the signal-to-noise ratio separately for the north and east components by comparing the amplitude $A_s$ within the group velocity window of 3–7 km s⁻¹ containing the surface wave train, with the amplitude $A_n$ in the group velocity window 2.5–3 km s⁻¹, for each period band. Waveforms with $A_s/A_n > 3$ for both horizontal components are accepted. If the waveform for one component has a lower signal-to-noise ratio for a period band, both components for that period band are rejected unless the high single-to-noise component is naturally rotated close to transverse. The gains of the two components are normalized and they are rotated to obtain the transverse component seismogram.

## 2.3 Extracting and fitting the secondary observables

Extracting the secondary observables (Cara & Lévêque 1987) and inverting them for the path average velocity structure follows closely the method employed by Debayle & Ricard (2012) but adapted to transverse component seismograms. The secondary observables are measured using cross-correlation. Two synthetic seismograms are calculated using an updated version of the DISPER80 code (Saito 1988). The first is a complete synthetic seismogram computed for a model which is then updated at each iteration of the inversion. This synthetic seismogram and the observed seismogram are cross-correlated with a single mode synthetic seismogram computed for a reference model (Fig. 2). The envelopes of these cross-correlograms are filtered in the frequency passbands determined suitable in the noise evaluation. The secondary variables picked are defined by three samples, one at the peak of the envelope and two on either side of the peak. The inversion minimizes the difference between these observables for the complete synthetic seismogram of the updated model and observed seismogram. The $\beta_h(z)$ model is updated by adjusting the current model using the partial derivatives (Fig. 3) with respect to the model parameter $\beta_h(z)$ and the attenuation calculated using the equations in Takeuchi & Saito (1972). The instantaneous phase of the cross-correlation is included in the inversion to adjust the waveform fit.

The starting wave-speed model for the inversion and the reference wave-speed model for the cross-correlation are identical. They consist of an *a priori* crustal model extracted from CRUST1.0 for each path overlying an upper-mantle wave-speed model based on PREM (Dziewoński & Anderson 1981) but with the 220 km discontinuities smoothed (Fig. 2). The reference model and the crust of the
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Locations of the stations (blue triangles) and the events (red stars) whose seismograms were used in the inversion for the $β_h(z)$ upper-mantle model.</figcaption></figure>

inversion model remain fixed for each iteration of the inversion but the mantle portion of the inversion model is updated at each stage of the inversion. The source excitation required for generating the synthetics is calculated using the expressions in Cara (1978) and the moment tensor source description is taken from the global CMT catalogue (Dziewoński et al. 1981; Ekström et al. 2012).

The criteria for successfully fitting the waveform are the same as given in Debayle & Ricard (2012). The waveform fit is successful if the secondary observables are fit, the inversion converges to a unique model, and the observed Love waveform is well matched by the final inversion synthetic waveform.

## 2.4 Extracting structure of the upper mantle

We extract information for the fundamental mode and the first five overtones from the seismograms. Including overtones in the tomography not only constrains deeper structure better than does the fundamental mode, the overtones improve resolution at shallower depths. An advantage in using the Debayle & Ricard (2012) implementation of the Cara & Lévêque (1987) method over other automated waveform inversion procedures is that it does not require separating the fundamental mode and overtones into unique time windows before making the wave speed measurements since the observed seismogram is modelled as a sum of interfering modes.

We evaluate the ability to extract information for the upper-mantle $β_h(z)$ structure by carrying out tests using synthetic seismograms. We compute sets of synthetic seismograms with a variety of sources in three depth ranges: 10–50, 50–350 and 350–600 km. We add real ground noise extracted at random time periods from recordings of the North American station CCM of the US national network. We invert the synthetic seismograms including noise once using the CMT fault description, source-time function and location, and once with the CMT fault description, source-time function but with the PDE location. In some cases the CMT location

and origin time can be significantly different from the PDE location and origin time. We test the inversion results with the CMT and PDE earthquake locations to evaluate the influence of potential source mislocation. The results for these tests are shown in Figs 4(a)–(f). For each set of tests we have recovered the $β_h(z)$ model which was used in computing the synthetic seismogram to at least 700 km depth. The difference in mantle structure between using the CMT and PDE locations in the inversion is small. Maggi & Priestley (2005) showed that artefacts due to event mislocation in the 3-D upper-mantle models were further minimized during the tomography.

In constructing the waveform inversion starting and reference models we have extracted the path average crustal model from CRUST1.0. Although CRUST1.0 contains more crustal constraints than older global crustal models, there are likely to be geographic regions where the crustal structure of CRUST1.0 is inaccurate. Bozdağ & Trampert (2008) have shown that errors in the crustal correction can have a detrimental effect on the $β_h(z)$ model. To assess the effect of errors in the *a priori* crust of the inversion starting and reference models on the final waveform inversion results, we re-invert the synthetic seismograms discussed in the proceeding paragraphs but choose the *a priori* crust with approximately 5 per cent difference in wave speed (Figs 4g–h) or up to 10 km difference in thickness from that used to compute the synthetic seismograms (Fig. 4i). These are extremes as it is unlikely that the crustal properties would be in error by as much as those assumed over the whole of the surface wave propagation path. Even for the most extreme crustal error in the inversion starting and reference models there is little effect on the resulting mantle structure below ~100 km depth. These tests demonstrate the level of error we need to be aware of in interpreting the $β_h(z)$ model, but they do not provide insight into lateral variations in resolution. In Appendix B, we discuss additional tests to evaluate the resolution of the $β_h(z)$ model.

Each waveform inversion yields a path-average wave speed model and the path-average phase and group velocity dispersion. In
<figure><img src="image_4.png" /><figcaption>**Figure 2.** The radially anisotropic reference wave-speed model used in the waveform inversion for the path average wave-speed models. This model is based on PREM (Dziewoński & Anderson 1981) but with the 220 km discontinuity smoothed. The 220 km discontinuity is smoothed in our reference model; if this was not the case this discontinuity would persist in the waveform inversion and be carried through to the tomographic model as has been done in some tomographic mantle models (e.g. Ritsema et al. 2011).</figcaption></figure>

addition, the waveform inversion provides an *a posteriori* covariance matrix $C_m$,

$$C_m = (I - R)C_{m0}$$

where $I$ is the identity matrix, $R$ is the resolution matrix and $C_{m0}$ is the *a priori* covariance matrix. The *a priori* and *a posteriori* errors are obtained from the square root of the diagonal terms of $C_{m0}$ and $C_m$ respectively. The *a posteriori* errors are close to the *a priori* errors when the resolution is poor ($R$ close to zero). The *a posteriori* errors decrease when resolution increases and tend towards zero for perfect resolution ($R = I$). If the *a posteriori* error is greater than 80 per cent of the *a priori* error (0.05 km s⁻¹ in this study), we reject that path for that depth as having poor resolution, thus ensuring that our model is not contaminated by low resolution paths.

For the assembled data set, more than 475 000 independent paths passed the waveform fitting and provide more than 17 000 000 constraints on upper-mantle structure. This is the total number of envelope picks used in the inversions (which consist of three picks per lobe on the cross-correlogram) in addition to the instantaneous phase information for all the paths.

## 3 THE $\beta_h(z)$ TOMOGRAPHIC MODEL

In the second stage of the analysis we combine the path-average wave speed models at each depth to produce a 3-D model of the wave speed and azimuthal anisotropy using the regionalization formulation of Debayle & Sambridge (2004) based on Montagner (1986). We refer to this model as CAM2016SH.

### 3.1 Clustering the paths

We first cluster paths which have epicentres within a defined clustering radius and with receivers within a defined clustering radius, assuming that these paths sample the same Earth structure. We use a cluster radius of 200 km for both sources and receivers. This is smaller for both events and stations than the size of the correlation length $L_c$ (explained below). From each cluster, we select the best-resolved path with the highest trace for the resolution matrix $R$ and use the *a posteriori* error to weight the path in the tomographic inversion.

### 3.2 The regionalization

In the regionalization, we invert for the 3-D wave speed and azimuthal anisotropy (2θ and 4θ terms) model consistent with the observed Love waveforms. For regularization in the tomographic inversion we apply a correlation length $L_c$ which imposes a correlation between neighbouring points using a Gaussian function. The azimuthal variation for Love waves depends strongly on sin 4θ and cos 4θ where θ is the azimuth of propagation (Smith & Dahlen 1973; Montagner & Nataf 1986). Good azimuthal path coverage is required to resolve the 4θ terms and this azimuthal variation will average out more easily when compared to Rayleigh waves which primarily have a 2θ dependence. In this study, we set $L_c$ to 400 km in the shallower part of the mantle, to 600 km in the mid-upper mantle and to 800 km in the deeper part of the mantle. A larger $L_c$ in the deeper mantle accounts for the decrease in path coverage with depth. We examine the effect of the $L_c$ length in Appendix A. We solve for azimuthal anisotropy but in this paper concentrate on discussing the isotropic parts of the $\beta_h(z)$ model.

## 4 THE HORIZONTAL SHEAR VELOCITY MODEL CAM2016SH

Fig. 5 shows maps of the $\beta_h(z)$ model at a number of upper-mantle depths. We discuss the resolution of CAM2016SH in Appendix B. Those tests and other tests not included show that the horizontal resolution in our model varies with geographic location but is on the order of 800–1000 km, meaning that, on average, two positive wave-speed features separated by 800–1000 km are resolved. The results shown in Fig. 5 need to be viewed alongside the resolution figures in Appendix B and sensitivity tests shown in Fig. 4. There is strong visual correlation between the variation in mantle wave speed and surface tectonics at upper-mantle depths (≈200 km), but this correlation decreases with increasing depth.

In the 100 km-depth map the Mid-Atlantic Ocean, Indian Ocean and southern Pacific Ocean spreading ridges appear as a restricted, slow wave-speed feature centred on the bathymetric expression of the ridges. These slow wave speeds give way to fast wave speeds away from the ridges over a few hundred kilometres as the plate forms (Parsons & Sclater 1977). The upper mantle beneath the East Pacific Rise shows as a broad, strong slow wave speed feature in the
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Partial derivatives calculated for the reference model (Fig. 2) for the different modes and periods: 50 s (blue), 100 s (green), 140 s (red), 200 s (turquoise) and 240 s (black).</figcaption></figure>

100 km depth map with a transition from slow to fast wave speeds toward the older portion of the Pacific Plate, similar to what has been observed in Rayleigh wave models of the Pacific Plate (Forsyth 1977; Zhang & Tanimoto 1991; Ritzwoller et al. 2004; Priestley & McKenzie 2006; Maggi et al. 2006; Priestley & McKenzie 2013). By ~200 km depth the low wave speeds associated with the Mid-Atlantic Ocean, Indian Ocean and southern Pacific Ocean spreading ridges largely disappear but the low wave speeds associated with the

East Pacific Rise persist to as much as 250 km depth. Portions of the Indian Ocean ridge fall in the poorest resolved part of CAM2016SH but the East Pacific Rise is in a well-resolved region (Figs B1 and B2) and the resolution tests shown in Figs B3 and B4 suggests that the deeper, slow wave speeds observed beneath the East Pacific Rise are not due solely to vertical smearing. There are some vestiges of the high wave-speed features beneath the northwest Pacific at deeper depths. The back arc basins of the western Pacific and the South
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Subplots (a)–(f) evaluate the ability of the Cara-Lèveûque waveform inversion procedure to extract $\beta_h(z)$ mantle structure. Synthetic seismograms were computed for various CMT source mechanisms and depths [(a,d) 10–50 km depth; (b,e) 50–300 km depth; (c,f) 350–600 km depth] for an input wave-speed model shown as the red line. Recorded ground noise was added to the synthetic seismograms and then the synthetic waveforms were inverted using the procedure outlined in Section 2.3 using a starting and reference model (black line) different from the wave-speed model used to compute the synthetic seismograms. The inversion results show that the method does satisfactorily retrieve structure to transition zone depths. The main discrepancies are in the region of the 410 and 670 km transition zone discontinuities. This is a result of the vertical correlation length imposed in the waveform inversion. Subplots (g)–(i) show the effect of using an incorrect *a priori* crustal correction. Subplot (g) shows the crust for the reference and starting model (black) with a different wave-speed structure from that used in computing the synthetic seismograms (red). The structure of the crust is fixed during the waveform inversion procedure. Subplot (h) shows the inversion result when the incorrect crustal wave speed is used for the reference and starting models. Subplot (i) shows the effect on the mantle structure when the crustal thickness is in error. The variation in crustal thickness is not shown as we have merely stretched the crustal layers and replaced the changes in the uppermost mantle layer with crustal material. In both cases, for an incorrect *a priori* crustal model, the main error in the inversion results is at depths shallower than 100 km. The test of the crustal correction is an extreme as it is unlikely that CRUST1.0 is in error by as much as we have used in these tests over the whole propagation path.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 5.** $\beta_h(z)$ model CAM2016SH plotted at eight depths. Velocities are plotted as perturbations to the average velocity at that depth (denoted above the individual maps)—±10 per cent in the upper 250 km and ±3.5 per cent for 350 and deeper depths. This change in scale is qualitatively justified as the resolution test shows that there is a decrease in the recovered amplitude with increasing depth. The maps are for the tomographic inversion with a correlation length of 400 km for depths 150 km or less, 600 km for the 200 and 250 km depth maps and 800 km at deeper depths. Plate margins are plotted as green lines and hot spots are denoted by red circles.</figcaption></figure>
Pacific super-swell (McNutt & Judge 1990) show as slow wave-speed features in the 100 and 150 km depth map and are still visible in the 200 km-depth map but disappear by 250 km depth.

The high wave-speed signature of the cratons is well defined in the shallow depth maps (Figs 5a and b) and for some cratons the high wave speeds persist down to ~350 km. The exceptions are the South Indian shield and the North China craton, neither of which has a deep, high wave-speed $\beta_h(z)$ root. The cratonic roots extend somewhat deeper in the CAM2016SH models than they do in the comparable $\beta_v(z)$ models (Debayle & Ricard 2012; Priestley & McKenzie 2013). Gung et al. (2003) suggest that roots of continental shields appear deeper in $\beta_h(z)$ models than they do in $\beta_v(z)$ models due to radial anisotropy arising from present-day asthenospheric flow aligning the mineral fabric leading to the strong fast $\beta_h(z)$ signature. Without the resolution tests (Figs B5–B8) we might draw the same conclusion. CAM2016SH is constrained by almost an order of magnitude more data than the model discussed in Gung et al. (2003) and the results of the resolution tests in Appendix B suggest that deeper high wave-speed features beneath some of the cratons may be an artefact and call into question this earlier idea that the deeper $\beta_h(z)$ signature compared to $\beta_v(z)$ is due to radial anisotropy resulting from present-day asthenospheric flow.

The underthrusting of the Indian Plate beneath Tibet in the Eurasian collision zone is distinguishable in the model at depths down to ~350 km. The high wave speeds of the South Indian Shield are clear at 100 km depth but are less clear by 150 km depth and have largely disappeared by 200 km depth, much shallower than most of the other cratons. This is consistent with Rayleigh wave studies which show that the South Indian Shield has a relatively shallow root (Mitra et al. 2006; Priestley et al. 2008). Interestingly, the synthetic test (Fig. B7) for India does not show as strong a high wave speed artefact in the upper mantle as observed in the tests for the other cratons. The alignment of the mineral fabric proposed by Gung et al. (2003) is related to the strains induced by the moving plates and this in turn is related to the translation speed of the plate. Australia and India are the two fastest moving plates containing continents and while the Australian craton does show the deep, high wave-speed root, the India plate – moving almost as fast – does not give further credence to the possibility that the high wave speeds seen at greater depths beneath most cratons in the $\beta_h(z)$ models might be related to model artefacts. North of Tibet high wave speeds occur at 100 km depth beneath the Tien Shan Mountains. The slow Tibetan wave speeds to the north of India seen in the 100 km-depth map become fast by 150 km depth and extend to at least 350 km depth, similar to what is observed in many Rayleigh wave models (e.g. Priestley et al. 2008; Ritsema et al. 2011).

At 100 km depth the region of the Afar plume and the region surrounding the Red Sea are slow and this slow feature extends northward to the Caucasus Mountains and then across the Anatolian Plateau. By 150 km depth the region beneath Anatolia is fast but the low wave-speed feature extending from Afar to the Caucasus persists to 200 km depth. At deeper depths this low velocity feature which appears to extend into the upper-mantle transition zone becomes more restricted to the region beneath the East African Rift. The most pronounced feature at deeper depths (>350 km) are the strong, low wave-speed anomalies roughly coinciding with the South Pacific Superswell and beneath the triple junction South of New Zealand. The low wave-speed feature spatially associated with the South Pacific Superswell extends through the upper-mantle transition zone but is more restricted in area as suggested by Niu et al. (2002), rather than a feature of broad lateral extent as postulated at shallower depths (McNutt & Judge 1990). There is a broad low

wave-speed feature within the transition zone which closely coincides with the geoid low between New Zealand and Antarctica (Sutherland et al. 2010). The other notable transition zone feature is the large low wave-speed region south of India coinciding with the Indian Ocean geoid low, the most prominent negative geoid feature on Earth.

## 5 CORRELATION WITH OTHER $\beta_h(z)$ MODELS

A number of radially anisotropic upper-mantle models have been developed (e.g. Zhou et al. 2005; Kustowski et al. 2008; Panning et al. 2010; Lekić & Romanowicz 2011; French et al. 2013; Auer et al. 2014; Chang et al. 2014; Moulik & Ekström 2014). Here, we compare our $\beta_h(z)$ upper-mantle model CAM2016SH with four recent $\beta_h(z)$ upper-mantle models: FFSW1 (Zhou et al. 2005), S362ANI (Kustowski et al. 2008), SAW642ANb (Panning et al. 2010) and SEMum2 (French et al. 2013). The first of these is a $\beta_h(z)$ model whereas the later three are radial anisotropic ($\xi$) model from which we have extracted the $\beta_h(z)$ component using the relationships between the $\xi$ model and $\beta_h(z)$ given in those papers. There have not been any detailed discussions of the resolution for these models. The resolution of CAM2016SH is discussed in Appendix B.

Model FFSW1 (Zhou et al. 2005) employs finite-frequency tomography using 3-D Born sensitivity kernels and inverts ~12 000 Love and Rayleigh dispersion fundamental measurements in the 5–15 mHz frequency range. The dispersion measurements were made for minor arc (G1,R1), major arc (G2,R2) as well as multiple orbits (G3,G4,R3,R4). Model S362ANI (Kustowski et al. 2008) is calculated from surface wave and body wave measurements. The upper mantle is constrained mainly by the surface wave information. Fundamental mode Love waves were inverted using a period range of 35–150 s. A total of ~55 000 paths were employed in their study. CRUST2.0 (Bassin et al. 2000) is used for crustal corrections. The model is expanded into 362 spherical splines, equivalent to a spherical harmonic expansion up to degree 18.

SAW642ANb (Panning et al. 2010) is derived using a non-linear asymptotic coupling theory (NACT; Li & Romanowicz 1995) including fundamental and higher mode information. This theory is based on normal modes which attempt to account for mode coupling and finite frequency effects of surface waves and body waves. Crustal properties are included in the inversion. The shortest period included in the SAW642ANb model is 60 s. The longest period is dependent on the magnitude of the earthquake and ranges from 220 to 3600 s. However, the path coverage in this model is low with only 21 101 Love wave paths. The model is expanded up to degree 24. Finally, model SEMum2 is a ‘hybrid’ spectral element method model. Unlike other spectral element models which use a full 3-D numerical calculation for the synthetic seismogram and Fréchet kernels (e.g. Fichtner et al. 2008; Bozdağ et al. 2011; Colli et al. 2013), SEMum2 uses a spectral element method to calculate synthetic seismograms and NACT to calculate the sensitivity kernels, making it faster than conventional full waveform inversion techniques. Short period group velocity dispersion maps are included to constrain the crustal structure. However, the method is still computationally intensive and a total of 13 192 Love wave fundamental mode wavepackets and ~14 000 overtone wavepackets were inverted in a period range from 60 s. The anisotropic part of SEMum2 is expanded up to degree 24.

Chang et al. (2014) found that there are more discrepancies between the published radial anisotropy models than there are
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Correlations computed between CAM2016SH and models FFSW1 (green), S362ANI (in red), SAW642ANb (in blue) and SEMum2 (in grey). Significance levels of 95 per cent (dotted) and 66 per cent (dashed) are plotted.</figcaption></figure>

similarities. As $\beta_v(z)$ upper-mantle models show a very good agreement, in some cases up to wavelengths of 650 km in the uppermost mantle (Debayle & Ricard 2012), it is likely that the observed discrepancies in radial anisotropy comes from differences in $\beta_h(z)$ models. We compare our $\beta_h(z)$ model (A) with published $\beta_h(z)$ models (B) using the following relationship:

$$C(l) = \sqrt{\frac{\sum_{m=-l}^{l} A_l^m B_l^{m*}}{S_A(l)S_B(l)}} \quad (1)$$

where $S_A(l) = \sum_{m=-l}^{l} A_l^m A_l^{m*}$ is the spectra for model A, $A_l^m$ is the spherical harmonic coefficient at degree $l$ and azimuthal order $m$ for model A, and likewise for model B. We calculate the correlation at depths of 100, 150, 200, 250, 350, 450, 550 and 610 km.

Fig. 6 shows the correlation between CAM2016SH and the published models. Due to the parametrization used in these models, we truncate the analysis to degree 20 for FFSW1, to degree 18 for S362ANI and degree 24 for SAW642ANb and SEMum2. At 100 km depth, CAM2016SH and FFSW1 correlate above the 95 per cent confidence level up to ~degree 15; CAM2016SH and S362ANI

correlate above the 95 per cent confidence level to ~degree 18; CAM2016SH and SAW642ANb correlate above the 95 per cent level to ~degree 20; and CAM2016SH and SEMum2 correlate above 95 per cent up to ~degree 24. At 150 km depth CAM2016SH correlates with all of the published models above 66 per cent below ~degree 15. SEMum2 and SAW642ANb remain mostly above the 66 per cent confidence level up to degree 24, S362ANI is correlated with CAM2016SH up to degree 18 and FFSW1 falls below the 66 per cent confidence level above degree 15. At these depths the structure in all of the models is primarily constrained by fundamental mode surface waves and it is not that surprising that the models correlate well with each other at lower degrees. At higher harmonic degree the spatial resolution depends on path coverage, the theory used, the parametrization and the regularization.

At 200 and 250 km the fundamental mode sensitivity drops off and the sensitivity of the higher modes start to dominate. At these depths CAM2016SH correlates with S362ANI above 95 per cent below degree 8-10. SAW642ANb and SEMum2 correlate with CAM2016SH above 66 per cent below degree 14. At 350 km, FFSW1 remains above the 66 per cent correlation level up to
degree 11. S362ANI, SAW682ANb and SEMum2 do not correlate well with CAM2016SH at this depth. Finally at depths greater than 400 km, no model correlates well with CAM2016SH above degree 5.

In Appendix C, we show the correlation between the four models FFSW1, S362ANI, SAW642ANb and SEMum2. All models have high correlation in the 100–150 km depth range. However, at depths greater than 200 km, this correlation drops significantly, far shallower than the correlation between $\beta_v(z)$ models (Debayle & Ricard 2012). Below 200 km, S363ANI correlates with SAW642ANb and SEMum2 above the 66 per cent confidence level up to degree ~5–10 at all depths. SAW642ANb correlates with SEMum2 up to ~degree 14. However, FFSW1 correlates poorly with all models below 250 km. The poor correlation below ~200 km depth suggests that it is the fundamental mode which primarily provides the model constraints and that $\beta_h(z)$ in these mantle models is still poorly constrained at depths greater than 200 km. By including the higher modes and longer period surface waves, we have greater sensitivity and resolution to the deeper upper-mantle structures. We show in Section 2.4 and Appendix B that the data employed in building CAM2016SH allow us to isolate seismic anomalies from the shallower and deeper structure, even when they are located within the transition zone.

# 6 SUMMARY AND CONCLUSIONS

This paper presents a new global $\beta_h(z)$ model CAM2016SH for the upper mantle obtained from a new data set analysed by a technique not previously applied to Love waves. To derived CAM2016SH we have adapted the procedure of Debayle & Ricard (2012), based on the waveform inversion technique of Cara & Lévêque (1987), so as to analyse horizontal component seismograms. We have applied this analysis method to a new, large, multimode Love wave data set. We have carried out extensive resolution tests for CAM2016SH and these tests suggest that we are able to resolve 800–100 km features to transition zone depths.

At shallow ($\lesssim 200$ km) depths the features in CAM2016SH match well with surface tectonics and geological features and correspond well with features seen in other $\beta_h(z)$ models and in $\beta_v(z)$ models— low wave speeds for the mid-ocean ridge system, backarc basins and tectonically active regions on the continents, and high wave speed areas of the continents and the old ocean basins. The low wave speeds beneath the mid-ocean ridges fade by ~150 km depth except for the upper mantle beneath the East Pacific Rise which remains slow to ~250 km depth. The resolution tests suggest that these low wave speeds beneath the East Pacific Rise are not solely due to smearing of shallow, low wave speeds to deeper depths.

The situation is different at deeper depths. A comparison of CAM2016SH with a number of published Love wave-based models and a comparison between the published models among themselves shows poor agreement in structures at depths greater than 250 km. In CAM2016SH at deeper depths, to ~350 km, we observe fast wave speeds persisting beneath many but not all of the continental shields. This is deeper than the high wave-speed signature of the shields observed in the $\beta_v(z)$ upper-mantle models. The deeper high wave-speed signature of the shields in CAM2016SH fits well with observations made by Gung et al. (2003) who suggest that this may be caused by radial anisotropy that develops from present day flow beneath these cratons. However, CAM2016SH is based on a much larger data density than the previously published upper-mantle models based on Love wave data, and resolution tests for CAM2016SH

show that the high wave-speed signature of the shields observed in the $\beta_h(z)$ models may be an artefact.

Four prominent, low wave-speed features occur in CAM2016SH at transition zone depths—one aligned along the East African Rift, one centred south of the Indian peninsula, one located south of New Zealand and one in the south Pacific Ocean coinciding with the South Pacific Superswell. The low wave-speed features south of New Zealand and south of the Indian peninsula correspond spatially with the two largest negative geoid lows on the Earth.

# ACKNOWLEDGEMENTS

This research was partially funded by an EPSRC studentship, with CASE funding from Weston Geophysical. This work was supported by the French ANR SEISBLOB no ANR-11-BLANC-SIMI5-6-016-01. The data used in this study was provided by the IRIS and DMC programs—networks: 1A 1D, 1E, 2C, 2D, 2G, 3A, 3E, 3F, 4A, 5A, 5C, 6C, 6E, 7A, 7B, 7D, 7E, 7F, 7J, 8A, 8B, 9B, 9D, AC, AE, AF, AG, AI, AK, AO, AR, AS, AT, AU, AV, AY, AZ, BC, BE, BK, BL, C, C1, CB, CC, CD, CH, CI, CK, CM, CN, CO, CT, CU, CW, CY, CZ, DK, DR, DW, EI, ER, ET, G, GB, GE, GH, GI, GO, GS, GT, GY, HG, HL, HT, HV, HW, IC, II, IM, IN, IO, IP, IU, IW, JM, JP, KC, KN, KO, KR, KY, KZ, LB, LD, LI, LX, MC, MI, MN, MX, MY, N4, NA, ND, NE, NI, NJ, NM, NN, NO, NR, NU, NV, NX, NZ, OE, OK, ON, OO, OV, PB, PE, PL, PM, PN, PO, PR, PS, PT, PY, RE, RM, RS, RV, SC, SH, SN, SP, SR, SS, SV, TA, TC, TD, TJ, TM, TO, TS, TT, TW, UK, UO, US, UU, UW, VU, WC, WI, WU, WY, X1, X4, X5, X6, X8, X9, XA, XB, XC, XD, XE, XF, XG, XH, XI, XJ, XK, XL, XM, XN, XO, XP, XQ, XR, XS, XT, XU, XV, XW, XY, XZ, Y1, Y2, Y3, Y6, Y7, YA, YB, YC, YD, YE, YF, YG, YH, YI, YJ, YK, YL, YM, YN, YO, YP, YQ, YR, YS, YT, YU, YV, YW, YX, YY, Z1, Z2, Z3, Z4, Z5, Z6, Z8, Z9, ZA, ZB, ZC, ZD, ZF, ZG, ZI, ZJ, ZK, ZL, ZM, ZN, ZO, ZP, ZR, ZT, ZU, ZV, ZW, ZX. Figures were prepared using the GMT software (Wessel & Smith 1995).

# REFERENCES

Aki, K. & Kaminuma, K., 1963. Phase velocity of low waves in Japan (Part 1): Love waves from the Aleutian Shock of March 9, 1957, *Bull. Earthq. Res. Inst. Univ. Tokyo*, **41**, 243–259.

Anderson, D.L., 1961. Elastic wave propagation in layered anisotropic media, *J geophys. Res.*, **66**, 2953–2963.

Auer, L., Boschi, L., Becker, T., Nissen-Meyer, T. & Giardini, D., 2014. Savani: a variable resolution whole-mantle model of anisotropic shear velocity variations based on multiple data sets, *J geophys. Res.*, **119**(4), 3006–3034.

Bassin, C., Laske, G. & Masters, G., 2000. The current limits of resolution for surface wave tomography in North America, *EOS, Trans. Am. geophys. Un.*, **81**, Fall Meet. Suppl., Abstract F897.

Becker, T.W., Kustowski, B. & Ekström, G., 2008. Radial seismic anisotropy as a constraint for upper mantle rheology, *Earth planet. Sci. Lett.*, **267**, 213–227.

Beucler, E., Stutzmann, E. & Montagner, J.-P., 2003. Surface-wave higher mode phase velocity measurements, using a roller-coaster type algorithm, *Geophys. J. Int.*, **155**, 289–307.

Bozdağ, E. & Trampert, J., 2008. On crustal corrections in surface wave tomography, *Geophys. J. Int.*, **172**, 1066–1082.

Bozdağ, E., Trampert, J. & Tromp, J., 2011. Misfit functions for full wave-form inversion based on instantaneous phase and envelope measurements, *Geophys. J. Int.*, **185**, 845–870.

Cara, M., 1978. Etude du manteau supérieur à partir des harmoniques des ondes de surface, *PhD thesis*, Université Pierre et Marie Curie.

Cara, M. & Lévêque, J., 1987. Waveform inversion using secondary observables, *Geophys. Res. Lett.*, **14**, 1046–1049.
Chang, S.-J., Ferreira, A.M.G., Ritsema, J., van Heijst, H.J. & Woodhouse, J.H., 2014. Global radially anisotropic mantle structure from multiple datasets: a review, current challenges, and outlook, *Tectonophysics*, **617**, 1–19.

Colli, L., Fichtner, A. & Bunge, H.-P., 2013. Full waveform tomography of the upper mantle in the south atlantic region: imaging a westward fluxing shallow asthenosphere, *Tectonophysics*, **604**, 26–40.

Debayle, E. & Ricard, Y., 2012. A global shear velocity model of the upper mantle from fundamental and higher Rayleigh mode measurements, *J. Geophys. Res.*, **117**, doi:10.1029/2012JB009288.

Debayle, E. & Sambridge, M., 2004. Inversion of massive surface wave data sets: Model construction and resolution assessment, *J. geophys. Res.*, **109**, doi:10.1029/2003JB002652.

Debayle, E., Kennett, B. & Priestley, K., 2005. Global azimuthal seismic anisotropy and the unique plate-motion deformation of Australia, *Nature*, **433**, 509–512.

Dziewoński, A., Chou, T. & Woodhouse, J., 1981. Determination of earthquake source parameters from waveform data for studies of global and regional seismicity, *J. geophys. Res.*, **86**, 2825–2852.

Dziewoński, A.M. & Anderson, D., 1981. Preliminary Reference Earth Model, *Phys. Earth planet. Inter.*, **25**, 297–356.

Ekström, G., 2011. A global model of Love and Rayleigh surface wave dispersion and anisotropy, 25–250 s, *Geophys. J Int.*, **187**, 1668–1686.

Ekström, G., Nettles, M. & Dziewoński, A.M., 2012. The global CMT project 2004–2010: centroid-moment tensors for 13,017 earthquakes, *Phys. Earth planet. Inter.*, **200**, 1–9.

Fichtner, A., Kennett, B., Igel, H. & Bunge, H.-P., 2008. Theoretical background for continental- and global-scale full-waveform inversion in the time-frequency domain, *Geophys. J Int.*, **175**, 665–685.

Forsyth, D.W., 1977. The evolution of the upper mantle beneath mid-ocean ridges, *Tectonophysics*, **38**, 89–118.

French, S., Lekic, V. & Romanowicz, B., 2013. Waveform tomography reveals channeled flow at the base of the oceanic asthenosphere, *Science*, **342**, 227–230.

Gung, Y., Panning, M. & Romanowicz, B., 2003. Global anisotropy and the thickness of continents, *Nature*, **422**, 707–711.

Kustowski, B., Ekström, G. & Dziewoński, A., 2008. Anisotropic shear-wave velocity structure of the Earth’s mantle: a global model, *J. geophys. Res.*, **113**, 1978–2012.

Larson, E. & Ekström, G., 2001. Global models of surface wave group velocity, *Pure appl. Geophys.*, **158**, 1377–1399.

Lekić, V. & Romanowicz, B., 2011. Inferring upper-mantle structure by full waveform tomography with the spectral element method, *Geophys. J Int.*, **185**, 799–831.

Li, X.-D. & Romanowicz, B., 1995. Comparison of global waveform inversions with and without considering cross-branch modal coupling, *Geophys. J Int.*, **121**, 695–709.

Li, X.-D. & Romanowicz, B., 1996. Global mantle shear velocity model developed using nonlinear asymptotic coupling theory, *J. geophys. Res.*, **101**, 2245–2272.

Maggi, A. & Priestley, K., 2005. Surface waveform tomography of the Turkish–Iranian Plateau, *Geophys. J Int.*, **160**, 1068–1080.

Maggi, A., Debayle, E., Priestley, K. & Barruol, G., 2006. Multimode surface waveform tomography of the Pacific Ocean: a closer look at the lithospheric cooling signature, *Geophys. J Int.*, **166**, 1384–1397.

Masters, G., Laske, G., Bolton, H. & Dziewoński, A., 2000. The relative behavior of shear velocity, bulk sound speed, and compressional velocity in the mantle: implications for chemical and thermal structure, in *Earth’s Deep Interior: Mineral Physics and Tomography From the Atomic to the Global Scale*, pp. 63–87, eds Karato, S.-I., Forte, A., Liebermann, R., Masters, G. & Stixrude, L., American Geophysical Union.

McEvilly, T.V., 1964. Central US crust–upper mantle structure from Love and Rayleigh wave phase velocity inversion, *Bull. seism. Soc. Am.*, **54**, 1997–2015.

McNutt, M.K. & Judge, A.V., 1990. The superswell and mantle dynamics beneath the south pacific, *Science*, **248**(4958), 969–975.

Mégnin, C. & Romanowicz, B., 2000. The three-dimensional shear velocity structure of the mantle from the inversion of body, surface and higher-mode waveforms, *Geophys. J Int.*, **143**, 709–728.

Mitra, S., Priestley, K., Gaur, V. & Rai, S., 2006. Shear-wave structure of the south Indian lithosphere from Rayleigh wave phase-velocity measurements, *Bull. seism. Soc. Am.*, **96**, 1551–1559.

Montagner, J.-P., 1986. Regional three-dimensional structures using long-period surface waves, *Ann. Geophys.*, **4**, 283–294.

Montagner, J.-P., 2007. Upper mantle structure: global isotropic and anisotropic tomography, in *Treatise on Geophysics*, Volume 1: *Seismology and Structure of the Earth*, pp. 559–590, eds Dziewonski, A.M. & Romanowicz, B., Elsevier.

Montagner, J.-P. & Nataf, H.-C., 1986. Simple method for inverting the azimuthal anisotropy of surface waves, *J. geophys. Res.*, **91**, 511–520.

Montagner, J.-P. & Tantimoto, T., 1991. Global upper mantle tomography of seismic velocity and anisotropies, *J. geophys. Res.*, **96**, 2037–2035.

Moulik, P. & Ekström, G., 2014. An anisotropic shear velocity model of the Earth’s mantle using normal modes, body waves, surface waves and long-period waveforms, *Geophys. J Int.*, **199**(3), 1713–1738.

Nettles, M. & Dziewoński, A.M., 2008. Radially anisotropic shear velocity structure of the upper mantle globally and beneath North America, *J. geophys. Res.*, **113**, 2303, doi:10.1029/2006JB004819.

Niu, F., Solomon, S.C., Silver, P.G., Suetsugu, D. & Inoue, H., 2002. Mantle transition-zone structure beneath the south pacific superswell and evidence for a mantle plume underlying the society hotspot, *Earth planet. Sci. Lett.*, **198**(3), 371–380.

Panning, M. & Romanowicz, B., 2006. A three-dimensional radially anisotropic model of shear velocity in the whole mantle, *Geophys. J Int.*, **361–379**.

Panning, M., Lekic, V. & Romanowicz, B., 2010. Importance of crustal corrections in the development of a new global model of radial anisotropy, *J. geophys. Res.*, **115**, doi:10.1029/2010JB007520.

Parsons, B. & Sclater, J., 1977. An analysis of the variation of ocean floor bathymetry and heat flow with age, *J. geophys. Res.*, **82**, 803–827.

Pasyanos, M.E., Masters, T.G., Laske, G. & Ma, Z., 2014. LITHO1.0: an updated crust and lithospheric model of the Earth, *J. geophys. Res.*, **119**, 2153–2173.

Priestley, K. & McKenzie, D., 2006. The thermal structure of the lithosphere from shear wave velocities, *Earth Planet. Sci. Lett.*, **244**, 285–301.

Priestley, K. & McKenzie, D., 2013. The relationship between shear wave velocity, temperature, attenuation and viscosity in the shallow part of the mantle, *Earth Planet. Sci. Lett.*, **381**, 78–91.

Priestley, K., Jackson, J. & McKenzie, D., 2008. Lithospheric structure and deep earthquakes beneath India, the Himalaya and southern Tibet, *Geophys. J Int.*, **172**, 345–362.

Ritsema, J., Deuss, A., van Heijst, H. & Woodhouse, J., 2011. S40RTS: a degree-40 shear-velocity model for the mantle from new Rayleigh wave dispersion, teleseismic traveltime and normal-mode splitting function measurements, *Geophys. J Int.*, **184**, 1223–1236.

Ritzwoller, M.H., Shapiro, N.M. & Zhong, S., 2004. Cooling history of the pacific lithosphere, *Earth Planet. Sci. Lett.*, **226**, 69–84.

Saito, M., 1988. Disper80: A subroutine package for the calculation of seismic normal-mode solutions, in *Seismological Algorithms*, pp. 293–319, ed. Doornbos, D.J., Academic Press.

Schaeffer, A.J. & Lebedev, S., 2013. Global shear speed structure of the upper mantle and transition zone, *Geophys. J Int.*, **194**, 417–449.

Shapiro, N. & Ritzwoller, M., 2002. Monte-carlo inversion for a global shear-velocity model of the crust and upper mantle, *Geophys. J Int.*, **151**, 88–1052.

Smith, M.L. & Dahlen, F.A., 1973. The azimuthal dependence of Love and Rayleigh wave propagation in a slightly anisotropic medium, *J. geophys. Res.*, **78**, 3321–3333.

Sutherland, R., Spasojevic, S. & Gurnis, M., 2010. Mantle upwelling after Gondwana subduction death explains anomalous topography and subsidence histories of eastern new Zealand and west antarctica, *Geology*, **38**(2), 155–158.
Takeuchi, H. & Saito, M., 1972. Seismic surface waves, in *Methods in Computational Physics*, 1st edn, pp. 217–295, ed. Bolt, B., Academic.

Trampert, J. & Woodhouse, J., 1995. Global phase velocity of Love and Rayleigh waves between 40 and 150 seconds, *Geophys. J. Int.*, **22**, 675–690.

Trampert, J. & Woodhouse, J., 2000. Assessment of global phase velocity models, *Geophys. J. Int.*, **144**, 165–174.

van Heijst, H. & Woodhouse, J., 1997. Measuring surface-wave overtone phase velocities using a mode-branch stripping technique, *Geophys. J. Int.*, **131**, 209–230.

van Heijst, H. & Woodhouse, J., 1999. Global high-resolution phase velocity distributions of overtone and fundamental-mode surface waves determined by mode branch stripping, *Geophys. J. Int.*, **137**, 601–620.

Visser, K., Lebedev, S., Trampert, J. & Kennett, B., 2007. Global Love wave overtone measurements, *Geophys. Res. Lett.*, **34**, doi:10.1029/2006GL028671.

Visser, K., Trampert, J. & Kennett, B., 2008. Global anisotropic phase velocity maps for higher mode Love and Rayleigh waves, *Geophys. J. Int.*, **172**, 1016–1032.

Wessel, P. & Smith, W., 1995. New version of the Generic Mapping Tool release, *EOS, Trans. Am. geophys. Un.*, **76**, 329.

Zhang, Y.-S. & Tanimoto, T., 1991. Global Love wave phase velocity variation and its significance to plate tectonics, *Phys. Earth planet. Inter.*, **66**, 160–202.

Zhou, Y., Dahlen, F., Nolet, G. & Laske, G., 2005. Finite-frequency effects in global surface-wave tomography, *Geophys. J. Int.*, **163**, 1087–1111.

## APPENDIX A: THE EFFECT OF CORRELATION LENGTH

The choice of correlation length $L_c$ depends on the overall path coverage of the model. In order to regularize regions where path coverage is poorer, a larger value of $L_c$ is necessary. However, large values of $L_c$ also tend to smooth out small-scale variations in structure. We have chosen to use an increasing correlation length with increasing depth, reflecting the decrease in path coverage with depth. We use $L_c = 400$ km at the shallower depths ($\le 150$ km) to avoid smoothing out smaller-scale structures at these depths, but $L_c = 600$ km at intermediate depths (200–250 km) and $L_c = 800$ km at deeper depths ($>350$ km). Fig. A1 shows the tomography models at three representative depths determined with different values of $L_c$.

## APPENDIX B: MODEL RESOLUTION TEST

The resolution of model CAM2016SH depends on the path density, frequency content, modal makeup of the measurements, regularization, parametrization and on the theory used in the modelling. Our waveform inversion approach assumes propagation along the great circle, without mode coupling. This is commonly assumed for long period surface waves and the uses of this simple theory provide the advantage of allowing the automated treatment of millions of seismograms, which would not be possible with more sophisticated theories. Because we assume great-circle propagation we do not include major arc or multi-orbit data. The waveform inversion provides a path average shear velocity model and a *posteriori* error for each seismogram and at each depth which allow us to determine if the inverted shear velocity is well resolved. Paths which poorly constrain the wave speed structure at a given depth are not included

in the tomographic inversion building CAM2016SH. Path density and azimuthal coverage for each depth are plotted in Figs B1 and B2. The number of paths gradually falls with increasing depth because sensitivity at these depths for the fundamental mode decreases and for deeper depths the model is largely constrained by the higher modes data. Coverage is good in the northern hemisphere, especially across Asia and North America. The poorest coverage is for the southwest Indian Ocean and Antarctica. Although the southern hemisphere is less well covered due to a smaller number of seismographs, we still achieved a path coverage of greater than 50 paths per $4^\circ \times 4^\circ$ cell for the southern Hemisphere.

We conduct synthetic tests to evaluate the resolution of CAM2016SH for specific features observed in the model. We compute synthetic seismograms using modal summation including the fundamental and up to the fifth overtone for all propagation paths used in constraining the model. The mantle model used in computing the synthetic seismograms for the resolution tests consists of the same laterally uniform reference mantle model used in the waveform inversion with a number of cylindrical perturbations of varying radii and thickness placed at various geographic locations corresponding to features seen in CAM2016SH and other $\beta_h(z)$ tomographic models. The radii, thickness and strength of the cylinders are chosen to roughly correspond to the features they are meant to represent. We have added real ground noise to the synthetic seismograms before inverting them.

We make measurements on the synthetic seismograms for the same frequencies and modes as for the actual data used in building model CAM2016SH. We then invert these measurements using the inversion procedure as outlined in Section 2. The information constraining the synthetic model will, therefore, be essentially the same as that constraining CAM2016SH. As in the development of model CAM2016SH, the horizontal correlation length is 400 km for the tests at 100–150 km depth, 600 km correlation for 200–250 km depth, and 800 km correlation for deeper depths. Tests such as these and other resolution tests not shown here give a realistic assessment of the vertical and horizontal resolution in various parts of the model. The tests suggests that the horizontal resolution of CAM2016SH is 800–1000 km even in areas with poor path coverage. In particular, these tests and those of Section 2.4 show that we do resolve features like the base of the cratons and structure within the mantle transition zone.

## APPENDIX C: MODEL COMPARISON

In Section 5, we compute the correlation between our model and four published $\beta_h(z)$ models for the upper mantle. That comparison showed relatively high correlation at shallow depths ($<200$ km) to about spherical degree 20 but that the correlation became significantly worse with increasing depth and higher spherical degree. Fig. C1 shows a similar comparison between three previously published $\beta_h(z)$ models. The correlation between these models is similar to that shown in Fig. 6; that is, relatively high correlation at shallow depths ($<200$ km) to about spherical degree 20 but poorer correlation with increasing depth and spherical degree. This lack of correlation between the various $\beta_h(z)$ models for the upper mantle suggests that the $\beta_h(z)$ structure between the shallow upper mantle and the transition zone is yet to be well constrained with confidence.
<figure><img src="image_9.png" /><figcaption>**Figure A1.** The tomography model at three depths, 100, 250 and 450 km, using $L_c = 400$, 600 and 800 km. The features observed in the 100 km deep, $L_c = 800$ km map are well defined and visually correlate with tectonic features—ridges, backarc basins, cratons, etc. These features become better defined and remain stable when $L_c$ is decreased to 600 and 400 km. This is not the case for the deeper maps. At 250 km depth, the features in the $L_c = 800$ km map remain stable in the $L_c = 600$ km map but begin to break up in the $L_c = 400$ km map. At 450 km depth, features look stable in the $L_c = 800$ km map but begin to break up in the $L_c = 600$ km depth map. In examining these and similar maps at intervening depths, we have chosen to use $L_c = 400$ for 100 and 150 km depth, $L_c = 600$ km for the 200 and 250 km maps, and $L_c = 800$ km for the maps at deeper depths.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Rays through a 4x4 cell</figcaption></figure>

**Figure B1.** Density for the paths per 4° × 4° cell constraining the tomographic model CAM2016SH at each depth. Coverage is good in the Northern Hemisphere; the Southern Hemisphere coverage is less good owing to the smaller number of seismographs but is still better than ~20 paths in the vicinity of the South Pole.
<figure><img src="image_11.png" /><figcaption>**Figure B2.** Voronoi maps showing the azimuthal coverage of the paths constraining the tomographic model CAM2016SH at each depth. The Voronoi cells in this figure display the smallest spatial region for which there is sufficient azimuthal coverage to resolve the 4θ term for Love waves. The azimuthal coverage is excellent for the Northern Hemisphere and much of the Southern Hemisphere. The azimuthal coverage is less good for Antarctica and poor for the southwestern Indian Ocean.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Figure B3.** Resolution tests for the East Pacific Rise. The spatial extent of the input cylinders is denoted by the red circles in the 100 km depth map and the red rectangles in the cross-sections.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure B4.** Resolution tests for the Indian Ocean Ridge. The format for this figure is the same as for Fig. B3.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Figure B5.** Resolution test for the African region. This test consists of three high-wave-speed cylinders representing the West African, Congo and Kalahari cratons and two low-wave-speed cylinders representing the Darfur and Afar hot spots. The format for this figure is the same as for Fig. B3.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Figure B6. Resolution tests for the Australian region. The format for this figure is the same as for Fig. B3.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure B7. Resolution tests for Asia. The format for this figure is the same as for Fig. B3.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>**Figure B8.** Resolution tests for the North America region. The format for this figure is the same as for Fig. B3.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>**Figure C1.** Correlations computed between S362ANI and SAW642ANb (black), between S362ANI and SEMum2 (red), between S362ANI and FFSW1 (green), between SAW642ANb and SEMum2 (blue), between SAW642ANb and FFSW1 (purple) and between SEMum2 and FFSW1 (grey). Correlations with S362ANI are truncated up to degree 18 and correlations with FFSW1 are truncated up to degree 20. Significance levels of 95 per cent (dotted) and 66 per cent (dashed) are plotted.</figcaption></figure>