Constraints on Bend‐Faulting and Mantle Hydration at the Marianas Trench From Seismic Anisotropy

Subduction zones are a key link between the surface water cycle and the solid Earth, as the incoming plate carries pore water and hydrous minerals into the subsurface. However, water fluxes from surface to subsurface reservoirs over geologic time are highly uncertain because the volume of water carried in hydrous minerals in the slab mantle is poorly constrained. Estimates of slab mantle hydration based on seismic tomography assume bulk serpentinization, representing an upper bound on water volume. We measure azimuthal seismic anisotropy near the Marianas Trench, use spatial variations in anisotropy to constrain the extent and geometry of bend‐related faulting, and place a lower bound on slab mantle water content for the case where serpentinization is confined within fault zones. The seismic observations can be explained by a minimum of ∼0.85 wt% water in the slab mantle, compared to the upper bound of ∼2 wt% obtained from tomography.

• We measure spatial variations in upper mantle anisotropy that indicate bend-faulting near the Marianas Trench • Hydration localized to bend-faults places a lower bound on the amount of water carried in the subducting slab mantle • Synthetic seismograms compared to the observed anisotropy indicate a minimum of 0.85 wt% water in the slab mantle

Supporting Information:
Supporting Information may be found in the online version of this article.
including seismic velocity and electrical conductivity (e.g., Christensen, 2004;Naif et al., 2015). However, hydration estimates from tomography models typically assume pervasive, bulk serpentinization of mantle material as regularized tomographic inversion is unable to resolve narrow fault zones (e.g., Van Avendonk et al., 2011). Bulk serpentinization is not necessarily the best physical model if water accesses the slab mantle through faults, and hydration estimates on this basis represent an upper bound. The true volume of water may be significantly less if serpentinized materials are localized in fault zones (Miller & Lizarralde, 2016;Obana et al., 2019).
Seismic anisotropy can be used to measure the effects of bend faulting and, by extension, mantle hydration at trenches. Azimuthal anisotropy in the oceanic upper mantle is typically attributed to either crystal preferred orientation (CPO) or fabric of mantle olivine generated by shearing during plate formation, or some form of thin layering such as aligned cracks, faults, or joints (Backus, 1962(Backus, , 1965Hess, 1964;Hudson, 1981;Skemer & Hansen, 2016;Zhang & Karato, 1995). Anisotropy in the incoming plate at a subduction zone is therefore expected to come from a combination of pre-existing olivine fabric and trench-parallel bend faulting, and spatial variations in anisotropy can be used to decompose the two components (e.g., Miller & Lizarralde, 2016;Miller et al., 2021). Olivine fabric formed in a nominally dry upper mantle is predicted to have 2θ periodicity with the anisotropic fast direction oriented along the direction of flow, and in oceanic lithosphere shearing associated with plate formation tends to form a CPO with the fast direction aligned with the paleospreading direction (Hess, 1964;Mainprice, 2015;Skemer & Hansen, 2016). According to analytic solutions, aligned cracks produce an azimuthal signal with strong 4θ periodicity, while the signal from faults or joints is predominantly 2θ with the fast direction oriented along strike of the faulting. If a plate is subducting at an angle with respect to the paleospreading direction from plate formation, the 2θ anisotropy from faults and fabric should be out of phase with each other.
In this study, we use active-source refraction data from the incoming plate at the Marianas Trench to measure azimuthal anisotropy in the upper mantle and investigate how anisotropy varies spatially. The incoming plate at Marianas is relatively old (∼150 Myr, Müller et al., 2008) and thus likely to carry more serpentinized material than a younger, hotter plate (Hacker, 2008). Normal faulting earthquakes on the incoming plate extend to 35 km depth below the seafloor , and both active-source studies (Eimer, 2020) and passive tomography (Cai et al., 2018) show decreases in seismic velocity interpreted as the effects of hydration along faults. The angle between the Marianas Trench and the orientation of the paleo-ridge where the incoming plate formed is ∼50°, indicating that anisotropy from faults and mineral fabric should have distinct orientations. The anisotropy measurement enables us to place bounds on bend-fault geometry, and thereby estimate a lower bound for slab mantle hydration at this subduction zone.

Data and Methods
We measure azimuthal anisotropy in the incoming plate upper mantle near the Marianas Trench using a delay time approach (Morris et al., 1969;Shearer & Orcutt, 1986). We obtain traveltime residuals, or delay times, by comparing picked traveltimes for mantle refractions to predicted times traced in a 3D isotropic reference velocity model. These delay times are inverted to obtain slowness variations as a function of azimuth. We partition the delay times based on source-receiver location and ray-bottoming depth, and invert subsets of the data to measure spatial variations in anisotropy.

The 2012-2013 Marianas Seismic Experiment
Data used in this study are from a 2012-2013 active and passive seismic experiment in the central Marianas Trench and forearc (Wiens, 2014). Four active-source reflection and refraction lines were collected over an array of short-period ocean bottom seismometers (OBS), moored hydrophones, and broadband OBS deployed across the incoming plate and forearc ( Figure 1). Instruments along the active-source lines were spaced ∼20 km apart and recorded shots fired every 500 m by the 36-element, 6,600 in 3 airgun array on the R/V Marcus G. Langseth. This study uses active-source refractions recorded by 38 instruments on the incoming plate. Shots occurring over the forearc and arrival times recorded by forearc instruments are excluded to avoid structural complications.
Previous work in Central Marianas has shown a region with reduced seismic velocity in the incoming plate upper mantle near the trench, interpreted as representing hydrous alteration of the slab mantle (Cai et al., 2018;Eimer, 2020). A low velocity anomaly imaged with Rayleigh wave tomography was interpreted as a 24 ± 5 km thick partially serpentinized slab mantle layer carrying ∼2 wt% water into the trench (Cai et al., 2018). The Rayleigh wave velocities showed significant anisotropy with N-S fast directions at periods sensing the uppermost mantle, interpreted as the effects of shape preferred orientation due to hydrated bending faults. A zone with uppermost mantle Vp as low as 7.3 km/s in the incoming plate bending region imaged with active-source refraction tomography was explained either by ∼20 vol% serpentinization of the mantle immediately below the Moho (corresponding to ∼2 wt% water), or by lesser amounts of serpentine confined to faults or joints (Eimer, 2020). Similar indications of mantle hydration via deep normal faults in the incoming plate have been seen in the Southern Marianas Trench (Zhu et al., 2021). Refraction tomography along the eastern trench-parallel active source line showed Vp ∼ 8.0 km/s immediately below the Moho, consistent with relatively unaltered uppermost mantle outboard of the onset of bending (Feng, 2016).

Data Reduction and Picking
The data are processed to enhance the airgun signals, attenuate noise, and facilitate picking of the mantle refraction Pn. The active-source shot lines were traversed multiple times, providing co-located traces which are stacked to increase the signal-to-noise ratio. The stacked traces are then filtered using a 4-pole, minimum phase Butterworth bandpass filter with corner frequencies of 4 and 15 Hz to enhance the predominantly ∼10 Hz signal from the airgun source. We calculate simple static time corrections for each source location and for moored hydrophones using a water velocity of 1.5 km/s and shift traces to correct for excess traveltime through deep water. The static shifts improve trace-to-trace coherence and clarify apparent velocity for Pn, and are removed after picking. The data are reduced at a velocity of 8 km/s for picking, and linear prediction deconvolution is applied to attenuate reverberations in shallow crustal structure (Henkart, 2003).
Traveltimes for Pn first arrivals are picked on the vertical channel where available; and in differential pressure gauge data for moored hydrophones and for OBS where the vertical channel data are of poor quality. Pn is identified based on increased apparent velocity relative to the crustal refraction Pg and, in places, by the triplication of the Moho reflection Pmp ( Figure S1 in Supporting Information S1). The resulting Pn traveltime picks for 12,888 unique source-receiver pairs provide good spatial coverage over the incoming plate and span all azimuths ( Figure S2 in Supporting Information S1).

Delay Times and Anisotropy Inversion
Delay times are obtained by comparing observed traveltimes to predicted times traced in a 3D isotropic velocity model. This implicitly assumes that our isotropic model is accurate; that crustal anisotropy has a negligible effect on Pn traveltimes due to the near-vertical propagation of Pn through the crust; that the traced mantle raypaths are not significantly different for isotropic versus anisotropic mantle; and that the traveltime residuals can therefore be attributed to mantle anisotropy. We construct a 3D isotropic Vp model for the study area using 2D refraction tomography models from previous studies (Eimer, 2020;Feng, 2016). The 3D model is built from east-west 2D slices, with each slice generated by (a) pulling bathymetry data for the profile, (b) determining the location of the trench axis along the profile, (c) interpolating the thicknesses of model layers (sediment, crust, subducting slab) from the active-source lines, (d) interpolating velocities from the sea surface down to the Moho from the active-source lines, and (e) incorporating prescribed mantle velocities beneath the Moho.
We prescribe mantle velocities as a function of distance from the trench axis using an east-to-west gradient because mantle velocities from the 2D refraction models are inherently not representative of the isotropic background velocity in the case of mantle azimuthal anisotropy and would bias the measurement. Velocities at the top of the mantle are set at the eastern and western edges of the model space, and values in between are interpolated using a piecewise linear gradient ( Figure S3 in Supporting Information S1). The gradient is larger west of 148.75°W, accounting for 60% of the total velocity difference between the endpoints, and smaller east of 148.75°W to make up the remaining 40%. The longitude for the breakpoint and the 60/40 balance are chosen to roughly mimic the trenchward velocity decrease at the top of the mantle from 2D refraction tomography (Eimer, 2020). The vertical velocity gradient is 0.023 km/s/km everywhere in the mantle.
We test several possible pairs of values for the eastern and western velocities, trace a random subsample of source-receiver pairs through each resulting 3D model, and choose a preferred mantle velocity model that minimizes the RMS traveltime residual for the Pn picks ( Figure S3 in Supporting Information S1). The RMS of the residuals is similar for several models with approximately the same average velocity at the top of the mantle. Our preferred isotropic model has Vp at the top of the mantle varying between 7.88 and 7.95 km/s from east to west across the study area.
The calculated delay times are inverted along with the lateral distances traveled in the mantle to obtain a model for slowness variations, dq, as a function of azimuth. We use the functional form for a transversely isotropic medium, a weighted sum of cosines and sines of 2θ and 4θ plus a constant term (Backus, 1965;Crampin, 1984;Thomsen, 1986). The mantle travel distances are calculated from the traced raypaths in the 3D isotropic model. We start by inverting all of the delay times together, and include a static term for each receiver in the inversion to capture potential contributions from unmodeled crustal structure ( Figures S4 and S5 in Supporting Information S1). These receiver statics are applied to the delay times for subsequent anisotropy inversions.
We measure how anisotropy varies spatially by partitioning delay times by longitude, and inverting subsets of the data for dq as a function of azimuth. The bins need to be narrow to measure lateral anisotropy variations, but must also each contain enough data for a robust inversion. We construct four overlapping longitudinal bins, numbered 0 (near trench) to 3 (far from trench), and a delay time is included in a bin if 90% of the source-receiver path is within the bin (Figure 1). Each bin contains a minimum of 1,600 delay times, and spans at least 0.7° of longitude. The overall conclusions described in the next section are not strongly dependent on the exact binning setup. The data in each bin are inverted for anisotropy (Figure 2, Figure S6 and Table S1 in Supporting Information S1). Coefficients for the slowness anisotropy models are converted to velocity parameters, referenced to the velocity at the top of the mantle in the center of each bin from the isotropic model.
Uncertainties are estimated for anisotropy parameters in each longitude bin using bootstrap resampling. The bootstrap scheme uses 200 random subsamples of 70% of the delay times in each bin, without replacement ( Figure S7 in Supporting Information S1). The fast direction azimuths are constrained to within ±1°.
We also measure depth variation in anisotropy near the trench. A depth is assigned to each delay time based on the maximum depth below the Moho reached by the corresponding raypath traced in the 3D isotropic model. The maximum traced depth is ∼14 km below the Moho, and 87% of paths remain above 6 km sub-Moho. Delay times within bins 0 and 1 are pooled and then divided into several depth ranges, and we invert for anisotropy in each depth range (Figure 2, Figure S8 and Table S1 in Supporting Information S1).

Constraints on Faults and Fabric From Anisotropy
Azimuthal anisotropy in the upper mantle near the Marianas Trench varies with distance from the trench axis and, near the trench, with depth below the Moho. The variations in anisotropic fast direction and in the amplitude of the 4θ component of anisotropy are indicative of bend-related faulting and cracking with orientations controlled by the trench (Miller et al., 2021). We use synthetic wavefield calculations to illustrate the combined effects of aligned faults and mineral fabrics on anisotropy for a range of faulting scenarios already determined to be consistent with 2D tomographic models (Eimer, 2020), and compare synthetics to the observed anisotropy to estimate a lower bound for the extent of fault-controlled slab mantle hydration at the Marianas Trench.

Spatial Variations in Anisotropy
The fast direction of anisotropy in the incoming plate upper mantle rotates from 124° ± 1°N in bin 3 to 153° ± 1°N in bin 0 ( Figure 2). Far from the trench, the fast direction roughly parallels the fossil spreading direction of 307°/127°N (Nakanishi et al., 1992;Seton et al., 2014), and the amplitude of the anisotropic velocity variation is ∼7.7% with respect to a mean of ∼8 km/s. This is consistent with the presence of olivine fabric from plate formation (Hess, 1964;Skemer & Hansen, 2016;Zhang & Karato, 1995) and the strength of anisotropy is similar to observations elsewhere in the Pacific (e.g., Mark et al., 2019;Morris et al., 1969) and ophiolite samples from fast-spreading environments (Ben Ismaïl & Mainprice, 1998). The near-trench fast direction is closer to the strike of the trench, indicating that trench-parallel aligned faults, cracks, and/or joints are partially overprinting the mineral fabric signal (Hudson, 1981;Miller & Lizarralde, 2016;Miller et al., 2021). The mean mantle velocity also decreases to 7.78 km/s near the trench, consistent with increasing serpentinization, and when anisotropy is factored in, the mantle velocity near the trench matches the 7.3 km/s observed along E-W refraction lines in previous work (Eimer, 2020).
Close to the trench, anisotropy also varies with depth, with a relatively large 4θ component present within 1 km of the Moho (Figure 2c). Anisotropy with strong 4θ periodicity is indicative of aligned thin cracks, in contrast to wider faults or joints which produce a predominantly 2θ signal (Hudson, 1981;Miller et al., 2021). The observed depth variation in anisotropy suggests that cracks are present in the uppermost mantle near the trench, closing within ∼1 km of the Moho, while faults or joints extend to greater depths. Farther from the trench, there is no significant depth variation in anisotropy over the depth range sampled by the data (Figure S9 in Supporting Information S1).

Synthetic Traveltimes Through Fabric and Faults
We use synthetic calculations to better understand how combining mineral fabrics and aligned faults might affect observations of anisotropy. Synthetic wavefields are propagated through 2D material models representing horizontal slices through the mantle (Cerjan et al., 1985;Kosloff & Baysal, 1982). The models have orthogonal fast and slow directions, and contain linear low-velocity fault zones oriented at an angle with respect to the background anisotropy. Fault-zone spacing is based on the observed spacing of fault scarps in bathymetry and multi-channel seismic (MCS) reflection data (Eimer, 2020). The distance between observable faults varies between 1.3 and 13.7 km, with an average spacing of ∼6 km. The background anisotropy mimics the observed anisotropy farthest from the trench. We vary fault-zone velocity and thickness to assess how faults contribute to anisotropy. Specifically, we start with three models for fault thickness and velocity that were previously found to be consistent with the 2D refraction tomography models (Table S2 and Figure S10 in Supporting Information S1) (Eimer, 2020). The fault-zone velocities are modeled after antigorite, lizardite, and serpentine gouge recovered from ODP leg 180 (Kopf, 2001). Variations on these models are used to assess the effects of fault thickness on anisotropy.
Since the synthetic calculations are for 2D representations of a horizontal slice through the mantle, they do not account for changes in fabric orientation due to plate bending. We also do not include fabric-type anisotropy within the fault zones, although serpentine minerals are highly anisotropic (Horn et al., 2020) and could potentially have a preferred orientation controlled either by preexisting olivine fabric (Wallis et al., 2011) or by fault orientation (O'Hanley, 1991). Rather, these idealized synthetic calculations enable us to explore the highly non-unique effects of fault parameters on anisotropy.
Hydration is estimated for each synthetic faulting scenario by calculating a weight percent of water at each point in the input velocity grid and averaging over the model domain. We take 8 km/s as the velocity for unaltered mantle and 5 km/s as the velocity for completely serpentinized but intact material. For points with velocities between 5 and 8 km/s we calculate a volume fraction of serpentine and corresponding weight percent of water assuming that serpentine has 13 wt% water (Carlson & Miller, 2003). Velocities less than 5 km/s are taken to indicate water-filled porosity within serpentinized material. The volume fraction of pore space is calculated from the velocity using the Wyllie time-averaged equation, which requires no assumptions about the aspect ratio of pore space (Wyllie et al., 1958).
The synthetic calculations are able to reproduce the observed rotation of the anisotropic fast direction from fossil spreading-parallel to approximately fault-parallel. Incorporating low-velocity fault zones also reduces the mean mantle velocity in the synthetics, consistent with observations (Eimer, 2020;Feng, 2016). All three initial models are able to explain the anisotropy observations well (Figure 3). While the models predict similar anisotropic velocity structure, the lizardite and gouge models are more physically likely as lizardite is favored over antigorite at the relatively lower temperatures expected in the slab (e.g., Boudier et al., 2010).
We use variations on the serpentine gouge model to look for a lower bound on hydration based on the anisotropy constraints. We test a suite of models with progressively thinner serpentinized fault zones, starting with the 500 m thickness used by Eimer (2020) to match the 2D refraction tomography. The fast azimuth is clearly rotated toward the fault strike azimuth down to 400 m thickness, but begins to shift back toward the background mineral fabric fast direction for faults 350 m or thinner (Figure 3, Table S2 in Supporting Information S1). The slowest velocities also increase to >7.4 km/s for 350 m thick faults, faster than the minimum velocities measured by refraction tomography (Eimer, 2020). This places a lower bound of ∼400 m on fault-zone thickness for faults filled with serpentine gouge, corresponding to ∼0.85 wt% water.
The anisotropic slow directions observed near the trench and predicted from synthetic calculations do not match, suggesting that the two-dimensional synthetic models do not fully capture the mantle structure. However, fast and slow direction orientations in the synthetics do match the anisotropy measured deeper than 1 km below the Moho near the trench (Figure 3a). This is consistent with the interpretation that strong 4θ anisotropy 0-1 km beneath the Moho is a signal from aligned cracks, since cracking is not included in the synthetics. The match between observations and synthetics deeper than 1 km below the Moho suggests that once cracks are closed due to increasing pressure, we are seeing anisotropy from aligned faults or joints.

Implications for Subduction Zone Water Flux
Any physical model for incoming plate hydration must simultaneously explain the spatially variable anisotropy and Vp from 2D refraction tomography, while respecting the Vsv model from surface wave tomography and constraints on fault spacing from MCS reflection imaging. Both the refraction and surface wave tomography models could be explained by either distributed mantle serpentinization or damaged and hydrated material localized in fault zones, but the anisotropy observations are best explained by faulting.
The refraction data sample only a few km beneath the Moho, and cannot provide information on how deep bend-related faulting and hydration extend. However, other observations, including Rayleigh wave dispersion and anisotropy (Cai et al., 2018) and earthquake hypocenters Emry et al., 2014), show that extensional faulting extends to greater depths than sampled in this study.
It is not possible to explain the anisotropy with bulk serpentinization by invoking serpentine fabric. The two most likely scenarios for fabric with bulk serpentinization are: (a) crystal orientations in serpentinized mantle peridotite follow those of relict olivine fabric (Boudier et al., 2010), or (b) serpentine orientations are random, destroying any pre-existing fabric (Rumori et al., 2004;Wallis et al., 2011). Neither scenario could produce the observed changes in anisotropy near the trench.
Under the constraint of the observed fault spacing, the faulting scenarios consistent with the anisotropy measurement imply that the incoming plate upper mantle at the Marianas Trench carries as little as ∼0.85 wt% water (Table S2 in Supporting Information S1) (Carlson & Miller, 2003;Wyllie et al., 1958). This represents a lower bound on water content, compared to the upper bound of ∼2 wt% water estimated from tomography models under the assumption of bulk serpentinization (Cai et al., 2018;Eimer, 2020). Global estimates of subduction zone water fluxes have typically had to rely on simplifying assumptions of bulk serpentinization in a mantle layer of prescribed thickness (e.g., Hacker, 2008). The evidence from anisotropy indicates that this global picture overestimates the mantle contribution to the subduction zone water budget at the Marianas Trench, and likely at other trenches with old, cold incoming plates. The degree of overestimation depends on several factors, including how deep serpentinization is able to extend in fault zones (e.g., Hatakeyama et al., 2017), effects of new versus reactivated faulting on the extent of hydration, and the potential effects of fracture-controlled serpentine CPO inside fault zones (O'Hanley, 1991).