Modification of crustal seismic anisotropy by geological structures (“structural geometric anisotropy”)

Seismic anisotropy is the cumulative interplay between propagating seismic waves and anisotropic earth material. Unraveling this effect in deformed crustal rocks is complex due to 3-D geological geometry and heterogeneity, seismic waves of limited bandwidth or selected travel paths, field experiments that may not offer full azimuth/inclination coverage, and the observation of anisotropy as often second-order waveform or traveltime features. Yet diverse studies of rock fabrics, metamorphic processes, and petrophysical properties associated with tectonic deformation reveal the common presence of strong elastic material anisotropy that will produce seismic anisotropy (e.g., Christensen, 1965; Kern and Wenk, 1990; Barruol and Mainprice, 1993; Christensen and Mooney, 1995; Godfrey et al., 2000; Cholach and Schmitt, 2006; Lloyd et al., 2009; Almqvist and Mainprice, 2017). Aligned cracks in the seismogenic crust have long been known to produce observable seismic anisotropy indicative of upper crustal stress fields near major faults or volcanic regions (e.g., Anderson et al., 1974; Hudson, 1981; Leary et al., 1990; Cochran et al., 2006; Savage et al., 2010; Almqvist and Mainprice, 2017).

Seismic anisotropy in the crust due to deformational and metamorphic fabrics has become increasingly observed using a variety of seismic body and surface wave phases (e.g., Brocher and Christensen, 1990; Carbonell and Smithson, 1991; Ozacar and Zandt, 2004; Sherrington et al., 2004; Readman et al., 2009; Lin et al., 2011; Bostock and Christensen, 2012; Okaya et al., 2016). More recently, studies that calculate seismic velocities from rock fabrics imaged with neutron or electron backscatter diffraction techniques suggest fabric anisotropy can have orthorhombic or lower symmetry (Ivankina et al., 2005; Lloyd et al., 2009, 2011). We show in this study that fabric-filled large-scale geological structures having 2-D or 3-D geometry can also produce lower orders of symmetry and hence can create seismic anisotropy signals that are inherently not two-cycle in azimuthal periodicity. Thus we conclude the starting assumption in seismic analysis methods that crustal anisotropy is slow-axis transverse isotropic is not always valid.

Recent studies in central Europe illustrate observations of two end-member scales of anisotropic structures: (1) bulk (uniform) regional medium and (2) structural modification of local rock anisotropy. For bulk medium, ´Sroda (2006) identified clear evidence of regional azimuthal P-wave anisotropy using CELEBRATION-2000 explosion first arrivals inverted for compressional velocity. This result was attributed to early Proterozoic meta-sedimentary layering tightly folded into a near-vertical orientation on a regional scale. Similarly, using 9-km-deep vertical seismic profiling data at the KTB deep drill hole site, Okaya et al. (2004) showed that despite significant interlayering (O’Brien et al., 1997), the high-grade metamorphic Bohemian massif in southern Germany can be characterized on a seismic scale as a regional bulk tilted transverse isotropic medium.

For the other end-member of structural modification, Bleibinhaus and Gebrande (2006) interpreted regional horizontal fast-axis symmetry from crustal seismic data in the Tauern Window, Swiss Alps. This symmetry was seemingly enigmatic because this region is geologically composed largely of (slow-axis symmetry) mica-rich phyllites and gneisses. These authors recognized that internal folding can redefine the regional-scale symmetry and thus alter the seismic response within their data. Figure 1A schematically illustrates how simple cylindrical folding of slow-axis symmetry fabric can become a net volume of fast-axis lower symmetry anisotropy. An implication is that the use of outcrop or hand sample-sized volumes to obtain seismic velocities via petrophysical or thin-section–based measurements may not be sufficient to characterize the regional anisotropy. The larger scale structural geometry of the anisotropic rocks needs to be assessed.

In this study, we show that the contributions of structural geometry to modify seismic anisotropy and as possible factors during seismic data analysis are significant. Here, we define “structural geometric anisotropy” (SGA) (see Table 1 for notation used within this study) as the alteration of seismic anisotropy that is produced by the 2-D or 3-D geometry of common crustal structures that contain internal rock fabric which locally is intrinsically anisotropic. Many common crustal rocks are elastically anisotropic owing to the inherent anisotropy of the constituent minerals as well as microstructural characteristics such as crystallographic or shape fabric. When these rocks are deformed into large structures such as folds, domes or shear zones, the anisotropic elastic properties of the rock are locally rotated across the structure. As a result, the structure’s larger scale volume can have its own bulk anisotropic properties. This effect is different from “structural anisotropy” that refers to seismic anisotropy produced by alignment of (isotropic) features or fine layering such as strata (e.g., Backus, 1962; Babuška and Cara, 1991).

Seismic waves are sensitive to volumes of rock at scales ranging from the 10s to 1000s of meters. For example, Figure 1B shows a km-scale antiform composed of Ragua schist in the Betic Mountains, southern Spain (Martínez-Martínez et al., 2002) that illustrates how structural geometry can affect the production of anisotropic seismic signals. Relatively short wavelength seismic waves may resolve individual parts of this structure. These waves may pick up different anisotropic characteristics depending on propagation direction across the fold limbs, parallel to the fold hinge (in-out of the photograph), or at an oblique angle to the limbs. However, an upcoming long-wavelength teleseismic wave will average through fine structure and just respond to the bulk structure. Representation of this fold volume as a long-wavelength effective medium will remove complexities of the finer scale details while still preserving the net anisotropic properties and resulting bulk seismic response.

We present in this study the role of large-scale geometry as applied to anisotropic fabrics and illustrate the SGA effects of several basic geometrical structures. We also use tensor algebra to show that under specific but major simplifying conditions, one can separate the geometrical orientation terms from a local already-anisotropic rock. This separation allows one to quantify the role of macroscale structure independent of the rock’s microscale material anisotropy when trying to understand the cause of observed seismic anisotropy signals. We demonstrate this separation using simple sinusoidal folds. Seismic compressional and shear wave velocities derived from SGA elasticity reveal complex angular variation in seismic anisotropy can exist. We illustrate these concepts with numerical (synthetic) anisotropic wave propagation examples. The main focus of our study is to understand the effect of structural geometry on seismic anisotropy; we are not solving for exact volume averaging or homogenization of folded structures.

Metamorphic and deformational processes produce fabrics such as schistosity that can be localized within deformational zones or form extensive penetrative regions of foliated rock. Laboratory investigations (e.g., Birch, 1958; Christensen, 1965; Fountain and Christensen, 1989; Weiss et al., 1999) show that most crustal metamorphic rocks are highly anisotropic due to preferred orientations of anisotropic silicate minerals. Material anisotropy can be severe; P-wave anisotropy in crustal foliated rocks (e.g., schists, gneisses, and amphibolites) may be as high as 17%–20% and S-wave anisotropy for phyllosilicate rocks as high as 25%–30% (Johnston and Christensen, 1995; Godfrey et al., 2000). In contrast, olivine-rich mantle rocks such as dunites average only 5%–6% (Christensen and Mooney, 1995). Crustal metamorphic rocks are often approximated as having slow-axis transverse isotropic symmetry.

A rock’s microstructure (or fabric) can be defined by the following characteristics of its constituent phases (minerals ± fluids): (1) spatial arrangement, (2) shapes, (3) modal abundances, (4) crystallographic preferred orientations and (5) shape preferred orientations. Each of these microstructural factors influences the bulk elastic properties of the rock and therefore must be taken into account to determine a precise rock stiffness tensor. These dependencies are caused by the cumulative effects of grain-scale mechanical interactions throughout a heterogeneous sample. We do not examine these microscale factors in this study and refer readers to studies of seismic anisotropy caused by microfabrics (e.g., Lloyd et al., 2009; Naus-Thijssen et al., 2011a, 2011b; Cyprych et al., 2017; summarized in Almqvist and Mainprice, 2017).

The alteration of a seismic signal that traveled through anisotropic earth can take different forms: traveltime advance or delay as a function of wavepath direction, alteration in particle motion, splitting of shear waves, or amplitude variation of a reflected or converted wave, all sensitive to wave direction. The wavelength of a propagating seismic wave compared to the length scale of the anisotropic earth feature directly affects the production of anisotropic signals. As is discussed below, the ratio of these length scales will determine if a wave accumulates signal along its path (e.g., raypath integration), responds to only a bulk average of fine-scale anisotropic earth, or may pass through and not “sense” any material anisotropy.

In addition, the amount of signal that a wave may pick up is related to the amount of exposure to anisotropic material in two ways: (a) the wave’s length of path in the material, or (b) the intensity of the material anisotropy. An azimuthal P-delay or shear wave splitting time can be represented as δt = Δpath*Δslowness. The longer the wave travels within anisotropic material, Δpath, the larger the produced seismic anisotropic signal. On the other hand, the strength of material anisotropy here is Δslowness = (1/V_slow–1/V_fast) which is related to anisotropy percentage. A split shear wave may be observable if path length or anisotropy percentage is large, even if the other is small. The relative raypath angles through the material are very important, particularly if the material changes strike or dip due to structure.

When a wave of length λ travels through material in which the scale of the anisotropic structure is L, then three cases can be considered (Fig. 2).

Case 1: L >> λ. When the anisotropic structure is much longer than the wave, the wave responds to each encountered segment of geology, and anisotropic signals will be accumulated along its path (e.g., path integration). The geology can be heterogeneous, the material anisotropy can change, and the seismic expression of anisotropy can change along the path.

Case 2: L ≈ λ. Anisotropic feature is about the same length of the wave. If the geological feature is isolated, then the seismic wave may pick up distinct identifiable phases, such as a reflection or transmitted conversion. This is similar to Case 1. If the feature is not isolated but is part of a more pervasive volume, the wave may respond to the net pattern. This is similar to case 3. An example of this latter situation is if a large layer contains local folding on the order of or shorter than the wave’s length. We also note that λ/4 to λ/8 is commonly used as the resolution limit for seismic reflections of thin layers (e.g., Widess, 1973).

Case 3: L << λ. the scale of anisotropic structure is much shorter than the wave (the long-wavelength case). The wave will be insensitive to individual small features but will respond to the net aggregate volume. This situation applies when an earth volume much larger than the wave is filled with small scale features. Complex earth volumes can be represented as effective media; simplified media that remove complexities but exhibit equivalent wave propagation behavior. Concepts of effective media for seismic wave propagation are well established with a wide range of applications. These applications include spatial averages of elastic heterogeneity using methods such as Voigt (1928), Reuss (1929), Hill (1952), or combinations thereof, or modern homogenization methods obtained from material sciences (e.g., Naus-Thijssen et al., 2011a, 2011b; Capdeville et al., 2010, 2013; Vel et al., 2016). Effective media theory is also applied to fractured rocks having crack alignment/density or fluids (e.g., O’Connell and Budiansky, 1974; Anderson et al., 1974; Budiansky and O’Connell, 1976, Crampin, 1981; Hudson, 1981; Douma, 1989); scattered media (Wu, 1989); and anisotropy due to subparallel thin layers that are isotropic (Backus, 1962) or anisotropic (Schoenberg and Muir, 1989). Equivalent media for aligned anisotropic mineral grains (crystallographic preferred orientations) can be determined using elastic tensor averaging methods (references within Mainprice, 2015); anisotropy of aligned isotropic bodies (shape preferred orientations) requires advanced methods of homogenization (e.g., Naus-Thijssen et al., 2011a, 2011b; Capdeville et al., 2010, 2013; Vel et al., 2016) or statistical analysis (Song and Jordan, 2017). Effective media derived from intrinsic anisotropic materials or from isotropic heterogeneity can produce similar seismic signals (Levshin and Ratnikova, 1984; Fichtner et al., 2013), highlighting the non-uniqueness of seismic waveform interpretation.

A primary condition for an effective medium is that the seismic wave under consideration is much longer than the scale of heterogeneity (i.e., Case 3). The heterogeneities can be compositional, and/or variations in degrees of anisotropy or structural orientation. Our definition of SGA applies only within this long-wavelength case. More specifically, we examine effective media for geometry of rocks that are intrinsically anisotropic. Our study here does not apply to structural geometry of isotropic rocks. We use established averaging and homogenization methods in this study in order to focus on structural orientation.

In this study we illustrate the effects of SGA as applied to local fabrics that are defined using real rock samples. While later in this study we use material symmetries that include full transverse isotropy and orthorhombic, we first employ elliptical transverse isotropy, one of the simplest symmetries, in order to demonstrate SGA effects without the complexities that lower order symmetries contribute. The rock samples we use represent three rock types: a schist (from the Orocopia Formation, southeast California), a slate (Central Range, Taiwan), and a gneiss (Nanga Parbat, Pakistan). Anisotropic seismic velocities and elastic stiffness tensors for these samples are listed in  Appendix 1, Tables A1 and A2, respectively. Laboratory acoustic measurements were obtained by N. Christensen (personal commun., 2017) for the schist and slate samples. Based on fabric and lineations, each has nine measured velocities that define orthorhombic symmetry (P-wave velocities in three axial and three diagonal directions plus S-wave velocities in three axial directions). As described in  Appendix 1, these orthorhombic measurements were simplified into full transverse isotropic (TI) symmetry by averaging the nine into five “measurement” velocities. This full TI symmetry was reduced to elliptical TI by further averaging from five to four velocities. The gneiss sample was measured by Meltzer and Christensen (2001) with five measurement velocities for full TI symmetry. We then defined its elliptical TI symmetry ( Appendix 1). Elastic tensors for all samples were calculated from the seismic velocities using Christoffel equations (Auld, 1973).

In the following sections we describe the effect of structural geometry on a volume filled with anisotropic rock fabric. For illustrative purposes we use the elliptical TI version of the schist as a representative fabric that becomes modified by structure. The seismic P-wave velocities (VP) of this schist measured at 200 MPa (∼7 km depth) are 5500 m/s normal to its foliation and 6865 m/s parallel to foliation; the corresponding S-wave velocities (VS) are 3160 and 4040 m/s, respectively (Fig. 3A;  Appendix 1). The diagonal VP value is set to the average of the two axial VP velocities based on the definition of elliptical TI symmetry (Auld, 1973). These velocities are used within Christoffel equations to obtain the velocities in all propagation directions (Auld, 1973). Figure 3A illustrates these velocity variations for propagation angles as referenced from the foliation normal which serves as the symmetry axis. The difference in the two VS curves indicates the possible production of shear wave splitting. The relative velocities in the axial directions define a characteristic of the anisotropy symmetry. When the velocities associated with the symmetry axis direction are slower than in the perpendicular direction (foliation plane), the symmetry is “slow-axis” (Figs. 3A, 3B). Levin and Park (1998) discussed the differences between slow- and fast-axis symmetries (e.g., their “pumpkins” and “melons”); Brownlee et al. (2017) provided additional characteristics of this terminology.

We replot these velocity curves in spherical coordinates using wave direction azimuth and inclination through the rock sample. The sample is oriented with its symmetry axis vertical and its foliation plane aligned horizontal with the “equator” of the sphere (“Orig” in Fig. 3B). Christoffel equations are used to solve for the P- and S-wave particle motions (e.g., Auld, 1973; Mainprice, 1990). We display the velocities and particle motions on the sphere surfaces by showing VP and translating the two VS into the equivalent amount of shear wave splitting (δt) per kilometer of propagation path. The S-wave particle motion bar (red) is the fast-S direction. Here, δt is an alternative to VS percentage as a way to illustrate the amount of VS anisotropy. The velocity patterns on the sphere surfaces are a direct expression of the underlying rock elastic tensor and easily convey the tensor’s anisotropy symmetry and tilt orientation.

With the rock sample oriented with its symmetry axis vertical, its velocity sphere patterns represent vertical-axis TI (VTI) symmetry or seismic radial anisotropy (Fig. 3B). The faster VP directions are clearly visible along the equator as are the larger S-wave splitting times. If the fabric and its elastic tensor are rotated by 90° so that the foliation plane is vertical, this orientation is horizontal-axis TI (HTI) symmetry or seismic azimuthal anisotropy. The 2θ azimuthal variations are noticeable not only along the equator but at oblique propagation directions. An inclined elastic tensor representing geological dip has tilted-axis TI (TTI) symmetry. Even though the fast zones are not aligned with the sphere axes, the characteristic of TI symmetry is easy to identify within the velocity patterns (Fig. 3B).

The elastic stiffness tensor for this representative fabric is: C₁₁ = 132.666, C₃₃ = 85.154, C₅₅ = 28.109, C₆₆ = 45.945, and C₁₃ = 46.372 (in GPa), using Voigt notation with indices mapped to coordinate axes using 1 = x, 2 = y, and 3 = z. Other elements are C₂₂ = C₁₁, C₄₄ = C₅₅, C₂₃ = C₁₃, and C₁₂ is calculated from C₁₁ and C₆₆ (see Appendix Table A2). All other elements are set to zero, and the tensor is diagonal symmetric. We apply elastic tensor decomposition in order to obtain a descriptive measure of the amounts of anisotropy within the tensor (Baerheim, 1993; Browaeys and Chevrot, 2004). This decomposition identifies tensor element patterns that match those of five orders of anisotropic symmetry (i.e., hexagonal, orthorhombic, tetragonal, monoclinic, and triclinic) after removing the isotropic contribution. Beneath each tensor sphere in Figure 3 and in the remaining figures, we provide three values that represent the TI, orthorhombic, and sum of all remaining lower symmetries as percentages of the original tensor. The net sum of these three values represents the total amount of anisotropy, and the net sum subtracted from 100% represents the isotropic component. The elastic tensor for the representative fabric shown in Figures 3A and 3B is 81.3% isotropic and 18.7% TI, with no amounts of orthorhombic or lower symmetries. The use of decomposition values when comparing two tensors allows for the difference in their anisotropy to be more easily identified.

In this section, we illustrate the modification to anisotropy that geometry can impart, using the example of a rock fabric that is sinusoidally folded (cylindrical fold). We define a layer filled with our representative fabric that starts out with horizontal foliation, then fold it and examine its bulk anisotropic properties. The sinusoidal folds are defined in the x-z cross-sectional plane using z(x) = Asin(2πx/L) where L is fold wavelength, A is amplitude, and its shape is characterized by the fold limb angle θ₀ = tan^–1(2πA/L). We assume the folds are infinitely long in the third dimension, y, defining cylindrical folds (Fig. 1A). We create a finely discretized numerical volume of length L with each subvolume containing our schist, tilted to conform to the shape of the fold. The bulk fold properties are represented by an averaged elastic tensor calculated by summing the subvolumes’ tilted tensors across one wavelength of fold.

For a fold described with fold limb angle θ₀ = 45°, we calculate its net volume tensor using the Voigt, Reuss, Hill, and asymptotic expansion homogenization (Naus-Thijssen et al. 2011a, 2011b; Vel et al., 2016) methods mentioned above. Figure 3C shows the seismic velocities obtained from these tensors. Symmetry decomposition percentage values (TI/orthorhombic/lower order) are also shown. There are three relevant observations. First, the anisotropy symmetry patterns for the fold are significantly different than that of the unfolded rock and its VTI-HTI-TTI orientations (Fig. 3B). The unfolded rock is 18.7% TI with no orthorhombic or lower contributions. Regardless of the averaging method, the fold has developed a larger component of orthorhombic symmetry (4.7%–6.1%) and a decreased amount of TI symmetry (6.3%–6.8%). This fold of unit wavelength has horizontal fast-axis orthorhombic symmetry, where one axis is much faster than the other two that are similar to each other, even though the internal rock is locally slow-axis TI. The sinusoidal tilting of the foliation in the x-z plane tends to moderate the averaged velocities in this cross-sectional plane, whereas the velocities in the y (fold hinge) direction are unchanged, producing the fast-axis orthorhombic symmetry. Second, the fold contains a lower amount of total anisotropy compared to the rock fabric that fills it. This is expressed in the anisotropic components of the symmetry decomposition. The rock fabric has net 18.7% anisotropic components. The fold has net 11.5%–12.5% anisotropy across the averaging methods. The structural folding has muted the anisotropy that would be measured from a local sample of the fabric. Third, comparison of tensor averaging methods reveals only second-order differences among the averaged results. These differences are less significant compared to the first-order change produced by the structural folding (SGA) of the original anisotropic rock.

We calculate the effective media for a series of folds of increasing steepness as expressed by fold limb angle. The numerical values of these effective medium tensors are provided in the Supplemental Material¹. We show the SGA velocities of this series in Figure 4. At a limb angle of θ₀ = 0° the medium is essentially the same as the original unfolded VTI rock (Fig. 3). The progression of fold limb angle from horizontal to vertical (Fig. 4C) correlates with the horizontal bulk schistosity becoming turned on-end. However, the fold elastic tensor does not simply rotate from horizontal to vertical slow-axis symmetry but changes from vertical slow-axis (θ₀ = 0°), weakly orthorhombic with vertical slow-axis as expressed in the 2.1%–4.8% orthorhombic component of symmetry decomposition (θ₀ = 30–45°), horizontal fast-axis TI (θ₀ = 60°; 0% orthorhombic component), orthorhombic with horizontal slow-axis (θ₀ = 75°; 4.2% orthorhombic component), to strong horizontal slow-axis (θ₀ = 90°). The overall amounts of anisotropy for the folds are less than that of the representative fabric that fills them. These systematic changes are well illustrated in Figure 4 and in an animation of velocity spheres for all fold limb angles between 0 and 90° (see Supplemental Material [^{footnote 1}]). Propagation directions of strongest shear wave splitting change as fold limb angle increases. As can be seen in Figure 4, seismic anisotropy signatures and azimuthal patterns differ significantly based on fold severity.

Elastic tensors and velocity spheres describe the bulk characteristics of a structure that is filled with anisotropic material, but seismic waves that travel through the structure offer a more direct demonstration of anisotropic effects. We illustrate the seismic response of the above fold cases using synthetic anisotropic wave propagation in an earth model that contains a horizontal layer internally filled with folds. In examples presented here, we compare azimuthal variation in seismograms through the folds and then the effects of fold wavelength (differing fold limb angles). We also show synthetic seismograms in unfolded VTI and HTI layers for calibration.

The earth model has map dimensions of 60 × 60 km plus 30 vertical km, discretized at 100 m grid spacing (Fig. 5). The model contains two layers, of which the shallower layer (0–10 km depth) contains anisotropic material that will be folded at different fold wavelengths. The lower background layer (10–30 km) has seismic velocities defined by a 1D isotropic P-wave gradient from 5600 to 6800 m/s with S-wave velocities defined using a Poisson’s ratio of 0.25 and density of 2670 kg/m³. We used a 3-D finite difference elastic wave propagation code that uses full elastic tensors and allows for full heterogeneity in composition and anisotropy symmetry/orientation (e.g., Okaya and McEvilly, 2003). A point source was positioned at the center in map view at a depth of 5 km (Fig. 5A). The source used a 3 s to 2 Hz narrow band waveform whose wavelength is ∼2.5 km at a propagation velocity of 6 km/s (Fig. 5D). This source was propagated first as an explosive P-wave for 20 s at 0.005 s sampling in the synthetic calculations, then as a radial SV source. For each source, three-component seismograms were collected every one arc-km (1°) along a circle array of radius 20.0 km centered over the source point (Fig. 5A). This source-to-array geometry produces oblique raypaths of nominal 75° incidence angle. The propagation code used an absorbing boundary condition at the top surface in order to suppress P-SV conversions and surface waves. Finally, we note that this synthetic modeling example was very idealized and resulted in exaggerated amounts of seismic anisotropy signals, such as arrival times and shear wave splitting delay times.

For calibration, the first earth model scenario has no folding within the anisotropic layer. The layer is filled with the representative fabric uniformally oriented with vertical symmetry axis so that the fabric is only horizontal. Figure 3B-VTI thus represents the local rock at each individual model node as well as the bulk properties of the entire layer; the layer has radial symmetry. The synthetic seismograms collected along the circle array are shown in Figures 6A and 6B. The three-component x-y-z seismograms are rotated into L-Q-T, a raypath coordinate frame (Plesinger et al., 1986) that maps compressional energy (P-wave) into the ray-parallel direction (L), radial shear wave energy (SV) into the sagittal plane (Q), and the second shear wave (SH) into the transverse direction (T). These phases appear in the synthetic seismograms as functions of source-to-receiver azimuth. For this VTI model, each phase has constant arrival time in agreement with the layer’s radial symmetry. The transverse SH (T) is the earlier arriving shear wave as was predicted in Figures 3A and 3B. The highly exaggerated split time delay of the slower-traveling SV (Q) phase is caused by the pure horizontal orientation of the anisotropic fabric.

In contrast, we made the opposite end member earth model that has the anisotropic layer filled with vertical fabric schist striking parallel to the Y-axis direction (Fig. 3B-HTI), representing azimuthal anisotropy. The synthetic seismograms produced by the P and SV sources are shown in Figures 6I and 6J. The P-wave arrival (L) is earlier in the 0° and 180° azimuthal directions in the circle array, forming a 2-cycle arrival pattern. Similarly, Q (qS1 in Fig. 3A) also shows a 2-cycle pattern with the same fast direction.

Figure 4C illustrates the geometry of folds with increasingly steep fold limb angle (increasingly large amplitude-to-wavelength ratio). We modeled three fold scenarios of θ₀ = 30°, 60°, and 75°, representing lateral wavelengths of 10.883, 3.628, and 1.684 km, respectively (Figs. 5A–5C). Fold amplitude was set at 1 km. Figures 6C–6H show the circle-array seismograms for the sequence of these three-fold scenarios. For θ₀ = 30°, the P-wave phase (L) shows little arrival time variation as a function of azimuth (Fig. 6C). As θ₀ increases across 60° to 75°, the P-wave phase picks up more arrival time variation (Figs. 6E and 6G). The variation becomes maximal when the fold limb angle is vertical and the anisotropy fabric possesses full azimuthal seismic anisotropy (Fig. 6I).

The length of the S-wave source is short with respect to the wavelength of the 30° folds (Fig. 6D) and responds to individual fold limbs, crests, and troughs. This fold model results in complex Q and T shear-wave arrival patterns. As θ₀ increases to 60°–75° (Figs. 6F and 6H), the fold wavelengths approach and become shorter than the source pulse length. The shear-wave arrivals remain complex. Scattering increases due to the compactness of the fold crests and troughs. From steep (θ₀ = 60°) to vertical folds, the fast shear wave switches from T to Q (SH to SV). The δt null directions change azimuthal directions.

The anisotropy behavior in these seismograms is in agreement with Christoffel solutions as shown in the effective medium velocity spheres (Fig. 4). The incidence angle obliqueness of the modeled source-to-circle array paths translates to a wave direction exiting the spheres at 15° north of the equator in all azimuthal directions. Careful examination of the fast-S particle motion bars along this “latitude” in Figure 4 help identify the Q-T and SV/SH wave behavior. For example, for θ₀ = 75° the particle motion bars indicate fast SV motion in 0° and 180° azimuthal directions and fast SH in narrow azimuthal windows centered at 90° and 270°. As a result, the first arriving S-wave does not have simple 2θ azimuthal behavior; in the synthetic seismograms the Q and T phases (that is, SV and SH) alternate as to which arrives first (Fig. 6H).

Because the production of seismic anisotropy is related to directional dependence of seismic waves with respect to rock properties, structural geometries can be particularly important in the analysis of seismic data. Mathematical functions have been suggested to describe idealized structural geometries. These functions can approximate shapes and provide tilt information for defining anisotropy-filled structures within discretized digital volumes. Figure 7 shows representative examples of structural geometries and descriptive equations. Functions that describe folds can be grouped into non-periodic functions such as polynomials, power functions and ellipse equations (Hudleston, 1973; De Paor, 1996; Bastida et al., 1999; Aller et al., 2004; Bastida et al., 2005) and periodic functions such as trigonometric functions represented as Fourier series and Bessel functions (Currie et al., 1962; Stabler, 1968; Hudleston, 1973; Bastida et al., 1999; Jeng et al., 2002; Bastida et al., 2005) (Fig. 7). The power function z = [1 - (1–4x/L)ⁿ]^1/a can describe folds such as chevron (n = 1, a = 1), parabolic (n = 2, a = 1), ellipsoidal (n = 2, a = 2) and box (n > 2, a = 1).

Other classes of structure that may commonly be encountered are doubly plunging folds, gneiss domes, salt domes and metamorphic core complexes (e.g., Whitney et al., 2004). Each of these can be treated by functions of the form z = f(x,y), for example Gaussian functions or 3-D refolding of the sinusoidal or power functions (Fig. 7). For azimuthally symmetrical shapes exhibited by many dome and basin structures, surfaces of revolution can be generated by rotating the two-dimensional curves about an axis of symmetry. Asymmetry, irregularity, or higher-order spatial patterns can be superimposed onto these functions. Scaling of functions can describe internal fabrics when concordant. Separate functions can be used when the trends of internal fabrics are different from the bounding structure. Spatial and volumetric analyses of fabric tilt patterns are possible because of the use of differentiable functions.

For each structure illustrated in Figure 7, we calculated its SGA effective medium for a unit shape based on a characteristic length or height, assuming the structure is filled with our representative fabric whose foliations conform to the shape. The corresponding anisotropic seismic velocities are obtained from the effective medium tensors and are also shown in Figure 7. The shear zone and monocline shown in Figure 7 have similar geometrical characteristics—a kink within an overall direction but with different orientations. The fabric in the shear zone is near vertical with a subhorizontal kink related to severity of shearing (Fig. 7A). The net fabric trend in the shear zone produces steeply tilted TTI symmetry with a trace of orthorhombic symmetry (1.3%). In contrast, the monocline is subhorizontal with a gentle-to-steep kink, and its bulk fabric-filled anisotropy is gently tilted TTI symmetry (Fig. 7B). The possible orthorhombic-symmetries of sinusoidal folds (Fig. 7C) were described in the previous sections. A plunging fold is elongate based on its lateral aspect ratio. Its bulk anisotropy exhibits orthorhombic symmetry (2.5%) with its fast direction parallel to the major fold axis and the slowest direction normal to the dome top (Fig. 7D). The largest δt is produced by horizontally propagating S-waves across the fold axis. S-waves propagating upward produce small amounts of splitting with fast direction parallel to the fold axes. In contrast to plunging folds, a dome of fabric has vertical fast-axis TI symmetry (Fig. 7E). Its uniaxial shape is expressed in its 6.2% TI and 0% orthorhombic symmetry components. Shear-wave splitting δt is largest for horizontal propagation directions and with vertical fast-directions. Null splitting is produced in the vertical propagation direction. These schematic structures in Figure 7 exhibit a wide range of SGA symmetries even though they are all filled with the same slow-axis TI rock. The effective medium tensors for these structures are provided in Supplemental Material (^{footnote 1}).

In the above structural examples, we illustrated SGA that is filled with rock that itself is already anisotropic. This question then arises: how much of the resulting seismic anisotropy is contributed by the geometry and how much is due to the local rock itself? If the rock anisotropy is changed within the structure (e.g., change in symmetry or percent intensity), does the bulk anisotropy change significantly? Different variations can be explicitly calculated using one of the elastic tensor volume averaging or homogenization methods. However, if we assume that the rock filling a structure is already anisotropic and can be treated as uniform throughout, then it is possible to separate in the volume averaging the structural orientation from this local rock. This assumption enables an examination of the changes in anisotropic velocities when one is changed while the other is held fixed.

If the anisotropic rocks within a structure are uniform or the seismic scale of interest allows for averaging of smaller scale heterogeneity (the long-wavelength case, L << λ, described in “Scales of Waves and Structures”), the various rocks within the structure can be replaced with a single representative rock that fills the structure. So far in our study, we have been using a representative fabric from a real rock sample, but this generalized rock could be defined in several ways: (a) be assigned from an individual rock sample as we have done, (b) combine petrophysics lab measurements or thin-section-based calculations of properties from several outcrops sites, (c) be a generalized rock, or (d) already be averaged from some other set of rock information. We denote this representative rock’s elastic stiffness as the tensor C^rep.

Because the volume-averaged tensor C* represents an effective medium that contains a structure, we name this to be a structural effective medium, EM_s. This formulation in Equation 2D allows the volume of bulk material within a geological terrane or feature to be separable into two components—a generalized type of rock that has its own elastic stiffness and the overall geometrical shape that systematically reorients the rock with respect to a geographical reference frame. This separation of volume averaged orientations and a constant but anisotropic rock tensor was similarly applied by Song and Jordan (2017) in order to study the effective media of mantle shape preferred orientations anisotropy.

This expression in Equation 2D has several important points. (1) Propagating seismic waves will respond to the structural effective medium (EM_s), which will have its own anisotropic seismic velocities. This EM_s is the stiffness tensor used directly in the elastic wave equation. (2) The anisotropic symmetry of the EM_s might not be the same as that of the filling representative rock C^rep. (3) Various geological structures have different SGO, and each will be a function of size or geometrical parameters such as spatial wavelength, amplitude, curvature, scaling, and self-similarity. Because of this, the EM_s of different structures will differ even if filled with the same C^rep rock. (4) For the same structure, a change in the rock type C^rep will affect the EM_s and subsequent seismic wave behavior. These changes pertain to the type of local rock anisotropy: elliptical or weak versus full anisotropy; TI or lower-order symmetry; and fast or slow-axis. Different patterns of azimuthal variation in seismic velocities will result. (5) In practice, the SGO may be efficiently evaluated by discretizing an earth volume into equal-volume voxels and averaging the local rotation operators. However, in selected cases the SGO can also be obtained analytically for simple geometries that are defined using elementary functions. (6) Because the SGO and C^rep are separable under the condition that the representative rock is the same throughout the volume of interest, their macroscale and mesoscale contributions to a bulk material anisotropy (EM_s) can be independently examined.

We note that the derivation of EM_s in Equation 2 is based on Voigt tensor averaging of stiffnesses. A similar derivation is possible using Reuss compliance averaging or Reuss-Voigt-Hill combinations. Because the Voigt and Reuss methods represent upper and lower bounds of averaging results, the Voigt-based and Reuss-based EM_s are also upper and lower bounds, respectively (Bunge, 1985; Vel et al., 2016). As illustrated in Figure 3C, anisotropic seismic velocities derived from the Reuss-based EM_s will be slightly slower than those derived from the Voigt-based EM_s.

We also note that the above points are valid for an anisotropic C^rep. When C^rep is isotropic, the SGO will be applied to this isotropic rock and the resulting EM_s will also be isotropic.

While effective media can be calculated by averaging tensors using 6 × 6 Voigt notation, the commutative law does not hold for this notation (Thomsen, 2002) and the rotation matrices M and M^T cannot be separated from C^rep in Equation 1B. As a result, SGO cannot be calculated using Voigt notation.

The SGO operator can be created via the summation of subvolume rotation matrices ROT_{(x,y,z) 21 × 21}. This 21 × 21 SGO rotation operator is large but not as unwieldy as the full (81 × 81) rotation matrix in ijkl notation. Numerical implementation is straightforward. We provide a table of all elements in the rotation matrix ROT_{(x,y,z) 21 × 21} as defined using directional cosines in Supplemental Material (^{footnote 1}). We also provide a set of software subroutines that can be used to calculate this rotation of a stiffness vector.

We demonstrate the equivalence of a structural effective medium to its equivalent but actually folded fabric. In Figures 6G and 6H we showed the synthetic seismograms produced by a layer of short wavelength folded fabric (λ = 1.684 km, θ₀ = 75°). Here we replaced the actual tilted nodes in the modeling earth volume with a uniform layer of this fold’s equivalent EM_s. The EM_s was obtained by applying Equation 5 to our original representative fabric tensor (Fig. 3A;  Appendix 1). The resulting nonzero EM_s stiffness coefficients for θ₀ = 75° are: EM_s11 = 95.842, EM_s12 = 44.934, EM_s13 = 42.213, EM_s22 = 132.670, EM_s23 = 42.213, EM_s33 = 118.940, EM_s44 = 41.363, EM_s55 = 29.628, EM_s66 = 32.692, in GPa.

The resulting synthetic seismograms through this EM_s layer are shown in Figure 8, next to the seismograms for the tilted tensor model (repeated from Figures 6G and 6H). The seismograms produced through the EM_s model in Figure 8B are simpler in that they do not contain scattered energy due to the actual fold crests and troughs of their counterparts in Figure 8A. The L-Q-T phases produced by the EM_s layer are identical to those from the tilted tensor version. The azimuthal anisotropy in both results is also similar.

Up to this point in our study we have used a representative fabric in order to illustrate folding effects. If the unfolded rock used in Figure 4 were different, its rock elastic tensor and velocity spheres would be different, but the same sequence of SGO rotations would be applied, with the result that similar EM_s tensor modification would occur during increased fold intensity. In Set 1 of Figure 9, fold effective media were recalculated after replacing the representative schist with a slate (“TAIELL”) and a gneiss (“NPELL” in  Appendix 1). We used the elliptical TI symmetry version of these two rock types. In both cases, folds filled with these rocks range between TI to orthorhombic symmetry with horizontal to vertical fast axis symmetry. Propagation zones of faster VP (Fig. 9A) and shear wave splitting directions (Fig. 9B) are similar to those of the representative schist fabric (Fig. 4 and “OROELL” in Fig. 9). For these three rock types, the symmetry decomposition values reveal the effects of folding. For the three types, folding produces the largest amounts of orthorhombic symmetries at θ₀ = 45° and 75°. Nearly pure TI symmetry is restored at θ₀ = 60° but with a symmetry axis that is horizontal, rotated 90° from the unfolded fabric. In all fold cases, the amounts of anisotropy are muted compared to the original fabrics. These changes are related to the geometry of cylindrical folds, not due to the rock anisotropy. The tensor values for these EM_s folds are available in the Supplemental Material (^{footnote 1}).

The role of geological structure on seismic anisotropy can be profound. While the local rock may exhibit fast or slow-axis behavior, we have shown that structure can map a slow-axis rock into a fast-axis regional material and even change the order of symmetry of the structure’s bulk EM_s. The effects of structural geometry are primarily at the macroscale whereas the specifics of rock stiffness are at the micro-to-mesoscale. Despite these differing scales, both the structure and the representative rock combine to produce the EM_s’s directional anisotropy, and the role of each must be studied.

The sinusoid folds presented above illustrate the effect that geometrical structures can impart on rock material anisotropy. An unfolded anisotropic layer possessing slow-axis transverse isotropic symmetry can be converted into moderated fast-axis orthorhombic symmetry with up to 90° reorientation of symmetry axis (e.g., Fig. 4 and the animation in the Supplemental Material [^{footnote 1}]). Inclusion of third dimensional structure onto these folds, such as doubly plunging synform/antiforms, doming, or other modification, has the potential to produce an even more varied bulk effective medium (Fig. 7). However, as the degree of 3-D structural complexity increases, the structural orientations can become so heterogeneous that the resulting effective medium becomes weakly anisotropic to isotropic (e.g., Naus-Thijssen et al., 2011a). Analytical or numerical construction of a structural geometry operator (SGO) facilitates the examination of the separate contributions of rock and geometry anisotropy.

In previous examples of structures filled with fabric, we have used rock material that has elliptical transverse isotropy. Real crustal rocks are more likely to be non-elliptical (e.g., full TI). This more complete form of TI anisotropy produces seismic wavefronts that in 3-D are not purely ellipsoidal. This results in directional variation in seismic velocities that are also less ellipsoidal with a component of 4θ cyclicity (Backus, 1965). The complete definition of transverse isotropic or hexagonal symmetry requires a fifth independent velocity that is commonly a diagonal P-wave velocity (Auld, 1973; Okaya and Christensen, 2002;  Appendix 1).

Figures 3 and 10 compare the differences between elliptical and full TI using the same representative (schist) fabric. A diagonal VP was made in the original laboratory measurements and is slower than the average of the two axial VP (N. Christensen, personal commun., 2017). This affects the angular variation of qVP and qVS₂; qVP remains slower over a wider range of propagation angles. The two qVS wave velocities (Fig. 10A) are quite different from the elliptical case (Fig. 3A). At ∼50°, these velocities “crossover” as to which is faster. A smaller angle range (∼65°–90°) exists with appreciable shear wave splitting (vertical blue lines). The spheres in Figure 10B show the smaller range of propagation angles for S-splitting in this case. The shear wave crossover produces a “halo” of null δt at mid-latitude incidence angles and marks the propagation directions of a major switch in which shear wave polarization is the “fast-S” arrival. Symmetry decomposition shows that the degree of TI symmetry is stronger than the elliptical TI version even though the four axial VP and VS velocities are the same. This non-elliptical anisotropy when mapped into a geological structure will produce a more complex bulk seismic anisotropy behavior.

The full, non-elliptical forms of our three rock types (schist, slate, and gneiss;  Appendix 1) were mapped into fold structures. The effective media of these folds are illustrated in Set 2 of Figures 9A and 9B. The P-wave anisotropy is vertical slow-axis for the unfolded forms (θ₀ = 0°) of all three rock types (OROFULL, TAIFULL, and NPFULL, respectively). Increased folding converts all three into weak horizontal fast-axis orthorhombic symmetry at 45° and 75° fold limb angle, as expressed in their symmetry decomposition values. These turn back into non-elliptical transverse isotropy but with horizontal slow-axis at 60° and when the layering becomes vertical. The S-wave splitting exhibits significant directional variation (Fig. 9B, Set 2). OROFULL and TAIFULL originally have null δt halos at mid-latitude incidence angles. These halos when mapped into the folds weaken the overall values of shear wave anisotropy compared to their elliptical anisotropy counterparts (Figs. 4 and 9B, Set 1). Sample NPFULL does not have a null δt halo but instead has a mid-latitude halo of maximum δt (Okaya and Christensen, 2002). This maximum δt halo when mapped into folds increases the overall value of shear wave splitting and degree of fast S-wave particle polarization. In all of these rock cases, structural folding modifies in a similar manner the rock tensors into their respective bulk fold EM_s.

In all prior sections above, fold examples have been filled with material anisotropy possessing transverse isotropy. Two of our rock samples, the schist and slate, were originally measured with orthorhombic symmetry (Set 3 in  Appendix 1). We recalculated the fold EM_s using orthorhombic symmetry; the resulting fold anisotropic velocities are shown in Set 3 of Figure 9. Orthorhombic behavior of VP is clearly visible for horizontal fabric of both samples. This results in weak horizontal fast axis TI symmetry for folding of 60° limb angle before becoming strong orthorhombic vertical fabric. The directional patterns of δt and fast direction are more complicated in comparison to their elliptical TI counterparts (Set 3 versus Set 1 in Fig. 9B). However, the change in anisotropy symmetry as related to increased folding (the SGO operator) is similar. Figure 9 illustrates that folding modifies the effective medium symmetries in similar ways regardless of the unfolded rock type or its original level of symmetry.

Many studies of crustal seismic anisotropy use analysis methods that assume that the anisotropic crust contains slow-axis elliptical transverse isotropy. As we demonstrated here with simple sinusoid folds and other schematic structures, this might not be a valid assumption if structure produced by tectonic deformation is present. When geological structures reorient rock anisotropy, seismic measurements of anisotropy may not necessarily fit simple 2θ-like patterns. These seismic data may have what appears to be significant scatter but might be better fit by considering lower symmetry and/or a switch in fast/slow axis behavior. The breadth of the geological patterns relative to the seismic wavelength and location of seismic stations needs to be taken into account—whether path integration or effective medium effects are influencing the seismic measurements.

In addition, the type of seismic phase being used and its associated wave propagation directionality (raypath back-azimuth/inclination) are important. Crustal seismic anisotropy studies typically use a single or limited number of phase types, ranging from P or S direct waves, Ps conversions (receiver functions), teleseismic SKS shear waves, and station-to-station correlated ambient noise for surface waves. These seismic phases propagate in fairly restricted directions in the context of the velocity spheres in our study. If the structural geometric anisotropy is quite varied in direction, an individual type of seismic wave reveals a limited amount of information about the full anisotropy provided by the structure. Future anisotropy studies and next-generation seismic field experiments may need to combine information using several wave types in order to piece together the full anisotropy beyond transverse isotropy assumptions.

We summarize this study with these points:

• Structures formed by tectonic deformation can contain a 3-D distribution of local anisotropic material producing structural geometric anisotropy (SGA). The net anisotropy of the structure can be different from the local rock’s anisotropy. Whereas we typically approximate rocks rich in aligned phyllosilicates as having slow-axis transverse isotropy, structure can alter this symmetry into regional media that can be fast-axis and/or orthorhombic or lower symmetry. Thus the origin of crustal seismic anisotropy signals may need to factor in larger scale SGA in addition to material crystallographic preferred orientations or shape preferred orientations.

• Seismic anisotropy produced by geological structures is sensitive to the length scale of the seismic wave relative to the size of the structure. A relatively short wavelength will respond to individual components of a relatively large structure (path integration) whereas a relatively long wavelength might only respond to the bulk average of all components of a relatively small structure (effective medium).

• For calculations of effective media, when larger volumes can be considered to contain the same rock elastic tensor (representative C^rep) tensor, the commutative law of tensor multiplication allows C^rep to be placed outside of the summation within the tensor averaging. When this operation is possible, the set of tilt operators can be volumetrically averaged by themselves. The average of the tilt operators becomes a bulk tilt operator (structural geometry operator or SGO) that maps the local rock C^rep into the associated structure’s effective medium (EM_s).

• An SGO is strictly related to the structural geometry and has no inherent information about the rock elasticity to be mapped into the structure. The modification to seismic anisotropy due to the structure can be examined independent of what rock fills the structure (that is, independent of the rock’s anisotropy intensity or symmetry).

• An SGO can sometimes be determined analytically from curve functions; otherwise, it can be derived numerically for more complex structures. Analytical calculations using the Voigt and Reuss approximations yield results close to those produced by recently developed physics-based tensor averaging methods.

• For cylindrical folds (not folded in the third dimension), their bulk volume anisotropy can be muted orthorhombic symmetry when filled with locally slow-axis transverse isotropic rocks, depending on fold steepness. Using a more complex local starting rock may produce an even lower ordered symmetry effective medium.

• Analysis of crustal seismic anisotropy data using the assumption of (slow-axis) elliptical transverse isotropy may produce overly simplistic interpretations when regional or pervasive fabric structure is present.

Seismic velocities for three rock samples used in this study are listed in Table A1. These velocities were used to construct elastic stiffness tensors of different levels of symmetry (Table A2). These tensors were used as representative tensors in the main text that were modified by geometry in order to produce structural effective media (e.g., Figs. 4, 7, and 9). These tensors were used within Christoffel equations to obtain phase velocities in all propagation directions that are displayed as spheres in the figures.

We obtained rock velocities from N. Christensen (personal commun., 2017) and Meltzer and Christensen (2001), who carried out petrophysical laboratory measurements. Two rock samples, the Orocopia schist (southeast California) and Taiwan Central Ranges slate, were measured for orthorhombic symmetry based on strong fabric and lineations (N. Christensen, personal commun., 2017). The Nanga Parbat gneiss was measured for full hexagonal symmetry but not for orthorhombic symmetry (Meltzer and Christensen, 2001). We used their values that were measured at 200 MPa confining pressure, representing a depth of several km where cracks are predominantly closed.

A full set of independent measurements defining an order of symmetry can be converted to higher order by averaging measurement directions. In this study we calculated the full transverse isotropic (TI) equivalent of the orthorhombic Orocopia schist and Central Range slate samples by averaging selected velocities using these procedures. (1) We identified a sample’s orthorhombic axes that defined its fabric plane (a-b in Table A1) and its normal (c), with a oriented parallel to lineation if present. (2) We assigned the c axis to become the TI symmetry axis (0° in Table A1). The a-b fabric plane became the perpendicular plane (90°) to the TI symmetry axis. For P-wave velocities (VP), orthorhombic VPc was assigned to TI’s VP_0°, and the average of VPa and VPb became TI’s VP_90°. For S-wave velocities, VSab was assigned to TI’s VS_90° and the average of VSac and VSbc became TI’s VS_0°. The fifth TI measurement, VP_45°, was obtained by averaging the orthorhombic VP45ac and VP45bc. These five TI velocities define full hexagonal/TI symmetry that might not have simple two-cycle variation with 360° direction (e.g., Auld, 1973; Okaya and Christensen, 2002).

We then converted the full TI symmetry velocities including the Nanga Parbat gneiss sample using this procedure. (1) All TI axial VP and VS velocities were preserved. (2) The fifth measurement, VP_45°, was replaced with the average of the two new axial VP velocities (VP_0°, VP_90°). This reassignment of the diagonal VP to become just the average of the axial VP velocities has the effect of simplifying this anisotropy into elliptical TI symmetry. In this situation, a point source will produce P- and S-waves that propagate outward with elliptical shaped wavefronts. In contrast, within a full TI symmetry medium the wavefronts will deviate from elliptical wavefront in all non-axial directions (Auld, 1973).

The elastic stiffness tensors use the Voigt notation C_ij where i,j = 1–6. Here, subscript 1 is aligned with symmetry a axis, 2 is aligned with b axis, and 3 is aligned with c axis. These tensors are decomposed into their symmetry components using the method of Browaeys and Chevrot (2004). This method identifies the percentage amounts of symmetry classes that sum into the tensor. In Table A2 we show the isotropic component, then list a three-value sequence of TI, orthorhombic, and all lower-order symmetries. This latter sequence represents the amount of anisotropy within the tensor.

For Figure 9 in main text and Figure A1 here, we group elliptical TI velocities as rock type Set 1, full TI as Set 2, and orthorhombic velocities as Set 3. Figure A1 shows a seismic velocity plot for each rock/symmetry combination in Table 1 for 360° propagation direction.

For elliptical TI symmetry (Set 1 “ELL” samples), the P-wave and qVS₁ (transverse motion or SH) velocities are elliptical over 360° directions (Fig. A1). The qVS₂ (sagittal equivalent to SV) is circular (constant in value). Shear wave splitting over 360° directions is indicated by the difference in the two shear wave curves.

The most noticeable feature of full TI symmetry (Set 2 “FULL” samples, Fig. A1) is that qVS₂ can either cross over with the other S-wave curve or may become slower in the diagonal directions. In either case, shear wave splitting is more varied over 360°. For the Set 3 orthorhombic versions (ORTH) of the schist and Taiwan slate, velocities in the three axial planes (a-b, a-c, and b-c) are shown in the bottom two rows of Figure A1. The 360° velocity curves in Figure A1 are related to the velocity spheres shown in Figures 3, 4, 9, and 10.

We also include in online Supplemental Material (^{footnote 1}) several Fortran code subroutines that define and perform the rotation of a 21 × 1 tensor given a set of directional cosines. These subroutines can be used to calculate the effective media of a geometrical structure filled with heterogeneous elastic tensors or modified into other programming languages. The codes can easily be converted to C, python, or Matlab scripts. The subroutine codes include:

• (a) Transfer between Voigt 6 × 6 elastic tensor and 21 × 1 stiffness vector:

•   subroutine sub_C6x6_to_C21x1(voigt6x6_in,c21x1_out)

•   subroutine sub_C21x1_to_C6x6(c21x1_in,voigt6x6_out)

• (b) Define 21 × 21 rotation matrix elements given 3 × 3 directional cosines.

•   subroutine sub_set_rotmatrix(dircos,rotmatrix_21x21)

• (c) Perform rotation by vector multiplication |R| = |ROTMAT| |C|

•   subroutine sub_perform_c21_rotation(c21_in,rotmat_21x21,r21_out)

The above subroutines can be packaged together in order to rotate a Voigt elastic tensor:

c umbrella subroutine to rotate a 6x6 elastic stiffness tensor using 21-element vectors.

• c (1) Define 21x21 rotation matrix given 3x3 directional cosines.

• c (2) Transfer 6x6 stiffness to 21x1 vector.

• c (3) Rotate stiffness vector, producing rotated R21x1 stiffness.

• c (4) Transfer rotated stiffness vector into 6x6 stiffness.

• c

•  subroutine rotate_stiffness_6x6(c6x6_in,dircos,r6x6_out)

•  implicit none

•  real*8 c6x6_in(6,6), r6x6_out(6,6), dircos(3,3)

•  real*8 c21_in(21),r21_out(21)

•  real*8 rotmatrix_21x21(21,21)

•  call sub_set_rotmatrix(dircos,rotmatrix_21x21)

•  call sub_C6x6_to_C21x1(c6x6_in,c21_in)

•  call sub_perform_c21_rotation(c21_in,rotmatrix_21x21,r21_out)

•  call sub_C21x1_toC6x6(r21_out,r6x6_out)

•  return

•  end

We thank Nik Christensen for petrophysical laboratory measurements of the schist and J.M. Martínez-Martínez for permission to use his photograph of the Betic Mountains. We thank Francis Wu for insightful discussions. We also thank two anonymous reviewers whose comments helped to improve this paper. The Nanga Parbat rock laboratory data used in this paper are available in Meltzer and Christensen (2001; https://doi.org/10.1029/2000GL012262). This study was funded by National Science Foundation EAR-1015599 (Okaya), EAR-1015349 (Vel, Johnson) and EAR-1118786 (Johnson, Vel). Figures were made using the Generic Mapping Tool software package (Wessel and Smith, 1998).