Inference of thermodynamic state in the asthenosphere from anelastic properties , with applications to North American upper mantle

5 Inference of thermodynamic state and full-spectrum mechanical behavior of the litho6 sphere and asthenosphere is a central problem in geophysics, implicating our understand7 ing of the convection patterns, transient responses and chemical composition of the planet. 8 Anelasticity is responsible for significant relaxation of stress associated with seismic wave 9 propagation in the asthenosphere, while irreversible transient creep may be important in the 10 lithosphere. This paper focuses on the processes that may act at the time scales of seismic 11 wave propagation, and current questions in the effort to determine the dependence of these 12 effects on thermodynamic state. We introduce a free code library, the “Very Broadband Rhe13 ology calculator” (VBRc), designed to calculate frequency-dependent mechanical properties 14 and easily compare different constitutive models favored by different laboratories. The meth15 ods operate only in the forward sense, starting with arrays of models of thermodynamic state, 16 proceeding to arrays of mechanical properties. These calculations are incorporated into a 17 Bayesian framework to infer variation in mantle thermodynamic state from Vs and Q, ap18 plied here to four locations in Western North America. The results demonstrate how well we 19 can constrain the state, given the input models and the measurements. Results for sites in the 20 Basin and Range, Colorado Plateau and interior craton east of the Rio Grande separate into 21 distinct state variable ranges consistent with their tectonic environments. 22 *corresponding author 23


24
The aim of this paper is to present an integrative framework for calculating the effects of anelas-25 ticity on seismic velocity and attenuation, and then apply a Bayesian inference framework to 26 several sites in western North America. Inference of thermodynamic state of the upper mantle 27 is central to understanding the mechanics of the lithosphere, the spatial variations in the degree 28 of mechanical coupling between plate motions and convection patterns, and any questions of 29 melting productivity and extraction physics. It is also critical for understanding surficial expres- geophysics, from completely unrelaxed to completely relaxed. In this paper, we only focus on 55 the inference of thermodynamic state from seismic measurements within the seismic band. The 56 code structure is also designed to be used to develop new constitutive models from laboratory 57 data, such that, by virtue of being a public code repository, new models can quickly be used by 58 geophysicists to interpret their measurements. 59 As the calculation is an entirely forward calculation, it can easily be used in the context of 60 a Bayesian inference approach to infer thermodynamic state over some representative volume contiguity formulations, the VBRc is implemented towards the aims described by Takei (2013). and observed crystal lattice preferred orientations and microstructures in xenoliths and ophiolites, 169 it may be that the subgrain size is the appropriate lengthscale for estimating the HTB attenuation, 170 with minor or significant additional effects coming from the grain boundary structure. 171 That said, the similitude (collapse of data onto the master curve by normalizing by the Maxwell  Below, we discuss various potential influences of melt, water and second phases on the HTB, 185 and then additional mechanisms that can elevate attenuation levels above the HTB.  T , and solidus, T s , the dramatic weakening attributed to melt begins at about T /T s = 0.95, or 201 95% of the melting temperature. This subsolidus weakening appears in the steady state viscosity 202 as well as the attenuation measurements. It is referred to as premelting and is attributed to in-203 creased grain boundary diffusivity due to a highly local increase in disorder rather than a direct 204 effect of a distinct melt phase (Takei, 2019), illustrated in Fig 1a.2. There is also good evidence pative mechanism, with the grain boundary's viscosity giving rise to a energy loss with a narrow 239 characteristic frequency. Transient diffusion creep invariably involves small-displacement sliding 240 on grain boundaries (hence the term "diffusionally assisted GBS") (e.g. Raj, 1975;Cooper, 2002;241 Morris and Jackson, 2009;Faul and Jackson, 2015)). In discussions of HTB attenuation, when 242 data exhibit only a Q ∝ f 1/3 behavior, this sliding is assumed to be on effectively inviscid grain 243 boundaries; the sliding dissipation is overwhelmed by the transient diffusion creep and is ignored.  attenuation in meso-scale structures such as organized melt networks ( Fig. 1c2-c4) has not yet 293 been quantitatively modeled, to our knowledge.

294
Another additional possible mechanism involves different local thermal changes in phases, 295 due to their different thermal expansion coefficients, driven by seismic strains (e.g. Budiansky 296 et al., 1983; Chrysochoos, 2012). The frequency dependence of the dissipation emerges due to 297 the relationship between the wave period and thermal exchange coefficients between adjacent 298 phases or anisotropic phases. This mechanism has been only minimally explored. In the lab, a single experiment can measure the elastic, transient and steady state viscous proper-301 ties. In practice, to calculate a full-spectrum mechanical response for rocks, an inherent source 302 of complexity and uncertainty is the extrapolation from lab-to-earth conditions, and to the wide 303 array of natural rock compositions and thermodynamic conditions that can only be minimally ex-304 plored in the lab. Self-consistent calculation across the wide array of parameters can be achieved 305 approximately, and we wish to point out here two aspects that emerged in the above discussions, the grain boundary structure, not the subgrain structure, to first order.

316
(2) While it is tempting to estimate effects of varying mantle composition by the simplest 317 route: varying anharmonic elastic moduli and density as functions of composition and then su-318 perimposing anelastic behavior calculated for olivine rocks, the effects of phase boundaries will 319 be ignored, which may be significantly incorrect. Composition effects need to be accounted for 320 across the entire spectrum, from unrelaxed to steady state creep, but may cause an increase in 321 attenuation that would not be captured by the end-member effects on the Maxwell time alone, 322 though that is the place to start.

323
Clearly, there are many open questions. We do not try to include these effects or calculate 324 any uncertainty from the constitutive models-doing so is beyond the scope of this paper. We 325 also remain agnostic to the different fitting and scaling models employed here, but in future work 326 intend to carry out comparisons of fitting models across laboratory data, towards convergence on 327 the many questions above.   we treat anharmonicity as simply as possible using a linear scaling of a generic modulus M from 356 reference pressure P R and temperature T R : The VBRc includes sets of anharmonic derivatives from Isaak (1992) and Cammarano et al.

358
(2003) for Fo90 olivine and in this paper, we employ those of Isaak (1992): dG/dT = −13.6 359 MPa • K and dG/dP = 1.8 with a value of G 0 = 80 at standard temperature and pressure.  kinetic parallel or mechanical series, so the total strain rate is then given byε = ε i , where an 377 individual mechanism's strain rate is given by (2) i = 1 for diffusion creep, i = 2 for dislocation creep and i = 3 for GBS creep. (sgn(σ) = sgn(ε) 379 is implied.) d is grain size, Q i is the thermal activation energy, σ is the differential stress, and 380 V * is the activation volume. Given the strain rate for each mechanism, the VBRc calculates the 381 viscosity of a single mechanism, η i = σ/ε i , as well as the total effective viscosity using the where η 0 is the flow law for the subsolidus viscosity and λ may vary for each deformation mech-389 anism, as melt will affect them differently.
where φ c is the critical melt fraction and x c is the weakening amplitude across φ c . It behaves In this section, we describe the approaches used in fitting and scaling anelastic models to ex-420 perimental data and emphasize those implemented in the VBRc. In general, the approach is to 421 fit an anelastic model to experimental data, which can then be used to scale to conditions in the 422 Earth. The models are generally different definitions of the "relaxation spectrum" as a function 423 of period, X(t) (Takei, 2013): where ∆ is the "relaxation strength" and X(τ ) is the relaxation spectrum. The methods described 425 below are all approximations of the full relaxation spectrum.

426
The VBRc currently implements four anelastic methods encompassing a range of models and where η ss is the steady state viscosity. Taking the Laplace transform of the creep function yields 443 the storage and loss compliances, J 1 and J 2 as a function of angular frequency ω: the Maxwell time and β * = β/J U .

446
Following Jackson and Faul (2010), the "pseudoperiod master variable" approach for scaling 447 from laboratory to earth conditions substitudes a master variable X a for the period in the angular 448 frequency, ω = 2π/τ = 2π/X a . The master variable X a is a function of the state variables 449 measured from a reference state: where τ 0 is the period of the oscillation.The VBRc implementation adds a dependence on melt 451 fraction using the diffusion creep values for melt weakening from section 2.2: where D(τ ) is the distribution of relaxation times of the series of Kelvin elements and ∆ is the other than the steady state process.

472
Multiple relaxation mechanisms may be superimposed with different relaxation strengths.

473
The distributions for the high temperature background attenuation of strength ∆ B and a dissipa-474 tion peak with relaxation strength ∆ P are given by where 0 < α < 1 when the relaxation time is between the high and low limits, τ L < τ < 476 τ H and the dissipation peak is a Gaussian distribution described by τ p and σ p . Including these 477 two distributions in the creep function and transforming to the frequency domain results in the 478 following storage and loss complicances: Scaling to other conditions is achieved through scaling the Maxwell time, τ M , lower and upper 480 integration limits τ L and τ H , and dissipation peak time τ P from reference values given activation 481 energy E and activation volume V : with i = L, H, P, M , grain size d, tempeature T , pressure P .

483
The VBRc includes four sets of fitting parameters that define the free parameters ∆ B , α, τ M R, and the relaxation spectrum is empirical is given as a piecewise function that depends on the 503 Maxwell-normalized period, τ : where τ = τ f M , with the Maxwell frequency f M given by where E u is the unrelaxed modulus and η dif f is the steady state diffusion creep viscosity of the 506 material.
where A p and σ p are both piecewise functions of T n : The above relaxation spectrum results in the following relationships for J 1 and J 2 : Constant parameters include: A B , α and τ P n .

527
The steady state Maxwell time is given by τ M = η/G U (T, P ) where η is the steady state 528 diffusion creep viscosity. Yamauchi and Takei (2016) introduce a scaling for the viscosity that 529 includes a dependence on T n : where η(T, P, d) is the viscosity at the current thermodynamic state (temperature, pressure, grain 531 size)neglecting any pre-melt and direct melt effects and A η (T n ) has the form where γ and T η n are fitting constants and λ is the steady state exponential melt dependence. Ya-  is calculated with appropriate constants for either composition and A η (T n ) is the same for both.
Note that in the remainder below, we frequently refer to the quality factor, Q, the inverse of 548 attenuation in lieu of Q −1 as Q varies with state variables in the same sense as V s : e.g., an 549 increase in temperature decreases both Q and V s .
The different anelastic models are influenced by melt fraction in different ways, as described in 572 Section 2.3 and summarized in Table 1.

573
From Fig. 4, it is clear that the poro-elastic effect must be included when interpreting observed 574 measurements. Fig. 4a and Fig. 4b     To better understand the grain size dependence of the various methods, we pick two tempera-   not only between location but also between methods.

650
Towards that end, we used the VBRc to vary φ,d and T between possible mantle values:

690
• p(s): the "prior models", the probability of sampling each set of state variables.

691
• p(m): the "measurement probability", the probability of observing the measurement itself.

692
In practice, p(m) is not known a priori but given that it is a normalizing constant, p(m) can 693 be neglected to calculate relative probabilities, or p(m) can be calculated by summing the 694 final relative probabilities.

695
The likelihood, p(m|s), comes from the χ 2 -misfit, which we use to construct a Gaussian

700
The prior probability, p(s), comes from our a priori knowledge of the state variables. In 701 the present study, we start with a uniform distribution across the φ, T and d parameter sweep.

702
However, given other constraints, e.g. geothermobarometry melt-and/or xenoliths, it is possible 703 to put in a more tightly constrained prior, with its own uncertainty, as applied below to the grain 704 size.
where we have written the p(V s , Q) as p o to emphasize its role as a normalization constant.

779
The present study is concerned with how the anelastic scalings may influence inferred state vari-780 ables rather than comparison of seismic models, of which several detailed studies exist (e.g.,

781
Cammarano and Guerri, 2017), or geodynamic interpretations. As such, we restrict our measure-  each anelastic method to grain size leads to larger differences at the smaller grain size. This is 912 visible in Fig. 11   (1) a propagating shear wave will have energy attenuated at the grain scale (2) at high temperature by transient diffusion creep, that occurs when sliding on a grain boundary causes peaks in traction at grain corners (red lines) that drive rapid local diffusion and diminish as the diffusion relaxes stress towards the traction profile of steady state creep (green lines). (3,4) As the wave arrives and passes, tractions develop on one set of grain edges and then switch; the total dissipation depends on frequency. (b) Properties that affect the HTB: (1) melt tubules (with topology determined by surface tension) at grain triple junctions can aid in rapid diffusion.
(2) premelting-or sub-solidus disordering of the grain boundary-leads to increased diffusivity which also relaxes traction peaks more quickly, as can (3) water-related defects in crystals or on grain boundaries and (4) phase boundaries. (c) secondary dissipative mechanisms include (1) elastically accommodated grain boundary sliding (eaGBS), (2) dislocation damping, (3) melt squirt, and (4) potential meso-scopic structures.      Table (LUT) slices showing tradeoffs between temperature, grain size, and melt fraction at a fixed pressure and frequency band for shear wave velocity V s (top row) and quality factor Q (bottom row). These calculations use the andrade psp anelastic method.  : Probability distributions of melt fraction φ, temperature T and grain/subgrain size d for the Basin and Range using the Andrade pseudoperiod scaling (andrade psp) for three cases: (top row) single parameter inference using shear wave velocity V s : p(φ, T, d|V s ), (middle row) single parameter inteference using quality factor Q p(φ, T, d|Q), and (bottom row) joint inference using V s and Q: p(φ, T, d|V s , Q). The 2D plots are the probability summed over the third variable that is not plotted and the 1D plots are the marginal probability of the that thrid variable. The bottom pressure-depth plots show the pressure range of the observation. Figure 10: Joint probability distribution p(φ, T, d|V s , Q) for the Basin and Range using the Andrade pseudoperiod scaling and assuming a log-normal distribution for the prior model of grain/subgrain size with median grain/subgrain size of (top row) 1 cm and (bottom row) 1 mm. See figure 9 for an explanation of the 2D and 1D plots. Figure 11: Likely melt fraction -temperature fields for each anelastic method and location for the two prior model cases in figure 10. Line color corresponds to location (abbreviations defined in Fig. 8), line thickness corresponds to the probability distribution interval (e.g., the thickest lines contain 70% of the distribution), and line style corresponds to the median grain/subgrain size of the prior model. Figure 12: Likely melt fraction -temperature fields for each location using an ensemble weighting of the probability distributions for each anelastic method. See the legend and caption of figure  11 for description of line colors and styles.