Periclase deforms more slowly than bridgmanite under mantle conditions

Transport of heat from the interior of the Earth drives convection in the mantle, which involves the deformation of solid rocks over billions of years. The lower mantle of the Earth is mostly composed of iron-bearing bridgmanite MgSiO3 and approximately 25% volume periclase MgO (also with some iron). It is commonly accepted that ferropericlase is weaker than bridgmanite1. Considerable progress has been made in recent years to study assemblages representative of the lower mantle under the relevant pressure and temperature conditions2,3. However, the natural strain rates are 8 to 10 orders of magnitude lower than in the laboratory, and are still inaccessible to us. Once the deformation mechanisms of rocks and their constituent minerals have been identified, it is possible to overcome this limitation thanks to multiscale numerical modelling, and to determine rheological properties for inaccessible strain rates. In this work we use 2.5-dimensional dislocation dynamics to model the low-stress creep of MgO periclase at lower mantle pressures and temperatures. We show that periclase deforms very slowly under these conditions, in particular, much more slowly than bridgmanite deforming by pure climb creep. This is due to slow diffusion of oxygen in periclase under pressure. In the assemblage, this secondary phase hardly participates in the deformation, so that the rheology of the lower mantle is very well described by that of bridgmanite. Our results show that drastic changes in deformation mechanisms can occur as a function of the strain rate.

Transport of heat from the interior of the Earth drives convection in the mantle, which involves the deformation of solid rocks over billions of years. The lower mantle of the Earth is mostly composed of iron-bearing bridgmanite MgSiO 3 and approximately 25% volume periclase MgO (also with some iron). It is commonly accepted that ferropericlase is weaker than bridgmanite 1 . Considerable progress has been made in recent years to study assemblages representative of the lower mantle under the relevant pressure and temperature conditions 2,3 . However, the natural strain rates are 8 to 10 orders of magnitude lower than in the laboratory, and are still inaccessible to us. Once the deformation mechanisms of rocks and their constituent minerals have been identified, it is possible to overcome this limitation thanks to multiscale numerical modelling, and to determine rheological properties for inaccessible strain rates. In this work we use 2.5-dimensional dislocation dynamics to model the low-stress creep of MgO periclase at lower mantle pressures and temperatures. We show that periclase deforms very slowly under these conditions, in particular, much more slowly than bridgmanite deforming by pure climb creep. This is due to slow diffusion of oxygen in periclase under pressure. In the assemblage, this secondary phase hardly participates in the deformation, so that the rheology of the lower mantle is very well described by that of bridgmanite. Our results show that drastic changes in deformation mechanisms can occur as a function of the strain rate.
This work aims to model steady-state creep of periclase (MgO) under lower mantle conditions, with an emphasis on the role of extremely low strain rates. It is based on the extensive knowledge we have of the deformation of this oxide following (1) the large number of experimental studies carried out on this ceramic since the 1970s, including the more recent mineral physics studies taking into account the effect of pressure, and (2) the important modelling effort carried out in recent years. At room pressure, periclase is easily deformed by glide of < 110> 1 2 dislocations on {110} planes. The thermal activation of this mechanism can now be modelled from the atomic scale 4 , even at low strain rates 5 . Pressure strongly affects this mechanism and induces a transition of glide towards the {100} planes above 50 GPa [6][7][8][9] . At high temperature, dislocations glide without friction (athermal regime) and the flow stress results from dislocation interactions, leading, at large strains and strain rates, to a pronounced hardening 10 . At low strain rates, ionic diffusion can enhance recovery mechanisms and allow steady-state creep 11 . In this study, we model the influence of mantle pressures on steady-state creep of periclase by 2.5-dimensional (2.5D) dislocation dynamics (DD). We focus on four pressures: 30, 60, 90 and 120 GPa, which, along the geotherm 12 , correspond to 2,000, 2,300, 2,500 and 2,800 K respectively, at around depths of 800, 1,500, 2,100 and 2,700 km respectively.

Diffusion coefficients
Dislocation climb in MgO is enabled by diffusion of Mg and O atoms, and the climb rate is controlled by the slower diffusing species. Diffusion in MgO has been studied for decades, and the mechanisms and rates of cation and anion diffusion are well established 13 . Mg and O diffuse using vacant sites on the cation and anion sublattices, respectively. Because the concentration of cation vacancies is much higher than that of anion vacancies, in both nominally pure synthetic and natural periclase samples, O diffusion is orders of magnitude slower than Mg diffusion, and is the rate-limiting step in dislocation climb. Whereas Mg diffusion is controlled by extrinsic vacancies, produced to charge-balance trivalent ( where D 0 = 7.2 × 10 −4 m 2 s −1 and ΔH sd are, respectively, the exponential prefactor and the activation enthalpy 11 . The activation enthalpies for each set of pressure/temperature and the diffusion coefficients are given in Table 2

Modelling dislocation creep of periclase
To model plasticity of periclase under mantle conditions we use a 2.5D DD simulation approach already benchmarked for periclase at ambient pressure 11 . DD simulation is a modelling tool that describes the collective motion and interactions of dislocations at the mesoscale. This technique is based on continuum elasticity theory, which provides the description of the elastic field induced by dislocations in a crystal, their interactions with each other and with respect to the stress field resulting from external loading. Within our DD code, dislocation glide is coupled with climb to investigate high-temperature creep (Methods).
Our creep results are presented in Fig. 2a, which shows the creep strain rates obtained from the 2.5D DD simulations as a function of the creep stress for the four (P, T) conditions selected. A constant slope is found for all these strain rate-stress curves, which gives the value of the power law stress exponent, n = 3.1. A strong pressure dependence is found between (30 GPa, 2,000 K) and (90 GPa, 2,500 K), whereas creep rates at (90 GPa, 2,500 K) and (120 GPa, 2,800 K) are very close due to very similar values for the diffusion coefficients ( Fig. 1 and Table 2). On Fig. 2a we compare the creep rates of periclase with those of bridgmanite at the same depths. For bridgmanite, we use the results of Reali et al. 21 for pure climb creep (PCC) with a vacancy concentration X v = 1 × 10 −5 ; these results were shown to be in good agreement with empirical geophysical constraints on the viscosity profile of Earth's lower mantle. Bridgmanite and periclase exhibit similar stress exponents. However, we see that these calculations lead to bridgmanite having higher creep rates than periclase at any depth. To better illustrate this difference, we report, on Fig. 2b, the viscosity contrast between those two phases, taken as the ratio η η / Per Bdm . We see that under lower mantle conditions bridgmanite creep rates are 10 3 to 10 7   and 120 GPa respectively. The range of temperatures considered in this study [2,000 K:2,800 K] is shown with the grey shading. The symbols (filled circles) correspond to the (P, T) conditions along the geotherm considered in the present study 12 (Table 2). greater than periclase creep rates. This result reflects the fact that in the lower mantle, anionic diffusion in periclase is much slower than cationic diffusion in bridgmanite (Extended Data Fig. 2). In the online Methods section and Extended Data Figs. 1-4 we consider other diffusion data from the literature for periclase and other vacancy concentrations in bridgmanite and show that the conclusions stated above remain valid in any case.

Laboratory versus nature
In this paragraph we address the influence of strain rate on deformation mechanisms. Our calculations showing that bridgmanite deforms faster than periclase under mantle conditions seem to contradict the common intuition that the silicate is stronger than the oxide. Beyond general considerations 1 or studies of analogous compounds 22,23 , Girard et al. 3 conducted experiments under mantle conditions of pressure and temperature that show very clearly that ferropericlase accommodates most of the strain compared to bridgmanite. We show here that there is no fundamental opposition between modelling and experiments and that our models are able to account equally for natural and experimental conditions. The experiments of Girard et al. 3 are conducted at a strain rate of 3 × 10 −5 s −1 and under very high stresses. Under these laboratory conditions, dislocation glide is expected as the dominant deformation mode. For bridgmanite, dislocation glide is thermally activated under these conditions and a model based on the Orowan equation can be built 24 . The stress values predicted by this model are reported on Fig. 3. They are close to the experimental values (slightly lower as this model was constructed for single crystals and does not involve grain boundaries). At 30 GPa and 2,000 K, dislocations in periclase glide in the athermal regime 6 . Several deformation experiments of periclase report that in the athermal regime, no steady state is reached and that strain hardening occurs [25][26][27] under the influence of dislocation-dislocation interactions. The origin of this strain hardening has been studied and modelled by Amodeo et al. 10 . It is well described by the Taylor mechanism, and the athermal stress τ f required to overcome the forest dislocation network is given by: where μ is the shear modulus, b the modulus of the Burgers vector, β is a forest strengthening coefficient calibrated by Amodeo et al. 10 , ρ f is the forest dislocation density and ρ ref is the reference dislocation density. Transmission electron microscopy observations 28 on the specimens of Girard et al. 3 confirm this view and show a very large dislocation density estimated at over 10 15 m −2 . Using this dislocation density and equation (3) leads to an excellent description of experimental values for the deformation of ferropericlase under shallow lower mantle P, T conditions at laboratory strain rates (Fig. 3). Our modelling approach allows us to reproduce very well the experimental data of Girard et al. 3 3 (symbols) can be compared to numerical models (lines). At laboratory stresses and strain rates numerical models based on dislocation glide are in fair agreement with experimental data and show that bridgmanite (Bdm) is stronger than periclase (Per). At mantle strain rates, creep (Weertman creep modelled by 2.5D DD in periclase (see text) and pure climb creep (PCC) in bridgmanite) is controlled by ionic diffusion, which leads to an inversion of rheology, with periclase deforming more slowly than bridgmanite. deformation conditions (that is, smaller strains or more recovery), the rheology contrast under laboratory conditions should be even larger. At 30 GPa, oxygen diffusion is already slow, and strain rates much lower than those obtainable in laboratory experiments would be required for climb to become active and promote the recovery necessary for steady-state creep.

Deformation of a bridgmanite-periclase aggregate
Our model places first-order constraints on the intrinsic properties of the major phases of the lower mantle. The rheological behaviour of rocks in the lower mantle certainly involves more complexity. A first step towards this direction is to use dedicated micromechanical models, to connect the viscoplastic behaviour of a polycrystalline aggregate representative of a lower mantle assemblage with those of the individual phases, bridgmanite and periclase, and, at the same time, quantify the contribution of both phases to the effective strain rate, in relation with the aggregate microstructure.
We have therefore applied a state-of-the-art scale transition homogenization method 29 (partially optimized second-order self-consistent procedure, denoted POSO-SC) adapted for multiphase aggregates with nonlinear rheologies, such as those considered here. The model considers an isotropic microstructure, that is, a random mixture of bridgmanite and periclase with no preferred elongated shape nor crystallographic texture for each phase. Details are given in the Methods. The effective viscoplastic behaviour can be expressed with a nonlinear creep law with a stress sensitivity n = 3.1 (see Methods). Figure 4 shows the results obtained for a vacancy concentration X v = 1 × 10 −5 in bridgmanite and a volume fraction of periclase of 25% (alternative figures considering different combinations of diffusion coefficients are presented in Extended Data Fig. 5). For pressures in the range 30-120 GPa, the effective rheology (blue line) lies very close to that of pure bridgmanite (red line), with the flow stress values differing by a factor of only around 1.65, whereas the flow stress of pure periclase (black line) is one to two orders of magnitude larger. Periclase deforms at least 100 times slower than bridgmanite in the aggregate and its contribution to the effective strain rate is negligible. In the online Methods section, we show that the volume fraction of periclase would have to exceed 40% to significantly affect the rheology of the assemblage. The presence of iron in MgO may affect anionic diffusivity but does not change our conclusions for the bulk lower mantle (see online Methods section).
Following the widespread acceptance that ferropericlase is the weakest phase 1 , and inspired by the experimental results of Girard et al. 3 , several studies 30-32 have discussed the implication of the morphology of ferropericlase (either forming a load-bearing framework or an interconnected weak layer 33 ) on the rheology of the lower mantle. Our study shows that under the conditions (mainly strain rate) of the lower mantle, the behaviour of periclase is radically different and that this fundamental assumption is now challenged. Periclase, minor in volume, hardly deforms and probably induces little evolution of the microstructure. Consequently, the rheology of the lower mantle is very well described by that of bridgmanite deforming by pure climb creep.

Online content
Any methods, additional references, Nature Portfolio reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41586-022-05410-9. Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

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 12 at 30, 60, 90 and 120 GPa, and use the periclase elastic properties from Karki et al. 34 (Extended Data Table 2).
To capture the statistically representative volume element of the microstructure, dislocations of Burgers vector b = 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 11 , 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. 11 .
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): v bτ B = g 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 35 and use as a starting point the experimental value recorded at 300 K 36 . 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 35 : where τ c is the climb stress, calculated similarly to τ but resolved here along the climb direction. α π R r = 2 /ln( / ) c is a geometrical factor that describes the cylindrical geometry of the vacancy flux field around the dislocation line, where R and r 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 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 / c ratio does not significantly affect the climb velocity values and here is taken as constant and equal to 100 37 . Ω is the formation volume of vacancies in the Mg and O sites and is calculated from the unit cell volume of MgO-that is, Ω a Z = / 3 where Z = 4 is the number of formula units per unit cell. c ∞ and c 0 are the vacancy concentrations far from the dislocation and the equilibrium vacancy concentration in the bulk volume, respectively. Following Boioli et al. 37 , 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 0 ). 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: ith ε 0 a reference strain rate, σ 0 a reference stress, n the stress sensitivity, ε eq the equivalent strain rate and σ eq the equivalent stress, defined as: with summation over repeated indices (Einstein convention), ε ij and σ′ ij components of the strain rate tensor and deviatoric stress tensor, respectively. We chose here to take ε = 1 × 10 s 0 −15 −1 , that is, a value typical for lower mantle convection. Doing so, the σ 0 values are of the order of magnitude of the equivalent stress encountered in situ during mantle flow. Values of σ 0 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 ε̇ and σ are the effective strain rate and stress tensors, respectively, given by the volume average (denoted • ⟨ ⟩) of local strain rate and stress fields: ε ε = ⟨ ⟩ and σ σ = ⟨ ⟩. ∼ 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 ∼ n n = =3.1. Finally, ∼ σ 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 29 .

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 3 bridgmanite with a vacancy concentration of 10 −5 . 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. 40 , which are faster than those of Ita and Cohen 14 and those of Yang and Flynn 15 . Fitting equation (1) on the data of Yoo et al. 40 (Extended Data Fig. 1), we determine two new constants for equation (1): E = 7200 × 10 −21 J and S k = 2 0 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. 40 are then slightly faster than with the diffusion coefficients of Ita and Cohen 14 , 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. 21 ). For this purpose, we vary the concentration of vacancies in bridgmanite X v between 10 −3 and 10 −6 . 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 3 , 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 21 , 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. 42 , 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. 42 . Considering the addition of up to 20 mol% FeO in periclase 43 , 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 σ 0 ∼ , respectively. In other words, σ 0 ∼ in the aggregate cannot be smaller than ∼ σ 0 Reuss nor larger than ∼ σ 0 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 29 , denoted HS+ in Extended Data Figs. 6 and 7, is another rigorous upper bound for σ 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 44 .
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 49 ; 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 −3 -10 −6 . 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 ∼ σ 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 σ ( ) n 0 ∼ . 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.

Data availability
The data of this manuscript are available at https://doi.org/10.5281/ zenodo.5970733. Source data are provided with this paper.