Yamauchi & Takei anelastic shear wave velocity-temperature conversion benchmark#

This section was contributed by Fred Richards.

This benchmark tests the implementation of the anelastic shear wave velocity-to-temperature conversion derived by Yamauchi and Takei [2016] based on forced-oscillation experiments conducted at seismic frequencies on polycrystalline borneol. This anelasticity parameterization has been calibrated against a range of observational constraints on upper mantle temperature, attenuation and viscosity structure, using the surface wave tomography model of Priestley et al. [2012] to constrain shear wave velocity (\(V_S\)) and the plate model of McKenzie et al. [2005] to estimate lithospheric temperature structure. The resulting \(V_S\)-to-temperature conversion accurately accounts for the strongly non-linear temperature dependence of \(V_S\) at near-solidus conditions and is therefore especially useful for initializing models with accurate temperature structure in the upper \(\sim\) 400 km of the mantle. This benchmark is located in the folder benchmarks/yamauchi_takei_2016_anelasticity.

../../../../../_images/YT16_benchmark.png

Fig. 193 V_S as a function of temperature in the oceanic lithosphere. Dotted lines: digitized results from Fig. 20 of Yamauchi and Takei [2016]; solid lines: results; red = 50 km; blue = 75 km. Temperatures are taken from the plate model of McKenzie et al. [2005] and V_{S} from the surface wave tomography model of Priestley et al. [2012].#

The parameterization of Yamauchi and Takei [2016] defines \(V_S\) as

(49)#\[V_S = \frac{1}{\sqrt{\rho J_1}} \left( \frac{1+\sqrt{1+(J_2/J_1)^2}}{2}\right)^{-\frac{1}{2}} \simeq \frac{1}{\sqrt{\rho J_1}}\]

where \(\rho\) is the density and \(J_1\) and \(J_2\) represent real and imaginary components of the complex compliance, \(J^*\), which is a quantity describing the sinusoidal strain resulting from the application of a unit sinusoidal stress. \(J_1\) represents the strain amplitude in phase with the driving stress, whilst the \(J_2\) component is \(\frac{\pi}{2}\) out of phase, resulting in dissipation. Density is calculated using the expression

(50)#\[\rho (P,T) = \rho_{0} \left\{ 1 - \left[\alpha(T - T_0)\right] + \frac{P}{K} \right\}\]

where \(\rho_{0} = 3291~\text{ kg m}^{-3}\) and \(\alpha = 3.59 \times 10^{-5}~\text{ K}^{-1}\) are the density and thermal expansivity corresponding to \(T_{0} = 873~\text{ K}\), \(P\) = pressure and \(K = 115.2~\text{ GPa}\) is the bulk modulus. \(J_1\) and \(J_2\) are expressed as

(51)#\[\begin{split}\begin{aligned} J_1(\tau_S^{\prime})= & J_U \left[ 1 + \frac{A_B[\tau_S^{\prime}]^{\alpha_B}}{\alpha_B} + \frac{\sqrt{2\pi}}{2} A_P~\sigma_P \left\{ 1-\text{erf}\left(\frac{\ln[\tau_P^{\prime}/\tau_S^{\prime}]}{\sqrt{2}\sigma_P}\right)\right\}\right] \\ J_2(\tau_S^{\prime}) = & J_U \frac{\pi}{2} \left[ A_B[\tau_S^{\prime}]^{\alpha_B} + A_P~\exp \left(-\frac{\ln^2[\tau_P^{\prime}/\tau_S^{\prime}]}{2\sigma_P^2}\right)\right] + J_U \tau_S^{\prime} \end{aligned}\end{split}\]

where \(A_B = 0.664\) and \(\alpha_B = 0.38\) represent the amplitude and slope of background stress relaxation and \(J_U\) is the unrelaxed compliance. Parameters \(A_P\) and \(\sigma_P\) represent the amplitude and width of a high frequency relaxation peak superimposed on this background trend such that

\[\begin{split}A_P(T^{\prime}) = \begin{cases} 0.01 & \text{for }T^{\prime} < 0.91 \\ 0.01+0.4(T^{\prime}-0.91) & \text{for }0.91\leq T^{\prime} < 0.96 \\ 0.03 & \text{for }0.96\leq T^{\prime} < 1 \\ 0.03+\beta(\phi_m) & \text{for }T^{\prime} \geq 1 \\ \end{cases}\end{split}\]

and

\[\begin{split}\sigma_P(T^{\prime}) = \begin{cases} 4 & \text{for }T^{\prime} < 0.92 \\ 4+37.5(T^{\prime}-0.92) & \text{for }0.92\leq T^{\prime} < 1 \\ 7& \text{for } T^{\prime} \geq 1 \\ \end{cases}\end{split}\]

where \(T^{\prime}\) is the homologous temperature (\(\frac{T}{T_{s}}\)) with \(T\) the temperature and \(T_s\) the solidus temperature, both in Kelvin. \(\phi_m\) is the melt fraction and \(\beta(\phi_m)\) describes the direct poroelastic effect of melt (assumed to be negligible under upper mantle conditions). For this case, \(J_U\) is the inverse of the unrelaxed shear modulus, \(\mu_{U}(P,T)\), such that

\[J_{U}(P,T)^{-1} = \mu_U(P,T) = \mu_U^0 + \frac{\partial{\mu_U}}{\partial{T}}(T -T_0)+ \frac{\partial{\mu_U}}{\partial{P}}(P-P_0)\]

where \(\mu_U^0\) is the unrelaxed shear modulus at surface pressure-temperature conditions, the differential terms are assumed to be constant and the pressure, \(P\), in GPa is linearly related to the depth, \(z\), in km by \(\frac{z}{30}\). The normalised shear wave period, \(\tau_S^{\prime}\), in (51) is equal to \(\frac{\tau_S}{2\pi\tau_M}\), where \(\tau_S = 100~\text{ s}\) is the shear wave period and \(\tau_M = \frac{\eta}{\mu_U}\) is the normalised Maxwell relaxation timescale. \(\tau_P^{\prime}\) represents the normalised shear-wave period associated with the centre of the high frequency relaxation peak, assumed to be \(6 \times 10^{-5}\). The shear viscosity, \(\eta\), is

(52)#\[\eta = \eta_r \left(\frac{d}{d_r}\right)^{m} \exp \left[ \frac{E_a}{R}\left(\frac{1}{T}-\frac{1}{T_r}\right) \right] \exp \left[ \frac{V_a}{R}\left(\frac{P}{T}-\frac{P_r}{T_r}\right) \right] A_{\eta}\]

where \(d\) is the grain size, \(m\) the grain size exponent (assumed to be 3), \(R\) the gas constant, \(E_a\) the activation energy and \(V_a\) the activation volume. Subscripts \([X]_r\) refer to reference values, assumed to be \(d_r = d = 1~\text{ mm}\), \(P_r = 1.5~\text{ GPa}\) and \(T_r = 1473~\text{ K}\) for the upper mantle. \(A_{\eta}\) represents the extra reduction of viscosity due to an increase in \(E_a\) near the solidus, expressed as

\[\begin{split}A_\eta(T^{\prime}) = \begin{cases} 1 & \text{for } T^{\prime} < T^{\prime}_{\eta} \\ \exp \left[-\frac{(T^{\prime}-T^{\prime}_{\eta})}{(T^{\prime}-T^{\prime}T^{\prime}_{\eta})} \ln(\gamma)\right]& \text{for }T^{\prime}_{\eta} \leq T^{\prime} < 1 \\ \gamma^{-1} \exp(\lambda\phi) & \text{for }T^{\prime} \geq 1\\ \end{cases}\end{split}\]

where \(T^{\prime}_\eta\) is the homologous temperature above which activation energy becomes \(E_a + \Delta E_a\), and \(\gamma = 5\) is the factor of additional reduction. \(\lambda\phi\) describes the direct effect of melt on viscosity, assumed to be negligible here. The solidus temperature, \(T_s\), is fixed to a value of 1599 K at 50 km equivalent to a dry peridotite solidus [Hirschmann, 2000] and linearly increases below this depth according to

(53)#\[T_{s}(z) = 1599 + \frac {\partial T_s}{\partial z} (z - 50000)\]

where \(\frac {\partial T_s}{\partial z}\) is the solidus gradient.

In this benchmark, a 2D input ASCII file containing \(V_S\) specified at 50 and 75 km depth is read in and converted to temperature using an initial temperature model that implements the \(V_{S}(P,T)\) formulation detailed above. The default parameters governing the relationship between \(V_S\) and temperature are set to the values calibrated by Yamauchi and Takei [2016], where \(\mu_U^0 = 72.45~\text{ GPa}\), \(\frac{\partial{\mu_U}}{\partial{T}} = -0.01094~\text{ GPa K}^{-1}\), \(\frac{\partial{\mu_U}}{\partial{P}} = 1.987\), \(\eta_r = 6.22 \times 10^{21}~\text{ Pa s}\), \(E_a = 452.5~\text{ kJ mol}^{-1}\), \(V_a = 7.913 \times 10^{-6}~\text{ m}^{3}~\text{ mol}^{-1}\) and \(\frac{\partial T_s}{\partial z} = 1.018~\text{ K km}^{-1}\). As \(V_S\) is a complex function of temperature, a Brent minimization algorithm is used to find optimal values. Fig. 193 shows that the implementation of this parameterization can accurately recreate the results shown by Yamauchi and Takei [2016] in their Fig. 20.

The parameter file and initial temperature model for this benchmark can be found at benchmarks/yamauchi_takei_2016_anelasticity/yamauchi_takei_2016_anelasticity.prm and benchmarks/yamauchi_takei_2016_anelasticity/anelasticity_temperature.cc. Code to recreate Fig. 193 is provided in benchmarks/yamauchi_takei_2016_anelasticity/plot_output.