1.1 Evidence for Hydration of the Downgoing Plate

We use P- and S-wave arrival times from regional earthquakes recorded during the Alaska Amphibious Community Seismic Experiment (AACSE) (Barcheck et al., 2020) (Figure 1) to develop three dimensional (3-D) compressional and shear-wave velocity (Vp and Vs) tomography models (Koulakov, 2009) (Text S1–S4 in Supporting Information S1). The models clearly delineate the transition from oceanic lithosphere (Vp > 7 km/s) to upper-plate lithosphere (Vp = ∼6–6.5 km/s) (Figure 2a), which is coincident with existing slab model JB2010 (Jadamec & Billen, 2010). Within the upper ∼25 km of the subducting plate, we image a trench subparallel region of elevated Vp/Vs $$\mathit{\approx }$$ 1.87 ± 0.08, a feature resulting primarily by a decrease in Vs (Figures 2c, 2e, 2g, and 3a). The high Vp/Vs shows along-strike amplitude variations that reach the highest values of ∼1.95 ± 0.10 in the Semidi segment. Resolution tests show: (a) that along-trench variations in amplitude of the high Vp/Vs feature are resolvable at lateral scales of 50 km or larger (Text S5, Figures S6 and S7 in Supporting Information S1); (b) that the feature is not an artifact arising from smearing of high Vp/Vs in the overriding accretionary prism (Text S8, Figure S16 in Supporting Information S1); and (c) the feature is not an artifact arising from using different P- and S-wave data quantities, as it persists in an inversion where only event-station pairs with both P- and S-wave picks are used, albeit with reduced magnitude (Text S11, Figure S20 in Supporting Information S1).

Figure 3 shows crossplots of Vs versus Vp datapoints extracted from our 3-D models along the top of the slab as well as 10 and 30 km below it. Due to model discretization and smoothing, slab-top datapoints reflect a mixture of oceanic crust and mantle, and bottom of the overriding plate, while those extracted 10 km below the plate interface are more representative of oceanic lower crust and mantle. Datapoints from the slab top with depths >∼30 km have Vp and Vs increasing with depth and match —within their uncertainties (Text S6 and S13; Figures S8–S10 in Supporting Information S1)— the expected velocities for subducting mid-ocean ridge basalt (MORB) and its prograde metamorphism (Hacker et al., 2003) (Figure 3b). A portion of datapoints from depths <∼30 km, mostly with Vp < 7 km/s, depart from this trend and cannot be explained by unmetamorphosed MORB. This general pattern is also observed when a more restricted set of P- and S-wave picks are inverted, albeit with less scatter (Figure S21b in Supporting Information S1). This departure can be explained by MORB with seawater-filled crack-like porosity, based on effective medium theory (Kuster & Toksöz, 1974; Toksöz et al., 1976) (i.e., thin crack-aspect-ratios of <0.01 and rock porosities of <1%, Figure 3b). However, the full scatter (i.e., Vs > 3.7, depths >20 km) requires porosities >3% which are unrealistic for those depths (Carlson & Herrick, 1990). Thus, this departure likely also reflects heterogenous lithologies near the plate interface, including basalts, subducted and accreted sediments, and forearc terranes.

In the slab mantle datapoints cluster within a broad range of Vp and Vs values (Vp $$\mathit{\approx }$$ 6.75–8.4 km/s, Vs $$\mathit{\approx }$$ 3.5–4.9 km/s) (Figure 3c). Serpentinization of mantle peridotites reduces Vp and Vs but at different rates, resulting in elevated Vp/Vs (Christensen, 2004; Grevemeyer et al., 2018; Miller et al., 2021). Nearly all datapoints, within their uncertainties, follow a serpentinization trend between unaltered peridotite and ∼30% serpentinization. Datapoints representing the along-trench high Vp/Vs feature depart from the serpentinization trend at Vs < 4 km/s and plot within a Vp of ∼7–7.5 km/s. This pattern is qualitatively explained by effective medium theory assuming seawater-filled cracks at low porosities of <∼0.5% with very small aspect ratios (<3·10⁻³) hosted in peridotite (Text S10 in Supporting Information S1). However, the test using a more restricted set of P- and S-wave picks reduces the scatter (Figure S21c in Supporting Information S1) and, considering the data uncertainties, both serpentinization or thin-crack porosity models can explain the observations. Thus, the high Vp/Vs feature requires the presence of subducting hydrated shallow mantle, although the mode of fluid storage cannot be uniquely determined from our models. Deeper in the oceanic mantle, at 30 km below the slab top, datapoints cluster around the reference peridotite model (Vp = 8.02 km/s, Vs = 4.53 km/s) (Figure 3d), indicating little serpentinization is present and that lithostatic pressure is high enough to close seismically detectable porosity (Cheng & Toksöz, 1979).

Previous Vs and Vp/Vs models for this region derived from AACSE data support the interpretation that the high Vp/Vs feature we image is in part caused by mantle hydration. Li et al. (2024) (cf. Figs. 7, 9) and Wang et al. (2024) (cf. Fig. 4d) show a similar high Vp/Vs feature. In the Shumagin segment, Li et al. (2024) interpret high Vp/Vs as serpentinization varying from 20% at the Moho to 0% at 18 km below, noting that large discrepancies between serpentinization estimates from Rayleigh wave Vs tomography versus active-source Vp tomography suggest the existence of water in joints and cracks. However, other models (Feng, 2021; Gou et al., 2022) do not include a trench-parallel reduction in incoming mantle Vs, or increase in Vp/Vs, comparable to our observations.

Despite the ambiguities in possible mechanisms, we suggest that the trench-parallel high Vp/Vs feature represents primarily a crustal and mantle pore-water reservoir, with serpentinization playing a secondary role and/or being confined to fault zones (e.g., Miller et al., 2021). The main argument for pore-water is the observation that the high Vp/Vs feature fades down-dip and disappears at depths >∼30 km (Figures 2g, 2h, and 3a). If the incoming mantle were pervasively serpentinized, its Vp/Vs signature should be observed deeper in the mantle because metamorphic dehydration of serpentine is not predicted to occur until 80 km depth (Abers et al., 2020). We therefore interpret the lack of high Vp/Vs deeper than ∼30 km as resulting from crack closure and pore-fluid expulsion from the mantle at ∼20–30 km depth.

1.2 Evidence for Fluid Release From the Downgoing Plate

We find evidence for slab fluid release by examining our models along the megathrust interface (Figure 4, Figures S17 and S18 in Supporting Information S1). Fracture zones (FZs) are high-porosity features capable of storing water at crustal and mantle depths (Roland et al., 2012; Wang et al., 2022). Off Kodiak Island, the incoming plate is traversed by a series of FZs (Figure 1a), and an active-source tomography model shows reduced crustal and mantle Vp in this region, suggesting increased hydration of the incoming plate (Acquisto et al., 2024). In our model, we find three distinct, elongated low velocity anomalies in the shallow part of the megathrust (A3, A4, A5, Figure 4) that spatially correlate, and have azimuths aligned with these FZs. Thus, anomalies A3–A5 most likely represent water released from these hydration anomalies, likely pore fluids, into the megathrust. Similar features, though with less distinct FZ geometry, are also present in the models of Gou et al. (2022) and Wang et al. (2024).

Beneath the Shumagin Islands between 162°W and 157°W, our models have two bands of trench-parallel negative anomalies (A1, A2) along the plate interface at 15–30 km and 40–60 km depth (Figure 4, Figures S17 and S18 in Supporting Information S1). Synthetic tests show they are well resolved by our data set (Figure S19 in Supporting Information S1). Similarly to anomalies A3–A5, we interpret A1 as representing entrapment of fluids at the plate interface sourced from the closing of pore spaces in the crustal and shallow mantle water reservoirs due to increasing confining pressure (Cai et al., 2018; Cheng & Toksöz, 1979). Anomaly A2 and a similar feature further northeast (A6), being distinctly located deeper down-dip, likely have a different fluid source such as metamorphic crustal dehydration reactions (Hacker et al., 2003; Peacock, 1993). Although significant metamorphic dehydration is not expected until 80 km depth (Abers et al., 2020), fluids may be transported up-dip as inferred at other subduction zones (Hyndman et al., 2015). Alternatively, anomalies A2 and A6 may partially result from an anomalously deep fluid source within the forearc slab mantle (Cordell et al., 2023) (Figure S25 in Supporting Information S1).

Anomalies A1 and A2 are not imaged east of ∼157°W, and only a small patch of negative δVp, δVs, and high Vp/Vs at 35 km depth is present. This could indicate reduced hydration and fluid release from the slab, consistent with a change in plate fabric orientation less favorable to bend-faulting induced hydration (Shillington et al., 2015) (Figure 1b). Or alternatively, it could reflect more efficient draining of the plate interface through the overriding plate in the Semidi than in the Shumagin segment. Backstop splay fault zones are present along the AP margin and are associated with seafloor fluid venting, thus providing a mechanism for fluid drainage from the plate interface. Fault zones in the Semidi segment are located ∼40–50 km landward from the trench and are therefore more favorably located for draining the megathrust than in the Shumagin segment, where the fault zone is located closer to the trench (∼10–15 km) (Bécel et al., 2017; von Huene et al., 2021). Lack of fluid-related anomalies in the Semidi segment, compared to the Shumagin, is a feature in all previous 3-D seismic models (Feng, 2021; Gou et al., 2022; Li et al., 2024; Wang et al., 2024).

1.3 Megathrust Structure and Implications for Earthquake Behavior

There is a striking correlation between the location of anomalies A1 and A2 and the up- and down-dip limits of the rupture zone of the 2020 M_w7.8 Simeonof earthquake (Crowell & Melgar, 2020; Liu et al., 2020; Xiao et al., 2021; Ye et al., 2021; Zhao et al., 2022) and the 2023 M7.2 (USGS, 2024), with a similar pattern in the vicinity of the 2021 M_w8.2 Chignik earthquake (Figure 4, Figure S24 in Supporting Information S1). Thus, we interpret anomalies A1 and A2 as rupture barriers for the Simeonof, Chignik, and 2023 M7.2 earthquakes. Although these patterns are observed in the δVp, δVs, and Vp/Vs models, we focus on the δVp anomalies because there is more P-wave data and the correlation is the clearest (Figures S17 and S18 in Supporting Information S1).

The Simeonof earthquake occurred ∼50 km east of the Shumagin islands and ruptured northwestward extending into the Shumagin segment (Crowell & Melgar, 2020). Published co-seismic slip models show there was no significant slip shallower than 15 km and that the largest along-strike slip was confined between ∼22 and 43 km depth (Figure 4) (Bai et al., 2022; Crowell & Melgar, 2020; Liu et al., 2020; Mulia et al., 2022; Xiao et al., 2021; Ye et al., 2021; Zhao et al., 2022). Understanding why the Simeonof earthquake was constrained to ∼22–43 km depth is challenging. Geodetically constrained locking models in this region are consistent with a partially locked trench, with coupling decreasing down-dip, or with a Gaussian locking distribution, with the highest locking up-dip from the maximum co-seismic slip (white rectangle, Figure 4) (Drooff & Freymueller, 2021; Li & Freymueller, 2018; Xiao et al., 2021). Higher shallow coupling (i.e., accumulated slip deficit) should favor rupture to the trench (e.g., Melgar et al., 2016), particularly if clay-rich subducted sediments, which are expected to be velocity-weakening at high slip velocities (Faulkner et al., 2011), are present at the plate interface (Creager et al., 1973; Li et al., 2018). However, it should be noted that coupling amount near the trench and whether the Simeonof event involved afterslip to the trench, are both poorly constrained due to modeling assumptions and lack of seafloor geodetic observations (Crowell & Melgar, 2020; Xiao et al., 2021; Zhao et al., 2022). One explanation for the Simeonof rupture at ∼22–43 km depth is that the largest slip occurred in a velocity-weakening asperity bounded by velocity-strengthening regions both down-dip from ∼40 km and up-dip from 30 km to the trench (Zhao et al., 2022). Based on the distribution of anomalies A1 and A2 along the slab interface, we argue that, in addition to the velocity-weakening to strengthening transitions proposed by Zhao et al. (2022), fluid release from the slab—and thus elevated pore pressure at the plate interface—is a major barrier that limited the rupture extent of the Simeonof earthquake.

Increased pore pressure—and thus reduced normal stress—in a velocity-weakening region will not stop rupture propagation, but it would make earthquake nucleation more difficult by increasing nucleation length (e.g., Rubin & Ampuero, 2005). If the nucleation length is larger than the region of elevated pore pressure, this region would be too stiff to slide unstably and thus could be freely slipping (Ruina, 1983). Assuming the up-dip velocity-strengthening region extends from the trench to ∼25 km depth, where it loads the velocity-weakening region below, we posit that anomaly A1 represents a region that is at least in part velocity-weakening but freely slipping due to a large nucleation length from increased pore pressure (Figure 5). Nucleation of an event such as the Simeonof earthquake would only occur once a creep front, extending down-dip from the creeping main zone of fluid accumulation into the region where nucleation length is shorter, becomes unstable. These same arguments apply to the down-dip rupture extent and anomaly A2, with the actual transition from velocity-weakening to strengthening friction, or from frictional to viscous deformation occurring at ∼50 km depth. The freely slipping regions with elevated pore pressure have no slip deficit and thus provide a rupture barrier, arresting co-seismic slip along the boundary of the low δVp anomalies, explaining the limited depth range of the Simeonof earthquake.

An exception to this pattern occurs locally east of the Shumagin Islands at 158.5°W, 55°N where the models of Liu et al. (2020), Ye et al. (2021), and Xiao et al. (2021) show slip extending up-dip and where the amplitude of anomaly A1 is diminished from negative ∼12% to ∼9% (although some of this amplitude variation may be a resolution effect, Figure S19 in Supporting Information S1), and therefore pore pressure is inferred to be lower (Figure 4, Figure S24 in Supporting Information S1). The patch with subdued amplitude in A1 hosts the Simeonof hypocenter, suggesting that local patches where lower pore pressure makes the megathrust more unstable, may serve as preferential nucleation regions.

The absence of pronounced negative δVp anomalies, and therefore inferred rupture barriers, in the region between 157.5°W and 155°W coincides with the area of largest slip during the M_w8.2 1938 earthquake (Estabrook et al., 1994; Freymueller et al., 2021; Johnson & Satake, 1994). Thus, efficient drainage of the plate interface and/or reduced fluid release in this region may have favored rupture at shallow depths, triggering a Pacific-wide tsunami (Freymueller et al., 2021; Johnson & Satake, 1994).

Fluid release in subduction zones is ubiquitous (van Keken et al., 2011) and thus linkages between fluid distribution along the plate interface and great earthquake slip extent are expected globally. Our observations and conceptual fluid-mediated rupture barrier model thus offer a new framework to better understand great earthquake rupture size and associated hazards across subduction zones, especially as more amphibious seismic and geodetic studies are conducted (Hilley et al., 2022).