Zones of Weakness Within the Juan de Fuca Plate Mapped From the Integration of Multiple Geophysical Data and Their Relation to Observed Seismicity

This study aims to explain the nonuniform earthquake pattern along the Cascadia Subduction Zone. In particular, we investigate the relationship between the tectonic features of the subducting oceanic Juan de Fuca slab and the onshore seismicity pattern. We have integrated multiple geophysical data sets toward three general objectives. The first study intends to study variations in physical properties along three 2‐dimensional models through regions of different seismicities that combine public gravity, magnetic, and seismic data sets. These models reveal multiple zones of decreased crustal density that we interpret as regions of weaker oceanic crust. The second objective is to delineate major tectonic features by performing spatial analysis of potential fields. The overall methodology comprises gravity and magnetic data filtering, followed by lineaments mapping and cross‐referencing interpretation with available seismic reflection data. This process allows delineating zones of crustal weakness by extrapolating outside our three 2‐D models. We also map multiple seamounts that appear to cluster along identified zones of weaker crust. Third, we investigate the relationship between the mapped tectonic elements, namely the zones of weak crust with accompanying seamounts, and the observed seismicity trends within the subducted slab. The alignment between those suggests that mapped weak crust zones and associated seamounts may have an influence on the overall subduction process. As more of these structures are heading toward the Washington portion of the margin than to the Oregon portion, more earthquakes are observed in the north than in the south.


10.1029/2023GC010943
2 of 18 want to map the major tectonic and topographic features of the subducting plate and explore its relationship to the seismicity pattern along the subduction margin.
The tectonic complexity of the CSZ is manifested in various geological structures.The young oceanic crust is complicated by known propagator wakes (PWs, shown as light blue polygons in Figure 1) and abundant seamounts.PWs are primarily mapped from the large offsets in sea floor magnetic anomalies and represent multiple broad-scale V-shaped features in map view (Hey, 1977;Wilson, 2002; Figure 1).These structures initiate at the mid-ocean spreading center and relate to the propagation of one ridge segment at the expense of the adjacent one (Harper et al., 2021;Hey, 2021).On the other hand, the onshore region is dominated by accreted terranes, such as the thick, basaltic Siletz terrain (shown by the dashed yellow outline in Figure 1).To unravel the interplay between these domains, it is necessary to establish the crustal heterogeneities of the Juan de Fuca (JdF) plate from the spreading center to the subduction zone.Therefore, our first objective is to develop three 2-D plate-scale geophysical models from multiple data sets over the existing seismic reflection profiles.The models incorporate seismic constraints for major tectonic elements to better understand the density variations based on observed potential field data within the JdF plate and the accreted Siletz terrane.Because our models utilize several data sets to determine subsurface structures and rock properties, we refer to them as integrated geophysical models in this paper.
Our second objective is to use gravity, magnetic, and seismic data to delineate major tectonic structures within the JdF plate through spatial analysis.This research builds on the established PWs from Wilson (2002).These features are referred to interchangeably in the literature as PWs or pseudofaults (Harper et al., 2021;Hasselgren et al., 1992;Hey, 1977Hey, , 2021;;Wilson, 1988Wilson, , 2002)).In our spatial analysis, we identify several new lineaments in magnetic data (in addition to previously established PWs), which we call pseudofault lineaments or PSFs.We refer to the previously established broad-scale structures as PWs, and newly mapped features as PSFs to avoid any ambiguity.In this study, we investigate the spatial relationship between the modeled lower density zones and the newly mapped PSFs.We extrapolated the mapped tectonic elements (i.e., modeled lower densities from the first objective) to the entire oceanic crust via lineaments observed in a filtered magnetic field.Another part of our second objective was to map buried seamounts from gravity and seismic data.In any subduction region, the topographic features on subduction interface, such as seamounts or depressions, greatly influence the overall rupture behavior (Ding & Lin, 2016;Hicks et al., 2014;Ruh et al., 2016;Wang & Bilek, 2014) by acting either as a facilitator to small and large earthquakes (Abercrombie et al., 2001;Bilek et al., 2003;von Huene et al., 2000;Husen et al., 2002;Mochizuki et al., 2008;Nakatani et al., 2015;Van Rijsingen et al., 2018;K. Wang & Bilek, 2011) or as a hindrance to any seismic slip (Kodaira et al., 2004) and sometimes lead to great earthquakes acting as strong asperity (Cloos, 1992;Scholz & Small, 1997).Moreover, there is evidence of subducted seamounts offshore north-central Oregon that can be spatially correlated with small earthquakes with magnitude less than 4.5 (Morton et al., 2018(Morton et al., , 2023;;Tréhu et al., 2012;Trehu et al., 2015), suggesting that they may be helping to release energy underneath the accretionary prism.Therefore, mapping buried seamounts in the non-subducting part of the JdF plate is crucial for understanding how these features can impact the overall seismicity pattern along the CSZ.
Our final objective is to evaluate the spatial relationship between mapped lineaments (PWs and PSFs) that correlate with lower density oceanic crust and seamounts, and the seismic activity observed in the continental domain of the CSZ.We limit our analysis to earthquakes within the subducting JdF slab, as we were concerned only with Figure 1.Map of the Cascadia Subduction Zone along with its tectonic elements.The inset map shows its location relative to the North American continent.The 2D integrated and gravity modeling profiles of this study are shown in this map (dark brown lines).Also, it has several public seismic profiles that are utilized in this research (shown with green, for reflection, and purple, for refraction).The modeling profiles are over the existing seismic reflection lines or align with them so that reflection data can be used as constraints in the integrated modeling.Seismicity from 1900 to present and above 4 magnitudes on the continental domain is represented by circles scattered from British Columbia in the north to Northern California in the south (University of Washington, 1963;U.S. Geological Survey & Earthquake Hazards Program, 2017;USGS Menlo Park, 1966).Several bathymetric seamounts are also pointed out in this map.These seamounts are more apparent on the Pacific plate; however, some are also evident on the Juan de Fuca plate.The boundary of the Siletz terrane is from Phillips et al. (2017).The propagator wakes are mapped on the basis of Wilson (2002).

Data and Methodology
To achieve our goals, we integrated multiple geophysical data sets in our analysis.In particular, bathymetry data (Figure 1) from Smith and Sandwell (1997), magnetic anomalies from the United States Geological Survey (USGS) magnetic data repository (Bankey et al., 2002;Figure 2a), gravity field from Sandwell et al. (2014) shown in Figure 2b (version 29.1), seismic reflection and refraction data from Han et al. (2016), Carpenter (1971), Parsons et al. (1999) and Gerdom et al. (2000), as well as earthquake hypocenter locations from the USGS catalog (University of Washington, 1963;U.S. Geological Survey & Earthquake Hazards Program, 2017;USGS Menlo Park, 1966).Sandwell et al. (2014) reported that the wavelength band of satellite gravity data that we used in our study allows for the investigation of structures as small as 6 km.The accuracy of this data set with respect to marine gravity was assessed by two recent studies.In particular, Sandwell et al. (2021) used a very precise marine gravity survey data set in the Gulf of Mexico to compute the median deviation between satellite and marine data sets of 1.33 mGal.Rouxel et al. (2023) compared satellite gravity with several marine gravity surveys over the Atlantic Ocean and concluded a dispersion between data sets within 1 mGal starting from satellite version 23 (we use version 29.1).Therefore, the observed local gravity anomalies such as local gravity lows up to ∼15 mGal and around 10-25 km wide over the known PWs (Figure 2a) are real and require geological explanations.Moreover, magnetic data (Figure 2b) clearly show magnetic stripes over the oceanic crust and apparent distortions associated with known PWs.
Our study utilizes three different methodologies, one per each objective.The first objective is achieved via integrated geophysical modeling along existing seismic reflection profiles.The primary goal of this modeling is to analyze the variations in the physical parameters of subsurface rocks (density and magnetic susceptibility) along the margin.We built the models in the GM-SYS module of the Geosoft software in 2D approximation that is sufficient for presented plate-scale modeling (Geosoft, 2021).Each model consists of the following geological layers: water, sediments, oceanic crust, upper mantle, accreted oceanic terrane, and the continental crust complicated by arc volcano lithologies.On the vertical axis, the models reach up to 90 km below the mean sea level.We use seismic reflection and refraction data (Carpenter, 1971;Gerdom et al., 2000;Han et al., 2016;Parsons et al., 1999) to constrain depths to the subsurface layers along the modeled profiles, while physical properties of subsurface layers (densities and magnetic susceptibilities) were modeled to fit both potential fields (i.e., gravity and magnetic data shown in Figure 2).In our modeling, we did not try to achieve a perfect fit between observed and computed potential fields as the data have some errors mentioned above, as well as the modeling assumption, such as 2D approximation and the uniform physical properties within each modeled block, leaving room for some uncertainty.Nevertheless, our target is to compute anomalies of comparable amplitude and with the same wavelengths and phases as the observed potential fields.In general, we try to keep the mismatch between observed and computed gravity within 2 mGal in the oceanic part of the model where we resolve 15-30 km long subsurface structures (e.g., PWs).In the continental part, we allow a larger error margin (5-25 mGal) because of the lack of dense seismic constraints.In this case, we primarily aim to fit the shape of observed anomalies and allow a larger mismatch in amplitude as we prefer to keep the physical properties of the corresponding layers the same for all three lines.For magnetic modeling, we try to explain the amplitudes and shapes of observed magnetic signals with the simplest possible model (i.e., a single magnetic susceptibility value for the entire oceanic crust).
The top layer in all models (Figure 3) is seawater with known physical properties, that is, 1.03 g/cm 3 density and zero magnetic susceptibility (Telford et al., 1990).Densities and magnetic susceptibilities of subsurface rocks were modeled to be within the range for the expected lithologies of each layer.For oceanic crust, our modeled density ranges from 2.8 g/cm 3 near the spreading center to 2.9 g/cm 3 in the subduction zone, which agrees with the average 2.85 g/cm 3 density for the entire ocean crust modeled by Blakely et al. (2005) and Finn (1990).In our models, this density range for oceanic crust represents a general increasing trend with occasional lower densities.We conduct a series of sensitivity studies to establish how gravity anomalies change with the density variations within the oceanic crust.In the test, we change the density value of the lower density region until its density matches the surrounding crust and then observe the mismatch between the observed calculated gravity.If the mismatch amount is significant (such as above 2 mGal to the west of the deformation front), it requires a geological explanation.For example, this sensitivity test over the westernmost lower density zone of the Oregon model (Figure 3a) reveals a 6.5 mGal mismatch if we make the density the same as its surrounding crust (Figure S1 in Supporting Information S1).As the geometries of our layers in this region are constrained by seismic refractions, the only parameter we can change is the sediment or crust density.It is very unlikely that densities of the sedimentary layers have vertical heterogeneities, as these would be evident in seismic data, but are not.The density distribution of other tectonic elements (e.g., accretionary prism, Siletz terrane, etc.) is also carried out on a trial and error basis, where their assigned densities generally agree with Blakely et al. (2005).The overall pattern in densities is also consistent with the seismic velocities derived from refraction data (Gerdom et al., 2000;Parsons et al., 1999).These seismic velocities frame the internal structure of the nearshore and onshore tectonic elements and we utilize them to make a more detailed model of each tectonic feature.The modeled densities of crustal blocks and the range of densities for the tectonic features within the subduction zone are in general agreement with Blakely et al. (2005).In the oceanic domain, the only magnetic layer is the JdF crust, which is modeled with a single value of magnetic susceptibility, but with different signs to account for positive and negative magnetic polarities of the ambient field.We keep the physical properties of the corresponding subsurface layers the same for all modeled profiles to ensure the consistency of the modeling.
For our Objective 2, we utilized the spatial analysis of potential fields where we filter gravity and magnetic fields to enhance lineaments that can be traced in the study area and correlated with geological structures.Prior to this analysis, potential fields were corrected (with Bouguer correction for gravity grid and differential reduction to pole for magnetic grid (Arkani-Hamed, 2007)), detrended, and filtered.For Bouguer correction (BC), we use the equation: where Δρ is the density difference between the first two layers (for ocean they are water and oceanic sediments, and for land, they are air and continental materials), G is the gravitational constant, and H represents the thickness of water.We assume Δρ as −0.97 g/cm 3 , considering the density contrast between water (density of 1.03 g/cm 3 ) and sea-floor sediments (assumed density is 2 g/cm 3 ).To reduce the magnetic anomaly to the pole, we computed magnetic inclination, declination, and total magnetic field using the International Geomagnetic Reference Field (Alken et al., 2021) with a magnetic epoch of 2000.
We calculate the regional trend for each potential field by continuing upward the data 1.25 km for gravity and 3 km for magnetic.These values are obtained from the trial-and-error exercise where we test the resulted grids against structures such as bathymetric seamounts and PWs and decide which residual fields highlight the shallow geology best.We then subtract the regional trend from the unfiltered potential field data (i.e., data containing deeper signals of lower frequency) to retrieve residual anomalies to represent the shallow subsurface structures (Figures 4 and 5).After testing many types of filters for the residual potential fields, we conclude that the tilt angle (referred to as tilt derivative in Oasis Montaj Geosoft software) is the most useful one for mapping the PWs and PSFs in the oceanic domain from magnetic data (Salem et al., 2008;Verduzco et al., 2004), as shown in the Figure S2 in Supporting Information S1, while the Bouguer residual map is the best for outlining the seamounts.
For our final objective, we correlated the interpreted tectonic structures of the JdF plate, namely the known PWs and our mapped PSFs, and seamounts, with the earthquakes associated with the subducting slab.We selected earthquakes with magnitude 2 and up from the USGS earthquake catalog (University of Washington, 1963) over the continental domain that occurred from 1900 to 2023.We then filter out the events with a location uncertainty of 2 km or greater in both vertical and horizontal directions.We isolate the remaining earthquakes to those within the Wadati-Benioff zone based on the slab geometry from McCrory et al. (2012).To understand their clustering, we performed two statistical analyses.The first method, Kernel density estimation, allows the recognition of broader earthquake patterns.The second method, K-means clustering, allows the grouping of earthquakes based on certain criteria.In our analysis, we group seismic events within certain latitude intervals; based on the number of groups matching the count of PWs and PSFs that end up within that interval during subduction.In particular, we assigned five different zones (between latitudes 43°-44°, 44°-46°, 46°-47°, 47°-48°, 48°-49°), and within these zones we specified the number of clusters as follows: 1, 2, 2, 3, and 2 clusters, respectively.Subsequently, regression analyses were performed within those grouped earthquakes to identify a statistical trend (i.e., direction) within each cluster.

Integrated Geophysical Modeling Along Three Plate-Scale Profiles (Objective 1)
We develop three regional models that extend from the JdF oceanic spreading center in the west to the continental portion of North America in the east (see location in Figures 1 and 2). Figure 3 shows the set of our integrated geophysical models starting from the southern one (Oregon model; Figure 3a), the central one (Central model; Figure 3b) and the northern one (Washington model; Figure 3c).In our 2-D plate scale models (Figure 3), the density of the undeformed sedimentary rocks overlying the JdF oceanic plate ranges from 2.00 to 2.40 g/cm 3 (corresponding to a P-wave velocity of 2.5-5 km/s from Gerdom et al. (2000) and Parsons et al. (1999)) and accounts for general compaction with depth trend within these sediments.In the accretionary prism, the density values range from 2.45 g/ cm 3 for the shallowest layer to 2.75 g/cm 3 for the deepest layer, which is required for matching the observed gravity data.This density range for the prism sediments also agrees with seismic studies from Gerdom et al. (2000) and Parsons et al. (1999) both in terms of internal structure and velocity range.The high density value (2.75 g/cm 3 ) for the deepest sedimentary layer relates to metamorphic changes at 10-15 km depth within the accretionary prism (Hyndman, 1988;Lewis et al., 1988;Nedimović et al., 2003).It also corresponds to gravity modeling by Finn (1990) over Washington, where the same density value (2.75 g/cm 3 ) was assigned for the deepest part of prism sediments.The density range within the accretionary prism is consistent for all three models presented in this study.
Based on the seismic velocity values from Gerdom et al. (2000), the Siletz terrane in the Oregon model has been divided into the shallower upper Siletz terrane (2.85 g/cm 3 and 0.005 cgs) and the deeper Lower Siletz Terrane (3.085 g/cm 3 and 0.002 cgs) (Figure 3a).These modeled physical properties of the accreted Siletz terrane are consistent with its lithology (Parsons et al., 1999;Phillips et al., 2017) and observed P-wave velocities from seismic refraction data (i.e., 6-7 km/s in the Siletz terrane region from Parsons et al. (2006)).Also, the 2.85 g/ cm 3 density of the shallower Siletz terrane in Oregon (Figure 3a) corresponds to the 2.92 g/cm 3 density of the upper Siletz terrane from Finn (1990).The Siletz terrane in the Washington and Central models has a density of 3.09 g/cm 3 (Figures 3b and 3c), which corresponds to the modeled 3.08 g/cm 3 density in the lower Siletz terrane  3 for each modeled block are shown with a light orange color in the subsurface model.Only for "a" and "c," the magnetic susceptibility values in cgs are also shown in the subsurface model with a light blue color.However, magnetic susceptibilities are not shown for every modeled element, especially where they are obviously zero, such as seawater.To account for the magnetic reversals in sea-floor, we use both positive and negative magnetic susceptibility interchangeably while modeling the JdF crust.The red star in the upper panels represents the tie point which shifts the calculated gravity and magnetic anomaly to the datum of the observed data.from Finn (1990), where they developed gravity models over Washington.In our plate-scale models (Figure 3), the easternmost part of the study area contains an active volcanic arc known as the Cascadia arc.Across all our models, the Cascadia arc consistently exhibits lower density compared to the Siletz terrane.In the Oregon model, the Cascadia arc has a density of 2.96 g/cm 3 , which is 0.125 g/cm 3 lower than the Siletz terrane (Figure 3a).Similarly, in the Central and Washington models, the densities of the Cascadia arc are 3.01 and 3.00 g/cm 3 , respectively, which are 0.08 and 0.09 g/cm 3 lower than the Siletz terrane (Figures 3b and 3c).This variation in density reflects the contrasting lithological composition and geological characteristics between the Cascadia arc and the Siletz terrane.Because this arc complex is in a subduction setting, the volcanic chambers of this region have a more intermediate composition than basalts (Green & Harry, 1999), which should exhibit lower density than the mafic Siletz terrane.In the Central and Washington models (Figures 3b and 3c), the transition line from the Siletz terrane to the Cascadia arc has a good correlation with the extent of the Siletz terrane (Figure 1), seismic tomography study by Parsons et al. (1999) and gravity modeling of Finn (1990).The modeled density transition also aligns with seismic velocities from Parsons et al. (1999).That study reported a significant depth increase in the 6.5 km/s contour line by ∼6 km on the transition line.This means that at an approximate depth range from 7 to 14 km, the P-wave velocity changes from 6.5 km/s in the Siletz terrane to 6 km/s in the Cascadia arc, which is indicative of lower density in the Cascade arc than in the Siletz terrane.Between around 5 to 15 km depth in the Cascadia arc region, Parsons et al. (1999) also noticed a vertical transition zone from 6 to 6.5 km/s, which they interpret as silicic intrusive rocks based on the laboratory analysis on plutonic inclusions of Mount St. Helens lava dome from Paine (1982).In our Oregon model (Figure 3a), we see that the Cascadia arc density has a lower value than the Washington region and a much higher difference from the Siletz terrane density.This finding is consistent with the observed lithology between Washington and Oregon.In Oregon, the arc-related lithologies mostly include younger basalts and rhyolites (0-8 million years ago), whereas in Washington, the materials include many older basaltic rocks (12-17 million years ago) that are mostly related to Columbia flood basalts (Christiansen et al., 2002).Furthermore, a low magnetic susceptibility in the basin overlying the Siletz terrane and Cascadia arc in the Washington model (Figure 3c) can be explained by the presence of the Columbia flood basalts in that area (Wolff et al., 2008).
The mantle has a constant assumed density (i.e., 3.3 g/cm 3 ) and no magnetic susceptibility (i.e., 0 cgs) values for all models (Telford et al., 1990) with the exception of the lower density of mantle rocks beneath the JdF ridge (i.e.,  (Carpenter, 1971).The profile for this seismic reflection line is also shown in Figure 1.(c) Image of the seismic reflection line 10 and few interpretations of its geological features.The abbreviated text "BT SM" and "BR SM" represent bathymetric and buried seamounts respectively.10.1029/2023GC010943 9 of 18 upwelling magma under the active seafloor spreading) that was required to fit the observed gravity anomaly.Over the JdF spreading ridge, the modeled crustal density is 2.8 g/cm 3 in all three models and magnetic susceptibility is 0.002 cgs.In the Oregon model (Figure 3a), there is an active volcanic feature, the Axial seamount, located over the spreading center (Chadwick et al., 2005).The seismic reflection profile of Han et al. (2016) provides a partial image of the magma chamber beneath the ridge.Therefore, we determined the geometry of the conical oceanic block under the spreading ridge from the modeling of gravity and magnetic anomalies.We model a similar decrease in mantle material beneath spreading centers in other models (Figures 3b and 3c).
The density of the JdF oceanic crust increases toward the subduction zone from 2.8 g/cm 3 in the spreading center and reaches 2.9 g/cm 3 under the Siletz terrane in all three models.This gradual increase in crustal density in the oceanic domain is due to the crust becoming older, colder and thus denser as it heads over to the subduction zone.This trend agrees with Carlson and Herrick (1990), who measured crustal densities in the Atlantic Ocean and found lower porosity and higher bulk density in the drilling location away from the spreading center compared to the drill site which is situated near the mid-oceanic ridge.In our modeling, this increase in oceanic crustal density toward the coastline has a few aberrations in all three models that are required by the gravity data.In particular, we must incorporate several blocks of relatively less dense oceanic crust in all three models to satisfy gravity anomalies (transparent gray boxes in Figure 3).The density contrasts with the adjacent crustal blocks vary from 0.02 to 0.06 g/cm 3 for all low-density regions.In order to assess the correlation between these zones, we perform the spatial analysis of potential fields described in the next section.

Mapping of Pseudofault Lineaments
We show two types of structures on the magnetic data (Figure 4).The first one consists of the previously known PWs that are traditionally mapped from significant offsets in magnetic stripes (Figures 4a and 4c) and represent zones of disturbed crust related to propagate spreading centers.In addition to those, we also identify a new lineament type from magnetic data that is associated with smaller-scale discontinuities (i.e., shorter offsets and smaller disturbances) in the sea-floor magnetic stripes.We name the newly identified smaller features pseudofault (PSF) lineaments, while we refer to the previously established larger structures as PWs.To select these lineaments, we consider the following types of disturbances in seafloor magnetic anomalies: gaps in the stripes, changes in the stripe orientation from ∼North-South to ∼East-West, and a considerable reduction of the stripe's width.Notably, we interpret these newly interpreted tectonic structures (PSFs) as line features.In reality, they represent zones with a certain width (smaller than the width of PWs); however, the data used in this study do not have the resolution to map them as zones, so they have been approximated as line features.For effective categorization, we numbered each newly identified PSF lineament from 1 to 13 starting from the north (Figure 4).We use Roman numerals to denote the previously interpreted PWs with the most northern one numbered as "i" to the "v" for the most southern one (Figure 4).Three of the interpreted lineaments are not crossed by any of our three plate scale models (i.e., PSF lineaments 1, 12, 13).Two important observations related to these lineaments are (a) the western end of an interpreted pseudofault lineament originates at the break or offset of either spreading centers or of the transform fault, and (b) they mostly correlate with the zones of lower densities determined from our 2-D modeling (shown as white bars in Figure 4).

Mapping of Seamounts
Seamounts represent another type of key tectonic structure that we investigate in this study.Some of the seamounts are evident in the bathymetry map (Figure 1); we refer to those as "bathymetric" seamounts.In addition to these exposed seamounts, sparse seismic reflection data, as shown in Figure 5, reveal several buried seamounts on the JdF plate.We used residual gravity data to map these buried seamounts.We did not utilize magnetic data because the strong magnetic reversals of the seafloor conceal the signals from seamounts.We have noticed that known bathymetric seamounts have a characteristic Bouguer gravity signal.Specifically, the peripheries of the bathymetric seamounts show relatively higher gravity values in the Bouguer anomaly map, while the central locations have a low gravity (Figure 5).We explain this with the density variations associated with the seamount.In particular, as the seamounts originate from mantle material that seeps through structurally weak crustal zones, their rocks are generally younger (i.e., less dense and hotter) than the surrounding oceanic crust, resulting in lower Bouguer gravity associated with the seamounts (Figures 5a and 5b).Moreover, this characteristic gravity signal was then justified with seamounts interpreted from seismic reflection images (Figure 5c).Based on this distinctive gravity signature (i.e., a gravity low surrounded by relatively higher gravity values), we interpret 32 seamounts from the residual Bouguer gravity anomaly (Figure 6).We categorize our identified seamounts into three categories based on the confidence of interpretation ranging from most confident (i.e., "bathymetric"), less constrained (i.e., "seismic," imaged in seismic data) and the least confident "gravity"-based seamounts that are not seen in bathymetry data (i.e., covered by sediments), and do not have seismic data to confirm their presence.All three types of seamounts have similar gravity signals, which are further illustrated in graphs in Figure S3 in Supporting Information S1.

Earthquake Clusters (Objective 3)
The earthquakes in Figure 6 show the 775 filtered events within the subducting slab.Figure 7 illustrates the results of two statistical analyses.The Kernel density (Figure 7a) distinctly delineates the high and low seismicity regions with their general boundaries being aligned with the projected subducted direction of identified PWs and PFS linements.In addition, K-means clustering and least-square linear regressions (Figure 7b) were employed to examine the statistical orientations of the earthquake clusters.The best fit lines for the majority of the earthquake clusters align with the projected subduction direction of the PSF lineaments and PWs.However, it is noteworthy that several trends in the northern part of the margin deviate from the subduction direction and instead coincide with the boundary of the Siletz terrane.We address this result in the next section.

Discussion
Our three 2-D subsurface models (Figure 3) reveal variations in the oceanic crust density.Previous density distributions for the JdF oceanic plate, such as the study of Romanyuk et al. (1998), lacked seismic constraints and resulted in rather questionable values in some regions, such as highly varying mantle densities (from 3.22 to 3.33 g/cm 3 ) without any supporting evidence.The 2-D plate-scale gravity model previously developed by Blakely et al. (2005) utilized a single density value for each tectonic element, such as the entire oceanic crust having a single density of 2.85 g/cm 3 .In our models, crustal densities gradually increase from 2.8 g/cm 3 at the spreading center to 2.9 g/cm 3 near the subduction zone (Figure 3) to account for age-related density variations within oceanic crust.Notably, the crustal density of Blakely et al. ( 2005) is an average value of our density range.It is also important to note that the density of a crustal block in our 2-D models represents the average density of different vertical layers of oceanic crust noted by Carlson and Herrick (1990).Horizontal extents of individual crustal blocks in our southern and northern models (namely Oregon and Washington models illustrated in Figure 3) are constrained with magnetic anomaly data, which provide a temporal subdivision of the JdF oceanic crust.These two models have reliable seismic reflection constraints from Han et al. (2016), allowing us to fix the geometry of the subsurface layers from the seafloor to the mantle.In contrast, the central model is based on a single channel vintage seismic data from Carpenter (1971) that does not image the Moho boundary.As we have to infer the crustal thickness in that model from regional seismic data, we choose not to do magnetic modeling for that line.Instead, we limit our analysis to gravity only to prove the presence of zones with low crustal densities and to determine their locations to constrain our spatial analysis.
Despite the increasing age-dependent density trend of oceanic crust in our 2-D models, several zones of lower crustal density are mandatory to fit gravity data (as illustrated in Figure 3).We explain these lower density zones by highlighting their alignment with known PWs and explaining their presence with abundant faulting and fracturing within PWs.The bookshelf-like faulting that occurs within a plate boundary zone allows for the movement of shear stress between propagating and failing spreading centers (Hey, 2021).This suggests that there is a greater occurrence of faulting within the PW, which potentially is the primary factor behind the observed decrease in density compared to the surrounding oceanic crust.The findings of this paper-PWs have lower densities than the surrounding crust-are also consistent with several studies relating the PWs of the JdF plate to observed seismicity clusters (Han et al., 2018;Merrill & Bostock, 2019;Nedimović et al., 2009), thus confirming crustal weaknesses along these structures.PWs around the world also show similar signatures of enhanced crustal fracturing and require lower densities during gravity modeling (Peirce et al., 2001;Wiedicke & Collier, 1993).Furthermore, multiple studies (e.g., Horning et al., 2016;McClymont & Clowes, 2005;Nedimović et al., 2008;Newman et al., 2011;Rathnayaka & Gao, 2017;Weekly et al., 2014) show that PWs in the JdF plate exhibit lower seismic velocities than the surrounding crust, supporting the interpretation that these regions represent zones of enhanced faulting and thus weaker crust.In addition to the large-scale disturbance in seafloor magnetic strips (i.e., PWs), we have also identified smaller scale offsets in sea-floor magnetic stripes that are mapped as PSF lineaments (Figure 4).This series of linear structures have certain spatial characteristics, namely (a) they mostly start at the breaks in the spreading centers (i.e., in between the ridge segments, and are roughly aligned with the subduction direction), and (b) they coincide with identified low-density zones from gravity analysis (Figures 3  and 4).
Figures 5 and 6 depict the seamounts mapped over the JdF plate through our spatial analysis.These seamounts have received limited attention in research, resulting in a lack of well-established information regarding their age and lithology in the existing literature.As a result, the origin of these seamounts remains a subject of ongoing debate.However, several studies indicate a potential association between PWs and seamounts, as they might originate from the same locations.Sreejith et al. (2016) noticed a potential correlation between mapped seamounts from gravity data in the Arabian Sea and their neighboring PWs.Ridge jumps and propagation in the South Atlantic Ocean may also be related to the formation of seamounts (Brozena & White, 1990).In the first seismic reflection study across a PW of the JdF plate, Calvert et al. (1990) suggested that seamounts in the northeast Pacific potentially formed as a result of rift propagation, which ties these seamounts to the same original location as the mapped PWs and PSFs of this study.In addition, a seismic reflection study conducted by Harris et al. (2020) revealed the presence of a buried seamount located above the PW "iii" region to the west of the deformation front offshore Washington.Similarly, Han et al. (2018) reported the identification of a buried seamount near the PW "iv" area in our study, situated just west of the deformation front offshore Oregon, through seismic reflection imaging.Moreover, the seamount formation requires significant magmatic supply; they most likely originated from the JdF spreading centers and Blanco transform fault zone.Since our mapped PSF lineaments and known PWs represent zones of weakness, it is natural for this magmatic supply to take advantage of these structures.
To evaluate this correlation, we count the number of seamounts within PWs and that are within a certain distance from our mapped PSF lineament.We utilize the ArcGIS geoprocessing tools, such as Buffer and Intersect, to carry out this analysis.For the distance parameter, we used 15 km, 10 km, 5 km and 3 km from the mapped lineaments (inset graph of Figure 6).All the percentage values reported in this section are rounded up to the nearest whole number.For distances of 15 and 10 km, respectively, all or part of the gravity and seismic seamounts fall within the distance limit in 100% and 88% of cases.Specifically, for a 15 km distance, 77% of gravity and seismic seamounts are completely within the limit, whereas for a 10 km distance, this figure is 49%.As for bathymetric seamounts, 95% of them are at least partially within both the 10 and 15 km distance limits.Among these, 76% of bathymetric seamounts fall entirely within the 15 km limit, while for the 10 km distance, this percentage is 59%.Even with stricter distance thresholds of 5 and 3 km, a significant proportion of gravity and seismic seamounts, namely 70% and 68%, respectively, are found within these limits, either fully or partially.Similarly, 65% and 60% of bathymetric seamounts are located within 5 and 3 km respectively of PWs and PSFs.
We interpret our mapped linear structures (PWs and PSFs) as zones of crustal weakness that formed during the propagation of spreading centers.Notably, some of these lineaments originate at the disruptions in the Blanco transform fault and are associated with reorganization between competing fault segments.Extensive faulting during rift propagation and transform fault reorganization resulted in reduced crustal densities, and consequently, has affected the crustal strength.The seamounts probably take advantage of these structures, which results in the potential clustering along these features, as documented in our study.We do not address the relative age of the seamounts and the source of magmatic material here; we simply document the fact of the formation of the seamounts along these zones of crustal weakness on the preexisting oceanic crust.
The crustal weakness zones, along with seamounts, predominantly extend toward the Washington margin (Figure 6).Our kernel density estimation analysis of onshore earthquakes within the subducting slab reveals a concentrated region of high earthquake density in Washington (Figure 7a).This region lies within the projected subduction direction between PSF lineament 9 and PSF lineament 11, or the PW "iv."This area is notable as at least four mapped oceanic structures are undergoing subduction.Conversely, one of the areas with a lower occurrence of earthquakes is observed between PW "iv" and PSF 12, where no mapped oceanic structures are traced during the subduction process.Moreover, the only high earthquake density region in Oregon can be spatially related to the projected subduction direction of the mapped PSF lineaments 12 and 13.Additionally, the few earthquake clusters in offshore central Oregon are around buried seamounts (Morton et al., 2018;Tréhu et al., 2012).This correlation of high earthquake density with the mapped oceanic structures, which are found to be a weaker oceanic crust in this study, suggests that these structures may have a significant contribution toward the broad seismicity pattern of CSZ.Our K-means clustering analysis of these subducting slab earthquakes also corroborates this relationship (Figure 7b).We perform regression analysis throughout these clusters to understand their general trend and the direction along which most earthquakes lie.In Washington, the best fit lines through groups 1 to 4 exhibit a spatial correlation with the Siletz terrane.This suggests that the terrane boundary has a higher degree of influence on the seismicity of northernmost Washington, which supports the findings of Merrill et al. (2020).Best fit lines through other clusters (groups 5 to 10) in southern Washington and Oregon do not show this stronger influence of the Siletz terrane boundary as they are more aligned to the projected subducted direction of mapped PSFs and PWs.These observations suggest that once these regions of weaker crust enter the subduction domain (Figure 8a), they are more susceptible to crustal rupture as earthquakes than the adjacent stronger oceanic crust.In contrast, fewer of the weaker crustal features subduct beneath the Oregon portion of the margin (Figure 8b), resulting in notably fewer earthquakes in that part of the CSZ.
A similar pattern can be observed for both the interplate and upper plate seismicity.It is not exactly clear how the structures on the subducting slab influence the upper plate.However, the presence of water content within the weaker portion of the subducting slab, such as the mapped PWs, PSFs, and seamounts, can impact the rupture behavior of both the upper crustal faults and the megathrust.Canales et al. (2017) showed that water is mostly concentrated in the upper part of the JdF subducting slab, leaving most of the crust dry.They report that two PWs of the JdF plate (PW "iii" and "iv" in our study) and a buried seamount that is located just west of the deformation front and near the PW "iv" (first noted in Han et al. (2018)) contribute significantly to the pore water content of the CSZ.In the northern CSZ, a zone of high Poisson's ratio (which is also electrically conductive) is noticed that may be caused by slab-derived fluids trapped at near-lithostatic pore pressure (Calvert et al., 2020).The location of this fluid rich crust is along the approximate subducted location of PW "ii."Furthermore, a teleseismic study by Abers et al. (2009) noted a layer of weak sediments or an overpressured fault zone in Washington.Because of the episodic tremor and slip signals in that region, they favor the model with near-lithostatic fluid pressure in that zone.These studies suggest that our mapped structures, especially PWs and PSFs, introduce hydration into the subduction zone and may have a significant impact on the characteristics of interplate seismicity.For instance, Schmalzle et al. (2014) explain the lack of earthquakes in the central portion of the margin as aseismic creep due to high fluid pressure below an impermeable thick part of the Siletz terrane.How this hydration of the subducting slab affects the upper crustal seismicity is not well established.Furthermore, the upper plate deformation can rely on the geology and the strength of the material (McKenzie et al., 2022).However, in Wells et al. (2017), 10.1029/2023GC010943 14 of 18 they show that low frequency tremor density (reported as highest over northern CSZ) statistically correlates with upper plate faults.They suggest that these faults extend to overpressured megathrusts and provide pathways for fluid escape.
The seamounts that we mapped can also play a significant role in the overall seismicity pattern of the Cascadia margin.In the north-central region of the CSZ, there is compelling evidence of seamounts that exhibit a spatial correlation with small-scale earthquakes (Morton et al., 2018(Morton et al., , 2023;;Tréhu et al., 2012;Trehu et al., 2015), meaning they are helping to release energy underneath the accretionary prism.Furthermore, Canales et al. (2017) suggested that buried seamounts are responsible for a significant portion of the total hydration in the CSZ.Based on the subduction direction, we can estimate that the majority of seamounts (∼80% of the total number, located north of the PW "iv") subduct underneath Washington, while the remaining ∼20%, located to the south of that PW, head toward Oregon.Therefore, our mapped seamounts are probably also responsible for the higher amount of interslab seismicity over the northern part of the Cascadia margin.
Our approximated region of fewer subducted features offshore Oregon (between 44°N to 46°N) suggests a lower slip rate deficit or plate coupling with the exception of the localized energy release from small earthquakes near subducted seamounts around 44.5°N (Morton et al., 2018(Morton et al., , 2023;;Tréhu et al., 2012;Trehu et al., 2015).This observation aligns with previous studies using geodynamic data to model slip rate deficits.

Conclusions
The present study relates the observed non-uniform seismicity along the CSZ to crustal inhomogeneities mapped within the incoming JdF oceanic plate.We identify a series of low-density crustal zones along three regional models that we interpret as evidence of weaker oceanic crust related to extensive faulting.Through the spatial analysis of magnetic data.We map a series of linear features through these weaker crustal regions that include both previously known PWs and newly interpreted pseudofault lineaments (PSFs).The latter are associated with smaller offsets in magnetic anomalies.Both PWs and PSFs represent linear regions of a less dense and weaker oceanic crust, and are generally aligned with the subduction direction.Notably, the buried seamounts that we map from residual gravity and seismic data appear to be clustered around these linear zones of weaker crust.Moreover, these zones are spatially aligned with observed earthquakes that occurred within the subducted JdF slab.We conclude that the higher seismicity recorded in the Washington portion of the margin relates to a greater number of these weaker crustal features in the subducted JdF slab in contrast to notably fewer of these heading toward the Oregon margin.Therefore, our findings explain the complex pattern of a subset of recorded onshore seismicity in the CSZ.
the JdF plate.We use the top of the subducting crust information from McCrory et al. (2012) to filter out the earthquakes that are above the Wadati-Benioff zone.

Figure 2 .
Figure 2. Potential field data used in this study are shown with propagator wakes and 2-D modeling profiles.(a) Total magnetic intensity data from Bankey et al. (2002) and (b) Free air gravity data from Sandwell et al. (2014).

Figure 3 .
Figure 3. 2-D integrated models of Cascadia Subduction Zone from the Juan de Fuca (JdF) spreading center to the continental domain.The Oregon (a) and Washington model (c) end in Oregon and Washington, respectively, while the Central model (b) ends in the boundary between Oregon and Washington.In each model, the top panels represent potential field data and the bottom panels show the subsurface model.Density values in g/cm3 for each modeled block are shown with a light orange color in the subsurface model.Only for "a" and "c," the magnetic susceptibility values in cgs are also shown in the subsurface model with a light blue color.However, magnetic susceptibilities are not shown for every modeled element, especially where they are obviously zero, such as seawater.To account for the magnetic reversals in sea-floor, we use both positive and negative magnetic susceptibility interchangeably while modeling the JdF crust.The red star in the upper panels represents the tie point which shifts the calculated gravity and magnetic anomaly to the datum of the observed data.

Figure 4 .
Figure 4. Mapping of PSFs on the Juan de Fuca (JdF) plate from the disturbances of sea-floor magnetic stripes.Figures show newly interpreted PSFs in the northern (a and b) and southern (c and d) parts of the JdF plate.The interpretation was made from both the reduced to pole total magnetic intensity map (a and c) and in its residual (upward continued by 3 km) tilt derivative (b and d).The tilt derivative is a unitless quantity.The white lines in the maps represent modeled the lower density zones shown in Figure 3.The roman numerals indicate propagator wakes, while the numbers represent PSF lineaments.

Figure 5 .
Figure 5. Mapping of seamounts using the gravity data.(a) Bouguer gravity anomaly and (b) residual Bouguer gravity anomaly map (see Section 2 for methodology) around the vintage Seismic reflection line 10 from cruise expedition RC1501(Carpenter, 1971).The profile for this seismic reflection line is also shown in Figure1.(c) Image of the seismic reflection line 10 and few interpretations of its geological features.The abbreviated text "BT SM" and "BR SM" represent bathymetric and buried seamounts respectively.

Figure 6 .
Figure 6.Diagram illustrating mapped tectonic features of the Juan de Fuca plate (i.e., previously established propagator wakes (PWs) and, newly mapped PSFs and seamounts) and their possible spatial correlation with the Wadati-Benioff zone earthquakes of Cascadia Subduction Zone.The graph in the inset shows the spatial relationship between the mapped linear tectonic features (PWs and PSFs) and seamounts.The bars in the graph show the percentage of seamounts that are within a certain distance of the PWs and PSFs.

Figure 7 .
Figure 7. (a) Kernel density estimation of onshore earthquakes within the subducting slab along Cascadia Subduction Zone, displaying varying densities with a color bar from highest (red) to lowest (blue).The mapped lineaments and oceanic structures, with similar symbology and numbering as Figure 6, are also shown in this figure.Approximate projection lines (blue dashed lines) from mapped lineaments are drawn in parallel with the subduction direction (yellow arrow) to showcase the correlation between oceanic structures and broader seismicity patterns.(b) K-means clustering of the onshore earthquakes (circles of different color) with corresponding best fit lines.Different colors represent distinct groups labeled Gr 1 to Gr 10.

Figure 8 .
Figure 8. Schematic block diagrams showing the lower density zones or weak crustal zones and the associated seamounts are subducting beneath Washington (a) in a higher amount than Oregon (b).Both the Washington and Oregon blocks are based on their respective models in Figure 3.
Schmalzle et al. (2014) found reduced interseismic vertical velocities, decreased slow slip earthquakes and significantly less locking fraction (in the Gaussian model) between 44°N and 46°N.Stress constrained geodetic modeling from Lindsey et al. (2021) and viscoelastic modeling from Li et al. (2018) also show less coupling in this region.