1.2 METHODS
1.2.1 Mantle Tomography
We compare the results of the geodynamic modeling with the characteristics of four recent global tomographic VS models that differ in data selection, theoretical framework, and model parameterization and regularization choices. SEMUCB‐WM1 (French and Romanowicz, 2014) is a global anisotropic shear‐wave tomographic model based on full‐waveform inversion of three‐component long‐period seismograms, divided into wavepackets containing body waves and surface waves including overtones. Accurate waveform forward modeling is achieved using the spectral element method (Capdeville et al., 2005), 2D finite frequency kernels are computed using nonlinear asymptotic coupling theory (Li and Romanowicz, 1995), and the tomographic model is constructed by iterating using the quasi‐Newton method. SEISGLOB2 (Durand et al., 2017) is an isotropic shear‐wave model constrained by data sensitive primarily to VSV structure–spheroidal mode splitting functions, Rayleigh wave phase velocity maps, and body wave traveltimes. The inverse problem is fully linearized, with sensitivity of normal mode and Rayleigh wave data computed using first order normal mode perturbation theory, and traveltime data using ray theory, both computed in a reference spherically symmetric model. GLAD‐M15 (Bozdağ et al., 2016) is a global tomographic model that, like SEMUCB‐WM1, aims to fit full, three‐component seismic waveforms. Once again, the spectral element method (Komatitsch and Tromp, 2002) enables accurate forward modeling, while 3D finite frequency kernels are computed using the adjoint method at each iteration. GLAD‐M15 uses S362ANI (Kustowski et al., 2008) as a starting model and carries out an iterative conjugate gradient refinement of the 3D velocity structure to match observations while accounting for the full nonlinearity of the sensitivity of the misfit function to earth structure. The fourth model, S362ANI+M (Moulik and Ekström, 2014), is a global model based on full‐spectrum tomography, which employs seismic waveforms and derived measurements of body waves (∼1–20s), surface waves (Love and Rayleigh, ∼20–300s), and normal modes (∼250–3300s) to constrain physical properties – seismic velocity, anisotropy, density, and the topography of discontinuities – at variable spatial resolution. The inverse problem is solved iteratively and jointly with seismic source mechanisms, using sensitivity kernels that account for finite‐frequency effects for the long‐period normal modes, along‐branch coupling for waveforms and ray theory for body‐wave travel times and surface‐wave dispersion. Model SEISGLOB2 was obtained as spatial expansions of Voigt VS in netcdf format from the IRIS Data Management Center (http://ds.iris.edu/ds/products/emc‐earthmodels/http://ds.iris.edu/ds/products/emc‐earthmodels/). S362ANI+M was provided in netcdf format on a two‐degree grid and at 25 depths. GLAD‐M15 was provided on a degree‐by‐degree grid and at 132 depths by the authors (Bozdagˇ, pers. comm.). Model SEMUCB‐WM1 was obtained from the authors’ web page and evaluated using the model evaluation tool a3d
on a degree‐by‐degree grid at 65 depths.
The tomographic models were expanded into spherical harmonics to calculate power spectra and to facilitate wavelength‐dependent comparisons among models. First, at each depth, models were resampled using linear interpolation onto 40,962 equispaced nodes, providing a uniform resolution equivalent to 1 × 1 equatorial degree. Spherical harmonic expansions were carried out using the slepian MATLAB routines (Simons et al., 2006). Spherical harmonic coefficients were computed using the 4π‐normalized convention, such that the spherical harmonic basis functions Ylm satisfy
where alm and blm are real spherical harmonic coefficients (and bl0 = 0 for all l ). We note that Equation (1.1) is a measure of power per unit area per spherical harmonic degree, and not per spherical harmonic coefficient. This choice differs by a factor of 1/(2l + 1) from the definitions adopted in Dahlen and Tromp (1998) and used in Becker and Boschi (2002), and our definition omitting the multiplicative prefactor gives the appearance of a flatter power spectrum. Because we use this convention uniformly here, the choice of convention for the normalization of power does not affect our interpretations. From the spherical harmonic power spectra, we computed a spectral slope, which contains information about the relative amounts of power present at different wavelengths. Multiple definitions of spectral slope have been used with spherical harmonic functions, depending on whether spectral fitting is carried out in log‐log or log‐linear space. Here, the spectral slope is defined based on a straight‐line fit to
Radial correlation functions (Jordan et al., 1993; Puster and Jordan, 1994; Puster et al., 1995) were calculated from the spherical harmonic expansions. The RCF measures the similarity of δV structures at depths z and z′ as
(1.2)
where θ and ϕ denote the polar angle and azimuthal angle and Ω refers to integration over θ and ϕ. When working with normalized spherical harmonic functions, the above expression is equivalent to the linear correlation coefficient of vectors of spherical harmonic coefficients representing the velocity variations. Because the denominator of the expression for radial correlation normalizes by the standard deviations of the fields at both depths, the RCF is sensitive only to the pattern and not to the amplitude of velocity variations.
1.2.2 Mantle Circulation Models
We carried out a suite of 3D geodynamic models in spherical geometry using CitcomS version 3.1.1 (Zhong et al., 2000, 2008) with modifications to impose time‐dependent plate motions as a surface boundary condition (Zhang et al., 2010). CitcomS solves the equations of mass, momentum, and energy conservation for incompressible creeping (zero Reynolds number) flow under the Boussinesq approximation in 3D spherical shell geometry. All of the models include a compositionally distinct layer (advected using tracers), meant to be analogous to the LLSVP material, which is assigned an excess density of 3.75%, equivalent to a buoyancy ratio of B = 0.5. The intrinsic density difference adopted here is chosen such that the compositionally distinct material remains stable against entrainment and is consistent with other geodynamic modeling studies (Mc‐Namara and Zhong, 2004, 2005), leading to a net buoyancy compatible with the available constraints from normal modes and solid earth tides Moulik and Ekström (2016); Lau et al. (2017). The compositionally