TY - JOUR
AU - Romanowicz,, Barbara
AB - SUMMARY We present a global upper-mantle shear wave attenuation model that is built through a hybrid full-waveform inversion algorithm applied to long-period waveforms, using the spectral element method for wavefield computations. Our inversion strategy is based on an iterative approach that involves the inversion for successive updates in the attenuation parameter (|$\delta Q^{-1}_\mu$|) and elastic parameters (isotropic velocity VS, and radial anisotropy parameter ξ) through a Gauss–Newton-type optimization scheme that employs envelope- and waveform-type misfit functionals for the two steps, respectively. We also include source and receiver terms in the inversion steps for attenuation structure. We conducted a total of eight iterations (six for attenuation and two for elastic structure), and one inversion for updates to source parameters. The starting model included the elastic part of the relatively high-resolution 3-D whole mantle seismic velocity model, SEMUCB-WM1, which served to account for elastic focusing effects. The data set is a subset of the three-component surface waveform data set, filtered between 400 and 60 s, that contributed to the construction of the whole-mantle tomographic model SEMUCB-WM1. We applied strict selection criteria to this data set for the attenuation iteration steps, and investigated the effect of attenuation crustal structure on the retrieved mantle attenuation structure. While a constant 1-D Qμ model with a constant value of 165 throughout the upper mantle was used as starting model for attenuation inversion, we were able to recover, in depth extent and strength, the high-attenuation zone present in the depth range 80–200 km. The final 3-D model, SEMUCB-UMQ, shows strong correlation with tectonic features down to 200–250 km depth, with low attenuation beneath the cratons, stable parts of continents and regions of old oceanic crust, and high attenuation along mid-ocean ridges and backarcs. Below 250 km, we observe strong attenuation in the southwestern Pacific and eastern Africa, while low attenuation zones fade beneath most of the cratons. The strong negative correlation of |$Q^{-1}_\mu$| and VS anomalies at shallow upper-mantle depths points to a common dominant origin for the two, likely due to variations in thermal structure. A comparison with two other global upper-mantle attenuation models shows promising consistency. As we updated the elastic 3-D model in alternate iterations, we found that the VS part of the model was stable, while the ξ structure evolution was more pronounced, indicating that it may be important to include 3-D attenuation effects when inverting for ξ, possibly due to the influence of dispersion corrections on this less well-constrained parameter. Waveform inversion, Seismic attenuation, Seismic tomography 1 INTRODUCTION When earthquake waves travel through the Earth’s mantle, they lose energy due to anelastic attenuation processes, which manifest themselves through amplitude decay. Mapping lateral variations in seismic attenuation holds great potential to provide constraints on the distribution of temperature anomalies, melt and water content in the mantle, and therefore to shed light on the mechanisms that drive mantle dynamics (e.g. Romanowicz 1990; Karato 1993; Jackson et al.2004; Dalton & Faul 2010). These constraints are complementary to those of elastic tomography, because of different sensitivities of velocity and attenuation to temperature and composition. For example, attenuation depends exponentially on temperature (Arrhenius law), whereas the dependence of velocity is more linear. Moreover, anelastic attenuation causes velocity dispersion, therefore constraining attenuation is important for improving the accuracy of elastic velocity models. However, extracting the anelastic signal from seismic data is a challenging task, because seismic wave amplitudes are also affected by scattering and focusing effects due to propagation in complex elastic structures (e.g. Woodhouse & Wong 1986; Romanowicz et al.1987; Durek et al.1993; Zhou 2009; Dalton et al.2014; Bao et al.2016), and by uncertainties in earthquake source parameters. This is why global attenuation tomography has been lagging behind elastic tomography in the last 20 yr. On the other hand, there are multiple studies on the global average shear attenuation profile with depth in the earth’s upper mantle, based on normal mode, surface wave and body wave measurements (e.g. Anderson et al.1965; Sailor & Dziewonski 1978; Anderson & Hart 1978; Widmer et al.1991; Durek & Ekström 1996; Resovsky et al.2005; Lawrence & Wysession 2006). From these studies, some robust features have emerged: Shear wave attenuation is low in the lithosphere. There exists a high-attenuation zone within the depth range of ∼80–220 km that approximately coincides with the low velocity channel in the upper mantle. Below 200 km, shear wave attenuation decreases sharply and remains low in the transition zone. These models also agree that bulk attenuation is much lower than shear wave attenuation in the mantle although there may be a need for a zone of high bulk attenuation to explain radial normal-mode data. There is no agreement on the location of this zone, which could be in the core (e.g. Anderson & Hart 1978; Dziewonski & Anderson 1981), in the upper mantle (e.g. Sailor & Dziewonski 1978) or in the asthenosphere (e.g. Durek & Ekström 1995). For further details and for a historical perspective on the development of various global average attenuation models, we refer to the reviews by Romanowicz & Durek (2000) and Romanowicz & Mitchell (2015). Existing 3-D global seismic attenuation models of the upper mantle have been constructed using either long-period surface wave amplitude or waveform data (e.g. Romanowicz 1990, 1995; Selby & Woodhouse 2002; Gung & Romanowicz 2004; Dalton et al.2008; Ma et al.2016), or differential amplitude measurements of body waves, such as S–SS (Bhattacharyya et al.1996) and P–PP (Warren & Shearer 2002). We here focus on attenuation modeling using long-period fundamental-mode and overtone surface waveforms. One of the central issues in long-period attenuation tomography is how to account for (i) uncertainties in the source contribution, and (ii) elastic focusing effects on seismic wave amplitudes. In early work, Romanowicz (1994, 1995) tried to minimize these effects by using consecutive Rayleigh wave trains (R1, R2 and R3) simultaneously for each record considered, and rejecting data visibly contaminated by focusing effects. She recovered a low-degree (up to the angular degree 5 in spherical harmonics expansion) shear attenuation upper-mantle model. As shown later by Selby & Woodhouse (2000), at long periods (e.g. 146 s) focusing effects do not significantly influence the seismic wave amplitudes. This makes the mapping of long-wavelength attenuation structure possible using amplitudes of the first arriving long-period surface wave trains (R1/R2, G1/G2), which should be the least contaminated by elastic heterogeneities because the paths traveled are shortest. Based on these results, Selby & Woodhouse (2002) and Gung & Romanowicz (2004) built seismic attenuation models with maximum angular order of 8, without including corrections for focusing effects. Gung & Romanowicz (2004) included constraints from both fundamental and overtone waveforms and checked that including focusing effects based on a degree-24 elastic tomographic model of the upper mantle did not affect their results. In order to obtain higher resolution attenuation images, focusing effects need to be accounted for. To that end, Dalton & Ekström (2006) developed an approach based on inversion of fundamental-mode Rayleigh wave amplitude data, correcting for focusing effects at each frequency using ray theory, and applied it to model both phase velocity and attenuation structures. They showed that it is also important to account for frequency-dependent source and receiver terms. This approach was used successively with different data sets, forward modeling approaches (great-circle ray approximation, exact ray theory and finite-frequency theory) and 2-D finite-frequency kernels (e.g. Dalton et al.2008; Ma et al.2016; Bao et al.2016) mapping consistent attenuation anomalies. Applying this technique, Dalton et al. (2017) showed that the path coverage and addressing the focusing effects are crucial for increasing the resolution of global attenuation models. Most recently, Adenis et al. (2017a,b) built a global upper-mantle shear attenuation model using a data set of fundamental Rayleigh wave attenuation curves measured by Debayle & Ricard (2012) with a method that relies on the secondary observables (Cara & Lévêque 1987). In these studies, the attenuation signal was extracted by considering frequency-dependent focusing effects following a similar approach to Dalton & Ekström (2006). One common feature of all the global attenuation models built so far is the use of approximate techniques for seismic wave propagation simulations. These techniques include approximations based on exact ray theory (e.g. Woodhouse & Wong 1986; Park 1987), asymptotic normal-mode perturbation theory (e.g. Romanowicz 1987; Li & Romanowicz 1995) or finite-frequency theory through single scattering approximation (Zhou et al.2004). All these approximate techniques introduce errors and uncertainty on the recovered seismic attenuation models, mostly because of their inaccuracy in modeling the influence of elastic heterogeneity on seismic wave amplitudes (e.g. Dalton & Ekström 2006; Bao et al.2016). Moreover, these corrections are generally based on smooth elastic models, in which lateral gradients may not be sufficiently well constrained for accurate focusing predictions. What these approximate techniques lack in accuracy can now be addressed by using the spectral element method (SEM) for seismic wavefield computations in 3-D models (Komatitsch & Vilotte 1998; Komatitsch & Tromp 2002). The accuracy of the SEM method extends to the modeling of seismic attenuation, whereby standard linear solids (Liu et al.1976) are used to model the Earth as a linearly anelastic medium. Modeling of the Earth’s mantle as a linear anelastic medium is considered as an acceptable assumption on the grounds of small deformations with respect to the wavelengths of the seismic waves being simulated at distances far enough from the earthquake sources for the non-linear effects to be negligible. Recently, Zhu et al. (2013) presented the first continental-scale attenuation model for Europe and North Atlantic, based on SEM simulations and an adjoint inversion approach. Since the development of the 3-D upper-mantle elastic velocity model SEMum (Lekić & Romanowicz 2011), our group has been using SEM for computing predicted teleseismic waveforms, in the context of a time-domain full-waveform inversion technique. This approach, which includes the computation of Fréchet derivatives using a normal-mode perturbation formalism, that is, non-linear asymptotic coupling theory (NACT, Li & Romanowicz 1995), has been used to build several global 3-D elastic, radially anisotropic models, first for the upper mantle, using long-period surface waveforms down to 60 s period (Lekić & Romanowicz 2011; French et al.2013), and later for the whole mantle, after inclusion of body waveforms down to 32 s period (French & Romanowicz 2014). In building these models, a smoothed version of the 1-D attenuation model QL6 (Durek & Ekström 1996) was considered, and kept fixed throughout the inversion process. We note that in the last 3-D upper-mantle shear attenuation model built by our group, Gung & Romanowicz (2004) employed an early version of our full-waveform inversion technique that relied on the computation of both synthetics and Fréchet derivatives using NACT for the seismic velocity and path average approximation (PAVA, Woodhouse & Dziewonski 1984; Romanowicz 1987) for the seismic attenuation inversions. Their data set included fundamental and overtone minor- and major-arc Rayleigh and Love waves in the period range of 400–60 s. In time-domain waveform inversions, it is important to only consider waveforms that are well aligned in phase with the synthetics of the starting elastic model, which implies careful data selection. For this reason and to try and minimize the contamination by elastic effects, a significant number of paths were eliminated that did not match phase alignment criteria or exhibited excessively strong amplitude anomalies. Gung & Romanowicz (2004) only inverted for spherical harmonic degrees up to l = 8 for which elastic anomalies do not influence the seismic attenuation measurements significantly (e.g. Selby & Woodhouse 2000). As the starting degree-zero (global average) attenuation structure, Gung & Romanowicz (2004) used various models, including the one proposed by Romanowicz (1995). They did not propose an updated 1-D model, because they found significant trade-offs with uncertainties in source moments. In this study, we present a global upper-mantle shear wave attenuation model built using a SEM-based hybrid-full-waveform inversion method. The salient features of the presented study include: Accurate computation of synthetics with SEM under the assumption of linear anelastic behaviour of mantle materials. Using the relatively high-resolution SEMUCB-WM1 (French & Romanowicz 2014), but with a different 1-D Qμ profile as the starting model. Update of source Centroid-Moment-Tensor (CMT) solutions using an NACT-based inversion technique. A stepwise inversion strategy that is based on successive updates of the 3-D shear attenuation model and 3-D radially anisotropic elastic model. Use of envelope waveforms to define misfit in the inversions for attenuation. In the following sections, we present the inversion methodology for elastic and anelastic structures, the data set, the process and results of inversions for source parameters, as well as the final model, which hereinafter is referred to as SEMUCB-UMQ, and compare it to two other existing models based on long-period surface wave data. Finally, we discuss the influence of a heterogeneous crustal attenuation model on SEMUCB-UMQ and present the results of resolution tests performed to quantify the limits on the scales of interpretable structures. 2 METHODOLOGY 2.1 Hybrid full-waveform inversion algorithm We have implemented a hybrid full-waveform inversion algorithm that combines the accuracy of the SEM (Komatitsch & Vilotte 1998) for the computation of the seismic wavefield, with the efficiency of the NACT (Li & Romanowicz 1995), used in the inversion step. SEM is currently the preferred method for the accurate computation of regional and global seismic wavefields in 3-D heterogeneous media. However, this advantage comes at the expense of increased numerical cost, especially for global simulations. One way to reduce the computational cost is to couple SEM with different, less demanding methods, such as normal-mode summation. In this work, as in our previous elastic tomographic models since the work of Lekić & Romanowicz (2011), we use the C-SEM code developed by Capdeville et al. (2003), which couples a 1-D normal-mode computation in the earth’s core with SEM in the 3-D heterogeneous mantle. To make sure that we include major-arc surface waves in our simulations, we computed the seismic wavefield for a time interval of 10 000 s after the origin time. This is long enough to recover (X)R2 and (X)G2 surface wave trains. For the inversion step, we compute the Fréchet derivatives and build the Hessian matrix by employing the NACT approximation. By including coupling between multiplets along and across normal-mode branches, NACT accounts for the variation of the Fréchet derivatives both in depth and along the source-to-receiver path. This better reflects the actual sensitivity of overtones and body waves to structure, compared to more conventional PAVA (Woodhouse & Dziewonski 1984; Romanowicz 1987), in which the Fréchet derivatives vary only in depth, and which is strictly valid only for isolated mode branches, such as the fundamental mode (Mégnin & Romanowicz 1999; Romanowicz et al.2008). The Hessian matrix is built from these 2-D Fréchet derivatives by combining tens of thousands of paths. For further details on NACT, we refer to Li & Romanowicz (1995). Having direct control over the NACT-based-Hessian matrix allows us to employ a rapidly converging Gauss–Newton optimization scheme, thus reducing the number of SEM-based forward wavefield simulations, which is computationally the most expensive step of our algorithm. In the inversion step, we work with time-windowed seismic records that correspond to significant seismic energy arrivals, or ‘wavepackets’ (Li & Romanowicz 1996). The data selection step is carried out through an automated process which compares synthetic waveforms for the current model with observed ones, and selects the wavepackets, following criteria such as signal-to-noise ratio and residual variance thresholds (Panning & Romanowicz 2006), as well as additional criteria in the case of inversion for attenuation, as described in the next section. At each iteration, once we have selected the data set and computed the synthetics and NACT kernels for the current model, we minimize the following misfit functional: \begin{eqnarray} {{\boldsymbol \Phi} (u,m)} &=& {\boldsymbol {\frac{1}{2}}} \left( {{\boldsymbol \Psi} (u,d) C_D^{-1} \Psi (u,d)} \right.\nonumber\\ &&\left. {+\, \alpha (m_k-m_0) C_M^{-1} (m_k-m_0)}\right) \end{eqnarray} (1) where |$\boldsymbol {m_k}$| is the current model that is recovered at the kth iteration, |$\boldsymbol {m_0}$| is the a priori model, which is a 1-D model that represents the degree zero structure of the current model, |$\boldsymbol {C_D}$| and |$\boldsymbol {C_M}$| are the data and model space covariance matrices, |$\boldsymbol {\alpha }$| is a regularization parameter that scales the a priori variance of the model and |$\boldsymbol {\Psi (u,d)}$| is the misfit quantifying the difference between the synthetic (|$\boldsymbol {u}$|) and data (|$\boldsymbol {d}$|) waveforms. We will further describe the covariance matrices in Section 2.6. The misfit functional (|$\boldsymbol {\Phi (u,m)}$|) is iteratively minimized through a Gauss–Newton optimization scheme that updates our model space with perturbations (|$\boldsymbol {\delta m_k}$|) computed as follows (e.g. Tarantola & Valette 1982): \begin{eqnarray} {\delta m_k}& =& - {({G_k^T C_D^{-1} G_k + \alpha C_M^{-1}})^{-1}}\nonumber\\ &&\times {\big ({G_k^T} {\boldsymbol \Psi} (u,d) + \alpha C_M^{-1} (m_k-m_0)\big )} \end{eqnarray} (2) \begin{eqnarray} {G_k} = {\partial {\boldsymbol \Psi} (u,d)}/{\partial m} \end{eqnarray} (3) where |$\boldsymbol {G_k}$| represents the Fréchet derivatives computed using NACT, and updated at each iteration. They will be discussed further in Section 2.2. Once |$\boldsymbol {\delta m_k}$| is computed, we remove the zero-degree structure from it and update our a priori model (|$\boldsymbol {m_0}$|). 2.2 Extension of NACT to attenuation Until now, NACT was used for the computation of the elastic gradient and the approximate Hessian in developing several recent global radially anisotropic velocity models (Lekić & Romanowicz 2011; French et al.2013; French & Romanowicz 2014). The independent parameters progressively updated in these models are the Voigt average isotropic S-wave velocity (VS) and radial anisotropy parameter (|$\xi = V^2_{SH}/V^2_{SV}$|). The perturbations in the remaining four parameters that define a vertically polarized anisotropic medium (VP, |$\phi = V^2_{PV}/V^2_{PH}$|, η and ρ) are scaled to VS and ξ perturbations through empirical scaling factors proposed for the upper mantle by Montagner & Anderson (1989). Further details on this approach and the physical reasoning can be found in appendix A of Panning & Romanowicz (2006). In this work, in addition to the two independent elastic parameters VS and ξ, we add the shear quality factor (Qμ), which represents the dominant contribution to intrinsic attenuation in the mantle compared to the bulk quality factor (Qκ). Further, we follow the assumption of frequency-independent Qμ as originally proposed by Liu et al. (1976) and Kanamori & Anderson (1977). Although there is evidence for frequency dependence of attenuation (e.g. Anderson & Hough 1984; Lekić et al.2009), it is weak for the relatively narrow period range considered in this study, and therefore, the frequency independent Qμ assumption is justified. We have extended NACT to incorporate attenuation perturbations in the computation of Fréchet derivatives. Further details on this can be found in Karaoğlu & Romanowicz (2018). 2.3 Envelope misfit functional In the construction of previous 3-D elastic models of the mantle in our group, we defined the misfit (|$ {{\boldsymbol \Psi} (u,m)}$|) as the point by point time-domain waveform difference (WF) between observed (d(t)) and synthetic (u(t)) windowed and weighted time-series (Gung et al.2003; Gung & Romanowicz 2004; Li & Romanowicz 1996; Panning & Romanowicz 2006): \begin{eqnarray} \Psi (u,d) = d(t)-u(t) \end{eqnarray} (4) Through some detailed synthetic tests, which we summarized in a recent paper (Karaoğlu & Romanowicz 2018), we confirmed the intuition that WF is not a good choice for imaging the attenuation structure in the framework of our hybrid full-waveform inversion technique. This is due to its high sensitivity to the phase differences between the observed and synthetic waveforms, compared to the differences in amplitude. In the same synthetic tests, we also assessed the performance of two other definitions of misfit based on differences between observed and predicted waveform envelopes (E-WF) and spectral amplitudes (SA), respectively. Compared to WF, the other two performed significantly better, showing less leakage from elastic heterogeneities to the attenuation model, leading us to use one of them for this study. After performing a series of tests in the early stages of model updates, we did not find any significant difference between the models recovered by using either of the two methods. In this study, we chose to use E-WF. The E-WF misfit is defined as: \begin{eqnarray} \Psi (u,d) &=& \text{env}(d(t))-\text{env}(u(t)) \nonumber\\ &=&\sqrt{d(t)^2 + \mathcal {H}(d(t))^2} - \sqrt{u(t)^2 + \mathcal {H}(u(t))^2} \end{eqnarray} (5) which requires the computation of the Hilbert transform (|$\mathcal {H}(t)$|) and the amplitude of the analytical function corresponding to the original waveform. The envelope waveform difference is computed at discrete time steps selected in sampling the waveforms. E-WF gives more weight to the amplitude compared to the phase difference at the expense of losing some of the higher frequency content in the waveforms. For the E-WF, the expression for the Fréchet derivative matrix (|$\boldsymbol {G}$| in eq. 2) is: \begin{eqnarray} \boldsymbol {G} = \frac{-1}{\text{env}(u(t))} \left ( u(t) \frac{\partial u(t)}{\partial m} + \mathcal {H}(u(t)) \frac{\partial \mathcal {H}(u(t))}{\partial m} \right ) \end{eqnarray} (6) where, according to our hybrid approach, we employ NACT for the computation of the partial derivatives (∂u(t)/∂m, |$\partial \mathcal {H}(u(t))/\partial m$|) and SEM synthetics for the other functions (u(t), |$\mathcal {H}(u(t)), \text{env}(u(t))$|). 2.4 Model-space parametrization We represent the physical model, defined by the 3-D perturbations in three independent parameters (VS, ξ and |$Q^{-1}_\mu$|) with respect to the global average structure, in terms of cubic b-splines radially, (Mégnin & Romanowicz 2000) and spherical splines laterally (Wang & Dahlen 1995). In this study, since we seek to build an attenuation model starting from SEMUCB-WM1, the b-spline and spherical spline parametrizations are kept the same as in the starting model for elastic parameters (VS and ξ). For Qμ, we considered the same b-splines as used for the elastic parameters, while spherical splines were defined on progressively finer tesselations, as iterations progressed and the bandpass was increased. Fig. 1 shows the distribution of the spline nodes used in mapping Qμ in space. We used 20 b-spline nodes in the vertical direction, distributed from the shallowest crustal depth (30 km) in SEMUCB-WM1 to the core–mantle boundary with the b-spline node locations at radii: 3480, 3600, 3775, 4000, 4275, 4550, 4850, 5150, 5375, 5575, 5750, 5900, 6050, 6100, 6150, 6200, 6250, 6300, 6346, and 6361 km. The nodes have a denser distribution closer to the top, where we expect higher radial resolution. Figure 1. Open in new tabDownload slide Model space parametrization in radial and lateral directions. We employ 20 b-spline basis functions (Mégnin & Romanowicz 2000) in the radial direction from core–mantle boundary to the shallowest Moho depth of 30 km in SEMUCB-WM1 model (French & Romanowicz 2014). The targeted 11 b-spline basis functions in this study are coloured red. The lateral model parametrization mesh for Qμ is progressively refined with increasing passband at subsequent iterations. For elastic parameters (VS, ξ), we used the same lateral model parametrization as in SEMUCB-WM1. The mean lateral grid node spacing (Δ) and the corresponding passbands targeted with the lateral mesh are indicated on each panel. Figure 1. Open in new tabDownload slide Model space parametrization in radial and lateral directions. We employ 20 b-spline basis functions (Mégnin & Romanowicz 2000) in the radial direction from core–mantle boundary to the shallowest Moho depth of 30 km in SEMUCB-WM1 model (French & Romanowicz 2014). The targeted 11 b-spline basis functions in this study are coloured red. The lateral model parametrization mesh for Qμ is progressively refined with increasing passband at subsequent iterations. For elastic parameters (VS, ξ), we used the same lateral model parametrization as in SEMUCB-WM1. The mean lateral grid node spacing (Δ) and the corresponding passbands targeted with the lateral mesh are indicated on each panel. To image the upper-mantle attenuation structure, we inverted only the top b-splines (coloured red in the figure), progressively increasing the number from 6 (down to 250 km depth) to 11 (down to 800 km depth) in the final iteration. The four spherical spline node distributions presented for Qμ from top to bottom were used progressively as we increased the period range considered, as described in Section 2.5. The final Qμ model is defined with a mean distance between nodes of 4°. The grids used to define VS and ξ perturbations presented in the figure are kept intact throughout the iterations as they were used with SEMUCB-WM1. Thus, in the last iteration, we inverted for 112 662 VS, 7062 ξ and 28 182 |$Q^{-1}_\mu$| spherical spline coefficients to update the model. For attenuation model updates, we inverted for perturbations in |$\delta Q^{-1}_\mu$|. This is a natural choice, as this term appears in the imaginary part of both the complex elastic moduli (δln(μ)) and the corresponding complex wave speed perturbations (δln(VS)) (Karaoğlu & Romanowicz 2018). To test the robustness of this choice, in the early stages of this study, we also inverted for |$\delta \text{ln}\left(Q^{-1}_\mu \right)$| with negligible differences in the recovered models. The advantage of inverting for |$\delta \text{ln}\left(Q^{-1}_\mu \right)$| is the automatic avoidance of recovering negative attenuation (Adenis et al.2017a), however by inverting for |$\delta Q^{-1}_\mu$|, one can adjust the scalar regularization parameter α so as to avoid any negative attenuation, which can be considered as a way of keeping the inverse problem well posed. 2.5 Inversion strategy Our hybrid full-waveform inversion algorithm suffers from many of the well-known limitations of waveform inversion techniques that arise from the non-linearity of the inverse problem. This implies that we need to be careful with the selection of the starting model as well as expand the data set incrementally to include shorter periods. In attenuation imaging, we also need to pay attention to how we respectively update the elastic and attenuation models, as they are linked to one another, and we need to address the uncertainties in the source parameters, instrument gain factors, and near-source and near-receiver structures. 2.5.1 Starting model SEMUCB-WM1 is a 3-D radially anisotropic whole-mantle model with a fixed 1-D attenuation structure, which is a smoothed version of QL6 of Durek & Ekström (1996). In this study, the starting model has the elastic structure of SEMUCB-WM1, and a simpler 1-D attenuation structure, with Qμ = 300 in the crust, Qμ = 165 in the upper mantle and transition zone. Below 670 km, the 1-D Qμ profile is the same as QL6. The value of Qμ in the upper mantle is that of QL6 for the depth range of 220–670 km. Starting with this 1-D attenuation model gives us more flexibility for our attenuation inversions, as there is no a priori high-attenuation zone in the upper mantle, which allows us to check how well we can recover the well-known features in the 1-D global Qμ profile, giving us a point of reference for the quality of the final 3-D model. Fig. 8(a) shows the evolution of the global average Qμ model through our model iterations. 2.5.2 Extending the bandpass to shorter periods We started our inversions by filtering waveforms in the period range of 400–120 s and incrementally included shorter periods in subsequent iterations of the inversion for 3-D attenuation structure. We considered the following criteria for extending the bandpass: (i) the increase in the number of selected wavepackets after each SEM simulation run in the updated model, and (ii) the stability of the model from one iteration to the next. Whenever the number of selected wavepackets and model perturbations seemed to stabilize for the inversions carried out in a given period range, we decreased the cut-off period by 20 s until we reached the bandpass of 400–60 s, in the fifth iteration. 2.5.3 Stepwise anelastic model inversion In this study, we follow a stepwise inversion method in which we iteratively update the attenuation model (Qμ) and the elastic model (VS, ξ). As described above, we use (E-WF) and WF misfit definitions respectively for Qμ and velocity/anisotropy. In the course of this study, we carried out eight iterations, updating the elastic model twice (second and fifth iterations) and the attenuation model six times (Table 1). Additionally, we carried out one inversion for the source parameters of our event catalogue. Table 1. Details of iterations, with the period bandpass (cut-off and corner periods), inverted parameters, included data type (F: fundamental-mode surface wave, O: overtone surface wave, M: mixed (fundamental + overtone) surface wave) and the inverted b-spline nodes counted from the top . Bandpass Updated Data Inv. b-spline period parameters type Node count Iter-1 400–250–150–120 s Qμ F 6 Iter-2 400–250–120–100 s Qμ F 7 Iter-3 400–250–80–60 s VS, ξ F, M, O 11 Iter-4 400–250–100–80 s Qμ F 8 Iter-5 400–250–80–60 s Qμ F 8 Iter-6 400–250–80–60 s VS, ξ F, M, O 11 Source 400–250–80–60 s MT F, O NA inversion* lat, lon, dep Iter-7 400–250–80–60 s Qμ F, O 11 Iter-8 400–250–80–60 s Qμ F, O 11 Bandpass Updated Data Inv. b-spline period parameters type Node count Iter-1 400–250–150–120 s Qμ F 6 Iter-2 400–250–120–100 s Qμ F 7 Iter-3 400–250–80–60 s VS, ξ F, M, O 11 Iter-4 400–250–100–80 s Qμ F 8 Iter-5 400–250–80–60 s Qμ F 8 Iter-6 400–250–80–60 s VS, ξ F, M, O 11 Source 400–250–80–60 s MT F, O NA inversion* lat, lon, dep Iter-7 400–250–80–60 s Qμ F, O 11 Iter-8 400–250–80–60 s Qμ F, O 11 *Source inversion updates the moment tensor (MT) and coordinates. See Section 4 for further details. Open in new tab Table 1. Details of iterations, with the period bandpass (cut-off and corner periods), inverted parameters, included data type (F: fundamental-mode surface wave, O: overtone surface wave, M: mixed (fundamental + overtone) surface wave) and the inverted b-spline nodes counted from the top . Bandpass Updated Data Inv. b-spline period parameters type Node count Iter-1 400–250–150–120 s Qμ F 6 Iter-2 400–250–120–100 s Qμ F 7 Iter-3 400–250–80–60 s VS, ξ F, M, O 11 Iter-4 400–250–100–80 s Qμ F 8 Iter-5 400–250–80–60 s Qμ F 8 Iter-6 400–250–80–60 s VS, ξ F, M, O 11 Source 400–250–80–60 s MT F, O NA inversion* lat, lon, dep Iter-7 400–250–80–60 s Qμ F, O 11 Iter-8 400–250–80–60 s Qμ F, O 11 Bandpass Updated Data Inv. b-spline period parameters type Node count Iter-1 400–250–150–120 s Qμ F 6 Iter-2 400–250–120–100 s Qμ F 7 Iter-3 400–250–80–60 s VS, ξ F, M, O 11 Iter-4 400–250–100–80 s Qμ F 8 Iter-5 400–250–80–60 s Qμ F 8 Iter-6 400–250–80–60 s VS, ξ F, M, O 11 Source 400–250–80–60 s MT F, O NA inversion* lat, lon, dep Iter-7 400–250–80–60 s Qμ F, O 11 Iter-8 400–250–80–60 s Qμ F, O 11 *Source inversion updates the moment tensor (MT) and coordinates. See Section 4 for further details. Open in new tab We updated the elastic model using data filtered in the period range of 400–60 s as in the work of Lekić & Romanowicz (2011) and French et al. (2013). Following their approach, we have included the same Rayleigh and Love wave group velocity dispersion data for periods from 25 to 150 s (Ritzwoller, private communication, 2009, Shapiro & Ritzwoller 2002), in order to constrain the shallow upper-mantle elastic structure. Additionally, we built a smooth, anisotropic and homogenized crustal model, which was updated after each iteration to match the surface wave velocity dispersion measurements from 60 to 25 s (for details see Lekić & Romanowicz 2011; French & Romanowicz 2014). 2.5.4 Source and receiver terms We have tried to minimize the effect of uncertainties in the seismic source parameters, site effects and instrument responses in three ways: (i) through a careful data selection as described in Section 3, (ii) by introducing component-dependent receiver terms and event-dependent source terms to be included as additional unknown parameters and (iii) by updating the source parameters (moment tensors, MT and hypocentral coordinates) after the fifth iteration, as described in more detail in Section 4. The source and receiver terms that we have included in our attenuation inversions are scaling parameters that absorb the bulk of the amplitude anomalies that cannot be explained by the recovered attenuation models and updated source parameters. 2.6 Prior information Assuming that we have a sufficiently good starting model and data coverage, we seek to make the inverse problem well posed through the data and model covariance matrices (|$\boldsymbol {C_D}$|, |$\boldsymbol {C_M}$|) and the regularization scaling parameter (|$\boldsymbol {\alpha }$|) shown in eq. (1). For |$\boldsymbol {C_D}$|, we consider a diagonal matrix that weighs the data in accordance with their quality and redundancy (Li & Romanowicz 1996,appendix A). |$\boldsymbol {C_M}$|, on the other hand, is built on the basis of permissible correlation lengths of heterogeneity (Lekić & Romanowicz 2011). In building SEMUCB-WM1, the latter were fixed in the radial direction, with values ranging from 50 km at the shallowest upper-mantle depths to 300 km in the mid-mantle, and in the horizontal direction, they were spatially adjusted according to the values on the diagonal of the Gauss–Newton NACT-based-Hessian matrix (|$\boldsymbol {G^T C^{-1}_D G}$|). This is done in order to account for data coverage density and quality. The correlation lengths were allowed to vary between 400 and 1200 km for VS and 1200 and 3600 km for ξ, following French & Romanowicz (2014). Further details on how we built |$\boldsymbol {C_M}$| can be found in Lekić & Romanowicz (2011). In attenuation imaging, we kept the correlation lengths in the radial direction the same as in SEMUCB-WM1 throughout all our iterations. For the horizontal direction, on the other hand, we varied the range of correlation lengths according to the mean distance between nodes in the spherical spline mesh considered at different iterations. Based on our preliminary synthetic tests, we concluded that the minimum anomaly size we can recover is twice the mean nodal distance of the tessellation employed. This value we chose as the lower limit of the correlation lengths. As the upper limit we chose six times the mean nodal distance. For the four spherical spline meshes used for Qμ perturbations, presented in Fig. 1, we assigned the following values as the limits of the range of correlation lengths: (i) 6400–12 800 km for the mesh with Δ = 32°, (ii) 3200–9600 km for the mesh with Δ = 16°, (iii) 1600–4800 km for the mesh with Δ = 8° and (iv) 800–2400 km for the mesh with Δ = 4°. 3 DATA SET The global upper-mantle shear wave attenuation model presented in this study is based on long-period, time-windowed seismic waveforms filtered between 400 and 60 s period, and group velocity dispersion data in the period range 25-150s, the latter used only for updating elastic parameters VS (isotropic shear velocity) and ξ = (VSH/VSV)2 (radial anisotropy parameter). The three-component waveforms were extracted from the data collection used in the development of the elastic 3-D whole mantle model SEMUCB-WM1 (French & Romanowicz 2014), consisting of 273 events with magnitudes ranging from 5.8 to 7.3, recorded at up to 524 broad-band seismic stations (Fig. 2). Figure 2. Open in new tabDownload slide Distribution of the 273 seismic sources (circles) and 524 receivers (green triangles) used in building the anelastic model (SEMUCB-UMQ). The seismic sources are coloured as a function of their hypocentre depth which varies from 12 to 648.3 km. Figure 2. Open in new tabDownload slide Distribution of the 273 seismic sources (circles) and 524 receivers (green triangles) used in building the anelastic model (SEMUCB-UMQ). The seismic sources are coloured as a function of their hypocentre depth which varies from 12 to 648.3 km. Our waveform data set comprises fundamental mode, overtone and mixed (fundamental and overtone energy together) Love and Rayleigh waves, recorded on all three components of motion, for both minor and major arcs. While we use the complete waveform data set obtained by an automated individual waveform selection routine, as in the construction of SEMUCB-WM1, for updating the elastic model, for the attenuation inversion steps, as described in the previous sections, we followed stricter selection rules, resulting in the elimination of a significant part of the original data set. These additional rules include: Unlike for the development of SEMUCB-WM1, we used a sampling rate of 0.1 Hz, as opposed to the Nyquist rate for the minimum period targeted (33 mHz). Since in pre-processing steps (instrument response removal and filtering), the seismograms are not perfectly filtered, they are susceptible to frequency components greater than the Nyquist frequency. Increasing the sampling rate helps to avoid these unwanted aliasing effects and waveforms. Compared to the earlier practice, the choice of higher than Nyquist sampling rate eliminates many originally selected waveforms. We imposed a minimum waveform and envelope correlation coefficient between the selected wavepackets and synthetics of 0.5. This rule was imposed to make sure that the phase difference between the synthetic and observed waveforms match was good enough for extraction of the anelastic attenuation signal from the amplitudes. We eliminated any mixed phase waveforms that include both fundamental mode and overtone energy, in order to work with clear signals. This also helps eliminate those wavepackets that contain overlapping fundamental-mode and overtone surface waves traveling, respectively, along minor and major arcs. For each station and each individual component of motion, we required at least 10 fundamental-mode and 10 overtone surface wave records that satisfy the criteria above. We kept only those events for which at least 50 fundamental-mode and 50 overtone surface wave records were kept after imposing the criteria above. Applying these selection criteria results in the loss of around 50 per cent of the original data set. This is a price to pay for working with a more reliable data set and in particular, reducing the uncertainty level in the source and receiver terms that are also inverted for, as explained in the previous section. In the first four attenuation iterations, we reduced the data set even further by using only the fundamental-mode Rayleigh and Love waves that are sensitive primarily to the top 250–300 km of the mantle. Including Love waves help better constrain the shallow upper-mantle attenuation structure. Next, we inverted for the source parameters (MT and hypocentre coordinates) of 193 events used in the last three attenuation model inversions. We excluded 80 events that failed to provide data satisfying our strict selection criteria for attenuation inversion. In the last two attenuation model iterations, we included overtones in the inversions as well. By doing so, we tried to extend our model to larger depths. The inversion steps are presented in a flowchart in Fig. 3 and discussed further in Section 5.1. Figure 3. Open in new tabDownload slide An overview of the stepwise inversion steps followed to build SEMUCB-UMQ model. Figure 3. Open in new tabDownload slide An overview of the stepwise inversion steps followed to build SEMUCB-UMQ model. Table 1 provides the details of the data set and bandpass used at each iteration. The data and synthetics are filtered with a cosine-taper window for varying cut-off and corner frequencies in accordance with our strategy to include progressively shorter periods into our inversions. The data set of French & Romanowicz (2014) comprised 168 986 fundamental and overtone surface wavepackets along with mixed surface and body wavepackets in the final iteration of the SEMUCB-WM1 construction. Once we imposed the additional selection criteria for attenuation imaging on that data set, the number of wavepackets dropped to 87 077, all three components included. With the same selection criteria imposed on the synthetics computed with SEMUCB-UMQ, 92 290 wavepackets were selected. In other words, SEMUCB-UMQ provides improved fits to a data set almost 6 per cent larger than the one that SEMUCB-WM1 can explain under the same selection criteria. Note that we changed the 1-D Qμ profile of SEMUCB-WM1 at the beginning, which partly explains the significant drop in wavepackets selected at the beginning of our inversion process. We consider the increase in the number of wavepackets accepted between the starting and final models as an indication that the model has evolved in the right direction. Further details on these figures are provided in Section 5.1. We note that major-arc phases make a significant contribution to the global ray density coverage we achieve with our data set. As shown in Fig. 4 for our final data set, the ray coverage in the Southern Hemisphere is mostly provided by major-arc phases. Still, major-arc phases come with a disadvantage for attenuation imaging. As they travel significantly longer distances, they become more prone to amplitude perturbations due to focusing from elastic heterogeneities, therefore, it is somewhat more difficult to extract a reliable attenuation signal from them. Figure 4. Open in new tabDownload slide Ray density in 5 × 5 degree cells for the data types considered, grouped with respect to the travel path (minor and major arcs) and receiver components. The maximum and minimum hit counts per cell are presented under the colour bars. Figure 4. Open in new tabDownload slide Ray density in 5 × 5 degree cells for the data types considered, grouped with respect to the travel path (minor and major arcs) and receiver components. The maximum and minimum hit counts per cell are presented under the colour bars. 4 SOURCE INVERSION In building SEMum, SEMum2 and SEMUCB-WM1, the MT and source locations were kept fixed at those provided by the Harvard Centroid Moment Tensor project (www.globalcmt.org). French & Romanowicz (2014) inverted for the source parameters of 20 events using surface wave and body wave data types to assess the effect of perturbing the CMT solutions on the variance reduction and concluded that source parameters can have a significant effect on the amplitude misfits especially for shorter period body waves (cut-off period of 32 s). For surface waves, on the other hand, the variations were not significant, at least for the subset of events they analysed. The methodology they used is similar to that of Liu et al. (2004), which is based on computing SEM-based 3-D Green’s functions and location partial derivatives with a 3-D heterogeneous model. The data set was built in a similar way to the structural inversion with the selected time-windowed waveforms, which were weighted according to their residual variance and redundancy with the latter defined on the basis of the source–receiver azimuth. Finally, the inversions were performed through a Gauss–Newton optimization scheme in which they observed that only the depth parameter needed to be damped to keep the inversion problem well posed. In the inversion step, the perturbations to the Global CMT solutions were enforced to have no isotropic component. The main obstacle with the source inversion methodology of French & Romanowicz (2014) is its computational cost. This approach requires nine simulations (six for MT components and three for hypocentre coordinates) per event. To overcome this obstacle, in this study, we computed the Green’s functions and location partial derivatives using NACT, similarly to our inversion for structural parameters. This reduces the cost significantly for building the Fréchet derivative matrix. To test this approach, we conducted a set of source inversions only for the focal mechanism keeping the hypocentre coordinates fixed as provided by the Global CMT solution. Inverting only for the MT is a linear problem, therefore, in this test, we did not need any starting model (zero MT). The goal was to see whether we can recover the MT using our data set consisting of fundamental and overtone waveforms filtered in the 400–60 s passband. Fig. 5 shows our solutions in comparison to the corresponding Global CMT solutions for six events at different depths. The details of the solutions for these events are presented in Table 2. Figure 5. Open in new tabDownload slide Seismic source inversion tests using NACT-based approach for six sample events. The geographic distribution of the events is shown with the corresponding focal mechanism beach-ball diagrams both for Global CMT (GCMT) and our solutions (NEW). The coordinates of the events are fixed to the values from the GCMT solution and we only invert for the MT. Further details are given in Table 2. Figure 5. Open in new tabDownload slide Seismic source inversion tests using NACT-based approach for six sample events. The geographic distribution of the events is shown with the corresponding focal mechanism beach-ball diagrams both for Global CMT (GCMT) and our solutions (NEW). The coordinates of the events are fixed to the values from the GCMT solution and we only invert for the MT. Further details are given in Table 2. Table 2. Comparison of the Global CMT and inverted source models for the six events presented in Fig. 5. Event ID Lat(°) Lon(°) Depth (km) GCMT NEW Berkeley conversion M0 (×1019 Nm) Mw M0 (×1019 Nm) Mw C032098D −50.06 162.89 15.0 1.31 6.68 1.26 6.67 C011197D 18.34 −102.58 40.0 6.09 7.13 4.70 7.05 C111202B −56.49 −26.89 116.3 0.25 6.21 0.29 6.24 C061603C 55.48 160.25 180.9 2.46 6.86 2.54 6.87 C031696C 29.12 139.12 477.9 1.06 6.62 1.07 6.62 C080700E −6.95 123.53 648.3 0.60 6.46 0.70 6.50 Event ID Lat(°) Lon(°) Depth (km) GCMT NEW Berkeley conversion M0 (×1019 Nm) Mw M0 (×1019 Nm) Mw C032098D −50.06 162.89 15.0 1.31 6.68 1.26 6.67 C011197D 18.34 −102.58 40.0 6.09 7.13 4.70 7.05 C111202B −56.49 −26.89 116.3 0.25 6.21 0.29 6.24 C061603C 55.48 160.25 180.9 2.46 6.86 2.54 6.87 C031696C 29.12 139.12 477.9 1.06 6.62 1.07 6.62 C080700E −6.95 123.53 648.3 0.60 6.46 0.70 6.50 Open in new tab Table 2. Comparison of the Global CMT and inverted source models for the six events presented in Fig. 5. Event ID Lat(°) Lon(°) Depth (km) GCMT NEW Berkeley conversion M0 (×1019 Nm) Mw M0 (×1019 Nm) Mw C032098D −50.06 162.89 15.0 1.31 6.68 1.26 6.67 C011197D 18.34 −102.58 40.0 6.09 7.13 4.70 7.05 C111202B −56.49 −26.89 116.3 0.25 6.21 0.29 6.24 C061603C 55.48 160.25 180.9 2.46 6.86 2.54 6.87 C031696C 29.12 139.12 477.9 1.06 6.62 1.07 6.62 C080700E −6.95 123.53 648.3 0.60 6.46 0.70 6.50 Event ID Lat(°) Lon(°) Depth (km) GCMT NEW Berkeley conversion M0 (×1019 Nm) Mw M0 (×1019 Nm) Mw C032098D −50.06 162.89 15.0 1.31 6.68 1.26 6.67 C011197D 18.34 −102.58 40.0 6.09 7.13 4.70 7.05 C111202B −56.49 −26.89 116.3 0.25 6.21 0.29 6.24 C061603C 55.48 160.25 180.9 2.46 6.86 2.54 6.87 C031696C 29.12 139.12 477.9 1.06 6.62 1.07 6.62 C080700E −6.95 123.53 648.3 0.60 6.46 0.70 6.50 Open in new tab Our preliminary tests on the non-linear source parameter inversion (MT and hypocentre coordinates), using the NACT-based approach showed negligible contributions in the bandpass 400–60 s using only the fundamental-mode surface waves. Therefore, we updated the Global CMT source parameters at the end of the sixth iteration for 193 events that provided data satisfying our strict selection criteria. The inversions were carried out using minor- and major-arc fundamental and overtone surface waves. In the inversions, the source depth parameter was damped more for the shallow events (<50 km) and less for the deep events following the observation we made in our preliminary tests. Similar to the work of French & Romanowicz (2014), our inversions led to negligible variations in the latitude and longitude of the sources. Source depths, on the other hand, show variations that vary with peak-to-peak perturbation values of −6.73 km to +2.60 km. Scalar moment magnitudes M0 also vary significantly with peak-to-peak values of |$-41\hbox{ per cent}$| to +26 per cent with respect to the Global CMT solutions. Fig. 6 shows M0 and depth perturbations we recovered both in histograms and in map view with their geographical distributions. Figure 6. Open in new tabDownload slide The scalar moment magnitude gain/loss (left) and depth perturbations (right) for the 193 inverted earthquake sources, with respect to the Global CMT solution. The magnitude and depth perturbations are presented in map form and in the form of histograms. Figure 6. Open in new tabDownload slide The scalar moment magnitude gain/loss (left) and depth perturbations (right) for the 193 inverted earthquake sources, with respect to the Global CMT solution. The magnitude and depth perturbations are presented in map form and in the form of histograms. The map distribution of M0 perturbations shows a visual correlation with the geographic settings of the events. In general, we see magnitude gain along the ridges and trenches and loss in the continental areas with some exceptions. This observation is in line with the study of Hjörleifsdóttir & Ekström (2010), who reported that Global CMT solutions might have a tendency to have higher M0 values in the areas with thick crust, like Tibet, and lower M0 values along oceanic ridges. They reported that for most of the 50 events simulated, the scalar moment magnitude perturbations were less than 10 per cent. In the present case, M0 variations have a mean of |$-2.8\hbox{ per cent}$| with a standard deviation of 10.8 per cent. As for depth perturbations, we observe larger values for deeper events which might be partially due to the damping we imposed. However, there is no clear correlation between the source depth and the depth perturbations as shown in Fig. S1 in the Supporting Information. 5 RESULTS 5.1 Inversion and model fits Our iterative inversion process is summarized in Fig. 3. The process involves a total of eight iterations, which include inversions for perturbations in |$Q^{-1}_{\mu }$| alternating with updates of the elastic structure (VS, ξ), and one iteration to update the source parameters. Starting with data filtered between 400 and 120 s, the bandpass is progressively increased to shorter periods (see Table 1). The starting model comprises the 3-D elastic structure of SEMUCB-WM1 and a simple 1-D Qμ model with constant value throughout the upper mantle (see Fig. 8). For the first iteration, in which only fundamental-mode three-component wavepackets were included, the number of fundamental-mode surface waveforms that satisfied our selection criteria was 12 595 of which 15 per cent was Love waves, and the total number of events was 123. We first carried out inversions to recover the long-wavelength structure in attenuation in the upper mantle, parametrizing the model laterally with 42 spherical splines (mean distance between nodes of 32°), and radially with six b-splines that extended from the Moho to ∼300 km depth. As we retrieved the long-wavelength structure, we refined our model space in the lateral direction and increased the number of b-spline nodes included in the inversion until we reached 2562 spherical spline nodes laterally (mean distance between nodes of 4°), and 11 b-spline nodes vertically, extending the model space to ∼800 km depth. At each iteration, our data set expanded, with additional waveforms satisfying our selection criteria. The fact that the data set increased at each iteration was a strong indication that the inversion process was working and the model improving. In the first four attenuation model inversions, we only included fundamental-mode surface waves, with the objective of retrieving the shallow upper-mantle attenuation structure before trying to map deeper attenuation heterogeneities. As long wavelength attenuation structure developed, we updated the elastic model (VS, ξ) twice. For the elastic model updates, we used the full long-period surface wave data set (period range 400–60 s) at our disposal from the construction of SEMUCB-WM1, with the same less strict wavepacket selection criteria as in French & Romanowicz (2014). However, different from the construction of SEMUCB-WM1, we did not include any shorter period body waveforms here. After updating the elastic model, a second time in the sixth iteration, we included overtone surface waveforms that satisfied our selection criteria, and performed one inversion for updates to the source parameters, following the technique explained in Section 4 for the 193 events that by then had accepted data according to our attenuation inversion criteria. In the last two iterations, we updated the attenuation model, using the new source parameters and including overtone surface waveforms in our data set to improve constraints on attenuation structure in the deeper part of our model. The conventional way of assessing the improvement in the fit to data, is to compare the variance reduction achieved with models obtained at successive iterations. Because the data set increased at each iteration, for this computation, we computed the variance reduction in the starting model and subsequent iterations using the data set obtained in the final iteration. Table 3 shows that the improvement in variance reduction from starting model to SEMUCB-UMQ is 32 per cent (from 55.4 per cent to 77.3 per cent), of which the 5 per cent comes from the source inversion. For the data set selected based on the starting model, the variance reduction increases by 21.4 per cent (from 71.4 per cent to 92.8 per cent). Table 3. Variance reduction (|$\text{VR} = 100x\left[1-||\mathbf {d-g(m)}||_2^2/||\mathbf {d}||_2^2\right]$|) calculated for the waveform windows selected in the last iteration for different models targeting the passband of 400–60 s, organized in receiver component and data type. Passband SEMUCB Start Iter-6 Iter-6 Iter-7 SEMUCB-UMQ 400–60 s (sr inv) (+overtones) (Iter-8) L (per cent) Fund 68.7 30.6 61.1 68.7 73.4 74.6 Over 79.0 73.5 81.6 82.1 82.7 83.8 T (per cent) Fund 77.1 31.1 72.8 75.5 76.4 76.1 Over 70.1 70.0 81.2 81.6 81.9 83.3 Z (per cent) Fund 71.6 31.1 61.4 69.5 74.7 76.7 Over 78.8 68.4 81.2 82.4 82.6 83.2 Total (per cent) 74.0 55.4 68.5 73.4 76.4 77.3 Passband SEMUCB Start Iter-6 Iter-6 Iter-7 SEMUCB-UMQ 400–60 s (sr inv) (+overtones) (Iter-8) L (per cent) Fund 68.7 30.6 61.1 68.7 73.4 74.6 Over 79.0 73.5 81.6 82.1 82.7 83.8 T (per cent) Fund 77.1 31.1 72.8 75.5 76.4 76.1 Over 70.1 70.0 81.2 81.6 81.9 83.3 Z (per cent) Fund 71.6 31.1 61.4 69.5 74.7 76.7 Over 78.8 68.4 81.2 82.4 82.6 83.2 Total (per cent) 74.0 55.4 68.5 73.4 76.4 77.3 Open in new tab Table 3. Variance reduction (|$\text{VR} = 100x\left[1-||\mathbf {d-g(m)}||_2^2/||\mathbf {d}||_2^2\right]$|) calculated for the waveform windows selected in the last iteration for different models targeting the passband of 400–60 s, organized in receiver component and data type. Passband SEMUCB Start Iter-6 Iter-6 Iter-7 SEMUCB-UMQ 400–60 s (sr inv) (+overtones) (Iter-8) L (per cent) Fund 68.7 30.6 61.1 68.7 73.4 74.6 Over 79.0 73.5 81.6 82.1 82.7 83.8 T (per cent) Fund 77.1 31.1 72.8 75.5 76.4 76.1 Over 70.1 70.0 81.2 81.6 81.9 83.3 Z (per cent) Fund 71.6 31.1 61.4 69.5 74.7 76.7 Over 78.8 68.4 81.2 82.4 82.6 83.2 Total (per cent) 74.0 55.4 68.5 73.4 76.4 77.3 Passband SEMUCB Start Iter-6 Iter-6 Iter-7 SEMUCB-UMQ 400–60 s (sr inv) (+overtones) (Iter-8) L (per cent) Fund 68.7 30.6 61.1 68.7 73.4 74.6 Over 79.0 73.5 81.6 82.1 82.7 83.8 T (per cent) Fund 77.1 31.1 72.8 75.5 76.4 76.1 Over 70.1 70.0 81.2 81.6 81.9 83.3 Z (per cent) Fund 71.6 31.1 61.4 69.5 74.7 76.7 Over 78.8 68.4 81.2 82.4 82.6 83.2 Total (per cent) 74.0 55.4 68.5 73.4 76.4 77.3 Open in new tab To better analyse the contribution of the 3-D Qμ structure to the improvement achieved, we simulated 10 randomly chosen earthquakes selected from our event catalogue in the SEMUCB-UMQ model, but including only the 1-D part of the Qμ model. For the selected 10 events, the variance reduction dropped by 4.2  per cent in total, compared to that in the 3-D model, with 7.1  per cent for fundamental-mode surface waves (L: 6.2 per cent, T: 7.7  per cent and Z: 6.4  per cent) and 2.9  per cent for overtone surface waves (L: 2.0 per cent, T: 3.0 per cent and Z: 3.2  per cent). This indicates that a large fraction of the variance reduction as computed above comes from improving the 1-D Q model. Fig. 7 shows the waveform and envelope fit improvement for some fundamental-mode minor- and major-arc surface waves recorded at two different receivers for one event. Figure 7. Open in new tabDownload slide Some examples of waveform fit improvement for different receiver components and seismic phases through attenuation imaging iterations. Compared are the synthetics computed for the starting model (blue), the recovered attenuation model using only the fundamental-mode data type at the sixth iteration (grey), after source inversion (green) and SEMUCB-UMQ. Waveforms are shown only for data and the synthetics computed with SEMUCB-UMQ for clarity. The bold dashed lines indicate the window which our automatic picking procedure selected. Figure 7. Open in new tabDownload slide Some examples of waveform fit improvement for different receiver components and seismic phases through attenuation imaging iterations. Compared are the synthetics computed for the starting model (blue), the recovered attenuation model using only the fundamental-mode data type at the sixth iteration (grey), after source inversion (green) and SEMUCB-UMQ. Waveforms are shown only for data and the synthetics computed with SEMUCB-UMQ for clarity. The bold dashed lines indicate the window which our automatic picking procedure selected. Another way of assessing the model improvement achieved through our iterative inversion process is to compare the number of wavepackets that satisfy our strict data selection criteria for attenuation imaging as a function of iteration. However, such a comparison is only valid if it is done within the same passband. To that end, we chose the period range of 400–60 s and compared the sizes of the corresponding data sets at different iterations. To compare the performances of SEMUCB-UMQ with SEMUCB-WM1 in terms of the number of wavepackets, we also imposed the strict data selection criteria for attenuation imaging on the data set used by French & Romanowicz (2014) in the same frequency band. The number of selected waveforms grouped by receiver component and data type are shown in Table 4. Note that, for the ‘Starting Model’, these numbers do not correspond to the data set used for the inversion, since initially we targeted the period range of 400–120 s. Comparing the starting model and SEMUCB-UMQ within the same period range (400–60 s), the increase we achieved in the number of selected wavepackets is around 80 per cent. The expansion in the data set comes mostly from fundamental-mode surface waveforms of which the number has increased by 360 per cent. The number of overtone surface waveforms, on the other hand, has increased by 22 per cent. The contribution of the seismic source updates to the number of selected wavepackets is less than 5 per cent. This is an indication of the robustness of the selection criteria which eliminates any earthquake sources and receivers with high uncertainty. A comparison between SEMUCB-UMQ and SEMUCB-WM1 shows an increase of 6 per cent in the size of the data set accepted by our ‘picking’ procedure (e.g. Table 4). Table 4. Waveform-window counts (Nwp) and number of data points (Nd), sampled at 0.1 Hz, satisfying our careful selection criteria for the passband of 400–60 s, organized in receiver component and data type. Bandpass SEMUCB-WM1 Start Iter-6 Iter-6 (sr inv) SEMUCB-UMQ (400–60 s) Nwp Nd Nwp Nd Nwp Nd Nwp Nd Nwp Nd L Fund 8 004 1 222 740 1780 257 000 7 207 1 082 245 8 154 1 243 657 8 826 1 349 951 Over 15 362 2 780 957 13 681 2 420 171 15 348 2 777 232 15 549 2 818 218 15 734 2 855 907 T Fund 13 724 1 511 338 3226 323 552 13 297 1 457 521 13 160 1 444 783 15 011 1 658 381 Over 11 427 1 180 053 10 540 1 092 688 11 581 1 194 996 12 035 1 240 287 12 307 1 267 326 Z Fund 15 394 2 389 901 3801 549 259 13 578 2 073 394 15 572 2 414 870 16 444 2 556 684 Over 23 166 4 010 440 18 498 3 043 013 23 032 3 978 762 23 664 4 106 618 23 968 4 165 485 Total 87 077 13 095 429 51 529 7 685 683 84 043 12 564 150 88 134 13 268 433 92 290 13 853 734 Bandpass SEMUCB-WM1 Start Iter-6 Iter-6 (sr inv) SEMUCB-UMQ (400–60 s) Nwp Nd Nwp Nd Nwp Nd Nwp Nd Nwp Nd L Fund 8 004 1 222 740 1780 257 000 7 207 1 082 245 8 154 1 243 657 8 826 1 349 951 Over 15 362 2 780 957 13 681 2 420 171 15 348 2 777 232 15 549 2 818 218 15 734 2 855 907 T Fund 13 724 1 511 338 3226 323 552 13 297 1 457 521 13 160 1 444 783 15 011 1 658 381 Over 11 427 1 180 053 10 540 1 092 688 11 581 1 194 996 12 035 1 240 287 12 307 1 267 326 Z Fund 15 394 2 389 901 3801 549 259 13 578 2 073 394 15 572 2 414 870 16 444 2 556 684 Over 23 166 4 010 440 18 498 3 043 013 23 032 3 978 762 23 664 4 106 618 23 968 4 165 485 Total 87 077 13 095 429 51 529 7 685 683 84 043 12 564 150 88 134 13 268 433 92 290 13 853 734 Open in new tab Table 4. Waveform-window counts (Nwp) and number of data points (Nd), sampled at 0.1 Hz, satisfying our careful selection criteria for the passband of 400–60 s, organized in receiver component and data type. Bandpass SEMUCB-WM1 Start Iter-6 Iter-6 (sr inv) SEMUCB-UMQ (400–60 s) Nwp Nd Nwp Nd Nwp Nd Nwp Nd Nwp Nd L Fund 8 004 1 222 740 1780 257 000 7 207 1 082 245 8 154 1 243 657 8 826 1 349 951 Over 15 362 2 780 957 13 681 2 420 171 15 348 2 777 232 15 549 2 818 218 15 734 2 855 907 T Fund 13 724 1 511 338 3226 323 552 13 297 1 457 521 13 160 1 444 783 15 011 1 658 381 Over 11 427 1 180 053 10 540 1 092 688 11 581 1 194 996 12 035 1 240 287 12 307 1 267 326 Z Fund 15 394 2 389 901 3801 549 259 13 578 2 073 394 15 572 2 414 870 16 444 2 556 684 Over 23 166 4 010 440 18 498 3 043 013 23 032 3 978 762 23 664 4 106 618 23 968 4 165 485 Total 87 077 13 095 429 51 529 7 685 683 84 043 12 564 150 88 134 13 268 433 92 290 13 853 734 Bandpass SEMUCB-WM1 Start Iter-6 Iter-6 (sr inv) SEMUCB-UMQ (400–60 s) Nwp Nd Nwp Nd Nwp Nd Nwp Nd Nwp Nd L Fund 8 004 1 222 740 1780 257 000 7 207 1 082 245 8 154 1 243 657 8 826 1 349 951 Over 15 362 2 780 957 13 681 2 420 171 15 348 2 777 232 15 549 2 818 218 15 734 2 855 907 T Fund 13 724 1 511 338 3226 323 552 13 297 1 457 521 13 160 1 444 783 15 011 1 658 381 Over 11 427 1 180 053 10 540 1 092 688 11 581 1 194 996 12 035 1 240 287 12 307 1 267 326 Z Fund 15 394 2 389 901 3801 549 259 13 578 2 073 394 15 572 2 414 870 16 444 2 556 684 Over 23 166 4 010 440 18 498 3 043 013 23 032 3 978 762 23 664 4 106 618 23 968 4 165 485 Total 87 077 13 095 429 51 529 7 685 683 84 043 12 564 150 88 134 13 268 433 92 290 13 853 734 Open in new tab The comparisons of the data sets at different iterations and with SEMUCB-WM1 within the period range of 400–60 s indicate: (i) the significance of the high-attenuation zone in the upper mantle as previously reported (e.g Gilbert & Dziewonski 1975; Dziewonski & Anderson 1981; Anderson & Hough 1984), (ii) stronger heterogeneity in the shallow upper mantle and (iii) improvements still needed to source parameters and elastic structure in order to exploit better the relatively low sensitivity of overtones to upper-mantle attenuation structure. 5.2 Upper-mantle attenuation structure 5.2.1 Global average Qμ profile Fig. 8(a) shows the evolution of the 1-D global average depth profile of our Qμ model in successive iterations of the inversion. While we did not include any high-attenuation zone in the uppermost mantle in the starting model, such a zone emerges in the 80–220 km depth range with an attenuation maximum at 129 km depth, where Qμ reaches a final value of 71.7. The depth range of the high-attenuation zone is comparable to that of previous studies (Fig. 8b), although there is some variability, due in particular to imposed layering in some models (e.g. Gung & Romanowicz 2004; Lawrence & Wysession 2006; Durek & Ekström 1996). The value of maximum in attenuation is comparable to QRFSI12 (Dalton et al.2008), while its depth is somewhat shallower (129 km compared to 165 km), in better agreement with that of the 1-D Qμ profile chosen in the development of SEMUCB-WM1 (French & Romanowicz 2014). The latter has a more pronounced and narrower maximum, due to the way the QL6 model was smoothed in its construction. The parametrization, in particular layering, also affects the recovery from the high-attenuation zone in the depth range 200–600 km. Our final model SEMUCB-UMQ exhibits a smooth profile with a slight minimum in attenuation around 300 km depth, lower attenuation than other models between 270 and 400 km depth, and slightly higher attenuation than most other models in the transition zone, except for QRLW8, which also included overtones in its construction. Figure 8. Open in new tabDownload slide (a) The evolution of the global average Qμ structure though iterations. (b) Comparison of different zero-degree Qμ models namely, QL6 of Durek & Ekström (1996), QLM9 of Lawrence & Wysession (2006), 1-D Qμ model of SEMUCB-WM1 by French & Romanowicz (2014) and the zero-degree structures of QRLW8 of Gung & Romanowicz (2004) and QRFSI12 of Dalton et al. (2008). (c) Normalized root mean square of the |$Q^{-1}_\mu$| perturbations in the compared 3-D heterogeneous models (QRLW8, QRFSI12 and SEMUCB-UMQ). The peak values are shown in the legend. Figure 8. Open in new tabDownload slide (a) The evolution of the global average Qμ structure though iterations. (b) Comparison of different zero-degree Qμ models namely, QL6 of Durek & Ekström (1996), QLM9 of Lawrence & Wysession (2006), 1-D Qμ model of SEMUCB-WM1 by French & Romanowicz (2014) and the zero-degree structures of QRLW8 of Gung & Romanowicz (2004) and QRFSI12 of Dalton et al. (2008). (c) Normalized root mean square of the |$Q^{-1}_\mu$| perturbations in the compared 3-D heterogeneous models (QRLW8, QRFSI12 and SEMUCB-UMQ). The peak values are shown in the legend. Since inverting for the 1-D shear attenuation profile from a starting constant Qμ model is the most non-linear part of our inversion process, the fact that the final model recovers the well-known main features of the upper-mantle profile encourages us to have confidence in the validity of our approach. Table 5 gives the values of the average 1-D Qμ structure of SEMUCB-UMQ. Table 5. Zero-degree Qμ profile of SEMUCB-UMQ model. Depth Qμ (km) 20 300.0 50 263.3 100 167.8 150 78.4 200 74.6 250 96.2 300 150.1 350 210.8 400 212.7 450 203.3 500 192.0 600 179.0 Depth Qμ (km) 20 300.0 50 263.3 100 167.8 150 78.4 200 74.6 250 96.2 300 150.1 350 210.8 400 212.7 450 203.3 500 192.0 600 179.0 Open in new tab Table 5. Zero-degree Qμ profile of SEMUCB-UMQ model. Depth Qμ (km) 20 300.0 50 263.3 100 167.8 150 78.4 200 74.6 250 96.2 300 150.1 350 210.8 400 212.7 450 203.3 500 192.0 600 179.0 Depth Qμ (km) 20 300.0 50 263.3 100 167.8 150 78.4 200 74.6 250 96.2 300 150.1 350 210.8 400 212.7 450 203.3 500 192.0 600 179.0 Open in new tab 5.2.2 3-D Qμ variations In this section, we describe the 3-D features of SEMUCB-UMQ and compare them to two other global upper-mantle attenuation (Qμ) models namely, QRLW8 of Gung & Romanowicz (2004) and QRFSI12 of Dalton et al. (2008). Model QRLW8 was also built from three-component fundamental-mode and overtone minor- and major-arc surface waveforms in the period range 400–60s. The differences in the construction of these two models are: (i) the starting elastic and 1-D Qμ models, (ii) the methods used for seismic wavefield simulations (PAVA for QRLW8 and SEM for SEMUCB-UMQ) and Hessian matrix computations (PAVA for QRLW8 and NACT for SEMUCB-UMQ) and (iii) the set of events and stations used to assemble the data set. Additionally, QRLW8 was built with a model parametrization laterally in spherical harmonics up to degree 8, and no correction for focusing was included. QRFSI12, on the other hand, was built from Z-component fundamental-mode Rayleigh wave records in the period range 250–50 s. Frequency-dependent focusing as well as source and station terms were considered in the inversion, and the model is parametrized laterally in spherical harmonics up to degree 12. Fig. 9 compares maps of lateral variations in Qμ at different upper-mantle depths for the three models. Down to 200-250 km depth, all three models show generally good correlation of the lateral variations in Qμ and tectonic provinces. At shallow depths, down to ∼100 km, the ridge signal dominates, most clearly in SEMUCB-UMQ and in QRFSI12, with clear low Qμ zones concentrated along the ridge system both in the Pacific, Indian and Atlantic oceans, somewhat more pronounced in the north Atlantic in SEMUCB-UMQ than in QRFSI12, while it is slightly less continuous in the former along the southeastern Indian ridge. Backarc regions are also prominent (high attenuation) in the Tonga Fiji area in both models, but in the Alaska/Aleutian area, this is less pronounced in SEMUCB-UMQ than in QRFSI12. A continuous band of slightly higher attenuation, probably associated with the Rea Sea/East Africa rift, extends into northern Africa, in all three models. Figure 9. Open in new tabDownload slide Map view of |$Q^{-1}_\mu$| variations with respect to the mean values at a range of depths for three global attenuation models built with surface wave data sets, namely QRLW8 of Gung & Romanowicz (2004), QRFSI12 of Dalton et al. (2008) and the model built in this study (SEMUCB-UMQ). The peak-to-peak variations are presented at the bottom of each panel both as the perturbation value (|$\delta Q^{-1}_\mu {\times} 10^3$|) and in percentage. The values in the parentheses shows the mean |$Q^{-1}_\mu {\times} 10^3$| values at the corresponding depths. The green circles denote the hotspot locations reported by Steinberger (2000). Figure 9. Open in new tabDownload slide Map view of |$Q^{-1}_\mu$| variations with respect to the mean values at a range of depths for three global attenuation models built with surface wave data sets, namely QRLW8 of Gung & Romanowicz (2004), QRFSI12 of Dalton et al. (2008) and the model built in this study (SEMUCB-UMQ). The peak-to-peak variations are presented at the bottom of each panel both as the perturbation value (|$\delta Q^{-1}_\mu {\times} 10^3$|) and in percentage. The values in the parentheses shows the mean |$Q^{-1}_\mu {\times} 10^3$| values at the corresponding depths. The green circles denote the hotspot locations reported by Steinberger (2000). Down to 200 km depth, the low-attenuation cratons stand out particularly well in the Canadian shield, Siberia, Australia, Antarctica and Africa, in both QRFSI12 and our new model SEMUCB-UMQ. In fact, in SEMUCB-UMQ, prominent high Qμ regions are clearly restricted to old continental areas in the depth range 150–200 km. Below 250 km, |$Q^{-1}_\mu$| perturbations become weaker (see change in colour scale) particularly so for models SEMUCB-UMQ and QRLW8 (which include constraints from overtones), and the pattern of heterogeneity progressively changes. We observe higher-than-average attenuation primarily in the southwestern Pacific and around the East African Rift extending towards south Africa. In SEMUCB-UMQ, the craton signal disappears, except perhaps in Australia, and is replaced by a low |$Q^{-1}_\mu$| band that extends from northern Siberia to the southern Atlantic. The long wavelength features, with high attenuation in the southwestern Pacific and low attenuation from northern Eurasia to the south Atlantic are present in all three models. A comparison of the root mean square (rms) of the |$Q^{-1}_\mu$| perturbations for the three models (Fig. 8c) shows that QRFSI12 has the strongest variations among the three models with a peak rms value of 3.90 in the high |$Q^{-1}_\mu$| zone, compared to 3.11 for SEMUCB-UMQ and 2.92 for QRLW8. The depth variation of rms values shows a similar trend for QRFSI12 and SEMUCB-UMQ, with the strongest heterogeneity located in the depth range 50–250 km. In SEMUCB-UMQ, the disappearance of the cratonic keels at larger depths seems to be a contribution of the overtone surface waves that we included in our data set in the last two iterations. Fig. 10 shows the distribution of |$Q^{-1}_\mu$| heterogeneity recovered below 250 km before and after adding the overtone surface waves in the inversions. Without the contribution of overtones, the high attenuation along the oceanic ridges and low attenuation in the shields and cratons seems to extend down to 450 km, manifesting poor depth resolution. Figure 10. Open in new tabDownload slide Map views of the global |$Q^{-1}_\mu$| variations in the recovered models at the end of the sixth iteration, which was built using only the fundamental-mode data type (R1, G1, R2, G2), and the final iteration, that was built with the inclusion of the overtone data type (XR1, XG1, XR2, XG2). The differences between the two models above 250 km are indistinguishable. Figure 10. Open in new tabDownload slide Map views of the global |$Q^{-1}_\mu$| variations in the recovered models at the end of the sixth iteration, which was built using only the fundamental-mode data type (R1, G1, R2, G2), and the final iteration, that was built with the inclusion of the overtone data type (XR1, XG1, XR2, XG2). The differences between the two models above 250 km are indistinguishable. 5.2.3 Spectral analysis and correlations Fig. 11 shows the spectrum of attenuation heterogeneity at different depths for the three models presented in Fig. 9. In contrast to the other two models, SEMUCB-UMQ shows strong power up to angular degree 16 at shallow upper-mantle depths, while the power becomes negligible below 450 km, even though we did not impose any stronger damping in the depth range extending down to 800 km. In all the attenuation models, the well-known degree-2 signal seen in elastic models is strong, reaching its peak between 150–200 km for both SEMUCB-UMQ and QRFSI12, and above 100 km depth for QRLW8. The signature of the continent–ocean distribution at degree 5 (e.g Meschede & Romanowicz 2015) shows up clearly for all the models down to 200–250 km. Degree 3 is also strong around 150 km depth in both SEMUCB-UMQ and QRFSI12. Figure 11. Open in new tabDownload slide Spherical-harmonic power spectral density comparison of the three global attenuation models shown in Fig. 9. Power spectral densities are computed for |$\delta Q^{-1}_\mu \times 10^3$|. Figure 11. Open in new tabDownload slide Spherical-harmonic power spectral density comparison of the three global attenuation models shown in Fig. 9. Power spectral densities are computed for |$\delta Q^{-1}_\mu \times 10^3$|. The correlations between the three attenuation models are shown in Fig. 12 for lmax = 8 and 12. In general, there is a better correlation between SEMUCB-UMQ and QRFSI12 compared to QRLW8, an older model which did not account for focusing effects either in the forward or in the inverse modeling stages. The correlation between SEMUCB-UMQ and QRFSI12 is particularly strong down to 150 km depth, but falls below 0.5 for depths greater than 250 km for lmax = 8 and 200 km for lmax = 12. The correlation between the three Qμ models at different angular orders is shown in Fig. S2 in the Supporting Information. Figure 12. Open in new tabDownload slide Correlation between the three global attenuation models in depth for |$\mathbf {l_{max} \in } \lbrace 8, 12\rbrace$|. The grey shaded area corresponds to correlation above 0.5. Figure 12. Open in new tabDownload slide Correlation between the three global attenuation models in depth for |$\mathbf {l_{max} \in } \lbrace 8, 12\rbrace$|. The grey shaded area corresponds to correlation above 0.5. In Fig. 13(a), we show the cross-correlation, as a function of depth, of 3-D structure in attenuation (Qμ) and in elastic parameters (VS, ξ) in model SEMUCB-UMQ, for lmax = 8, 12 and 16. Figure 13. Open in new tabDownload slide (a) Cross-correlation of |$Q^{-1}_\mu$| with the elastic parameters (VSiso, ξ) of SEMUCB-UMQ in depth. (b) Depth correlation of |$Q^{-1}_\mu$|. Figure 13. Open in new tabDownload slide (a) Cross-correlation of |$Q^{-1}_\mu$| with the elastic parameters (VSiso, ξ) of SEMUCB-UMQ in depth. (b) Depth correlation of |$Q^{-1}_\mu$|. There is overall a strong negative correlation between the velocity and attenuation models down to 250–300 km, with values of the correlation coefficient, in absolute value, larger than 0.5, that drop slightly with increasing |$\mathbf {l_{max}}$| and depth. The lateral variations in the radial anisotropy parameter (ξ) are positively correlated with those in |$Q^{-1}_\mu$|, but not very strongly. In Fig. 13(b), we show the depth correlation function for the |$Q^{-1}_\mu$| structure in SEMUCB-UMQ, which confirms the change in pattern between the upper 250 km of the mantle and deeper depths, as already noted in Gung & Romanowicz (2004). 5.2.4 Update of 3-D elastic structure Among the eight iterations, we carried out building model SEMUCB-UMQ, two were for the elastic parameters [(VSiso, ξ), Fig. 3]. In these elastic model updates, we relaxed the data selection criteria to those of French & Romanowicz (2014). However, our present data set is composed only of long-period surface waves in the period range of 400–60 s and lacks the body wave phases that were used in building SEMUCB-WM1. We have not included body waves in our elastic model updates, as we are focusing on the upper mantle and because of computational cost. Also, in their study, French & Romanowicz (2014) observed that in the upper mantle, the contribution of body waves was not significant beyond the first two iterations, which led them to fix the upper-mantle elastic model for their last few iterations and invert surface and body waveforms jointly only for the depths below 300 km. Figs 14 and 15 show a comparison of maps of lateral variations in VS and ξ, respectively, for SEMUCB-WM1 and SEMUCB-UMQ, as well as the differences between the two models, global average depth profiles and rms of the perturbations as a function of depth. Comparing the VS models, we see that our elastic model updates with only long-period surface waveforms tend to produce a smoother model by filtering out the short wavelength features. The strength of heterogeneity remains at the same level as in SEMUCB-WM1 down to 300 km and decreases below that, as shown in the rms depth variations. Although the general features seem to be stable, below 300 km, the difference in rms is more pronounced, which we attribute primarily to the lack of shorter period body wave constraints in our current data set. The global average VS profile remains stable. Figure 14. Open in new tabDownload slide (a) Map view of VS variations in per cent with respect to the mean VS value at a range of depths for the starting (SEMUCB-WM1) and final (SEMUCB-UMQ) models. The difference refers to |$100 \times (\delta \text{ ln}(V_S^{\text {SEMUCB-UMQ}})-\delta \text{ ln}(V_S^{\text {SEMUCB-WM1}}))$|. (b) The mean VS structures. (c) The normalized root-mean-square variation of δVS in depth with the peak values shown in the legend. Figure 14. Open in new tabDownload slide (a) Map view of VS variations in per cent with respect to the mean VS value at a range of depths for the starting (SEMUCB-WM1) and final (SEMUCB-UMQ) models. The difference refers to |$100 \times (\delta \text{ ln}(V_S^{\text {SEMUCB-UMQ}})-\delta \text{ ln}(V_S^{\text {SEMUCB-WM1}}))$|. (b) The mean VS structures. (c) The normalized root-mean-square variation of δVS in depth with the peak values shown in the legend. Figure 15. Open in new tabDownload slide (a) Map view of ξ variations expressed in per cent relative to isotropy (ξ = 1), at a range of depths in comparison to the starting elastic model (SEMUCB-WM1). The difference refers to |$100 \times (\xi ^{\text {SEMUCB-UMQ}} - \xi ^{\text {SEMUCB-WM1}})$|. (b) The mean ξ structures. (c) The normalized root-mean square variation of 100 × δξ with the peak value shown in the legend. Figure 15. Open in new tabDownload slide (a) Map view of ξ variations expressed in per cent relative to isotropy (ξ = 1), at a range of depths in comparison to the starting elastic model (SEMUCB-WM1). The difference refers to |$100 \times (\xi ^{\text {SEMUCB-UMQ}} - \xi ^{\text {SEMUCB-WM1}})$|. (b) The mean ξ structures. (c) The normalized root-mean square variation of 100 × δξ with the peak value shown in the legend. The changes in the ξ structure seem to be more pronounced than for VS, not only in the lateral distribution of structure but also in the global average. This may indicate stronger trade-offs between Qμ and ξ structure, which is not surprising given that both parameters are more poorly resolved than VS. Still, the main characteristics of the patterns in ξ are preserved, with predominantly positive anomalies (VSH > VSV) in the top 200 km of the mantle, and more isotropic continents than oceans (e.g Babuska & Plomerova 1999). Below 250–300 km, the signature of the western Pacific subduction zone (VSV > VSH) and to a lesser extent south America (e.g Panning & Romanowicz 2006; Moulik & Ekström 2014) is present in both models, more focused in SEMUCB-UMQ than in SEMUCB-WM1. An intriguing zone of (VSV > VSH) appears in SEMUCB-UMQ at 350 km depth under some hotspots in the south and west Pacific, also indicating the presence of vertical flow, while the ξ < 1 signal present in SEMUCB-WM1 under the east Pacific rise vanishes at shallower depths in SEMUCB-WMQ. Because our focus here is not on the elastic structure, the robustness of these new features in the ξ model will need to be confirmed in later studies, however, this also points to the possible influence of lateral variations in Qμ on elastic structure, which may manifest itself more significantly when inverting for ξ, possibly due to the additional dispersion effects on a poorly constrained elastic parameter. 5.3 Test of the effect of the crustal |$\mathbf {Q_\mu }$| model Mapping mantle structure accurately requires accounting for the effects of crustal structure on the waveforms. To that end, our group has been using a Monte Carlo-type inversion algorithm to build smooth, radially anisotropic crustal models with spatially varying thickness, that are equivalent to the real crust in the bandpass considered, using surface wave group velocity dispersion data. Further details on this technique and the rationale behind it can be found in Lekić & Romanowicz (2011) and French & Romanowicz (2014). Once the crustal model is constructed, it is used to calculate the spatially varying perturbations on the normal-mode frequencies through effective discontinuity kernels, as described in Lekić et al. (2010). Each time we update the upper-mantle structure, we adjust the crustal model, since even when only inverting for perturbations in the attenuation structure, the dispersion corrections thus introduced may change the uppermost mantle elastic structure, to which the dispersion data are sensitive. We fixed the crustal Qμ to a value of 300 as in QL6 (Durek & Ekström 1996). There is a general consensus that the mean attenuation in the crust is low. However, spatial variations of attenuation can be strong as shown by Mitchell (1995) and reviewed in Romanowicz & Mitchell (2015). We expect this crustal attenuation heterogeneity to affect long-period surface waves, especially Love waves, given their high sensitivity to near-surface structure. Here, we explore the potential influence of crustal attenuation 3-D structure on the recovered global upper-mantle Qμ model. For our test, we built a synthetic crustal attenuation model based on our crustal velocity model with Qμ values ranging from 80 to 2000, and a mean value of 290. Shown in Fig. 16, the test model follows the tectonic features closely, with high attenuation in regions of young, mid-age oceanic crust and backarcs, and low attenuation in cratons, stable parts of continents and regions of old oceanic crust. Figure 16. Open in new tabDownload slide (a) Heterogeneous crustal Qμ model used for evaluating the influence of the crustal Qμ structure on the synthetics and resulting mantle Qμ model. (b) T and Z component synthetics for the time windows corresponding to the arrival of different surface wave trains, computed for a path shown on the map for which details are given in the grey box, with and without the heterogeneous crustal Qμ model with SEMUCB-UMQ. Synthetics are filtered in the period range of 400–60 s. The minor and major arcs are indicated on the map by solid and dashed lines, respectively. Figure 16. Open in new tabDownload slide (a) Heterogeneous crustal Qμ model used for evaluating the influence of the crustal Qμ structure on the synthetics and resulting mantle Qμ model. (b) T and Z component synthetics for the time windows corresponding to the arrival of different surface wave trains, computed for a path shown on the map for which details are given in the grey box, with and without the heterogeneous crustal Qμ model with SEMUCB-UMQ. Synthetics are filtered in the period range of 400–60 s. The minor and major arcs are indicated on the map by solid and dashed lines, respectively. We computed synthetic seismograms using SEM in SEMUCB-UMQ, substituting this crustal Qμ model to the original 1-D Qμ crust, for the same data set as used in the construction of SEMUCB-UMQ, filtered in the same period range of 400–60s. Comparing these synthetic waveforms to the data, the heterogeneous crustal Qμ model had negligible effect on the number of selected wavepackets and variance reduction, showing that crustal Qμ model has no significant influence on this measure of fit, for the targeted period range. Presented in Fig. 16 are the time-windowed synthetic waveforms computed for a record for which the Love wave train G1 shows the highest waveform amplitude residual. For this great-circle path, the surface waves traveling along the minor-arc (G1, R1) sample areas in the crust with less-than-average attenuation whereas the ones traveling along the major-arc (G2, R2) sample regions with larger-than-average attenuation. The influence of the heterogeneous crustal attenuation model is clearly higher for the Love waves than the Rayleigh waves, as expected. To evaluate the influence of a heterogeneous crustal Qμ model on the recovered upper-mantle attenuation structure, we repeated our last iteration using the synthetics computed with the crustal test model. Fig. 17 shows the difference in upper-mantle attenuation structure between SEMUCB-UMQ and the model recovered. The crustal Qμ model slightly affects the upper-mantle attenuation model with decaying influence in depth. Wherever the crust has a higher/lower-than-average attenuation, we see an increase/drop in the attenuation perturbations in the upper mantle, respectively. These observations are in agreement with the tests of Dalton et al. (2008). A comparison of the attenuation perturbation differences presented in Fig. 17 with the absolute perturbation amplitudes in Fig. 9 shows that the effect of the crustal Qμ model is an order of magnitude smaller than the recovered heterogeneity within the period range considered. Figure 17. Open in new tabDownload slide Differences in upper-mantle Qμ maps for the recovered models with and without the 3-D crustal Qμ test model shown at the centre. The colour bar is saturated to make the small variations visible. Below 150 km, the perturbation differences are even smaller, and not visible with the present colour scale. Figure 17. Open in new tabDownload slide Differences in upper-mantle Qμ maps for the recovered models with and without the 3-D crustal Qμ test model shown at the centre. The colour bar is saturated to make the small variations visible. Below 150 km, the perturbation differences are even smaller, and not visible with the present colour scale. 5.4 Linear resolution analysis Here, we use the resolution matrix constructed from the approximate NACT-based-Hessian matrix obtained in the last iteration of our inversion, applied to two checkerboard patterns with block widths of 30° and 15°. Known to be valid only for linear problems, the type of resolution analysis presented here is commonly used for seismic tomography models under the assumption of being in the vicinity of global optimum model. Under this assumption, resolution analysis is considered to be sufficiently valid even in the case of non-linear inversion (Tarantola 2005). Although the analysis itself does not perfectly fit with our hybrid approach and tends to underpredict input model amplitudes (French & Romanowicz 2014; Lekić & Romanowicz 2011), it is a useful tool for assessing the resolving power of the model given the data quality, ray density coverage and a priori information. However, it does not address the accuracy of elastic corrections (focusing). Because model SEMUCB-WM1 is still relatively smooth (lateral heterogeneity resolved down to a scale of ∼600 km, (Meschede & Romanowicz 2015)), and because focusing effects depend on the second transverse derivatives of structure along the wavepath, and so are very sensitive to details in the elastic model, there may still be considerable room for improvement of the elastic model resolution. Figs 18 and 19 summarize the results of our resolution tests. In general, the input pattern is recovered fairly well down to 300 km depth, with degrading resolution and stronger lateral and in depth smearing at larger depths. The input model anomaly amplitudes, on the other hand, are not well recovered with the best peak-to-peak amplitude recovery of +80 per cent to 70 per cent at 80 km depth for the coarse checkerboard pattern (Fig. 18). Not presented here, the checkerboard test with block width of 12.5° failed to be recovered sufficiently well, which leads us to conclude that the lateral scales of interpretable structures in SEMUCB-UMQ are on the order of 15°. Figure 18. Open in new tabDownload slide Linear resolution analysis for the global upper-mantle Qμ structure in SEMUCB-UMQ at a range of depths, for two different checkerboard patterns with block widths of 30° and 15°. The maps are presented with multiple orientations, so that the structure at the poles can be seen clearly. The input models are shown in the grey box with the corresponding equatorial block widths (°). Figure 18. Open in new tabDownload slide Linear resolution analysis for the global upper-mantle Qμ structure in SEMUCB-UMQ at a range of depths, for two different checkerboard patterns with block widths of 30° and 15°. The maps are presented with multiple orientations, so that the structure at the poles can be seen clearly. The input models are shown in the grey box with the corresponding equatorial block widths (°). Figure 19. Open in new tabDownload slide Smearing in depth through linear resolution analysis for the global upper-mantle Qμ structure. Presented are the depth variations of the recovered checkerboard patterns at 120 km (left-hand column) and 300 km (right-hand column) depth. Figure 19. Open in new tabDownload slide Smearing in depth through linear resolution analysis for the global upper-mantle Qμ structure. Presented are the depth variations of the recovered checkerboard patterns at 120 km (left-hand column) and 300 km (right-hand column) depth. 6 DISCUSSION AND CONCLUSIONS We have presented a SEM-based upper-mantle anelastic model built through a hybrid full-waveform inversion method, starting from the 3-D radially anisotropic model SEMUCB-WM1 (French & Romanowicz 2014). Compared to previous global attenuation models built using surface waves (with or without overtones), presented by our or other groups, the new aspects in this study include: The use of synthetic seismograms computed using SEM in a relatively high-resolution 3-D elastic model, assuming that the mantle behaves as a linearly anelastic solid. The use of NACT-based Fréchet derivatives, which is important, compared to the standard PAVA approximation, for the modeling and inversion of overtone surface and body waves (e.g Mégnin & Romanowicz 1999). The use of an envelope full-waveform inversion method through a Gauss–Newton optimization scheme for the attenuation model updates. Being aware of the well-known limitations of the full-waveform inversion methods that stem from the non-linear nature of the inverse problem, we followed a careful strategy through an iterative inversion of the attenuation and elastic models, as well as the progressive expansion of the data set to include shorter periods and, in the last two iterations only, overtones. By the end of six attenuation and two elastic model updates, as well as one inversion iteration for source parameters, we obtained SEMUCB-UMQ, a 3-D Qmu model that extends down to 450–500 km depth. Our final average 1-D shear attenuation structure shows high attenuation in the depth range 80–220 km, in agreement with previous models (e.g. Anderson & Hough 1984; Durek & Ekström 1996; Lawrence & Wysession 2006). As our starting model did not include a high-attenuation zone, and recovering 1-D structure is the most non-linear part of long-wavelength seismic tomography, we see this recovery as an indication of the validity of our inversion approach. To first order, SEMUCB-UMQ shows similarities in its attenuation structure with QRLW8 of Gung & Romanowicz (2004) and QRFSI12 of Dalton et al. (2008) with patterns correlating with tectonic features down to 250 km, with low attenuation under cratons, shields and old oceanic crust and high attenuation along mid-ocean ridges and backarcs and a change in pattern at greater depths. A quantitative comparison between the latter three attenuation models through spectral analysis shows a better correlation between our final model SEMUCB-UMQ and QRFSI12 compared to QRLW8 with values of the correlation coefficient above 0.5 in the shallow upper mantle. In fact, the correlation coefficient reaches 0.8 around a depth of 100 km, which is remarkable, given that these two models have been constructed using a very different data set and a different theoretical approach. This is a promising step towards the recovery of robust 3-D global shear attenuation structure in the upper mantle. Below 250–300 km, SEMUCB-UMQ recovers two dominantly strong attenuation zones, one in the southwestern Pacific and other in eastern Africa, which also exist in QRFSI12 and QRLW8 with the latter showing the high-attenuation zone in Africa more to the south. However, the correlation between the recovered attenuation heterogeneity patterns in those three models is poor below 250 km. A comparison of the spherical harmonics degrees-2 and 3 components of the Qμ models and of the VS model of SEMUCB-UMQ is presented in Fig. 20. The geographic distributions of the peak positive Qμ anomalies at degree 2 appear slightly shifted to the west in the Pacific and Africa for QRLW8 and SEMUCB-UMQ, compared to VS and QRFSI12. A comparison of the degree-3 patterns reveals a closer correlation between VS of SEMUCB-UMQ and Qμ of QRLW8. Note, however, the generally good correlation between the degree-2+3 geographical pattern of SEMUCB-UMQ and QRFSI12. Figure 20. Open in new tabDownload slide Comparison of Qμ structure of SEMUCB-UMQ, QRLW8 and QRFSI12 and VS structure of SEMUCB-UMQ at 150 km depth for spherical harmonic degrees 2 (first row) and 3 (second row). The last row shows the pattern recovered by the addition of angular degrees 2 and 3. The perturbations are expressed in terms of |$\delta Q_\mu ^{-1}$|. They are normalized with respect to the absolute peak values and peak-to-peak values are shown on each panel. Figure 20. Open in new tabDownload slide Comparison of Qμ structure of SEMUCB-UMQ, QRLW8 and QRFSI12 and VS structure of SEMUCB-UMQ at 150 km depth for spherical harmonic degrees 2 (first row) and 3 (second row). The last row shows the pattern recovered by the addition of angular degrees 2 and 3. The perturbations are expressed in terms of |$\delta Q_\mu ^{-1}$|. They are normalized with respect to the absolute peak values and peak-to-peak values are shown on each panel. In SEMUCB-UMQ, we observe a strong negative correlation between VS and Qμ anomalies in the shallow upper mantle down to 250 km depth, which we interpret as indicating a common relation to surface tectonics, and the dominance of temperature effects. Including overtone surface waves into our data set, we tried to extend the sensitivity of our inversions to transition zone depths. The last two iterations that include overtone surface waves led to changes in the updated model below 250 km, bringing out a high-attenuation zone in eastern Africa and replacing low attenuation zones beneath some cratons and shields (e.g Canadian shield and Siberian craton), which were possibly smeared from shallower mantle structures due to low sensitivity of the fundamental modes below 300 km, by relatively higher attenuation. However, isolating the anelastic attenuation signal in overtone surface waveforms, which usually have lower amplitudes than fundamental-mode surface waves and are more sensitive to earthquake hypocentre depths, is a challenging task. Additionally, a reliable recovery of the deeper attenuation anomalies will require first, an improvement in the elastic model at those depths by including shorter period body waves, as in the work of French & Romanowicz (2014). As the elastic structure does not seem to be improving below 300–400 km using only surface waveforms, we did not extend the attenuation model beyond 450 km depth. Aside from starting our inversions with SEMUCB-WM1, we tried to improve the quality of our inversions and reduce the uncertainty inherent to the interpretation of seismic wave amplitudes through: (i) a careful data selection, (ii) joint inversion of |$\delta Q^{-1}_\mu$| and source and receiver terms, assigning a source term to each event and receiver terms to each component of ground motion and, (iii) updating the CMT source parameters, including both MT and hypocentre locations, using an NACT-based inversion algorithm, in the recovered 3-D heterogeneous anelastic models. The source inversions showed a correlation between the scalar moment magnitude variations and the geographic settings of the earthquake locations with magnitude loss in the continental areas and gain along the plate boundaries, in agreement with the Global CMT solution magnitude biases reported by Hjörleifsdóttir & Ekström (2010). As another source of uncertainty, we tested a heterogeneous crustal Qμ model and showed that within the period range of 400–60 s, the influence of crustal attenuation perturbations, even with Qμ values ranging from 80 to 2000, is low. Still, no matter how low they seem to be, some of the differences we observe between the compared attenuation models might be due to the crustal models considered, that differ both in elastic and anelastic parameters. Our final model’s VS structure is smoother than that of the starting model (SEMUCB-WM1) and shows weaker strength of heterogeneity below 300 km, due to the removal of short-wavelength anomalies that are not constrained by long-period surface waves. However, we do not consider this model as an update of SEMUCB-WM1 because of the lack of constraints from body waves. The changes in ξ are more pronounced, with visible changes even in the global average structure. While part of this difference may be due to the lack of constraints from body waves, it may also indicate that 3-D attenuation heterogeneity may trade-off with the less well-resolved elastic parameter ξ, through the indirect effect of dispersion corrections, and therefore, it may be important to consider lateral variations in global attenuation structure when inverting for radial anisotropy. Although updating elastic parameters leads to a smoother model, we have observed through our inversion procedure that it is a necessary step to increase the resolution of the attenuation model. This, we attribute primarily to the consideration of Qμ anomalies that leaked into the starting elastic model in its construction, since in their work, French & Romanowicz (2014) used only a fixed 1-D Qμ model. Therefore, by updating the elastic model, we achieved higher resolution in the attenuation model at the expense of smoothing the elastic model due to the lack of body waves in our data set. The resolution tests performed on SEMUCB-UMQ show that the size of interpretable attenuation anomalies in this model is about 15° in the shallow upper mantle and increases with depth. This might be partially due to the inaccuracies in the elastic model, the ray density coverage (data set size), as well as the limitation of the inversion method. Indeed, the fact that we do not include focusing in building the Fréchet derivatives matrix is a source of error which may be more significant with attenuation imaging than with the elastic model inversions. We consider SEMUCB-UMQ as a complement to SEMUCB-WM1 in the upper mantle and intend to improve it further. In future work, we plan to expand our data set by increasing the number of events in our catalogue, extending the bandpass to include body waveforms down to 30 s, as in the construction of SEMUCB-WM1, and extending the depth range of the Qμ inversion to lower mantle depths. This in turn would require a reliable crustal Qμ model, as the sensitivity of shorter period seismic waves to crustal heterogeneities increases. To that end, we intend to work towards building a crustal Qμ model which will serve initially to update our mantle elastic model and to construct a 3-D attenuation model that includes body waveforms. ACKNOWLEDGEMENTS The authors thank Vedran Lekić and Scott French who provided code and insights in the early stages of this project, and Wolfgang Friederich and two reviewers who contributed to improving the original manuscript. This work was supported by the European Research Council under the EC’s 7th Framework Programme (FP7-IDEAS-ERC)/ERC Advanced Grant WAVETOMO. BR acknowledges support from NSF grant EAR-1417229. The computations were performed on OCCIGEN of the HPC national facility CINES, in France, on EDISON and CORI at the National Energy Research Scientific Computing Center, USA, and the local cluster of Institut de Physique du Globe de Paris. Spectral analysis presented in this work relies on the SHTOOLS tool set built by Wieczorek et al. (2017). All waveform data used in this study are available from the Incorporated Research Institutions for Seismology (http://www.iris.edu). REFERENCES Adenis A. , Debayle E. , Ricard Y. , 2017a . Seismic evidence for broad attenuation anomalies in the asthenosphere beneath the pacific ocean , Geophys. J. Int. , 209 ( 3 ), 1677 – 1698 . https://doi.org/10.1093/gji/ggx117 Google Scholar Crossref Search ADS WorldCat Adenis A. , Debayle E. , Ricard Y. , 2017b . Attenuation tomography of the upper mantle , Geophys. Res. Lett. , 44 ( 15 ), 7715 – 7724 . https://doi.org/10.1002/2017GL073751 Google Scholar Crossref Search ADS WorldCat Anderson D.L. , Hart R. , 1978 . Attenuation models of the earth , Phys. Earth planet. Inter. , 16 ( 4 ), 289 – 306 . https://doi.org/10.1016/0031-9201(78)90068-7 Google Scholar Crossref Search ADS WorldCat Anderson D.L. , Ben-Menahem A. , Archambeau C.B. , 1965 . Attenuation of seismic energy in the upper mantle , J. geophys. Res. , 70 ( 6 ), 1441 – 1448 . https://doi.org/10.1029/JZ070i006p01441 Google Scholar Crossref Search ADS WorldCat Anderson J.G. , Hough S.E. , 1984 . A model for the shape of the fourier amplitude spectrum of acceleration at high frequencies , Bull. seism. Soc. Am. , 74 ( 5 ), 1969 – 1993 . WorldCat Babuska V. , Plomerova J. , 1999 . Seismic anisotropy: a method for studying the fabric of deep continental lithosphere , Global Tectonics Metallogeny , 7 ( 1 ), 3 . WorldCat Bao X. , Dalton C.A. , Ritsema J. , 2016 . Effects of elastic focusing on global models of rayleigh wave attenuation , Geophys. J. Int. , 207 ( 2 ), 1062 – 1079 . https://doi.org/10.1093/gji/ggw322 Google Scholar Crossref Search ADS WorldCat Bhattacharyya J. , Masters G. , Shearer P. , 1996 . Global lateral variations of shear wave attenuation in the upper mantle , J. geophys. Res. , 101 ( B10 ), 22 273 – 22 289 . https://doi.org/10.1029/96JB01782 Google Scholar Crossref Search ADS WorldCat Capdeville Y. , To A. , Romanowicz B. , 2003 . Coupling spectral elements and modes in a spherical earth: an extension to the sandwich case , Geophys. J. Int. , 154 ( 1 ), 44 – 57 . https://doi.org/10.1046/j.1365-246X.2003.01959.x Google Scholar Crossref Search ADS WorldCat Cara M. , Lévêque J.J. , 1987 . Waveform inversion using secondary observables , Geophys. Res. Lett. , 14 ( 10 ), 1046 – 1049 . https://doi.org/10.1029/GL014i010p01046 Google Scholar Crossref Search ADS WorldCat Dalton C.A. , Ekström G. , 2006 . Global models of surface wave attenuation , J. geophys. Res. , 111 ( B5 ), B05317 , doi:10.1029/2005JB003997 . https://doi.org/10.1029/2005JB003997 Google Scholar Crossref Search ADS WorldCat Dalton C.A. , Faul U.H. , 2010 . The oceanic and cratonic upper mantle: clues from joint interpretation of global velocity and attenuation models , Lithos , 120 ( 1 ), 160 – 172 . https://doi.org/10.1016/j.lithos.2010.08.020 Google Scholar Crossref Search ADS WorldCat Dalton C.A. , Ekström G. , Dziewoński A.M. , 2008 . The global attenuation structure of the upper mantle , J. geophys. Res. , 113 ( B9 ), doi:10.1029/2007JB005429 . https://doi.org/10.1029/2007JB005429 WorldCat 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 . https://doi.org/10.1093/gji/ggt365 Google Scholar Crossref Search ADS WorldCat Dalton C.A. , Bao X. , Ma Z. , 2017 . The thermal structure of cratonic lithosphere from global Rayleigh wave attenuation , Earth planet. Sci. Lett. , 457 , 250 – 262 . https://doi.org/10.1016/j.epsl.2016.10.014 Google Scholar Crossref Search ADS WorldCat 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 ( B10 ), B10308 , doi:10.1029/2012JB009288 . https://doi.org/10.1029/2012JB009288 Google Scholar Crossref Search ADS WorldCat Durek J.J. , Ekström G. , 1995 . Evidence of bulk attenuation in the asthenosphere from recordings of the bolivia earthquake , Geophys. Res. Lett. , 22 ( 16 ), 2309 – 2312 . https://doi.org/10.1029/95GL01434 Google Scholar Crossref Search ADS WorldCat Durek J.J. , Ekström G. , 1996 . A radial model of anelasticity consistent with long-period surface-wave attenuation , Bull. seism. Soc. Am. , 86 ( 1A ), 144 – 158 . WorldCat Durek J.J. , Ritzwoller M.H. , Woodhouse J.H. , 1993 . Constraining upper mantle anelasticity using surface wave amplitude anomalies , Geophys. J. Int. , 114 ( 2 ), 249 – 272 . https://doi.org/10.1111/j.1365-246X.1993.tb03914.x Google Scholar Crossref Search ADS WorldCat Dziewonski A.M. , Anderson D.L. , 1981 . Preliminary reference earth model , Phys. Earth planet. Inter. , 25 ( 4 ), 297 – 356 . https://doi.org/10.1016/0031-9201(81)90046-7 Google Scholar Crossref Search ADS WorldCat French S. , Romanowicz B. , 2014 . Whole-mantle radially anisotropic shear velocity structure from spectral-element waveform tomography , Geophys. J. Int. , 199 ( 3 ), 1303 – 1327 . https://doi.org/10.1093/gji/ggu334 Google Scholar Crossref Search ADS WorldCat French S. , Lekic V. , Romanowicz B. , 2013 . Waveform tomography reveals channeled flow at the base of the oceanic asthenosphere , Science , 342 ( 6155 ), 227 – 230 . https://doi.org/10.1126/science.1241514 Google Scholar Crossref Search ADS PubMed WorldCat Gilbert F. , Dziewonski A.M. , 1975 . An application of normal mode theory to the retrieval of structural parameters and source mechanisms from seismic spectra , Phil. Trans. R. Soc. Lond. A , 278 ( 1280 ), 187 – 269 . https://doi.org/10.1098/rsta.1975.0025 Google Scholar Crossref Search ADS WorldCat Gung Y. , Romanowicz B. , 2004 . Q tomography of the upper mantle using three-component long-period waveforms , Geophys. J. Int. , 157 ( 2 ), 813 – 830 . https://doi.org/10.1111/j.1365-246X.2004.02265.x Google Scholar Crossref Search ADS WorldCat Gung Y. , Panning M. , Romanowicz B. , 2003 . Global anisotropy and the thickness of continents , Nature , 422 ( 6933 ), 707 – 711 . https://doi.org/10.1038/nature01559 Google Scholar Crossref Search ADS PubMed WorldCat Hjörleifsdóttir V. , Ekström G. , 2010 . Effects of three-dimensional earth structure on cmt earthquake parameters , Phys. Earth planet. Inter. , 179 ( 3 ), 178 – 190 . https://doi.org/10.1016/j.pepi.2009.11.003 Google Scholar Crossref Search ADS WorldCat Jackson I. , Faul U.H. , Gerald F. , John D. , Tan B.H. , 2004 . Shear wave attenuation and dispersion in melt-bearing olivine polycrystals: 1. specimen fabrication and mechanical testing , J. geophys. Res. , 109 ( B6 ), B06201 , doi:10.1029/2003JB002406 . https://doi.org/10.1029/2003JB002406 Google Scholar Crossref Search ADS WorldCat Kanamori H. , Anderson D.L. , 1977 . Importance of physical dispersion in surface wave and free oscillation problems: review , Rev. Geophys. , 15 ( 1 ), 105 – 112 . https://doi.org/10.1029/RG015i001p00105 Google Scholar Crossref Search ADS WorldCat Karaoğlu H. , Romanowicz B. , 2018 . Global seismic attenuation imaging using full-waveform inversion: a comparative assessment of different choices of misfit functionals , Geophys. J. Int. , 212 ( 2 ), 807 – 826 . https://doi.org/10.1093/gji/ggx442 Google Preview WorldCat COPAC Karato S.-i. , 1993 . Importance of anelasticity in the interpretation of seismic tomography , Geophys. Res. Lett. , 20 ( 15 ), 1623 – 1626 . https://doi.org/10.1029/93GL01767 Google Scholar Crossref Search ADS WorldCat Komatitsch D. , Tromp J. , 2002 . Spectral-element simulations of global seismic wave propagation. -II. Three-dimensional models, oceans, rotation and self-gravitation , Geophys. J. Int. , 150 ( 1 ), 303 – 318 . https://doi.org/10.1046/j.1365-246X.2002.01716.x Google Scholar Crossref Search ADS WorldCat Komatitsch D. , Vilotte J.-P. , 1998 . The spectral element method: an efficient tool to simulate the seismic response of 2d and 3d geological structures , Bull. seism. Soc. Am. , 88 ( 2 ), 368 – 392 . WorldCat Lawrence J.F. , Wysession M.E. , 2006 . Qlm9: a new radial quality factor (q μ) model for the lower mantle , Earth planet. Sci. Lett. , 241 ( 3 ), 962 – 971 . https://doi.org/10.1016/j.epsl.2005.10.030 Google Scholar Crossref Search ADS WorldCat Lekić V. , Romanowicz B. , 2011 . Inferring upper-mantle structure by full waveform tomography with the spectral element method , Geophys. J. Int. , 185 ( 2 ), 799 – 831 . https://doi.org/10.1111/j.1365-246X.2011.04969.x Google Scholar Crossref Search ADS WorldCat Lekić V. , Matas J. , Panning M. , Romanowicz B. , 2009 . Measurement and implications of frequency dependence of attenuation , Earth planet. Sci. Lett. , 282 ( 1 ), 285 – 293 . https://doi.org/10.1016/j.epsl.2009.03.030 Google Scholar Crossref Search ADS WorldCat Lekić V. , Panning M. , Romanowicz B. , 2010 . A simple method for improving crustal corrections in waveform tomography , Geophys. J. Int. , 182 ( 1 ), 265 – 278 . WorldCat Li X.-D. , Romanowicz B. , 1995 . Comparison of global waveform inversions with and without considering cross-branch modal coupling , Geophys. J. Int. , 121 ( 3 ), 695 – 709 . https://doi.org/10.1111/j.1365-246X.1995.tb06432.x Google Scholar Crossref Search ADS WorldCat Li X.-D. , Romanowicz B. , 1996 . Global mantle shear velocity model developed using nonlinear asymptotic coupling theory , J. geophys. Res. , 101 ( B10 ), 22245 – 22272 . https://doi.org/10.1029/96JB01306 Google Scholar Crossref Search ADS WorldCat Liu H.-P. , Anderson D.L. , Kanamori H. , 1976 . Velocity dispersion due to anelasticity; implications for seismology and mantle composition , Geophys. J. Int. , 47 ( 1 ), 41 – 58 . https://doi.org/10.1111/j.1365-246X.1976.tb01261.x Google Scholar Crossref Search ADS WorldCat Liu Q. , Polet J. , Komatitsch D. , Tromp J. , 2004 . Spectral-element moment tensor inversions for earthquakes in southern california , Bull. seism. Soc. Am. , 94 ( 5 ), 1748 – 1761 . https://doi.org/10.1785/012004038 Google Scholar Crossref Search ADS WorldCat 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 . https://doi.org/10.1093/gji/ggv480 Google Scholar Crossref Search ADS WorldCat Mégnin C. , Romanowicz B. , 1999 . The effects of the theoretical formalism and data selection on mantle models derived from waveform tomography , Geophys. J. Int. , 138 ( 2 ), 366 – 380 . https://doi.org/10.1046/j.1365-246X.1999.00869.x Google Scholar Crossref Search ADS WorldCat 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 ( 3 ), 709 – 728 . https://doi.org/10.1046/j.1365-246X.2000.00298.x Google Scholar Crossref Search ADS WorldCat Meschede M. , Romanowicz B. , 2015 . Non-stationary spherical random media and their effect on long-period mantle waves , Mon. Not. R. Astron. Soc. , 203 ( 3 ), 1605 – 1625 . https://doi.org/10.1093/gji/ggv356 Google Scholar Crossref Search ADS WorldCat Mitchell B.J. , 1995 . Anelastic structure and evolution of the continental crust and upper mantle from seismic surface wave attenuation , Rev. Geophys. , 33 ( 4 ), 441 – 462 . https://doi.org/10.1029/95RG02074 Google Scholar Crossref Search ADS WorldCat Montagner J.-P. , Anderson D.L. , 1989 . Petrological constraints on seismic anisotropy , Phys. Earth planet. Inter. , 54 ( 1 ), 82 – 105 . https://doi.org/10.1016/0031-9201(89)90189-1 Google Scholar Crossref Search ADS WorldCat 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 . https://doi.org/10.1093/gji/ggu356 Google Scholar Crossref Search ADS WorldCat Panning M. , Romanowicz B. , 2006 . A three-dimensional radially anisotropic model of shear velocity in the whole mantle , Geophys. J. Int. , 167 ( 1 ), 361 – 379 . https://doi.org/10.1111/j.1365-246X.2006.03100.x Google Scholar Crossref Search ADS WorldCat Park J. , 1987 . Asymptotic coupled-mode expressions for multiplet amplitude anomalies and frequency shifts on an aspherical earth , Geophys. J. Int. , 90 ( 1 ), 129 – 169 . https://doi.org/10.1111/j.1365-246X.1987.tb00679.x Google Scholar Crossref Search ADS WorldCat Resovsky J. , Trampert J. , Van der Hilst R. , 2005 . Error bars for the global seismic q profile , Earth planet. Sci. Lett. , 230 ( 3 ), 413 – 423 . https://doi.org/10.1016/j.epsl.2004.12.008 Google Scholar Crossref Search ADS WorldCat Romanowicz B. , 1987 . Multiplet-multiplet coupling due to lateral heterogeneity: asymptotic effects on the amplitude and frequency of the earth’s normal modes , Geophys. J. Int. , 90 ( 1 ), 75 – 100 . https://doi.org/10.1111/j.1365-246X.1987.tb00676.x Google Scholar Crossref Search ADS WorldCat Romanowicz B. , 1990 . The upper mantle degree 2: constraints and inferences from global mantle wave attenuation measurements , J. geophys. Res. , 95 ( B7 ), 11051 – 11071 . https://doi.org/10.1029/JB095iB07p11051 Google Scholar Crossref Search ADS WorldCat Romanowicz B. , 1994 . On the measurement of anelastic attenuation using amplitudes of low-frequency surface waves , Phys. Earth planet. Inter. , 84 ( 1–4 ), 179 – 191 . https://doi.org/10.1016/0031-9201(94)90040-X Google Scholar Crossref Search ADS WorldCat Romanowicz B. , 1995 . A global tomographic model of shear attenuation in the upper mantle , J. geophys. Res. , 100 ( B7 ), 12375 – 12394 . https://doi.org/10.1029/95JB00957 Google Scholar Crossref Search ADS WorldCat Romanowicz B. , Durek J.J. , 2000 . Seismological constraints on attenuation in the earth: a review , in Earth’s Deep Interior: Mineral Physics and Tomography from the Atomic to the Global Scale , pp. 161 – 179 . Google Scholar Crossref Search ADS Google Preview WorldCat COPAC Romanowicz B. , Mitchell B.J. , 2015 . Deep earth structure–qof the earth from crust to core-1.21. COPAC Romanowicz B. , Roult G. , Kohl T. , 1987 . The upper mantle degree two pattern: constraints from geoscope fundamental spheroidal mode eigenfrequency and attenuation measurements , Geophys. Res. Lett. , 14 ( 12 ), 1219 – 1222 . https://doi.org/10.1029/GL014i012p01219 Google Scholar Crossref Search ADS WorldCat Romanowicz B.A. , Panning M.P. , Gung Y. , Capdeville Y. , 2008 . On the computation of long period seismograms in a 3-d earth using normal mode based approximations , Geophys. J. Int. , 175 ( 2 ), 520 – 536 . https://doi.org/10.1111/j.1365-246X.2008.03914.x Google Scholar Crossref Search ADS WorldCat Sailor R.V. , Dziewonski A.M. , 1978 . Measurements and interpretation of normal mode attenuation , Geophys. J. Int. , 53 ( 3 ), 559 – 581 . https://doi.org/10.1111/j.1365-246X.1978.tb03760.x Google Scholar Crossref Search ADS WorldCat Selby N.D. , Woodhouse J.H. , 2000 . Controls on rayleigh wave amplitudes: attenuation and focusing , Geophys. J. Int. , 142 ( 3 ), 933 – 940 . https://doi.org/10.1046/j.1365-246x.2000.00209.x Google Scholar Crossref Search ADS WorldCat 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 . https://doi.org/10.1029/2001JB000257 Google Scholar Crossref Search ADS WorldCat Shapiro N. , Ritzwoller M. , 2002 . Monte-Carlo inversion for a global shear-velocity model of the crust and upper mantle , Geophys. J. Int. , 151 ( 1 ), 88 – 105 . https://doi.org/10.1046/j.1365-246X.2002.01742.x Google Scholar Crossref Search ADS WorldCat Steinberger B. , 2000 . Plumes in a convecting mantle: models and observations for individual hotspots , J. geophys. Res. , 105 ( B5 ), 11127 – 11152 . https://doi.org/10.1029/1999JB900398 Google Scholar Crossref Search ADS WorldCat Tarantola A. , 2005 . Inverse Problem Theory and Methods for Model Parameter Estimation , SIAM . Google Scholar Crossref Search ADS Google Preview WorldCat COPAC Tarantola A. , Valette B. , 1982 . Generalized nonlinear inverse problems solved using the least squares criterion , Rev. Geophys. , 20 ( 2 ), 219 – 232 . https://doi.org/10.1029/RG020i002p00219 Google Scholar Crossref Search ADS WorldCat Wang Z. , Dahlen F. , 1995 . Spherical-spline parameterization of three-dimensional earth models , Geophys. Res. Lett. , 22 ( 22 ), 3099 – 3102 . https://doi.org/10.1029/95GL03080 Google Scholar Crossref Search ADS WorldCat 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 ), ESE 6-1–ESE 6-11, 2342 . https://doi.org/10.1029/2001JC001123 Google Scholar Crossref Search ADS WorldCat Widmer R. , Masters G. , Gilbert F. , 1991 . Spherically symmetric attenuation within the earth from normal mode data , Geophys. J. Int. , 104 ( 3 ), 541 – 553 . https://doi.org/10.1111/j.1365-246X.1991.tb05700.x Google Scholar Crossref Search ADS WorldCat Wieczorek M.A. , Meschede M. , Sales de Andrade E. , Oshchepkov I. , Xu B. , Walker A. , 2017 . SHTOOLS—Tools for working with spherical harmonics (v4.2) , doi:10.5281/zenodo.1067108 . Google Preview WorldCat COPAC Woodhouse J. , Wong Y. , 1986 . Amplitude, phase and path anomalies of mantle waves , Geophys. J. Int. , 87 ( 3 ), 753 – 773 . https://doi.org/10.1111/j.1365-246X.1986.tb01970.x Google Scholar Crossref Search ADS WorldCat Woodhouse J.H. , Dziewonski A.M. , 1984 . Mapping the upper mantle: three-dimensional modeling of earth structure by inversion of seismic waveforms , J. geophys. Res. , 89 ( B7 ), 5953 – 5986 . https://doi.org/10.1029/JB089iB07p05953 Google Scholar Crossref Search ADS WorldCat Zhou Y. , 2009 . Surface-wave sensitivity to 3-d anelasticity , Geophys. J. Int. , 178 ( 3 ), 1403 – 1410 . https://doi.org/10.1111/j.1365-246X.2009.04230.x Google Scholar Crossref Search ADS WorldCat Zhou Y. , Dahlen F. , Nolet G. , 2004 . Three-dimensional sensitivity kernels for surface wave observables , Geophys. J. Int. , 158 ( 1 ), 142 – 168 . https://doi.org/10.1111/j.1365-246X.2004.02324.x Google Scholar Crossref Search ADS WorldCat Zhu H. , Bozdağ E. , Duffy T.S. , Tromp J. , 2013 . Seismic attenuation beneath europe and the north atlantic: implications for water in the mantle , Earth planet. Sci. Lett. , 381 , 1 – 11 . https://doi.org/10.1016/j.epsl.2013.08.030 Google Scholar Crossref Search ADS WorldCat SUPPORTING INFORMATION Supplementary data are available at GJI online. Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper. Author notes Now at: Department of Earth Sciences, Universiteit Utrecht, Utrecht, the Netherlands. © The Author(s) 2018. Published by Oxford University Press on behalf of The Royal Astronomical Society.
TI - Inferring global upper-mantle shear attenuation structure by waveform tomography using the spectral element method
JO - Geophysical Journal International
DO - 10.1093/gji/ggy030
DA - 2018-06-01
UR - https://www.deepdyve.com/lp/oxford-university-press/inferring-global-upper-mantle-shear-attenuation-structure-by-waveform-1UWaCfCt1n
SP - 1536
VL - 213
IS - 3
DP - DeepDyve
ER -