Periclase deforms more slowly than bridgmanite under mantle conditions – Nature

Dislocation dynamics calculations

The 2.5D DD simulation method is used throughout this study to model the creep properties of periclase under lower mantle conditions under a constant creep stress ranging from 10 to 150 MPa. We consider the P, T conditions on the geotherm of Stacey and Davis¹² at 30, 60, 90 and 120 GPa, and use the periclase elastic properties from Karki et al.³⁴ (Extended Data Table 2).

To capture the statistically representative volume element of the microstructure, dislocations of Burgers vector b = (frac{1}{2})<110> are introduced edge-on in two intersecting slip systems within a square element of typical cell size L_x = L_y between 8 to 60 µm (function of the desired dislocation densities) and periodic boundary conditions are applied. As extensively described in some of our previous works¹¹, specific characteristics of a realistic MgO microstructure are introduced as local rules, including junction formation strengthening and a dislocation multiplication rate to allow random introduction of fresh dislocations of opposite signs on all slip systems. All local rules are parameterized according to experimental observations or three-dimensional DD simulations as described in Reali et al.¹¹.

For a realistic simulation of the creep behaviour, we use a combination of glide and climb events to control the displacement of dislocations in the simulation volume. The Burgers vector contained in the reference plane defines the slip direction, whereas the climb direction is defined as the one in the reference plane, orthogonal to the Burgers vector. The positive orientation of the climb direction is taken along the vacancy emission direction. During a given timestep of the simulation, dislocations are displaced along the glide or climb direction according to the velocity formulation presented below.

In our simulations, dislocations glide in the athermal regime according to a « free-flight » velocity v_g, which is a function of the effective resolved shear stress τ = τ_app+ τ_int (that is, including the applied resolved shear stress and the elastic interaction stress):

As, in the range of P and T conditions investigated here, the lattice friction of MgO vanishes^5,6, the glide velocity is thus expressed with a viscous drag coefficient B. Such a coefficient is known to be temperature and material dependent. Thus, we consider a linear scaling of B with temperature³⁵ and use as a starting point the experimental value recorded at 300 K³⁶.

Climb events occur through the self-diffusion of vacancies adsorbed or emitted by the dislocation line. To describe the climb rate in steady-state creep equilibrium conditions, we use a dislocation climb velocity v_c expressed as follows³⁵:

where τ_c is the climb stress, calculated similarly to τ but resolved here along the climb direction. (alpha =2pi /{rm{ln}}(R/{r}_{{rm{c}}})) is a geometrical factor that describes the cylindrical geometry of the vacancy flux field around the dislocation line, where R and ({r}_{{rm{c}}}) represent the two radii of the cylindrical surfaces through which the vacancy flux is calculated. Again, in the range of temperature investigated, the dislocation line must contain a high density of jogs, and thus be saturated with jogs, enabling vacancies to be absorbed or emitted instantaneously by the dislocation. Here, ({r}_{{rm{c}}}) is the dislocation core capture radius and R is taken as a fraction of the average dislocation distance. Being within the logarithmic term, the (R/{r}_{{rm{c}}}) ratio does not significantly affect the climb velocity values and here is taken as constant and equal to 100³⁷. Ω is the formation volume of vacancies in the Mg and O sites and is calculated from the unit cell volume of MgO—that is, (varOmega ={a}^{3}/Z) where Z = 4 is the number of formula units per unit cell. c_∞ and c₀ are the vacancy concentrations far from the dislocation and the equilibrium vacancy concentration in the bulk volume, respectively. Following Boioli et al.³⁷, we assume that far away from the dislocations (that is, on the external limit R of the cylindrical flux) the vacancy concentration is constant and equal to the equilibrium concentration in the bulk volume (c_∞= c₀). Finally, D_sd is the vacancy self-diffusion coefficient, which controls the flow of atomic species from and to the climbing dislocation.

Micromechanical model of the rheology of the aggregate

Local behaviour

We apply state-of-the-art formulations of nonlinear homogenization methods to estimate and bound the effective viscosity of a two-phase material constituted by bridgmanite and periclase with various volume fractions, at pressures and temperatures relevant for the lowermost part of Earth’s mantle. These phases are supposed to be randomly mixed so that the effective mechanical behaviour of the whole aggregate can be considered isotropic. Bridgmanite and periclase are considered to be two isotropic viscous phases too. In reality, they are both made of grains, but we consider that these grains are equiaxed and randomly oriented so that, on average, the polycrystalline aggregates can be replaced by two homogenous and isotropic phases. The role of grain or phase boundaries is not explicitly taken into account. We are interested in the effective rheological behaviour in the permanent (that is, secondary) creep regime. Thus, the elastic behaviour of both phases does not come into play. At the scale of each phase, the viscoplastic behaviour can be described by the following nonlinear creep law:

with ({dot{varepsilon }}_{0}) a reference strain rate, σ₀ a reference stress, n the stress sensitivity, ({dot{varepsilon }}_{{rm{eq}}}) the equivalent strain rate and σ_eq the equivalent stress, defined as:

with summation over repeated indices (Einstein convention), ({dot{varepsilon }}_{{ij}}) and ({sigma {prime} }_{{ij}}) components of the strain rate tensor and deviatoric stress tensor, respectively. We chose here to take ({dot{varepsilon }}_{0}=1times {10}^{-15}{{rm{s}}}^{-1}), that is, a value typical for lower mantle convection. Doing so, the σ₀ values are of the order of magnitude of the equivalent stress encountered in situ during mantle flow. Values of σ₀ for bridgmanite and periclase are fitted from the creep laws of Fig. 2 and available from the data deposit. The stress sensitivity n is the same for both phases, n = 3.1.

Effective behaviour

The effective behaviour reads in a similar way as the local behaviour:

where (dot{bar{varepsilon }}) and (bar{sigma }) are the effective strain rate and stress tensors, respectively, given by the volume average (denoted (leftlangle bullet rightrangle )) of local strain rate and stress fields: (dot{bar{varepsilon }}=leftlangle dot{varepsilon }rightrangle ) and (bar{sigma }=leftlangle sigma rightrangle ). (widetilde{n}) is the effective stress sensitivity, which can, in general, be inferred from the homogenization procedure. In the case considered here, things are simpler as bridgmanite and periclase exhibit the same n value. Therefore (widetilde{n}=n=3.1). Finally, ({widetilde{sigma }}_{0}) is the effective flow stress associated with the (nonlinear) viscosity of the aggregate in situ. It is determined by the homogenization procedure detailed below.

Homogenization procedure

To describe the possible rheology of an assemblage representative of the lower mantle, one can reasonably assume that bridgmanite and periclase phases are randomly mixed. In that case, a good estimation of the effective behaviour is provided by the Self-Consistent scheme (denoted SC below), based on a fully disordered microstructure^38,39. To extend the SC scheme to nonlinear rheologies, as here, we have made use of the very efficient and accurate linearization procedure (so-called partially optimized second-order procedure (POSO)) of Ponte Castaneda²⁹.

Sensitivity to calculation parameters

The calculations presented in the main text are made based on the diffusion coefficients for pure MgO periclase and for MgSiO₃ bridgmanite with a vacancy concentration of 10⁻⁵. For the calculation of the aggregate rheology, we considered a periclase proportion of 25% by volume. In the following, reasonable alternatives to these choices are discussed to assess the robustness of our conclusions.

Diffusion coefficients

As the creep rate critically depends on the diffusion coefficients, we have re-calculated the creep rates of periclase using another dataset of diffusion coefficients for oxygen in periclase, those of Yoo et al.⁴⁰, which are faster than those of Ita and Cohen¹⁴ and those of Yang and Flynn¹⁵. Fitting equation (1) on the data of Yoo et al.⁴⁰ (Extended Data Fig. 1), we determine two new constants for equation (1): E = 7200 × 10⁻²¹ J and ({S}_{0}=2{k}_{{rm{B}}}). The activation enthalpies obtained with these new constants are given in Extended Data Table 1 and the evolution of the two sets of diffusion coefficients are compared in Extended Data Fig. 2. The creep rates of periclase with the diffusion coefficients of Yoo et al.⁴⁰ are then slightly faster than with the diffusion coefficients of Ita and Cohen¹⁴, but still significantly lower than those of bridgmanite (Extended Data Fig. 3). For this comparison, we have also taken into account the uncertainty on the rheology of bridgmanite owing to the lack of constraints on diffusion coefficients (this is discussed in detail by Reali et al.²¹). For this purpose, we vary the concentration of vacancies in bridgmanite X_v between 10⁻³ and 10⁻⁶. Extended Data Fig. 2 shows that no matter which combination of data we may choose, Mg (or Si) diffusion in bridgmanite is always significantly faster than oxygen diffusion in periclase. The contrast is particularly marked at high pressures. Extended Data Figs. 3–7 show that even taking into account these uncertainties in the most conservative way does not change any of our conclusions.

Influence of iron

Our study is based on the end-member compositions MgO and MgSiO₃, whereas in the mantle these minerals contain iron (and aluminium for bridgmanite). These differences in composition are likely to affect the values of the diffusion coefficients on which our calculations are based. In bridgmanite, this is taken into account by the range of vacancy concentrations considered²¹, see discussion above and Extended Data Figs. 3, 4, 6 and 7. In ferropericlase, it has been shown that oxygen diffusion coefficients are unaffected by the presence of trivalent cations^13,41, because these do not alter the concentration of neutral cation–anion vacancy pairs that are responsible for oxygen transport. The addition of significant amounts of FeO, however, is likely to alter oxygen diffusion coefficients. This effect was considered by Reali et al.⁴², who showed that anion diffusion in MgO and in other oxides and halides with rock salt structures is well described by a homologous temperature scaling, wherein the diffusion coefficients as a function of temperature collapse onto a common curve when normalized to the melting temperature. The addition of FeO lowers the melting temperature of periclase, and hence is expected to enhance oxygen diffusivity. We have thus calculated oxygen diffusion coefficients under lower mantle conditions for ferropericlase as a function of its Fe content, under the same assumptions made by Reali et al.⁴². Considering the addition of up to 20 mol% FeO in periclase⁴³, we do find that the enhancement of oxygen diffusion is up by a factor of approximately 2000 relative to pure MgO. However, this enhancement does not change the conclusions of this study for the bulk lower mantle, as oxygen diffusion in ferropericlase remains significantly slower than Mg or Si in bridgmanite.

Influence of periclase volume fraction

We have computed the rheology of bridgmanite–periclase aggregates for various volume fractions of periclase based on four homogenization procedures (Extended Data Figs. 6 and 7), the difference between them being related to the knowledge on the specimen microstructure that is taken into account.

First, Reuss and Voigt bounds are valid for any microstructures. They only depend on the volume fraction of both phases, not on their geometrical arrangement. They provide rigorous lower and upper bounds for ({widetilde{sigma }}_{0}), respectively. In other words, ({widetilde{sigma }}_{0}) in the aggregate cannot be smaller than ({{widetilde{sigma }}_{0}}^{{rm{Reuss}}}) nor larger than ({{widetilde{sigma }}_{0}}^{{rm{Voigt}}}). The Reuss bound is obtained when assuming that both phases are submitted to the same stress, whereas the Voigt bound is obtained when assuming a uniform strain rate in both phases.

The Hashin-Strikman variational upper bound²⁹, denoted HS+ in Extended Data Figs. 6 and 7, is another rigorous upper bound for ({widetilde{sigma }}_{0}). It is more stringent than the Voigt upper bound as is based on an extra assumption concerning the specimen microstructure, that is, that the microstructure is isotropic. It is obtained when assuming that the softer phase (here, bridgmanite) is spread as inclusions in a stiffer matrix (periclase) while the associated nonlinear homogenization problem is addressed using the variational procedure introduced by Ponte Castaneda, also known as ‘modified secant’ method⁴⁴.

To stick more closely to the possible rheology of an assemblage representative of the lower mantle, one can reasonably assume that bridgmanite and periclase phases are randomly mixed, that is, that contrary to what is assumed within the HS+ bound, both phases are on the same footing, with neither of them playing the role of a matrix and the other one being inclusions. In that case, a good estimation of the effective behaviour is provided by the POSO extension of the Self-Consistent scheme (SC in Extended Data Figs. 6 and 7), introduced in the previous section. Here, unlike previous applications to polycrystalline aggregates, where individual crystal orientations were taken into account^45,46, a two-phase aggregate was considered. A fully optimized linearization procedure for the SC scheme has been proposed recently^47,48 and applied to minerals of the mantle transition zone⁴⁹; results are close to the POSO version when the mechanical contrast between the phases and the nonlinearity of their behaviour is not loo large, as here.

Typical results are shown in Extended Data Figs. 6 and 7, calculated here for pressures between 30 and 120 GPa and X_v = 10⁻³–10⁻⁶. Due to the large mechanical contrast between periclase and bridgmanite, the effective rheology of the aggregate, predicted by the POSO-SC model, lies very close to the Reuss lower bound for volume fractions lower than 30%, whereas it is close to the Voigt upper bound for a volume fraction larger than 70%. In between, the effective stress ({widetilde{sigma }}_{0}) evolves significantly, by more than one order of magnitude. This would correspond to a change in viscosity of more than three orders of magnitude, the viscosity being proportional to ({left({widetilde{sigma }}_{0}right)}^{n}). One other specific feature of these results is illustrated in the right panels of Extended Data Figs. 6 and 7, which show the mean equivalent strain rate in the two phases within the aggregate, normalized by the effective strain rate, as predicted by the SC model. It is obtained that the periclase strain rate is really small at volume fractions lower than 40%. This volume fraction corresponds to a ‘mechanical percolation threshold’ observed in previous work^50,51. At small volume fractions of the highly viscous periclase phase, the aggregate can deform with almost no plastic strain in that phase.