Analytical and numerical models of viscous anisotropy: A toolset to 1 constrain the role of mechanical anisotropy for regional tectonics and 2 fault loading

mechanical anisotropy


Introduction
Mechanical anisotropy can refer to either elastic moduli or creep viscosities depending on the style and orientation of deformation.The former is important for seismic wave propagation, but the viscous, long-term deformation type of mechanical anisotropy may be important for geodynamic processes, which is the focus of this study.
The resulting mechanical anisotropy can be preserved at distributed lithospheric scale within presently inactive, formerly deformed suture, i.e., tectonic inheritance, or concentrated into narrow shear zones within active plate boundaries (Vauchez et al., 1998;Mühlhaus et al., 2004).Spatial variations in mechanical anisotropy may result in strain localization in plate interiors that may affect flexural strength (e.g., Simons and van der Hilst, 2003) or play a role for intraplate seismicity (Mameri et al., 2021).
Olivine-aggregate deformation experiments show textures with significant viscous anisotropy (e.g., Hansen et al., 2016).Mechanical anisotropy is thus expected as a result of CPOs, and the development of the latter is explored widely in the context of connecting mantle flow and seismic anisotropy (e.g., Becker and Lebedev, 2021).Any feedback between mechanical anisotropy and convection may then affect the predictions for seismic anisotropy, for example (e.g., Chastel et al., 1993;Blackman et al., 2017).However, at least within an instantaneous mantle flow or lithospheric deformation scenario, mechanical anisotropy can be hard to distinguish from isotropic weakening (Becker andKawakatsu, 2011, Ghosh et al., 2013).Time-dependent scenarios of deformation are expected to be more modified by mechanical anisotropy compared to isotropic zones of weakness, e.g. for lithospheric instabilities and shear zones (Mühlhaus et al., 2004, Lev and Hager, 2008, 2011, Perry-Houts and Karlstrom, 2019), for post-glacial rebound (Schmeling, 1985, Han andWahr, 1997), or on plate scales (Honda, 1986, Christensen, 1987, Király et al., 2021).
It is thus important to further constrain the role of mechanical anisotropy for the lithosphere, and observations from tectonically well constrained regional settings provide an opportunity to explore complementary strain and stress sensitive data (e.g., Mameri et al., 2021, Schulte-Pelkum et al., 2021).In turn, mechanical anisotropy may affect some of the methods used to infer stress or stressing rate close to faults, such as inversion of focal mechanisms (e.g.Kaven et al., 2011).In Southern California, for example, inherited CPOs and alignment of weak layers through SPO could both be a source of mechanical anisotropy.This could possibly explain some of the mismatch between geodetically inferred strainrates and focal-mechanism derived stress close to faults, and the reactivation of preexisting fault structures may affect the tectonic deformation response and local fault loading (Schulte-Pelkum et al., 2021 and references therein).
Studies that explore the effects of mechanical anisotropy on regional scales for Southern California are, however, still limited.Ghosh et al. (2013) implemented an anisotropic San Andreas Fault (SAF) as a shear zone in a 3-D global, viscous deformation model but failed to identify robust indicators of mechanical anisotropy on regional scales.However, if mechanical anisotropy is considered in a regional scale model, it may be easier to assess the documented non-coaxiality between stress and strain (Schulte-Pelkum et al., 2021), and to eventually incorporate time dependence in a field-observation validated way.This suggests an opportunity to develop new methods for inferring mechanical anisotropy from field observations and further constrain fault loading.
In this study, we work toward a theoretical framework and first solve analytically the deformation of a simple 2-D model with a viscously anisotropic layer which highlights some of its fundamental mechanical behavior.The solution shows stress heterogeneity, strain-rate enhancement, and non-coaxial principal stress and strain-rates inside the anisotropic layer and reveals the mechanics behind such heterogeneity.We explore how the orientation and strength of mechanical anisotropy affect the non-coaxiality, stress heterogeneity, and strain rate enhancement.Second, we present a new, open-source finiteelement tool, its validation against the analytical solution, and applications to more complex 3-D scenarios.Lastly, we discuss the implications and potential applications of the method and tools.
2 The 1-D analytical solution of a viscously anisotropic layer subjected to simple shearing Motivated by the not necessarily intuitive solutions produced by earlier numerical tests for mechanical anisotropy, e.g., based on our implementations (Moresi et al., 2003, Becker andKawakatsu, 2011), we proceed to solve analytically the incompressible Stokes flow equation for a layered model subjected to simple shearing over the thickness, where a central viscously anisotropic layer is sandwiched between two isotropic layers (Figure 1).

Governing equations and rheology
The general boundary-value problem of incompressible Stokes flow equation is described by the force balance for a continuum (eq. 1) and the incompressible fluid assumption (eq.2) at any point in a domain Ω, where  is the stress tensor,  is the body force, and  is the velocity field.We use an incompressible, Newtonian flow constitutive law such that where p is pressure,  the deviatoric stress tensor, D the 4 th -order viscosity tensor, I the identity matrix, and ̇ the strain-rate tensor.
For isotropic and anisotropic domains, the viscosity  will be  iso and  ani , respectively.In the isotropic domains, with scalar dynamic viscosity .In the anisotropic domains, Here we solve a system with the hexagonal anisotropy following formulations in Mühlhaus et al. (2002) and Moresi and Mühlhaus (2006) (MM hexagonal anisotropy) with  the "director" of the weak viscous direction.Following eq. ( 3) in Mühlhaus et al. (2002), where in   ( = , ) is the components of the normal "director",  is the 'normal' shear viscosity, and   is the weak shear viscosity along the weak layer., , ,  = , .As shown in Figure 1,  is the angle between  and axis y, and then   = − (),   =  () (cf.Christensen, 1985).
A general set of boundary conditions on the boundary Ω = Γ  ∪ Γ  is given by where Γ  and Γ  stand for Dirichlet boundary and Newmann boundary, respectively, and   is the normal to Γ  .

Solution specifics
For our example problem, we chose as boundary conditions where a horizontal velocity   0 is applied to the top side, no velocity at the bottom, and periodic velocity and pressure on the west and east sides.Given the symmetry of model geometry and boundary conditions along x, the velocity, pressure, and stress are invariant along x, and vertical velocity is zero, which give   = 0;  , = 0;  , = 0;  , = 0 (11) where, for example,  , stands for

𝜕𝑥
, and ,  = , .Therefore, we solve the 1-D analytical solution of velocity, pressure, and stress along the vertical thickness (y axis).
Substituting eq. ( 11) into eq.( 5), we get ẋ  =  , = 0 (12a) In the isotropic layer, the deviatoric stress components follow as In the anisotropic layer, following eq.( 8), the deviatoric stress components are The task now is to find solutions of velocity gradients  , in the isotropic ( 1 ) and anisotropic ( 2 ) layers.Eq. ( 12) gives   =   = 0,   =  1 (15) and eq.( 13) yields The continuity condition for shear stress   and normal stress   +  on the interfaces between the isotropic and anisotropic layers require where  iso and  aniso are pressures inside the isotropic and anisotropic layers, respectively.
The boundary condition for   ( = 0) =   0 and   ( = −) = 0 and the integration of  , over the entire thickness w can be expressed as  17) and ( 20), we get Substituting  1 and  2 to eqs.(15,16,18), we get solutions for velocities, stresses, and pressure as a function of thickness y.Substituting  1 and  2 to eq. ( 11), we get the expressions for shear strain-rate in the isotropic and anisotropic layers as We use the square root of the  2 , deviatoric invariant of strain-rate tensor to measure the deformation, and in 2-D Then, in the isotropic and anisotropic layers, We define the ratio between square root of  2 invariant of the strain-rate tensor in anisotropic and isotropic layers  as strain-rate enhancement to measure the heterogeneity of deformation caused by mechanical anisotropy, and

The character of the analytical solution
We compute a scenario with  = 1,  = 1,   0 = 1, and  = 0.4 (thickness between −0.1 and −0.5) with variables defined as in Figure 1.We change the director  of the weak viscous direction by varying  from 0° to 90°, and the viscosity contrast  = /  in the anisotropic layer to explore their effects on stress and strain-rate.We first set  = 10.  Figure 3a shows the angles between  1 , 1, and  as a function of  in the anisotropic layer for  of 2, 10, and 100, respectively. 1 and  2 are angles between  1 and n, and between 1 and n, respectively.The mismatch  =  1 −  2 .For all s,  increases with increasing starting from 0°, reaches to a maximum, and then decreases to 0° when  reaches 45°.
The maximum  depends on viscosity contrast .With the larger  of 100, the maximum The maximum  for a wider range of  and the corresponding  that this maximum  is achieved is shown in Figure 4.If  is close to 1,  will approach to zero and the model recovers the isotropic scenario.If  increases,  will increase to the maximum 45° when  approaches to zero, akin to deformation along the weak anisotropic direction being a stressfree boundary.For  = 10, perhaps appropriate for olivine CPOs (Hansen et al., 2012), the maximum angular mismatch  could be as large as about 27.45° when  = 8.8°.Figures 3b and c show the maximum shear stress   max and pressure  in the anisotropic layer and the isotropic layer, respectively, as a function of  and . Figure 3d shows the difference between Figures 3b and c, and the difference shows similar trends as to the mismatch  that increases to a maximum and then decreases to zero when  varies from 0° to 45°.For  = 2, 10, and 100, the difference of   max is 0.05, 0.31, and 0.45, which occur when  = 18.8°, 13.5°, and 11.6°, respectively.
Figure 4. Maximum angular mismatch  between principal stress  1 and principal strain rate ̇1 as a function of viscosity contrast .For each ,  defines the normal vector of weak anisotropic direction at which the maximum  occurs.
The weak viscous anisotropy enhances strain-rate in the anisotropic layer.The enhancement can be measured by , the strain-rate enhancement as defined in eq. ( 27).
Figure 5 shows the normalized strain-rate enhancement /, caused by various viscosity contrast s as a function of .The maximum strain-rate enhancement occurs when  = 0° with a normalized value of unity, i.e., the enhancement  = .The strain-rate enhancement decreases with increasing  until there is no strain-rate enhancement with  = 1 when  = 45°.
where  and  are bilinear and linear terms of the variational formulation,  is the flux on the Newmann boundary.
Following the Stokes tutorial, the sign of pressure is flipped from the strong form given above.The purpose is to have a symmetric but not positive-definite system of equations in the finite-element implementation, which can be solved iteratively after properly preconditioning of the system.We precondition the linear system of equations with the preconditioner defined as Viscous anisotropy can be decomposed into components with different symmetries, e.g., similarly to what was explored by Browaeys and Chevrot (2004) for elastic anisotropy in the Voigt approximation.Here we derive and compare 3-D mathematical formulations of hexagonal anisotropy, which describe physical structures with a weak plane as shown in MM hexagonal anisotropy, and orthorhombic anisotropy, which is a closer approximation to full crystal structure of olivine that dominates the upper mantle, here modeled under the incompressible fluid assumption.
We define local material coordinate system with axes 1, 2, 3, and finite-element coordinate system with axes x, y, z.To simplify the structure of the 4 th order viscosity tensor expressed as a 6 × 6 Voigt matrix form, axes to symmetry planes in viscosity are aligned with axes 1,2,3.Different formulations for hexagonal viscous anisotropy are in use.With the deviatoric stress vector and strain rate tensor defined as  = ( 11 ,  22 ,  33 ,  23 ,  13 ,  12 ) and ̇= (̇1 1 , ̇2 2 , ̇3 3 , 2̇2 3 , 2̇1 3 , 2̇1 2 ), following eq.( 8), the Voigt form viscosity matrix  MM of MM hexagonal anisotropy is where  is a reference shear viscosity and   is the weak anisotropic viscosity.
Han and Wahr (1997) derive a hexagonal viscous anisotropy from a different method, and the Voigt form viscosity matrix  HW is where  1 ,  2 are isotropic shear viscosity, and weak shear anisotropic viscosity, respectively.And  1 (or  2 ) corresponds to 'normal' anisotropic viscosity (see, e.g., Christensen 1987).Not all four non-zero parameters are independent.Following the derivations in Han and Wahr (1997),  =  HW ̇ gives 22 = ( 2 + 2 2 )̇2 2 (32b) The incompressible fluid assumption is, and zero of the trace of deviatoric stress tensor gives Substituting eqs.(32) to eq. ( 34), we get To ensure eq. ( 33) is satisfied for any strain-rate tensor, eq. ( 35) gives The difference between  MM and  HW are the off-diagonal terms  13 HW and  31 HW .If  1 = 0,  HW collapses to  MM , that is MM hexagonal anisotropy is a simplified version of HW hexagonal without the correlation of deformation of normal strain-rates inside the weak plane.
For orthorhombic anisotropy, we add on top of  HW an additional orthorhombic component inferred from analogy to the orthorhombic elastic tensor in Browaeys and Chevrot (2004), which we define as where , , ,  are non-zero parameters.
Then, the orthorhombic viscosity matrix  ORTHOR is  ORTHOR =  HW +  ORTHOR (38a) The four non-zero parameters are not all independent given the incompressible fluid assumption.Following the same method above,  =  ORTHO ̇ gives Substituting eqs.( 39) into (34), we get To ensure eq. ( 33) is satisfied for any strain-rate tensor, and combining eq. ( 36),  =  = −.Therefore, of the four non-zero parameters, only  and  are independent.
Rotations of 4 th -order viscosity tensor are required to translate viscosity matrix from the material coordinate system to the finite-element one, and vice versa.In later 3-D scenarios with the anisotropic shear zone under simple shearing, we consider two elementary rotations of material coordinate system relative to the finite-element coordinate system, as shown in Figure 9. Axes 1, 2, 3 are originally aligned with axes x, y, z.For hexagonal anisotropy, axis 2 is the normal director to the weak viscosity plane.
For the first elementary rotation, axis 2 is rotated counterclockwise away from axis y around axis z(3) for an angle of .This rotation is similar to the rotation of  in the 2-D analytical model.For the second elementary rotation, axes 1 and 3 are further rotated around axis 2 counterclockwise for an angle of .
In the following sections, we first verify the finite-element implementation against the analytical solution by modeling the same problem presented in Section 2. We then increase the complexity slightly by introducing a Gaussian distribution of weak anisotropy across the thickness of the anisotropic layer.We next simulate a set of 2-D models inspired by a vertical fossil shear zone subjected to misoriented shortening to explore the strain-rate enhancement caused by the mechanical anisotropy.Then, 3-D shear zones with orthorhombic and two forms of hexagonal anisotropy subjected to simple shearing are simulated.Lastly, we present results from a 3-D model inspired by the Leech River Schist above the Cascadia subduction zone (Bostock and Christensen, 2012, and references therein) under convergent margin loading conditions.

Verification of the FEniCS code against the analytical solution
We simulate the 2-D model in Figure 1 with our FEniCS code and verify the implementation against analytical solutions derived in Section 2. Figure 6 shows matching FEniCS and analytical solutions for velocity, strain-rate enhancement, effective stress, and pressure over the whole thickness of the model, indicating that the code correctly implements this case of anisotropy.

"Fossil mantle" shear zone subjected to misoriented shortening.
We now consider strain-rate enhancements from a set of models with anisotropic shear zones subjected to misoriented shortening, partially inspired by the work of Mameri et al. (2021) and our earlier exploration of potential signals of mechanical anisotropy in southern California (Schulte-Pelkum et al., 2021).
The anisotropic shear zone is characterized by MM hexagonal anisotropy with the weak plane aligned with the strike of the shear zone.We simulate the deformation and stress/pressure from 2-D models of 2.5 by 1 along  and  directions, respectively, with viscosity contrast  = 10.The shear zone width is 0.1 and it is striking at an angle of  to the unit shortening (  = 1) along  on the west side (Figure 7a).The east side is free slip.
For the north and south sides, two scenarios are considered.In the Free Sides scenario, both sides are free, which simulates the extreme condition that the interacting blocks  The weak viscosity in the shear zone enhances strain-rates.The enhancement depends on the style of rheology and boundary conditions.Figure 8 shows strain-rate enhancement caused by the weak shear zone for various s, the angle between the normal to the shear zone strike and the horizontal shortening.The strain-rate enhancement is calculated by the average of square root of J2 invariant of the strain rate tensor in the shear zone divided the average outside of the shear zone along a horizontal profile.For Free Sides scenarios, if the shear zone is MM hexagonal anisotropy, the maximum strain-rate enhancement reaches 10, the same as the viscosity contrast  = 10 given, when  = 45°.If the shear zone is isotropic weak  iso = 0.1 , the maximum strain-rate enhancement is ~5.4.Either by increasing or decreasing  away from 45°, strain-rate enhancement decreases.
The maximum strain-rate enhancement with the isotropic weak shear zone is lower than for the MM hexagonal anisotropy due to lower shear stress along the inclined shear zone.
The driving force is normal stress   , which mainly affects flow ẋ  through the corresponding normal viscosity.In the isotropic weak shear zone, not only the shear viscosity is lower than the isotropic surrounding, as in the MM hexagonal anisotropic shear zone, but also the normal viscosities are lower than those in both the isotropic surrounding and MM shear zone.As a result, stresses and pressure are heterogeneous across the shear zone in the isotropic weak scenario while they are homogenous for MM scenario.In particular,   is lower inside the isotropic shear zone, which leads to lower shear stress along the inclined shear zone.The boundary conditions also matter.Mameri et al. (2021) discussed the effect of boundary conditions with free slip/lithospheric pressure conditions given their viscoelastic rheology.
In our models, the north and south sides in Pure Shear scenarios are more restricted compared to Free Sides scenarios where material is free to flow along the shear zone and outwards the north and south sides.As shown in Figure 8, for either anisotropic or isotropic weak shear zone, Pure Shear scenarios give less strain-rate enhancement compared to Free Sides scenarios.The maximum strain-rate enhancement occurs when  = 65° and it decreases with decreasing .Pure Shear isotropic weak shear zone produces less strainrate enhancement compared to anisotropic scenarios.

3-D shear zone with hexagonal and orthorhombic anisotropy under simple shearing
We simulate 3-D shear zones with MM and HW hexagonal anisotropy and orthorhombic anisotropy under simple shearing.Figure 9   Following the decomposition method in Browaeys and Chevrot (2004), we can compute the contributions to viscosity from isotropic, hexagonal, and orthorhombic symmetries.
Tetragonal and other lower symmetries such triclinic and monoclinic in the viscosity are not included in this study.As a demonstration, we choose  = 1,   = 0.1,  1 = 0.3,  = 0.6, and  = 0, which parameters give ~76% isotropic and ~24% hexagonal component weights for MM hexagonal anisotropy, and ~70% isotropic and ~21% hexagonal and 9% orthorhombic component weights for ORTHOR anisotropy, analogous to the composition of elasticity tensor of olivine.
We simulate models with  from 0° to 90° at 10° step size, and  from 0° to 90° at 15° step size.Figure 10a and b show the mismatch of principal stress and strain rate axes at the center (x = 0.5, y = 0.5, z = 0.5) of the anisotropic zone for ORTHOR and MM anisotropy, respectively.For  = 0° or 90°, the mismatch is zero for both ORTHO and MM anisotropy,

Leech River Schist above the Cascadia subduction zone
We expect that viscous anisotropy may arise from structural anisotropy like schist, rocks that has highly developed layered textures, which are generally exposed and associated with subduction zone environments (e.g., Chapman et al., 2010, Bostock and Christensen, 2012, Chapman, 2016, Xia and Platt, 2017).It appears the schist may overlap on top of the subducting oceanic plate as reconstructed geologically in the southern California case (Xia and Platt, 2017), though the schists were transferred to shallow depth in subsequent geologic episodes.If viscous anisotropy may cause non-coaxial stress/strain-rate axes and significant stress heterogeneity and enhance strain-rates as we demonstrate in previous theoretical setups, the migration of schist and its close relation to subduction zones may play an important role in the tectonic deformation of the lithosphere.Here, we focus our attention to the non-coaxially of stress strain-rate axes from a regional wedge-shaped schist structure subjected to subducting loading.
In Cascadia between southern Puget Sound and central Vancouver Island, the Leech River Schist (LRS), which is bounded by two north dipping thrusts forming a wedge (Bostock and Christensen, 2012, and references therein).The LRS rides on top of the subducting Juan de Fuca plate relative to North America.The schistosity, which is the parallel alignment of platy mineral constituents that reflects a considerable intensity of metamorphism, is generally west-east and vertically dipping and the relative plate motion direction is N56°E (Bostock and Christensen, 2012).Figure 12 shows a finite-element model and boundary conditions inspired by the LRS.The model domain is dimensionless and 10 by 10 by 3 along x, y, and z, respectively.The grid size inside the schist wedge is 0.1, which gradually increases to 1 near the model boundaries.The schist wedge is 2 by 1 on the free surface and vanishes at depth of −1.The schist is assumed to be with MM hexagonal anisotropy and the weak viscosity is aligned with the general strike of the schist, which is ~60° relative to the y axis.The viscosity contrast is 10. Figure 13 presents the principal stress and strain-rate axes on three orthogonal crosssections, x-y plane at z = −0.5, y-z plane at x = 5, and x-z plane at y = 5, that cut through the schist, respectively.The subducting loading and the wedge shape of the anisotropic regime are different from previous models and produce different stress and strain-rate axes patterns.
In map view (Figure 13a), the whole schist shows non-coaxial stress and strain-rate axes with mismatch angles about 27 − 30°.Strain-rate axes inside the anisotropic zone are largely aligned with those in the isotropic regime.The stress axes, on the other hand, are rotated away from those in the isotropic regime.The side view on the yz plane (Figure 13b) also shows significant stress and strain-rate non-coaxiality with mismatch angles increase with depth.The mismatch could reach a notable 90° near the sharp wedge bottom.The other side view on xz plane (Figure 13c) shows very limited angular mismatch of just a few degrees, when the subduction is near parallel to the weak direction of the anisotropy.In addition, the stress and strain-rate axes dip out of the horizontal plane.The implication is that loading style and the shape of anisotropic structure could be important in producing mismatch between principal stress and strain-rate axes, and dipping principal axes.The results assume that the schist can be approximated with hexagonal viscous anisotropy and the deformation and stress features reflect the current loading condition.The schist may, of course, carry stress and strain signatures inherited from previous tectonic episodes and is subjected to temporal change depending on the viscosity of the structure and the time length scale of interest.Further exploration of observations of stress and strain-rate orientations associated with the structure and a suite of models that have various viscosity contrasts would be helpful to differentiate signatures from present and inherited.

An approach to constrain viscous anisotropy
The difference of stress and pressure between the isotropic and anisotropic layers could influence mechanical processes in such a system like a fault zone (e.g., Hardebeck andMichael, 2004, Hirano andYamashita, 2011).Non-coaxiality between principal stress and strain-rate axes from viscous anisotropy, such as due to SPOs and CPOs, could be assessed quantitatively, and they can infer stress and pressure heterogeneity.This motivates reassessment of independent measures for inferring stress or stressing-rates (e.g., Michael, 1984) and strain-rates derived from geodetic constraints (e.g., Smith-Konter and Sandwell, 2009).Close to faults in southern California, the two fields match in their alignment on broad scales, but there are also significant local deviations (Becker et al., 2005, Yang and Hauksson, 2013, Schulte-Pelkum et al., 2021, Johnson, 2024) which are expected to be of relevance for long-term tectonics as well as setting local stress conditions for earthquake rupture.One interpretation was deformation memory from the Farallon subduction and subsequent extension.

Schulte
Notably, a comparison of focal mechanism-based principal stress axes (Yang and Hauksson, 2013) with GNSS-derived principal strain rates (Sandwell et al., 2016) shows an angular mismatch with a peaked distribution centered on an azimuth (CW from N) of −6° with a standard deviation of 19° (Schulte-Pelkum et al., 2021).Based on our results (Figure 3a), the observations may indicate mild mechanical anisotropy of viscosity contrast of 2 to 10 in the region for nearly all the  if we assume the weak anisotropy were parallel to the simple shearing loading.The higher viscosity contrast of 100 is also possible if 20°<  < 70°.It could be also possible that the anisotropic structure is subjected to misoriented shortening or additional factors should be considered such as more complex loading conditions, special shapes of structures, inheritance from previous geodynamical processes, and combinations of any few.For misoriented orthorhombic anisotropy or the case of Leech River Schist where the loading is oblique to anisotropic regime with special shape, dips of principal axes could be used to infer mechanical anisotropy if they were measurable.Alternative sources that can help narrow down candidate scenarios are helpful.
The non-coaxiality of principal stress and strain-rate is more visible if the loading direction is misoriented from the weak anisotropic direction (cf.Ghosh et al., 2013).The case of Leech River Schist and the structure in southern California illustrate that the combining condition of misoriented loading and weak anisotropy (such as schistosity) may be common in nature.In addition to non-coaxial principal axes, heterogeneity of stress and pressure, and enhanced strain-rate may occur as well.For example, using teleseismic receiver functions, Audet (2015) finds that the plane of fast velocity strikes parallel to the San Andreas fault while dipping mildly throughout the crust near Parkfield.He interprets the mid-crustal anisotropy as fossilized fabric within fluid-rich foliated mica schists.Our results suggest that heterogeneity of stress and pressure might indeed be induced by the mechanical anisotropy of the schist, which could influence the stress distribution in the region and nearby earthquakes.

Conclusion
We present a 1-D analytical solution to a viscously anisotropic layer subjected to simple shearing which predicts significant stress heterogeneity and non-coaxial stress and strain rates.Observations of the non-coaxiality and dips of principal axes could give us constraints on mechanical anisotropy in nature.Such analysis may be possible, e.g., by comparing stress inversions from focal mechanisms, surface strain-rates from geodetic measurements, and integrated strain from seismic anisotropy (Schulte-Pelkum et al., 2021, and references therein).
To accelerate such studies, we develop an open-source finite-element code using FEniCS, verify the 2-D version of the code against the analytical solution, and explore a number of 2-D and 3-D illustrative cases with various loading styles, hexagonal and orthorhombic anisotropy, and the wedged shape Leech River Schist above the Cascadia subduction zone.
We hope that this exploration of mechanical anisotropy for tectonic problems and our new implementation will help advance model and verification of mechanically anisotropic lithospheric models, and their implications, from long-term plate boundary evolution to fault loading and rupture propagation.

Figure 1 .
Figure 1.Schematic diagram of the 2-D layered model with a viscously anisotropic layer subjected to simple shearing. is the "director" of weak viscous ( weak ) direction.The viscosity of the strong direction in the anisotropic layer and the isotropic viscosity are  strong and  iso , respectively.The model domain is L by w with the anisotropic layer with a thickness of d.The angle  is counted counterclockwise from the y axis to .The bottom of the model is no slip, zero velocity.The top of the model shears horizontally with a velocity of   0 .Velocity and pressure on the west and east boundaries are periodic, and the 1-D analytical solution applies with thickness.

Figure 2
Figure 2 shows the maximum principal stress  1 (white bars) and maximum principal strain rate ̇1 (red bars) between −0.45 and −0.55 thickness, and the maximum shear stress   max (background) between −0.4 and −0.6 thickness, for various s.Sharp changes of physical quantities occur at the isotropic-anisotropic interface at −0.5 thickness.In the anisotropic layer, principal stress axes are mismatched at an angle  to the principal strainrate axes, which are always at 45° to the horizontal axis.The mismatch occurs for a wide range of  and the magnitude of  depends on  .The maximum  is ~27.45°.With increasing  from 0°,  increases from 0° to the peak of ~27.45° when  = 8.8°, and then decreases to 0° when  reaches 45°.When  further increases from 45°,  increases from 0° again to ~27.45° but with sign reversed until  = 81.2°,then decreases to 0° when  reaches 90°.

Figure 2 .
Figure 2. Principal stress  1 (white bars), principal strain rate 1 (red bars), and maximum shear stress   max (background with the colorbar) as a function of  with viscosity contrast of 10.The isotropic-anisotropic interface is at -0.5 thickness, and the domain above is anisotropic and below is isotropic, as indicated by 'ani' and 'iso', respectively.

Figure 3 .
Figure 3. (a) Angular relations between principal stress  1 , principal strain rate 1, and the normal director n of the weak anisotropic viscosity for three viscosity contrasts  s.Maximum shear stress and pressure as a function of  in the anisotropic (b) and isotropic layer (c) for three  values.(d) The difference between (b) and (c).

Figure 5 .
Figure5.Normalized strain-rate enhancement / for various  s and  s.Srain-rate enhancement  and viscosity contrast  are defined in eq.(27).3Numerical solutions for 2-D and 3-D problems3.1 Overview of the finite-element method and formulations of various viscous anisotropyFor increased transparency, accessibility, and expandability for more complicated 2-D and 3-D scenarios, including for regional settings, we develop a new finite-element code using the open-source computing platform FEniCS with a user-friendly Python interface(Logg et al., 2012, Logg andWells, 2010) (https://fenicsproject.org/) to simulate incompressible Stokes flow with viscous anisotropy.The finite-element implementation follows the FEniCS Stokes tutorial (link provided in the Data and Software Statement).The material matrix for viscous anisotropy is fully expressed by 4 th -order tensors through a set of Python functions, which currently support hexagonal and orthorhombic anisotropy, and can be readily expanded to anisotropy with more general symmetries.For the choices of function spaces, we use second-order Continuous Galerkin (CG2) elements for velocity, and first-order Continuous Galerkin (CG1) elements for pressure in 2-D.For 3-D problems, we use third-order Continuous Galerkin (CG3) elements for velocity, and second-order Discontinuous Galerkin (DG2) elements for pressure.The choices of the function space pairs satisfy the Ladyzhenskaya-Babuška-Brezzi (or inf-sup) compatibility condition (seeBrezzi and Fortin (1991) for more details).The theoretical considerations behind the choices are described in Chapter 20 inLogg et al. (2011) and references therein.We use built-in mesh generator of FEniCS with triangles in 2-D and tetrahedrals in 3-D for simple model geometries, and the open-source mesh generator Gmsh(Geuzaine and Remacle, 2009) (https://gmsh.info/) for more complicated model geometries.FEniCS provides API to Gmsh for a seamless integration of the two tools.

Figure 6
Figure 6 also shows results of a scenario with Gaussian distribution of weak anisotropy where   = 1 − (1 − 1  )  (− ( −  ℎ ) 2 ), perhaps closer to what might be expected in a natural shear zone.Here,   = −0.7 is the thickness at the center of the anisotropic layer, ℎ = 0.1, and  is 10.  is 1  = 0.1 at   , and about unity, i.e. the isotropic shear viscosity,when  approaches the edges of the anisotropic layer ( = −0.5 and  = −0.9).The   in the Gaussian scenario is mostly larger than the constant 0.1 in the analytical solution over the anisotropic layer.Therefore, amplitudes of heterogeneities of strain-rate enhancement, stress and pressure are less pronounced compared to the analytical solution and the peaks occur within a narrower thickness.

Figure 6 .
Figure 6.Verification of FEniCS finite-element solution against and analytical solution for horizontal velocity,   , strain-rate enhancement, effective stress   + , and pressure , over thickness.Results with weak anisotropy following a Gaussian distribution in the anisotropic layer are in red lines. denotes the orientation of weak anisotropy director defined in Figure 1.
outside of the north and south of the domain are extremely weak.In the Pure Shear scenario, the north and south sides extrude at absolute velocities of |  | = 0.2, simulating the other extreme condition that the interacting blocks are sufficiently strong compared to the simulated domain.Because we are solving incompressible Stokes flow, the extruding velocity of 0.2 is calculated by conserving the total volume.We vary  from 5° to 65° in 5° step size.We also consider scenarios with the shear zone to be isotropic but with weaker viscosity 1  = 0.1 than the surroundings.

Figure 7 .
Figure 7. (a) Schematic diagram of 2-D shear zone subjected to misoriented shortening.The west side has a unit shortening of   = 1 and the east side is free slip.The north and south sides are either free or extruding at a fixed velocity.The shear zone is at an angle of  to the unit shortening., the normal director to the weak anisotropy, is always normal to the shear zone strike.

Figure 8 .
Figure 8. Strain-rate enhancement caused by 2-D weak viscous shear zone subjected to misoriented shortening.
shows the unit box that has the anisotropic zone enclosed by isotropic layers.The north side has a unit velocity along x.The top, bottom, and south sides are free slip, and the east and west sides are periodic for both velocity and pressure.The volume of the model does not change, compatible to the incompressible fluid assumption.

Figure 9 .
Figure 9. Diagram of 3-D anisotropic shear zone under simple shearing.Two elementary rotations from local material coordinate system 1,2,3 that define the Voigt form of viscosity matrix, to finite-element coordinate system x, y, z are shown.
consistent with results from 2-D models.For other s but same , mismatch peaks at  = 10° or 80° and decreases when  changes toward 45°.The mismatch for MM anisotropy does not depend on , as expected from the fact that hexagonal anisotropy is isotropic inside the weak plane.The mismatch angles are the same as the 1-D analytical solutions for same s in Figure3a.In contrast, the mismatch for ORTHOR anisotropy depends on  and increases when  increases from 0° to 90° ( 33 ORTHOR <  22 ORTHOR <  11 ORTHOR ) for most s except for  = 40° or 50°.For one , the spread of mismatch for different s ranges from ~5° ( = 10° or 80°) to ~2°.HW hexagonal anisotropy gives the same mismatch angle results to MM anisotropy.

Figure 10 .
Figure 10.Angular mismatch of principal stress and strain-rate axes for orthorhombic (a) and Mühlhaus and Moresi hexagonal anisotropy (b) at the center of the anisotropic zone in the 3-D model subjected to simple shearing.In addition to the -dependence of mismatch for ORTHOR anisotropy, it tilts the principal stress and strain rate axes out of the horizontal x-y plane.Figure11aand b show the dip angles of axes of principal stress (a), and strain-rate (b) at the center of the ORTHOR anisotropic zone for s and s.The axes of principal stress do not dip much.Larger dips occur with  > 40°.The peak dip is ~2° when  = 80° and  = 30°/45° (Fig 12a).The dips of axes of principal strain rates show higher values when  < 60° with peak value at ~7° when  = 20° and  = 45° (Fig 12b).For hexagonal anisotropy, the principal axes all stay inside the horizontal x-y plane.

Figure 11 .
Figure 11.Dips of axes of principal stress (a), and strain-rate (b) at the center of the orthorhombic anisotropic zone for different s and s.

Figure 12 .
Figure 12.(a) Finite-element model of the Leech River Schist model.The schist is at the center of the model with west-east trending and vertically dipping schistosity.East is indicated.Dashed lines show the subducting of the Juan de Fuca plate.Except for the free surface, other boundaries are free slip.(b) Tetrahedral finite-element mesh generated by the open-source mesh generator Gmsh with refined mesh inside the schist.

Figure 13 .
Figure 13.Principal stress (black) and principal strain-rate (red) axes of a horizontal crosssection (a) at z = −0.5, of two vertical cross-sections (b) at x=5 and (c) at y =5 that cut through the Leech River Schist.
-Pelkum et al. (2021) discussed a wider range of deformation indicators for southern California from the surface to the asthenosphere mantle.They found general consistency with N-S compression and E-W extension near the surface and in the asthenospheric mantle, but all lithospheric anisotropy indicators deviate from such patterns.