the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
ECOMAN: an opensource package for geodynamic and seismological modelling of mechanical anisotropy
Manuele Faccenda
Brandon P. VanderBeek
Albert de Montserrat
Jianfeng Yang
Francesco Rappisi
Neil Ribe
Mechanical anisotropy related to rock fabrics is a proxy for constraining the Earth's deformation patterns. However, the forward and inverse modelling of mechanical anisotropy in 3D largescale domains has been traditionally hampered by the intensive computational cost and the lack of a dedicated, opensource computational framework. Here we introduce ECOMAN (Exploring the COnsequences of Mechanical ANisotropy), a software package for modelling strain and stressinduced rock fabrics and testing the effects of the resulting elastic and viscous anisotropy on seismic imaging and mantle convection patterns.
Differently from existing analogous software, ECOMAN can model straininduced fabrics across all mantle levels and is optimised to run efficiently on multiple CPUs. It also enables modelling of shape preferred orientation (SPO)related structures that can be superimposed over lattice/crystallographic preferred orientation (LPO/CPO) fabrics, which allows the consideration of the mechanical effects of fluidfilled cracks, foliated and lineated grainscale fabrics, and rockscale layering.
One of the most important innovations is the Platform for Seismic Imaging (PSI), a set of programs for performing forward and inverse seismic modelling in isotropic–anisotropic media using real or synthetic seismic datasets. The anisotropic inversion strategy is capable of recovering parameters describing a tilted transversely isotropic (TTI) medium, which is required to reconstruct 3D structures and mantle strain patterns and to validate geodynamic models.
 Article
(10328 KB)  Fulltext XML

Supplement
(4395 KB)  BibTeX
 EndNote
The study of the Earth's interior has been traditionally based on seismological and geodynamic modelling, the former providing important information about its presentday structure, composition, and state (Chang et al., 2015; French and Romanowicz, 2015; Schaeffer et al., 2016; Debayle et al., 2020) and the latter about its dynamics and compositional evolution (Davies et al., 2012; Crameri and Tackley, 2014; Müller et al., 2022). Seismological modelling and geodynamical modelling are very often conducted independently, which creates mechanical and geometrical inconsistencies across the models, hampers the interpretation of seismic observations in terms of geodynamic processes, and exacerbates the nonuniqueness of geodynamic model predictions. An alternative approach is combining computational seismology and geodynamics with mineral physics, which provides a comprehensive understanding of the Earth's interior processes, seismic behaviour, and material properties. This multidisciplinary methodology has been used in previous studies to postprocess geodynamic flow calculations with thermodynamically selfconsistent models of mantle mineralogy and convert thermal structure into isotropic elastic parameters. The obtained seismic mantle structure can then be used in simulations of global wave propagation such that specific hypotheses on mantle dynamics can be tested directly against seismic data (Styles et al., 2011; Schuberth et al., 2012; Maguire et al., 2018). The inverse procedure consists of converting seismic anomalies into density anomalies driving mantle flow models and is typically employed to quantify dynamic topography and mantle viscosity structure, as well as to reproduce largescale mantle flow patterns (e.g. Bunge et al., 2003; Steinberger and Calderwood, 2006; Rudolph et al., 2015). However, isotropic seismic imaging provides limited information regarding local and regionalscale dynamical processes (Fraters and Billen, 2021), and a better way to couple the geodynamic evolution and seismological structure of the Earth's interior is by accounting for the strain and stressinduced mechanical anisotropy of crustal and mantle rocks.
Mechanical anisotropy refers to the directional dependence of mechanical properties in a material and is well known to affect both elastic and viscous deformational behaviour. Mechanical anisotropy depends on several factors, including the lattice or crystal preferred orientation (LPO, CPO) of intrinsically anisotropic minerals and extrinsic mechanisms related to the shape preferential orientation (SPO) of melt, fluid, or airfilled fractures and nonspherical pores, as well as grain or rockscale compositional layering. Most micro and macroscale fabrics are acquired as a function of the cumulative deformation and material mechanical properties, and as such they constitute an important source of information about the Earth's dynamical behaviour.
Elastic anisotropy is directly connected to seismic anisotropy, which is a phenomenon in which the seismic wave speed varies as a function of the propagation direction. It is mainly observed in the crust, mantle boundary layers, and inner core (Almqvist and Mainprice, 2017; Karato, 1998; Kendall, 2000; Long and Becker, 2010; Deuss, 2014). Understanding and modelling seismic anisotropy are crucial for determining longterm deformational patterns and the presentday stress field in the crust, as well as constraining geodynamic modelling predictions (Jadamec and Billen, 2010; Hu et al., 2017; Zhou et al., 2018; Lo Bue et al., 2022). In addition, the ability to account for anisotropic effects can improve the quality of subsurface imaging. Indeed, it has been demonstrated that, because of uneven seismic ray coverage, failing to account for seismic anisotropy may generate strong artefacts that substantially bias our understanding of mantle structures and dynamics in different tectonic settings (Bezada et al., 2016; VanderBeek and Faccenda, 2021; VanderBeek et al., 2023; Faccenda and VanderBeek, 2023). Considering the widespread presence of seismic anisotropy, anisotropic seismic models provide a more realistic representation of the Earth's subsurface compared to isotropic models.
Viscous anisotropy modelling refers to the study and simulation of materials that exhibit varying degrees of viscosity (resistance to viscous deformation) in different directions. Although viscous anisotropy has been traditionally associated with the mechanical behaviour of multilayered media (Mühlhaus et al., 2002; Kocher et al., 2006), it has been also observed in experimentally deformed micarich and olivine crystal aggregates (Shea and Kronenberg, 1993; Hansen et al., 2014). Previous numerical studies demonstrated that viscous anisotropy can potentially stabilise longwavelength convective patterns (Christensen, 1987; Mühlhaus et al., 2004) and more generally affect processes such as plate motion (Király et al., 2020), postglacial rebound (Han and Wahr, 1997), lithospheric shear zone reactivation (Tommasi et al., 2009), and dripping (Lev and Hager, 2008a).
Over the last few years a few attempts have been made to integrate micromechanical modelling of fabric evolution with largescale geodynamic models using either directors (Lev and Hager, 2008b; Halter et al., 2022), the CPO model DREX (Kaminski et al., 2004; Becker et al., 2006; Jadamec and Billen, 2010; Faccenda and Capitanio, 2013; Faccenda, 2014; Ito et al., 2014; Hu et al., 2017; Zhou et al., 2018; Fraters and Billen, 2021), or the CPO model viscoplastic selfconsistent (VPSC; Tommasi et al., 2009; Li et al., 2014). However, each of these methodologies has its own limitations mainly associated with the accuracy of the estimates, the large computational burden, or software accessibility, which have impeded a more widespread diffusion in the geodynamic community. At the same time, the recovery of 3D seismic anisotropy patterns has been traditionally considered intractable due to the highly underdetermined nature of the inverse problem, and, although it has been addressed by several previous studies (e.g. Debayle et al., 2005; Panning and Romanowicz, 2006; Abt and Fischer, 2008; Long et al., 2008; Wookey, 2012; Mondal and Long, 2019), only recently have a theoretical background and computational algorithms been developed to simultaneously invert for P and/or Swave isotropic velocity anomalies and anisotropy resulting from arbitrarily oriented structures (Munzarova et al., 2018; VanderBeek and Faccenda, 2021; Rappisi et al., 2022; Wang and Zhao, 2021; VanderBeek et al., 2023; Del Piccolo et al., 2024). Yet, to date there is no freely available software capable of modelling seismic anisotropy related to arbitrarily oriented structures.
In this contribution we present the new and opensource software package ECOMAN (Exploring the COnsequences of Mechanical ANisotropy) that enables linking seismology and geodynamics by providing a set of computationally optimised programs for (i) estimating rock mechanical anisotropy as a function of the geodynamic model deformation history and compositional, rheological, stress, pressure, and temperature fields and (ii) solving forward and inverse seismological problems accounting for seismic anisotropy. In the next sections, we first describe the different ECOMAN modules, after which we discuss the advantages and limitations of the software package and the roadmap for future developments.
ECOMAN includes several programs that are complementary and can be grouped into three main categories (Fig. 1):

programs that estimate strain and stressinduced rock fabrics (LPO and SPO) and their elastic and viscous anisotropic mechanical properties (DREX_S, DREX_M, EXEV);

programs that postprocess the simulated rock fabrics for visualisation of their isotropic–anisotropic mechanical properties and deformational history (VIZTOMO, VIZVISC) and format the elastic tensors generating input files for seismological synthetics (VIZTOMO); and

programs that test the elastic response of anisotropic media by performing seismological forward and inverse modelling and, in particular, isotropic and anisotropic seismic tomographies on synthetic and real seismic datasets (SKSSPLIT, PSI).
Most of the code modules are written in the Fortran programming language, except for the PSI program that is written in Julia. Visualisation of the output is done through the MATLAB MTEX toolbox (Mainprice et al., 2011) for singleaggregate fabrics and ParaView (Ahrens et al., 2005) for 2D and 3D simulations.
2.1 Rock fabric and mechanical property simulations
The evolution of the straininduced LPO can be simulated with the DREX model (Kaminski et al., 2004), which has been adapted to reproduce the fabric evolution of single (DREX_S) or multiple (DREX_M) twomineralphase aggregates representative of the entire Earth's mantle. The accuracy of the DREX model has been tested against analytical solutions derived by Fraters and Billen (2021) (Appendix A). Five sets of twomineralphase mantle aggregates can be defined as a function of depth or density ρ (Fig. 2; see Table 1 for a list of abbreviations and physical properties):

olivine + enstatite for the upper mantle (UM: 0–410 km or $\mathrm{3000}<\mathit{\rho}\le \mathrm{3650}$ kg m^{−3});

wadsleyite + majoritic garnet for the upper mantle transition zone (UTZ: 410–520 km or $\mathrm{3650}<\mathit{\rho}\le \mathrm{3870}$ kg m^{−3});

ringwoodite + majoritic garnet for the lower mantle transition zone (LTZ: 520–660 km or $\mathrm{3870}<\mathit{\rho}\le \mathrm{4150}$ kg m^{−3});

bridgmanite + ferropericlase for the lower mantle (LM: 660–2900 km or ρ>4150 kg m^{−3});

postperovskite + ferropericlase for the bottom of the lower mantle according to the parameterised phase boundary $P\left(\mathrm{GPa}\right)=\mathrm{98.7}+T\left(\mathrm{K}\right)\cdot \mathrm{0.00956}$ (Oganov and Ono, 2004).
The straininduced LPO is computed for the phases in bold, while other major phases such as garnet, ringwoodite, and ferropericlase are considered to be isotropic and their distribution is set to be random. Thus, no LPO is calculated for the LTZ such that (minor) anisotropy arises only when SPO modelling due to compositional layering is active (Sect. 2.1.3; Faccenda et al., 2019). The LPOs of wadsleyite, bridgmanite, and postperovskite are computed according to the same DREX methodology originally applied to uppermantle phases and described in detail in Kaminski and Ribe (2001) and Kaminski et al. (2004) but including additional slip systems as indicated in Table A1. For each of these new phases it is possible to define in the DREX_S and DREX_M input files the same free parameters of the DREX model (the dimensionless grain boundary mobility M^{*}, nucleation parameter λ^{*}, threshold volume fraction χ, and in addition the powerlaw exponent for dislocation creep deformation within the grains) as for uppermantle aggregates. The full elastic tensor is then calculated according to the crystal orientation, volume fraction, phase abundance, P–T conditions, and bulk rock composition, as well as using Voigt–Reuss–Hill averaging schemes (see Sect. 2.1.2 and Appendix B for more details).
The elastic properties related to strain and stressinduced SPO fabrics can instead be calculated at the grain or rock scale and for layered or twophase (matrixellipsoidal inclusions) systems using the isotropic elastic moduli of the different (fluid, mineral, rock) components (EXEV).
The elastic properties and density of the aggregates characterised by LPO and/or SPO fabrics are estimated at relevant mantle P–T conditions using the singlecrystal elastic moduli and their P–T derivatives for the main mineral phases and compiled from different mineral physics studies, together with lookup tables of the isotropic elastic moduli, density, and mineralphase volume fraction generated with MMA_EoS (Chust et al., 2017) for five bulk rock compositions (dunite, harzburgite, pyrolite, MORB, pyroxenite) at temperature and pressure ranges of $T=\mathrm{300}:\mathrm{50}:\mathrm{4500}$ K and $P=\mathrm{0}:\mathrm{0.1}:\mathrm{140}$ GPa, respectively (see Appendix B and Tables B1 and B2).
2.1.1 DREX_S
DREX_S is a program designed for modelling the evolution of straininduced LPO fabrics and related elastic properties of a single, twomineralphase mantle aggregate as a function of the imposed flow field, amount of strain, crystal plasticity, P–T conditions, and additional effects related to SPO fabrics. It builds on the original DREX software (Kaminski et al., 2004) for modelling the straininduced LPO, and it includes MATLAB scripts to generate pole figures of the LPO and isotropic–anisotropic seismic properties with the MTEX software (Mainprice et al., 2011) (Fig. 3a–d), together with the possibility to display the evolving fabric strength (M index, J index) and the fraction of different anisotropy components (obtained via tensor decomposition; Browaeys and Chevrot, 2004; Fig. 3e–h).
DREX_S is particularly useful for users who are not familiar with LPO modelling and, more in general, anyone interested in performing parameter sensitivity tests on different mantle mineral aggregates before launching largescale simulations. In addition, the microstructures generated with DREX_S can be used in the DREX_M 2D–3D simulations to impose preexisting (e.g. fossil) fabrics on multiplecrystal aggregates located within a specific subdomain (see Sect. 2.1.2).
As an application and with the aim of debunking the common misconception that the DREX model can only produce LPO fabrics that are too strong, here we show that with a suitable choice of DREX free parameters it is possible to reproduce both weak and strong olivine fabrics as recorded in natural samples (Warren et al., 2008) and highstrain laboratory experiments (Hansen et al., 2014; Tasaka et al., 2017) (Fig. 4).
2.1.2 DREX_M
DREX_M is a program that computes the evolution of the LPO and related elastic properties of multiple, twomineralphase mantle aggregates as a function of the singlecrystal plastic and elastic properties and of the flow field, deformation mechanisms, and P–T conditions resulting from 2D–3D geodynamic simulations. It builds on the original DREX software, which includes routines for estimating the straininduced LPO and elastic properties (i) for uppermantle polycrystalline aggregates only, (ii) using singlecrystal elastic tensors derived at room P–T conditions and averaged using a Voigt scheme, (iii) in a 2D Cartesian domain, (iv) for a single (steadystate) flow field, and (v) whereby the whole deformation is considered to be accommodated by dislocation creep assisted by grain boundary sliding (Kaminski et al., 2004). DREX_M additionally models the following.

Fabrics relevant to the middle and lowermost mantle (Faccenda, 2014), including those with postperovskite. Phase transitions can be set to occur at predefined depths (e.g. 410, 660 km), density crossovers (which allow modelling the deflection of phase boundaries with nonzero Clapeyron slopes; Fig. 7), and parameterised phase boundaries as for the case of postperovskite (Oganov and Ono, 2004). When a particle enters the stability field of another anisotropic phase, the LPO can either be reset or retained (which implies axisymmetric topotactical growth).

Elastic properties and density as a function of the bulk rock composition and local P–T conditions (Faccenda, 2014; Chang et al., 2016; Ferreira et al., 2019; Sturgeon et al., 2019) (Appendix B). The isotropic component of the elastic tensors and density are taken from the lookup tables generated by MMAEoS for a given mantle lithology (Table A1), as well as the anisotropic component from the mineral singlecrystal elastic moduli and their pressure and temperature derivatives listed in Table B2. This strategy ensures a gradual transition of the seismic properties at phase boundaries where phase transformations occur. Voigt–Reuss–Hill averaging schemes of the elastic moduli are included.

Nonsteadystate flows in 2D–3D Cartesian and polar grids (Faccenda and Capitanio, 2012, 2013; Hu et al., 2017; Zhou et al., 2018; Lo Bue et al., 2022; Faccenda and VanderBeek, 2023; Rappisi et al., 2024). In polar coordinates the velocity gradient tensor must be computed in the external Cartesian reference frame, as described in Appendix C. The globalscale models are spatially discretised using socalled Yin–Yang grids (Kageyama and Sato, 2004). Several examples (cookbooks) are provided on how to use the software in different coordinate systems and in steadystate or timedependent flow conditions.

Fabric evolution in the presence of multiple creep mechanisms. At any time step, the fraction of deformation accommodated by dislocation creep in a given point of the geodynamic model defines the fraction of time spent for intracrystalline deformation assisted by grain boundary sliding according to the DREX model. The remaining time is used to apply, when present, fluid deformation rotation to the whole crystal aggregate (e.g. Hedjazian et al., 2017) (Appendix D).

A preexisting (fossil) fabric (precomputed with DREX_S) within a subdomain, typically the lithosphere. This is often the case for geodynamic models where the lithosphere accretion is not modelled, and its geometry is initially prescribed.
The DREX_M input files should contain information about the geodynamic model evolution. Critical information includes the components of the velocity vector V field and, for timedependent flow models, the total elapsed time and time step. These are used to compute the velocity gradient tensor, as well as the LPO evolution and advection of the crystal aggregates. Additional fields defined on the Cartesian–spherical grid that can be included in the input files are the following.

When the geodynamic model is thermomechanical, the temperature T and total pressure P fields, which can be used to compute the singlecrystal elastic tensor and aggregate phase transitions as a function of the local P–T conditions.

When the rheological model of the geodynamic simulation is based on multiple viscoplastic deformation mechanisms, the fraction of deformation accommodated by dislocation creep ${F}_{\mathrm{d}}={\mathit{\eta}}_{\mathrm{disl}}/{\mathit{\eta}}_{\mathrm{eff}}$, where $\mathrm{0}\le {F}_{\mathrm{d}}\le \mathrm{1}$, η_{disl} is the viscosity calculated with the dislocation creep flow law, and η_{eff} is the effective viscosity calculated with the harmonic average of each of the viscosities representing a different deformation mechanism.
In summary, the time and the velocity field V are essential pieces of information, while the P, T, and F_{d} fields are optional and depend on the type of (mechanical vs. thermomechanical) geodynamic and (single vs. multiple viscoplastic deformation mechanisms) rheological models.
While the V, P, T, and F_{d} fields are defined on the Eulerian grid, the distribution, size, modal composition, and mechanical properties of mantle aggregates are defined on the Lagrangian particles. After initialising the Eulerian grid and Lagrangian particles, the entire run consists of three main steps: (i) backward advection of the particles for a given time span or up to a maximum cumulated strain, (ii) forward advection and update of the LPO and deformation gradient tensor F, and (iii) computation of the full elastic tensor and creation of the output file. This strategy ensures a homogeneous final distribution of the aggregates (coincident with the initial one) which is desired for visualisation of the aggregates' mechanical and physical properties and for running seismological synthetics. The DREX_M output file includes, for each mineral aggregate, the elastic tensor C, density, and deformation gradient tensor F. These properties can then be processed and visualised with the software VIZTOMO (Sect. 2.2.1).
DREX_M is parallelised using a hybrid MPI and OpenMP scheme to take advantage of multipleCPU architectures of modern HPC clusters. Runtime is obviously affected by the number of crystal aggregates, number of crystals per aggregate, number of active slip systems per crystal, and number of time steps (Fig. 5a). The parallel efficiency is 70 %–90 %, and the update of the LPO and F tensor during forward advection is the most timeconsuming part of the run (Fig. 5b and c). The small performance degradation is due to the initialisation of the Eulerian–Lagrangian grids and arrays, as well as to I/O operations which are executed serially within each process. As a result, the efficiency of the timedependent flow models is lower than that of steadystate models, as the latter only require a single velocity (and P, T, and F_{d}) field to be loaded and processed.
2.1.3 EXEV
EXEV includes routines to compute the extrinsic elastic and viscous anisotropy using effective medium theory modelling for a multicomponent layered system (smoothed transversely isotropic longwavelength equivalent, STILWE; Backus, 1962) or a twocomponent system with similar ellipsoidal inclusions in a uniform background matrix (differential effective medium, DEM; e.g. Mainprice, 1997).
The elastic tensor C due to SPO fabrics can either be estimated independently or, when using the DEM approach, superimposed on that obtained from the straininduced LPO modelling. SPO fabrics that can be modelled are those related to rock or grainscale layering (e.g. Faccenda et al., 2019) or to the presence of preferentially aligned ellipsoidal inclusions (e.g. melt and/or fluidfilled cracks). The user then needs to specify the following.

For grainscale layered fabrics, a dominant ultramafic or mafic lithology. In this case the mineralphase proportions from the MMaEoS lookup tables define the mixture for the STILWE model.

For rockscale layered fabrics, the relative abundance of the five available ultramafic–mafic lithologies (dunite, harzburgite, pyrolite, basalt, and pyroxenite) defining the mixture for the STILWE model.

For matrixinclusion fabrics, the elastic tensors of the two components and the inclusion's shape and volume fraction as required by the DEM modelling. The matrix elastic tensor can be replaced with that from the LPO modelling in order to estimate the combined effect of LPO and SPO fabrics.
The SPO fabrics can then be oriented at any angle relative to the principal finitestrain ellipsoid (FSE) axis or, in the case of cracks, the local principal stress obtained from the “presentday” (i.e. last) velocity field.
Summarising, SPO fabric modelling requires one or more of the following: F and/or C obtained from DREX_M; P, T, and/or V and/or fluid and melt fraction fields for the “presentday” state of the geodynamic model. Consequently, this modelling is performed when postprocessing the DREX_M output with the software VIZTOMO (see Sect. 2.2.1).
The total or deviatoric component of the viscous tensor η due to SPO fabrics is estimated for twophase systems with ellipsoidal weak or hard inclusions using the DEM theory, as well as the parameterisation of the viscous tensor evolution and orientation as a function of the cumulated deformation (F) obtained following de Montserrat et al. (2021). Indeed, most, if not all, mantle levels are composed of two main mineral phases that control both the elastic and viscous properties. The case of a multicomponent layered medium is not considered because its viscous tensor can either be approximated with flat inclusions or more simply computed using the Voigt and Reuss averages of the layers' isotropic viscosity.
The first modelling phase requires running the subprogram DEMviscous to generate a database of viscous tensors for a range of inclusion shapes and volume fractions, as well as inclusionmatrix viscosity contrasts. The latter implies that the viscous moduli are dimensionless and can therefore be interpreted as scaling factors with respect to an isotropic effective viscosity of the bulk rock or most abundant mineral phase. Subsequently, the database can be exploited by largescale geodynamic simulations to either (i) return the viscosity tensor η from a lookup table for every point of the computational domain, which can be superimposed on the isotropic effective viscosity computed from flow laws (coupled mechanical simulations), or (ii) estimate η (uncoupled mechanical simulation) and/or visualise its anisotropic viscous properties for the “presentday” state of the geodynamic model with the software VIZVISC (Sect. 2.2.2).
2.2 Visualisation of the mechanical properties and data formatting in preparation for seismological synthetics
2.2.1 VIZTOMO
VIZTOMO processes the DREX_M output for the visualisation of the aggregates' elastic and deformational history properties. Estimation of extrinsic anisotropy effects via the EXEV routines is possible at this stage. The properties of the C and F tensors can either be determined at the position of the Lagrangian mantle aggregates or, in the case of the elastic properties and density, interpolated to a structured (tomographic) grid. In the latter case the grid can be saved in a format suitable for generating synthetic seismic datasets via the PSI_D package (see Sect. 2.3.2) or for 3D waveform simulations in SPECFEM (Komatitsch and Tromp, 1999).
Several properties of the elastic tensor C can be visualised:

isotropic or raypathdependent velocity anomalies and anisotropic elastic properties of body waves (i.e. Pwave anisotropy and direction of maximum Pwave velocity, direction and magnitude of maximum Swave splitting delay time, and Swave radial and azimuthal anisotropy);

reflection and/or transmission energies resulting from the whole range of P–S conversions occurring at discontinuities (useful for studies based on receiver function analysis);

the fraction of the elastic tensor anisotropic component relative to the total and relative contributions of five different anisotropy classes (hexagonal, orthorhombic, tetragonal, monoclinic, triclinic) obtained through elastic tensor decomposition (Browaeys and Chevrot, 2004);

the orientation of the hexagonal symmetry axis (already present in the original DREX).
The deformation history stored in F can be visualised in terms of the FSE shape and orientation and/or length or orientation of its minimum and maximum semiaxes.
The different fields are saved in specific file formats which can be imported by the opensource ParaView software (Ahrens et al., 2005) for visualisation. Figs. 6–9 display some of these fields computed for different geodynamic models.
2.2.2 VIZVISC
VIZVISC processes the DREX_M output for the visualisation of the aggregate properties such as the viscous anisotropy and related deformational history (in terms of the FSE) in ParaView. The deviation from isotropy is evaluated by computing the radial and azimuthal components of viscous anisotropy in a similar way as for the elastic tensor. More in detail, radial viscous anisotropy is defined as $\mathit{\xi}=N/L$, while azimuthal viscous anisotropy is defined by the magnitude $G=\sqrt{{G}_{\mathrm{c}}^{\mathrm{2}}+{G}_{\mathrm{s}}^{\mathrm{2}}}$ and azimuth $\mathit{\varphi}={\mathrm{tan}}^{\mathrm{1}}\left({G}_{\mathrm{s}}{G}_{\mathrm{c}}\right)$, where $N=\mathrm{1}/\mathrm{8}({\mathit{\eta}}_{\mathrm{11}}+{\mathit{\eta}}_{\mathrm{22}})\mathrm{1}/\mathrm{4}{\mathit{\eta}}_{\mathrm{12}}+\mathrm{1}/\mathrm{2}{\mathit{\eta}}_{\mathrm{66}}$, $L=\mathrm{1}/\mathrm{2}({\mathit{\eta}}_{\mathrm{44}}+{\mathit{\eta}}_{\mathrm{55}})$, ${G}_{\mathrm{c}}=\mathrm{1}/\mathrm{2}({\mathit{\eta}}_{\mathrm{55}}{\mathit{\eta}}_{\mathrm{44}})$, and G_{s}=η_{45}. Radial anisotropy and azimuthal viscous anisotropy are evaluated in the FSE (and thus inclusions) reference frame, whereby the minor semiaxis is oriented along the vertical direction, and the intermediate and major semiaxes are in the horizontal plane. As such, radial anisotropy is always ≥1. Figure 10 displays some of these fields computed for a 2D steadystate model of mantle convection with periodic upwellings and downwellings.
2.3 Seismic forward and inverse modelling
2.3.1 SKSSPLIT
SKSSPLIT estimates the SKS splitting at a grid of virtual seismic stations placed at the top of the DREX_M model as a function of the back azimuth using the Fortran routines included in FSTRACK (SchultePelkum and Blackman, 2003; Becker et al., 2006). These routines use the reflectivity method (for details, see Booth and Crampin, 1983; Chapman and Shearer, 1989) to compute synthetic seismograms through layered anisotropic media traversed by an incident plane wave (5° for typical SKS arrivals) over a range of frequencies (0–25 Hz). Bandpass filters from 0.1 to 0.3 Hz are then applied to construct synthetic seismograms in the SKS band (3.3 to 10 s). Successively, the splitting is determined with the crosscorrelation method of Menke and Levin (2003). The routines have been adapted to load the DREX_M output, stack the elastic tensors and densities in an uppermantle rock column beneath each virtual seismic station (Faccenda and Capitanio, 2013), and run in parallel using MPI. The averaged fast azimuth scaled by the delay time can then be visualised in ParaView as shown in Figs. 7 and 8a.
2.3.2 PSI
PSI (Platform for Seismic Imaging) is a Julia package for performing tomographic inversion of both real and synthetic seismic datasets. In the context of the ECOMAN software, it is a useful tool for exploring how features in highresolution geodynamic simulations are mapped into lowerresolution seismic tomography models; this is a critical step if one aims to evaluate geodynamic results against existing tomographic images.
PSI can be used to forwardmodel P and S travel times, splitting intensity (Chevrot, 2000), and shear wave splitting parameters (i.e. the fastpolarisation azimuth and the delay time between the fast and slowpolarised waveforms). Forward modelling of these seismic observables is supported for three different model parameterisations: (1) isotropic P and Swave speeds, (2) hexagonally anisotropic media defined by the five Thomsen parameters and the azimuth and elevation of the symmetry axis, and (3) fully anisotropic models defined by the densitynormalised 21component elastic tensor which can be generated from DREX_M + VIZTOMO (Sect. 2.1.2 and 2.2.1). Seismic phase velocities are computed following Thomsen (1986) for hexagonally anisotropic models, while the Christoffel equations are solved when the elastic tensor parameterisation is used. Anisotropic travel times and splitting intensities for arbitrarily polarised Swaves are computed following VanderBeek et al. (2023) using a longwavelength approximation in which the accumulated delay time between fast and slowpolarised qS waves is less than their dominant period. Lastly, splitting parameters are predicted using the matrix propagation method of Rumpker and Silver (1998). In this initial release of PSI, all predictions are made via integration along ray paths traced through a userdefined 1D reference velocity model using the TauP Toolkit (Crotwell et al., 1999). The 3D model properties are subsequently interpolated to these paths before computing the seismic observables. We anticipate releasing an update to the package that includes both 3D anisotropic ray tracing and finitefrequency kernels that are currently under development.
Travel time and splitting intensity datasets can be inverted individually and jointly either for isotropic (V_{p} and V_{s}) or hexagonally anisotropic model parameters. Hexagonal anisotropy is defined by up to five free parameters, the isotropic (1) P and (2) S velocity, (3) anisotropic magnitude, and the (4) azimuth and (5) elevation of the hexagonal symmetry axis. Only a single anisotropic magnitude parameter is required because the strength of P and S anisotropy is strongly correlated (e.g. Becker et al., 2006). Consequently, the ratios between the three Thomsen parameters (ε, δ, γ; Thomsen, 1986) and the inverted anisotropic magnitude parameter must be chosen a priori and can be spatially variable. Source and receiver statics for each observation and seismic phase type may also be included as free parameters. The tomographic model is obtained by iteratively solving a system of linearised equations relating the perturbations in the inversion parameters to the data residuals augmented with damping and smoothing constraints. The solution is obtained via the LSQR algorithm (Paige and Saunders, 1982). Full details on the tomographic method can be found in VanderBeek and Faccenda (2021) and VanderBeek et al. (2023). The final tomographic solution is written to a VTK file to be visualised in ParaView. Tomographic inversions can be run from a workstation. For the problem shown in Fig. 11 consisting of 12 320 observations and 415 044 free parameters, solutions can be obtained within ∼10 min (less for isotropic inversions) using six cores and 1.8 GB peak RAM.
The PSI inversion methodology has been tested on synthetic models of oceanic plate subduction, intraoceanic upwelling plume, spreading oceanic ridge, and central–eastern Mediterranean subduction (VanderBeek and Faccenda, 2021; Lo Bue et al., 2021, 2022; VanderBeek et al., 2023; Faccenda and VanderBeek, 2023) and applied to the isotropic and anisotropic imaging of the central Mediterranean (Rappisi et al., 2022) and of the Mt. Etna volcanic field (Lo Bue et al., 2024). In Fig. 11, we illustrate results obtained from a synthetic inversion of direct teleseismic P and Swave relative travel times computed through the subduction zone shown in Fig. 7. Synthetic data were computed using the full 21component elastic tensor, while the inversion was performed for the bestfit hexagonal anisotropic parameters. We considered an array of 770 receivers extending from ±7.5° in longitude and ±11.5° in latitude (∼75 km spacing) that recorded 16 events: 8 at 50° and another 8 at 80° from the origin of the subduction zone model and equally distributed in the back azimuth. This example is included in the PSI package.
3.1 Software advantages
When compared to other similar software, ECOMAN (i) is a more versatile package suitable for any geodynamic simulations (2D and 3D; Cartesian and polar coordinate systems; regional and global settings); (ii) allows for the possibility to account for the timedependent deformational history of the mantle (which is usually not steadystate, especially close to plate boundaries and hotspot settings); (iii) computes the straininduced fabric and elastic tensor of different mantle layers (i.e. not only the upper mantle); (iv) includes effective medium theories (STILWE, DEM) and a parameterisation of the fabric evolution of twophase composites to predict elastic and viscous extrinsic anisotropy; (v) generates realistic grid structure distributions of mantle elastic properties to be used for forward and inverse seismological modelling (e.g. PSI, SPECFEM3D); and (vi) performs synthetic seismic inversions (e.g. P and Swave travel time tomographies, Swave splitting intensities) within the computational domain, which facilitates the comparison with other tomographic models and the estimation of apparent anomalies (artefacts) due to, for example, unaccounted for elastic anisotropy (Bezada et al., 2016; VanderBeek and Faccenda, 2021; VanderBeek et al., 2023) and/or regularisation.
3.2 Software (real or potential) current limitations
DREX_S and DREX_M compute the straininduced LPO for the most abundant and highly anisotropic mineral phases, while assuming random orientation for several secondary phases that instead could contribute substantially to the aggregate anisotropic properties. For example, cubic ferropericlase is relatively abundant (<20 %) and becomes highly anisotropic in the lower half of the lower mantle (Avs = 30 %–50 %) such that it could dominate seismic anisotropy at these depths (Marquardt et al., 2009). Davemaoite (Caperovskite) is also a highly anisotropic mineral with cubic symmetry (Avs = 25 %–15 %; Kawai and Tsuchiya, 2015) but its model abundance is quite low (<10 %). Recently, micromechanical simulations of straininduced LPO in aggregates with a pyrolite mantle composition have shown that the bulk aggregate seismic anisotropy is controlled by bridgmanite and postperovskite, while the cubic secondary phases appear to only slightly reduce the amplitude of anisotropy (Chandler et al., 2021). The latter effect can thus be well approximated by randomising the orientation of the cubic phases as assumed here.
A main limitation of ECOMAN is that it does not yet include the modelling of viscous anisotropy due to the intrinsic mechanical anisotropy of crystals. Tommasi et al. (2009) used VPSC simulations to show that olivine CPO in the lithosphere can control the reactivation of fossil faults misoriented with respect to the stress field. Király et al. (2020) employed the modified director method (MDM) and estimated up to 1 order of magnitude of intrinsic viscous anisotropy in olivine aggregates, which can control the kinematics and dynamics of tectonic plates. Both the VPSC and MDM approaches are computationally expensive and appear to be prohibitive for the number of mantle aggregates required to discretise the mantle domain of largescale 3D simulations. An alternative approach which minimises the computational time is therefore desired.
Fraters and Billen (2021) released a version of the DREX model that is embedded into the ASPECT geodynamic modelling software (Bangerth et al., 2020) and that takes into account changes in olivine fabric as a function of the water content and deviatoric stress (Karato et al., 2008). Although at present DREX_M only allows the definition of a single olivine fabric per run, changes in crystal fabrics can be easily implemented in the future by loading the water content and deviatoric stress field from the geodynamic models and by parameterising the boundaries of the olivine fabric domains.
Tape and Tape (2024) recently reformulated the elastic decomposition method proposed by Browaeys and Chevrot (2004) that they argued is not entirely accurate. Although this will only affect the visualisation of the tensor components as performed by DREX_S and VIZTOMO, it is planned to replace the existing elastic tensor decomposition methodology with the one proposed by Tape and Tape (2024).
A minor limitation is that the different code modules are based on different programming languages (Fortran, Julia), libraries (OpenMP, MPI, HDF5), and software (MATLAB, ParaView), whose installation on local devices might discourage potential users. There are a few main reasons for this. First of all, the original DREX software was written in Fortran, and thus its modification in DREX_S and DREX_M software was more straightforward by maintaining the same programming language. Fortran has high performance standards on HPC clusters that are often if not always higher than other interpreted programming languages (e.g. MATLAB or Python). Given the largescale computational power needed for the DREX_M simulations, especially those with hundreds of thousands or millions of crystal aggregates, and considering that the required compilers and libraries are routinely installed in most (if not all) HPC clusters, the usage of the ECOMAN's Fortranbased applications in highperformance environments is warranted. In contrast, the visualisation of the Fortranbased applications' output in MATLAB and ParaView can be performed on local devices. In particular, the MTEX software is a MATLAB toolbox for visualisation of the LPO and elastic tensor pole figures which should be downloaded (https://mtextoolbox.github.io/download, last access: 1 October 2024) and installed along with MATLAB only when using D_REX_S. Seismic forward and inverse modelling with PSI can be performed on either an HPC cluster or local devices upon installation of the Julia package. Julia was chosen because it is an opensource and highlevel language that offers a number of performance benefits over other popular scientific computing languages such as Python and MATLAB. It is important to stress that for any of these applications the user only needs to modify input text files, and thus no particular prerequisite or computational skill is required.
Finally, running DREX_M requires allocating ∼160 GB of memory per million crystal aggregates with 1000×2 crystals each. This potential problem can be addressed by distributing the computational and memory load over several nodes, which is possible thanks to the hybrid MPI and OpenMP parallelisation scheme.
ECOMAN is an opensource software package for estimating strain and stressinduced fabrics in mantle aggregates, their mechanical properties, and how mechanical anisotropy affects the geodynamic evolution and seismic imaging of the Earth's interior. Programs included in ECOMAN are portable across different HPC and local device systems (provided the Julia package and Fortran compilers are available) and are applicable to any 2D–3D geodynamic simulation. Computationally expensive programs such as DREX_M are parallelised, offering a nearly perfect scaling with an increasing number of cores. As a result, the straininduced fabrics of millions of mantle aggregates can now be estimated with a reasonable quantity of time and computational resources.
As ongoing developments, we are seeking to include in ECOMAN micromechanical modelling methods that are capable of estimating the straininduced LPO and/or the intrinsic viscous anisotropy at computational speeds that are orders of magnitude faster than current ones. For instance, Ribe et al. (2019) proposed an analytical finitestrain parameterisation for texture evolution in deforming olivine polycrystals that is ≈10^{7} times faster than full homogenisation approaches such as the secondorder selfconsistent model. When implemented in ECOMAN, preliminary tests indicate that this new method outperforms DREX by 1–2 orders of magnitude (Ribe et al., 2023).
In addition, in the near future the PSI software will be updated to include transdimensional Bayesian Monte Carlo sampling methods that, in contrast to deterministic approaches, address the consequences of underdetermination in seismic imaging by constraining the uncertainty (Del Piccolo et al., 2024). Isotropic and anisotropic seismic imaging with PSI is currently feasible using body wave information such as travel times and Swave splitting intensity. However, when local deep seismicity is absent, as is the case in warm subduction zones or at spreading ridges and intraplate settings, the retrieved isotropic and anisotropic mantle structures are only partially recovered and often affected by smearing (e.g. Faccenda and VanderBeek, 2023). Consequently, we are planning to complement body wave information with surface wave data to improve the seismic ray coverage and the resolving power of tomographic models. Lastly, to improve the prediction of seismic observables, 3D anisotropic ray tracing and finitefrequency kernels are planned for a future release.
We compare the numerical solution of the DREX model included in DREX_S and DREX_M with the analytical solution provided by Fraters and Billen (2021, Supplement) for a test with an uppermantle aggregate ($Ol:Ens=\mathrm{50}:\mathrm{50})$ composed of five crystals for each phase, ${M}^{*}=\mathrm{50}$, χ=0.3, and $\mathrm{nCRSS}=(\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{\infty})$ for the four olivine slip systems [100](010), [100](001), [001](010), and [001](100) (Table A1). For the first crystal, the time derivative of the cosine direction matrix (CDM) is $\frac{\mathit{\delta}{a}_{ij}^{\mathrm{v}}}{\mathit{\delta}t}=\mathrm{0}$ and the vorticity matrix is ${\stackrel{\mathrm{\u203e}}{\mathit{\omega}}}_{ij}^{\mathrm{v}}={\mathit{\epsilon}}_{jki}{\mathit{\omega}}_{k}^{v}=\mathrm{0}$, while for the other four crystals
with a maximum error of $\mathrm{3.1}\times {\mathrm{10}}^{\mathrm{3}}$ %, which we attribute to truncation errors, and
The changes in volume fractions are $\frac{\mathit{\delta}{f}_{\mathrm{m}}^{\mathrm{v}}}{\mathit{\delta}t}=\mathrm{0.630116863029}$ for grain 1 and $\frac{\mathit{\delta}{f}_{\mathrm{m}}^{\mathrm{v}}}{\mathit{\delta}t}=\mathrm{0.157529215757}$ for the other four crystals, resulting in a misfit of $\mathrm{6.5244}\times {\mathrm{10}}^{\mathrm{4}}$ % relative to the analytical solution.
We provide the DREX_S_test program for reproducing this simple test in the Supplement.
At any P–T condition, the singlecrystal elastic moduli ${C}_{ij}^{\mathrm{m}}$ are (in Voigt notation)
where the partial derivatives are listed in Table B2 and ΔP and ΔT are deviations from the room conditions (${P}_{\mathrm{0}}={\mathrm{10}}^{\mathrm{4}}$ GPa, T_{0}=298 K; Mainprice, 2007). The compliance tensor is the inverse of the stiffness tensor, i.e. ${\mathbf{S}}^{\mathrm{m}}={\mathbf{C}}^{{\mathrm{m}}^{\mathrm{1}}}$. These tensors can be rotated as
where a is the CDM defining the crystal orientation with respect to the external reference frame. The stiffness and compliance tensors of the crystal aggregate (i.e. rock) are then obtained by arithmetic averaging:
where M=2 is the number of minerals in the twophase mantle aggregates modelled in DREX_M and DREX_S, X_{m} their relative abundance, N the number of crystals of each phase, and X_{v} the normalised crystal volume fraction ($\sum _{v=\mathrm{1}}^{N}{X}_{\mathrm{v}}=\mathrm{1}$). C^{r} and (S^{r})^{−1} correspond to the Voigt and Reuss average stiffness tensor of the rock, i.e. ${\mathbf{C}}_{\mathrm{Voigt}}^{\mathrm{r}}$ and ${\mathbf{C}}_{\mathrm{Reuss}}^{\mathrm{r}}$. More in general, the average stiffness tensor resulting from the linear combination of the two end members is found as
where X_{Voigt} is a free parameter defining the fraction of ${\mathbf{C}}_{\mathrm{Voigt}}^{\mathrm{r}}$ relative to ${\mathbf{C}}_{\mathrm{Reuss}}^{\mathrm{r}}$. It is clear that the Hill average stiffness tensor ${\mathbf{C}}_{\mathrm{Hill}}^{\mathrm{r}}$ is obtained when X_{Voigt}=0.5.
The isotropic component of the average stiffness tensor ${\stackrel{\mathrm{\u203e}}{\mathbf{C}}}_{\mathrm{ISO}}^{\mathrm{r}}$ can be assembled as ${\stackrel{\mathrm{\u203e}}{C}}_{\mathrm{ISO}\mathrm{\_}\mathrm{11}}^{\mathrm{r}}={\stackrel{\mathrm{\u203e}}{C}}_{\mathrm{ISO}\mathrm{\_}\mathrm{22}}^{\mathrm{r}}={\stackrel{\mathrm{\u203e}}{C}}_{\mathrm{ISO}\mathrm{\_}\mathrm{33}}^{\mathrm{r}}=K+\frac{\mathrm{4}}{\mathrm{3}}G$, ${\stackrel{\mathrm{\u203e}}{C}}_{\mathrm{ISO}\mathrm{\_}\mathrm{12}}^{\mathrm{r}}={\stackrel{\mathrm{\u203e}}{C}}_{\mathrm{ISO}\mathrm{\_}\mathrm{13}}^{\mathrm{r}}={\stackrel{\mathrm{\u203e}}{C}}_{\mathrm{ISO}\mathrm{\_}\mathrm{23}}^{\mathrm{r}}=K\frac{\mathrm{2}}{\mathrm{3}}G$, and ${\stackrel{\mathrm{\u203e}}{C}}_{\mathrm{ISO}\mathrm{\_}\mathrm{44}}^{\mathrm{r}}={\stackrel{\mathrm{\u203e}}{C}}_{\mathrm{ISO}\mathrm{\_}\mathrm{55}}^{\mathrm{r}}={\stackrel{\mathrm{\u203e}}{C}}_{\mathrm{ISO}\mathrm{\_}\mathrm{66}}^{\mathrm{r}}=G$, knowing that:=
The isotropic component of the aggregate stiffness tensor ${\stackrel{\mathrm{\u203e}}{\mathbf{C}}}_{\mathrm{ISO}}^{\mathrm{r}}$ is then replaced with the isotropic stiffness tensor ${\mathbf{C}}_{\mathrm{ISO}}^{\mathrm{MMA}\phantom{\rule{0.25em}{0ex}}\mathrm{EoS}}$ obtained using the bulk and shear moduli from the lookup tables computed with the MMA_EoS software (Chust et al., 2017) as
^{*} (1) Kelemen and Ghiorso (1986); (2) Chemia et al. (2015); (3) Workman and Hart (2005); (4) Hirschmann et al. (2003).
^{*} (1) Isaak (1992); (2) Zha et al. (1998); (3) Chai et al. (1992); (4) Jackson et al. (2007); (5) Zha et al. (1997); (6) Zhou et al. (2022); (7) Sinogeikin et al. (2003); (8) Chai and Brown (1997); (9) Sinogeikin and Bass (2002); (10) Wentzcovitch et al. (2004). Derivatives at 100 GPa and 300 K: (11) Zhang et al. (2013). Derivatives extrapolated by interpolation: (12) Karki et al. (2000); (13) Zhang et al. (2016). Derivatives extrapolated by interpolation.
In polar grids the velocity field is typically defined by the longitudinal, radial, and (in 3D) colatitudinal components ${\mathit{V}}_{(\mathit{\varphi},\mathrm{r},\mathit{\theta})}=\left({V}_{\mathit{\varphi}}{V}_{\mathrm{r}}{V}_{\mathit{\theta}}\right)$. However, the crystals' orientation and their rotation are defined relative to the external Cartesian reference frame. As such, the velocity gradient tensor must be computed in a Cartesian coordinate system. This is done by rotating the velocity field from spherical to Cartesian coordinates as
Using the finitedifference approach, the spatial gradients of the Cartesian velocity components are computed on the spherical grid as
Defining J as the Jacobian matrix,
the velocity gradient tensor in Cartesian coordinates is ${L}_{\frac{(x,y,z)}{(x,y,z)}}={L}_{\frac{(x,y,z)}{(\mathit{\varphi},r,\mathit{\theta})}}\mathbf{J}$:
When multiple creep mechanisms are active, we assume that only the fraction of deformation accommodated by dislocation creep contributes to the LPO development. This fraction of the total deformation is defined as ${F}_{\mathrm{d}}={\mathit{\eta}}_{\mathrm{eff}}/{\mathit{\eta}}_{\mathrm{disl}}$, where $\mathrm{0}\le {F}_{\mathrm{d}}\le \mathrm{1}$ is interpolated from the largescale geodynamic model, η_{eff} is the effective viscosity resulting from the harmonic average of the viscosities from the flow law of different creep mechanisms, and η_{disl} is the viscosity from the dislocation creep flow law. At a given time step Δt, crystal rotation according to the DREX model is applied for the fraction $\mathrm{\Delta}{t}_{\mathrm{disl}}=\mathrm{\Delta}t\cdot {F}_{\mathrm{d}}$, while for the remaining time $\mathrm{\Delta}{t}_{\mathrm{diff}}=\mathrm{\Delta}t\cdot (\mathrm{1}{F}_{\mathrm{d}})$ we follow the numerical study of Hedjazian et al. (2017) and apply fluid body rotation by multiplying the CDM of the aggregates by the rotation of the FSE. In this way, we preserve and do not alter the strength of the LPO.
Given the orthogonal CDMs defining the original and final crystal orientation a^{t} and ${\mathbf{a}}^{t+\mathrm{\Delta}{t}_{\mathrm{diff}}}$ and those defining the original and final FSE semiaxis orientation v^{t} and ${\mathit{v}}^{t+\mathrm{\Delta}{t}_{\mathrm{diff}}}$ (i.e. the eigenvectors of, respectively, the left stretch tensor LS^{t}=F^{t}(F^{t})^{T} and ${\mathbf{LS}}^{t+\mathrm{\Delta}{t}_{\mathrm{diff}}}={\mathbf{F}}^{t+\mathrm{\Delta}{t}_{\mathrm{diff}}}({\mathbf{F}}^{t+\mathrm{\Delta}{t}_{\mathrm{diff}}}{)}^{T}$), the following relation implies that after fluid body rotation the orientation of the crystal relative to FSE semiaxes must remain constant:
The new crystal orientation is then
The version of ECOMAN described in this article is 2.0 and is freely available via GitHub (https://github.com/ecomangeos, ECOMAN project, 2024) and Zenodo (ECOMAN2.0geodynamics: https://doi.org/10.5281/zenodo.10599735, Faccenda, 2024a; ECOMAN2.0seismology.SKSsplit: https://doi.org/10.5281/zenodo.11173260, Faccenda, 2024b; ECOMAN2.0seismology.PSI_D: https://doi.org/10.5281/zenodo.11186805, VanderBeek,, 2024).
The software package contains the input files to generate the synthetic models and datasets discussed in this manuscript. No data have been produced for this work.
The supplement related to this article is available online at: https://doi.org/10.5194/se1512412024supplement.
MF and BPV developed the software package, with important contributions from AdM, JY, and FR. MF wrote the manuscript draft. MF acquired the funding to support the research activities. All authors have contributed to the discussion and manuscript editing.
The contact author has declared that none of the authors has any competing interests.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.
Manuele Faccenda is indebted to Fabio Antonio Capitanio for supporting the beginning of this journey, Thorsten Becker for providing routines to compute shear wave splitting, Edouard Kaminski for fruitful discussions on straininduced LPO modelling, and Ana Ferreira for the continuous and stimulating feedback. Two anonymous reviewers helped improve the article.
This research has been supported by the European Research Council, H2020 European Research Council (grant no. 758199).
This paper was edited by Juliane Dannberg and reviewed by two anonymous referees.
Abt, D. L. and Fischer, K.: Resolving threedimensional anisotropic structure with shear wave splitting tomography, Geophys. J. Int., 173, 859–886, 2008.
Ahrens, J., Geveci, B., and Law, C.: ParaView: An EndUser Tool for Large Data Visualization, in: Visualization Handbook, edited by: Hansen, C. D. and Johnson, C. R., Elsevier, 717–731, https://doi.org/10.1016/B9780123875822/500381, 2005.
Almqvist, B. S. G. and Mainprice, D.: Seismic properties and anisotropy of the continental crust: Predictions based on mineral texture and rock microstructure, Rev. Geophys., 55, 367–433, https://doi.org/10.1002/2016RG000552, 2017.
Backus, G. E.: Long‐wave elastic anisotropy produced by horizontal layering. J. Geophys. Res., 67, 4427–4440, https://doi.org/10.1029/JZ067i011p04427, 1962.
Bangerth, W., Dannberg, J., Gassmoeller, R., and Heister, T.: Aspect v2.2.0, Zenodo [code], https://doi.org/10.5281/zenodo.3924604, 2020.
Becker, T. W., Chevrot, S., SchultePelkum, V., and Blackman, D. K.: Statistical properties of seismic anisotropy predicted by upper mantle geodynamic models, J. Geophys. Res., 111, B08309, https://doi.org/10.1029/2005JB004095, 2006.
Bezada, M., Faccenda, M., and Toomey, D. R: Representing anisotropic subduction zones with isotropic velocity models: A characterization of the problem and some steps on a possible path forward, Geochem. Geophy. Geosy., 17, 3164–3189, https://doi.org/10.1002/2016GC006507, 2016.
Booth, D. C. and Crampin, S.: The anisotropic reflectivity technique: Theory, Geophys. J. Roy. Astron. Soc., 72, 31–45, 1983.
Browaeys, J. T. and Chevrot, S.: Decomposition of the elastic tensor and geophysical applications, Geophys. J. Int., 159, 667–678, 2004.
Bunge, H.P., Hagelberg, C. R., and Travis, B.J.: Mantle circulation models with variational data assimilation: inferring past mantle flow and structure from plate motion histories and seismic tomography, Geophys. J. Int., 152, 280–301, 2003.
Chai, M. and Brown, J.: The elastic constants of a pyropegrossularalmandine garnet to 20 GPa, Geophys. Res. Lett., 24, 523–526, 1997.
Chai, M., Brown, J., and Slutski, L.: The elastic constants of an aluminous orthopyroxene to 12.5 GPa, J. Geophys. Res., 102, 779–785, 1992.
Chandler, B. C., Chen, L.W., Li, M., Romanowicz, B., and Wenk, H.R.: Seismic anisotropy, dominant slip systems and phase transitions in the lowermost mantle, Geophys. J. Int., 227, 1665–1681, https://doi.org/10.1093/gji/ggab278, 2021.
Chang, S.J., Ferreira, A. M. G., Ritsema, J., van Heijst, H. J., and Woodhouse, J. H.: Joint inversion for global isotropic and radially anisotropic mantle structure including crustal thickness perturbations, J. Geophys. Res.Solid, 120, 4278–4300, https://doi.org/10.1002/2014JB011824, 2015.
Chang, S.J., Ferreira, A. M. G., and Faccenda, M.: Upper and midmantle interaction between the Samoan plume and the TongaKermadec slabs, Nat. Commun., 7, 10799, https://doi.org/10.1038/ncomms10799, 2016.
Chapman, C. H. and Shearer, P. M.: Ray tracing in azimuthally anisotropic media: II. Quasishear wave coupling, Geophys. J. Int., 96, 65–83, 1989.
Chemia, Z., Dolejš, D., and SteinleNeumann, G.: Thermal effects of variable material properties and metamorphic reactions in a threecomponent subducting slab, J. Geophys. Res., 120, 6823–6845, 2015.
Chevrot, S.: Multichannel analysis of shear wave splitting, J. Geophys. Res., 105, 21579–21590, 2000.
Christensen, U. C.: Some geodynamical effects of anisotropic viscosity, Geophys. J. Roy. Astron. Soc., 91, 711–736, 1987.
Chust, T. C., SteinleNeumann, G., Dolejs, D., Schuberth, B. S., and Bunge, H. P.: MMA EoS: a computational framework for mineralogical thermodynamics, J. Geophys. Res., 122, 9881–9920, 2017.
Crameri, F. and Tackley, P. J.: Spontaneous development of arcuate singlesided subduction in global 3Dmantle convection models with a free surface, J. Geophys. Res.Solid, 119, 921–5942, https://doi.org/10.1002/2014JB010939, 2014.
Crotwell, H. P., Owens, T. J., and Ritsema, J.: The TauP Toolkit: flexible seismic traveltime and raypath utilities, Seismol. Res. Lett., 70, 154–160, 1999.
Davies, D. R., Goes, S., Davies, J. H., Schuberth, B. S. A., Bunge, H.P., and Ritsema, J.: Reconciling dynamic and seismic models of Earth's lower mantle: The dominant role of thermal heterogeneity, Earth Planet. Sc. Lett., 353–354, 253–269, 2012.
Debayle, E., Kennet, B., and Priiestley, K.: Global azimuthal seismic anisotropy and the unique platemotion deformation of Australia, Nature, 433, 509–512, 2005.
Debayle, E., Bodin, T., Durand, S., and Ricard, Y.: Seismic evidence for partial melt below tectonic plates, Nature, 586, 555–559, https://doi.org/10.1038/s4158602028094, 2020.
Del Piccolo, G., VanderBeek, B., Faccenda, M., Morelli, A., and Byrnes, J.: Imaging upper mantle anisotropy with transdimensional Bayesian Monte Carlo sampling, Bull. Seismol. Soc. Am., 114, 1214–1226, https://doi.org/10.1785/0120230233, 2024.
de Montserrat, A., Faccenda, M., and Pennacchioni, G.: Extrinsic anisotropy of twophase Newtonian aggregates: Fabric characterization and parameterization, J. Geophys. Res.Solid, 126, e2021JB022232, https://doi.org/10.1029/2021JB022232, 2021.
Deuss, A: Heterogeneity and anisotropy of Earth's inner core, Annu. Rev. Earth. Planet. Sci., 42, 103–126, https://doi.org/10.1146/annurevearth060313054658, 2014.
ECOMAN project: ecomangeos, GitHub, https://github.com/ecomangeos (last access: 1 October 2024), 2024.
Faccenda, M.: Mid mantle seismic anisotropy around subduction zones, Phys. Earth Planet. Int., 227, 1–19, 2014.
Faccenda, M.: ECOMAN2.0geodynamics, Zenodo [code], https://doi.org/10.5281/zenodo.10599735, 2024a.
Faccenda, M.: ECOMAN2.0seismology. SKSsplit, Zenodo [code], https://doi.org/10.5281/zenodo.11173260, 2024b.
Faccenda, M. and Capitanio, F. A.: Development of mantle seismic anisotropy during subductioninduced 3D flow, Geophys. Res. Lett., 39, L11305, https://doi.org/10.1029/2012GL051988, 2012.
Faccenda, M. and Capitanio, F. A.: Seismic anisotropy around subduction zones: insights from threedimensional modeling of upper mantle deformation and SKS splitting calculations, Geochem. Geophy. Geosy., 14, 243–262, https://doi.org/10.1002/ggge.20055, 2013.
Faccenda, M. and VanDerBeek, B. P.: On constraining 3D seismic anisotropy in subduction, midoceanridge, and plume environments with teleseismic body wave data, J. Geodyn., 158, 102003, https://doi.org/10.1016/j.jog.2023.102003, 2023.
Faccenda, M., Ferreira, A. M. G., Tisato, N., LithgowBertelloni, C., Stixrude, L., and Pennacchioni, G.: Extrinsic anisotropy in a compositionally heterogeneous mantle, J. Geophys. Res., 124, 1671–1687, https://doi.org/10.1029/2018JB016482, 2019.
Ferreira, A. M. G., Faccenda, M., Sturgeon, W., Chang, S.J., and Schardong, L.: Ubiquitous lowermantle anisotropy beneath subduction zone, Nat. Geosci., 12, 301–306, https://doi.org/10.1038/s4156101903257, 2019.
Fraters, M. R. T. and Billen, M. I.: On the implementation and usability of crystal preferred orientation evolution in geodynamic modelling, Geochem. Geophy., Geosy., 22, e2021GC009846, https://doi.org/10.1029/2021GC009846, 2021.
French, S. and Romanowicz, B.: Broad plumes rooted at the base of the Earth's mantle beneath major hotspots, Nature, 525, 95–99, https://doi.org/10.1038/nature14876, 2015.
Halter, W., Macherel, E., Duretz, T., and Schmalholz, S. M.: Numerical modelling of strain localization by anisotropy evolution during 2D viscous simple shearing, EGU General Assembly 2022, Vienna, Austria, 23–27 May 2022, EGU2211438, https://doi.org/10.5194/egusphereegu2211438, 2022.
Han, D. and Wahr, J.: An analysis of anisotropic mantle viscosity, and its possible effects on postglacial rebound, Phys. Earth Planet. Int., 102, 33–50, 1997.
Hansen, L. N., Zhao, Y.H., Zimmerman, M. E., and Kohlstedt, D. L.: Protracted fabric evolution in olivine: Implications for the relationship among strain, crystallographic fabric, and seismic anisotropy, Earth Planet. Sc. Lett., 387, 157–168, https://doi.org/10.1016/j.epsl.2013.11.009, 2014.
Hedjazian, N., Garel, F., Rhodri Davies, D., and Kaminski, E.: Ageindependent seismic anisotropy under oceanic plates explained by strain history in the asthenosphere, Earth Planet. Sc. Lett., 460, 135–142, 2017.
Hirschmann, M., Kogiso, T., Baker, M., and Stopler, E.: Alkalic magmas generated by partial melting of garnet pyroxenite, Geology, 31, 481–484, 2003.
Hu, J., Faccenda, M., and Liu, L.: Subductioncontrolled mantle flow and seismic anisotropy in South America, Earth Planet. Sc. Lett., 470, 13–24, https://doi.org/10.1016/j.epsl.2017.04.027, 2017.
Isaak, D.: Hightemperature elasticity of ironbearing olivines, J. Geophys. Res., 97, 1871–1885, 1992.
Ito, G., Dunn, R., Li, A., Wolfe, C. J., Gallego, A., and Fu, Y.: Seismic anisotropy and shear wave splitting associated with mantle plumeplate interaction, J. Geophys. Res.Solid, 119, 4923–4937, https://doi.org/10.1002/2013JB010735, 2014.
Jackson, J., Sinogeikin, S., and Bass, J.: Sound velocities and singlecrystal elasticity of orthoenstatite to 1073 K at ambient pressure, Phys. Earth Planet. Int., 171, 1–12, 2007.
Jadamec, M. and Billen, M.: Reconciling surface plate motions with rapid threedimensional mantle flow around a slab edge, Nature, 465, 338–341, https://doi.org/10.1038/nature09053, 2010.
Kaminski, E. and Ribe N. M.: A kinematic model for recrystallization and texture development in olivine polycrystals, Earth Planet. Sc. Lett., 189, 253–267, 2001.
Kageyama, A. and Sato, T: “YinYang grid”: An overset grid in spherical geometry, Geochem. Geophy. Geosy., 5, Q09005, https://doi.org/10.1029/2004GC000734, 2004.
Kaminski, E., Ribe, N. M., and Browaeys, J. T.: DRex, a program for calculation of seismic anisotropy due to crystal lattice preferred orientation in the convective upper mantle, Geophys. J. Int., 158, 744–752, https://doi.org/10.1111/j.1365246x.2004.02308.x, 2004.
Karato, S.I.: Seismic anisotropy in the deep mantle, boundary layers and the geometry of mantle convection, Pure Appl. Geophys., 151, 24, 565–587, https://doi.org/10.1007/s000240050130, 1998.
Karato, S.I., Jung, H., Katayama, I., and Skemer, P.: Geodynamic significance of seismic anisotropy of the upper mantle: New insights from laboratory studies, Annu. Rev. Earth Planet. Sci., 36, 59–95, https://doi.org/10.1146/annurev.earth.36.031207.124120, 2008.
Karki, B., Wentzcovitch, R., de Gironcoli, S., and Baroni, S.: Highpressure lattice dynamics and thermoelasticity of MgO, Phys. Rev. B, 61, 8793–8800, 2000.
Kawai, K. and Tsuchiya, T.: Small shear modulus of cubic CaSiO_{3} perovskite, Geophys. Res. Lett., 42, 2718–2726, https://doi.org/10.1002/2015GL063446, 2015.
Kelemen, P. and Ghiorso, M.: Assimilation of peridotite in zoned calcalkaline plutonic complexes: evidence from the Big Jim complex, Washington Cascades, Contrib. Min., 94, 12–28, 1986.
Kendall, J.M.: Seismic anisotropy in the boundary layers of the mantle, in: Earth's deep interior: Mineral physics and tomography from the atomic to the global scale, Geophysical Monograph Series 117, American Geophysical Union, Washington, DC, 133–159, https://doi.org/10.1029/GM117p0133, 2000.
Király, Á., Conrad, C. P., and Hansen, L. N.: Evolving viscous anisotropy in the upper mantle and its geodynamic implications, Geochem. Geophy. Geosy., 21, e2020GC009159, https://doi.org/10.1029/2020GC009159, 2020.
Kocher, T., Schmalholz, S. M., and Mancktelow, N. S.: Impact of mechanical anisotropy and powerlaw rheology on single layer folding, Tectonophysics, 421, 71–87, 2006.
Komatitsch, D. and Tromp, J.: Introduction to the spectral element method for threedimensional seismic wave propagation, Geophys. J. Int., 139, 806–822, 1999.
Lev, E. and Hager, B. H.: Rayleigh–Taylor instabilities with anisotropic lithospheric viscosity, Geophys. J. Int., 173, 806–814, 2008a.
Lev, E. and Hager, B. H.: Prediction of anisotropy from flow models: A comparison of three methods, Geochem. Geophy. Geosy., 9, Q07014, https://doi.org/10.1029/2008GC002032, 2008b.
Li, Z.H., Di Leo, J. F., and Ribe, N. M.: Subductioninduced mantle flow, finite strain, and seismic anisotropy: Numerical modeling, J. Geophys. Res.Solid, 119, 5052–5076, https://doi.org/10.1002/2014JB010996, 2014.
Lo Bue, R., Faccenda, M., and Yang, J.: The role of Adria Plate Lithospheric Structures on the Recent Dynamics of the Central Mediterranean Region, J. Geophys. Res., 126, e2021JB022377, https://doi.org/10.1029/2021JB022377, 2021.
Lo Bue, R., Rappisi, F., VanderBeek, B. P., and Faccenda, M.: Tomographic Image Interpretation and CentralWestern MediterraneanLike Upper Mantle Dynamics From Coupled Seismological and Geodynamic Modeling Approach, Front. Earth Sci., 10, 884100, https://doi.org/10.3389/feart.2022.884100, 2022.
Lo Bue, R., Rappisi, F., Firetto Carlino, M., Giampiccolo, E., Cocina, O., Vanderbeek, B. P., and Faccenda, M.: Crustal structure of Etna volcano (Italy) from Pwave anisotropic tomography, Geophys. Res. Lett., 51, e2024GL108733, https://doi.org/10.1029/2024GL108733, 2024.
Long, M. D. and Becker, T. W.: Mantle dynamics and seismic anisotropy, Earth Planet. Sc. Lett., 297, 341–354, 2010.
Long, M. D., Maarten, V., and Van Der Hilst, R. D.: Waveequation shear wave splitting tomography, Geophys. J. Int., 172, 311–330, https://doi.org/10.1111/j.1365246X.2007.03632.x, 2008.
Maguire, R., Ritsema, J., Bonnin, M., van Keken, P. E., and Goes, S.: Evaluating the resolution of deep mantle plumes in teleseismic traveltime tomography, J. Geophys. Res., 123, 384–400, https://doi.org/10.1002/2017JB014730, 2018.
Mainprice, D.: Modelling the anisotropic seismic properties of partially molten rocks found at midocean ridges, Tectonophysics, 279, 161–179, 1997.
Mainprice, D.: Seismic Anisotropy of the Deep Earth from a Mineral and Rock Physics Perspective, in: Treatise on Geophysics, 2, edited by: Schubert, G., Elsevier, 437–491, ISBN 9780444538031, 2007.
Mainprice, D., Hielscher, R., and Schaeben, H.: Calculating anisotropic physical properties from texture data using the MTEX open source package, in: Deformation Mechanisms, Rheology and Tectonics: Microstructures, Mechanics and Anisotropy, edited by: Prior, D. J., Rutter, E. H., and Tatham, D. J., Geol. Soc. Lond. Spec. Publ., 360, 175–192, 2011.
Marquardt, H., Speziale, S., Reichmann, H. J., Frost, D. J., Schilling, F. R., and Garnero, E. J.: Elastic Shear Anisotropy of Ferropericlase in Earth's Lower Mantle, Science, 324, 224–226, https://doi.org/10.1126/science.1169365, 2009.
Menke, W. and Levin, V.: The crossconvolution method for interpreting SKS splitting observations, with application to one and two layer anisotropic earth models, Geophys. J. Int., 154, 379–392, 2003.
Mondal, P. and Long, M. D.: A model space search approach to finitefrequency SKS splitting intensity tomography in a reduced parameter space, Geophys. J. Int., 217, 238–256, https://doi.org/10.1093/gji/ggz016, 2019.
Mühlhaus, H.B., Moresi, L., Hobbs, B. and Dufour, F.: Large amplitude folding in finely layered viscoelastic rock structures, Pure Appl. Geophys., 159, 2311–2333, 2002.
Mühlhaus, H.B., Moresi, L., and Cada, M.: Emergent anisotropy and flow alignment in viscous rock, Pure Appl. Geophys., 161, 2451–2463, 2004.
Müller, R. D., Flament, N., Cannon, J., Tetley, M. G., Williams, S. E., Cao, X., Bodur, Ö. F., Zahirovic, S., and Merdith, A.: A tectonicrulesbased mantle reference frame since 1 billion years ago – implications for supercontinent cycles and plate–mantle system evolution, Solid Earth, 13, 1127–1159, https://doi.org/10.5194/se1311272022, 2022.
Munzarova, H., Plomerova, J., and Kissling, E.: Novel anisotropic teleseismic body wave tomography code AniTomo to illuminate heterogeneous anisotropic upper mantle: Part I – theory and inversion tuning with realistic synthetic data, Geophys. J. Int., 215, 524–545, 2018.
Oganov, A. R. and Ono, S.: Theoretical and experimental evidence for a postperovskite phase of MgSiO_{3} in Earth's D^{′′} layer, Nature, 430, 445–448, 2004.
Paige, C. C. and Saunders, M. A.: LSQR: An algorithm for sparse linear equations and sparse least squares, ACM Trans. Math. Softw., 8, 43–71, 1982.
Panning, M. and Romanowicz, B.: A threedimensional radially anisotropic model of shear velocity in the whole mantle, Geophys. J. Int., 167, 361–379, 2006.
Rappisi, F., VanderBeek, B., Faccenda, M., Morelli, A., and Molinari, I: Slab geometry and upper mantle flow patterns in the Central Mediterranean from 3D anisotropic Pwave tomography, J. Geophys. Res., 127, e2021JB023488, https://doi.org/10.1029/2021JB023488, 2022.
Rappisi, F., Witek, M., Faccenda, M., Ferreira, A. M. G., and Chang, S.J.: Artificial ageindependent seismic anisotropy, slab thickening and shallowing due to limited resolving power of (an)isotropic tomography, Geophys. J. Int., 237, 217–234, 2024.
Ribe, N. M., Hielscher, R. H., and Castelnau, O.: An analytical finitestrain parameterization for texture evolution in deforming olivine polycrystals, Geophys. J. Int., 216, 486–514, https://doi.org/10.1093/gji/ggy442, 2019.
Ribe, N. M., Faccenda, M., and Hielscher, R. H.: SBFTEX: An Analytical Parameterization for Finite StrainInduced upper mantle Anisotropy, in: AGU Fall Meeting 2023, San Francisco, USA, DI42A08, https://agu.confex.com/agu/fm23/meetingapp.cgi/Paper/1239978 (last access: 1 October 2024), 2023.
Rudolph, M. L., Lekić, V., and LithgowBertelloni, C.: Viscosity jump in Earth's midmantle, Science, 350, 1349–1352, https://doi.org/10.1126/science.aad1929, 2015.
Rumpker, G. and Silver, P. G.: Apparent shearwave splitting parameters in the presence of vertically varying anisotropy, Geophys. J. Int., 135, 790–800, 1998.
Schaeffer, A. J., Lebedev, S., and Becker, T. W.: Azimuthal seismic anisotropy in the Earth's upper mantle and the thickness of tectonic plates. Geophys. J. Int., 207, 901–933, https://doi.org/10.1093/gji/ggw309, 2016.
Schuberth, B. S. A., Zaroli, C., and Nolet, G.: Synthetic seismograms for a synthetic Earth: longperiod P and Swave traveltime variations can be explained by temperature alone, Geophys. J. Int., 188, 1393–1412, https://doi.org/10.1111/j.1365246X.2011.05333.x, 2012.
SchultePelkum, V. and Blackman, D. K.: A synthesis of seismic P and S anisotropy, J. Geophys. Int., 154, 166–178, 2003.
Shea Jr., W. T. and Kronenberg, A. K.: Strength and anisotropy of foliated rocks with varied mica contents, J. Struct. Geol., 15, 1097–1121, https://doi.org/10.1016/01918141(93)901587, 1993.
Sinogeikin, S. and Bass, J: Elasticity of pyrope and majorite – pyrope solid solutions to high temperatures, Earth Planet. Sc. Lett., 203, 549–555, 2002.
Sinogeikin, S., Bass, J., and Katsura, K.: Singlecrystal elasticity of ringwoodite to high pressures and high temperatures: implications for 520 km seismic discontinuity, Phys. Earth Planet. Int., 136, 41–66, 2003.
Steinberger, B. and Calderwood, A.: Models of largescale viscous flow in the Earth's mantle with constraints from mineral physics and surface observations, Geophys. J. Int., 167, 1461–1481, https://doi.org/10.1111/j.1365246X.2006.03131.x, 2006.
Sturgeon, W., Ferreira, A. M. G., Faccenda, M., Chang, S.J., and Schardong, L.: On the origin of radial anisotropy near subducted slabs in the midmantle, Geochem. Geophy. Geosy., 20, 5105–5125, https://doi.org/10.1029/2019GC008462, 2019.
Styles, E., Davies, D. R., and Goes, S.: Mapping spherical seismic into physical structure: Biases from 3D phasetransition and thermal boundarylayer heterogeneity, Geophys. J. Int., 184, 1371–1378, https://doi.org/10.1111/j.1365246X.2010.04914.x , 2011.
Tape, W. and Tape, C.: A reformulation of the Browaeys and Chevrot decomposition of elastic maps, J. Elast., 156, 415–454, https://doi.org/10.1007/s1065902410056x, 2024.
Tasaka, M., Zimmerman, M. E., Kohlstedt, D. L., Stünitz, H., and Heilbronner, R.: Rheological weakening of olivine + orthopyroxene aggregates due to phase mixing: Part 2. Microstructural development, J. Geophys. Res.Solid, 122, 7597–7612, https://doi.org/10.1002/2017JB014311, 2017.
Thomsen, L.: Weak elastic anisotropy, Geophysics, 51, 1954–1966, 1986.
Tommasi, A., Knoll, M., Vauchez, A., Signorelli, J. W., Thoraval, C., and Logé, R.: Structural reactivation in plate tectonics controlled by olivine crystal anisotropy, Nat. Geosci., 2, 423–427, https://doi.org/10.1038/ngeo528, 2009.
VanderBeek, B. P.: ECOMAN2.045 seismology PSI_D, Zenodo [code], https://doi.org/10.5281/zenodo.11186805, 2024.
VanderBeek, B. P. and Faccenda, M.: Imaging upper mantle anisotropy with teleseismic Pwave delays: insights from tomographic reconstructions of subduction simulations, Geophys. J. Int., 225, 2097–2119, 2021.
VanderBeek, B. P., Lo Bue, R., Rappisi, F., and Faccenda, M.: Imaging upper mantle anisotropy with traveltime and splitting intensity observations from teleseismic shear waves: Insights from tomographic reconstructions of subduction simulations, Geophys. J. Int., 235, 26402670, https://doi.org/10.1093/gji/ggad389, 2023.
Wang, Z. and Zhao, D.: 3D anisotropic structure of the Japan subduction zone, Sci. Adv., 7, 4, https://doi.org/10.1126/sciadv.abc9620, 2021.
Warren, J. M., Hirth, G., and Kelemen, P. B.: Evolution of olivine preferred orientation during simple shear in the mantle, Earth Planet. Sc. Lett. 272, 501–512, https://doi.org/10.1016/j.epsl.2008.03.063, 2008.
Wentzcovitch, R., Karki, B., Cococcioni, M., and de Gironcoli, S.: Properties of MgSiO_{3}perovskite: Insights on the nature of the earth's lower mantle, Phys. Rev. Lett., 92, 018501, https://doi.org/10.1103/PhysRevLett.92.018501, 2004.
Wookey, J.: Direct probabilistic inversion of shear wave data for seismic anisotropy, Geophys. J. Int., 189, 1025–1037, https://doi.org/10.1111/j.1365246X.2012.05405.x, 2012.
Workman, R. and Hart, S.: Major and trace element composition of the depleted morb mantle, Earth Planet. Sc. Lett., 231, 53–72, 2005.
Zha, C., Duffy, T., Mao, H., Downs, R., Hemley, R., and Weidner, D.: Singlecrystal elasticity of βMg_{2}SiO_{4} seismic discontinuity in the earth's mantle, Earth Planet. Sc. Lett., 147, E9–E15, 1997.
Zha, C., Duffy, T., Downs, R., Mao, H., and Hemley, R.: Brillouin scattering and xray diffraction of San Carlos olivine: direct pressure determination to 32 GPa, Earth Planet. Sc. Lett., 159, 25–33, 1998.
Zhang, Z., Stixrude, L., and Brodholt, J.: Elastic properties of MgSiO_{3}perovskite under lower mantle conditions and the composition of the deep earth, Earth Planet. Sc. Lett., 379, 1–12, 2013.
Zhang, Z., Cottaar, S., Tao, L., Stackhouse, S., and Militzer, B.: Highpressure, temperature elasticity of Fe and Albearing MgSiO_{3}: Implications for the earth's lower mantle, Earth Planet. Sc. Lett., 434, 264–273, 2016.
Zhou, Q., Hu, J., Liu, L., Chaparro, T., Stegman, D. R., and Faccenda, M.: Western U.S. seismic anisotropy revealing complex mantle dynamics, Earth Planet. Sc. Lett. 500, 156167, https://doi.org/10.1016/j.epsl.2018.08.015, 2018.
Zhou, W.Y., Hao, M., Zhang, J., Chen, B., Wang, R., and Schamndt, B.: Constraining composition and temperature variations in the mantle transition zone, Nat. Commun., 13, 1094, https://doi.org/10.1038/s41467022287097, 2022.
 Abstract
 Introduction
 Software package structure
 Discussion
 Conclusions and outlook
 Appendix A: DREX model vs. analytical solution
 Appendix B: Aggregate elastic stiffness tensor
 Appendix C: Velocity gradient tensor in polar coordinates
 Appendix D: Aggregate rotation due to fluid body rotation
 Code availability
 Data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Supplement
 Abstract
 Introduction
 Software package structure
 Discussion
 Conclusions and outlook
 Appendix A: DREX model vs. analytical solution
 Appendix B: Aggregate elastic stiffness tensor
 Appendix C: Velocity gradient tensor in polar coordinates
 Appendix D: Aggregate rotation due to fluid body rotation
 Code availability
 Data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Supplement