Deformation of the Juan de Fuca Plate Beneath the Central Cascadia Continental Margin (44°-45°N) in Response to an Upper Plate Load

Introduction

Many investigators have attempted to estimate the strength of the subducting plate based on the flexural response of the plate to subduction as measured from topography as the plate descends into the trench (e.g., Caldwell and Turcotte, 1979; Forsyth, 1980; Bry and White, 2007; Contreras-Reyes and Osses, 2010; Craig and Copley, 2014). An expected relationship between plate age and strength at the trench, however, has been elusive. Bry and White (2007) discerned a possible linear increase in apparent elastic thickness from 5–10 km for ages of 0–20 Ma from a global survey of subduction zones, whereas Contreras-Reyes and Osses (2010) argued for an elastic thickness of 13–15 km for a plate age of ∼5 Ma offshore Chile. Craig and Copley (2014) found no consistent relationship between plate age and strength and attributed this to bending beyond the elastic limit, which weakens the plate and results in outer rise seismicity (Forsyth, 1980).

Other studies have highlighted the impact of upper plate structure on subduction zone earthquakes (e.g., Song and Simons, 2003; Wells et al., 2003; Collot et al., 2004; Lotto et al., 2017; Salleres and Ranero, 2019; Egbert et al., 2022). Recently Arnulf et al. (2022) and Kimura et al. (2022) showed that the presence of a large mafic intrusion in the forearc of the Nankai Tough beneath the Kii peninsula deflects the subducting plate. Arnulf et al. (2022) argued that the high-density intrusion, which has an area of ∼8,000 km², influences upper mantle hydration, earthquake nucleation, and plate boundary segmentation. Kimura et al. (2022) argue that it acts as a barrier to rupture propagation, with large plate boundary earthquakes nucleating near the edges of the intrusion and propagating away from it. We discuss a new three-dimensional (3D) model of the P-wave velocity (Vp) structure of the central Cascadia margin from 44° to 45°N and suggest that similar processes may be occurring on the central Cascadia margin.

The model was constructed from amphibious, controlled-source seismic data obtained through several piggyback programs over the past three decades and includes a high velocity (and by inference density) anomaly in the upper plate that deflects the lower plate. In contrast to the Kii peninsula region of the Nankai Trough, in which an ∼100 km wide pluton loads a 27 million-year-old plate that is subducting at a rate of ∼47 mm/yr, the central Cascadia upper plate load is a narrow (<10 km) vertical slab that loads a 10 million-year-old plate (Wilson, 2002) subducting at a rate of ∼31 mm/yr (Wilson and McCrory, 2022). We use the geometry of the lower plate in response to the load to estimate plate strength. The relict nature and compact geometry of the upper plate load and its direct correlation with a deflection in the Moho of the subducted plate provide a unique opportunity to derive an in situ estimate of the elastic strength of the subducted plate. We discuss the results in the context of models for the thermal structure of the plate (Cozzens and Spinelli, 2012), contemporary interplate locking (DeSanto et al., 2022) and the extent of sediment subduction derived from a 3D electrical resistivity model (Egbert et al., 2022). Because the frictional behavior of a fault and the amount of strain energy that can be stored in the adjacent rock depends on the balance between the strength of the fault surface and the strength of the adjacent rock (Cook, 1981), our model provides new insights into the geodynamic response of the Cascadia megathrust.

Background

The Cascadia Subduction zone is exceptional in its low level of instrumentally recorded seismicity, although the presence of an active volcanic arc and historical and paleoseismic evidence for large plate boundary earthquakes indicate that it remains an active subduction zone and that large earthquakes are likely in the future (Walton, Staisch, et al., 2021). Geologic and historic tsunami data indicate that the most recent subduction-zone-wide great earthquake occurred in 1700 (Atwater et al., 2005), and paleoseismic data indicate an average, but highly variable, recurrence interval of ∼500 years for such events. Figure 1A shows that the offshore forearc between ∼43° and 47°N has been nearly devoid of earthquakes with M >3 since at least 1989, except between 44° and 45°N (Figure 1A).

FIGURE 1. (A) Central Cascadia subduction zone showing earthquakes with magnitude greater than 3 in the ANSS Comprehensive Catalog from June 1989 through April 2023 (red dots). Violet lines offshore show locations of 3 strike-slip faults from Goldfinger et al. (1997): AC, Alvin Canyon fault; DB, Daisy Bank fault; WF, Wecoma fault. PO shows location of Pythia’s Oasis (Philip et al., 2023), an unusually warm and vigorous fluid seep. Black dashed lines offshore show the outline of basins from Wells et al. (2003): WB, Willapa Basin; AB, Astoria Basin; NB, Newport Basin bifurcated by a structural high known as the Stonewall Bank (SB), a NW-trending anticline with Miocene sediments exposed at the seafloor that has been active in the Holocene (Yeats et al., 1998); CB, Coos Bay Basin. White dashed line in (A) and (B) shows the seaward edge of the early Eocene Siletz terrane (Wells et al., 1998). Dashed lines onshore show late Eocene volcanics that penetrated through the underlying Siletz terrane (Wells et al., 2014): YB, Yachats Basalt; CH, Cascade Head; TV, Tillamook Volcanics; GH, Greys River Volcanics. White rectangle shows location of inset. Black rectangle shows location of (B). (B) Aeromagnetic map of the study region with ANSS-catalog earthquakes relocated by Tréhu et al. (2008, 2015, 2018) and Williams et al. (2011) in red. The depth distribution of these events, relocated relative to the M4.7 and M4.8 earthquakes, are shown in the depth section. Additional earthquakes located by Stone et al. (2018) and Morton et al. (2018) in dark and light blue, respectively. Depths of these earthquakes are scattered and not reliable (Tréhu et al., in prep). Regional moment tensor solutions are shown for the four largest events since 2004 (Tréhu et al., 2008, 2015, 2018). Black dashed line shows the extent of the magnetic anomaly indicative of the Yachats basalts. Grey dashed circles labeled “smt” show the location of subducted seamounts inferred from the magnetic anomalies (Tréhu et al., 2012). Cross-section shows the major structural boundaries interpreted from Gerdom et al. (2000) and relocated earthquakes. JdF, Juan de Fuca plate. Moho (solid grey lie) and the contact between the accretionary complex, Newport Basin and Siletz terrane (long grey dashed line) in the upper plate are well constrained by the Vp model; the plate boundary (short dashed grey line) is inferred based on the thickness of oceanic crust prior to subduction.

In 2004, two moderate earthquakes (M4.9, M4.7), occurred approximately 1 month apart and were followed by several aftershocks (Figure 1B). Based on regional moment tensor (RMT) analysis, high frequency pP observations and raytracing of secondary phases through a 2D velocity model of the source region (Gerdom et al., 2000), Tréhu et al. (2008) argued that these events were low-angle thrusts located on or near the plate boundary with slip vectors approximately aligned with the predicted plate motion vector (Figure 1B). Since 2004, smaller earthquakes have been reported by the PNSN (Pacific Northwest Seismic Network) at a rate of ∼6/year in clusters around these larger events (e.g., Tréhu et al., 2015; Tréhu et al., 2018), and inclusion of ocean bottom seismometer data from the Cascadia Initiative (Toomey et al., 2014), which provided offshore data from 2011 and 2015, has resulted in identification of many additional earthquakes in these clusters (Morton et al., 2018; 2023; Stone et al., 2018). Locations of these events using velocity models that approximate the structure beneath the continental shelf rather than the simpler models used by the PNSN suggest that hypocenters in the PNSN catalog are systematically overestimated, and relative relocation of clusters of earthquakes suggest that many occur in a narrow depth range near the plate boundary (Williams et al., 2011; Tréhu et al., 2015).

The relationship between the seismicity and forearc crustal structure is illustrated in Figure 1B based on potential field data calibrated by 2D seismic transects (Tréhu et al., 2012). The seaward edge of the Siletz terrane, an oceanic plateau that was accreted to North America at ∼50 Ma (see Wells et al., 2014 for a comprehensive review), is clearly defined by high-resolution aeromagnetic data (US Geological Survey, 2002), as are several subducted seamounts. A cluster of earthquakes is associated with a subducted seamount that appears to be impinging on the seaward edge of the Siletz terrane (Figure 1B). RMT solutions for 2 M >3.8 earthquakes in this cluster that occurred in 2012 and 2017, show thrust mechanisms that are steeper than those of the 2004 earthquake and may reflect the impact of a subducted seamount on the leading edge of Siletzia (Tréhu et al., 2015; Tréhu et al., 2018). A second cluster is located several kilometers landward of the seaward edge of the Siletz terrane and coincident with a local gravity low that reflects an anomalously deep part of the Newport basin (Tréhu et al., 2012). It was interpreted to indicate erosion of the Siletz terrane from below, resulting in the relatively recent formation of a basin within the larger and older Newport basin (Figure 1A), which was initiated in the late Miocene (McNeill et al., 2000).

The distinctive northwest-trending lineations in the aeromagnetic data (outlined and labeled YB in Figure 1B) coincide with coastal exposures of the Yachats Basalt (Figure 1A insert) and extend ∼17 km offshore (Figure 1B). The Yachats Basalt consists of seaward-dipping basalt flows about 37 m.y. old that were erupted from a west-northwest-oriented dike swarm intruding the accreted Siletz terrane onshore (Snavely and MacLeod, 1974; Snavely et al., 1976; Wells et al., 2014; Wells and Blakely, 2019). Wells et al. (2014) have argued that the Yachats Basalt, as well as several other middle-to-late Eocene volcanic units (Figure 1A), formed in a margin parallel extensional stress field and were subsequently rotated approximately 50° clockwise along with the entire Oregon block of the Siletz terrane. The linear topographic ridges and magnetic anomalies formed by the Yachats Basalt dikes are unique along the Oregon and Washington coastline (Figure 1A), and they project toward the anomaly offshore (Wells and Blakely, 2019).

The forearc between 44° and 45°N is also characterized by several other along-strike changes in geologic and geophysical parameters interpreted as potential indicators of interplate dynamics (Walton et al., 2021). These include: 1) a trio of strike-slip faults that offset the deformation front and can be traced across the outer wedge (Figure 1A; Goldfinger et al., 1997); 2) a lower degree of locking than north of 46°N or south of 43.5°N (Schmalzle et al., 2014); 3) a low slip region between two high slip patches in the 1700 earthquake at 44.4°N (Wang et al., 2013); 4) an offshore gravity high interpreted to be a potential rupture boundary (Wells et al., 2003); 5) segment boundaries at 44.2 and 45°N in the paleo-earthquake record derived from marine turbidites (Goldfinger et al., 2012); and 6) an abrupt southward increase in the slope of the outer wedge at 45°N (Watt and Brothers, 2020).

Data Acquisition

The Cascadia subduction zone between 44° and 45°N has been the focus of several controlled source seismic experiments in the past several decades (Figure 2). The primary dataset for this study was acquired as a piggyback project on the 2012 Ridge-to-Trench project (e.g., Han et al., 2016; Canales et al., 2017; Han et al., 2018). This included, in addition to the primary Ridge-to-Trench data, four coast-parallel lines of airgun shots with shorter connecting lines that were located to complement 2 coast-parallel shot lines acquired in 1996. The airgun sources from the R/V Marcus Langseth (cruise MGL1211) were recorded on seven additional ocean bottom seismometers (OBSs) and 35 temporary short period, continuously-recording seismometers in the Coast Range (Figure 2). No seismometers were deployed in the Willamette Valley because seismometers in the Willamette Valley installed in 1996 had background noise levels in the frequency range of interest that were ∼5 times higher than those in the Coast Range.

FIGURE 2. Map showing location of sources (lines) and receivers (inverted triangles) from the 1989 (orange), 1996 (yellow) and 2012 (red) data acquisition programs that were used for this analysis. The three sites and lines of sources for which data are shown in Figure 3 are highlighted and labeled. Source line T07 is also labeled because it is mentioned in the text. The white rectangle shows the limits of Figure 1B. The dashed blue line shows the location of CASIE21 line PS01.

Data from two prior experiments are also included in the dataset for the model. In 1989, as a piggy-back on an ODP seismic reflection site survey of the deformation front (MacKay, 1995), we acquired a single multichannel seismic reflection (MCS) transect from the Juan de Fuca plate, across the deformation front, to the coast near 44.8°N (Tréhu et al., 1995). The seismic sources for the MCS profile were recorded on a linear array of eight ocean bottom and ten onshore seismometers, which detected energy from diving waves and wide-angle reflections to offsets of up to 160 km (USGS Open-File reports 93-317, 93-318). The data were modeled using MacRay, a 2D forward modeling approach based on user-specified layers and lateral velocity variations with the layers (Tréhu et al., 1994). The result showed that thickened Siletz terrane rocks extended 30 km offshore, forming a nearly vertical backstop to the accretionary prism and restricting sediment subduction to at most a 2 km thick channel.

In 1996, a second amphibious, large-aperture transect was acquired near 44.6°N (Figure 1B). Airgun shots were recorded on 21 ocean bottom seismometers or hydrophones and on 29 onshore stations on a line that extended from the coast to the Cascade foothills (Gerdom et al., 2000). The 1996 experiment also included two north-south profiles on the continental margin (Gerdom et al., 2000), which were used to develop velocity/depth models for locating offshore earthquakes in this region (Williams et al., 2011; Morton et al., 2018; Stone et al., 2018). These margin-parallel shot lines were recorded as “fan” shots on the transect and on five additional seismometers that had been installed north and south of the primary transect to inform on 3D crustal structure (Figure 2). As in the 1989 data set, strong secondary arrivals interpreted to be reflections from the Moho of the subducted plate were observed, indicating that PmP observations could be used to map the Moho of the subducted plate in 3D with a suitable distribution of sources and receivers (as had been done in 1998 beneath the Olympic Peninsula; Tréhu et al., 2002; Preston et al., 2004), providing the impetus for the 2012 onshore experiment. A simplified version of the 2-dimensional Vp model derived from this experiment (Gerdom et al., 2000) served as a starting model for this study. A subset of the 1996 dataset was incorporated into the database for this study (Figure 2).

A much larger network of 755 vertical component nodal seismometers was installed in 2021 with an interstation spacing similar to that in the prior experiments but with a much larger footprint that extended from the California/Oregon border into southwest Washington (Cascadia 2021 Field Team, 2022). The nodal seismometers recorded shots from the 2021 CASIE21 multichannel seismic (Carbotte et al., 2021) and OBS (Canales et al., 2023) experiment and will enable extension of 3D velocity modeling to the margin from the Oregon/California border to the Olympic peninsula. Because only one margin-parallel line of shots was acquired along the forearc and many sites along this segment of the margin were reoccupations sites installed in 2012, the new large aperture data were not used in this study since they provide limited additional information on the structures of interest for this study. The CASIE21 reflection data along line PD01 (Figure 2), however, are important for validating our assumption that the Moho of the subducting plate is a proxy for the position of the plate boundary.

Data Processing and Interpretation

The integrated data set from 1989, 1996, and 2012 used to develop the velocity model presented here includes 57 onshore and 35 ocean bottom recorders. Station coordinates and elevations are given in Table 1 of Kenyon (2016). The continuous data recorded on the onshore and ocean bottom seismometers were merged with the shot instant data to generate SEGY-format data and record sections for each station/line pair. Record sections were examined with different bandpass filters and reduction velocities. Results from 2D transects were used to identify the offset range at which to expect high amplitude secondary arrivals representing reflections from the base of the subducted crust and potentially other distinctive secondary arrivals on fan shots, defined as source to receiver geometries in which the shots and receiver are not aligned. To illustrate the magnitude of this task, we note that each of the stations installed in 2012 resulted in 18 record sections from the vertical component alone; in total, ∼1,000 record sections were considered, with the majority of these showing signals from the offshore shots.

P-wave arrival times were picked from the vertical component of the onshore stations. Either the vertical geophone or hydrophone was used to pick arrivals on ocean bottom seismometers, depending on which component had the greater signal-to-noise ratio. Although some of the stations included horizontal geophones, no horizontal component data were used for this analysis. Pg, PmP, and Pn phases, defined as diving and refracted waves through the crust, wide-angle reflections from the base of the crust, and refracted wave in the uppermost mantle, were picked manually using tlPicker software developed at the University of Washington. Given the large amount of data, we were conservative and picked only arrivals for which the phase identification was unambiguous, resulting in 53,855 Pg, 17,922 PmP, and 13,143 Pn arrival time picks.

Three examples of record sections showing data from a single line of shots recorded at a single receiver are shown in Figure 3, with the location of the shot lines and stations shown in Figure 2. These examples are typical of the data and were chosen to illustrate aspects of the data that control the features of the model discussed in this paper. Figure 3A shows a rapid travel time advance in Pg between offsets of 70 and 80 km, suggesting that crustal P-waves from sources south of this transition are traveling through material with a much faster velocity than waves from sources to the north and that the crust contains a strong horizontal velocity gradient. This record section also displays strong PmP arrivals indicative of reflections from the base of the crust with angles of incidence close to the critical angle. The source-receiver distance range in which PmP is a strong arrival depends on the depth to Moho and the velocity structure and is a∼57–110 km in Figure 3A.

FIGURE 3. Examples of vertical component data recorded onshore that illustrate the impact of 3D crustal structure on the data. The data are displayed as equally spaced traces (an approximately linear function of distance along the shot line given) with the X-axis labeled by the distance between source and receiver, which is a non-linear function. Travel time picks for first arrivals interpreted to be Pg and secondary arrival interpreted to be PmP are overlain on the data. (A) Station S04, Line T03. (B) Station S16, Line T05. (C) Station 21, Line T05.

In Figure 3B, the shortest distance between the line of sources and the receiver is 49 km, near the middle of the record section. A clear difference in the travel time for sources from the north compared to shots from the south is observed, and the shortest time does not correspond to the shortest source/receiver distance. Water depths are nearly constant along most of the line (Figure 2) and cannot explain this asymmetry, which must result from along-strike variations in upper plate velocity. Strong secondary arrivals consistent with PmP are also observed for source/receiver offsets >50 km. No picks were made on the southern end of this section because of uncertainty about how to classify these arrivals.

Figure 3C shows arrivals from the same shots as in Figure 3B recorded on a seismometer located only ∼15 km farther south. As in Figure 3B, for the same source/receiver distance, ray paths from the north are much faster than ray paths from the south. PmP is observed along the entire line, and the time difference between Pg and PmP for a particular source/station distance varies with the azimuth as well as the distance between the source and receiver. Note also strong scattering of Pg on Figure 3B, recorded on Yachats Basalt, compared to Figure 3C.

Inversions were run using Stingray/Tomolab software developed at the University of Oregon (Toomey et al., 1994; Dunn et al., 2005). Crustal Pg arrivals were inverted with a variety of starting models, including a 1D velocity versus depth function fixed to the seafloor surface and an east-west oriented 2D velocity model representing a smoothed version of the model of Gerdom et al. (2000). Within the volume that is well sampled by the ray paths provided by our experimental geometry these two starting models yield similar results, with a slightly higher final misfit for the 1D model (see Figure 14 in Kenyon, 2016). Because of the overall eastward thickening of the crust, however, the 2D starting model results in fewer artifacts around the edges of the well sampled volume and produces a model that can more easily be used for studies for which we must extrapolate from the well-resolved volume (e.g., for earthquake relocations or for depth migration of multichannel seismic reflection data). Different smoothing and regularization approaches were also tested and compared as well as synthetic exercises to test the resolution provided by the experimental geometry. Examples of resolution tests are given in Supplementary Figure S1, with additional tests in Kenyon (2016).

Once the preferred crustal Vp model based on first arrivals had been obtained, PmP and Pn arrival times were included in the arrival time dataset for the inversion, which included a surface corresponding to the base of the subducted crust. Because Pn arrivals were primarily observed from shots located seaward of the deformation, they constrain the upper mantle velocity west of ∼124.7°W. Beneath the continental shelf and Coast Range, the position of the Moho surface is based almost entirely on PmP reflections. The sensitivity of the Moho geometry to different starting models and smoothing parameters is discussed in the next section.

Results

The inversion of travel times converged rapidly and smoothly with most of the decrease in the residuals obtained after 3 iterations and no significant improvement between iterations 5 and 6. The root mean squared misfit between the picked and predicted travel times for the initial 2D starting model was 1.1, 0.73, and 0.36 for Pg, PmP, and Pn, respectively, and decreased to 0.058, 0.072 and 0.04, respectively, after 6 iterations (close to within the estimated pick uncertainty of 0.04 s for Pg and Pn and 0.08 s for PmP). The residual for each iteration and the residuals for individual observations for the sixth iteration as a function of offset are shown in Supplementary Figure S2. Figure 4 shows horizontal and vertical slices spaced at 2 km intervals in depth and at 0.1° (∼8 km) intervals in longitude.

FIGURE 4. Selected equally spaced slices from the Vp model. Left side shows horizontal slices referenced to the seafloor. Right side shows vertical slices along lines of longitude. The thick grey line shows where the Moho surface intersects the slice. Bright green points highlighted by horizontal arrows on the N-S slice along −124.5° are top-of-crust and Moho picks from CASIE21 line PS01. Arrows on slices at 8 km depth and at −124.3° show the high-Vp, near-vertical, NW-trending slab discussed in the text.

The most striking and well-resolved feature of the Vp model is a northwest trending, near vertical, high Vp anomaly in the upper plate. This feature of the model is clearest in horizontal slices at depths of 6–10 km below the seafloor and is shown to extend through the crust in the north-south slices. It is highlighted by arrows in the horizontal slice at 8 km and in the vertical slice at −124.3°. Resolution tests (Supplementary Figures S1A, B) indicate that the shape and internal velocity of this body are well resolved by the experiment geometry, although we are not able to resolve whether its base is in contact with the crust of the subducting plate (possibly even penetrating the upper crust of the subducted plate) or whether a thin, but unresolved, low-velocity zone separates the two. Another well-resolved large anomaly is an oval, low-velocity region located north of the northwest trending high velocity anomaly that extends from the surface to 8 km depth.

The primary feature of the Moho surface model (the thick grey line on the NS slices in Figure 4 and the surface shown in Figure 5A) is an apparent downward deflection of the Moho beneath the high-density crustal slab. Because of the experimental geometry, resolution is best east of −124.7°, which represents the approximate position of reflection points for PmP arrivals from line T7 into stations near the coast. Because of the slab dip, reflection points for PmP are shifted towards the source rather than being midway between source-receiver pairs.

FIGURE 5. (A) Result of inversion of PmP arrivals for the Moho surface corresponding to the model of Figure 4 (model V2). Red line is the coast. Black star shows the apex of the deflection. (B) Cross-section through three 2D starting models used to test the resolution of the Moho surface. (C) Profiles of the Moho surface after six iterations along longitude −124.3°. The result for starting model V2 is shown with two different values of smoothing. The part of the surface sampled by PmP observations is shown as solid lines and the unsampled region is shown by dotted lines. The depth of the starting surface at this longitude is shown by the horizontal dashed lines. Our qualitative estimate of the range of possible values for the amplitude and width of the deflection are also indicated. (D) Profiles of the Moho surface after 6 iterations along latitude 44.35. Starting models V1 and V3 are also shown. (E) Moho surface obtained after six iterations for starting models V3 (i) and V2 with additional smoothing (ii). Red line is the coast; dashed black line indicates the region within which the difference after six iterations for all starting models tested are the same to within 1 km (see Supplementary Figure S4 for difference plots). Discontinuous offsets smaller than ∼1 km may be present but are not resolved by the model.

To test whether this deflection is required by the data, we inverted for the Moho surface for three different starting models. V2 is the starting model that resulted in the surface shown in Figures 4, 5A. V3 was generated by projecting the maximum deflection to the north and south to generate a linear trough rather than a “dimple.” V4 is parallel to V2 but 3 km deeper. Results for the three initial models, including two different smoothing functions for V2, are shown for a north-to-south cross-section at −124.29° in Figure 5C and for a west-to-east cross-section at 44.35° in Figure 5D. The resulting surfaces for V2 with greater smoothing and for V3 are shown in Figure 5E. We conclude that the initial deflection and downdip flattening of the Moho beneath the high velocity slab are robust features of the model that do not depend on the starting model.

Having confirmed that the data require a depression in the Moho beneath the high velocity/density anomaly, we consider how well we can estimate its apex, amplitude and width. The apex of the deflection (shown by a star on Figure 5A) is near 44.35°N, 124.6°W. Estimates for the amplitude and width are correlated, with a width of ∼70 km and amplitude of ∼2 km for the solution for starting models V3 and V4 and a width of ∼120 km and amplitude of ∼4 km for starting model V2 (Figure 5C). Although estimating the width and amplitude of the deflection in the down-dip direction is complicated by an overall dip of ∼14° to the east and limited spatial coverage in the down-dip direction, we estimate a width of ∼50 km and amplitude of 2–3 km, similar in magnitude to the estimate in the strike direction.

Discussion

Geologic Interpretation of the Seismic Model

Figure 6A shows the velocity model at 8 km depth overlain by selected other geological and geophysical features of the continental margin in this region. Figure 6B shows a cross-section through the model that samples the near-vertical, high-velocity anomaly and shows that the contrast in velocity between this feature and the adjacent rocks of the Siletz terrane extends through the upper plate wedge (Figure 6C). The average velocity above the plate boundary (defined as a surface parallel to and 6 km above the Moho surface of Figure 5A for reasons discussed in this section) is shown in Figure 7, highlighting the magnitude and focused nature of the velocity contrast.

FIGURE 6. (A) Velocity slice at 8 km depth with geological and geophysical features discussed in the text overlain. The locations of the vertical slices in (B) and (C) are shown. Thick grey lines labeled in km are contours of depth to the plate boundary for model V2 of in Figure 5 assuming the plate boundary is 6 km above the Moho. Bars above the map show the transition zone from 90 or 100% interplate locking to the zone of episodic tremor and slip (ETS) beneath the Coast Range (from DeSanto et al., 2022). ACF, Alvin Canyon fault; DBF, Daisy Bank Fault; PO, Pythia’s Oasis. (B) Vertical slice through the high Vp slab. Arrows along the top of the cross-section show the seaward edge of Siletz based on magnetic anomaly data and the coastline, respectively. Thin dashed grey line is the interpreted top of Siletzia and a subducted seamount located near the line of this section and associated with seismicity. The outline of the high velocity slab and the Moho and top of crust for model V2 are shown by thick grey lines and labeled. (C) Vertical slice ∼30 km north of (B). Annotations the same as for (B). (D) Temperature 3 km beneath the top of the subducted crust from the model of Cozzens and Spinelli (2012). The red line shows results for a purely conductive model; the blue line shows the high flow (coolest) model of Cozzens and Spinelli (2012).

FIGURE 7. Average velocity above the plate boundary, defined as a surface 6 km above the Moho surface of Figure 5A. Star is the apex of the deflection. Bold dashed line is the region within which the Moho surface inversion is independent of the starting model. Light dashed line is the line along which estimates for the deflection amplitude and width were made. Red line is the coast. Contours are at 0.5 km/s intervals.

The location and orientation with respect to the Yachats Basalt and its offshore extension, as indicated by the outline of their magnetic signature (Figure 6A), suggest that the high-velocity, near vertical anomaly in the upper plate represents the dike complex that fed the near-surface basalt flows, consistent with geologic observations (Snavely and MacLeod, 1974). Geochronologic and paleomagnetic studies indicate that the flows were emplaced within and through the Siletz terrane between 36 and 37.5 million years ago in an extensional forearc environment and have rotated clockwise with Siletz Oregon Block since then (Wells et al., 2014). Since higher Vp implies higher density, the dike complex represents a static load embedded within the upper plate.

The low velocity zone in the upper 8 km north of the dike complex is interpreted to be a basin within the larger Newport Basin and corresponds to topography on a mid-Miocene unconformity mapped throughout the region by McNeill et al. (2000). The model confirms that the earthquakes in the northern cluster are located on or near the plate boundary consistent with prior suggestions that the basin is the results of deep subduction of a seamount beneath the Siletz terrane and erosion of the base of the overlying plate (Wells et al., 2003; Tréhu et al., 2012).

The seaward edge of the Siletz terrane (SES) appears as a north-south ridge of relatively high velocity near −124.5° that clearly bounds the seaward edge of the basin in the slices at 6 and 8 km depth (Figure 4). The cross-sections in Figures 6B, C show an interpretation of the upper and western edges of the Siletz terrane based on the 3D seismic velocity model. This interpretation takes into account the depth-dependent increase in velocity within the terrane (e.g., Tréhu et al., 1994). In Figure 6B, the seaward edge appears to be nearly vertical. The apparent seamount inferred from aeromagnetic anomaly data and associated with the southern cluster of seismicity is not well resolved in this model, although the model suggests that the seismicity is located primarily within the upper plate rather than on the plate boundary, consistent with the conclusion of Tréhu et al. (2012); Tréhu et al. (2015) and Morton et al. (2018) that they result from the interaction between a subducted seamount and the SES. In contrast to the near-vertical SES in Figure 6B, the small velocity inversion along that boundary in Figure 6C suggests the presence of a poorly resolved channel of underthrust sediment (see Supplementary Figure S1).

Figures 6B, C also show the position of the Moho surface from Figure 5A. Based on the Juan de Fuca plate crustal thickness from seismic refraction data from the incoming plate (Tréhu et al., 1994; Gerdom et al., 2000; Canales et al., 2017) and on the multichannel seismic reflection data along line PS01 acquired during cruise MGL2104 (Carbotte et al., 2021) we interpret the plate boundary to be surface located 6 km above the Moho. CASIE21 line PS01 (location shown on Figure 2) shows a strong and nearly continuous reflection, interpreted to be the top of the subducted crust, that overlies an intermittent reflection that is coincident with the Moho surface in the 3D model (Figure 4 on the slice along −124.5°). The agreement between the CASIE21 seismic reflection data and the 3D Vp model validates this assumption. The CASIE21 observations also justify an assumption that the top of the subducted crust and plate boundary in this region are approximately coincident because there is no shallower strong and continuous reflection from the lower crust between 44° and 45°N that would indicate the presence of a layer of subducted sediment between the plate boundary and the top of the subducted crust thick enough to be resolved by our data. The average dip of the subducting plate increases beneath the SES and may flatten toward the coast, and the local deflection beneath the dike complex is superimposed on this regional trend.

We consider it unlikely that the local deflection is a pre-existing feature of the subducting slab that is coincidentally located beneath the crustal slab because we would expect the topography of the overlying plate to sag above such a depression if the depression had not been filled prior to subduction; instead, the Yachats Basalt is associated with topographic highs (Figure 1). On the other hand, if the depression had been filled prior to subduction, we would expect to see a plate boundary reflection above the reflection interpreted to be the top of the subduction crust–a scenario that is not compatible with the CASIE21 reflection data.

Deformation Modeling of the Juan de Fuca Plate

Having ruled out a preexisting depression in the subducted plate, we consider two possible endmember explanations for the deflection: elastic or plastic deformation. If the deformation is purely elastic, the depression should disappear once the slab has subducted beneath the region where the slab is present. If it is plastic, a trough parallel to the convergence vector should remain in the slab at greater depth. The tests of the sensitivity of the Moho surface to the starting model in Figure 5 indicate that the initial depression flattens at ∼30 km depth, consistent with an elastic rebound of a loaded plate. We recognize, however, that this is close to the edge of the resolved model space and cannot resolve whether the rebound is partial (elasto-plastic) or complete (purely elastic).

Assuming that the observed deformation is purely elastic, we can estimate the flexural rigidity and elastic plate thickness of the subducted Juan de Fuca plate by estimating the load resulting from the dike complex and modeling the observed amplitude and wavelength of deflection. We model the deflection as a point load (Brotchie and Silvester, 1969) on an unbroken plate and consider two options for estimating the load. These are depicted in Figure 8A. In the first, the point source is estimated by collapsing the excess load of a column that is 20 km high with a 10 × 10 km cross-section. In the second, the entire mass excess of the dike complex (approximated by a 20 km high prism with a 10 × 40 km cross-section) is collapsed into the point load. This second scenario attempts to qualitatively capture the dynamic effect of a dipping slab moving toward the load. As a third option, we calculate the impact of a line load perpendicular to the deflection (Turcotte and Schubert, 1982).

FIGURE 8. (A) Schematic scenario for flexural modeling. The load was calculated for three scenarios: 1) a point source representing the load of a 10 × 10 × 20 km column; a 10 × 40 × 20 km rectangular slab; and a line source based on a 10 × 20 km cross-section. (B) Average P-wave velocity and corresponding density in the upper 20 km of the model along a north-south transect at −124.29° (light blue line) and a simplified model of a density contrast between the high-velocity dike complex and the adjacent Siletz terrane. (C) Overview of the equations and procedure for calculating the elastic thickness of the subducting plate given a load and the width and amplitude of the depression resulting from the load. Calculations in Table 1 assume E = 80 GPa, n = 0.25 and r of 130, 160 and 190 kg/m³. See Figure 5C for estimates for the amplitude and width of the deflection.

To estimate the density contrast between the slab and the adjacent rock, we convert the average velocity above 20 km depth along a north-south slice of the Vp model at −124.29° and convert Vp to density assuming the Nafe-Drake curve (Ludwig et al., 1970; Figure 8B). This exercise suggests that the average density contrast between the dike complex and the surrounding rock is ∼160 kg/m³.

The equations used for calculating the elastic thickness are shown in Figure 8C. Observables are the width and depth of the observed depression along a north-south profile at longitude −124.29, as shown in Figure 5C, avoiding the eastward increase in dip of the slab. Calculations for the flexural rigidity and elastic plate thickness for the three source scenarios for assumed density contrasts of 130, 160 and 190 kg/m³ and different possible combinations of the width and amplitude of the deflection are given in Table 1 assuming a Young’s modulus of 80 GPa and a Poisson’s ratio of 0.25. Although the resolution tests indicate that observations of the minimum deflection and width are not independent, we included an extreme case of minimum deflection and maximum width to illustrate the sensitivity of the calculation to the observations.

TABLE 1. Summary of results of flexural modeling.

All combinations of parameters yield very low elastic plate thicknesses, ranging from 1.3 to 6.4 km, with a thickness >5 km only achieved with the extreme endmember line-load model of 2 km of deflection over a width of 120 km. We recognize that these simple static elastic calculations are an oversimplification of reality, in which the load is partially supported by the adjacent crust of the upper plate, the motion of the slab beneath the load is time dependent, and viscosity of material both below and above the elastic plate may affect the rates of deflection and rebound. Including these effects as well as the superimposed effect of subduction to the east beneath a thickening forearc wedge of increasing and heterogeneous density is beyond the scope of this study.

These estimates are on the low end of the range of observations reported for the flexural strength of young oceanic lithosphere based on seamount loading of oceanic lithosphere (Watts, 2001). In many cases in nature, estimates of the flexural strength of oceanic lithosphere are complicated by time dependence of the load generation and/or plastic deformation. We have argued that the load in the scenario we model here results from an old magmatic event and therefore is independent of time on a scale relevant to this calculation (although off course, it is susceptible to erosion on a longer time scale) and that the apparent rebound of the deflection argues for a primarily elastic response. Moreover, partial support of the load by the upper plate or delayed deflection because of mantle viscosity would imply an even weaker plate given the observed amplitude and width of the deflection.

The lack of seismicity associated with the lines of maximum curvature of the deflection is consistent with our interpretation of dominantly elastic deformation, although, given the limited length of the seismic record and the ∼1 km resolution of the smoothed Moho surface obtained from the inversion, we cannot rule out the presence of small fault offsets. Our conclusion of dominantly elastic deformation in response to the focused load of the dike complex also does not preclude the presence of larger amplitude elastic and plastic deformation of the subducted plate elsewhere in Cascadia as it encounters large scale variations in the upper plate load.

Comparison to Geodetic and Thermal Models

We now evaluate whether elastic deformation is compatible with geodetic and thermal models. Although geodetic models based only on onshore GPS measurements cannot resolve whether interplate locking is restricted to a narrow band of shore or distributed more widely across a zone of partial locking (e.g., Schmalzle et al., 2014), and increasing database of offshore geodetic measurements supports the former model (DeSanto et al., 2022). Bars above the map in Figure 6A illustrate two models that fit the combined onshore/offshore data from this segment of the Cascadia margin (DeSanto et al., 2022; D. Schmidt, personal communication, 2023). Locking degree is shown as a function of distance from the deformation with black bars representing strong locking and the variable grey shade representing an exponential decrease to 0% locking 165 km east of the deformation front. Along this segment of the margin episodic tremor and slip (ETS), which is thought to result from slip on the nearly unlocked portion of the plate boundary (e.g., Walton, Staisch, et al., 2021) is observed in a band ∼25–65 km from the coast (e.g., Boyarko et al., 2015).

Figure 6D shows predicted temperatures 3 km beneath the plate boundary (i.e., at the midline of the subducted oceanic crust) from the model of Cozzens and Spinelli (2012) for a transect across the Cascadia forearc at 44°N. Results are shown for their endmember models of purely conductive heat transport and compared to crust that has been cooled by vigorous fluid circulation. Although Cozzens and Spinelli (2012) were not able to resolve the impact of fluid circulation on the temperature of the plate boundary beneath the central Cascadia margin from heat flow measurements available at the time, recent heat flow measurements over two buried seamounts located ∼25 km west of the deformation front indicate that vigorous hydrothermal circulation in the upper crust of the Juan de Fuca plate persists to the deformation front (Norvell et al., 2023), suggesting that the in situ temperature is between these two endmembers. The temperature in the subducted crust where the deflection is observed is, therefore, likely well below the brittle/ductile transition for basalt (Violay et al., 2012).

We reconcile the apparently elastic behavior of the subducting plate in this region with the GPS-based estimates of interplate locking, which indicate that the plate boundary where we image the depression in the subducting plate is downdip of the zone of strong locking, by calling on different brittle/ductile transitions for the plate boundary and for the subducted crust. A temperature of 350°C is widely cited as marking the transition between locking and transitional behavior of the plate boundary (Hyndman and Wang, 1993), with 450°C representing the boundary between transitional and freely slipping interplate behavior. These temperatures are thought to represent the behavior of a silica-rich sediment and gouge-filled plate boundary shear zone. In contrast, Violay et al. (2012) have argued, based on laboratory measurements extrapolated to in situ conditions, that the brittle/ductile transition for basalt is ∼550°C. The geodetic data suggest that the Cozzens and Spinelli (2012) endmember of rapid fluid circulation is too cold to be consistent with the narrow band of strong locking indicated by the GPS data and that the conductive model is too hot. Stronger constraints on the temperature structure of the incoming plate are needed to better constrain the amount of cooling provided by fluid circulation (Norvell et al., 2023). Regardless of the amount of cooling due to circulation, however, the subducted crust where the deflection is observed is likely cold enough to be dominated by elastic and brittle deformation.

Potential Implications for Seismic Hazard Evaluation

The amount of elastic energy that can be stored in the rocks adjacent to a locked fault, to be released when the locked fault ruptures in an earthquake, depends on the elastic strength of those rocks as well as on the frictional properties of the fault surface. An elastically weak rock can store more elastic energy for a given fault strength than a stiff rock (Cook, 1981). While Salleres and Ranero (2019) recently explored this phenomenon using Vp at the base of the upper plate as a proxy for the elastic modulus of the upper plate wedge, they discounted the impact of the lower plate by concluding that it does not affect potential down-dip variations in strain accumulation across the megathrust.

Our results indicate that, at least in Cascadia, the elastic properties of the lower plate should not be ignored when modeling subduction zone dynamics. The elastic strain implied by the imaged deflection is added to the strain resulting from plate bending due to subduction. We speculate that the transient deformation responsible for this added lower plate strain might play in role in triggering slip on the megathrust. As elastic strain energy accumulates over time in both the upper and lower plates due to locking on the plate boundary, the additional strain in the lower plate due to the ongoing deflection of the plate beneath the upper load may indirectly trigger slip on the plate boundary, and we pose this scenario as one worth exploring with more sophisticated strain accumulation and earthquake triggering models. As Arnulf et al. (2022) and Kimura et al. (2022) noted, the Nankai earthquakes of 1944 and 1946 initiated on either side of large upper plate pluton and propagated away from the pluton. This speculation is also consistent with the observation that great subduction zone earthquakes around the globe more often than not initiate at basement highs on the edge of forearc basins and propagate beneath the basins. (Song and Simons, 2003; Wells et al., 2003; Llenos and McGuire, 2007). The Wang et al. (2013) model for patchy slip during the 1700 Cascadia earthquake may also be consistent with this scenario because two patches of high slip flank the zone of lower plate deformation that we image. By analogy with Nankai, we speculate that rupture may have initiated here and propagated bilaterally in 1700.

We also speculate that our local result can provide an insight into the mechanical properties of the rest of the Cascadia subduction zone. The plate age in our study area is similar to the age of the subducted crust beneath the continental shelf along much of the Cascadia margin (Wilson, 2002). The lack of similar examples of focused lower plate flexure may be due in part to limited data coverage as well as to the fact that other sources of heterogeneous forearc loading have a more complex temporal history or are located farther from the deformation front, further complicating modeling to derive an estimate of effective elastic strength. Ultimately the 2021 onshore nodal (Cascadia 2021 Field Team, 2022), OBS data (Canales et al., 2023), integrated with the CASIE21 seismic reflection images (Carbotte et al., in review), will lead to extension of our 3D model up-dip (Nolan, 2023) and to the north (Ashraf et al., 2022) and south (work in progress).

Our results may also contribute insights into the interpretation of a 3D model for the electrical conductivity structure of Cascadia (Egbert et al., 2022). In Figure 6A, we show the outline of the low conductivity “g” anomaly, which represents a discontinuity in a high conductivity region immediately above the subducting plate that is interpreted to represent sediment subducted beneath the edge of the Siletz terrane (Egbert et al., 2022). We speculate that the load represented by the dike complex has pinched off the subducted sediment channel beneath the margin here. Future analyses can test this hypothesis by increasing the volume imaged seismically to increase overlap with the conductivity model and, potentially, through a joint inversion of seismic and electrical conductivity that uses the different sensitivity of these two imaging approaches to better constrain the spatial resolution of both.

Geodynamic modeling of plate interactions that incorporate the geometric relationship between the upper and lower plates imaged by this study as well as its rheological implications may also provide new insights into the driving forces that have caused brittle deformation of the upper and lower plates farther up-dip in the region characterized by the Alvin Canyon, Daisy Bank and Wecoma faults (Figure 1A; Goldfinger et al., 1997) and on possible connections between these faults, the deeply sourced seafloor seep known as Pythia’s Oasis (Philip et al., 2023), Holocene upper plate deformation that formed Stonewall Bank (Yeats et al., 1998), and clusters of recent seismicity (Morton et al., 2023).

Conclusion

A 3D Vp model of the central Cascadia margin between 44 and 45°N that includes a model of the Moho surface based on wide-angle reflections reveals that the forearc wedge contains a northwest-trending, near-vertical high-velocity anomaly that extends through the upper plate wedge. The anomaly is interpreted to be the dike complex that fed the mid-to-late Eocene Yachats Basalt, which is exposed in the Coast Range. This dike complex acts as a local, static load that deflections the plate boundary. The imaged scenario provides the opportunity to estimate the effective elastic thickness of the subducted Juan de Fuca plate beneath the margin. We conclude that:

• The experiment geometry, in which six parallel lines of offshore sources spaced 6–12 km apart and recorded by a grid of stations onshore spaced ∼8 km apart, complemented by ocean bottom seismometers, provides good resolution of the crustal velocity structure and the geometry of the Moho surface beneath the continental shelf and coastline.

• The velocity anomaly is ∼10 km wide, much narrower than the region of extrusive basalt that was mapped using field samples and aeromagnetic data.

• The Moho is depressed beneath the dike complex. The amplitude and width of the deflection, measured perpendicular to the curvature of the subduction zone, are 2–4 km and 70–120 km, respectively. A down-dip profile suggests at least partial elastic rebound.

• Multichannel seismic reflection data from the 2021 CASIE21 experiment as well as older data from the oceanic plate indicate that the plate boundary is located ∼6 km above Moho in this region.

• Assuming the Nafe-Drake relationship between velocity and density, the contrast between the dike complex and the adjacent mafic rocks of the Siletz terrane is ∼160 kg/m³.

• An apparent elastic thickness of 2–6 km is obtained from simple analytic elastic plate models for three different assumed load configurations, density contrasts of 130–190 kg/m³ and a range of estimates for the amplitude and width of the deflection, implying that the subducted Juan de Fuca plate is very weak.

• Although we assume pure elastic behavior, we cannot resolve whether the observed deformation is purely elastic or elasto-plastic (i.e., only partially recoverable). Elastic behavior is compatible with thermal models and observed seismicity but small offset faults may be present but are not resolved by our analysis.

• Elastic energy stored in the lower plate should be considered in models of interplate dynamics.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://ds.iris.edu/ds/nodes/dmc/data/types/time-series-data/.

Author Contributions

AT: conceived of study, obtained funding, organized field work, supervised analysis, drafted manuscript. KD: modeled Moho surface using PmP arrivals, generated figures, discussed interpretation, edited manuscript. CK: modeled Pg arrivals, generated figures, discussed interpretation, edited manuscript. SC: provided source for 2012 study, discussed interpretation, edited manuscript. JN: assisted with field work, discussed interpretation, edited manuscript. DT: developed inversion software, discussed interpretation, edited manuscript. WW: developed analysis software, discussed interpretation, edited manuscript. All authors contributed to the article and approved the submitted version.

Funding

This work was supported by NSF grant EAR1147975 and USGS grant G17AP00046 to Oregon State University.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We thank the many individuals who helped to acquire the field data in 1989, 1996 and 2012 and the public land managers and private landowners who provided access to their land. We thank Kelin Wang for interesting discussions. Ray Wells and Nick Harmon provided thorough and constructive reviews.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.escubed.org/articles/10.3389/esss.2023.10085/full#supplementary-material

References