Testing the D-Rex model of flow-induced preferred orientation

Last modified on March 10, 2023

The D-Rex model is a kinematic model for flow-induced CPO (crystalline preferred orientation) in minerals. It has been used to simulate CPO development in upper mantle olivine by e.g. Kaminski & Ribe (2002), Hedjazian et al. (2017) and Fraters & Billen (2021). In this article I document some preliminary testing of the model. In particular, I derive some simple analytical results for comparison with simulation outputs. I am using an unfinished Python rewrite based on the reference implementation, which I plan to develop further as a standalone plugin for geodynamics simulations.

This article assumes a fair amount of background in crystallography and geodynamics. I will link to related introductory articles when I have written them. I am also attempting to use an unambiguous mathematical notation for this website.

Contents:

---

Overview of the D-Rex model

D-Rex is a kinematic model for flow-induced texture development in the Earth's mantle. It was first developed by Kaminski & Ribe (2001) for simulations of pure olivine minerals, but later expanded to handle two-phase olivine-enstatite aggregates. For an expanded overview, see "Kinematic models of flow-induced preferred orientation". Kinematic constraint CPO models are generally the simplest of CPO models. This also makes them the most computationally efficient.

The D-Rex model considers aggregate minerals, which are composed of a fixed number of "grains". However, these grains are not used to track grain boundaries or model realistic but computationally expensive microscopic stress variations. Rather, they are used as simulation nodes that store an orientation matrix and a volumetric weighting factor (called the grain "size").

Macroscopic deformation is accommodated by rigid-body rotation and dislocation creep. Dislocation creep requires a certain minimum grain size, as small grains are assumed to migrate via grain-boundary sliding. While the number of grains is fixed during the simulation, their sizes will change depending on the density of accumulated dislocations. The D-Rex model describes two processes that can cause grain size changes.

Grains with low dislocation density gradually annex parts of neighbouring grains if they have a significantly higher dislocation density. This is known as grain-boundary migration. Grain sizes also change based on nucleation of internal subgrains. Subgrain boundaries are caused by heterogeneity in the internal dislocation density, which lead to subgrain rotation. Both of these processes are particular cases of dynamic recrystallisation (Kaminski et al., 2004).

Because the relative positions of grains within the aggregate is not tracked in D-Rex, grain size changes are calculated by comparing a residual strain energy of each grain to the mean value for the aggregate. The strain energy formula includes the effects of both recrystallisation processes simultaneously. Also, because the number of grains is fixed, new subgrains simply inherit the rotation of their parent and only contribute by enforcing an upper bound on grain size.

Slip activation in olivine

Olivine crystallises with dipyramidal orthorhombic symmetry, which means that it has three symmetry planes. Slip may activate on any of these planes in any direction, however slip directions parallel to one of the symmetry axes tend to be more favourable. At high temperatures and stresses, activation of four slip systems has been experimentally demonstrated (Karato et al., 2008). These are depicted along with an illustration of the symmetry class by the octahedron below. Note that olivine crystals exhibit a more complex but symmetrically equivalent geometry.

Olivine may exhibit five different kinds of deformation response. The fabric of a particular olivine crystal depends on a variety of factors including hydration. Fabric types can be distinguished by their (dimensionless) critical resolved shear stress (CRSS)^[1] values which are experimentally determined (Kaminski et al., 2004, Karato et al., 2008):

These values define the threshold of shear stress, relative to the softest system, which is required to initiate slip on a particular system. A value of ∞ indicates a permanently inactive system. A-type is the most commonly observed of the listed fabrics (Karato et al., 2008).

Solution for a single grain under simple shear conditions

I will not reproduce the general theory from the papers here, and will focus directly on simplified analytical solutions to certain test cases.

In order to validate the outputs of the implementation, I will be adding a suite of tests to the code. First, I am checking that the rotation rates make sense for single grains. For these tests, the macroscopic deformation is set to a simple shear. The code does not yet handle transient CPO, so these tests are being written for steady-state deformation gradients. I am only considering olivine grains for now.

The macroscopic deformation is induced by a velocity gradient \(\bm{L}=∇\bm{v}\), which can be decomposed into symmetric and antisymmetric components,

\[ \dot{ε}_{0}\bm{L} = \dot{ε}_{0}(\dot{\bm{ε}} + \bm{Ω}) \]

where \(\dot{ε}_{0}\) is the strain rate scale, which is removed to achieve a convenient dimensionless representation. Let \(\begin{Bmatrix} \bm{\hat{x}}_{i}\end{Bmatrix} \in ℝ^3\) denote an orthonormal basis. Simple shear in the plane spanned by \(\begin{Bmatrix} \bm{\hat{x}}_{1}, \bm{\hat{x}}_{2}\end{Bmatrix}\) can be defined as

\[ \bm{L} = \begin{bmatrix} 0 & 0 & 2 \\ 0 & 0 & 0 \\ 0 & 0 & 0\end{bmatrix} ,\quad \bm{\dot{ε}} = \begin{bmatrix} 0 & 0 & 1 \\ 0 & 0 & 0 \\ 1 & 0 & 0\end{bmatrix} ,\quad \bm{Ω} = \begin{bmatrix} 0 & 0 & 1 \\ 0 & 0 & 0 \\ -1 & 0 & 0\end{bmatrix} \]

where the matrices \(\dot{\bm{ε}}\) and \(\bm{Ω}\) represent the strain rate and vorticity, respectively.

I will first derive a solution for rotation of the grain around a single axis, which is a simple and intuitive case. I choose an initial rotation around \({\bm{\hat{x}}}_{2}\) so that the grain will be subjected to the macroscopic shear. The orientation matrix in the global, i.e. macroscopic frame is^[2]

\[ \bm{a} = \begin{bmatrix} \cosθ & 0 & \sinθ \\ 0 & 1 & 0 \\ -\sinθ & 0 & \cosθ\end{bmatrix} \]

First, we extract the strain rate invariants for the four slip systems of olivine. I will call these the slip invariants. They are extracted from the strain rate tensor by calculating

\[ I = l^{i}n^{j} \dot{ε}_{ij} \]

for each slip system, where \(\bm{l}\) and \(\bm{n}\) are unit vectors parallel and normal to the slip plane, respectively. These are found by extracting the appropriate rows from the orientation matrix \(\bm{a}\) by multiplying the transpose of the orientation matrix by the basis vector, e.g. \(\bm{n} = \bm{a}^{⊺}\bm{\tilde{n}}\), as shown in the slip vector table:

For the strain rate from (2), only two terms will contribute to the calculation of (4), namely the \(\bm{l}_{1}\bm{n}_{3}\) and \(\bm{l}_{3}\bm{n}_{1}\) terms. Considering A-type olivine, and using superscripts \(s = \{a,b,c,d\}\) to represent the slip systems as ordered in the slip vector table, we have

\[ \begin{align*} I^{a} &= \bm{l}^{a}_{1} \bm{n}^{a}_{3} + \bm{l}^{a}_{3} \bm{n}^{a}_{1} = \bm{a}_{11} \bm{a}_{23} + \bm{a}_{13} \bm{a}_{21} = 0 \\ I^{b} &= \bm{l}^{b}_{1} \bm{n}^{b}_{3} + \bm{l}^{b}_{3} \bm{n}^{b}_{1} = \bm{a}_{11} \bm{a}_{33} + \bm{a}_{13} \bm{a}_{31} = \cos^{2}θ - \sin^{2}θ = \cos2θ \\ I^{c} &= \bm{l}^{c}_{1} \bm{n}^{c}_{3} + \bm{l}^{c}_{3} \bm{n}^{c}_{1} = \bm{a}_{31} \bm{a}_{23} + \bm{a}_{33} \bm{a}_{21} = 0 \\ I^{d} &= - I^{b} \end{align*} \]

The activity of each slip system is proportional to the ratio of the slip invariant \(I^{s}\) to the CRSS value \(τ^{s}\). In this case, only the second slip system is active, because the invariant of the fourth system is irrelevant (divided by infinity) for A-type olivine.

The microscopic velocity gradient, or deformation rate \(\bm{d}\), can be decomposed using an equation analogous to (1),

\[ \bm{d} = \bm{G}γ + \bm{ω} \]

where \(\bm{G}\) is a tensor describing the relative slip rates, which is multiplied by the slip rate on the most active slip system, \(γ\). The second quantity, \(\bm{ω}\), describes the rotation rate of the crystallographic axes (Kaminski & Ribe, 2001). Knowing the relative slip activities \(β^{s}\), we can construct the tensor

\[ G_{ij} = 2 \sum_{s} β^{s} l^{s}_{i} n^{s}_{j} \]

where in this case we don't need the sum over \(s\) and \(β=1\) because only one system is active, i.e. the \((001)[100]\) system. The only nonzero components are

\[ \begin{align*} G_{11} &= 2 \bm{l}_{1} \bm{n}_{1} = 2 \bm{a}_{11} \bm{a}_{31} = - 2 \sinθ\cosθ \\ G_{13} &= 2 \bm{l}_{1} \bm{n}_{3} = 2 \bm{a}_{11} \bm{a}_{33} = 2 \cos^{2}θ \\ G_{31} &= 2 \bm{l}_{3} \bm{n}_{1} = 2 \bm{a}_{13} \bm{a}_{31} = - 2 \sin^{2}θ \\ G_{33} &= 2 \bm{l}_{3} \bm{n}_{3} = 2 \bm{a}_{13} \bm{a}_{33} = 2 \cosθ\sinθ = -G_{11} \\ \end{align*} \]

Derivation of the general expression for the slip rate on the most active system, \(γ\), is omitted here for the sake of brevity. Note that the version in Kaminski & Ribe (2001) contains an error, which has been corrected in equation 4 of Fraters & Billen (2021). For simple shear along one of the global axes, the expression simplifies to

\[ γ = \frac{ 2 \bm{G}_{ij} \bm{L}_{ij} - \left(\bm{L}_{k,k+1} - \bm{L}_{k+1,k}\right) \left(\bm{G}_{k,k+1} - \bm{G}_{k+1,k}\right) }{ 2 (\bm{G}_{ii}^2 + \bm{G}_{ij}^2 + \bm{G}_{ji}^2 + \bm{G}_{jj}^2) - \left(\bm{G}_{k,k+1} - \bm{G}_{k+1,k}\right)^{2} } \]

for the pair of indices \((i,j)\) that denote the location of the only nonzero velocity gradient entry, and the indices \((k,k+1)\) that satisfy \(k + 1 = a + 2 \mod 3\) for the rotation axis \(\bm{\hat{x}}_a\) in the right-handed coordinate system. Evaluating the individual terms for \(\bm{\hat{x}}_{2}\), we obtain

\[ \begin{align*} 2 \bm{G}_{13} \bm{L}_{13} &= 2 ⋅ 2\cos^{2}θ ⋅ 2 = 8\cos^{2}θ, \\ (0 - \bm{L}_{13})(\bm{G}_{31} - \bm{G}_{13}) &= -2\left(-2\cos^{2}θ - 2\sin^{2}θ\right) = 4, \\ 2\left(\bm{G}_{11}^2 + \bm{G}_{13}^2 + \bm{G}_{31}^2 + \bm{G}_{33}^2\right) &= 2\left(2\bm{G}_{11}^2 + \bm{G}_{13}^2 + \bm{G}_{31}^2\right) \\ &= 2\left(2 ⋅ 4\sin^{2}θ\cos^{2}θ + 4\cos^{4}θ + 4\sin^{4}θ\right) \\ &= 2\left(2\cos^{2}θ + 2\sin^{2}θ\right)^{2} = 8, \\ \left(\bm{G}_{31} - \bm{G}_{13}\right)^{2} &= \left(-2\cos^{2}θ - 2\sin^{2}θ\right)^{2} = 4 \\ \end{align*} \]

which gives

\[ γ = \frac{8\cos^{2}θ - 4}{8 - 4} = 2\cos^{2}θ - 1 = \cos2θ \]

The rotation rate of the crystalline lattice, \(\bm{ω}\), can be obtained from (Kaminski & Ribe, 2001)

\[ \begin{align*} \bm{\bar{ω}}_{i} &= \sum_{i=1}^{3} \frac{\bm{L}_{i+2,i+1} - \bm{L}_{i+1,i+2}}{2} - \frac{γ\left(\bm{G}_{i+2,i+1} - \bm{G}_{i+1,i+2}\right)}{2} \\ ω^{ij} &= ε^{i\phantom{k}j}_{\phantom{i}k\phantom{j}}\bar{ω}^{k} \\ \end{align*} \]

where the indices of \(\bm{L}\) and \(\bm{G}\) are cyclical from 1 to 3, and \(\bm{ε}\) is the permutation symbol,

\[ ε_{ijk} = \begin{cases} +1 &\text{if } (i,j,k) \text{ is } (1,2,3), (2,3,1), \text{ or } (3,1,2), \\ -1 &\text{if } (i,j,k) \text{ is } (3,2,1), (1,3,2), \text{ or } (2,1,3), \\ 0 &\text{if } i = j, \text{ or } j = k, \text{ or } k = i \\ \end{cases} \]

We only need to evaluate one entry of the \(\bm{\bar{ω}}\) spin vector, because there is only one nonzero element of the velocity gradient matrix:

\[ \begin{align*} \bm{\bar{ω}}_{2} &= \frac{\bm{L}_{13} - \bm{L}_{31}}{2} - \frac{γ\left(\bm{G}_{13} - \bm{G}_{31}\right)}{2} \\ &= \frac{2 - 0}{2} - \frac{\cos2θ\left(2\cos^{2}θ + 2\sin^{2}θ\right)}{2} \\ &= 1 - \cos2θ \end{align*} \]

which yields the rotation rate

\[ \bm{ω} = \begin{bmatrix} 0 & 0 & 1 - \cos2θ \\ 0 & 0 & 0 \\ \cos2θ - 1 & 0 & 0\end{bmatrix} \]

The change in crystalline axes directions is given by

\[ d_{t^{\ast}}\bm{a} = \bm{ω}^{⊺} \bm{a} \]

for the dimensionless timescale \(t^{\ast} = t \dot{ε}_{0}\), with \(\dot{ε}_{0}\) defined in (1). The matrix transpose is needed once again due to using the active rotation convention. Rotation of coordinate axes is the inverse operation of vector rotation. The matrix entries evaluate to

\[ \begin{align*} d_{t^{\ast}}\bm{a} &= \begin{bmatrix} 0 & 0 & \cos2θ - 1 \\ 0 & 0 & 0 \\ 1 - \cos2θ & 0 & 0 \\\end{bmatrix} \begin{bmatrix} \cosθ & 0 & \sinθ \\ 0 & 1 & 0 \\ -\sinθ & 0 & \cosθ \\\end{bmatrix} \\ &= \begin{bmatrix} \sinθ (1 - \cos2θ) & 0 & \cosθ (\cos2θ - 1) \\ 0 & 0 & 0 \\ \cosθ (1 - \cos2θ) & 0 & \sinθ (1 - \cos2θ) \\\end{bmatrix} \end{align*} \]

To verify that the analytical result makes sense, let us consider some characteristics of the solution.

Firstly, \(θ = ±45° ⟹ \cos2θ = 0\) and the rotation rate from (15) will be equal to the macroscopic vorticity from (2). This means that no dislocation creep is occurring, and the grain will simply be advected within the viscous flow. None of the slip systems are active (all of the slip invariants from (5) are zero). This makes sense, because for these orientations the only slip system that can activate for this rotational axis (the \((001)[100]\) system) is aligned with either purely compressional or purely extensional macroscopic forces, i.e. one of the principal axes of the finite strain ellipse.

For \(θ = 0° ⟹ \cos2θ = 1\), the rotation rate from (15) will be a null matrix, which means that the dislocation creep is completely accommodating the imposed vorticity. The grain is not rotating at all, because the \((001)[100]\) slip system is at the most favourable orientation, i.e. aligned to the shear plane. Dislocation planes will allow for the individual grain to undergo shear. For an analogy, consider how a deck of cards would respond to simple shear imposed at the top.

When \(θ = ±90° ⟹ \cos2θ = -1\), the norm of the rotation rate from (15) will be at its maximum. The slip plane is furthest from equilibrium and therefore rotates fastest, and realignment due to dislocation creep acts in phase with the rigid-body rotation.

In conclusion, the overall expected physical behaviour is reproduced by the model. Derivations for other configurations of simple shear and single-axis rotation are very similar. For example, a simple shear defined by

\[ \bm{L} = \begin{bmatrix} 0 & 0 & 0 \\ 2 & 0 & 0 \\ 0 & 0 & 0\end{bmatrix} ,\quad \bm{\dot{ε}} = \begin{bmatrix} 0 & 1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 0\end{bmatrix} ,\quad \bm{Ω} = \begin{bmatrix} 0 & -1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 0\end{bmatrix} \]

and initial grain orientation of

\[ \bm{a} = \begin{bmatrix} \cosθ & -\sinθ & 0 \\ \sinθ & \cosθ & 0 \\ 0 & 0 & 1 \\\end{bmatrix} \]

result in the rotation rate

\[ \bm{ω} = \begin{bmatrix} 0 & -1 - \cos2θ & 0 \\ 1 + \cos2θ & 0 & 0 \\ 0 & 0 & 0 \\\end{bmatrix} \]

and a change in orientation given by

\[ d_{t^{\ast}}\bm{a} = \begin{bmatrix} \sinθ (1 + \cos2θ) & \cosθ (1 + \cos2θ) & 0 \\ -\cosθ (1 + \cos2θ) & \sinθ (1 + \cos2θ) & 0 \\ 0 & 0 & 0 \\\end{bmatrix} \]

Using these results, I have been able to implement analytical tests for the code. Statistical tests for mineral aggregates will follow. This article was necessarily technical, and is primarily intended as a supplementary reference for my future work. Documentation of the work-in-progress software may at some point refer to or incorporate this information.

---

[1]	In some of the referenced literature this is called reference resolved shear stress (RRSS).

[2]

I have found that the reference implementation and existing literature does not clearly qualify its rotation conventions. For this article and the code I am developing I will be adopting active rotations relative to the right-handed global reference frame, using pre-multiplication (\(\bm{v}'=R\bm{v}\)). This matches the conventions used by the scipy.spatial.transform.Rotation class. For the Euler decomposition, I will prefer the extrinsic 3-1-3 (ZXZ) convention, although the as_euler method allows for any choice of decomposition.

References

• Fraters, M. R. T., & Billen, M. I. (2021). On the Implementation and Usability of Crystal Preferred Orientation Evolution in Geodynamic Modeling. Geochemistry, Geophysics, Geosystems, 22(10). https://doi.org/10.1029/2021gc009846
• Hedjazian, N., Garel, F., Davies, D. R., & Kaminski, É. (2017). Age-independent seismic anisotropy under oceanic plates explained by strain history in the asthenosphere. Earth and Planetary Science Letters, 460, 135–142. https://doi.org/10.1016/j.epsl.2016.12.004
• Kaminski, É., & Ribe, N. M. (2001). A kinematic model for recrystallization and texture development in olivine polycrystals. Earth and Planetary Science Letters, 189(3–4), 253–267. https://doi.org/10.1016/s0012-821x(01)00356-9
• Kaminski, É., & Ribe, N. M. (2002). Timescales for the evolution of seismic anisotropy in mantle flow. Geochemistry, Geophysics, Geosystems, 3(8), 1–17. https://doi.org/10.1029/2001gc000222
• Kaminski, É., Ribe, N. M., & Browaeys, J. T. (2004). D-Rex, a program for calculation of seismic anisotropy due to crystal lattice preferred orientation in the convective upper mantle. Geophysical Journal International, 158(2), 744–752. https://doi.org/10.1111/j.1365-246x.2004.02308.x
• Karato, S., Jung, H., Katayama, I., & Skemer, P. (2008). Geodynamic Significance of Seismic Anisotropy of the Upper Mantle: New Insights from Laboratory Studies. Annual Review of Earth and Planetary Sciences, 36(1), 59–95. https://doi.org/10.1146/annurev.earth.36.031207.124120

Bibliography

• Karato, S. (2008). Deformation of Earth Materials. https://doi.org/10.1017/cbo9780511804892