# Seismic evidence for broad attenuation anomalies in the asthenosphere beneath the Pacific Ocean

Alice Adenis, Eric Debayle, Yanick Ricard

► To cite this version:

Alice Adenis, Eric Debayle, Yanick Ricard. Seismic evidence for broad attenuation anomalies in the asthenosphere beneath the Pacific Ocean. Geophysical Journal International, 2017, 209 (3), pp.1677–1698. 10.1093/gji/ggx117. hal-02328803

HAL Id: hal-02328803

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

Submitted on 26 Aug 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>

# Seismic evidence for broad attenuation anomalies in the asthenosphere beneath the Pacific Ocean

Alice Adenis, Eric Debayle and Yanick Ricard

Laboratoire de Géologie de Lyon: Terre, Planètes et Environnement, CNRS, Université de Lyon 1, Ecole Normale Supérieure de Lyon, Villeurbanne, France.

E-mail: aliceadenis@gmail.com

Accepted 2017 March 20. Received 2017 January 31; in original form 2016 October 18

## SUMMARY

We present QADR17, a global model of Rayleigh-wave attenuation based on a massive surface wave data set (372 629 frequency-dependent attenuation curves in the period range 50–260 s). We correct for focusing-defocusing effects and geometrical spreading, and perform a stringent selection to only keep robust observations. Then, data with close epicentres recorded at the same station are clustered, as they sample the same Earth’s structure. After this pre-selection, our data set consists of about 35 000 curves that constrain the Rayleigh-wave intrinsic attenuation in the upper mantle. The logarithms of the attenuation along the individual rays are then inverted to obtain global maps of the logarithm of the local attenuation. After a first inversion, outliers are rejected and a second inversion yields a variance reduction of about 45 per cent. Our attenuation maps present strong agreement with surface tectonics at periods lower than 200 s, with low attenuation under continents and high attenuation under oceans. Over oceans, attenuation decreases with increasing crustal ages, but at periods sensitive to the uppermost 150 km, mid-ocean ridges are not characterized by a very localized anomaly, in contrast to what is commonly observed for seismic velocity models. Attenuation is rather well correlated with hotspots, especially in the Pacific ocean, where a strong attenuating anomaly is observed in the long wavelength component of our signal at periods sampling the oceanic asthenosphere. We suggest that this anomaly results from the horizontal spreading of several thermal plumes within the asthenosphere. Strong velocity reductions associated with high attenuation anomalies of moderate amplitudes beneath the East Pacific Rise, the Red Sea and the eastern part of Asia may require additional mechanisms, such as partial melting.

**Key words:** Composition of the mantle; Surface waves and free oscillations; Seismic attenuation; Seismic tomography; Pacific Ocean.

## 1 INTRODUCTION

Most global tomographic studies of the upper mantle focus on the 3-D distribution of seismic velocities. These studies are generally based on surface waves, which provide global coverage of the upper mantle and have a strong sensitivity to the shear wave (S-wave) velocity. Over the last decade, there have been important improvements in mapping 3-D shear-velocity heterogeneities of the upper mantle (e.g. Ritsema et al. 2011; Debayle et al. 2016). Recent global S-wave tomographic models are obtained from the inversion of massive data sets and show robust patterns, especially in the top of the upper mantle where they agree for horizontal wavelengths smaller than 1000 km (Meschede & Romanowicz 2015; Debayle et al. 2016). However, the Earth’s mantle does not behave as a perfectly elastic body and it is in principle possible to extract information on its anelastic properties from the decay of the amplitude of seismic waves.

Such attenuation studies are difficult, because the wave amplitude is influenced by a number of mechanisms. It is affected by uncertainties in the source excitation (Um & Dahlen 1992) including the scalar seismic moment $M_0$, by the geometrical spreading of the wave front, by propagation effects such as focusing and defocusing (Lay & Kanamori 1985; Woodhouse & Wong 1986), by short wavelength scattering (e.g. Ricard et al. 2014), by local site response and by the calibration of the measuring devices (Dalton et al. 2014). Amplitude is also affected by various intrinsic anelastic mechanisms converting elastic energy into heat, mechanisms such as interaction of the waves with phase changes (Durand et al. 2012), with crystal dislocations or with partial melting among other mechanisms (see Jackson 2007, for a review). Attenuation studies aim to correct the measured amplitudes from all the elastic effects, in order to obtain the intrinsic attenuation. Because of the difficulties of measurements and of the numerous mechanisms that need to be accounted for, less work has been done on attenuation than on
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Up: average $\ln(Q)$ curves of the entire data set, with the root-mean-square deviation represented as error bar. Dashed curve corresponds to the average of the 372 629 $\ln(Q_i)$ measurements before pre-processing while the continuous line is the average of the selected measurements used for the construction of the final model (between 25 971 and 33 180 paths). Bottom: radial sensitivity kernel of $\ln(Q)$ with respect to shear attenuation for Rayleigh-wave fundamental model.</figcaption></figure>

velocity, and the agreement between recent global S-wave attenuation models is limited to very long wavelengths, greater than ~5000 km (Selby & Woodhouse 2002; Gung & Romanowicz 2004; Dalton et al. 2008).

Attenuation and velocity have different sensitivities to the physical properties of the Earth such as temperature, composition, and the presence of partial melt or water content. Improving our knowledge of attenuation would bring information complementary to velocity about Earth's interior. Furthermore, seismic velocities are affected by dispersion effects caused by attenuation, and a precise knowledge of the attenuation would improve the resolution of the velocity structure of the Earth (Romanowicz 1990; Karato 1993).

Most surface wave attenuation studies are based on Rayleigh-wave amplitudes. As shown in Fig. 1, the fundamental mode of Rayleigh wave for periods up to 250 s provides sensitivity to the whole upper mantle. This sensitivity can be improved by including overtones in the modelling (Gung & Romanowicz 2004). In previous surface-wave studies, the attenuation is generally found to be coherent with large-scale surface tectonics down to about 200–250 km depth, with low attenuation beneath continents and high attenuation under oceans. Some studies (Billien et al. 2000; Dalton & Ekström 2006; Dalton et al. 2008; Ma et al. 2016) suggest a good correlation with the velocity models. Ma et al. (2016), Dalton et al. (2008) and Warren & Shearer (2002) find a strong dependence on seafloor age whereas Gung & Romanowicz (2004) and Romanowicz (1995) observe high-attenuation anomalies under the southern Pacific and Africa, correlated with the hotspot distribution. Below 200–250 km, a change in the pattern is observed: high-attenuation regions are found under the southeastern Pacific and beneath eastern Africa and low-attenuation regions seem to be associated with subduction zones.

In this study we map upper-mantle Rayleigh-wave attenuation at global scale using a data set of 372 629 attenuation curves measured by Debayle & Ricard (2012). This data set, which has never been used in an attenuation study, provides global coverage with a large redundancy. We have made a major effort to reject measurements that are likely to be affected by mechanisms that are not accounted for in our modelling. We also correct measurements for focusing-defocusing and exploit the large redundancy of our data set to minimize errors. Our final Rayleigh-wave attenuation model presents a strong correlation with surface tectonics, (high-attenuating oceans and low-attenuating continents), and a decrease of attenuation with the age of the ocean floor. We show that mid-ocean ridges are less prominent than in velocity models and that attenuating regions are located around hotspot in oceanic regions. A high-attenuation anomaly is found in the middle of the Pacific in the period range 50–100 s, this anomaly is not observed in velocity models. We suggest that attenuation and velocity models are compatible with a thermal anomaly in the central Pacific resulting from several plumes, and with the presence of partial melt beneath the East Pacific Rise, the Red Sea and the eastern part of Asia.

## 2 DATA AND METHOD

The attenuation model is built using a two-step procedure. The first step provides a set of path-average attenuation curves and is summarized in Section 2.1. An important pre-processing is made in order to remove attenuating effects other than the intrinsic attenuation from our measurements. This pre-processing is described in Section 2.2. The last step, presented in Section 2.3, involves the regionalization of the attenuation curves by a modified version of the continuous regionalization approach used previously for velocities (Montagner 1986; Ricard et al. 1996; Debayle & Sambridge 2004).
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Attenuation maps at 100 s (top) and 200 s (bottom) periods obtained after the first inversion (left column); after rejecting data for which the misfit increases (final inversion, middle column); from the inversion of the rejected data (right column). The global correlation between the first and final inversion is above 0.91. The data set used to obtain maps (a) contains both (b) and (c) data sets. Maps built with the rejected data (right column) do not display a pattern coherent with large scale tectonics. Hotspot locations, according to Müller *et al.* (1993), are indicated with blue circles. Plate contours are in green.</figcaption></figure>

## 2.1 Data

Our data set consists of 372 629 attenuation curves for the fundamental mode of Rayleigh waves measured by Debayle & Ricard (2012) from the automation of a waveform inversion approach initially developed by Cara & Lévêque (1987). We summarize here their approach for extracting an attenuation curve from a Rayleigh-wave seismogram and refer to sections 2, 3 and 4.1 of Debayle & Ricard (2012) for a detailed description. This waveform inversion is applied to every single seismogram.

First, a synthetic seismogram is calculated for the epicentre-station path corresponding to the recorded waveform. This synthetic is computed in a 1-D model which includes a path-averaged crust structure estimated from 3SMAC (Nataf & Ricard 1996) and a radially anisotropic mantle, close to PREM (Dziewonski & Anderson 1981). Debayle & Ricard (2012) noted that the attenuating layer located between 80 km and 200 km depths in PREM is not adapted to continental paths for which the attenuating layer is less pronounced or located deeper. For this reason, they replace as a starting model, the quality factor $Q_μ(z)$ of PREM in the whole upper mantle by the uniform value of 200. Fig. 2 of Debayle & Ricard (2012) displays their starting model. Phase velocities are corrected from physical dispersion using Kanamori & Anderson (1977) assuming a reference period of 100 s.

Second, the difference between the synthetic seismogram and the actual waveform is inverted for a set of fundamental and higher modes dispersion $c(T)$ and attenuation curves $Q(T)$, where $T$ is the period. In order to account for period-independent amplitude differences between the synthetic and recorded waveforms, Debayle & Ricard (2012) also invert for the scalar seismic moment $M_0$.

We first describe our processing of this data set, prior to build attenuation maps and then discuss in details their geodynamical implications. We focus on the fundamental mode which is easier to interpret and leave the higher modes and the inversion with depth to a further study.

## 2.2 Pre-processing

The amplitude of a fundamental mode Rayleigh wave measured at a seismic station can be expressed as

$$A = A_S \cdot A_I \cdot A_G \cdot A_{\text{Foc}} \cdot A_{\text{Diff}} \cdot A_{\text{Int}} \quad (1)$$

where $A_S$ is due to excitation at the source, $A_I$ is due to the instrumental response, $A_G$ is the geometrical spreading factor, $A_{\text{Foc}}$ is produced by the focusing-defocusing of seismic waves due to 3-D heterogeneities, $A_{\text{Diff}}$ is the effect of scattering and diffraction on amplitudes and $A_{\text{Int}}$ describes the amplitude decay due to the intrinsic attenuation of the Earth. These amplitude terms are generally functions of the period of the wave in addition to various other quantities (epicentre-station distance, source radiation pattern, quality of the seismometer...). As we are interested in the effect of the intrinsic attenuation $A_{\text{Int}}$, we need to make sure that our measurements are not biased by a poor knowledge of $A_S$, $A_I$ and $A_{\text{Diff}}$ ($A_G$ is well known), and that we properly account for the important propagation effect $A_{\text{Foc}}$. This requires some pre-processing before the tomographic inversion.

The computations of the source, the instrument and the geometric effects, $A_S$, $A_I$ and $A_G$ are included in the synthetic calculation of the automated waveform inversion (see eq. 1 of Debayle & Ricard 2012). Possible errors in $A_I$ are also accounted for in the waveform inversion. First, Debayle & Ricard (2012) rejected waveforms for which the amplitude of the synthetic before inversion differs by a factor greater than 10 from the amplitude of the actual data. In this study, we use a more drastic criterion and we keep only data for which the amplitudes of the synthetic and actual waveforms differ by a factor less than 2. This reduces the 372 629 paths of Debayle & Ricard (2012) to 359 627 paths (Table 1). Second, Debayle & Ricard (2012) invert for the scalar seismic moment $M_0$, parametrized by $log_{10}(M_0)$ and use a large *a priori* standard deviation ($\sigma_{log_{10}(M_0)} = 0.5$). If required by the data, this absorbs in the inverted $M_0$ any difference between the synthetic and the observation which is frequency independent.
<table><caption>**Table 1.** Number of data at each step of the selection process, from the pre-processing step up to the construction of the final model.</caption><thead><tr><th>Processing step</th><th>Number of data</th></tr></thead><tbody><tr><td>Initial data set from Debayle & Ricard (2012)</td><td>372 629</td></tr><tr><td>Ratio between synthetic and real seismogram ≤ 2</td><td>359 627</td></tr><tr><td>Focusing effect - Path length ≤ 110°</td><td>332 354</td></tr><tr><td>Source sensitivity ≤ 20 per cent</td><td>186 024–194 113</td></tr><tr><td>Clustering step</td><td>65 277–66 605</td></tr><tr><td>Residual based selection</td><td>25 971–33 180</td></tr></tbody></table>

$A_{Diff}$ represents the effects of scattering and diffraction which are difficult to estimate. Ricard *et al.* (2014) indicate that, for a realistic power-law distribution of heterogeneities, scattering attenuation should be frequency independent. At long period (≥50 s), fundamental Rayleigh waveforms are in general simple and various studies suggest that the scattering effects are negligible (Ferreira & Woodhouse 2007; Dalton *et al.* 2014; Meschede & Romanowicz 2015). For these reasons, we consider that $A_{Diff} \sim 1$, but our final attenuation maps may include a small contribution from the local density of scatters.

Our modelling approach (described in Section 2.3) uses the great-circle ray approximation (GCRA) which considers that surface waves propagate along the source-station great-circle and that they are only sensitive to the structure along a zero-width ray. This theory has proved to be of great efficiency in mapping velocities (e.g. Ritsema *et al.* 2011; Debayle *et al.* 2016) and it has recently been shown that GCRA accurately predicts phases of long period Rayleigh waves in the period range 50–160 s for distances shorter than 110° (Parisi & Ferreira 2016). However, GCRA does not account for focusing-defocusing effects which are known to affect the amplitude of surface waves (Lay & Kanamori 1985). We improve our amplitude estimate using GCRA by computing $A_{Foc}$ along each ray using the formalism of Woodhouse & Wong (1986) and the phase velocity maps of Durand *et al.* (2015). According to Woodhouse & Wong (1986),

$$ \ln(A_{Foc}) = \frac{1}{2 \sin(\Delta)} \int_{0}^{\Delta} \sin(\Delta - \phi) \left[ \sin \phi \frac{\partial^2}{\partial \theta^2} - \cos \phi \frac{\partial}{\partial \phi} \right] \frac{\delta c}{c} d\phi, \quad (2) $$

where $\Delta$ is the epicentre-station distance, $\phi$ is the along-path coordinate away from the epicentre, $\theta$ is the path-perpendicular coordinate, $\delta c/c$ is the relative perturbation in phase velocity at a given period. $A_{Foc}$ depends primarily on the second derivative of the relative perturbation in phase velocity perpendicular to the path. The phase velocity maps have been obtained from the regionalization of Debayle & Ricard’s (2012) data set extended with measurements at longer periods (100–350 s). The Woodhouse & Wong’s (1986) method gives a fair approximation of focusing effects for path lengths between 70° and 110°. For paths shorter than 70°, focusing is slightly underestimated, while it can be strongly overestimated above 110° (see fig. 11(a) in Dalton *et al.* 2014). As the errors remain small for paths shorter than 70° (Dalton *et al.* 2014), we retain these data in order to preserve an optimal coverage, especially in oceanic regions. This allows us to keep 332 354 paths for next step (Table 1). We checked that removing the paths shorter than 70° does not affect the conclusions of this study. Examples of attenuation maps obtained with and without focusing correction for our final data set are displayed in Fig. C1. Fig. C2 illustrates the effect of our focusing corrections beneath the East Pacific Rise and North America. These figures show a good agreement with similar

estimates made by previous authors (compare with e.g. Dalton & Ekström 2006, figs 11 and 12).

The source excitation $A_S$ is computed following Cara (1979) using the global Centroid-Moment Tensors (CMT) solution (Dziewonski *et al.* 1981; Ekström *et al.* 2012). To improve accuracy, we use a specific 1-D model for the source by extracting density, seismic velocities and attenuation from 3SMAC (Nataf & Ricard 1996) beneath the epicentre. It is well known that deviation in take-off azimuth of the ray at the source can produce an amplitude anomaly of the same order than produced by focusing-defocusing (Um & Dahlen 1992). We therefore eliminated all paths close to a nodal plane of the source radiation pattern (we remove all paths for which the computed amplitudes vary by more than 1 per cent for a take-off azimuth variation of 1°). Fig. C3 shows that removing these paths has a moderate effect on the inverted model. It is likely that the requirement of an amplitude factor <2 between synthetic and actual waveforms has already allowed us to reject most data close to a nodal plane in the source radiation pattern. However, this second criteria guarantees that all paths for which a small error in the take-off azimuth or in the nodal plane orientation produces a large amplitude error are rejected. This selection is applied for every frequency, and the number of remaining data after this selection varies between 186 024 at 260 s of period and 194 113 at 90 s of period (Table 1).

## 2.3 Tomographic inversion

Assuming that the pre-processing described in Section 2.2 has allowed us to reject all paths for which $A_S$ and $A_I$ are likely to suffer from significant errors, that we know $A_G$, $A_I$ and $A_S$, that we correctly account for $A_{Foc}$ and can reasonably neglect $A_{Diff}$, we can extract the quality factor $Q_{Int}(T)$ for each ray $i$, hereafter called $Q_i$, related to $A_{Int}$ by

$$ \ln[A_{Int}(T, \Delta)] = - \frac{\pi \Delta}{T c(T) Q_{Int}(T)}. \quad (3) $$

Note that we use the phase velocity in the definition of the quality factor in agreement with what Aki & Richards (2002) called spatial quality factor. Another definition, using the group velocity is possible (the temporal quality factor of Aki & Richards 2002). As the quality factor varies laterally by one to two orders of magnitude while the phase and group velocities only differ by a few 10 per cent, the difference in the definitions is not crucial.

We then combine the $Q_i$ measurements in a tomographic inversion using a continuous regionalization scheme (Montagner 1986; Ricard *et al.* 1996; Debayle & Sambridge 2004): we invert the path-average $\ln(Q_i)$ for the local quality factor $\ln(Q(r))$ at position $r$. The use of logarithms brings our data set closer to a Gaussian distribution, whose average and root-mean-square deviation are presented for each period in Fig. 1. It also guarantees to avoid negative values for the attenuation in the inverted model.

For a source-receiver path $i$, the forward problem can be written at a given period $T$:

$$ \ln(Q_i^{-1}) = \ln \left( \frac{1}{\Delta_i} \int_i \exp[\ln(Q^{-1}(s))] ds \right), \quad (4) $$

where $\Delta_i$ is the epicentre-receiver distance, $Q_i$ is the path-average quality factor for path $i$, $s$ is the abscissa along the path and $Q(s)$ is the local quality factor at abscissa $s$ (i.e. at position $r(s)$). This problem can be considered as a non-linear relationship $d = g(m)$ between a data vector $d$ that contains $\ln(Q_i^{-1})$ along each path $i$ and a parameter vector $m$ that contains $\ln(Q^{-1}(r))$. The inversion is
performed using the iterative least-square solution for continuous non-linear inverse problem proposed by Tarantola & Valette (1982). The developments of the equations are detailed in Appendix A.

In order to obtain a smooth model, we assume a Gaussian *a priori* covariance function for the model:

$$C_{m_0}(\mathbf{r}, \mathbf{r}') = \sigma_m^2 \exp \left( -\frac{\Delta^2(\mathbf{r}, \mathbf{r}')}{2L_{\text{corr}}^2} \right) \quad (5)$$

where $\mathbf{r}$ and $\mathbf{r}'$ are the coordinates of two geographical points separated by an angular distance $\Delta(\mathbf{r}, \mathbf{r}')$, the angular correlation length $L_{\text{corr}}$ controls the lateral degree of smoothing of the final model, and the *a priori* error $\sigma_m$ the allowed variation of amplitude at every point of the map. At each period $T$, the initial model $m_0$ is the average of the observed $\ln(Q_i)$ and the *a priori* error on the model $\sigma_m$ is set to 20 per cent of $m_0$.

After several tests, we choose a correlation length $L_{\text{corr}} = 10^\circ$. This choice is conservative and we favour a smooth model. We tried other correlation lengths (see Fig. C5) that show that the variance reduction does not significantly improve for shorter correlation lengths. Our choice is also justified by the fact that focusing effects are corrected with phase velocity maps built using $L_{\text{corr}} = 3.6^\circ$ (Durand et al. 2015). As the focusing correction depends on the derivatives of the phase velocity anomaly, we expect that corrections performed with such phase velocity maps are valid for smoother $Q$ models. A synthetic experiment (see Section 3.2) suggests that with a $10^\circ$ correlation length, we resolve the first 12 spherical harmonic degrees very accurately and from degrees 12 to ~16 with a damped amplitude.

The *a priori* covariance on the data $C_{d_0}$ is diagonal, and its diagonal terms are $\sigma_{d_i}^2$ where $\sigma_{d_i}$ is the *a priori* data error that controls the contribution of path $i$ in the final model. We chose $\sigma_{d_i}$ constant, and equal to the root-mean-square deviation of the data set at each period. To avoid redundant data, we cluster paths corresponding to events with close epicentres recorded at a same station, assuming the waves propagate through the same attenuation structure. We use a cluster radius of $2^\circ$, small compared the horizontal degree of smoothing used in the inversion. We reject outliers in each cluster (see Appendix B) and compute an average $\ln(Q_i)$ per cluster at each period $T$. Fig. C4 compares maps with (right column) and without (left column) clustering at 70, 140 and 200 s. Results are very close for all degrees, with correlation above the 95 per cent significance level. This confirms that grouping nearby paths has no effect on the long wavelength pattern discussed in this study. Without grouping, each inversion requires more than one week of computation per period compared to one day with grouping. We therefore cluster nearby paths, which gives us more flexibility to test the model. After the clustering step, the number of paths reduces to values ranging from 65 227 to 66 605, depending on the period (Table 1).

A first inversion is performed using eqs (A7), (A8) and (A9). The residual between a data $i$ and the model at iteration $k$ is defined as $r_{i,k} = |d_i - g_i(\hat{m}_k)|$ and the variance reduction for the $n$ observations as:

$$V_R^k = 1 - \frac{\sum_{i=1}^n r_{i,k}^2}{\sum_{i=1}^n r_{i,0}^2} \quad (6)$$

where $r_{i,0}$ and $r_{i,k}$ are initial ($k=0$) and final ($k$ = last iteration) residuals for each data $i$. In all the inversion that we have performed, we observe very little evolution of the model after the first iteration. For this reason, all inversions presented in this paper are stopped after the first iteration.

After the first inversion, the variance reduction reaches hardly 10 per cent suggesting that significant noise remains in our data set.

We compare the residual $r_i$ before and after inversion for every path $i$. Data for which $r_i$ increases are better explained by our starting 1-D model than the inverted model. These data are rejected; they are likely associated with cases where the ray path inversion did not resolve the attenuation that was therefore maintained close to the *a priori* uniform guess. Typically, 50 per cent of the path clusters are rejected and the number of selected data after the first inversion ranges between 25 971 and 33 180 depending on period (Table 1).

The second and final inversion is then conducted with the remaining data set and using the same *a priori* information as in the first inversion. The second inversion significantly improves the variance reduction from 10 per cent up to 48 per cent. Fig. 2 displays the maps at 100 and 200 s periods obtained after the first and second inversions and after an inversion of the rejected data. The maps after the first and second inversions display very similar patterns, although contrasts are greater after the second inversion, due to the removing of inconsistent data. The inversion of the rejected data yields maps with a weaker amplitude close to the *a priori* model, that do not display a pattern coherent with large scale tectonics. The average $\ln(Q)$ of our final data set and its standard deviation are presented for each period in Fig. 1. The large standard deviations reflect the broad range of attenuation variations observed in the upper mantle.

Ray density maps corresponding to the final data set are presented at different periods on Fig. 3. We compute the ray density $D(\mathbf{r})$ at each geographical point $\mathbf{r}$ using

$$D(\mathbf{r}) = \sum_i \exp \left( -\frac{\Delta^2(\mathbf{r}, \mathbf{r}_i)}{2L_{\text{corr}}^2} \right) \quad (7)$$

where $\mathbf{r}_i$ is the point of path $i$ which is the closest to $\mathbf{r}$, and $L_{\text{corr}}$ is set to 10 degrees as in eq. (5). Therefore, at each geographical point $\mathbf{r}$, the contribution of rays is weighted by a Gaussian having the same width as that controlling the horizontal smoothing in our model. $D(\mathbf{r})$ can therefore be seen as a measure of the number of rays in a surface of radius $L_{\text{corr}}$. Fig. 3 shows that density maps are typical of global tomography, with the highest coverage beneath the continents of the northern hemisphere. However, even in oceanic regions of the southern hemisphere where coverage is the poorest, there are hundreds of paths within a distance $L_{\text{corr}}$ at each geographical points $D(\mathbf{r})$, which confirms that our data set provides global coverage with large redundancy.

# 3 RESULTS

## 3.1 QADR17 model

Fig. 4 presents our Rayleigh-wave attenuation model QADR17 at different periods between 50 and 240 s. High $Q$ values (low attenuation) are in blue, low $Q$ values (high attenuation) in red. The sensitivity kernels are shown in Fig. 1. In the period range 50–240 s, our data set provides sensitivity to the upper mantle, with maxima varying between ~70 km at 50 s and ~350 km at 240 s.

QADR17 shows a strong correlation with surface tectonics at periods lower than 140 s, mostly sensitive to the uppermost 200 km. The dominant signal is the difference between continents and oceans, associated with anomalously low and high attenuation, respectively. Beneath continents, weak attenuation is found under old continental roots, such as the African, North American, Amazonian, Siberian and Australian cratons, the Russian platform and Antarctica. A stronger attenuation is observed under Phanerozoic
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Ray density maps (≈ number of rays within a distance of 10°) obtained for our final data set at 50, 70, 100, 140, 200 and 240 s period.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Our model QADR17 of lateral variations in Rayleigh-wave attenuation at different periods. We plot values of *Q* with a logarithmic scale with the geometrical average of *Q* above the colour scale.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 5.** (a) Map of the tectonic provinces according to their age. (b) Average *Q* of our model QADR17 on every tectonic provinces as a function of the period.</figcaption></figure>

lithosphere and in active tectonic areas such as western North America, the Afar province, the Bénoué rift and the Andes. Mid-ocean ridges are associated with higher than average attenuation but interestingly, their signature does not dominate the attenuation pattern of oceanic regions. This contrasts with velocity models where the ridge signature is very strong at short periods (see Fig. C6). The location of hotspots is in general well correlated with attenuating regions and some of them are associated with a very strong signal. This is particularly clear for the Hawaii, Marquesas, Society and Galapagos hotspots in the Pacific, Meteor, Tristan da Cunha-Gough and Cap Verde hotspots in the Atlantic, and Reunion, Kerguelen and Marion hotspots in the southern Indian ocean. A strong attenuation is observed between periods of 50 and 140 s in the Wharton basin, which is bounded to the west and south by the Ninetyeast and Broken Ridges. At periods greater than 140 s, both the amplitude of attenuation variations and the correlation with surface tectonics decrease. However, broad attenuating regions remain beneath most Pacific, north Africa, south India and Atlantic ocean hotspots up to 240 s.

We display in Fig. 5 the *Q* values as a function of period, averaged for various tectonic provinces extracted from the *a priori* model 3SMAC (Nataf & Ricard 1996). To first order, *Q* varies with lithospheric age (young oceans have higher attenuation than old continents). The amplitude of variations peaks at 70 s periods, corresponding to depths of 100–150 km and then decreases with the period. Under continents, *Q* varies from about 250 beneath Paleozoic regions up to about 600 under Archean regions at 70 s. Under oceans, *Q* increases from an average value of ~130 under 0-5 Ma-old lithosphere to ~250 under 135–150 Ma lithosphere at 70 s period. Fig. 5 and Fig. C7 suggest that *Q* increases with seafloor age at periods lower than 100 s, although it does not strike the eye from Fig. 4. In the following subsection, we per-

form a synthetic experiment in order to test the robustness of our observations.

## 3.2 Synthetic experiment

We perform a synthetic experiment to explore the limits of path coverage and damping, where the input model is the pattern of attenuation extracted from the 3SMAC *a priori* model of Nataf & Ricard (1996). 3SMAC provides $Q_μ(z)$ maps, that we converted in $Q(T)$ using the sensitivity kernels. After this conversion $Q(T)$ varies from 25 beneath mid-ocean ridges up to 225 beneath cratons. A number of studies argue for higher $Q(T)$ values beneath cratons. Rayleigh-wave attenuation models of Dalton & Ekström (2006) and Ma *et al.* (2016) found that the Rayleigh-wave quality factor can reach values of ~2300 at 50 s of period beneath cratons. Our own model has variations between ~30 and ~1,800 at 50 s of period. We rescaled the attenuation variations of 3SMAC by a constant factor in order to fit the extreme $Q(T)$ values of QADR17. From the input model, a synthetic data set is calculated, using the actual ray density. This data set is then inverted using the same inversion parameters as in our real inversion ($L_{corr} = 10°, σ_{di}$ is the root-mean-square deviation of the data set and $σ_m = 20$ per cent of $m_0$).

Fig. 6 shows the input model (left column) and the output of the inversion (right column). The general pattern of the input model is well retrieved although variations are smoothed out and amplitudes are damped. It is for example difficult to distinguish on the inverted maps at 100 and 140 s the difference between Archean and Proterozoic lithospheres. The attenuation signature of mid-ocean ridges is recovered at every period. Note that in 3SMAC, mantle plumes are characterized by a high-attenuation anomaly located in narrow $2° × 2°$ conduits and that our inversion is unable to recover these narrow anomalies. This implies that the high-attenuation anomalies observed in the vicinity of some hotspots (Fig. 4) are produced by much broader structures, probably due to the spreading of plumes at the bottom of the oceanic lithosphere. For completeness, we show in Fig. 7 the correlation between the input model and the inverted model as well as their spectral ratio (output over input). The spectral amplitude of the input model is very well recovered up to degree 12, while the lateral variations are damped by the inversion process at higher degrees. The correlation between output and input model is very high, especially at low spherical harmonic degrees. Fig. 7 demonstrates that the long wavelength component of the input model is very well recovered, especially at spherical harmonic degrees $≤6$.

Fig. 8 shows the average $Q(T)$ values for the tectonic provinces of Fig. 5(a). Fig. 8 confirms that the inversion retrieves the pattern of age variation present in 3SMAC, although amplitudes are reduced by the inversion process. This damping is present at all periods and is slightly stronger at shorter periods. We observe a slight decrease of the global average after inversion (black dashed line in Fig. 8). This small bias can be attributed to uneven sampling and suggests that QADR17 slightly underestimates the global average of *Q*.

Our synthetic test shows that our data coverage and inversion process allow to recover the large scale Rayleigh-wave attenuation pattern, with only a moderate damping of the lateral variations. The mid-ocean ridge signature dominates the oceanic attenuation signal in 3SMAC and this pattern is well recovered. This comforts us with the idea that, although present in our model, the mid-ocean ridge signature is not the dominant signal in oceanic regions, contrary to what is suggested by 3SMAC.
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Synthetic test: left column is the input model 3SMAC (rescaled to match the observed amplitudes) from Nataf & Ricard (1996) and right column the retrieved model. The average values of *Q* are reported above the logarithmic scales. Although some smoothing and damping are observed, note the high correlation between the input model and the inverted model.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Top: correlation between the input and output synthetic models in Fig. 6; bottom: their spectral ratio (output/input). On the top figure the dashed line indicates the 95 per cent significance level (i.e. there is 95 per cent of chance that the correlation is not due to chance).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Average *Q* for every tectonic provinces, according to our resolution test, as a function of period. Continuous curves represent average values for the input model, dashed curves show the output values after the synthetic experiment.</figcaption></figure>

# 4 DISCUSSION

## 4.1 Comparison with other models

We first compare our attenuation maps with those published by Ma *et al.* (2016) (hereafter referred to as QMA16) and Dalton & Ekström (2006) (hereafter referred to as DEA06a). It is worth noting that these models are very similar, which may be due to a fairly similar approach for extracting the signal of attenuation from the amplitude data. The good correlation of these two models for *Q*(T) is quite different from what happens with the attenuation models for *Q*(z) (e.g. Selby & Woodhouse 2002; Gung & Romanowicz 2004; Dalton *et al.* 2008), which are very different. Fig. 9 shows that all maps share the ocean-continents differences, with low attenuation beneath most cratons and a higher attenuation in younger regions. There are more contrasts in QADR17, especially in oceanic regions that display strong attenuation variations. The signature of mid-ocean ridges is in general weak at 50 s period in QMA16 and DEA06a, except for the East Pacific Rise. It is stronger beneath the Indian and Atlantic ocean in QADR17. At 100 s and 140 s a stronger signature of mid-ocean ridges is observed in QMA16 and DEA06a, especially in the Pacific ocean. In addition to the stronger contrasts,

a significant difference between QADR17 and other models is the presence of stronger attenuating regions in old oceanic basins, some of them correlated with the presence of hotspots.

In Fig. 10, we display the correlation between QADR17 and the two other models. The correlation is fair up to degree 12 at 50 and 100 s, but low at 140 s between degree 3 and 6. The difference between models is partly due to the presence of broad attenuation anomalies in the central Pacific. We believe from our synthetic tests (Fig. 7) that these broad anomalies are well resolved in our model.

We attribute the stronger contrasts of our model with respect to the two other models displayed in Fig. 9 to our choice to invert for ln(*Q*) rather than for *Q* or *Q*⁻¹. This brings the data close to a Gaussian distribution and only data with Gaussian distributions are appropriate for a least-square algorithm. ln(*Q*) is also better suited to recover the large variations of attenuation and avoids negative values of *Q*. As an example, we consider that *Q* can vary between 70 and 600 from an average value of *Q*₀ = 200, which is chosen as the starting model of the inversion. An inversion of *Q* that would allow perturbations of 65 per cent in *Q*₀ (i.e. *Q* = *Q*₀(1 ± 0.65) = 200 ± 130) might recover the low *Q* values but not the large quality factor. An inversion allowing perturbations in *Q*₀ of ±400 could easily lead to spurious negative attenuations (the same difficulties obviously occur for an inversion on terms of *Q*⁻¹). On the contrary, no such problems occur when using ln(*Q*), and changing the starting model by only 20 per cent is sufficient to retrieve the full range of variation (i.e. ln(*Q*) = ln(*Q*₀)(1 ± 0.2) with *Q*₀ = 200 allows *Q* to be in between ~70 and ~600). We have also been very careful in our selection process to remove data with likely source or station errors. In addition, our synthetic tests (Fig. 7) show that long wavelengths are very well retrieved. We are therefore confident in the long wavelength pattern of our model.

## 4.2 Correlation with the velocity

We compare our attenuation maps with phase velocity maps built using the same data set. This comparison can bring insights on temperature and composition as attenuation and velocities have different sensitivities to temperature, water content and other parameters such as melt fraction, major elements chemistry or grain size (e.g. Shito *et al.* 2006).

We use the dispersion curves of the fundamental Rayleigh waves which have been selected for QADR17, and build a dispersion model using the same ray coverage and the same correlation length ($L_{corr} = 10^7$). As for QADR17, we use a constant *a priori* error for the data, equal to the root-mean-square deviation of the data set. A description of the inversion method can be found in Durand *et al.* (2015). The obtained dispersion maps are shown at periods between 50 and 240 s in Fig. C6.

To compare Rayleigh-wave attenuation and dispersion, we expand our attenuation and dispersion maps in spherical-harmonics. We note a global correlation (for all degrees and orders together) between seismic velocity and attenuation for all periods. The global correlation coefficients are above 0.34 for all periods between 50 and 200 s. This suggests a common origin to both attenuation and velocity perturbations. For the period of 100 s, Fig. 11 depicts the spherical harmonic coefficients of the attenuation $\tilde{Q} = \delta \ln(Q)/\ln(Q)$ as a function of those of the phase velocity $\tilde{c} = \delta c/c$. The correlation between these two quantities is $C = 0.52$ and their variances $\sigma_{\tilde{Q}}^2 = 1.21 \times 10^{-2}$ and $\sigma_{\tilde{c}}^2 = 1.27 \times 10^{-4}$. The proportionality factor between these two quantities should be in between $C\sigma_{\tilde{Q}}/\sigma_{\tilde{c}} = 5.11$ (least-square estimate assuming that the velocities
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Comparison of QADR17 (left column) with two sets of attenuation maps: QMA16 (middle column, Ma *et al.* 2016) and DEA06a (right column, Dalton & Ekström 2006). At each period we display the perturbation in $\delta Q^{-1}/Q^{-1}$ in per cent according to a reference $Q^{-1}$ value which is indicated on each map. The three rows correspond to observations at periods of 50 s, 100 s and 150 s, respectively.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Correlations as a function of spherical harmonic degree between QADR17, QMA16 (Ma *et al.* 2016) and DEA06a (Dalton & Ekström 2006). DEA06a is only given up to degree 12.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Coefficients of $\delta\ln(Q)/\ln(Q)$ as a function of those of $\delta c/c$, measured at a period of 100 s. The regression suggests a slope of 9.74 between these two quantities (solid line). The dotted lines have slopes of 5.11 and 18.56. The correlation between the variables is 0.52.</figcaption></figure>

are perfectly known) and $\sigma_{\tilde{Q}}/(\sigma_c C) = 18.56$ (least-square estimate assuming that the attenuations are perfectly known). As we do not properly estimate the uncertainties of the two tomographic models to perform a sophisticated statistical adjustment (Deming 1964), we can roughly consider that the proportionality factor between the two quantities is the geometric average of the two previous estimates:

$$ \frac{d \ln(Q)}{d c/c} = \frac{\sigma_{\tilde{Q}}}{\sigma_c} \approx 9.74. \qquad (8) $$

This ratio is similar to what can be estimated assuming that both velocity and attenuation are functions of temperature (Nataf & Ricard 1996). In the same way that scaling factors between P and S velocities or between density and S velocity are often imposed in 3-D tomographic models (e.g. Ishii & Tromp 2001; Durand et al. 2016), this ratio between quality factor and velocities could also be imposed. The inversion of the attenuation with depth will likely provide a more precise indication of the attenuation-velocity relations.

We can then compute at each period the correlation between QADR17 and the corresponding phase velocity map as a function of spherical-harmonic degree (see Fig. 12). Correlations are above the 95 per cent confidence level (i.e. there is 95 per cent probability that these correlations are not due to chance) for degrees $\ge 4$ at periods smaller than 100 s and $\ge 6$ at longer periods. To first order, this means that velocity and attenuation are consistent for wavelengths smaller than 10 000 and 6600 km, respectively. However, the models are not correlated between degrees 2 and 4 or 6 depending on the period. This indicates a strong disagreement between velocity and attenuation at very long wavelengths.

To identify the regions responsible for this disagreement, we present in Fig. 13 the attenuation and dispersion maps at 100 s.

Panels (a) and (b) display unfiltered maps, panels (c) and (d) only include degrees 1, 2 and 3 and panels (e) and (f) are for degrees $\ge 4$. At first sight, the unfiltered attenuation and dispersion maps show similar patterns with high velocity and low attenuation under continents, and low velocity and high attenuation under oceans. There are however significant differences. Beneath continents, some Phanerozoic regions (eastern Asia, western north-America and the northeast of Africa) have a broader signature in velocity than in attenuation. In oceanic regions, a clear low-velocity anomaly follows the mid-ocean ridges and the velocity increases with seafloor age. This is for example very clear in the Pacific ocean (Fig. 13b). This evolution with age is less clear in the attenuation map and the attenuation in the central part of the Pacific, in the southern Indian ocean or in the Wharton basin is as strong as under ridges (Fig. 13a). The filtered maps (panels (c) and (d)) highlight these discrepancies. Fig. 13(c) reveals a strong attenuating anomaly in the middle of the Pacific ocean which encompasses most Pacific hotspots, whereas the strongest low-velocity anomaly is centred beneath the Pacific ridge. There is a difference in the sign of the perturbation (high Q and low velocity) under a broad region extending from the northern part of Africa to eastern Asia, Indonesia and the Philippines. This region covers entirely the southern part of Asia. The differences observed under the Pacific ocean and from Africa to eastern Asia explain the poor correlations observed for degrees $\le 3$ in Fig. 12. Figs 13(e) and (f) shows similar velocity and attenuation patterns for spherical harmonic degrees $\ge 4$, in agreement with correlations obtained in Fig. 12.

The differences at very long wavelengths (degrees $\le 3$) are observed at all periods although the strength of the anomalies is maximum at 50 and 100 s (see Fig. C8). From the sensitivity kernels (Fig. 1) the physical processes responsible for these differences probably originate in the depth range 70–200 km, which corresponds to the oceanic asthenosphere and the base of the continental lithosphere. We are confident that these differences are not biased by ray density, because the very long wavelengths of our models are well constrained by our data set, and because the Rayleigh-wave attenuation and dispersion maps are built using the same ray coverage. In the following section, we discuss some possible origins of these differences.

## 4.3 Interpretation

Our filtered attenuation map (Fig. 13c) suggests that the Pacific high-attenuation anomaly has a very broad component within the oceanic asthenosphere. It could be rooted deeper, as we still observe the high attenuation anomaly at 200 s north of Fiji-Tonga subduction (Fig. C8). This is in good agreement with the attenuation maps of Romanowicz & Gung (2002) at 300 km, where they interpret a

<figure><img src="image_14.png" /><figcaption>**Figure 12.** Correlation for every spherical-harmonic degree between our model QADR17 and a dispersion model built with the same ray density.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Figure 13.** Maps of lateral variations of *Q* (left) and of phase velocity (right) at a period of 100 s. We plot raw maps in panels (a) and (b), we filter out spherical harmonic degrees larger than 3 in panels (c) and (d) and we filter out spherical harmonic degrees lower than 4 in panels (e) and (f). Averaged attenuation and velocity (in km s⁻¹) are indicated above the scales.</figcaption></figure>

region of high attenuation in the central Pacific extending from south of the equator to Hawaii, as the deflection under the lithosphere of a thermal plume.

Thermal plumes are expected to produce anomalies at least 1000 km in diameter, when they spread horizontally in the low-viscosity asthenosphere (Davies & Richards 1992). A number of broad plumes have recently been detected in global S-wave tomography beneath some major hotspots in the Pacific ocean (French & Romanowicz 2015). The very broad attenuation signature observed in the Pacific ocean (Fig. 13c) could therefore be produced by several plume anomalies in the asthenosphere (in Figs 13a and e, most Pacific hotspots are located in attenuating regions). In addition, recent shear-velocity models show a strong low-velocity layer associated with the oceanic asthenosphere (e.g. Debayle & Ricard 2012) and this layer is particularly pronounced in the Pacific ocean (e.g. Ritzwoller et al. 2004; Maggi et al. 2006). At periods smaller than 200 s, we observe on the dispersion maps that seismic velocities gently increase with seafloor age in the Pacific ocean (Fig. C6). However, at periods greater than 70 s (Fig. C6), the Pacific hotspots are generally located in low-velocity regions.

As attenuation and velocity are at first order related to temperature (see eq. 8), it is likely that temperature contributes to the positive correlation between attenuation and velocity maps for degrees higher than 6 (Figs 12, and 13e and f). Temperature could for example explain a number of hotspot anomalies (Afar, Reunion, Marion, Kerguelen, Tristan da Cuhna, Cape Verde, etc.) which are associated with high attenuation and low velocities.

In addition to temperature, seismic velocities and attenuation are sensitive to water content, melt fraction, major element chemistry and grain size. Experimental studies show that water is present in the mantle as hydrogen incorporated in anhydrous mineral (Karato 2003). In upper mantle minerals, an increase in water content produces an increase in attenuation and a decrease in velocity

(Karato 1995, 2003). Water may therefore also contribute to explaining the strong attenuation observed in the central Pacific.

The interpretation of the Pacific attenuation anomaly as resulting from several thermal plumes is appealing, because it would also explain the strong radial anisotropy which has been observed in the Pacific ocean by Ekström & Dziewonski (1998) with a lateral extension similar to what we observe in Fig. 13(c). Radial anisotropy means that S-waves polarized in the horizontal direction (SH-waves) and in the vertical direction (SV-waves) propagate with different velocities. Ekström & Dziewonski (1998) find that SH-waves are faster than SV-waves in the Pacific ocean, an observation which is commonly interpreted as due to the horizontal alignment of anisotropic crystals in an horizontal shear-flow (Montagner 1994). A confinement of the horizontal shearing might be due to the addition of low viscosity material brought by the hotspots. This would explain both the strong radial anisotropy and the azimuthal anisotropy, remarkably well aligned with plate motion beneath the fast-moving Pacific plate (Debayle & Ricard 2013).

If the upwelling of hot and wet material is a plausible explanation for the central Pacific attenuation anomaly, the seismic signatures of most ocean ridges remain surprising. This is particularly true for the East-Pacific Rise which is associated with very low velocities but not with a major high-attenuation signal. Our findings are however in agreement with the MELT experiment which reports in the same area, much less attenuation than what is predicted by thermal models (Yang et al. 2007). The presence of melt under ridges may therefore play an important role, but this would imply that melt affects the seismic velocity more than the attenuation which is not what is observed experimentally (Faul & Jackson 2015). However, the hypothesis that melt affects the seismic velocities more that the attenuation has been advocated on the base of seismic observations (Forsyth et al. 1998) and mechanical models (Hammond & Humphreys 2000a,b).
Note finally that partial-melt is also a plausible explanation for other regions where low velocities are associated with high *Q*. This is for example observed beneath the Red Sea and in East Asia (Figs C8a and b), where low-velocity zones in the upper mantle can be explained by melting of wet plumes at the base of lithosphere (Richard & Iwamori 2010).

## 5 CONCLUSIONS

We have presented QADR17, a set of global attenuation maps obtained in the period range 50–240 s from the analysis of 372 629 fundamental Rayleigh-wave attenuation curves published by Debayle & Ricard (2012). We invert for the logarithm of *Q*, a parametrization which brings our data set closer to a Gaussian distribution, avoids negative values and allows for the known large range of variations of the quality factor. This parametrization may explain part of the differences observed between QADR17 and some other recent published models.

To first order, we find that *Q* varies with lithospheric age (young oceans have higher attenuation than old continents). Contrary to what is commonly observed for velocity models, the mid-ocean ridge signature does not dominate the attenuation pattern of oceanic regions. Attenuating regions are rather well correlated with hotspots and some of them are associated with a very strong attenuation signal which persists at periods greater than 140 s.

We compare QADR17 maps with a set of phase velocity maps built using the same approach and the same Rayleigh-wave data set and find no correlation at degrees 3, corresponding to wavelengths of ~13 000 km. The long wavelength component of our maps reveal a broad attenuating anomaly located in the middle of the Pacific ocean and which encompasses most Pacific hotspots, whereas the strongest low velocity anomaly is centred beneath the Pacific ridge. High *Q* and low velocities are found under a broad region extending from the northern part of Africa to eastern Asia, Indonesia and the Philippines.

Sensitivity kernels suggest that the physical processes responsible for these differences originate in the depth range 70–200 km, which corresponds to the oceanic asthenosphere and the base of the continental lithosphere. The most likely explanation for the broad attenuating anomaly located in the middle of the Pacific ocean is the existence of several thermal plume anomalies that would pond in the asthenosphere. This interpretation provides an explanation for the broad radial anisotropy anomaly which has been observed in the Pacific ocean (Ekström & Dziewonski 1998) and is consistent with a previous study by Romanowicz & Gung (2002). The strong low velocity signature of the East Pacific Rise is difficult to reconcile with the sole effects of temperature. The presence of partial melt has been invoked to explain a strong velocity reduction associated with a moderate attenuation anomaly (Yang *et al.* 2007) beneath a small portion of the East Pacific Rise. Our observations of a strong reduction in phase velocities combined with a relatively weak attenuation suggest that partial melting occurs at larger scale along the entire East Pacific Rise, and possibly in other regions beneath the Red Sea and the eastern part of Asia.

## ACKNOWLEDGEMENTS

This work was supported by the French ANR SEISGLOB ANR-11-BLANC-SIMI5-6-016-01. We thank the IRIS and RESIF data centres for providing broadband data. We thank Jean-Jacques Lévêque and Luis Rivera for helpful discussions, Gaby Laske and two anonymous reviewers for thorough and constructive reviews that helped to improve an earlier version of the manuscript.

## REFERENCES

Aki, K. & Richards, P.G., 2002. *Quantitative Seismology*, 2nd edn, University Science Books.

Billien, M., Leveque, J.J. & Trampert, J., 2000. Global maps of rayleigh wave attenuation for periods between 40 and 150 seconds, *Geophys. Res. Lett.*, **27**(22), 3619–3622.

Cara, M., 1979. Lateral variations of S-velocity in the upper mantle from higher Rayleigh modes, *Geophys. J. R. astr. Soc.*, **57**(3), 649–670.

Cara, M. & Lévêque, J.J., 1987. Waveform inversion using secondary observables, *Geophys. Res. Lett.*, **14**(10), 1046–1049.

Dalton, C.A. & Ekström, G., 2006. Global models of surface wave attenuation, *J. geophys. Res.*, **111**, B05317, doi:10.1029/2005JB003997.

Dalton, C.A., Ekström, G. & Dziewonski, A.M., 2008. The global attenuation structure of the upper mantle, *J. geophys. Res.*, **113**, B09303, doi:10.1029/2007JB005429.

Dalton, C.A., Hjörleifsdóttir, V. & Ekström, G., 2014. A comparison of approaches to the prediction of surface wave amplitude, *Geophys. J. Int.*, **196**(1), 386–404.

Davies, G.F. & Richards, M.A., 1992. Mantle convection, *J. Geol.*, **100**(2), 151–206.

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**, B10308, doi:10.1029/2012JB009288.

Debayle, E. & Ricard, Y., 2013. Seismic observations of large-scale deformation at the bottom of fast-moving plates, *Earth planet. Sci. Lett.*, **376**, 165–177.

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

Debayle, E., Dubuffet, F. & Durand, S., 2016. An automatically updated S-wave model of the upper mantle and the depth extent of azimuthal anisotropy, *Geophys. Res. Lett.*, **43**, 674–682.

Deming, W.E., 1964. *Statistical Adjustment of Data*, Dover Publications.

Durand, S., Chambat, F., Matas, J. & Ricard, Y., 2012. Constraining the kinetics of mantle phase changes with seismic data, *Geophys. J. Int.*, **189**(3), 1557–1564.

Durand, S., Debayle, E. & Ricard, Y., 2015. Rayleigh wave phase velocity and error maps up to the fifth overtone, *Geophys. Res. Lett.*, **42**(9), 3266–3272.

Durand, S., Debayle, E., Ricard, Y. & Lambotte, S., 2016. Seismic evidence for a change in the large-scale tomographic pattern across the D'' layer, *Geophys. Res. Lett.*, **43**(15), 7928–7936.

Dziewonski, A. & Anderson, D., 1981. Preliminary reference Earth model, *Phys. Earth planet. Inter.*, **25**(4), 297–356.

Dziewonski, A.M., Chou, T.-A. & Woodhouse, J.H., 1981. Determination of earthquake source parameters from waveform data for studies of global and regional seismicity, *J. geophys. Res.*, **86**(B4), 2825–2852.

Ekström, G. & Dziewonski, A.M., 1998. The unique anisotropy of the Pacific upper mantle, *Nature*, **394**(6689), 168–172.

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

Faul, U. & Jackson, I., 2015. Transient creep and strain energy dissipation: an experimental perspective, *Annu. Rev. Earth Planet. Sci.*, **43**, 541–569.

Ferreira, A.M.G. & Woodhouse, J.H., 2007. Source, path and receiver effects on seismic surface waves, *Geophys. J. Int.*, **168**(1), 109–132.

Forsyth, D.W., Webb, S.C., Dorman, L.M. & Shen, Y., 1998. Phase velocities of Rayleigh waves in the MELT experiment on the East Pacific Rise, *Science*, **280**(5367), 1235–1238.

French, S.W. & Romanowicz, B., 2015. Broad plumes rooted at the base of the Earth's mantle beneath major hotspots, *Nature*, **525**(7567), 95–99.

Gung, Y. & Romanowicz, B., 2004. *Q* tomography of the upper mantle using three-component long-period waveforms, *Geophys. J. Int.*, **157**(2), 813–830.
Hammond, W.C. & Humphreys, E.D., 2000a. Upper mantle seismic wave attenuation: effects of realistic partial melt distribution, *J. geophys. Res.*, **105**(B5), 10 987–10 999.

Hammond, W.C. & Humphreys, E.D., 2000b. Upper mantle seismic wave velocity: effects of realistic partial melt geometries, *J. geophys. Res.*, **105**(B5), 10 975–10 986.

Ishii, M. & Tromp, J., 2001. Even-degree lateral variations in the Earth's mantle constrained by free oscillations and the free-air gravity anomaly, *Geophys. J. Int.*, **145**(1), 77–96.

Jackson, I., 2007. Properties of Rocks and Minerals – Physical Origins of Anelasticity and Attenuation in Rock, in *Treatise on Geophysics*, ed. Schubert, G., pp. 493–525, Elsevier.

Kanamori, H. & Anderson, D.L., 1977. Importance of physical dispersion in surface wave and free oscillation problems: review, *Rev. Geophys.*, **15**(1), 105–112.

Karato, S.-I., 1993. Importance of anelasticity in the interpretation of seismic tomography, *Geophys. Res. Lett.*, **20**(15), 1623–1626.

Karato, S.-I., 1995. Effects of Water on Seismic Wave Velocities in the Upper Mantle, *Proc. Japan Acad. B*, **71**(2), 61–66.

Karato, S.-I., 2003. Mapping water content in the upper mantle, in, *Inside the Subduction Factory*, Geophysical Monograph Series, vol. 138, pp. 135–152, ed. Eiler, J., American Geophysical Union.

Lay, T. & Kanamori, H., 1985. Geometric effects of global lateral heterogeneity on long-period surface wave propagation, *J. geophys. Res.*, **90**(B1), 605–621.

Ma, Z., Masters, G. & Mancinelli, N., 2016. Two-dimensional global Rayleigh wave attenuation model by accounting for finite-frequency focusing and defocusing effect, *Geophys. J. Int.*, **204**(1), 631–649.

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**(3), 1384–1397.

Meschede, M. & Romanowicz, B., 2015. Lateral heterogeneity scales in regional and global upper mantle shear velocity models, *Geophys. J. Int.*, **200**(2), 1076–1093.

Montagner, J., 1986. Regional 3-Dimensional Structures Using Long-Period Surface-Waves, *Ann. Geophys. B*, **4**(3), 283–294.

Montagner, J.-P., 1994. Can seismology tell us anything about convection in the mantle?, *Rev. Geophys.*, **32**(2), 115–137.

Müller, R.D., Royer, J.-Y. & Lawver, L.A., 1993. Revised plate motions relative to the hotspots from combined Atlantic and Indian Ocean hotspot tracks, *Geology*, **21**(3), 275–278.

Nataf, H.-C. & Ricard, Y., 1996. 3SMAC: ana priori tomographic model of the upper mantle based on geophysical modeling, *Phys. Earth planet. Inter.*, **95**(1), 101–122.

Parisli, L. & Ferreira, A. M.G., 2016. Empirical assessment of the validity limits of the surface wave full ray theory using realistic 3-D Earth models, *Geophys. J. Int.*, **205**(1), 146–159.

Ricard, Y., Nataf, H.-C. & Montagner, J.-P., 1996. The three-dimensional seismological model a priori constrained: confrontation with seismic data, *J. geophys. Res.*, **101**(B4), 8457–8472.

Ricard, Y., Durand, S., Montagner, J.P. & Chambat, F., 2014. Is there seismic attenuation in the mantle?, *Earth planet. Sci. Lett.*, **388**, 257–264.

Richard, G.C. & Iwamori, H., 2010. Stagnant slab, wet plumes and Cenozoic volcanism in East Asia, *Phys. Earth planet. Inter.*, **183**(1–2), 280–287.

Ritsema, J., Deuss, A., van Heijst, H.J. & Woodhouse, J.H., 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**(3), 1223–1236.

Ritzwoller, M.H., Shapiro, N.M. & Zhong, S.-J., 2004. Cooling history of the Pacific lithosphere, *Earth planet. Sci. Lett.*, **226**(1–2), 69–84.

Romanowicz, B., 1990. The upper mantle degree 2: Constraints and inferences from global mantle wave attenuation measurements, *J. geophys. Res.*, **95**(B7), 11 051–11 071.

Romanowicz, B., 1995. A global tomographic model of shear attenuation in the upper mantle, *J. geophys. Res.*, **100**(B7), 12 375–12 394.

Romanowicz, B. & Gung, Y., 2002. Superplumes from the core-mantle boundary to the lithosphere: implications for heat flux, *Science*, **296**(5567), 513–516.

Selby, N.D. & Woodhouse, J.H., 2002. The *Q* structure of the upper mantle: constraints from Rayleigh wave amplitudes, *J. geophys. Res.*, **107**(B5), ESE 5-1-ESE 5-11.

Shito, A., Karato, S.-I., Matsukage, K. & Nishihara, Y., 2006. Towards mapping the three-dimensional distribution of water in the upper mantle from velocity and attenuation tomography, in *Earth's Deep Water Cycle*, Geophysical Monograph Series, vol. 168, pp. 225–236, eds Jacobsen, S.D. & Van Der Lee, S., American Geophysical Union, Washington, DC.

Tarantola, A. & Valette, B., 1982. Generalized nonlinear inverse problems solved using the least squares criterion, *Rev. Geophys.*, **20**(2), 219–232.

Um, J. & Dahlen, F.A., 1992. Source phase and amplitude anomalies of long-period surface waves, *Geophys. Res. Lett.*, **19**(15), 1575–1578.

Warren, L.M. & Shearer, P.M., 2002. Mapping lateral variations in upper mantle attenuation by stacking P and PP spectra, *J. geophys. Res.*, **107**(B12), 2342, doi:10.1029/2001JB001195.

Woodhouse, J.H. & Wong, Y.K., 1986. Amplitude, phase and path anomalies of mantle waves, *Geophys. J. Int.*, **87**(3), 753–773.

Yang, Y., Forsyth, D.W. & Weeraratne, D.S., 2007. Seismic attenuation near the East Pacific Rise and the origin of the low-velocity zone, *Earth planet. Sci. Lett.*, **258**(1–2), 260–268.

## APPENDIX A: NON-LINEAR INVERSION

The problem to solve is to find the local attenuation $Q(r)$ when the integrated attenuation $Q_i^{-1}$ is known for a number of ray paths $i$. With the great circle approximation, the integrated attenuation is

$$Q_i^{-1} = \int_i Q^{-1}(s) \frac{ds}{\Delta_i}, \quad (A1)$$

where $s$ is the abscissa along the ray of length $\Delta_i$ and $Q^{-1}(s)$ is a shortcut for $Q^{-1}(r(s))$, the line integral is along the path $i$.

To work with a data distribution closer to a Gaussian, we use logarithms and define the data $d_i$ and the parameters $\hat{m}(r)$ as

$$d_i = \ln Q_i^{-1} \text{ and } \hat{m}(r) = \ln \hat{Q}^{-1}(r). \quad (A2)$$

Here, $\hat{m}$ is the model retrieved from our data set, which is not necessarily identical to the real attenuation of the Earth. The problem is to find $\hat{m}$ such that the vector $d$, containing the observations $d_i$, is obtained from the model $\hat{m}$. Practically, solving the problem consists in resolving $d = g(\hat{m})$, which minimizes the distance between the observation $d_i$ and the prediction $d_i^p$ containing the predicted data along the ray $i$ computed from the model $\hat{Q}(r)$:

$$ \begin{aligned} d_i^p &= \ln((Q_i^p)^{-1}) = \ln \left[ \int_i \hat{Q}^{-1}(s) \frac{ds}{\Delta_i} \right] \\ &= \ln \left[ \int_i \exp(\hat{m}(s)) \frac{ds}{\Delta_i} \right] = g_i(\hat{m}). \end{aligned} \quad (A3) $$

Following Tarantola & Valette (1982) we impose a smoothing, regularizing condition and we minimize

$$ S = \sum_i \left[ \frac{d_i^p - d_i}{\sigma_i} \right]^2 + \iint [\hat{m}(\mathbf{r}) - m_0(\mathbf{r})] \cdot C_{m_0}^{-1}(\mathbf{r}, \mathbf{r}') \cdot [\hat{m}(\mathbf{r}') - m_0(\mathbf{r}')] \, d\omega \, d\omega', \quad (A4) $$

that is, we want the prediction $d_i^p$ to fit the data $d_i$ within an uncertainty $\sigma_i$, while the model $\hat{m}(r)$ remains close to a reference model $m_0(r)$ (here a uniform model of quality factor $Q_0$) through a weighting covariance function of the form

$$ C_{m_0}(\mathbf{r}, \mathbf{r}') = \sigma_m^2 \exp \left( - \frac{\Delta(\mathbf{r}, \mathbf{r}')^2}{2L_{\text{corr}}^2} \right). \quad (A5) $$

In (A4), $d\omega$ and $d\omega'$ are two surface elements (solid angles) around the positions $\mathbf{r}$ and $\mathbf{r}'$ separated by the angular distance $\Delta(\mathbf{r}, \mathbf{r}')$, the
two integrals are over the whole sphere. In term of attenuations, the quantity to minimize (A4) can also be written as

$$S = \sum_i \left[ \frac{\ln(Q_i / Q_i^p)}{\sigma_i} \right]^2 + \iint \ln(Q_0 / \hat{Q}(\mathbf{r})) \cdot C_{m_0}^{-1}(\mathbf{r}, \mathbf{r}') \cdot \ln(Q_0 / \hat{Q}(\mathbf{r}')) \, d\omega \, d\omega'. \quad (A6)$$

Following Tarantola & Valette (1982), the solution of this inverse problem is

$$\hat{m} = m_0 + \mathbf{C}_{m_0} \cdot \mathbf{G}^T \cdot \mathbf{S}^{-1} \mathbf{D} \quad (A7)$$

with

$$\mathbf{S} = \mathbf{C}_{d_0} + \mathbf{G} \cdot \mathbf{C}_{m_0} \cdot \mathbf{G}^T \quad (A8)$$

and

$$D = d - g(\hat{m}) + \mathbf{G} \cdot (\hat{m} - m_0), \quad (A9)$$

where $T$ denotes the transpose, $m_0$ is the a priori model, and the matrix $\mathbf{G}$ (formally $\partial g / \partial \hat{m}$) contains the partial derivatives of the predictions with respect to the model.

The implementation of an apparently discrete formalism (A7)–(A9) is somewhat tricky for continuous functions. For our problem the $\mathbf{G}$ matrix is in fact the operator with lines:

$$\mathbf{G}_i = \int_i \frac{Q_i^p}{\hat{Q}(s)} \cdot \frac{ds}{\Delta_i} \quad (A10)$$

(the • remains that applied to a model $m(\mathbf{r})$, with means multiplied by $m(s)$, it leads to $i$ values). Eqs (A9) and (A8) can therefore be written as

$$D_i = - \ln \left[ \int_i \frac{Q_i}{\hat{Q}(s)} \frac{ds}{\Delta_i} \right] + \int_i \frac{Q_i^p}{\hat{Q}(s)} \ln \left[ \frac{Q_0}{\hat{Q}(s)} \right] \frac{ds}{\Delta_i} \quad (A11)$$

$$S_{ij} = \sigma_i^2 \delta_{ij} + \sigma_m^2 \int_i \int_j \frac{Q_i Q_j}{\hat{Q}(s) \hat{Q}(s')} \times \exp\left(-\frac{\Delta(\mathbf{r}(s), \mathbf{r}'(s'))^2}{2L_{\text{corr}}^2}\right) \frac{ds}{\Delta_i} \frac{ds'}{\Delta_j}, \quad (A12)$$

Finally $\mathbf{C}_{m_0} \cdot \mathbf{G}^T$, is a smoothing operator, which applies to a vector $\mathbf{X}$ with dimension the number of rays (i.e. $\mathbf{X} = \mathbf{S}^{-1} \mathbf{D}$), to

generate a continuous function on the sphere, that is,

$$\mathbf{C}_{m_0} \cdot \mathbf{G}^T \cdot \mathbf{X} = \sigma^2 \sum_i X_i \int_i \frac{Q_i^p}{\hat{Q}(s')} \exp\left(-\frac{\Delta(\mathbf{r}, \mathbf{r}'(s'))^2}{2L_{\text{corr}}^2}\right) \frac{ds'}{\Delta_i} \quad (A13)$$

This operator gives the estimate for the attenuation at position $\mathbf{r}$. Classically, this inversion is proceeded iteratively where an estimate $\hat{m}$ (i.e. a model $Q(\mathbf{r}))$, is used to compute the operators of the right side of (A7) which provides an improved estimate $\hat{m}$ (left side of (A7).

## APPENDIX B: REJECTING OUTLIERS IN EVERY CLUSTER

This data selection is included in the clustering step performed before the regionalization. We cluster paths corresponding to epi-centres close together recorded at a single station, assuming waves propagate in the same attenuation structure. If an attenuation curve differs too much from other curves in the cluster, we suspect that this observation is problematic. The aim of this data selec-tion is to identify these outliers and to exclude them from each cluster.

First, we compute at each period $T_k$ the average quality factor for the cluster $\overline{\ln Q(T_k)}$ and its variance $\sigma_g(T_k)$. Second, we average $\sigma_g(T_k)$ over the periods $T_k$:

$$\sigma^2 = \frac{1}{p} \sum_{k=1}^{p} \sigma_g^2(T_k) \quad (B1)$$

where $p$ is the number of periods. Third, we compute a distance between each attenuation curve $i$ of the cluster and the average attenuation curve:

$$R_i^2 = \frac{1}{p} \sum_{k=1}^{p} \left( \ln(Q_i(T_k)) - \overline{\ln(Q(T_k))} \right)^2. \quad (B2)$$

If $R_i^2$ is higher than $\sigma^2$, we consider the attenuation curve $i$ as an outlier. The curve $i$ is discarded if the remaining number of paths in the cluster is larger than 2.
APPENDIX C: FIGURES

<figure><img src="image_16.png" /><figcaption>**Figure C1.** Maps of lateral variations of *Q*. Left column: maps are corrected for all effects except focusing (see eq. 1); right column: maps also include the focussing correction. Although the focusing corrections are visible on the maps, these corrections remain moderate. Hotspot locations are in blue (Müller *et al.* 1993), plates in green.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>**Figure C2.** Maps of attenuation for 70 s Rayleigh waves beneath the East Pacific Rise (top row) and for 140 s Rayleigh waves beneath North America (bottom row). The colour scale shows the perturbation in attenuation $\delta Q^{-1}$ according to the reference indicated above the colour scale (similar plots can be found in Dalton & Ekström 2006). In each row, the left column displays maps without the focusing correction; the middle column, with the focussing correction; the right column displays the effect of focusing (i.e. left column minus middle column).</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>**Figure C3.** QADR17 (right column) is compared with the result of an inversion (left column) in which we keep paths for which the computed amplitudes vary rapidly with azimuths (i.e. by more than 1 per cent for a take-off azimuth variation of 1°). There are only minor differences between the two inversions, suggesting that most data with large CMT errors have been rejected by the requirement of an amplitude factor < 2 between synthetic and actual waveforms.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>**Figure C4.** Top: maps of lateral variations of *Q* at 70, 140 and 200 s. Left panel displays the maps obtained without clustering the paths and the right panel displays QADR17, where nearby paths are clustered. Bottom: correlation between the maps above, at 70, 140 and 200 s period. The dashed lines indicate the 95 per cent significance levels.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>**Figure C5.** Lateral variations in Rayleigh-wave attenuation at a period of 100 s obtained using different correlation lengths $L_{corr}$. Other *a priori* informations are identical to our preferred model QADR17. We plot values of *Q* with a logarithmic scale.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>**Figure C6.** Lateral variations of Rayleigh phase velocities at different periods. Perturbation in percent are shown according to a reference value (in km s⁻¹) which is indicated under each map.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>**Figure C7.** Same as Fig. 5 except that we plot attenuation in ocean basins as a function of the square root of sea floor age. Up to 100 s periods, the attenuation decreases rather linearly with the root square of age.</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>**Figure C8.** Maps of lateral variations of *Q* (left) and phase velocities (right) plotted after filtering out spherical harmonic degrees ≥ 4. At each period we plot *Q* with a logarithmic scale and velocity perturbation in per cent from the globally average phase velocity.</figcaption></figure>