Seismic Anisotropy in Subduction Zones: Evaluating the Role of Chloritoid

Subduction zones are often characterized by the presence of strong trench-parallel seismic anisotropy and large delay times. Hydrous minerals, owing to their large elastic anisotropy and strong lattice preferred orientations (LPOs), are often invoked to explain these observations. However, the elasticity and the LPO of chloritoid, which is one of such hydrous phases relevant in subduction zone settings, are poorly understood. In this study, we measured the LPO of polycrystalline chloritoid in natural rock samples, obtained the LPO-induced seismic anisotropy, and evaluated the thermodynamic stability field of chloritoid in subduction zones. The LPO of chloritoid aggregates displayed a strong alignment of the [001] axes subnormal to the rock foliation, with a girdle distribution of the [100] axes and the (010) poles subparallel to the foliation. New elasticity data of single-crystal chloritoid showed a strong elastic anisotropy of chloritoid with 47% for S-waves (VS) and 22% for P-waves (VP), respectively. The combination of the LPO and the elastic anisotropy of the chloritoid aggregates produced a strong S-wave anisotropy with a maximum AVS of 18% and a P-wave anisotropy with an AVP of 10%. The role of chloritoid LPO in seismic anisotropy was evaluated in natural rock samples and a hypothetical blueschist. Our results indicate that the strong LPO of chloritoid along the subduction interface and in subducting slabs can influence the trench-parallel seismic anisotropy in subduction zones with “cold” geotherms.


INTRODUCTION
Solid-state mantle flow is inferred from seismically observed anisotropy in the upper mantle, in subducting slabs, and in many forearc zones and arc regions. Seismic anisotropy is intricately related to the lattice preferred orientation (LPO) of major mantle minerals (Nicolas and Christensen, 1987;Ben Ismaïl and Mainprice, 1998;Savage, 1999;Karato et al., 2008;Long and Silver, 2008;Jung, 2017). Due to its volumetric abundance in the upper mantle, the LPO of olivine has been extensively studied (Jung and Karato, 2001;Ohuchi et al., 2011;Hansen et al., 2014;Tommasi and Vauchez, 2015;Michibayashi et al., 2016;Skemer and Hansen, 2016;Boneh et al., 2017;Cao et al., 2017;Précigout et al., 2017;Soustelle and Manthilake, 2017). It is well-known that the subducting mantle is variably hydrated, and although olivine is volumetrically the most abundant mineral in subducting slabs, olivine-dominant lithologies alone are often insufficient to provide adequate explanations for the large delay times observed by teleseismic studies (e.g., Katayama et al., 2009;Nagaya et al., 2016). Possible alternative mechanisms suggest that volumetrically less abundant hydrous phases are significantly more elastically anisotropic than major mantle phases (Faccenda et al., 2008;Katayama et al., 2009;Jung, 2011;Mookherjee and Capitani, 2011;Mookherjee and Mainprice, 2014;Kim and Jung, 2015). However, in contrast to olivine, the elasticity and LPOs of volumetrically less abundant hydrous mineral phases are incompletely understood and require further study. The few hydrous minerals like serpentine, amphibole, and chlorite, which have been recently examined, exhibit extreme elastic anisotropy and often a strong LPO (Katayama et al., 2009;Mainprice and Ildefonse, 2009;Jung, 2011Jung, , 2017Puelles et al., 2012;Mookherjee and Mainprice, 2014;Kim and Jung, 2015;Nagaya et al., 2016;Almqvist and Mainprice, 2017;Kang and Jung, 2019;Horn et al., 2020). Chloritoid is such a hydrous phase whose elasticity and LPO remain unknown. Since chloritoid is often found in metamorphosed rocks relevant for subduction zone settings, constraining the elasticity and LPO of chloritoid, it is likely to shed valuable insight on the seismic anisotropy of subduction zones.
In this study, we measured the LPO of chloritoid in schistose rock samples by EBSD analysis and used the results to explore the influence of chloritoid-bearing lithologies on the seismic anisotropy in subduction zone settings. We further reevaluated the pressure-temperature stability field of chloritoid in subduction zones. In addition, the role of chloritoid in the generation of seismic anisotropy of the slab was examined by comparing its seismic signature with that of other previously studied hydrous minerals, which were found to be stable with chloritoid.

Chloritoid-Bearing UHP Rock Samples
The chloritoid-bearing schists studied here come from the Makbal Complex, which is an ultrahigh-pressure (UHP) metamorphic complex located in Kyrgyzstan's Northern Tianshan orogenic belt (Figure 1 and Supplementary Figure 1). The Makbal Complex is a part of a tectonic mélange and comprises an area of ∼12 × 25 km, in which abundant (U)HP talc-bearing schists are exposed as up to 200 m thick layers in the central zone within the UHP area (ca. 10 km in diameter; Meyer et al., 2014;Lee et al., 2020 and references there in). Previous studies indicate that the Makbal complex protolith represents highly altered oceanic crust, and the equilibrium peak mineral assemblage of chloritoid + garnet + talc + phengite + coesite + rutile revealed peak metamorphic conditions of ∼ 2.85 GPa and 580 • C (Meyer et al., 2014). These results imply a subduction up to a depth of 110 km along a cold geotherm, followed by rapid exhumation (Beaumont et al., 2009;Chen et al., 2013;Klemd et al., 2015).

Measurement of the LPO
We determined the foliation by examining the compositional layering of talc and phengite, and the lineation by the elongation of chloritoid and layered minerals along the foliation plane. Thin sections were prepared along the X-Z plane, where X is parallel to the lineation and Z is normal to the foliation, and mechanically and chemically polished using Syton colloidal silica. The LPO of chloritoid was measured using electron-backscattered diffraction (EBSD) with a JEOL JSM-6380 scanning electron microscope (SEM) located in the School of Earth and Environmental Sciences at Seoul National University. The SEM was operated using an accelerating voltage of 20 kV, a vacuum level of chamber of 9.6 × 10 −5 Pa, 1-nA current, 5-µm beam size, and a working distance of 15 mm. All thin sections were carbon coated. To avoid repeated data collection, we indexed the EBSD patterns manually for every single grain. For a better data quality, the LPO of chloritoid was measured with a mean angular deviation of less than 1.0 • and a minimum of 11 diffraction Kikuchi bands for each point measurement. A monoclinic instead of a triclinic space group symmetry of the chloritoid was used for the analysis, because the detected Kikuchi bands closely matched the simulated Kikuchi bands of the monoclinic chloritoid.
The chloritoid LPOs in the Grt-Cld-Tlc schist samples were plotted using the MATLAB toolbox MTEX, version 5.4.0 (Mainprice et al., 2011). For plotting the LPO, the orientation distribution functions of the chloritoids were calculated using a  (Bunge, 1982) and M-index (Skemer et al., 2005).

Elasticity of Chloritoid
We performed first principles simulations based on the density functional theory (Hohenberg and Kohn, 1964;Kohn and Sham, 1965). The chloritoid was investigated with a widely used approximation of the exchange-correlation functionals-the local density approximation (LDA) and the generalized gradient approximation (GGA) with the highly accurate projector augmented wave (PAW) method as implemented in the Vienna ab initio simulation package (VASP) (Kresse and Joubert, 1999). We used the monoclinic (C2/c space group) crystal structure (Comodi et al., 1992) with a Mg 2 Al 4 Si 2 O 10 (OH) 4 stoichiometry and performed a series of convergence tests by varying the energy cut-off and k-points. An observed energy cut-off of E cut −off = 900 eV and a k-point mesh of the 2 × 4 × 1 Monkhorst-Pack grid (Monkhorst and Pack, 1976) with six irreducible k-points were sufficient for describing the energetics of chloritoid. In order to determine the full elastic constant tensor, the lattice parameters were strained by 0.5, 1.0, and 1.5% for one static volume, and we observed that an appropriate ∼1% strain is within the elastic limit (e.g., Chheda et al., 2014). The elastic stiffness constants (C ij s) for monoclinic magnesiochloritoid, calculated using the GGA and LDA methods, are reported in Table 2.

Calculation of the Seismic Velocity and Anisotropy
The seismic anisotropy of single-crystal chloritoid and that of polycrystalline chloritoid were obtained with the MSAT toolkit (Walker and Wookey, 2012), using elastic constants determined in this study. We determined the elastic constants at ambient pressures (P = 0 GPa) using first principles with both the GGA and the LDA. The zero-pressure density predicted using the LDA of ∼ 3.30 gcm −3 is larger than that using the GGA of ∼ 3.10 gcm −3 . Also, the elastic constants predicted using the LDA are stiffer than that of using the GGA ( Table 2). This difference between the GGA and the LDA results is expected and has already been reported for other hydrous mineral assemblages (Mookherjee and Bezacier, 2012;Peng and Mookherjee, 2020). We observed that the aggregate bulk moduli for the GGA are in better agreement with the experimental estimates of the bulk modulus (Grevel et al., 2005). The seismic velocities (V P and V S ) were calculated with the Christoffel equation (e.g., Mainprice et al., 2008). A plot of the elastic anisotropies of the single-crystal chloritoid is presented in Figure 4. The P-wave anisotropy (AV P ) for the mineral aggregates was calculated as: where V Pmax refers to the maximum P-wave velocity (V P ) and V Pmin refers to the minimum P-wave velocity (V P ).
The S-wave anisotropy for the mineral aggregates (AV S ) is defined as using the Voigt-Reuss-Hill averaging scheme, where the V S1 and V S2 represent the fast and slow S-wave velocity, respectively (Mainprice, 1990). The seismic anisotropies of the chloritoid aggregates were calculated with the C ij values determined in this study using the GGA (Figures 4, 5). We calculated the seismic properties of other hydrous minerals in the Grt-Cld-Tlc schists, using the single-crystal elastic constants (C ij ) of chloritoid (this study), chlorite , and for phengitic mica, we used the elastic constants of muscovite (Vaughan and Guggenheim, 1986). For consistency, ambientpressure elastic constants were used for all minerals. To calculate the bulk seismic anisotropy of the Grt-Cld-Tlc schist and a hypothetical blueschist (detailed in section "Effect of Hydrous Minerals in Blueschist-Facies Rock on Seismic Anisotropy"), the elastic constants of the composite minerals were averaged using 2 | Elastic stiffness tensor (C ij in GPa) of single crystal chloritoid using the GGA method (P = 0 GPa, density = 3.10 g/cm 3 ) and the LDA method (P = 0 GPa, density = 3.30 g/cm 3 ).  the Voigt-Reuss-Hill scheme and the modal compositions (see Table 1 and captions of Table 3).

Calculation of the Thermodynamic Stability and Modal Abundance of Chloritoid
We used thermodynamic calculations to obtain the equilibrium pressure-temperature (P-T) stability field for chloritoid, commonly referred to as pseudosections. In addition, we determined the modal abundance (i.e., volume percentage of chloritoid as a function of pressure and temperature) in the different pseudosections calculated in this study. The bulk compositions of three different rock types were considered: (1) Grt-Cld-Tlc schist (this study), (2) average pelite (Shaw, 1956) (Supplementary Figure 3), and (3) aluminous metagabbro (Messiga et al., 1999) (Supplementary Figure 4). For the Grt-Cld-Tlc schist, the effective bulk rock composition of sample #10-16 was used (Meyer et al., 2014). The Theriak-Domino thermodynamic software (de Capitani and Petrakakis, 2010) and the recent version of the Holland and Powell (2011) thermodynamic dataset was used to calculate the P-T pseudosections in the CaO- NCKFMASH) model systems with activity models for chloritoid (White et al., 2000), garnet (White et al., 2007), clinopyroxene (Green et al., 2007), amphibole (Diener and Powell, 2012), talc , white mica (Coggon and Holland, 2002), chlorite , and feldspar (Baldwin et al., 2005). Pure end-member mineral phases include lawsonite, kyanite, quartz/coesite, and water, which was set to be in excess, i.e., all the reactions were considered at water saturated conditions. The P-T range of each pseudosection covers 1.0-3.5 GPa and 450-650 • C.

LPOs of Chloritoid
The chloritoid LPOs are all characterized by a strong alignment of the [001] axes subnormal to the foliation (Figure 3). The strongest chloritoid LPO in sample #15R is characterized by the concentration maxima of the [100] axes subparallel to the foliation and subnormal to the lineation, and by the (010)

Elastic Constant and Seismic Anisotropy of Chloritoid
The elastic anisotropy of the single-crystal chloritoid was calculated by using the elastic constant (C ij ) obtained by first principles simulation at ambient pressure conditions (P = 0 GPa) ( Table 2 and Figure 4). The elastic constants show a similar anisotropy of stiffness components between the directions normal to the c-axis (C 11 , C 22 ) and the c-axis direction (C 33 ), and a strong anisotropy observed among the shear elastic components, i.e., C 66 >> C 44 ∼ C 55 ( Table 2). The seismic velocities and seismic anisotropies of the single-crystal chloritoid for a given propagation direction were calculated from the C ij , the contoured stereograms of the P-wave velocity (V P ), the percentage shear wave anisotropy (AV S ), and the polarization directions of the fastest S-wave (V S1 pol.) for ambient pressures (Figure 4). The pattern of the V P showed high velocities (maximum V P = 9.8 km/s) normal to the c-axis, and low velocities (minimum V P = 7.9 km/s) oblique to the c-axis, displaying a high seismic P-wave anisotropy (AV P = 21.9 and 19.4% for the GGA and LDA, respectively). The AV S normal to the c-axis exhibits maximum AV S = 46.5 and 49.2% for the GGA and LDA, respectively. The AV S parallel to the c-axis exhibits a minimum AV S of 0.6 and 0.1% for the GGA and LDA, respectively. The fast S-wave has a polarization in the basal plane, which is typical for layered silicates (e.g., Morales et al., 2013). However, the elastic constants predicted by the LDA are stiffer than those using the GGA, while the elastic anisotropies predicted by the GGA and LDA are similar (Figure 4 and Table 2). The seismic anisotropies of the P-and S-waves of the polycrystalline chloritoid were calculated (Figure 5) by using the elastic stiffness constant ( Table 2) and the LPOs of chloritoid (Figure 3). The magnitudes of the seismic anisotropies of the polycrystalline chloritoid are large (up to AV P = 10.3%) for the P-waves (sample #15R) and larger (up to AV S = 18.1%) for the S-waves (sample #15R) (Figure 5). The seismic anisotropy of the S-waves is high at angles shallow to the foliation and at a maximum subparallel to the lineation (flow direction) (samples #15R and #10-16), where the polarization direction of the fast S-wave (V S1 ) is nearly perpendicular to the lineation (X-direction).
The fast S-wave polarization direction of the chloritoid strongly depends on the wave propagation direction ( Figure 5). As the ray path proceeds far-off the Z-direction and close to the direction parallel to the foliation, the polarization direction varies from oblique to nearly perpendicular to the rock lineation (X) (Figure 5).
The bulk rock seismic anisotropies of the three Grt-Cld-Tlc schist samples were estimated (Supplementary Figure 5) by averaging the elastic constant of the mineral aggregate, considering the LPOs of the chloritoid and other composite minerals (such as garnet and talc), their individual single-crystal elastic constant, the densities of the individual phases, and the modal volume percentages of the sample ( Table 1). The magnitudes of the whole-rock seismic anisotropies were found to be very large for the P-waves (up to AV P = 28.1%) and large for the S-waves (up to AV S = 15.5%). For all three Grt-Cld-Tlc schist samples, the entire rock AV P (AV P = 14.7-28.1%) is higher than that of the chloritoid aggregates (AV P = 5.3-10.3%). On the other hand, the maximum whole rock AV S of samples #15R and #12-52 (max AV S = 15.5 and 7.0%, respectively) is lower than that of the chloritoid aggregates (max AV S = 18.1 and 9.7%, respectively). FIGURE 4 | Elastic anisotropy of single-crystal monoclinic chloritoid at a pressure of 0 GPa. The P-wave anisotropy was high (AV P = 21.9 and 19.4% for the GGA PAW and LDA methods, respectively) and the shear wave anisotropy was very high (max AV S = 46.5 and 49.2% for the GGA PAW and LDA methods, respectively). The polarization directions of the fast shear wave (V S1 polarized) are shown as black bars. A reference frame for monoclinic symmetry is shown on the left and depicts three orthogonal axes X1, X2, and X3, where X2 is parallel to the crystallographic b-axis and the plane normal to the b-axis is a mirror plane containing X1 and X3 (where X1 corresponds to the reciprocal a*-axis and X3 is parallel to the crystallographic c-axis). V P , P-wave velocity; AV S , shear wave anisotropy; V S1 pol, polarization direction of the fast S-wave; V S1 , fast S-wave velocity; V S2 , slow S-wave velocity.
The maximum whole rock AV S (max AV S = 14.5%) of sample #10-16 is higher than that of the chloritoid aggregates (max AV S = 12.6%). The pattern of the V S1 polarization direction of the bulk rocks is very similar to that of the chloritoid aggregates.

Thermodynamic Stability and Modal Abundance of Chloritoid
The P-T pseudosection of the Grt-Cld-Tlc schist (sample #10-16) calculated in this study shows that chloritoid is thermodynamically stable in the temperature range of 450-600 • C ( Figure 6A). The modal abundance of chloritoid was calculated for the three different bulk compositions ( Figure 6B and Supplementary Figures 3B, 4B). The calculations for the Grt-Cld-Tlc schist revealed that the highest amount of chloritoid was 28 vol.% within the given P-T conditions, which is in good agreement with the observed amount of chloritoid in the samples studied here (20-25 vol.% in samples #10-16 and #12-52) ( Figure 6B). The possible modal amounts of chloritoid in the average pelite (Shaw, 1956) and the Mg-Al metagabbro (Messiga et al., 1999) are up to 20 and 22 vol.%, respectively. The variation of the chloritoid modal volume shows the development of a complex pattern during pressure increase in all three rock types. Along cold subduction geotherms (5-10 • C/km), the maximum chloritoid amount reaches 26 vol.% ( Figure 6B) and is highest above the chlorite-out reactions in the Grt-Cld-Tlc schist (Figure 6). In the Grt-Cld-Tlc schist, chloritoid is stable at <600 • C, and its amount generally decreases with increasing P-T conditions (Figure 6). Other hydrous minerals stable with chloritoid are chlorite, lawsonite, talc, mica, Mg-carpholite, and zoisite. When considering the average pelitic composition, the stability field of chloritoid is constrained to be at >1.5 GPa and <580 • C (Supplementary  Figure 3). Hydrous minerals stable with chloritoid are chlorite, glaucophane, lawsonite, mica, and carpholite. When considering the Mg-Al gabbro composition, chloritoid occurs at relatively high P-T conditions (>2 GPa and >570 • C) (Supplementary Figure 4). Other hydrous minerals stable with chloritoid are chlorite, glaucophane, lawsonite, and zoisite.

LPO Development of Chloritoid
A strong alignment of the [001] axes of the chloritoid was observed in the UHP schist samples from the Makbal Complex. This is consistent with findings of Haerinck et al. (2015), although they used a different analyzing method for the fabric measurement (synchroton X-ray diffraction) on only one chloritoid-bearing sample.
The flattening type of fabric is prevalent in previous studies of layered silicates (e.g., Lloyd et al., 2009;Wenk et al., 2010;Dempsey et al., 2011). The basal glide is a dominant deformation mechanism for most phyllosilicates, allowing the basal (001) planes to occur parallel to the rock foliation or the cleavage planes (Wenk et al., 2010). Although chloritoid FIGURE 5 | Seismic anisotropies of polycrystalline chloritoid in the UHP Makbal schist samples. The elastic anisotropy using GGA method (Figure 4) was applied to the calculation (see "Materials and Methods" section). The P-wave anisotropy is up to AV P = 10.3% (sample #15R), and the S-wave anisotropy is up to AV S = 18.1% (sample #15R). The anisotropy is shown on a stereonet, in which the center of a plot (Z) represents the direction normal to the foliation (X-Y plane), and X corresponds to the lineation. V P , P-wave velocity; AV S , shear wave anisotropy; V S1 pol, polarization direction of the fast S-wave; V S1 , fast S-wave velocity; V S2 , slow S-wave velocity.
is classified as orthosilicate, it commonly exhibits the physical and structural properties of phyllosilicates due to the layered structure, which consists of octahedral sheets bonded by layers of isolated silicate tetrahedra (Klein, 2002;Nesse, 2009). The basal glide could, therefore, be a plausible deformation mechanism of the chloritoid. An alternative explanation is the formation of the chloritoid LPO by dislocation creep, manifested by undulose extinction of the chloritoid in sample #10-16 (Supplementary Figure 2).

Chloritoid Stability and Its Implication for the Seismic Anisotropy in Subduction Zones
A strong trench-parallel seismic anisotropy of S-waves has been observed in many subduction zones such as Ryukyu, Izu-Bonin, Mariana, and Tonga (Long, 2013 and references there in). An anisotropic forearc mantle has been proposed as the major source in terms of the B-type LPO of olivine in the hydrated mantle (Jung and Karato, 2001;Kneller et al., 2005;Jung et al., 2006;Karato et al., 2008;Nagaya et al., 2014), the trench-parallel flow in the sub-slab mantle (Russo and Silver, 1994), the pressureinduced olivine LPO (Jung et al., 2009;Ohuchi et al., 2011;Lee and Jung, 2015), and the strong LPO of hydro-phyllosilicates (Katayama et al., 2009;Jung, 2011Jung, , 2017Nagaya et al., 2016;Ha et al., 2018;Lee et al., 2020). However, the strong seismic anisotropy in some subduction zones cannot be solely due to the LPO of olivine (Wirth and Long, 2012). Recent geodynamic studies using realistic slab parameters showed that pure trenchparallel flow is unlikely to be dominant in the sub-slab mantle (e.g., Alisic et al., 2012). The LPO of anisotropic hydrous minerals in subduction zones may, therefore, be an important manifestation in the trench-parallel seismic anisotropy.
Thermodynamic phase stability indicates that chloritoid is likely to be stable between a depth range of 80 and 120 km in subducting hydrated MORB (Schmidt and Poli, 1998; Figure 7A). The slab geotherms in steeply dipping subduction zones such as the Mariana (slab-dip angle = 57-62 • ), the Izu-Bonin (slab-dip angle = 46-63 • ), the Tonga-Kermadec (slabdip angle = 52-56 • ), and the Nicaragua (slab-dip angle = 62 • ) (Syracuse et al., 2010) settings, which display a strong trenchparallel seismic anisotropy, suggest that chloritoid is likely to be stable along the interface of the overlying mantle wedge and the hydrated oceanic slab at a depth range between 80 and 120 km ( Figure 7B) (Schmidt and Poli, 1998;Syracuse et al., 2010). Recent studies suggest that the subducting oceanic crust and the slab-mantle interface are pervasively hydrated (Abers et al., 2017) and thus might be an important source of the trench-parallel anisotropy (Huang et al., 2011).
To evaluate the effect of chloritoid on the seismic anisotropy observed in subduction zones, the three-dimensional seismic anisotropy pattern of chloritoid aggregates was calculated using the chloritoid LPO of sample #15R (Figure 7C), which showed the strongest AV S (up to ∼18%) among the samples analyzed here. In the present study, the subduction zone setting was simplified by using a 2-D corner (slab-parallel) flow model that resulted from viscous coupling between the downgoing slab and the overlying mantle (van Keken et al., 2002). Thus, the respective slab-dip angle was applied to the chloritoid LPO and the seismic anisotropy pattern by rotation around the axis parallel to the Y-direction, which is normal to the lineation parallel to the downgoing flow direction (Figure 7C). The higher the slab-dip angle, the higher is the strength of the AV S using the chloritoid aggregate, which most likely also causes a longer delay time >0.2 s at a dip-angle θ >45 • . Subduction zones that show strong trench-parallel S-wave anisotropies (with a long delay time of up to 1-2 s) are also characterized by large dipping slab angles of 40-60 • as displayed by the Ryukyu, Mariana, Izu-Bonin, and Tonga subduction zone settings (Long and Silver, 2008). For example, the Tonga and Mariana slabs have steep subducting angles (50-70 • ) and show a strong trench-parallel seismic anisotropy under the forearc and arc regions (Smith et al., 2001;Pozgay et al., 2007). Since the AV S pattern of the polycrystalline chloritoid was controlled by the chloritoid LPO and implied an importance of high-angle slab dip geometry (Figure 7C), the effect of the chloritoid LPO on the seismic anisotropy is thought to be important in cold and high-angle subduction zones.

Effect of Chloritoid LPO on Seismic Anisotropy of the Grt-Cld-Tlc Schist
Garnet is one of the major minerals in the Grt-Cld-Tlc schists, and is seismically nearly isotropic (AV P = 0.5%, max AV S = 1.2% of single garnet crystal) (Babuška et al., 1978). Chloritoid and talc, which are the other major constituent minerals in the Grt-Cld-Tlc schist, are hydrous layered silicates and are elastically anisotropic. Talc in particular has been shown to be an important mineral that influences the strong trench-parallel seismic anisotropy of hydrated slab-mantle interfaces (Lee et al., 2020;Nagaya et al., 2020). The bulk seismic anisotropies of the three Grt-Cld-Tlc schist samples studied here (#15R, #10-16, and #12-52) were determined in an earlier study (Lee et al., 2020), which, however, did not consider the chloritoid LPO due to the lack of elasticity data of single chloritoid crystal at that time. Using the new elasticity data of single-crystal FIGURE 7 | (A) Schematic diagram illustrating a subducting slab with a dipping angle of 60 • . Also shown are the thermodynamic stability fields of hydrous phases in the hydrated peridotitic layer, overlying the oceanic crystal layer and the mantle wedge region (modified after Schmidt and Poli, 1998). Stippled lines are isotherms and arrows indicate flow lines in the mantle wedge. (B) Pressure-temperature stability fields of hydrous phases in hydrated MORB under subduction zone conditions (modified after Schmidt and Poli, 1998). Dashed lines are the P-T paths of the slab surface for the Mariana, Izu-Bonin, Tonga-Kermadec, and Nicaragua settings (Syracuse et al., 2010). Yellow shaded region indicates the stability field of chloritoid based on the pseudosection modeling conducted here ( Figure 6 and Supplementary Figures 3, 4) and also on previous estimates (Schmidt and Poli, 1998). (C) Schematic diagram of a subduction zone illustrates the influence of the chloritoid LPO on the seismic anisotropy in the subduction zone. The applied slab-dip is 60 • (rotation around the axis parallel to the Y-direction; see section "Discussion") and the lineation (X-direction) is assumed to be parallel to the subduction direction. Red bars represent the polarization direction of the fast S-wave. Inset figure shows shear wave anisotropy (AV S ) and the polarization directions of fast S-waves (red and black bars) going through the chloritoid aggregates (sample #15R). X and Z represent the lineation (flow direction) and the direction perpendicular to the foliation (flow plane), respectively. chloritoid and the LPO of chloritoid aggregates, we calculated the complete bulk seismic anisotropy of the Grt-Cld-Tlc schist samples in this study (Supplementary Figure 5). The integration of the chloritoid LPO and its volume fraction revealed that the chloritoid LPOs caused a weakening of the bulk P-wave seismic anisotropies of the Grt-Cld-Tlc schist (AV P = 14.6-28.1%) when compared to those determined in the former "chloritoid-absent" whole rock (AV P = 19.1-31.2%) study of Lee et al. (2020). This is thought to be associated with the significantly lower P-wave seismic anisotropies of the chloritoid aggregates (AV P = 5.3-10.3%), seeing the much higher P-wave seismic anisotropies (AV P = 67.3-72.3%) of the abundant talc aggregates (Lee et al., 2020). On the other hand, the chloritoid LPO slightly increased the maximum AV S of the entire rock samples #15R and #12-52, and decreased that of sample #10-16 (Lee et al., 2020). This observation indicates that the whole rock mineral assemblage is important to determine the role of chloritoid in the whole rock S-wave anisotropy.
We used sample #15R, which displayed the strongest chloritoid LPO (Figure 3), to calculate the effect of the slab-dip angle and the sample geometries on the ray path of a vertically incident S-wave through chloritoid and the entire rock specimen (Figure 8). The sample geometry is defined with respect to the sample foliation (X-Y plane normal to Z) and sample lineation (X). The values approaching 0 • correspond to the trench-parallel V S1 polarization direction, and those approaching 90 • to the trench-normal V S1 polarization direction (see right column of Figure 8). These calculations revealed that under broad sample geometry conditions, the talc LPO produced a longer delay time than the chloritoid LPO. The maximum delay time generated FIGURE 8 | Illustration of the 3-D effect of chloritoid, talc, and the Grt-Cld-Tlc schist on a vertically incident S-wave. LPO data of garnet and talc (Lee et al., 2020) and chloritoid (this study) of sample #15R were used for the calculation. An anisotropic layer with a thickness of 10 km was assumed. The rotation around X and Y from 0 • to 90 • was illustrated to display the effect of the dipping foliation (X-Y plane) where X is parallel to the shear direction and Y is parallel to the strike of the hypothetical trench. The left column shows contour plots displaying the effect of sample orientation on the delay time (dt), with blue (short delay times) to red (long delay times) colors. The right column shows contour plots displaying the variation of the angle between the fast S-wave polarization direction (V S1 pol. direction) and the hypothetical strike of the trench (Y) as the sample is rotated around X and Y. Colors range from blue (trench-parallel fast directions) to red (trench-normal fast directions). In the uppermost left, the schematic sample geometry is illustrated to show how the sample (X-Y plane = foliation, dots, and lines = lineation) is rotated with respect to the incoming S-waves (red line).
by chloritoid is ∼70% of that generated by talc. In the Grt-Cld-Tlc schist, a trench-parallel V S1 polarization direction is achieved when the shear plane dips at an angle >30 • from the horizontal plane (i.e., rotation around the Y direction is greater than 30 • ). When rotating the shear plane around the shear direction at an angle lower than 30 • , the V S1 polarization direction is normal to the trench. The maximum delay time (∼0.4 s) of the Grt-Cld-Tlc schist is generated when the shear plane is dipping < 50 • and its rotation around the shear direction is >60 • (Figure 8), resulting in a trench-normal V S1 polarization direction of the Grt-Cld-Tlc schist. However, When the shear plane dipped at an angle >30 • (i.e., rotation around the Y-direction >30 • ), the variation of the whole rock V S1 polarization direction is similar to that of the talc aggregate when using a shear plane dip angle >30 • (i.e., rotation around the Y-direction > 30 • ). Overall, it was observed that the geometries of the S-wave anisotropy of the Grt-Cld-Tlc schist are significantly more influenced by the talc LPO than by the chloritoid LPO.

Effect of Hydrous Minerals in Blueschist-Facies Rock on Seismic Anisotropy
Considering previous geothermobarometrical studies (see Meyer et al., 2014, and referenced therein) and the pseudosections calculated here, chloritoid is found to be stable under highpressure conditions ranging from the blueschist to the eclogite facies. The here calculated pseudosections for metapelites and metagabbros revealed that a number of hydrous minerals such as glaucophane (Gln), lawsonite (Lws), chlorite (Chl), and phengitic white mica (Ms) are stable with chloritoid under high-pressure conditions (Supplementary Figures 3A, 4A). In order to compare the contributions of these minerals to the seismic anisotropy of a subducted slab, we calculated the delay time and the polarization directions of the fast S-waves, with a vertically incident S-wave ray path traveling though the hydrous minerals, and various sample geometries and slab-dip angles (Figure 9). In addition, we compared these data with the whole-rock seismic anisotropy of a hypothetical blueschist, comprising the respective hydrous mineral assemblages (Figure 9) with a modal composition of Gln:Lws:Chl:Phg:Cld = 30:25:15:15:15 (vol.%). For sample geometry consistency, the foliations of the different rock samples were assumed to be parallel, and, furthermore, we followed the previously published original 2-D plots of the seismic anisotropy of both glaucophane and lawsonite (cf. Cao et al., 2013;Choi et al., in unpublished). The magnitude of the seismic anisotropies and the S-wave velocities are summarized in Table 3.
Our calculations demonstrated that the ray path (a, b, c, d, and e) and the sample geometries, which can generate a long delay time with the trench-parallel S-wave anisotropy, were different for the five hydrous minerals. The chloritoid and phengite LPOs produced their maximum delay time (∼ 0.3 and ∼1.0 s, respectively) with a trench-parallel or -subparallel V S1 polarization direction, a steep dipping shear plane > 60 • (i.e., rotation around Y direction was greater than 60 • ), and a shear plane rotation around the shear direction (i.e., rotation around the X-direction) at an angle lower than 30 • (ray path "a" in Figure 9). The chlorite LPO (cf. sample #15R) displayed its maximum delay time (∼0.4 s) with a trench-normal V S1 polarization direction, a shear plane dip angle < 40 • , and a rotation angle of < 10 • or > 70 • around the shear direction (ray path near "c" and "d, " respectively in Figure 9). The glaucophane LPO generated its maximum delay time (∼0.2 s) with a trench-normal V S1 polarization direction, a shear plane dip-angle below ∼40 • and a rotation angle > 70 • around the shear direction (ray path near "d" in Figure 9). The lawsonite LPO generated its maximum delay time (∼0.4 s) with a steep shear plane dip-angle greater than ∼70 • (ray path from "a" to "b" in Figure 9), or with a shear plane rotation around the shear direction at an angle of >70 • (ray path from "b" to "d" in Figure 9). The latter case corresponds to a highangle fault present in subducting slabs. The lawsonite LPO also produced a trench-parallel anisotropy at a shear plane rotation angle >50 • around the shear direction. Among the different minerals considered here, the chloritoid and phengite LPOs mainly contribute to the trench-parallel polarization of the V S1 with a long delay time (near the ray path "a" in a wide field of blue color in the right column plots in Figure 9) at a large slab dip angle, corresponding to a cold and old subducted slab.
The hypothetical blueschist displayed the maximum delay time of ∼0.2 s, which is the shorter than that of all the composite minerals studied here (see bottom of Figure 9). This is thought to be due to the following two characteristics: (1) the presence of large amounts of glaucophane, which displays the shortest maximum delay time of the minerals studied here, in the blueschist (Figure 9), and (2) the low absolute values of the S-wave velocities of the polycrystalline phengite, lawsonite, and chlorite ( Table 3). The presence of "slow" phengite (i.e., minimum V S1 = 3.2 km/s and minimum V S2 = 2.9 km/s) seems to significantly decrease the maximum S-wave velocity of the blueschist (Table 3), thereby causing a decrease of the maximum S-wave anisotropy and the shorter delay time of the blueschist (Figure 9). Thus, the S-wave velocities of seismically anisotropic minerals such as chlorite and phengite with different velocity ranges may reduce the entire rock's seismic anisotropy. In the present study, chloritoid is the only hydrous mineral (V S1 = 4.7-5.4 km/s, V S2 = 4.5-5.2 km/s) that displayed higher S-wave velocities than glaucophane (V S1 = 4.4-4.9 km/s, V S2 = 4.3-4.5 km/s) ( Table 3).
The delay time pattern of the hypothetical blueschist was mostly similar to that of the chloritoid compared to the other hydrous minerals, i.e., long delay times were generated by the ray paths "a" and "b, " and short delay times by the ray paths "c" and "e, " with the shortest delay time along the ray path "e." The pattern of the V S1 polarization direction of the blueschist was also similar to that of the chloritoid (Figure 9). In order to generate trench-parallel polarization of the blueschist with a measurable delay time of >0.1 s near the ray path of "a, " a shear plane dip-angle of >∼50 • to a horizontal surface (i.e., rotation around the Y direction > 50 • ) was required. Similar conditions were required to generate trench-parallel anisotropy in the LPOs of chloritoid, phengite, and glaucophane (Figure 9). On the other hand, the trench-normal S-wave anisotropy with a measurable delay time of the blueschist is mostly affected by the LPO of the chlorite, phengite, and glaucophane near the ray path of "b." Chloritoid and phengite mostly contribute to the S-wave anisotropy of blueschists and can influence the trench-parallel anisotropy with a long delay time in cold subducting slabs. ). An anisotropic layer thickness of 10 km was assumed for consistency. The reference frame for the rotation and the sample orientation (schematics in the upper left) is the same as in Figure 8. Color ranges describe the same parameters as in Figure 8 for both delay time and fast polarization direction. Symbols marked with "a"-"e" indicate five different ray paths and sample geometries.

CONCLUSION
The elastic anisotropy of chloritoid single-crystals was calculated for the first time and used to examine the seismic anisotropy for chloritoid-bearing rocks in subducting oceanic crust. The elastic stiffness tensor of single-crystal monoclinic chloritoid showed a high elastic anisotropy of chloritoid reflected by AV P = 22% and max AV S = 47%. The LPO of the polycrystalline chloritoid from the UHP schists showed a strong alignment of the [001] axes of chloritoid subnormal to the foliation. The strongest chloritoid LPO showed a concentration or girdle distribution of the [100] axes and the (010) poles subparallel to the foliation. The re-evaluated stability field of the chloritoid revealed the stability of a significant amount of chloritoid in high-pressure and low-temperature subduction zones, particularly in Al-rich metapelites and metagabbros. The LPOs of the polycrystalline chloritoid produced a strong trench-parallel seismic anisotropy. The chloritoid LPOs tend to reduce the bulk seismic anisotropy of the Grt-Cld-Tlc schists due to their lower volume fraction and lower seismic anisotropy compared to other abundant hydrous minerals such as talc. However, depending on the mineral assemblage and slab geometry, chloritoid, in association with glaucophane and phengite, can influence the seismic anisotropy of blueschist-facies rocks in cold subduction zones.

DATA AVAILABILITY STATEMENT
All datasets generated for this study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

AUTHOR CONTRIBUTIONS
JL measured the LPO and calculated the seismic velocity and anisotropy. MM calculated the elastic constants. TK assisted the use of thermodynamic code. HJ guided the project. RK collected the rock samples. JL, MM, and HJ analyzed the data. JL, MM, TK, HJ, and RK wrote the manuscript. All authors contributed to the article and approved the submitted version.