Mantle Convection and Surface Expressions

(1.1)

where a lmand b lmare real spherical harmonic coefficients (and b l0= 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/(2 l + 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 vs. l in log‐log coordinates, for spherical harmonic degrees 2–20. A slope of zero would indicate uniform power at all spherical harmonic degrees (the power spectrum of a delta function on the sphere (e.g., Dahlen and Tromp, 1998)), and increasingly negative slopes indicate more rapidly decreasing power at shorter wavelengths.

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 δ Vstructures 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.

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 distinct material is initially present as a uniform layer of 250 km thickness. All of the models include time‐dependent prescribed surface plate motions, which shape the large‐scale structure of mantle flow. We adopt plate motions from a recent paleogeographic reconstruction by Matthews et al. (2016), which spans 410 Ma‐present, although some of the calculations do not include the entire plate motion history. All of the models except Case 40 ( Table 1.1) impose the initial plate motions for a period of 150 Myr to spin‐up the model and initialize large‐scale structure following Zhang et al. (2010). The mechanical boundary conditions at the core‐mantle boundary are free‐slip, and the temperature boundary conditions are isothermal with a nondimensional temperature of 0 at the surface and 1 at the core–mantle boundary. We use a temperature‐ and depth‐dependent viscosity with the form η ( z ) = η z( z ) exp [ E (0.5 − T )], where η z( z ) is a depth‐dependent viscosity prefactor and E = 9.21 is a dimensionless activation energy, which gives rise to relative viscosity variations of 10 4due to temperature variations. The models are heated by a combination of basal and internal heating, with a dimensionless internal heating rate Q = 100.

We include the effects of a phase transition at 660 km depth in some of the models. Phase transitions are implemented in CitcomS using a phase function approach (Christensen and Yuen, 1985). We adopt a density increase across 660 km of 8%, a reference depth of 660 km, a reference temperature of 1573 K, and a phase change width of 40 km. We assume a Clapeyron slope of –2 MPa/K. Recent experimental work favors a range of –2 to –0.4 MPa/K (Fei et al., 2004; Katsura et al., 2003), considerably less negative than values employed in earlier geodynamical modeling studies that produced layered convection (Christensen and Yuen, 1985). The models shown in the present work are a subset of a more exhaustive suite of models from Lourenço and Rudolph (in review), which consider a broader range of convective vigor and additional viscosity structures. We list the parameters that are varied between the five models in Table 1.1and the radial viscosity profiles used in all of the models are shown in Figure 1.4.

We carried out inversions for the mantle viscosity profile constrained by the long‐wavelength nonhydrostatic geoid. The amplitude and sign of geoid anomalies depend on the internal mantle buoyancy structure as well as the deflection of the free surface and core‐mantle boundary, which, in turn, are sensitive to the relative viscosity variations with depth (Richards and Hager, 1984; Hager et al., 1985). Because the long‐wavelength geoid is not very sensitive to lateral viscosity variations (e.g., Richards and Hager, 1989; Ghosh et al., 2010), we neglect these, solving only for the radial viscosity profile. The geoid is not sensitive to absolute variations in viscosity, so the profiles determined here show only relative variations in viscosity, and absolute viscosities could be constrained using a joint inversion that includes additional constraints such as those offered by observations related to glacial isostatic adjustment. In order to estimate the viscosity profile, we first convert buoyancy anomalies from mantle tomographic models into density anomalies and then carry out a forward model to generate model geoid coefficients. We then compare the modeled and observed geoids using the Mahalanobis distance

(1.3)

where denotes a vector of geoid spherical harmonic coefficients calculated from the viscosity model with parameters , is the data‐plus‐forward‐modeling covariance matrix. The Mahalanobis distance is an L 2‐norm weighted by an estimate of data+forward modeling uncertainty, and is sensitive to both the pattern and amplitude of misfit.