Fluid transport and storage in the Cascadia forearc influenced by overriding plate lithology

Subduction of hydrated oceanic lithosphere can carry water deep into the Earth, with consequences for a range of tectonic and magmatic processes. Most of the fluid is released in the forearc where it plays a critical role in controlling the mechanical properties and seismic behaviour of the subduction megathrust. Here we present results from three-dimensional inversions of data from nearly 400 long-period magnetotelluric sites, including 64 offshore, to provide insights into the distribution of fluids in the forearc of the Cascadia subduction zone. We constrain the geometry of the electrically resistive Siletz terrane, a thickened section of oceanic crust accreted to North America in the Eocene, and the conductive accretionary complex underthrust along the margin. We find that fluids accumulate over timescales exceeding 1 My above the plate in metasedimentary units, while the mafic rocks of Siletzia remain dry. Fluid concentrations tend to peak at slab depths of 17.5 and 30 km, suggesting control by metamorphic processes, but also concentrate around the edges of Siletzia, suggesting that this mafic block is impermeable, with dehydration fluids escaping up-dip along the megathrust. Our results demonstrate that the lithology of the overriding crust can play a critical role in controlling fluid transport in a subduction zone. The lithology of the overriding plate plays a critical role in determining fluid transport in subduction zones, according to magnetotelluric imaging of the impact of the dry, mafic Siletzia terrane on fluids in the Cascadia subduction zone, North America.

and spatial context for the features seen in a single profile. The resolution of the offshore structure is also substantially improved owing to the inclusion of seafloor data. We focus on 3D variations in large-scale, deep crustal resistivity of the forearc, which are well resolved with our dataset. Figure 2 reveals a conductive body (C1) extending from ~50-70 km offshore to an abrupt landward termination at ~20 km slab depth beneath a massive resistive body (R1). We interpret C1 as representing fluid-rich sedimentary and metasedimentary rocks of the accretionary complex, and R1 as representing comparatively dry mafic rocks of Siletzia. Indeed, the geometry of R1 is broadly consistent with geological and geophysical constraints on Siletzia, including its spatial extent 36 , greater thickness in Oregon 5 and substantial disruption and unroofing of the mafic section in the Olympic Peninsula 6,37 (R2). In Oregon and southern Washington, Siletzia is very thick in the east, thinning abruptly to ~10 km in a layer that extends offshore. The conductive accretionary complex is thrusted ~50 km (more in places) under this thin layer, terminating at the nearly vertical edge of R1. In western Washington, the magnetotelluric images are in good agreement with the geometry of subducted sediments inferred from seismic tomography 24 (Extended Data Fig. 6). Offshore, the conductivity of C1 is highest from 44-47° N; however, this may reflect the limited offshore data coverage outside this area (Fig. 1). A deeper, higher-amplitude peak in conductivity (C2) occurs in all sections shown at a slab interface depth of 30-35 km, sometimes shallower. C2 commonly appears as a second peak within a broader conductive layer encompassing C1. Finally, there are conductive upwellings (C3) that rise towards the surface beneath the volcanic arc, beyond the eastern edge of Fig. 2. These arc conductive features (which can be seen more clearly in Supplementary Fig. 1) are well documented in Cascadia 32 and many other subduction zones 38 but will not be discussed further here.
Resistive material extends beyond the southern edge of Siletzia into the Klamath terranes but changes to a relatively thin upper crustal layer (R3). Resistivity variations in this area are complex ( Supplementary Fig. 1), with interlaced zones of high and low resistivity, consistent with the Klamath terrane's origin as an amalgam of pre-Tertiary oceanic units with plutonic intrusions 3 . Although the outboard edge is still dominated by conductive Franciscan material extending to the surface near the coast, the coherent accretionary complex seen along most of Siletzia was not imaged here. However, land magnetotelluric data are relatively limited in this area, and there are no offshore sites.
To provide a more complete view of the geometry of key features, we derived two plan-view images from the 3D model (Fig. 3a,b). First, we computed the depth to the bottom of the resistive bodies, defined by the 300 Ω m isosurface (Fig. 3a). This revealed a series of 30-40-km-thick resistive blocks (a-e), distributed from Puget Sound to the inferred southern boundary 36 , where a resistive curtain 'f ' extends to the plate interface. This refines previous models of the deep geometry of Siletzia based on a single seismic refraction profile 5 . These thick resistive blocks coincide with magnetic highs (Extended Data Fig. 7a), and they tend to align with high Vs above the plate interface (Extended Data Fig. 8). There is a sharp transition from the thin outboard edge of Siletzia to the much thicker blocks to the east, marking the abrupt termination of the underthrust accretionary complex (C2). One exception to this picture occurs in central Oregon, where 'g' (R4; Fig. 2) extends to the plate interface near the coast, with thin Siletzia to the east underlain by a shallow conductor (C2).
The geometry of Siletzia revealed here is consistent with its hypothesized origin as a chain of seamounts 39 61 . Coastlines and state borders are shown with white lines, and the outline of Siletzia 36 by the green line, both at z = 0. Principal conductive and resistive features discussed in text are labelled C1-C3, and R1−R4, respectively. An animation showing this 3D view from different angles is provided in Supplementary Fig. 2. plateau 40 . It is plausible that the deep resistive patches represent roots of individual eruptive centres. Some features (for example, g) may be due to post-accretion basaltic magmatism 2 . There is a strong correlation between the imaged thickness of Siletzia and modern tectonic and magmatic activity. Faults in the forearc crust 18 (Fig. 3a) do not cut the thickest resistive blocks, which also exhibit reduced recent seismicity (Extended Data Fig. 7b). Furthermore, the most westerly eruptive vents in the arc occur where Siletzia is thin or absent-for example, Mount St Helens lies between blocks c and d, and the most extreme westerly vents are in the gap between blocks b and c-suggesting that pre-existing crustal structures control modern magmatic processes 41 .
A second view of the 3D model is provided by a map of conductance of the 10 km layer above the plate interface ( Fig. 3b). High-conductance features parallel to slab contours, and aligned with the outer edge of Siletzia, extend from 45-48° N (A-C in Fig. 3b). The northern terminus coincides with the limit of densest data coverage, and is not well constrained. The band continues to the southern edge of Siletzia, with decreased conductance except for the peak at 'I' . These conductors are abruptly truncated down-dip by thick resistive blocks at slab depths ranging from 35 km in the Olympic Peninsula to 20 km in central Oregon, where the structure is more complex. Here, conductive feature 'D' is deeper (35 km slab depth) and sits in the 'cavity' formed by deep resistive features a, b, and g. Along the southern edge of Siletzia a thin band of high conductance (E) is located just up-dip and north of the resistive curtain f. Farther south, in the Klamath terrane, conductor 'F' is at 35 km slab depth.
As summarized in Fig. 4a, conductance above the plate interface is systematically enhanced at three slab depths: (1) offshore (G-I, ~17.5 km), (2) near the FMC (A-D, F; ~30 km) and (3) beneath the arc (~65 km). We focused on the forearc peaks (1 and 2), which are plotted as a function of latitude on a map of ETS density in Fig. 3c. In the Klamath terrane (north of 42° N where data coverage is good) and in the Olympic Peninsula these peaks occur at nominal slab depths of 17.5 and 30 km, with the deeper peak having a higher amplitude. Depths, and relative amplitudes, are more variable beneath Siletzia. Between 45° N and the southern Siletz margin, the seaward edge of resistive block a extends up-dip only to 40 km, and peak depths are similar. However, where resistive blocks b-d extend farther up-dip, the inboard peak is deflected seawards to slab depths as shallow as 20 km.
As in previous works 32 we interpreted conductive features as reflecting fluids released from the subducting plate through sediment compaction, dehydration of clay minerals, pore collapse in metabasalts and a series of prograde metamorphic reactions 27,42,43 . Although these processes overlap in pressure-temperature space, pulses of fluid release are expected as the mechanical thresholds and stability fields for mineral assemblages are crossed 44,45 . The observed characteristic depths for enhanced conductivity probably reflect these facies transitions. However, along-margin variations in lithology will also control the transport and accumulation of fluids, and we argue that these variations can explain the observed spatial patterns.

Implications for fluid transport and storage
Conductivity in the forearc crust varies by over three orders of magnitude, often abruptly. We believe these differences are most reasonably explained by variations in lithology and their corresponding changes in fluid content, especially given substantial seismic and geodynamic evidence for fluids in the megathrust 12,29,30 and in the overriding crust 24,25 . Assuming that resistivity can be taken as an indicator of lithology, the magnetotelluric images provide important constraints on the structure, creation and evolution of the forearc crust. In particular, our results support previous suggestions 46 that the entire Cascadia margin is being underthrust by accretionary complex material, driving uplift of the Coast Ranges 47 , and provide constraints on the 3D geometry of Siletzia. Fluid inputs at the trench include pore water and hydrous minerals, both in the sediment layer and in the underlying oceanic crust. Free water in sediments is lost through compaction at depths <10 km (refs. 27,48 ), and is not likely to be the main source for the patches of high conductivity, which are all deeper. Similarly, dehydration of clay minerals should be largely complete at temperatures of ~150 °C, which in Cascadia is well above the depth of the imaged conductors. We thus suggest that the peak centred at ~17.5 km primarily reflects pore collapse in the metabasalts, which is believed to occur at temperatures of 200-400 °C (refs. 27,44 ). Some of the fluid may also result from sediment dewatering and metamorphism and, as we argue below, in some cases may be more deeply sourced by fluids flowing up the megathrust.
The highest peaks in above-slab conductance at 20-35 km depth almost certainly reflect fluids released from metabasalt dehydration.
The subducting oceanic crust is relatively dry in Cascadia, with at most 5 vol% structural water 20 . We show in the Methods that the implied rate of fluid input, together with estimates of the fluid volume required to explain observed conductance (Fig. 4a), places a lower bound on fluid residence times of ~1 Myr. Thus, the main conductive features in Fig. 3b represent areas where fluids are sequestered, rather than source regions.
Within Siletzia the main areas of fluid storage are just up-dip of thickened resistive crust that extends to near the plate interface. We suggest that these fluids are sourced from dehydration reactions occurring at a relatively constant range of slab depths, near the FMC 49 , and that Siletzia is relatively impermeable. Where these rocks extend to near the plate interface, fluids are transported up-dip along the megathrust, ultimately leaking out along the outer edge where they accumulate in metasediments (Fig. 5). Support for this conclusion is provided in Fig. 4b, where we plot above-slab conductance (Fig. 3b) averaged in wide strips along the margin. At latitudes where conductance (and fluid accumulation) is low near the source (27.5-37.5 km), more fluid accumulates up-dip (17.5-27.5 km). Even where the roots of Siletzia do not extend to the plate interface, peak conductance always occurs up-dip of the FMC, and the peak in tremor intensity (Figs. 3b,c and 4). In fact, metabasalt dehydration reactions occur over a range of depths, near and down-dip from the ETS zone 23,43,49 . In our model (Fig. 5) dehydration fluids escaping into the crust move more rapidly up-dip than vertically owing to strong permeability anisotropy along and near the plate interface. As fluid accumulation volumes reflect an integral over source depths and time, conductance peaks are shifted systematically up-dip relative to the source region-by the greatest amount where impermeable thick Siletzia extends landwards into the source region.
Seismically imaged low-velocity layers near the plate interface in Cascadia 12,14,23,24,29,50 and elsewhere 30,51 have been interpreted either as a thin layer of fluid trapped in the upper oceanic crust 51 or a thicker zone of fluid-saturated metasediments in the overriding crust 24,52,53 ; for example, see Extended Data Fig. 6. The magnetotelluric results are more consistent with the latter view. As illustrated in Extended Data Fig. 9, a 5-km-thick layer with a porosity of 2.7-4% as estimated for the low-velocity layer beneath Vancouver Island 29,51 would have a conductance of 300-1,000 S, and cannot be present everywhere in the Cascadia forearc, as some have argued 13,14 . However, such a layer could be present in places, and a thinner or less-conductive layer could be present over much of the forearc without contradicting the magnetotelluric data.
Seismic reflection images 52,54 from southern Vancouver Island may provide insight into the nature of the conductive underthrust rocks. A zone of high reflectivity (E), interpreted as a zone of strongly sheared metasediments with sub-horizontal fluid-filled cracks, thickens from 2 km offshore where the megathrust is locked to a 5-7 km layer inland where deformation becomes aseismic 52 . The conductive zones in Fig. 3b are probably similar, representing thick layered metasedimentary packages that are ductile, strongly sheared, laminated and filled with fluids. Indeed, the E layer beneath Vancouver Island is conductive 53,55 .
As others have noted 15,18,25 , the reduction in ETS frequency between 43 and 47° N (Fig. 3c) occurs beneath the thickest parts of Siletzia (Fig. 3a). This has been ascribed to reduced fluid pressures, due either to reduced hydration 22 or permeability 16 of the subducting plate or enhanced permeability of the overriding crust 18 . However, the small seismically inferred variations in plate hydration 20,22 or sediment water content 21 (in both cases greater off Oregon) are inconsistent with increased fluid storage to the north, as we infer. It thus seems likely that fluid inputs are fairly uniform, and our inference that Siletzia is impermeable would actually suggest higher fluid pressures near the FMC in central Oregon.
A source of fluid is a necessary condition for ETS, if not sufficient alone. ETS most probably results from complex interactions in the stressed rock-fluid system 56 that involve fracturing, fluid movement and crack sealing by mineral precipitation 57,58 . Variations in lithology (that is, Siletzia versus accretionary complex), with corresponding variations in rheology, rock strength and mineral dissolution and precipitation kinetics, may play an important role in controlling ETS intensity in Cascadia. There is some evidence that fault zone heterogeneity, with stronger asperities embedded in a more easily deformed matrix, may be a requirement for ETS 56,59 . Thus, one possibility is that where Siletzia extends to near the plate interface, the subduction channel is narrower and more strongly sheared, resulting in fewer large asperities near the megathrust. It is also possible that increased transverse permeability of a more focused subduction channel allows more rapid up-dip fluid escape, and reduced fluid pressures in the ETS zone.
To make further progress, the magnetotelluric images should be more fully integrated with other geophysical constraints. The accretionary complex (mostly, but not completely, metasediments) and the competent rocks of Siletzia seem to have vastly different capacities for fluid storage, allowing magnetotelluric data to effectively map lithology, and hence to constrain spatial variations in rheology. Incorporating this information into geodynamic models 60 might provide insights into seismicity, ETS and plate coupling in Cascadia and other subduction zones.

online content
Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/ s41561-022-00981-8.

Methods
Magnetotelluric data and inversion. Our magnetotelluric dataset includes 165 new sites (including 64 offshore) from the Magnetotelluric Observations of Cascadia with a Huge Array (MOCHA) project, straddling the continental margin in the central part of the Cascadia subduction zone (Fig. 1). These were supplemented with 142 EarthScope transportable array magnetotelluric sites 35 and legacy long-period data from six 2D profiles (13 from EMSLAB 31 , 18 from CAFE-MT 34 , 13 from SWORMT 64 , 15 from KLMD 32 , 17 from POM 65 and 10 from ABCS 66 ). Most of the profile data are discussed in ref. 32 . For the MOCHA onshore sites and the five profiles, site spacing was typically 10-20 km, whereas for the offshore sites, profiles were ~40 km apart with ~20 km site spacing, except for the central profile at 44.25° N, which had a denser site spacing of 5 km. In total there were 393 long-period magnetotelluric sites, with data covering the period ranging from 7 s to 14,678 s. All onshore sites measured vertical magnetic fields; offshore sites did not. Data were processed with standard robust processing algorithms 67,68 .
For the 3D inversions we used the ModEM code 69 . This inversion code allows specification of a prior model, deviations from which are penalized by the regularization. For results reported here the prior model consisted of a layered 1D Earth with a resistivity of 100 Ω m to a depth of 410 km, 10 Ω m from there to 660 km and 3 Ω m below this (Extended Data Fig. 1a). In the first layers of the starting model we included known bathymetry with the resistivity of seawater set to 0.3 Ω m (ref. 70 ). We also included a realistic ocean sediment distribution 71 (Extended Data Fig. 1d), with the electrical resistivity of the saturated sediments set to 1 Ω m.
Finally, we included the subducting Juan de Fuca plate in the prior model as a dipping resistive layer (Extended Data Fig. 1b), as in ref. 34 . There is substantial evidence that oceanic lithosphere is generally relatively resistive 72 . Although resistive bodies are typically observed in inversions of magnetotelluric data from subduction zones 38 without imposing this structure, magnetotelluric data cannot vertically resolve resistive bodies. Imposing the seismically constrained geometry of the plate in the prior model can improve the resolution of conductive features above the plate 34 . We took the depth to the top of the shallower parts of the plate (above 90 km) from ref. 61 and extended to greater depths using seismic tomographic models 73 . The thickness of the plate was set to 50 km on the basis of published models of resistive ocean lithosphere 74 .
The model grid (Extended Data Fig. 1c) had a resolution of 11 × 9.5 km 2 horizontally, with vertical layers logarithmically spaced, starting from very fine layers (20 m) to discretize the bathymetry and the seafloor sediments. The resulting core grid was 120 × 108 × 57 cells in latitude, longitude and depth, respectively. Edges were padded with seven logarithmically increasing grid cells on all sides to extend the boundary of the study area ~600 km in all directions. We ran the inversion more than 30 times to test the effects of various parameters (for example, the thickness of the slab, the effect of offshore sites, the effects of seafloor sediments and covariance parameters) and inversion strategies. Selected results are shown in Extended Data Fig. 5.
The seafloor data presented some challenges to the 3D inversions. The lack of short-period data made it difficult to constrain shallow structures such as the seafloor sediments. Inversion tests demonstrated that including a realistic distribution of conductive ocean sediments a priori (Extended Data Fig. 1d) was essential; anomalous phases of the nominal transverse electric (TE) mode (coast-parallel electric field) impedances, both for onshore and offshore sites near the coast could not be fit without imposing this structure.
Note that including a sediment layer greatly reduced the resistivity contrast at the seafloor, which is probably an important factor in improving the performance of the inversion. Even with the sediment layer imposed, fitting the offshore TE mode data was a challenge that we found could be mitigated with a two-stage inversion strategy: we first ran the inversion with large error floors on the nominal transverse magnetic (TM) mode (10% of the off-diagonal impedance), forcing the inversion to concentrate on fitting TE mode data, for which we used a smaller (5%) error floor. In the second stage the inversion was restarted, using 5% error floors for both modes. The final normalized root-mean-square misfit with 5% error floors was 2.74, with the largest misfits for TE mode data localized on the seafloor and near the coast (Extended Data Fig. 4).
Responses for the full dataset (observed and computed) are shown as maps of apparent resistivity and phase in Extended Data Figs. 2 and 3, along with site-by-site normalized root-mean-square for each impedance component in Extended Data Fig. 4.
Residence time for fluids. Assuming a 4-km-thick hydrated layer with 5 vol% water (ref. 20 ) and plate velocity of .03 m yr −1 , the flux of bound water per metre width of plate was estimated as F = 0.05 × 4,000 × 0.03 = 6 m 2 yr −1 . To convert imaged 'patch conductance' C (conductivity integrated over the cross-sectional area) to crustal water volume V f for the same patch we used Archie's law 75 , assuming a fluid conductivity of σ f = 30 S m −1 (the salinity of seawater at ambient temperature 76 ). A simple calculation showed that for a cross-sectional patch of area A (in m 2 ) and total conductance C, , where m is the exponent in Archie's law. Assuming steady state, the minimum residence time of fluids in the high-conductivity patches is T = V f /F yr. For the main conductive features in our inverse solution typical values for cross-sectional areas were A = 4 × 10 8 m 2 and for total patch conductance C = 5 × 10 7 S m −1 . For m = 1.3, 1.5 and 2.0 we found T = 1.0, 1.9 and 4.8 Myr. Although these calculations were very rough, the conclusion that fluids must be stored for >1 Myr was robust, especially given that some fluids are subducted to greater depths to drive arc magmatism (refs. 34,77 ) or to serpentinize the forearc mantle (refs. 26,78 ).

Data availability
The data that support the findings of this study are publicly available online at http://ds.iris.edu/ds/products/emtf/. The electrical resistivity model file can be accessed online at https://doi.org/10.5281/zenodo.6303537 Source data are provided with this paper. Fig. 2 | observed (columns 1 and 3) and computed (columns 2 and 4) apparent resistivities for 4 representative periods. Left two columns are for the nominal tE mode, with current flow parallel to the coastline. Right columns are for nominal tM mode, with current flow perpendicular to coast. Plotted values are interpolated from values at sites using natural neighbor scheme. Fig. 3 | observed (columns 1 and 3) and computed (columns 2 and 4) Fig. 1). Note that for the Vp plots cool colors are low velocities (a-b), and are interpreted as subducted sediments and metasediments. these should have low resistivities (hot colors in (c-d). there is a very good agreement between geometries imaged by the two geophysical methods.

Extended Data
Extended Data Fig. 9 | Conductance (C = σ bulk H) of a fluid layer, as a function of porosity (ϕ) and layer thickness (H). We assume a fluid of conductivity σ fluid = 30 S/m (salinity of seawater at ambient temperatures 76 ) and compute bulk resistivity using Archie's law (σ bulk = σ fluid ϕ m ) for m = 1.5, 1.75, 2 for the three panels (a, b and c). Colormap for conductance is identical to that used for Fig. 3b. Even with an Archie's law exponent of m = 2 (panel c) conductance exceeds ∼ 300 S for layer thickness 51 and porosities 29 previously postulated. the Mt conductance maps (Fig. 3b) show that such a layer is not present everywhere, but a thinner layer, or one with lower σ fluid could be. Note also that conductance of observed anomalies integrated over a 10 km layer, exceed peak values shown here, requiring either higher σ fluid , lower values of m or some combination.