Ferropericlase Control of Lower Mantle Rheology: Impact of Phase Morphology

The rheological properties of Earth's lower mantle play a key role for global mantle dynamics. The mineralogy of the lower mantle can be approximated as a bridgmanite‐ferropericlase mixture. Previous work has suggested that the deformation of this mixture might be dramatically affected by the large differences in viscosity between bridgmanite and ferropericlase. Here, we employ numerical models to establish a connection between ferropericlase morphology and the effective rheology of the Earth's lower mantle using a numerical‐statistical approach. Using this approach, we link the statistical properties of the two‐phase composite to its effective viscosity tensor using analytical approximations. We find that ferropericlase develops elongated structures within the bridgmanite matrix that result in significantly lowered viscosity. While our findings confirm previous endmember models that suggested a change of mantle viscosity due to the formation of interconnected weak layers, we show that significant rheological weakening can thus be already achieved even when ferropericlase does not form an interconnected network. Additionally, the alignment of weak ferropericlase leads to a pronounced viscous anisotropy that develops with total strain, which may have implications for understanding the viscosity structure of Earth's lower mantle as well as for modeling the behavior of subducting slabs. We show that to capture the effect of ferropericlase elongation on the effective viscosity tensor (and its anisotropy) in large‐scale mantle convection models, the analytical approximations that have been derived to describe the evolution of the effective viscosity of a two‐phase medium with aligned elliptical inclusions can be used.


Introduction
The Earth's lower mantle links the hot, vigorously convecting outer core to the upper mantle and ultimately the Earth's surface. It therefore has a strong impact on the dynamics of outer core and mantle convection and consequently also on Earth's thermal evolution. Despite its importance, the rheological properties of the Earth's lower mantle are not very well constrained. Seismic tomography provides a glimpse into the present-day structure of the Earth's mantle. Although the resolution becomes worse with depth, these tomographies have revealed a number of large-and smaller-scale heterogeneities in the lower mantle (e.g., Garnero & McNamara, 2008;Hedlin et al., 1997;Kaneshima & Helffrich, 2010). In particular, it has been shown that only a fraction of subducting slabs directly sink into the lower mantle but that they rather stagnate at 660-km depth (Fukao et al., 2009). This stagnation is most likely related to the phase transition of ringwoodite to bridgmanite and ferropericlase. Additionally, a large number of slabs seems to stagnate at 1,000-km depth (e.g., Fukao et al., 2009;Goes et al., 2017). The cause of this stagnation has been attributed to a viscosity increase occurring at this depth (Marquardt & Miyagi, 2015;Rudolph et al., 2015).
In a lower mantle of pyrolitic composition, bridgmanite is expected to compose the largest part of the lower mantle (with estimates ranging from approximately 70-90%), followed by ferropericlase and calcium silicate perovskite (e.g., Hirose et al., 2017;Stixrude & Lithgow-Bertelloni, 2012). It has been suggested that the rheology of the Earth's lower mantle is governed by the rheology of its two main constituents bridgmanite and ferropericlase; however, both rheologies are not very well known at lower mantle conditions and are subject of active research (e.g., Girard et al., 2016;Immoor et al., 2018;Kaercher et al., 2016;Merkel et al., 2003;Miyagi & Wenk, 2016;Reali et al., 2019). The rheology of calcium silicate perovskite is even less known due to experimental difficulties (Miyagi et al., 2009;Nestola et al., 2018). It has been argued that due to the absence of seismic anisotropy, the lower mantle should be deforming in diffusion creep (Karato & Li, 1992). The viscosity of bridgmanite is estimated to be significantly larger than that of ferropericlase, thus resulting in significant viscosity variations depending on the quantity and degree of interconnectivity of In this study, we extend previous studies by employing a combined numerical-statistical approach to establish a connection between the phase morphology of ferropericlase and the effective rheology of the bridgmanite-ferropericlase mixtureby using statistical descriptors. We then show how a randomly distributed weak phase affects the effective viscosity even without interconnection and that elongation may result in significant weakening. We then compare the numerical results to analytical approximations for two-phase composites, which are used in conjunction with the statistical description of the weak phase morphology. Finally, we study the formation of elongated structures due to simple shear deformation and their corresponding effective rheology using numerical models as well as analytical approximations.

Methodology
As mentioned above, we follow a statistical-numerical approach to investigate the effective rheology of bridgmanite-ferropericlase mixtures. This approach consists of the following steps: (1) characterization of ferropericlase morphology, (2) artificial sample generation for direct numerical modeling, (3) direct numerical modeling, and (4) the prediction of numerical results via an analytical approximation. In the following sections, each step will be described in detail. In Section 2.1, we show how a random heterogeneous two-phase medium can be characterized in a statistical manner and how key statistical properties can be extracted from this description. In Section 2.2, we explain how these statistical properties can be used to create artificial samples that will then be later used for direct numerical modeling as described in Section 2.3. Finally, in Section 2.4, we shortly discuss different analytical approximations that have been proposed to describe the effective rheology of two-phase media.

Characterization of Ferropericlase Morphology
The microstructure of bridgmanite-ferropericlase mixtures has been investigated by several researchers. Yamazaki et al. (2009) synthesized (Mg,Fe)SiO 3 perovskite and ferropericlase under hydrostatic conditions. The observed texture shows that ferropericlase does not take a specific shape, but is rather randomly distributed (see Figure 1). As the sample has been synthesized under hydrostatic conditions, ferropericlase and bridgmanite grains exhibit rather isotropic shapes with no visible shape preferred orientation. On the other hand, the grains in the deformed sample in Girard et al. (2016) exhibit an elongated structure. In both cases, the observed ferropericlase structures are rather randomly distributed. It is therefore necessary to characterize these structures using more general statistical descriptors.
Two-phase media can be characterized by a number of statistical descriptors (e.g., Torquato, 2002), some of which will be described below. A two-phase medium is defined by the following function (e.g., Jiao & Chawla, 2014): Several questions arise when the microstructure of a two-phase medium is investigated. One of these questions is: Choosing an arbitrary point x 1 within this medium, what is the probability of finding phase i at this point? The answer to this question is given by the so-called 1-point probability function S (i) 1 (x), which describes the probability of finding phase i at any given point x and is given by the average of I (i) (x): which is equivalent to the volume fraction. A second question is: What is the probability of finding two arbitrary points x 1 and x 2 within the same phase i? The 2-point probability function S (i) 2 (x 1 , x 2 ) gives the probability of finding phase i at both points x 1 and x 2 : The 2-point probability function is also known as the autocorrelation function (e.g., Torquato, 2002) and has been used to analyze the fabrics of geological materials, in particular to determine the grain size and shape preferred orientation of thin sections (e.g., Heilbronner, 1992).
While the 2-point probability function does contain some information about the morphology of phase i, it does not contain information on the connectivity between different clusters of this phase, where a cluster is defined as a region of phase i that can be reached from a point in that phase without having to go through another phase (Torquato, 2002). This information is contained in a subset of the 2-point probability function, the 2-point cluster function C (i) 2 (x 1 , x 2 ) (Jiao et al., 2009;Torquato et al., 1998), which answers the question: What is the probability of finding two arbitrary points x 1 and x 2 within the same cluster of phase i? As we only investigate the effect of a single additional phase (ferropericlase) on the effective rheology in this study, we will drop the index i in the following and refer to the 2-point probability and the 2-point cluster functions as S 2 and C 2 , which will be used to describe the morphology of ferropericlase inclusions.
For a statistically homogeneous and isotropic medium, both the 2-point probability and the 2-point cluster function only depend on the distance r = |r| = |x 2 −x 1 | between the points x 1 and x 2 . However, if a medium is anisotropic, both probability functions are not only a function of the length r but also of the orientation of r. For this reason, 2-point correlation functions are often computed in multiple directions (e.g., Jiao & Chawla, 2014;Sheehan & Torquato, 2000).
Here, we compute the 2-point probability and cluster functions not only for a certain set of directions but obtain a complete two-dimensional image of S 2 and C 2 . For the two samples presented in Yamazaki et al. (2009) and Girard et al. (2016) these functions are shown in Figure 1 together with the original scanning electron microscope backscatter images. The isotropic nature of the undeformed sample is directly visible in the circular shape of the 2-point cluster function, whereas the anisotropic structure of the deformed sample is reflected by the tilted ellipsoidal shape of the 2-point cluster correlation function.
As the 2-point cluster function contains significantly more information than the 2-point probability function, it is generally considered to be superior to other 2-point statistical descriptors (Jiao et al., 2009). The 2-point cluster functions have been frequently employed to determine the effective mechanical properties of two-phase composites in order to generate materials with desired properties (e.g., Torquato, 2010). Additionally, these correlation functions have not only been used to quantify mechanical but also hydraulic and electrical properties (e.g., Torquato, 2002). In this study, we aim to connect the properties of the 2-point cluster function of the ferropericlase inclusions to the rheological properties of the ferropericlase-bridgmanite mixture.
The 2-point probability and 2-point cluster functions can be used to statistically determine the morphology (size of inclusions and their anisotropy) of the random heterogeneous two-phase mixture given by the bridgmanite/ferropericlase mixtures. As we can see in Figure 1, both probability functions capture the isotropic and anisotropic nature of the two-phase mixture. To quantify the degree of anisotropy, we can now compute the correlation lengths c * i of the two probability functions, which are defined as the distance over which C (i) 2 (x 1 , x 2 ) (or S (i) 2 (x 1 , x 2 )) decreases by 1∕e (e.g., Taconet & Ciarletti, 2007). A similar procedure was employed in Heilbronner (1992) to determine the grain size in thin sections. If the fabric is isotropic, the 1∕e-isoline exhibits a circular shape (see gray line in inset in the top right image of Figure 1) and the correlation length can be obtained by fitting a circle to this isoline to obtain a single value for c * . In the anisotropic case, the 1∕e-isoline exhibits an elliptical shape, which can be fitted with an ellipse whose semimajor and semiminor axes correspond the the correlation lengths c * x and c * . Additionally, the angle of the semimajor axes with the horizontal yields information about the average orientation of ferropericlase inclusions.

Artificial Generation of Heterogeneous Random Media
As can be seen in Figure 1, the distribution of ferropericlase in the ferropericlase-bridgmanite composite is rather random. This calls for a statistical approach to determine the effect of weak phase elongation on the effective viscosity. We therefore have to create a large number of artificial samples with prescribed statistical properties. The distribution of the two phases is set via the following procedure (Cahn, 1965;Roberts & Teubner, 1995;Vel & Goupee, 2010): For an artificial sample with the desired statistical properties (given by the correlation lengths c x , c and , which are input parameters), a correlated random surface is generated by convoluting an uncorrelated Gaussian distribution of random numbers z u (x, ) with a two-dimensional Gaussian autocovariance function F(x, ): where z u (x, ) and F(x, ) are given by The term in the exponent in equation (6) is the general equation of a tilted ellipse with its center at the origin. Consequently, the factors e 1 , e 2 , and e 3 can be computed from the semimajor and semiminor axes (given here by input parameters c x and c ) as well as the inclination angle : THIELMANN ET AL. The prefactor A in equation (4) is then chosen such that the values of z(x, ) lie in the interval [0, 1]. The random surface is then converted into a binary distribution by choosing a threshold level p below which z(x, ) corresponds to Phase 2, otherwise to Phase 1. p can then be iteratively adapted until the desired fraction of Phase 2 is reached (see Figure 2 for an example of a random surface and its conversion into binary distributions).
The main purpose of this study consists of linking the topology of media with different amounts of layering to their effective rheology. The layering (i.e., anisotropy) can be controlled by varying the value of c x . However, it turns out that simply increasing c x results in an undesired increase of weak phase cluster size and a decrease in cluster number (see Figure 3). A remedy to this issue consists of additionally rescaling c as follows: where c x,iso = c ,iso are the correlation lengths for the isotropic case where the product of c x,iso c ,iso equals the product c x c . Due to this rescaling, the number of clusters as well as their area is preserved in the resulting anisotropic random structure ( Figure 3).

Distribution Properties
To assess the effect of the binarization process on the statistical properties of the artificially generated random media, we computed the effective correlation lengths c * x and c * for samples with given correlation lengths ranging from c x ∈ [0.05, 0.5]. c was computed according to equation (8). For each value of c x , we created 2,500 samples and computed the effective correlation lengths c * x and c * as described in Section 2.1. Not surprisingly, our procedure to generate heterogeneous media with given correlation lengths c x and c produces binary distributions with a range of values for c * x and c * , as can be seen in Figure 4. As the size of the weak phase clusters depends on their desired volume fraction, we do not expect c * x and c * to equal c x and c . However, it turns out that c * x is a linear function of c x and that the scaling relationship (8) is not affected by the binarization process. The effective correlation length c * is therefore given by where c * x,iso = c * ,iso is the effective correlation length in the corresponding isotropic case. x . Values for c * were obtained in the same way as for the plot above. The solid line denotes a fit to the data using a function of the same form as (8). The obtained constant approximately corresponds to the product of c * x,iso c * ,iso , which shows that the scaling relationship (8) between both correlation lengths is not affected by the binarization process.
As weak phase clusters become increasingly elongated, they may at some point span the domain. The formation of such spanning clusters is very important in two-phase media, as it has a strong impact on mechanical properties (e.g., Roberts & Garboczi, 2002). However, in the samples shown here, the occurrence of a spanning cluster may be related to sample size and correlation lengths of the clusters and as such exhibit a finite size effect. Consequently, it is important to distinguish samples with spanning clusters from those without spanning clusters, as they may exhibit different mechanical properties. For large elongations, the 1∕e contour line of the cluster correlation function C (i) 2 (x 1 , x 2 ) may also not form an ellipse but span the entire domain and a determination of the effective correlation lengths is not possible.

Direct Numerical Modeling 2.3.1. Governing Equations
In its simplest form, the creeping flow of mantle rocks can be described with the incompressible Stokes equations: where v denotes the velocity vector, P pressure, and the viscosity (the forcing term due to gravity has been neglected here). Recent studies (Girard et al., 2016;Nzogang et al., 2018) indicate that the deformation of bridgmanite and ferropericlase may be accomodated via dislocation activity and thus potentially giving rise to nonlinear behavior. However, given the uncertainties related to the extrapolation of laboratory data to the Earth, it is unclear which deformation mechanism is active in the lower mantle. In this study, we thus focus on linear viscous rheologies for simplicity. In the case of an isotropic, linear viscous material, i can be related to the strain rate tensor where . i is given by We can now nondimensionalize these equations using characteristic scales for viscosity, length, and time which are given by where we choose the characteristic scales as the viscosity of the host material 1 , the size of the model domain L, and the background strain rate . BG . This procedure significantly reduces the number of free parameters, as host viscosity and background strain rate take the nondimensional value of 1. The computed effective viscosities from these models are then also normalized by the host viscosity.

Numerical Solution of the Governing Equations
The system of equations given in (10) is solved in primitive variable formulation using the hybrid finite element particle-in-cell method by employing an optimized version of the open-source code MVEP2 (e.g., Thielmann et al., 2014). The model domain consists of a 1 × 1 quadratic box. In order to properly resolve the given structures, the model domain is discretized using a regular mesh of 1, 024 × 1, 024 Q 1 P 0 elements (see also appendix A). For advection we employ passive markers using a Runge-Kutta fourth-order scheme. In time-dependent simulations, we employ a constant nondimensional time step of 10 −4 , which implies that the system experiences the same amount of overall strain per time step. We use initially 16 markers per cell to track the different phases, which amounts to ≈16.8 million markers. If the number of markers per cell falls below six markers, 10 new markers are inserted with their properties being determined by the nearest neighboring marker. Material properties are interpolated from markers to nodes via harmonic interpolation to the integration points. To obtain the different components of the viscosity tensor, two deformation steps are performed, one in pure shear and one in simple shear. Coordinate updates (in time-dependent simulations) are only done after a simple shear step with a constant nondimensional time step Δt = 5 · 10 −4 .

Postprocessing
As already discussed above, even though we assume the viscosity of the individual phases to be isotropic, the effective viscosity of the two-phase medium can be anisotropic due to the morphology of the weak phase, which implies that viscosity is not a scalar, but rather a tensor. To determine the different components of the viscosity tensor i kl , we extract the spatial averages of the stress tensor ⟨ i ⟩ and ⟨ . i ⟩ from the numerical simulations. Using equation (21), we then compute the respective components of the viscosity tensor, keeping in mind that ⟨ Simple shear boundary conditions are employed by setting the horizontal velocity of the top boundary condition to a nondimensional value of 2 and its vertical velocity to 0 to obtain a nondimensional strain rate of 1. Bottom boundary conditions are zero velocity for both horizontal and vertical velocity components. Side boundaries are periodic. In the case of pure shear deformation, vertical velocities at top and bottom boundaries are set to −1 and 1, respectively. Finally, we obtain the principal components n and s of the viscosity tensor and its orientation through principal component analysis.
In addition to computing these components directly from the spatial averages of stress and deformation rate tensors, we also characterize each sample using the 2-point cluster correlation function C (i) 2 (x 1 , x 2 ), where we obtain the correlation lengths c * x and c * as well as the inclination by fitting an ellipse to the 1∕e contour line.

Analytical Approximations
The effective rheology of two-phase aggregates has been the subject of many studies in the past. The bounds for the effective viscosity of a two-phase medium with Newtonian rheologies are given by the Voigt and Reuss bounds (e.g., Handy, 1990;Takeda, 1998): where 1,2 are the viscosities of the two phases, is the volume fraction of Phase 2, and Δ = 2 ∕ 1 is the viscosity contrast between Phase 2 and Phase 1. The two bounds represent the upper and lower bounds on the viscosity of a two-phase mixture. The Voigt bound (or the arithmetic mean) describes the volume-averaged effective viscosity of a layered material, which is deformed in pure shear, with the layering being orthogonal to the compression direction. The Reuss bound (harmonic mean) describes the volume-averaged effective viscosity of a layered medium being deformed in simple shear with the layering being parallel to the shear direction (e.g., Schmeling et al., 2008). With increasing viscosity contrast between the two deforming phases, these bounds diverge and predict significantly different values for the effective viscosity. However, both endmember models assume a layered structure of the two phases. This severely limits the applicability of these bounds (e.g., in geodynamical models). As a potential remedy, several researchers have proposed different approaches to overcome this issue (e.g., Dabrowski et al., 2012;Fletcher, 2009;Treagus, 2002;Vel et al., 2016). However, each of these approaches either only deals with very specific inclusion shapes (Dabrowski et al., 2012;Fletcher, 2004;Treagus, 2002) or only considers statistically isotropic two-phase random media (Vel et al., 2016). Here, we combine these previous works to study the effect of weak phase morphology using random heterogeneous media. Based on the observation that the 2-point cluster function of the deformed aggregate from Girard et al. (2016) exhibits an elliptical shape, we turn to analytical solutions that address the effective rheology of a two-phase medium containing aligned elliptical inclusions. In particular, we employ the anisotropic differential effective medium (aDEM) approach presented in Dabrowski et al. (2012) as well as the self-consistent approximation (SCA) developed by Treagus (2003) and Fletcher (2004), which are summarized below.

The aDEM Approach
The aDEM approach extends the differential effective medium approach (Bruggeman, 1935) (based on the assumption of circular inclusions) to systems with aligned elliptical inclusions. Contrary to circular inclusions, the presence of aligned elliptical inclusions gives rise to a nonisotropic viscosity tensor. In plane strain deformation in the x plane, this tensor can be described by a normal viscosity n , a shear viscosity s , and a misalignment angle . Based on the differential effective medium approach, Dabrowski et al. (2012) derived a system of coupled differential equations that can be used to compute the normalized (with the host viscosity 1 ) effective normal and shear viscosity n = n ∕ 1 and s = s ∕ 1 of a composite with aligned elliptical inclusions: where = n ∕ s and p = 0.5 (a∕b + b∕a) are anisotropy factor and shape factor, with a and b being the semimajor and semiminor axis of the elliptical inclusion Δ denotes the viscosity contrast between both phases. The coupled system of equations (16) and (17) can be integrated numerically to obtain the desired values for n and s .

The SCA Approach
Another commonly employed approximation to compute effective medium properties is the SCA (Budiansky, 1965;Hill, 1965). For the same system of aligned ellipses as described in the section above, the first expression to compute n and s based on SCA was proposed by Treagus (2003): where may be either the normal or shear viscosity. To obtain the normal viscosity n , the shape factor p T equals shape factor p mentioned above, whereas it equals p −1 to compute the shear viscosity s .
However, this approximation did not consider the effective host to become anisotropic as soon as an inclusion was added and thus yields too low values for s . This was incorporated by Fletcher (2004), who took into account the host anisotropy. To compute the effective normal and shear viscosities n and s , the following system of nonlinear equations has to be solved for n and (Fletcher, 2004): where p and are the same shape factor and anisotropy factor as described above.
THIELMANN ET AL.

Deformation and Rotation of an Elliptical Viscous Inclusion
The theoretical considerations described above provide estimates for the normalized effective viscosity (given by n and s in a reference frame that is aligned with the principal axes of the aligned inclusions). Generally, the viscosity tensor relating the spatial averages of the deviatoric stress tensor components ⟨ xx ⟩ and ⟨ x ⟩ to the deformation rate tensor components ⟨ . xx ⟩ and ⟨ .
x ⟩ is given by (e.g., Weijermars, 1992) where the components of i kl in the x plane can be computed from the effective normal and shear viscosity n and s as well as the orientation of the major ellipse axes to the horizontal (e.g., Weijermars, 1992): cos 2 sin 2 x x = n sin 2 2 + s cos 2 2 When in addition the weak phase morphology of a random heterogeneous medium is modified during deformation, its viscosity tensor changes continuously. In simple shear, weak phase clusters become more elongated and rotate in the shear direction. The rotation and stretching of ellipses during deformation has been extensively studied (e.g., Bilby & Kolbuszewski, 1977;Mulchrone & Walsh, 2006) and also applied to systems of aligned ellipses (Dabrowski et al., 2012;Fletcher, 2009). The change of aspect ratio r and orientation of an elliptical inclusion is given by where . incl xx , . incl x and W incl x are deformation and vorticity inside the inclusion. These can be related to the far field deformation rates and vorticity as follows (Dabrowski et al., 2012;Fletcher, 2009): where Δ n = incl n = incl n is the ratio between inclusion viscosity and the normal host viscosity (or their respective normalized values).
Comparing the analytical prediction to the results from numerical models, Dabrowski et al. (2012) showed that the analytical solution does predict the evolution of a two-phase medium with initially spherical inclusions up to a shear strain of = 1. At larger strains, predictions and model results start to diverge, the main reason being the interaction between different inclusions which result in nonelliptical inclusion shapes.

Instantaneous Simulations
To determine the impact of weak phase topology, in particular of the elongation in shear direction, we first conducted a series of simulations where we prescribed the degree of layering via the correlation length c x . For the isotropic case, the correlation lengths c x,iso and c ,iso were set to a value of 0.05. The result of two such simulations with a weak phase fraction of = 0.2, a viscosity contrast of Δ = 10 −2 , and different elongations c x = 0.05 and c x = 0.4 (the respective values for c ) is shown in Figure 5, where the distribution Figure 5. Stress fields of a statistically isotropic (left) and anisotropic (right) random heterogeneous medium with a weak phase fraction of = 0.2. The weak phase is 100 times weaker than the surrounding medium. The presence of the weaker phase results in significant stress heterogeneity. Due to the viscosity contrast, stresses in the weaker phase are significantly reduced. It can be seen that the overall stress level is lower in the medium with strongly elongated inclusions. of the effective stress II within both samples is shown. Due to the viscosity contrast between both materials, stresses within the sample are highly heterogeneous, with stresses in the weak phase being significantly lower. Stress levels in the strong phase are generally higher in the isotropic case than in the anisotropic (elongated) case.
For each combination of c x , and Δ , we conducted 2,500 realizations per combination (427,500 simulations in total). We then separated the data set into two parts, one containing the samples without a spanning weak phase cluster and one containing the samples with one or more spanning weak phase cluster(s). For each realization, we also computed the shear and normal viscosities s and n according to (21) as well as the effective correlation lengths c * x and c * . In general, the existence of one or more spanning clusters does not significantly affect the value of n and the analytical predictions do fit quite well. In the case of s , results are more scattered and realizations with one or more spanning clusters exhibit significantly lower values for s , in particular for higher viscosity contrasts. At a viscosity contrast of Δ = 10 −1 , the aDEM and the SCA approximations yield similar results, with the aDEM approximation performing slightly better than the SCA approximations. At larger viscosity contrasts, the viscosity predictions of the different approximations  (14) and (15). Solid lines denote the aDEM approximation by Dabrowski et al. (2012) and the SCA approximations by Fletcher (2004) and Treagus (2003). diverge. The best viscosity prediction is made by the SCA approximation of Fletcher (2004). The aDEM approach tends to overestimate the shear viscosity for larger viscosity contrasts, whereas the SCA approximation by Treagus (2003) underestimates it. Interestingly, these two approximations seem to form an envelope of the numerically determined effective viscosities, in particular at larger viscosity contrasts.
In summary, both numerical results and analytical approximations show that the existence of a weak phase does have an effect on effective composite rheology, even though the weak phase does not form an interconnected network. As expected, this effect becomes more pronounced for larger elongations. The effect of random weak phase morphology can be quantified with a statistical description employing the 2-point cluster function C 2 and analytical approximations.   (14) and (15). Solid colored lines denote the aDEM approximation by Dabrowski et al. (2012) and the SCA approximations by Fletcher (2004) and Treagus (2003). Gray areas indicate where numerical results become less reliable due to loss of resolution.

Microstructure Evolution During Simple Shear Deformation
The results of the previous section have shown that it is possible to quantify the effect of weak phase morphology on the effective rheology of a random heterogeneous medium by using a combination of its statistical description, namely, the 2-point cluster function together with an analytical approximation. However, the simulations shown were only "instantaneous" simulations and did not address the formation of elongated structures during deformation. To further explore the potential of the combined statistical-analytical approach to predict not only effective viscosities for a given microstructure but also its evolution during deformation, we conducted an additional set of experiments for the same volume fractions and viscosity contrasts shown in Figure 6. For each parameter combination, we deformed 10 different random samples up to a strain of 5 (which amounts to 10 4 time steps per sample). The initial phase morphology was chosen to be statistically isotropic.
The temporal evolution of a sample with a volume fraction of 20% and a viscosity contrast Δ = 10 −2 is shown in the top row of Figure 7. As expected the weak phase becomes increasingly elongated during deformation and eventually forms thin layers along which deformation is localized. This increasing localization can be seen in the strain rate distributions within both the hard and weak phase shown in the bottom row of Figure 7. Here, the normalized probability distribution of the second invariant of the strain rate tensor . II is shown for both the hard and the weak phase. The increased localization within the weak phase causes a shift of the maximum of the weak strain rate distribution to significantly larger values.
The morphology of the weak phase can be described using the 2-point cluster function C 2 (x). The shape of C 2 (x) reflects the increasing weak phase elongation and aligned layer formation during simple shear Figure 10. Comparison of numerically obtained inclination angle (left) and aspect ratio (right) to analytical predictions. Shown are the results from the simulation set with Δ = 10 −2 and Φ = 0.2 as an example (all other simulation sets show the same behavior). The dashed black line denotes the inclination angle obtained from principal component analysis of the viscosity tensor, whereas the solid black line denotes the inclination angle (left) and aspect ratio (right) obtained from analysis of the 2-point cluster function. Gray areas indicate strains at which results become unreliable due to loss of resolution. deformation (Figure 8). Employing the same procedure as in the single step simulations, we extract the effective normal and shear viscosity of the deformed material using equation (21) and (22). As expected, the normal and shear viscosities n and s exhibit an increase respectively a decrease corresponding to weak layer formation. We observe this behavior for all samples.
In Figure 9, we show the average evolution of the effective normal and shear viscosities for phase fractions Φ of 0.05, 0.1, and 0.2 and viscosity contrasts Δ of 10 −1 , 10 −2 , and 10 −3 . The evolution of normal viscosities (upper lines) is predicted fairly well by all approximations at smaller viscosity contrasts but shows increasing discrepancies at larger viscosity contrasts. Here, the analytical approximations tend to overestimate normal viscosities. All analytical approximations predict shear viscosities decreasing faster than what is observed in the numerical models and also predict considerably smaller viscosities at large strains. At large strains, some numerical results even indicate an increase in shear viscosity. However, this increase is related to a loss of resolution where the thickness of elongated structures becomes less than the grid resolution. For this reason, numerical results at strains above 3 become inaccurate.
The discrepancy between analytical predictions and numerical results can be explained with the differences between the assumptions made for the analytical approximations and the actual numerical model setup. All analytical approximations assume that the weak phase is distributed as elliptical inclusions that do not interact with each other. The assumption of elliptical inclusions also implies that the deformation gradient within each inclusion is constant (Eshelby, 1957), which is not the case for randomly heterogeneous media. The heterogeneity of the deformation field within the weak phase can also be seen in Figure 7, where we observe a large range of strain rates in the weak phase.
These discrepancies in strain rates do not seem to have such a large effect on the prediction of the inclination angle. Both the inclination angle obtained from principal component analysis of the viscosity tensor and the angle extracted from the 2-point cluster correlation function C2 provide very similar results and are predicted very well by the analytical approximations, as can be seen in Figure 10. However, the aspect ratio is predicted to increase much faster than what we observe in the numerical models, which also saturate at larger strains. The saturation behavior of the numerically determined aspect ratio can be explained with the finite size effect of the model domain that limits the horizontal extent of the 2-point cluster correlation function, while the lesser increase in aspect ratio may be related to the fact that the analytical predictions assume a constant strain rate field within an ellipse to determine the change in inclination angle. As already stated above, this assumption does not represent the actually observed strain rate field within the weak phase ( Figure 7) and is likely the cause for the discrepancy between analytical predictions and numerical results.

Discussion
In geodynamics, the rheology of the lower mantle is commonly assumed to be linear viscous and isotropic. Based on a large set of numerical simulations, we have shown here that the elongation of ferropericlase can have a significant impact on the effective viscosity and produces viscous anisotropy.
Notably, the viscosity of the lower mantle will already be significantly lowered upon deformation even before the weak ferropericlase develops an interconnected weak layer network. This is contrary to previous works that suggested a binary viscosity scenario, where either bridgmanite dominates lower mantle viscosity through the formation of a load-bearing framework or ferropericlase dominates viscosity once an interconnected weak layer network forms (Ballmer et al., 2017;Marquardt & Miyagi, 2015). Our modeling shows that there is a gradual transition between these two endmember scenarios. This implies that strain localization might become self-accelerating in lower mantle regions where initial deformation occurred, such as around subducting slabs that enter the lower mantle.
However, the resulting viscous anisotropy can reach several orders of magnitude between normal and shear direction. Importantly, the direction of lowest viscosity is not exactly aligned with the shear direction. Upon initial deformation, the angle between viscosity minimum and shear direction is about 45 • , but continuously rotates toward the shear direction during further deformation. The development of viscous anisotropy in the vicinity of subducting slabs where deformation will be localized is likely to affect the slab's descent into the lower mantle. One can envisage a scenario where the slab rotates in the lower mantle to move along the direction of minimal viscosity.
The presence of viscous anisotropy in the lower mantle is of relevance for modeling a wide range of geophysical measurements, including the overall viscosity profile of the lower mantle or postglacial rebound (Christensen, 1987). The here-predicted viscous anisotropy in the lower mantle is significantly larger than the viscous anisotropy predicted for olivine in the upper mantle that has been invoked to explain several large-scale geophysical phenomena (Hansen et al., 2012).
The development of a preferred orientation of ferropericlase in the bridgmanite matrix might also contribute to seismic anisotropy observations in the uppermost lower mantle  and possibly the D" layer (Kendall & Silver, 1996) given that the velocity contrast between bridgmanite and ferropericlase is sufficiently strong. Qualitatively, we would expect a V SH > V SV seismic shear wave splitting below regions where slabs move horizontally in the uppermost mantle, consistent with recent observations , which would require significantly differing seismic properties of ferropericlase and bridgmanite. The seismic anisotropy originating from crystallographic preferred orientation (CPO) alignment, however, may be reduced due to the heterogeneous deformation field in ferropericlase (e.g., Miyagi & Wenk, 2016). This may (at least partly) explain the absence of significant seismic anisotropy in most of the lower mantle. In the lowermost mantle, however, the transition of bridgmanite to postperovskite potentially results in a significant viscosity decrease of the matrix surrounding the ferropericlase inclusions (e.g., Ammann et al., 2010), which may result in a less heterogeneous deformation field and thus allow for consistent CPO evolution resulting in a stronger seismic anisotropy. However, this strongly depends on the active deformation mechanisms at this depth, which are yet unclear (e.g., Reali et al., 2019). Further studies are required to better quantify the effect of rheology on ferropericlase shape preferred orientation formation, CPO evolution, and the related seismic anisotropy.
Our results show that the prediction of the mechanical properties of a random heterogeneous bridgmanite-ferropericlase mixture can be predicted via a combination of statistical descriptors, namely, the 2-point cluster correlation function and analytical approximations for simpler systems of aligned ellipses. While the analytical predictions for the effective viscosity at given aspect ratios do successfully predict the numerical results, predictions for the evolution of aspect ratio and inclination angle predict a faster increase in aspect ratio than what is observed in numerical models. The most likely explanation for this discrepancy lies within the overestimation of elongation. This issue will have to be addressed in future studies.
The simulations performed in this paper as well as the analytical predictions only consider two-dimensional structures and deformations. However, both deformation and structures in nature represent 3-D problems. Faccenda et al. (2019) performed 3-D models for the evolution of initially spherical inclusions during simple shear for viscosity contrasts of 0.1, 1, and 10. For weak inclusions, they observe flattening of the inclusions and out of plane spreading. The 2-D models assume the weak ferropericlase inclusions to be infinitely stretched in this direction, which is why they can be seen as an extreme endmember of the 3-D case. The results presented in this study therefore have to be seen as conservative estimates for specific geometries. In his review on the subject, Jiang (2016) presents a comprehensive overview of the theoretical developments of the deformation of viscous inclusions and also provides general expressions of inclusion deformation for the 3-D case. However, "The behaviors of viscous inclusions in general anisotropic materials are terra incognita" (Jiang, 2016), and research in this field is still ongoing. Extension of the numerical-statistical approach presented in this study to 3-D is in principle straightforward. At present, it is still hampered by the significant computational expense that comes with modeling the deformation of random heterogeneous media in 3-D with a satisfactory resolution, which should be remedied in the future.
Lastly, the rheology of the Earth's lower mantle is determined not only by the interaction between bridgmanite and ferropericlase but by a large number of other parameters. Christensen (1987) already stated that "unfortunately, at the moment it is not clear if any of the predicted effects can be linked to observation in a definite manner, in order to decide whether or not viscous anisotropy plays a role in the mantle. In most cases there will be a trade-off with the influence of temperature, pressure, or stress on the rheology." This statement is still true as of today and all these effects have to be quantified in order to render a judgment on the importance of anisotropy development due to ferropericlase elongation. Figure A1. Effect of numerical resolution (n el denotes the number of elements used in each direction) on model results for (left) normal viscosity and (right) shear viscosity for different amounts of elongation (see legend). The weak phase is 100 times weaker than the strong phase.

Conclusions
In this study, we present a numerical-statistical approach to determine the effective viscosity of random heterogeneous two-phase media to estimate the impact of ferropericlase elongation on lower mantle rheology. We used the 2-point cluster correlation function to statistically describe microstructures that have been observed in experiments. These statistical properties were then used to artificially generate samples with equivalent properties. For a range of weak phase (ferropericlase) fractions and viscosity contrasts between the strong (bridgmanite) and weak phase (ferropericlase), the effective viscosity tensor was then numerically computed. Based on the observation that the 2-point cluster correlation function of the observed media exhibits an elliptical shape, we used analytical approximations for two-phase media with aligned elliptical inclusions to predict the effective viscosity tensor depending on inclusion elongation, phase fraction, and viscosity contrast. Numerical and analytical values were shown to agree very well, in most cases.
To determine the capability of the numerical approximations to predict not only the viscosity tensor of a given microstructure but also its evolution during simple shear deformation, we also performed time-dependent numerical simulations. For each simulation, we extracted the evolution of the viscosity tensor as well as the evolution of weak phase morphology as described by the 2-point cluster correlation function. Using the same connection between 2-point cluster correlation function, its elliptical shape, and the effective viscosity tensor, we used analytical solutions for the evolution of a two-phase composite with aligned elliptical inclusions to predict the weak phase morphology during simple shear deformation. These predictions do fit the numerical results to first order but tend to overestimate the stretching of the weak phase and thus predict stronger weakening over time compared to the numerical results. Further research is thus needed to better understand the differences between general random heterogeneous media and systems consisting of aligned elliptical inclusions. The extension of the 2-D results to 3-D likewise needs further research. Nevertheless, the present study shows that existing analytical predictions for the effective viscosity of two-phase composites can be indeed used to first order to predict the viscosity evolution of the bridgmanite-ferropericlase mixture in the lower mantle due to ferropericlase elongation.

Appendix A: Resolution Study
To determine the reliability of our model results with respect to model resolution, we performed a series of resolution tests for varying amounts of weak phase elongation. As can be seen in Figure A1, results for the effective normal viscosity n converge at moderate resolutions, but a much higher resolution is needed to resolve the effective shear viscosity s . For strongly elongated structures, results do not converge even at the highest resolution of 1,024 elements but show starting asymptotic behavior, as it is not possible anymore to resolve the vertical dimensions of some of the weak phase clusters. It is thus important to keep in mind that numerical results tend to be less accurate in case of strongly elongated clusters and likely overestimate Figure B1. Same as Figure 4 but for = 0.05. Simulations were performed on cluster btrzx2, University of Bayreuth. M. T. thanks M. Dabrowski for discussions and for sharing the software to compute the aDEM approximation. Figures shown in this manuscript are available from Thielmann (2019b). The numerical code to generate the data used in this manuscript is available from Thielmann (2019a). We want to thank two anonymous reviewers whose comments helped to improve and clarify the manuscript.