On the Implementation and Usability of Crystal Preferred Orientation Evolution in Geodynamic Modeling

Flow in the Earth's mantle causes the preferred orientation of crystals called lattice/crystal preferred orientation (LPO or CPO). This preferred orientation is one of the main reasons why seismic anisotropy is observed. Seismic anisotropy observations could therefore be used to constrain the mantle flow in geodynamic models through tracking CPO evolution, and computing the resulting elastic tensor and the anisotropy predicted at the surface. Even though there are many types of CPO models, only a few studies include CPO calculations due to the complexity and computational cost. Here, we implemented an extended version of the CPO model D‐Rex into the open‐source community geodynamics code ASPECT. We show that the implementation is correct, how to use it, and that it is feasible and important to use in large 3D models. We also show that it is important to calculate CPO, especially for models focusing on plate boundary of smaller scale flow because the resulting fast axis directions can greatly deviate from the flow direction. The added infrastructure will also allow for future enhancement, testing, and even replacement of the CPO model.

2 of 27 randomly distributed, no seismic anisotropy will be observed, meaning that the mantle will appear isotropic (Karato et al., 2008). Conversely, when mineral grains are deformed by mantle flow the crystal lattice aligns in a preferred orientation determined by the shear strain of the flow. This alignment is referred to as lattice preferred orientation or crystallographic preferred orientation (CPO). When the shear strain causes strong alignment over a large enough volume, the CPO creates an observable anisotropy in seismic observations made at stations on the Earth's surface. Relating observations of seismic anisotropy to the fast axes direction of the CPO of olivine is difficult in practice. This is due to several possible origins of the seismic anisotropy besides CPO (including alignment of melt-filled lenses and channels [Holtzman et al., 2003], and alignment of fractures filled with fluids or distinct minerals [Faccenda et al., 2008;Ranero et al., 2003]) and because relating measurements of the fast axis direction of the CPO to mantle flow is not a simple one-to-one relationship either. The relationship between the mantle flow direction and the CPO depends on several factors. First, the CPO alignment is caused by shear strain, which is usually not aligned directly with the flow, except in cases like simple shear. Second, the development of CPO is not an instantaneous process, even on geologically relevant timescales, therefore the history of deformation will affect the present-day alignment (Wenk et al., 1991;Zhang & Karato, 1995). Furthermore, the direction of the CPO fast-axis can change up to 90°with respect to the shear strain direction due to changes in water content and stress applied to the minerals (Karato et al., 2008;Kneller & van Keken, 2008). Understanding mantle flow through seismic anisotropy and CPO will therefore require good constraints on both conditions in the mantle, such as the current and past mantle flow, water content and stress in the area of study, and a good model of how crystal lattices align at these conditions. Geodynamic models, used in conjunction with a range of geophysical observations, are a powerful tool to constrain unknown physical parameters (e.g., using the gravity field to constrain the mantle viscosity structure). However, despite the unique ability of seismic anisotropy to shed light on the details of flow in the mantle, the number of studies using this approach has been limited (e.g., Blackman et al., 1996;Conrad et al., 2007;Di Leo et al., 2014;Faccenda & Capitanio, 2013;Forte et al., 2010;Ito et al., 2014;Jadamec & Billen, 2010;Long & Becker, 2010;McNamara et al., 2001;Miller & Becker, 2012;Tommasi et al., 2004Tommasi et al., , 2009. The main reasons that there are a limited number of such studies are the high computational cost of most CPO models and the lack of integration between CPO modeling software and geodynamic modeling software. To overcome this barrier to scientific advancement, we have integrated a CPO model with geodynamic modeling software. There are many numerical models in the geodynamic modeling community, each built for a different purpose. In this study we will be implementing a CPO model using the D-Rex method (Kaminski & Ribe, 2001;Kaminski et al., 2004) into a geodynamic modeling software package called Aspect (Bangerth, Dannberg, et al., 2020;Heister et al., 2017;Kronbichler et al., 2012). The main advantage of using Aspect is that it is an actively developed open source community code built to solve geodynamic problems that are highly scalable (from a laptop up to tens of thousands of cores), and it is widely used (e.g., CIG tracking of Aspect publications lists about 100 publications between 2012 and 2021, CIG, 2021). Furthermore, it contains many advanced methods that are needed or useful for geodynamic CPO modeling such as a scalable and flexible Particle-in-Cell method that works with adaptive mesh refinement (Gassmöller et al., 2018), advanced nonlinear rheologies (Glerum et al., 2018), and a stabilized Newton solver to accurately solve for the flow field with nonlinear rheology . Similarly, there are many different CPO modeling algorithms available, with each having its advantages and disadvantages. The CPO models can be divided into two broad groups: models based on the finite strain ellipsoid (FSE) and polycrystal models. The FSE based models use the strain at a location to compute the CPO of the material and assume that the CPO is independent of the deformation path. Using the FSE to compute the CPO is based on the hypothesis that the CPO is a unique function of the finite strain (McKenzie, 1979). This hypothesis was tested against the more accurate, but computationally expensive, visco-plastic self consistent (VPSC) type of CPO models in simulations of varying complexity and found acceptably similar results in most flow fields (Blackman & Kendall, 2002;Li et al., 2014;Ribe, 1992). FSE models have been used by many authors to interpret seismic anisotropy measurements (e.g., Becker et al., 2003;Long et al., 2007;Morgan, 1987). However, because FSE models do not track grains, they can not account for an initial CPO or CPO evolution. This makes this method much less accurate in settings where strain rates can change significantly both spatially and temporally, such as in subduction zones. A special case of FSE models is the Infinite Strain Axis (ISA) model. The ISA is defined as the long axis of finite strain ellipsoid in the limit of an infinite strain (Kaminski & Ribe, 2002). Jadamec and Billen (2010) showed that the results from this 3 of 27 method were generally consistent with the SKS splitting measurements for the eastern Alaska subduction zone. This agreement was likely the result of a high strain rate immediately surrounding the slab, such that the instantaneous flow and ISA approximation captured the full re-alignment of the CPO. This approach has also been used computing the seismic anisotropy for global mantle flow (e.g., Conrad et al., 2007). Polycrystal type CPO methods are models that track the evolution of individual grains to arrive at the CPO of the material. In the polycrystal group, there are three classes: full-field models, mean-field homogenization models, and kinematic models. Full field models (e.g., Lebensohn, 2001;Moulinec & Suquet, 1998;Sarma & Dawson, 1996) treat both the spatial distribution of continuous grains and the stresses and allow stress and strain to vary both among and within individual grains in a physically realistic way. This realism comes at the expense of being computationally many orders of magnitude too slow for routine use in geodynamic models (Ribe et al., 2019). Mean-field homogenization models are computationally faster than full-field models because they do not use spatial information of the individual grains. Compatibility of stress and strain is not enforced between spatially contiguous grains, but between each phase and a homogeneous effective medium defined by the average of all the other phases. However, this method is four orders of magnitude slower than the kinematic type of models (Ribe et al., 2019) although faster versions are being developed (Signorelli et al., 2020). Examples of mean-field models are the VPSC model (Lebensohn & Tomé, 1993;Molinari et al., 1987) and the second-order self-consistent model (Liu & Ponte Castañeda, 2004;Ponte Castañeda, 2002).
The kinematic class of models is based on a quantitative description of the deformation-induced rate of crystallographic rotation as a function of grain orientation. Kinematic models do not account explicitly for stress compatibility among grains and therefore can not be used to predict rheological properties of a deforming aggregate (Ribe et al., 2019). However, the effects of secondary phases, such as enstatite (Kaminski et al., 2004), water (Kaminski, 2002), and melt (Kaminski, 2006), can be included in the method. In this paper, we show how we implemented a modified version of the polycrystal kinematic CPO model D-Rex (Kaminski & Ribe, 2001;Kaminski et al., 2004) into the geodynamics code Aspect. The choice for a polycrystal type model comes from our interest in having initial and history-based evolution of CPO in our models, which is not possible in FSE models. The choice for a kinematic type model is based on its low computational cost while still providing accurate fast axis directions under most conditions. A final note of the advantages of using D-Rex is that it has the ability to incorporate enstatite and dynamic recrystallization into the calculation. We show that the implementation is fast enough for routine usage in geodynamic models, how to use the new implementation, and some of the limitations with the current implementation.

Geodynamics Simulation Software: Aspect
Aspect is an open source finite element model geodynamics code (Bangerth, Dannberg, et al., 2020;;Heister et al., 2017;Kronbichler et al., 2012), which is built on Deal.ii 9.2.0 (Bangerth et al., 2011;Bangerth, Maier, et al., 2020), p4est (Burstedde et al., 2011), and trilinos (Trilinos Project Team, n.d.). It contains many advanced features such as two-phase flow for melts (Dannberg & Heister, 2016), a free surface (Rose et al., 2017), nonlinear viscoplasticity (Glerum et al., 2018), a Newton solver for nonlinear Stokes problems , a scalable particle-in-cell method for massively parallel computation (Gassmöller et al., 2018) and Matrix free solvers (Clevenger & Heister, 2019). Aspect also has an integration with the Geodynamic World Builder (Fraters, 2020; which allows for the easy and reproducible setup of models, and which we will use in the subduction example shown here. The D-Rex kinematic CPO algorithm operates on an aggregate of grains. In this study, we have chosen to use the particles in Aspect (Gassmöller et al., 2018) to store and advect these properties. The main advantages of using particles (and not a field variable) are that they can be distributed independently of the mesh, so that only areas of interest have the particles, and that the implementation can efficiently handle a large amount of data the particles need to carry. To fully integrate the D-Rex algorithm in Aspect, it is necessary to: 1. define input variables that will allow the user to specify the initial CPO (random or with a preferred alignment), and volume fraction distribution (random, uniform, or a function of the orientation), parameters defining the mineral composition (volume fraction of olivine and enstatite), and parameters that control the recrystallization behavior 4 of 27 2. create the initial conditions for the aggregates of grains that will evolve in time as the particles are advected with the flow 3. compute the evolution of the CPO and grain size using the flow solution and other environmental variables (e.g., water content) at each time step, and 4. post process the CPO solution for output and analysis Our implementation is set up in the following ways: 1. each particle can contain several minerals (e.g., olivine, enstatite) 2. each mineral contains a number of grains 3. each grain contains a (a) nondimensional size where all grain sizes of a mineral sum to one, and (b) rotation matrix, which tracks the CPO.
Because it is useful to create initial conditions for the CPO rotations based on Euler angles and also output Euler angles from general rotation matrices, we can use a direction cosine matrix (DCM) to construct and derive Euler angles. DCM is a general rotation matrix formed by the multiplication of three specific (or basic) successive rotation matrices, where each matrix rotates about a separate axis. The angles specifying each of these rotations are the Euler angles. When Euler angles are used, the implementation uses the rotations around the Z-X-Z axes respectively, which is the same as in D-Rex (Ribe & Yu, 1991).

CPO Evolution Algorithm: D-Rex
To implement the D-Rex algorithm in Aspect, we followed the algorithm descriptions provided in several publications (Kaminski & Ribe, 2001;Kaminski et al., 2004;Ribe & Yu, 1991). In some instances, there are discrepancies or incomplete information in the publications. In these cases, we consulted the Fortran90 D-Rex code (Kaminski, 2003) and a separate Matlab version (Thissen, 2015). Therefore, in this section, we provide a complete description of the algorithm as we have implemented it in Aspect. We also include some extra information on differences between the provided equations and the actual implementation in the D-Rex code and this code and explain why they are equivalent or justified. All variables and definitions related to the D-Rex algorithm can be found in Appendix A. The definitions in Table A2 are ordered such that they can be read as an algorithm. The particles containing the grains and their orientations are advected with the velocity field. This means that even if the orientations of the grains would remain the same with respect to the particle, they would rotate with the particle (see Appendix B on the passive rotation benchmark). To ensure that the grain properties are properly advected with the flow, and advection scheme is required. In order of rising stability, accuracy, and computational cost, we have implemented the For- The D-Rex kinematic CPO model captures four of the main process that determine the evolution of CPO in an aggregate. For more information on the processes see Tullis (2002). Two processes determine the rotation and two the grain size. The rotation is determined by: (a) rotation of the crystal lattice due to dislocation creep and (b) grain boundary sliding for small grains. First, the orientation of the crystal lattice changes due to the motion of dislocations on slip planes in the crystal lattice. The rotation of the crystal lattice depends on the orientation of the slip planes relative to the deformation field (the strain-rate tensor), and the amount of slip on each slip system as determined by the relative strength of each slip system. Second, very small grains will not deform by dislocation creep, and therefore their lattice orientations will not evolve in time.
The change in grain size is determined by (a) change in grain size due to grain boundary migration and 5 of 27 (b) subgrain rotation, which forms small strain-free grains. First, changes in the grain size due to grain boundary migration are calculated. Grain boundary migration occurs due to differences in the internal strain energy of adjacent grains caused by differences in the number of dislocations in the grain. The number of dislocations tends to increase with continued deformation. Grains with lower strain energy tend to grow into adjacent higher strain energy grains, reducing the dislocation density (and the strain energy) in the larger grains. Second, new strain-free grains form from grains with high strain energy (i.e., high dislocation density). Note that these new small strain-free grains inherit the rotation of the old grain because there is a fixed number of grains and the accounting of small grains is the same for both grain boundary sliding and subgrain rotation.

Computation of Olivine Fabric
Before the CPO update can be computed, first the relative activities of the different slip systems must be determined. This is controlled by the fabric of the mineral. For olivine, the fabric type depends on both the water content and the state of stress, while for enstatite there is in practice only one fabric. The fabric is determined by the activity of a slip system, which is primarily controlled by the length of the Burgers vector, favoring the shortest vector (Karato et al., 2008). For olivine, both Burgers vectors  have only a small difference in length and therefore changes in the physical and/or chemical environment influence the dominant slip direction, leading to different fabrics. The fabric for olivine is highly dependent on stress and water content. CPO in orthopyroxenes such as enstatite has not been studied as well as for olivine. Although orthopyroxene CPO has been observed to be able to have slip-on both the [100] and [001] planes (Carter et al., 2013), the slip system (100) [001] E is seen as the dominant slip system (Karato et al., 2008). This can be explained by the length of the  [001] E b Burgers vector, which is much shorter than the length of the other Burgers vectors (Karato et al., 2008). The fabric properties are summarized in Table 1. To determine the deformation type for olivine in our implementation in Aspect, olivine fabric is determined from the current water content and magnitude of stress, based on information from the mesh, at the location of the particle and are by default compared against a digitized version of Figure 4 in Karato et al. (2008) to determine the correct fabric (see Appendix C). The boundaries between the different fabrics are not well constrained and more data is needed (Ohuchi et al., 2012).

Change in Orientation of the Crystals
The change in orientation of a grain is given by Equations 1 and 2 (adapted from Equation 9 in É. Kaminski & Ribe, 2001). The difference between Equation 9 in Kaminski and Ribe (2001) and Equations 1 and 2 is that ij E a , which is the current orientation of the crystallographic axis, has been separated from the rest of the equation, which is placed in Equation 2. This is because ij E a is added separately in the advection step, which is required by some of the advection methods. This representation follows our implementation in Aspect.
Dominant slip systems for olivine are based on Table 1 in Karato et al. (2008) and Table 1 in Kaminski et al. (2004). The value between the parentheses is normal to the slip plane, while the value between the square brackets is the slip direction (i.e., the direction of the Burgers vector). The Reference Resolved Shear Stress (RRSS) shows the dimensionless RRSS for each slip system in the order of the slip systems as shown in Table 2. ∞ indicates an inactive shear plane. The values in the RRSS column are a combination of values from Table 1 in Kaminski (2002) and Table 1 in Kaminski et al. (2004). is defined as a combination between rotation and deformation of the grain: In Equation 3, index i takes the values of 1, 2, and 3, * ij E L is the dimensionless velocity gradient tensor, defined as L L ij ij Kaminski & Ribe, 2001, Equation 7, see the Supporting Information S1): quantifies the relative activities of the slip systems for each mineral phase (Kaminski et al., 2004). This means that the strength of each slip system has to be computed. This can be done by first computing  sv E I , where  E s is the same as E s but with an arbitrary order of the slip systems (we use the first column of Table 2 but the order is not important since we are going to reorder it). The strength of each slip system can then be computed with: where  sv E I is defined as (below Equation 5 in É. Kaminski & Ribe, 2001): Note. The rotation matrix entries of v E n and v E l can be easily computed by multiplying n * and l * respectively by the inverse rotation matrix (e.g., n a n a n v t   1 * * ). The values in the first column are only used for the index of the rrss vector and on their own do not indicate anything about the strength of the slip system!

Change in Volume Fraction of Grains
Grain size changes are modeled by grain boundary migration, which is driven by the differences in strain energy within the grains and subgrain rotation, which drives the nucleation of new strain-free grains. In the D-Rex algorithm, these changes in grain size are determined by comparing the strain energy of grain to the mean strain energy of the aggregate (adapted from Equation 12; Kaminski et al., 2004): For enstatite there is no experimental data for the mobility. Therefore, as in Kaminski et al. (2004), we use the same value of mobility for enstatite as we use for olivine.
The strain energy of a grain, v E E depends on the strain-rate acting on each slip system. The dimensionless strain energy of the grain *v E E , and can be computed by (Kaminski et al., 2004;Equation 10): where v E E is the strain energy of the crystal in phase E m , ski et al., 2004, Equation 9),  * m E is the dimensionless nucleation rate (Kaminski et al., 2004, Equation 8): where  m E is the nucleation parameter, which determines the rate at which new strain-free grains form from high strain energy grains due to subgrain rotation.  * m E is an input parameter.  *sv E is the dimensionless dislocation density per slip system based on the definition in Equation 11, Kaminski et al. (2004): . This does not make sense since the values are already nondimensional and this is also not what the D-Rex code implements. Also, as has been done in Thissen (2015), we fixed an indexing issue in which summation corresponding to  4v E was not added in the computation of the dislocation density in the D-Rex code. This indexing issue affects the grain size, because the computed dislocation density can be too small in some cases. Note that, as stated in Kaminski et al. (2004), the relationships to compute the 8 of 27 strain energy break down if the density of dislocations is not at a steady state, and therefore D-Rex should not be used for models in which static recrystallization is important.
As stated in Section 2.2.3, small grains (  gbs v E f f ) are assumed to deform by grain boundary sliding and they are free of dislocations. Therefore, the strain energy of these grains is set to zero:  a (9 values) and the volume fraction, v E f , (1 value) for each grain and mineral, for each particle. Therefore, if there are 1,000 grains per particle and two minerals, there are 20,000 output values per particle at each time step. These are the outputs used to create the pole figures. It is also possible to output Euler angles (3 values a , but this is generally not recommended since it is computationally more expensive and can suffer from issues like the gimbal lock (Hemingway & O'Reilly, 2018). At this time the post-processing is limited to the creation of pole figures (see the Supporting Information S1) for comparison with previous results and laboratory experiments, and two different methods for determining the orientation of the fast axis direction: the Bingham average and the maximum eigenvector from the elastic tensor. The user can select when and whether or not to write out the raw information.

Bingham Average
The DCM/rotation matrix for each grain corresponds with the three directions of a rotated coordinate system. These direction are used to create pole figures, but for spatial visualization purposes in the models themselves, a representative direction of all those grains can be useful. To this end we use Bingham statistics (Bingham, 1974;Onstott, 1980) for each of the orthogonal directions of a mineral, which can estimate the center of the densest cluster, similar to a mean direction. This method makes use of the fact that rotation matrices are orthogonal matrices, which means that the rows and columns form a set of orthogonal directions. Using these directions a 3×3 matrix x is one of those vectors and E v is a vector related to a rotation matrix of grain E v of a mineral. Note that the i E p from Equation 4 in Onstott (1980) is set to 1, since we use the same approach to volume weighting as with the pole figures (see Supporting Information S1), so we do not account for that here. The resulting matrix per direction is: The eigenvector corresponding to the largest eigenvalue of this matrix is the estimate for the densest cluster and tends to coincide with the mean vector (Onstott, 1980). This is computed separately for all three orthogonal directions in the rotation matrix providing estimates of the orientation of the crystal axis (a, b, and c) for olivine (and enstatite if included).

Elastic Tensor and Decomposition
We also compute and decompose the elastic tensor resulting from all the enstatite and olivine grains present in the particle. To compute the elastic tensor we use the reference elastic tensor as listed in the D-Rex code 9 of 27 (note this is slightly different than that given in Browaeys and Chevrot (2004), and so we list it in the Supporting Information S1. With the reference elastic tensor for each mineral m pgrs E C , we can rotate each tensor with the rotation matrix E M and compute the Voigt average elastic stiffness tensor (Mainprice, 1990): These computations can also be performed entirely in Voigt notation (Carcione, 2015) which is computationally more efficient: The symmetry cartesian coordinate system (SCCS) of this elastic tensor can be computed. With the three SCCS directions, the elastic tensor can be decomposed into the different symmetries in those three SCCS direction, that is, triclinic, monoclinic, orthorhombic, tetragonal, hexagonal, and isotropic (Browaeys & Chevrot, 2004). The computation of these different symmetries have been implemented and checked against the results in Browaeys and Chevrot (2004).

Results
In this section we present three different types of numerical experiments in which we test the capabilities of the implementation of the CPO algorithm in Aspect. The first set of experiments uses a 3D Cartesian box undergoing simple shear. This set of shearbox experiments is divided in two parts. One part uses a box size and strain-rates on the order of that expected in laboratory experiments, the second part uses a unit box size and strain-rates more typical for a geodynamic setting. These experiments from the second part can be compared to the results in Kaminski et al. (2004) and demonstrate some of the additions we have made compared to what is possible with D-Rex such as computation of the direction of the largest eigenvector of the elastic tensor. The second experiment is a mid-ocean ridge corner flow model, similar to what has been presented in Kaminski et al. (2004) except that it is in 3D. We also use this experiment to demonstrate the fabric selector for stress/water dependence of the CPO fabric. The final experiment is a 3D subduction case in which we show that the CPO implementation is usable in practice for complex geodynamic simulations. For each experiment we note the number of grains (per mineral per particle) used in each experiment, as well as the input parameters used. All the experiments are performed with a 70% olivine and 30% enstatite mix.
Besides these experiments, we have also implemented unit tests, which check the produced derivatives of the code against the analytical result as computed in the Supporting Information S1 and the advection test as presented in Appendix B. All figures use scientific color scales (Crameri, 2018(Crameri, , 2021Crameri et al., 2020).

3D Shear-Box Experiments
We present four different 3D shear box experiments which provide a basic test of the implementation of CPO development allowing for comparison to the previous results from the D-Rex code and expected 10 of 27 results from laboratory experiments. The model domain is a 3D box being sheared in four different tests. All the particles containing the CPO are located in the center of the box (x = 0, y = 0, z = 0). These tests are similar to the tests in Figure 2b in the appendix of MacDougall et al. (2017). The first test is designed to be compared with experimental results in for example Figure 3 in Kaminski et al. (2004). The second test is designed to show that the different olivine fabrics produce the expected CPO. The third test is designed to establish the sensitivity of the results to the number of grains per particle. The last shear-box test demonstrates the sensitivity of the CPO to a change in direction in the flow. s , for the last test.

Tests 1 and 2: Comparing to D-Rex and Experimental Results
This test uses the shear box model to replicate some of the experiments as presented in É. Kaminski et al. (2004) and also compares against laboratory experiments. The flow field is defined as: , creating a simple shear strain field in which the X-Y plane is the shear plane and X is the shear direction ( Figure 1). The test uses a particle in the location where the velocity is zero    ( 0, 0, 0) E x y z . Figure 2 shows the pole figures of both the olivine and enstatite fractions of the CPO separately at  5 10 E t s. The test was run with 5,000 grains and has the same parameters as used in Figure 3 of Kaminski et al. (2004). Both the olivine and enstatite pole figures have the primary and secondary peaks at the correct locations relative to the shear direction.
The second test is designed to show that the implemented algorithm produces reasonable results for each fabric compared to experimental results (e.g., Figure 7 Jung et al., 2006). Therefore, we have repeated the first shear-box test for each of the fabrics, while placing 500 particles at the same location. Using this many particles at a single location (same velocity) allows us to evaluate the range of evolution paths given different random starting orientations. The pole figures are presented for time 5 10 E s, when the integrated strain is 0.5. They show approximately the expected directions for each fabric (right side of Figure 3, compare to Figure 3 from Karato et al. [2008]). Note that we did not reproduce those experiments directly, since we wanted to keep the experiment parameters the same for all fabrics. Because the box is 3D, we can consider how the initially random CPO orientations evolve within the planes that are parallel (X-Y plane) to, and perpendicular (Y-Z plane) to, the shear plane (X-Y) (left side of Figure 3). As expected, the component of the CPO in the parallel plane aligns with the shear direction more quickly (0.5 s) than the component in the perpendicular plane. Furthermore, the direction of the olivine A-axis Bingham average (brown line in Figure 3) is almost indistinguishable from the direction of the largest eigenvector (i.e., the eigenvector associated with the largest eigenvalue) of the elastic tensor (cyan line in Figure 3) in these and the next experiments. This is probably due to a combination of the higher percentage of olivine in each particle and the olivine stiffness matrix is much more anisotropic (Browaeys & Chevrot, 2004). This means that in practice either method can be used to determine the fast direction.

Test 3: Grains per Particle
In this test, we establish some guidelines for how many grains per particle are needed for reliable results. Note that the same number of grains is used for each mineral, therefore a value of 5,000 means 5,000 grains for both olivine and enstatite and the total number of grains being tracked for each particle is 10,000. This test uses the same flow field, particle location, and number of particles as in the fabric test discussed above. The results are shown in Figure 4. For 50 grains, after a run time of  5 3 10 E s, the average directions converge

of 27
to nearly the same values retrieved with more grains, however, there are some grains that are nearly 45° away from the expected direction. Also note that with fewer grains, it takes longer for the range of directions to decrease. For example, with 50 grains, the range in directions is 90° up to nearly  5 1 10 E s, while with 5,000 grains, the range in directions is only about 10° after  5 1 10 E s.
It is difficult to recommend an exact number of grains to use since the potential applications can vary wildly.
If you have a model with a high particle density and you only need a general sense of the direction, then 500 grains may be sufficient. However, in most cases, based on the results shown here, we recommend using at least 1,000 grains. The pole figure of the model with 4 10 E grains is much sharper than the other plots: this is probably both because the data are more aligned and because the way the pole figures are constructed is not well calibrated for so many grains.

Test 4: Change in Shear Direction
Geodynamic flows are expected to be dynamic with rapid changes in the direction of shear especially around slabs (e.g., Jadamec & Billen, 2010). Therefore, we test the influence of a change in shear direction on the CPO evolution and analyze the spread in directions from 500 particles, like in the previous experiments. To make the results of this experiment more relevant for geodynamic simulations, we changed the box size to 100 × 100 × 100 km and lowered the velocity. The model is run for 50 Ma, resulting in an accumulated strain invariant of 2.5. To begin the shear direction is in the X direction in the X-Y plane, similar    Figure 3 in Kaminski et al. [2004]). Top row are the pole figures for olivine. Bottom row are the pole figures for enstatite. Single particle with 5,000 grains. Arrows indicate shear direction (X) in the X-Y plane. Colors are multiples of a uniform distribution plotted as described in the Supporting Information S1.

of 27
The first thing to note is that the orientation of the olivine and enstatite a-axis do not respond to the change in shear direction in the same way. Second, the change in orientation within the X-Y plane is different than the change in orientation in the Y-Z plane. The olivine a-axis rotates relatively quickly with a small spread in directions to align with the new direction, while the enstatite direction takes longer to align to the new direction and shows a much larger spread of orientations in the process. This is because enstatite only has one slip system, and therefore it is more difficult to re-align. When there is a switch in shear directions, the orientation of the enstatite a axis in the X-Y plane and the olivine a-axis in the Y-Z plane actually shows a reversal away from the shear direction, before rotating back toward an angle of 0°. Also, although most of the time the olivine a-axis and the orientation of the largest elastic eigenvector are aligned, they separate by a few degrees following the change in direction of the flow field. Since the largest elastic eigenvector represents the resulting fast direction from both the olivine and enstatite, we will preferentially use this in the next experiments. This experiment also demonstrate the importance of tracking the CPO in 3D, since the out-of-plane rotation affects the time-scale for re-aligning the CPO as the shear direction changes. Similar behavior has been shown in experiments (e.g., Boneh et al., 2015). . The plots are made with 500 particles each containing 5,000 grains starting with an initial random orientation and placed in the same position in space. Right side: Pole figures from a single particle at time 5 1 E . Since the strain-rate is experienced by the particles is    6 1 5 10 s E , the integrated strain at this time is 0.5. Each row is for a different crystal preferred orientation fabric (A-E). These can be loosely compared with Figure 3 of Karato et al. (2008) or Figure 7 in Jung et al. (2006). For details on how the pole figure is created, see the Supporting Information S1. Left side: two graphs are showing the CPO evolution over time using the angle on the shear plane between the a-axis direction and the shear direction X, and the angle between the a-axis direction and the shear direction (X) in the plane perpendicular to the shear plane (the X-Z plane). The a axis direction for each particle for both the olivine (brown) and enstatite (magenta) is determined using the Bingham average. A third curve is the orientation of the largest eigenvector of the symmetry cartesian coordinate system elastic tensor. The solid lines give the average value for all the particles, while the shaded area shows the spread (the area containing all the values).

Mid-Ocean Ridge/Corner Flow
This test emulates the mid-ocean ridge spreading example in (Kaminski et al., 2004), with the addition of the automatic fabric selector and two regions with different water content. The model domain is 300 × 300 × 25 km and contains 750 randomly distributed particles that each contain 1,000 grains. The boundary at 0 km is the location of the ridge axis and this boundary and the front and back boundaries are free-slip boundaries. The top boundary has a fixed velocity of 1 cm/yr, while the remaining two sides have open boundaries. The temperature is constant (i.e., this is a purely mechanical model, with no thermally defined lithosphere) with an isoviscous viscosity of 5 22 E e . The water is distributed such that above the line . The conditions and plots are the same as described in Figure 3 except each row is a simulation using a different number of grains per particle, as indicated.
14 of 27    0.5 * 300000 E z x the water content is 25 ppm H/Si, while below it the water content is 800 ppm H/Si. The results can be seen in Figure 6. The CPO directions that we refer to in this section are always the eigenvector belonging to the largest eigenvalue of the SCCS elastic matrix.
Although the values in this experiment have not been chosen to be representative of a natural ridge system, it is important to note the dynamic nature of the CPO in this case. For example, at 10 Ma in the top right corner (furthest from the ridge axis) the patterns of the CPO are almost dipping a 45° angle in the vertical plane (side view), while in the horizontal plane (top view) most of the CPO aligns very quickly with the flow. Likewise, with increasing distance from the ridge axis, fewer of the CPO at the top of the model align at a 90° angle from the flow. This can be explained by activation of the B-type fabric (high stress and moderate water content) as the flow at the top of the model leaves the ridge axis. Over the course of the model run (e.g., 20 to 30 Ma), this B-type fabric becomes more pronounced, i.e., more grains in each particle have a CPO in this preferred direction. This is because these particles have stayed longer in the B-type fabric stability field than the particles that were at the corner at the start of the model. This makes it progressively harder for the CPO to be reset to the A and D-type fabric even though the stress is lower further from the ridge axis.

of 27
We want to emphasize that these results should not be interpreted to mean that his alignment is predicted for real mid-ocean ridge settings: we have not made any effort here to model the thermal structure, true rheology, hydration, or stress state of a mid-ocean ridge. Rather, this model is used to illustrate how the fabric selector works within a region of changing flow direction and water content and that switching alignment once there is a strong CPO depends on the strain rate.

Subduction
To test the feasibility of using the CPO calculation in real geodynamic models, we added the CPO calculation to a more complex and realistic geodynamic setting than a corner flow model: a simple 3D subduction zone. The model uses many advanced features such as open boundary conditions, complex rheologies (Glerum et al., 2018), and nonlinear materials and solvers (e.g., Fraters, 2019), but the tectonic setting is relatively simple: a subduction zone that is orthogonal to the convergence direction. The model box has a width of 2,500 km (x), a length of 2,000 km (y), and is 800 km deep (z). The top boundary has a free slip boundary condition and the bottom boundary has a zero velocity boundary condition. The sides are all divided vertically where the top part has a velocity boundary condition where some directions are left free and the bottom has a traction boundary condition which is set to the initial pressure at the sides. For further  17 of 27 details, see the input files provided in the Supporting Information S1. The trench is placed in the center of the y axis and is 1,000 km wide. On both sides of the slab and overriding plate there are lithospheric weak zones. In this experiment we use 13,500 randomly distributed particles (on average 1 per The fabric selector is used to determine the fabric of the olivine. Since the fabric selector requires water to determine the fabric type, the model contains a wet mantle, a dry lithosphere (top 100 km and slab), and a very wet top layer (top 30 km). The exact amount of water at each location is not important to showcase the CPO implementation and would require further investigation before using models to simulate real settings. The model parameters have also not been tuned to represent a specific subduction zone or to match any observables. The main goal is to show that a relatively simple 3D subduction zone setting can create a complex 3D field of CPO directions. For more details on the model setup and parameters see the Supporting Information S1. Even with these precautions about interpreting the results for this model, some interesting patterns arise illustrating how the complex dynamics around a slab might manifest in terms of the CPO pattern.  figure). The CPO directions that we refer to in this section are always the eigenvector belonging to the largest eigenvalue of the SCCS elastic matrix. This means that it takes into account both the olivine and enstatite grains in the particle and their stiffness matrices. Figure 8 combines the velocity field and CPO directions at two different depth slices at 30 Ma. The top part of the figure represents the model results at a depth between 200 and 300 km and the bottom part of the figure represents the model results at a depth between 300 and 400 km. The size of the CPO cylinders is determined by the integrated strain invariant, with a max at 0.5. The integrated strain invariant is computed by multiplying the second invariant of the strain-rate tensor by the timestep size and adding the result to the currently integrated strain invariant (starting a 0 at the start of the model). This approach is not very accurate, but we only use it to make sure that we only plot directions that have actually experienced enough strain to get a reliable direction.
At 10 Ma the CPO directions as seen from the top view in Figure 7 are mostly in the direction of the flow. The only obvious exception is the CPO directions around the edges of the slab, which show the effect of toroidal flow around the slab edge driven by the sinking slab. The CPO in this region curves around the slab edges, while the flow is pulled away from the edges by the motion of the unsubducted parts of the subducting plate adjacent to the overriding plate. This is a consistent pattern that can be seen during the whole model run and the difference between the CPO and flows direction is especially visible at a depth of 300-400 km (see slab edges in Figures 8a and 8c).
In the cross section of Figure 7 at 10 Ma, we see that the CPO orientation behind the slab is in the direction of the flow, but in front of the slab, it is following the sinking slab instead of the flow direction. The fabric is mostly A-type, with some E-type, and C-type fabric in the trench. As noted before, we have not used realistic water distribution in this model, therefore the fabric types should not be over-interpreted, but rather illustrate the range of a variable fabric possible. At 20 Ma a new CPO direction starts to become visible below the overriding plate and keeps growing at least up to 40 Ma. This new CPO direction is parallel to the trench and is formed in particles where the olivine A-fabric is active. The location of these trench parallel CPO directions is in the cross section at 40 Ma, located between 400 and 1,000 km from the edge of the model at the overriding plate side and at a depth of between 200 and about 350 km (see the black box in Figure 7).
With the precautions noted above, we can make some interpretations. First, the CPO around the edges of the slab is interpreted as a consequence of a combination of the toroidal flow around the edges of the slab, driven by the sinking slab, and dragging of the mantle by the un-subducted portions of the

of 27
subducting plate adjacent to the overriding plate. These two components of the flow field cause shear strains adjacent to the slab edges, which are dominated by the toroidal flow on the back side of the slab, but not on the front side of the slab. Second, the CPO directions below the overriding plate, which are aligned parallel to the trench, are interpreted as a consequence of the component the shear strain in the mantle caused by flow entrained by the moving portions of the subducting plate adjacent to the fixed overriding plate. In the top view, the flow beneath this region flows away from the center at an angle of around 45°, which results in a shear direction parallel to the trench. Third, the CPO near the trench is strongly aligned with the sinking direction of the slab. The dipping orientation of this CPO would strongly influence the hypothetical seismic anisotropy measurements made at the surface. Specifically, when modeling seismic anisotropy measurements, it is usually assumed that the fast axis is aligned in the horizontal plane and is not dipping (Savage, 1999). Therefore, using the output of these kinds of models to forward predict the expected anisotropic signal will greatly improve interpretations of seismic anisotropy data. Finally, note that there are no B and D-type fabrics present in the model. This is either due to the stresses not being high enough to generate these fabrics and/or due to a lack of resolution in the mantle wedge to resolve this.

Discussion
By adding the CPO module to the open-source code Aspect, we have now made it possible for any Aspect user to add the CPO calculation to a geodynamic simulation. While a few studies have used codes with an integrated CPO calculation (e.g., Faccenda & Capitanio, 2012;Faccenda & Capitanio, 2013), many studies have either used some other approximation of CPO (e.g., Jadamec & Billen, 2010;Kneller & van Keken, 2008;Lev & Hager, 2011;Lin, 2014;Lin & Kuo, 2013) or analyzed CPO as a post-processing step (e.g., Di Lin, 2014;Lassak et al., 2006;Morishige & Honda, 2010). While calculating CPO as a post-processing step may be cost-effective for steady-state models in which only a single CPO calculation is needed, the cost quickly becomes prohibitive for time-dependent models in which the evolution of the CPO is of interest.
The addition of the CPO calculation to a fully parallelized geodynamic code greatly reduces the overhead cost of including CPO analysis in geodynamic models. However, adding the CPO calculation to models in Aspect can still add significantly to the overall computation cost. For example, for the 40 Ma subduction example, updating and advecting of the 13,500 particles, which contained the CPO (each having 1,000 olivine and 1,000 enstatite grains), the CPO Bingham average, the elastic tensor and its decomposition, and the integrated strain invariant, took about 20% of the total computation time. Although we are confident that CPO calculation cost can be improved with further development, it is also true that users can reduce this cost by carefully choosing which output values will be needed for analysis. In some cases, only the olivine and enstatite directions are of interest, or in other cases, there is only interest in CPO in a certain area of the model: then overhead could be reduced significantly. Nonetheless, while adding CPO does add a significant overhead, the added benefit of being able to use seismic anisotropy as a constraint on the simulations for many more studies will outweigh this added cost.
Besides the addition of passive CPO to many geodynamic models, one other major future research avenue would be to actively link the CPO to the viscosity calculation in the model. Although Aspect has capabilities to model anisotropic viscosities, more development remains to actively link these two physical processes in the code. Another future development that could be added as a post-processing step is the direct computation of seismic anisotropy from the modeled CPO. The current implementation allows for the output of both the raw and volume-weighted CPO data, which could be used for this purpose in a post-processing step. We were purposeful in our implementation to not make any assumptions about crystal symmetry (Savage, 1999) so that this choice (or multiple analyses with different choices) could be made independent of the CPO calculation.
The CPO calculation itself also still has room for improvement and research potential. The fabric model we have used for the fabric selector is simple digitization (see the figure in Appendix C) of the plot of fabric dependence on water and stress from Karato et al. (2008). Because the experiments in that figure do not fully constrain the boundaries between different fabrics, one should assume that there are large 20 of 27 error bars on the transition region from one fabric to another. Furthermore, the water field in Aspect is a passively advected field. To make it more realistic, and thereby making the CPO calculation more realistic, a model for the distribution and evolution of water in the mantle needs to be developed. However, once that is done, no additional change to the CPO module is needed to account for the effects of water.
Finally, the CPO framework now included in Aspect greatly facilitates future improvements to the D-Rex algorithm or the addition of a completely different CPO algorithm. This is because there is a strict separation between several parts of the code. The D-Rex algorithm is only responsible to compute per grain  E and  E (see Equations 2 and 11) and the infrastructure allows to easily exchange what code computes those values through the input file. Besides D-Rex there is currently also the option to let the CPO passively rotate with the flow. With this flexibility, future work could also link with grain size evolution models, as for example has been done in Dannberg et al. (2017) using the paleowattmeter (Austin & Evans, 2007) or (Rozel et al., 2011), and anisotropic rheology (e.g., Lev & Hager, 2011). There are two main drawbacks of the grain size evolution in D-REX that make linking the grain size distribution in D-Rex to other material properties such as viscosity challenging: the grain size is non-dimensional (with no clear way to dimensionalize) and the grain size evolution does not depend on stress (or mechanical work; i.e., Austin & Evans, 2007). Future work would need to assess whether the grain size distribution from D-REX can be scaled to dimensional values and whether the resulting distribution (and mean grain size) was similar to other grain size evolution models (e.g., Austin & Evans, 2007;Rozel et al., 2011). The way the code is set up also allows for easy output of the CPO data (directions and sizes) and the full elastic tensor. This allows it to be used in future research for computing an "observed" seismic anisotropy based on the model results.

Conclusions
Seismic anisotropy is a powerful and yet under-utilized observational constraint on geodynamic models. By implementing the D-Rex algorithm in the open-source, fully parallelized geodynamic simulation code Aspect, we have removed two of the major barriers to utilizing seismic anisotropy constraints more often: access and computational cost. This should facilitate more studies on the development of anisotropic structure in different tectonic environments and parts of the mantle, and the interpretation of observations of seismic anisotropy in terms of mantle flow. Similarly, we expect that improvements to the D-Rex algorithm and calibration of the algorithm against experimental data will be facilitated by its integration into Aspect. Finally, we have laid the groundwork for future work that could integrate the CPO evolution with grain size evolution models and anisotropic rheology. The implementation in Aspect has been carefully tested, and we have presented several cases to show that it works and that it is cost-effective to use in realistic geodynamic models on a regular basis.

Appendix A: Tables of D-Rex Related Symbols
This appendix contains two tables. The first table list the indices used in the equations. The second table lists all the variables related to the D-Rex calculation with both their representation in this paper and in the original D-Rex papers (É. Kaminski & Ribe, 2001;É. Kaminski et al., 2004) (Table A1).
x dt x x (B6) After the advection step we normalize the grain sizes back to a sum of one and project ij E a back to an orthogonal tensor following the implementation in Thissen (2015). Note. Sb is the symbol used in this paper, DR is the equivalent symbol in the É. Kaminski et al. (2004) or É. Kaminski and Ribe (2001). The variables are ordered such that the variables in the definition are defined above the current entry. Some values are not defined, because they only appear in the definitions of input parameters. Note that if a variable is dependent on index E v it is automatically assumed to also be dependent on index E m . This table is ordered such that is can be read as an algorithm.

of 27
To test the advection schemes we put a particle at position (0,1) with a rotation matrix of and taking the square root of the squared sum of the elements of that result. For the particle advection, we have chosen to use the fourth order Runga-Kutta scheme. This prevents the size of the error from the particle advection being dominant. The results are shown in Figure B1. The convergence order for the Crank-Nicolson scheme is second-order, as expected. The convergence order of the Forward and Backward Euler schemes is also second-order: one order higher than expected. This is probably due to the test being relatively simple given that the spin tensor is constant in time Figure C1. Figure B1. convergence rates (rotation matrix with respect to timestep size) of the different implemented advection algorithms.