The recycling of water into the Earth’s mantle via hydrated oceanic lithosphere is believed to have an important role in subduction zone seismicity at intermediate depths. Hydration of oceanic lithosphere has been shown to drive double planes of intermediate-depth, Wadati-Benioff zone seismicity at subduction zones. However, observations from trenches show that pervasive normal faulting causes hydration ~25 km into the lithosphere and can explain neither locations where separations of 25–40 km between Wadati-Benioff zone planes are observed nor the spatial variability of the lower plane in these locations, which suggests that an additional mechanism of hydration exists. We suggest that intraplate deformation of >50-m.y.-old lithosphere, an uncommon and localized process, drives deeper hydration. To test this, we relocated the 25 November 2018 6.0 MW Providencia, Colombia, earthquake mainshock and 575 associated fore- and aftershocks within the interior of the Caribbean oceanic plate and compared these with receiver functions (RF) that sampled the fault at its intersection with the Mohorovičić discontinuity. We examined possible effects of velocity model, initial locations of the earthquakes, and seismicphase arrival uncertainty to identify robust features for comparison with the RF results. We found that the lithosphere ruptured from its surface to a depth of ~40 km along a vertical fault and an intersecting, reactivated normal fault. We also found RF evidence for hydration of the mantle affected by this fault. Deeply penetrating deformation of lithosphere like that we observe in the Providencia region provides fluid pathways necessary to hydrate oceanic lithosphere to depths consistent with the lower plane of Wadati-Benioff zones.

The long-term evolution of both Earth’s mantle and hydrosphere is strongly influenced by the incorporation of water into the lithosphere of oceanic plates and its eventual subduction (e.g., Rüpke et al., 2006; Peslier et al., 2017, and references therein; Ohtani, 2021, and references therein). Understanding the mechanisms that drive this hydration and its extent requires inferences about the relationship between seismogenic behavior at depth within subducted slabs and observations of plate deformation or the seismic velocity structure within an oceanic plate’s mantle at or near the surface (e.g., Ranero et al., 2003; Emry and Wiens, 2015; Cai et al., 2018; Boneh et al., 2019; Obana et al., 2019). Intermediate-depth (60–300 km) Wadati-Benioff zone seismicity within subducted oceanic plates is often taken as a proxy for the extent of hydration of oceanic lithosphere (e.g., Peacock, 2001; Hacker et al., 2003; Ferrand et al., 2017), though not all proposed mechanisms for this seismicity require hydration (Kelemen and Hirth, 2007; Reynard et al., 2010; Prieto et al., 2013). Recent observations from the Northwest Pacific (e.g., Hasegawa and Nakajima, 2017, and references therein; Kita and Ferrand, 2018) and Southeast Pacific (e.g., Kumar et al., 2016; Bloch et al., 2018; Linkimer et al., 2020; Wagner et al., 2020) strongly suggest that intermediate-depth seismicity is associated with the release of fluids from subducting, hydrated oceanic lithosphere.

A clear association between the amount of intermediate seismicity within a subduction zone and the extent of outer-rise normal faulting is well documented both for specific subduction zones (e.g., Ranero et al., 2005) and globally (e.g., Boneh et al., 2019). The amount of intermediate-depth seismicity appears to increase with greater fault throw at the outer rise (Boneh et al., 2019) and decrease with greater thickness of sediments deposited on the subducting plate (Grevemeyer et al., 2005). The former relationship suggests that fault damage zones act as pathways for fluids to enter the plate (Grevemeyer et al., 2005; Boneh et al., 2019), while the latter suggests that thick sedimentary layers inhibit fluid entry (Grevemeyer et al., 2005).

Extensive efforts have been made to document hydration of oceanic lithosphere related to outer-rise faulting near trenches; however, significant hydration has been identified only to ~10–25 km below the oceanic Mohorovičić discontinuity (Moho) (e.g., Ranero et al., 2005; Grevemeyer et al., 2005; Contreras-Reyes et al., 2008; Emry and Wiens, 2015; Cai et al., 2018; Obana et al., 2019; Zhu et al., 2021). This contrasts with numerical models, which predict that deeper hydration of the plate should occur at the outer rise in older oceanic lithosphere (e.g., Iyer et al., 2012). These observations have left the cause and site of deeper hydration of the oceanic lithosphere contested (see review in Hasegawa and Nakajima, 2017, and references therein). Recent studies also indicate that the lower-plane Wadati-Benioff zone seismicity may require a more limited degree of hydration than that of the upper plane (e.g., Ferrand et al., 2017; Florez and Prieto, 2019), which suggests that other mechanisms of hydration may be responsible. Heterogeneities within the overriding plate’s forearc and a resulting localized increase in deformation of the subducting plate have been used to explain localized hydration up to 30 km below the oceanic Moho (Arnulf et al., 2022). Models focusing on sites of deep, limited lithospheric hydration away from trenches include thermal cracking as an oceanic plate ages (Korenaga, 2017) or the influence of plume magmatism (Seno and Yamanaka, 1996). Consideration of global Wadati-Benioff zone seismicity, discussed in detail in the next paragraph, suggests that thermal cracking of old oceanic plates is an unlikely cause of hydration associated with the lower plane of seismicity. Deep hydration and a lower plane of Wadati-Benioff zone seismicity have not been observed in association with the plume-derived Nazca Ridge (Kumar et al., 2016; Wagner et al., 2020) and O’Higgins Seamount Group/Juan Fernández Ridge (Kopp et al., 2004; Linkimer et al., 2020), which indicates that plume magmatism is unlikely to be a driver of deep oceanic lithosphere hydration.

Patterns in global Wadati-Benioff zone seismicity indicate that the cause of deep hydration associated with the lower plane of seismicity differs from that responsible for upper-plane seismicity. While global compilations of Wadati-Benioff zone seismicity plot broad, trench-parallel swaths of events to show a first-order trend of widening separation between the upper and lower planes of Wadati-Benioff zones with increasing age (or the thermal parameter derived from age) in the ~10–160-m.y.-old subducting oceanic lithosphere of the Indo-Pacific basins (Brudzinski et al., 2007; Florez and Prieto, 2019), other observations suggest that additional factors strongly affect this relationship. Most significantly, the 220–340-m.y.-old oceanic lithosphere (Speranza et al., 2012; Granot, 2016; Dannowski et al., 2019; El-Sharkawy et al., 2021; Segev et al., 2018, and references therein) subducting beneath the Hellenic arc shows only a single Wadati-Benioff zone plane (Halpaap et al., 2019). Likewise, the 90–110-m.y.-old lithosphere (Müller et al., 2008) subducting beneath the Lesser Antilles arc exhibits a single Wadati-Benioff zone plane (Bie et al., 2019). In addition, Wadati-Benioff zone seismicity within the lithosphere subducted by the Japan-Kuril subduction zones, when examined in detail, shows differing spatial characteristics between the Wadati-Benioff zone’s upper and lower planes: the upper plane is present throughout the region, while the lower plane exhibits a variable, highly localized patchiness (Kita et al., 2010). Taken together, these observations suggest that there may be two differing pathways for hydrating oceanic plates: (1) a general pathway common to all subducting plates that is responsible for the Wadati-Benioff zone’s upper planes as well as lower planes separated from upper planes by ≤25 km; (2) a much more localized pathway that is responsible for the lower planes separated from the upper planes by >25 km.

Here, we focus on a possible mechanism for generating localized pathways capable of hydrating oceanic lithosphere at depths sufficient to explain the deep lower plane of Wadati-Benioff zones: intraplate faulting within comparatively old oceanic lithosphere. This mechanism has previously been proposed based on activesource seismic survey observations of sub-mantle reflectors in the Wharton Basin of the Indian Ocean (Qin and Singh, 2015). These reflectors extend down to 30 km below the oceanic crust (Qin and Singh, 2015). In the same region, teleseismic observations have indicated that large earthquakes rupture to >45 km depth (Wei et al., 2013; Hill et al., 2015) and aftershocks, lacking formal depth uncertainties, may occur to at least 40 km depth (Kwong et al., 2019). However, this region lacks local to regional observations capable of providing well-constrained location uncertainties for earthquake sequences that may define possible fluid pathways or that are capable of detecting evidence for hydration along the segment of the fault ruptured by these earthquakes.

In the present study, we use a unique data set to examine an analogous sequence of earthquakes occurring within oceanic lithosphere that was observed with a local seismic network. Observation of intraplate oceanic settings is typically extremely challenging. Intraplate earthquakes observed in this setting are rare and have long recurrence times; sparse station coverage limits observations to relatively high-magnitude events (e.g., the 64 events over an ~50 year period compiled by Okal [1984] for the entirety of the South Pacific with a minimum magnitude of ~4.2 Mb). The rareness of these intraplate events makes local observations with oceanbottom seismometers difficult, unlike at oceanic plate boundary settings such as transform faults, which may exhibit quasiperiodic earthquake recurrence (e.g., McGuire et al., 2012). While teleseismic recordings may be used to detect and locate events greater than ~3–4 MW occurring in oceanic intraplate settings (e.g., Kwong et al., 2019), detecting small aftershocks and accurately relocating both small and large events requires stations within ~1.4× the depth of the events (Gomberg et al., 1990; Waldhauser, 2001). This requirement is difficult to meet in an oceanic intraplate setting; however, this was achieved with data recorded by the permanent broadband seismic stations of the Servicio Geológico Colombiano’s Red Sismológica Nacional de Colombia (RSNC) for the 25 November 2018 6.0 MW Providencia, Colombia, earthquake and its aftershocks (Fig. 1; Servicio Geológico Colombiano, 2018a, 2018b; U.S. Geological Survey, 2017).

We use arrival times from the Providencia earthquake’s mainshock along with 575 fore- and aftershocks to this event reported by RSNC to provide high-precision relocations of these events. While most stations in the region are too far from the fault segment ruptured by the Providencia sequence to be used in receiver function (RF) analysis, the presence of the RSNC’s station PRV on Isla Providencia directly above the ruptured segment provides an unusual opportunity to test for the effects of hydration at depth with RFs. We utilized seven years of data from station PRV, archived with the Incorporated Research Institutions for Seismology and the associated Earthscope Automated Receiver Survey (EARS) tool (IRIS DMC, 2010; Crotwell and Owens, 2005), to calculate RFs and tie observations of the Providencia earthquake sequence to hydration at >20 km depth.

The relocated events show that the rupture from the earthquake sequence extends beyond the oceanic crust of the Caribbean plate’s Lower Nicaragua Rise to reach a depth of ~40 km, which indicates that even moderately sized oceanic intraplate events can produce a pathway for fluids from the surface to the depths required to hydrate the mantle lithosphere and produce a lower plane of Wadati-Benioff zones. Our RF results show that the contrast in seismic velocities across the Moho in the vicinity of the fault are reduced relative to the contrast away from the fault, which indicates the presence of fluids or resulting serpentinization of the mantle around the fault and confirms its role as a fluid pathway. We then compare our findings in the context of prior observations of oceanic lithosphere both in intraplate settings and near trenches to propose a possible explanation for the spatially limited extent of lower-plane Wadati-Benioff zones.

The 25 November 2018 Providencia earthquake and its associated fore- and aftershocks occurred within the Caribbean plate in Colombia’s San Andres and Providencia Archipelago (Fig. 1). These islands and associated banks lie >~500 km from the boundaries of the Caribbean plate and near three significant tectonic features: the Pedro Bank Fault Zone, the San Andres Rift, and the Hess Escarpment. These features bound a distinct region of the Caribbean plate known as the Lower Nicaragua Rise.

Plate Boundaries and Major Intraplate Faults of the Western Caribbean Plate

The northwest Caribbean plate is separated from the Cocos plate by the Middle America Trench (Mann et al., 2006; Linkimer et al., 2010; see Fig. 1A) and from the North American plate by a variably transpressional to transtensional strike-slip boundary that includes the Motagua Suture Zone and Cayman Trough (Mann et al., 2006; Ratschbacher et al., 2009). The southern boundary is somewhat more complicated by incipient subduction of the Caribbean plate beneath the North Panama Deformed Belt (Camacho et al., 2010) and the transition to the Caribbean–South American plate boundary along the South Caribbean Deformation Front (Mantilla-Pimiento et al., 2009). The plate boundaries all lie ~500 km from the location of the Providencia sequence of earthquakes, while active deformation associated with these boundaries occurs at least 300 km from the sequence. This separation highlights the intraplate nature of these events.

Subplate-scale tectonic features exert a greater influence on the San Andres and Providencia Archipelago than distant plate boundaries. The most significant of these features is the Pedro Bank Fault Zone and the associated San Andres Rift. The Pedro Bank Fault Zone is a northeast-striking, left-lateral strike-slip fault separating the Upper (or northern) and Lower (or southern) segments of the Nicaragua Rise (Case et al., 1990; Holcombe et al., 1990; Ott, 2015; see overview of the nature of the Lower Nicaragua Rise below). Directly north of Isla Providencia, the Pedro Bank Fault Zone intersects the San Andres Rift, a north–south-striking extensional feature that joins the Pedro Bank Fault Zone and the Hess Escarpment (Carvajal-Arenas and Mann, 2018). Both the Pedro Bank Fault Zone and the San Andres Rift are presently active (Ott, 2015; Carvajal-Arenas and Mann, 2018) and represent pre-existing tectonic boundaries that have reactivated since the Miocene (Ott, 2015; Carvajal-Arenas and Mann, 2018; Torrado et al., 2019). The Hess Escarpment is an enigmatic and largely inactive structure (Holcombe et al., 1990) that parallels the Pedro Bank Fault Zone to the north. Only the southwesternmost segment of this feature, west of the San Andres Rift, has remained active into the late Cenozoic to present as a left-lateral, strike-slip fault (Bowland, 1993; Carvajal-Arenas and Mann, 2018). Previously, the full length of the Hess Escarpment has been interpreted as either a left-lateral fault active from the Cretaceous to early Cenozoic (Burke et al., 1984; Mauffret and Leroy, 1997) or a right-lateral fault active from the early to mid-Cenozoic (Boschman et al., 2014).

Smaller-scale faulting in the region (see Fig. 1B), like that active during the Providencia sequence, is subparallel to either the San Andres Rift or to the Pedro Bank Fault Zone/Hess Escarpment (see regional faults in Carvajal-Arenas and Mann, 2018; Idárraga-García and León, 2019; Idárraga-García et al., 2021). The proximity of the Providencia mainshock to a mapped fault trace southeast of Isla Providencia, the mainshock’s fault-plane orientation, and the fault-plane orientation of fore- or aftershocks (Idárraga-García and León, 2019; Idárraga-García et al., 2021; see Fig. 1C) strongly suggest that the Providencia sequence occurred along the northeast-trending Southern Roncador fault (alternately named the Serrana Fault System; Idárraga-García et al., 2021) ~20 km south of the island.

Nature of the Lower Nicaragua Rise and the San Andres and Providencia Archipelago

The origin of the Lower Nicaragua Rise remains contested. There is general agreement that a combination of continental crust (e.g., Case et al., 1990; Mann et al., 2006; Sanchez et al., 2019;) and oceanic arc crust (e.g., Holcombe et al., 1990; Mann et al., 2006; Lewis et al., 2011; Sanchez et al., 2019) makes up the Upper Nicaragua Rise north of the Pedro Bank Fault Zone, and that this crust differs significantly from that of the Lower Nicaragua Rise to the south of the fault zone (Holcombe et al., 1990; Mauffret and Leroy, 1997; Sanchez et al., 2019). The Lower Nicaragua Rise has been traditionally interpreted as a part of the overthickened oceanic crust of the Caribbean Large Igneous Province on the basis of seismic (e.g., Case et al., 1990; Holcombe et al., 1990; Mauffret and Leroy, 1997) and petrological (e.g., Feigenson et al., 2004; Gazel et al., 2011; Dürkefälden et al., 2019; Smith et al., 2021) data. However, recent geophysical modeling has suggested the Lower Nicaragua Rise may represent thinned continental or transitional crust, partly on the basis of isostatic and Bouguer gravity estimates of crustal structure (Carvajal-Arenas and Mann, 2018). In addition, 1661–454 Ma zircon xenocrysts recently identified in rocks from Isla Providencia may have been inherited from an underlying continental crust or crustal fragment (Smith et al., 2021), supporting a continental origin for the rise. However, these zircons may instead be detrital material derived from older neighboring terranes or contamination from paleo-Pacific material subducted during the Permian to Triassic, as has been proposed for other sites within the Caribbean (Smith et al., 2021, and references therein), consistent with an oceanic origin for the rise. The presence of these zircons complicates interpretation of the region. However, Srisotope values of the lavas hosting these xenocrysts are inconsistent with a continental crustal origin (Smith et al., 2021), again suggesting an oceanic origin for the rise.

Here, we favor the interpretation of the Lower Nicaragua Rise as a part of the Caribbean Large Igneous Province. Rock samples recovered from the Lower Nicaragua Rise are inconsistent with continental affinities and are geochemically similar to rocks recovered from elsewhere within the Caribbean Large Igneous Province (Dürkefälden et al., 2019). In addition, Miocene to Pliocene leadisotope data from volcanic rocks recovered from Isla Providencia are most consistent with the source mantle material beneath the island having been influenced by the Galapagos hot spot, which is the Caribbean Large Igneous Province’s likely origin (Feigenson et al., 2004). These rocks predate any potential slab tears or slab windows that have been suggested to allow interaction of hot spot material with the presently active Nicaraguan–Costa Rican volcanic arc (Feigenson et al., 2004; Gazel et al., 2011). Extensive petrological and geochemical investigation of the volcanic units making up Isla Providencia have further demonstrated that the lithospheric mantle in the region has a close affinity to both the Caribbean Large Igneous Province and the Galapagos hot spot (Smith et al., 2021). The seismic velocity structure of the Lower Nicaragua Rise east of our study area shows a crustal structure with a thin (~2–4 km) oceanic layer 2 and thick (14 km) oceanic layer 3 (Mauffret and Leroy, 1997). Both of these values are typical for overthickened oceanic crust formed by hot spots (e.g., Mutter and Mutter, 1993; Hampel et al., 2004). S-wave receiver function results calculated for the seismic station located on Isla Providencia indicate that the crust in our study area is ~20 km thick (Blanco et al., 2017) and comparable to the crust east of our study area. The P-wave velocities reported by Mauffret and Leroy (1997) for oceanic layer 2 (4.5–6.0 km/s) and oceanic layer 3 (6.0–8.0 km/s) are likewise typical for overthickened oceanic crust (e.g., Mutter and Mutter, 1993; Contreras-Reyes et al., 2019). While these velocity values may overlap with those of continental crust (Carvajal-Arenas and Mann, 2018), the overall structure and reported geochemistry are more consistent with the Lower Nicaragua Rise being part of the Caribbean Large Igneous Province or a closely related oceanic crustal feature.

We also caution that efforts to use Bouguer gravity or isostasy as constraints to determine the Lower Nicaragua Rise’s origin (e.g., Carvajal-Arenas and Mann, 2018) must consider positive buoyancy related to the proportionately thicker harzburgitic mantle layer associated with thicker oceanic crust (e.g., van Hunen et al., 2002) and the long-lasting positive buoyancy effects associated with hot spot-affected oceanic plates (e.g., Contreras-Reyes et al., 2019). Both features allow for a dense, relatively thin oceanic crust to support shallower bathymetry than is assumed in Bouguer gravity or isostasy modeling.

Isla Providencia is composed of a core of Miocene–Plioceneage volcanic materials (Concha and Macia, 1995; Pacheco-Sintura et al., 2014) surrounded by carbonate deposits (Geister, 1992). Exposed rocks on San Andres and other shallow bathymetric features in the region, including the Roncador and Serrana Banks, are composed primarily of carbonates (Geister, 1992; Carvajal-Arenas and Mann, 2018). These carbonate features have been presumed to be deposited on volcanic materials similar to those on Providencia (Geister, 1992). However, analysis of seismic reflection data indicates that at least San Andres represents an uplifted rift shoulder and may not lie on top of significant Cenozoic volcanic material (Carvajal-Arenas and Mann, 2018).

We use the compositions, layer thicknesses, seismic velocities, and uncertainties from these previous studies to construct and test possible regional velocity models for this part of the Caribbean Large Igneous Province and examine their effects on our results (see Effects of Velocity Model on Event Relocation section and Effects of Moho Depth section).

While the hot spot origin and resulting thick oceanic crust of the Lower Nicaragua Rise produce a petrologic structure that differs from that of ordinary oceanic plate material, the rheologic structure and patterns in seismicity of the rise are comparable to that of an ordinary oceanic plate. Both diabase/gabbro (Ranalli, 1997; He et al., 2003) and peridotite (Ranalli, 1997; Ueda et al., 2020) show brittle-to-ductile or semibrittle behavior to temperatures of 700–800 °C, which indicates that the lower part of the rise’s thick oceanic crust should behave in a rheologically similar way to the uppermost mantle of ordinary oceanic plates at comparable temperatures. As further explored below, the maximum depth of brittle deformation associated with the seismicity we observe for the Providencia sequence appears to be similar to thermally comparable oceanic lithosphere elsewhere (e.g., Hill et al., 2015; Kwong et al., 2019; and compilation in McKenzie et al., 2005). This similarity in brittle behavior indicates that processes inferred from our observations may be generalized to other intraplate oceanic settings.

Local to Regional Earthquake Relocation Data Set

We use the reported P-and S-wave arrival times, along with the hypocenter locations, provided by the RSNC catalog (Servicio Geológico Colombiano, 1993; for the region around Isla Providencia (see Fig. 1B), as well as the magnitudes (local magnitude [Ml], moment magnitude [Mw], and body wave magnitude [Mb]) calculated by the RSNC to better understand the Providencia sequence. The RSNC utilizes and reports arrival times from both Colombian and neighboring national and international networks (see Figs. 1A, 1B, and Fig. S11) to calculate hypocenter locations, and here we use reported arrival times from four RSNC stations (Servicio Geológico Colombiano, 1993), three Red Sísmica Nacional de Nicaragua (RSNN) stations (Instituto Nicaragüense de Estudios Territoriales, 1975), one station from the Caribbean U.S. Geological Survey Seismic Network (C-USGS, Albuquerque Seismological Laboratory/USGS, 2006), and one station from the Global Seismograph Network (GSN, Scripps Institution of Oceanography, 1986). Details about these stations are provided in Table 1.

To calculate our relocations, we extracted arrival times and origins from the RSNC catalog for all events within the square region from 12°N to 14°N and from 80°W to 82°W. We extracted events with reported depths greater than the approximate maximum depth to the seabed (≥2 km) occurring between 25 February 2018 and 4 July 2019 (the last available date at the time of extraction). Events with fewer than four stations reporting arrivals were not included in our data sets. This restriction effectively excludes the initial foreshock events that occurred before the deployment of the Roncador Bank (RNCC), Serrana Bank (SERC), and Isla San Andres (SAIC) stations, leaving us with a total of 493 potentially relocatable events out of 576 total events reported in the region. The first occurred on 27 March 2018, and the last occurred on 30 June 2019. We created two data sets by extracting arrival times for stations within a 300 km radius of each event to test possible velocity models focused on the four RSNC stations closest to the events and for a less restrictive 450 km radius to compensate for a temporary loss of station RNCC.

The 300-km-radius data set contains 490 events with 1871 P-phase and 1862 S-phase arrivals at the four RSNC stations and RSNN stations BLUN and ESPN. We applied a distance-based weighting to each arrival of 0.9 for station-event distances of less than 100 km, 0.8 for distances between 100 km and 150 km, 0.7 for distances between 150 km and 200 km, and 0.5 for distances greater than 200 km to emphasize observations from stations closest to the events. Double-difference relocation algorithms like the one that we use for relocation utilize the differential times between pairs of P-phases or pairs of S-phases for two or more events recorded at a station rather than the absolute travel times of individual phases. These differential times are the actual inputs used to solve for earthquake locations. Further data selection occurs automatically during the calculation of these differential times, which reduces the size of the actual data set used. These automatic criteria reduced the 300-km-radius data set to 406 events (83% of initial events) with a total of 63,699 P-phase differential times and 63,702 S-phase differential times for input to the relocation algorithm.

The 450-km-radius data set contains 493 events with 1986 P-phase and 1950 S-phase arrivals at the four RSNC stations and all five of the other stations. After applying the same distance weighting scheme to these events, the automatic criteria of the relocation algorithm reduced the data set to 407 events (83% of initial events) with 65,279 P-phase differential times and 65,233 S-phase differential times for input to the relocation algorithm.

Teleseismic Receiver Function Data Set

The RF data set for station PRV obtained from the EARS tool contains 127 pairs of radial and tangential RFs calculated using a low-pass Gaussian filter with a filter-width of 2.5 through the iterative time-domain deconvolution method (Crotwell and Owens, 2005; Ligorría and Ammon, 1999). The RFs were calculated from direct P-wave arrivals for teleseismic earthquakes with magnitudes of any type >5.0 at distances of between 30° and 100° between 31 October 2013 and 15 August 2018 (covering the full range of data accessible through EARS) as reported by the USGS National Earthquake Information Center Preliminary Determination of Epicenters (PDE) catalog. Ninety-five pairs of RFs failed the EARS internal quality checks. We then visually inspected the remaining 32 pairs of RFs to ensure no systematic errors had occurred before stacking the pairs by back azimuth for analysis.

We defined a magnitude of completeness (Mc) for the Providencia sequence using all 576 events reported in the region by the RSNC for the time period between 27 March 2018 and 30 Jun 2019. The Mc allows us to better understand the effects of the loss of station RNCC for the first ~3 months of 2019 on the network’s ability to detect events at the start of 2019 and to better characterize the sequence. To determine Mc, we followed Zúñiga and Wyss (1995) in assuming that seismicity in the region follows the Gutenberg-Richter relationship between event magnitude and frequency. In this case, Mc may be determined by plotting the log of the number of events with magnitudes greater than or equal to any particular magnitude value; Mc will mark the location where the plot shifts from an approximately linear trend.

Our calculated Mc of 2.2 (Fig. 2) for the Providencia sequence compared with the Mc of ~3 for the foreshocks alone provides some constraints for interpreting patterns in the foreshock seismicity of the sequence (Fig. 3). Small-magnitude events below the calculated foreshock Mc appear only after the installation of stations RNCC, SAIC, and SERC between 17 October 2018 and 29 October 2018, which is consistent with Mc being higher prior to this time and indicates that small foreshocks were not detected. MW magnitudes from the RSNC show a difference in magnitude of 0.9 between the mainshock and largest foreshock and 0.3 between the mainshock and largest aftershock, which suggests that the foreshocks may have been generally smaller than the aftershocks. The relatively small magnitude of these foreshocks again suggests that a number of small foreshocks may have been unrecorded. Large foreshocks (magnitude >4.0) began on 31 August 2018 and increased in magnitude until 4 October 2018 before a period of relative quiescence in the region that lasted until the mainshock occurred on 25 November. This period of quiescence is clear for foreshocks larger than the foreshock Mc of ~3, but it may be an artifact for events below this value.

The relatively long-lasting large-magnitude foreshock sequence contrasts with the aftershock portion of the sequence. Relatively large aftershocks essentially ceased by 5 December 2018, and aftershock productivity noticeably decreased by 14 January 2019. Unrest continued through the end of our obser-vational period, with events of Ml ~3 (comparable to the magnitude of events consistently recorded prior to the installation of additional stations) occurring as late as 12 June 2019. Beyond the period examined here, the RSNC reported 54 events of ≤Ml 4 up to 7 April 2021 in the area of the Providencia sequence.

Prior to the onset of the Providencia sequence, the RSNC reported only six events with magnitudes of between 2.8 and 4.0 at the southern margin of our study region and spatially distinct from the sequence. The spatial separation indicates that the events we analyzed as part of the Providencia sequence are not overprinted by background seismic activity in the region.

Relocation Using the HypoDD Method

HypoDD (Waldhauser and Ellsworth, 2000) is an approach to earthquake relocation that significantly reduces the influence of seismic velocity heterogeneities on the determination of location by comparing the differential arrival times reported at individual stations for events initially located within an arbitrary, relatively small distance of one another. This algorithm allows for improved recovery of spatial patterns in clustered seismicity, and by inference, the configuration of faults on which these clusters occur.

This approach requires eight observations for each pair of events to constrain the location and origin time. Event location and uncertainty estimation can be further improved by requiring that each event be linked to more than one neighboring event, which results in a network or constellation of events in which each event acts as a node. Here, we required that each event have a minimum of eight links with neighbors. We required that the reported RSNC catalog locations of neighboring events had to be within 40 km of each other prior to running the HypoDD inversion, and we found that, on average, neigh-boring events in our data set are within 9.6 km of one another. We ran each HypoDD inversion for a total of eight successful (meaning no earthquakes with negative depth were present) iterations using the singular value decomposition approach to obtain meaningful uncertainty estimates on event location (Waldhauser and Ellsworth, 2000). The first four iterations assigned P-arrival data a weighting of 1 and S-arrival data a weighting of 0.5, while the final four iterations assigned both arrival types a weighting of 1.

HypoDD’s use of both P- and S-wave arrival times can significantly improve event location and uncertainty estimates (Waldhauser and Ellsworth, 2000). The inclusion of S-wave data in event location has been shown to provide information that can greatly improve the estimation of event depth when at least one station is within a distance of <1.4× hypocenter depth (Gomberg et al., 1990), and station PRV is within this distance range for a large portion of the Providencia sequence. The incorporation of both P- and S-wave data has also been shown to greatly enhance the ability of small networks to accurately determine the location and depth of events (Chiu et al., 1997), as long as the velocity model and P-wave velocity to S-wave velocity ratio (Vp/Vs) are well constrained.

The P-wave velocity and largescale structure of the Lower Nicaragua Rise have been reasonably well constrained by prior studies, and we used these results to define a range of possible velocity models (see Table 2). These models include: (1) the sedimentary section based on thicknesses reported for the Isla Providencia region (Carvajal-Arenas et al., 2015) and average sedimentary P-wave velocities for the Lower Nicaragua Rise (Mauffret and Leroy, 1997); (2) the oceanic layer 2 and layer 3 velocities and thickness reported for the Lower Nicaragua Rise (Mauffret and Leroy, 1997); and (3) the Moho depth reported for the Lower Nicaragua Rise by Mauffret and Leroy (1997) and the S-wave receiver function equivalent below the station PRV identifiable in Blanco et al. (2017).

We follow Chiu et al. (1997) in using Wadati plots to calculate the regional Vp/Vs but utilize a bootstrapping approach to constrain uncertainty. We find that either using only the four closest stations or all stations in our data set yields a Vp/Vs of 1.72 and a 2σ uncertainty of ± 0.01. This ratio is consistent with calculated Vp/Vs values for mantle materials at 10 km to 30 km depth with a temperature of 400 °C to 500 °C (Abers and Hacker, 2016), which suggests that the majority of ray paths from the Providencia sequence largely pass through the uppermost mantle. The predominance of uppermost mantle paths is consistent with both the regional crustal thickness and relatively deep hypocenter locations reported by the RSNC. This interpretation is reinforced by calculating a modified Wadati plot to estimate the Vp/Vs of the crust (e.g., Chatterjee et al., 1985; Poveda et al., 2015) beneath the station closest to the Providencia sequence (PRV) and using a bootstrapping approach to constrain uncertainty. We find that the Vp/Vs of the crust beneath station PRV is 1.87 with a 2σ uncertainty of ± 0.03. This value is comparable to the intrinsic Vp/Vs value of 1.84 for oceanic crustal gabbro (Saito et al., 2015), as is expected for crust in this region.

In addition to testing alternate velocity models, we added perturbations to the RSNC arrival pick times and hypocenter locations to better understand possible uncertainties in our results. We tested uniform offsets of 10 km in each cardinal direction, uniform offsets in depths of 5 km and 10 km, and three sets of semirandom offsets (preventing offsets that would bring an event above our 2-km-depth cut off) capped at 5 km, 10 km, and 20 km distance from the RSNC’s reported hypocenters. The uniform offsets can be interpreted as testing systematic errors in event location, and the semirandom offsets as testing increased uncertainty in events’ initial locations. We also tested random perturbations to each reported arrival time capped at 0.1 s, 0.5 s, 1.0 s, and 2.0 s, as well as random perturbations applied uniformly to arrival times of each event using the same increments. These perturbations can be interpreted as tests of uncertainty in the choice of arrival time and event origin time, respectively.

Possible Sources of Relocation Uncertainty

We focus here on results obtained using our preferred velocity model and alternate depths for the Moho, as the Moho’s position has the greatest influence on our interpretation. Results from alternate velocity models are only summarized here, and plots of these results are available in the Supplemental Material (see footnote 1). We also examine possible effects of uncertainty in input event starting location and potential errors in arrival pick time or event origin time on our results. We find that while individual event locations may be variable under alternate velocity models or assumed input uncertainties, the broad structural features highlighted by the Providencia sequence are consistent and robust: the sequence occurred on a nearly vertical fault to ~15 km depth before continuing onto a potentially connecting fault dipping at ~50° to a depth of ~40 km, well into the mantle lithosphere.

Effects of Velocity Model on Event Relocation

We used the velocity models reported in Table 2 to relocate events in the Providencia region contained in our 300-km-radius data set. All models over-lie a half-space with a P-wave velocity of 8.0 km/s. These models represent variations on our preferred model determined from prior studies of the Lower Nicaragua Rise as described above in our discussion of the tectonic background of the region and were tested using only the four stations closest to the Providencia sequence to avoid possible biases from more distant stations. When these more distant stations were included for the final inversion, we were able to test the effects of changing the configuration of stations on our relocations.

Variation in the thickness of the sediment layer must account for the ~2 km water depth in the region. HypoDD is unable to account for negative topography, and its limitation to 1-D velocity models prevented us from more accurately accounting for water depth. Models with a 2.5-km-thick top layer (~0.5-km-thick sediment) or a 3.0-km-thick top layer (~1-km-thick sediment) both relocate only two events to a depth of <2 km, and the 3.5-km-thick top layer (1.5-km-thick sediment) results in only one event at <2 km depth, which suggests that our inability to better account for water depth has little effect on our ability to relocate events. All three models with varying upper-layer thickness successfully relocate essentially the same number of events (see Table 2), while average uncertainty (see Table 2) for the preferred 3.5-km-thick top layer model is notably better than that of the 2.5 km model and slightly better than that of the 3.0 km model. As seen in Figure S2 (footnote 1), the configurations of events for the preferred model and the 3.0 km model have no notable differences that would affect our interpretation of a nearly vertical feature to ~15 km depth rapidly transitioning to an ~50° dipping feature extending to 35–40 km depth. The 2.5 km model shifts all events slightly to the southeast and steepens the dipping feature to ~55°, which again has little effect on our overall interpretation.

These results suggest that our relocations are largely insensitive to the thickness of sediments in the region. To check the sensitivity to sediment thickness, we relocated the Providencia sequence using two models that take alternate approaches to removing the sediment layer while maintaining a 21.5 km depth to the Moho (see Fig. S3, footnote 1). The first thickens oceanic layer 2 by 3.5 km, while the second thickens oceanic layer 3 by the same amount. The thick layer 2 model relocates as many events as our preferred model with slightly greater average uncertainty, while the thick layer 3 model is able to relocate ~8% fewer events with notably greater average uncertainty. These results indicate some sensitivity to a lowvelocity top layer for our data set. Little difference is apparent in the relocations made using our preferred model or the thick layer 2 model results (Figs. S3A and S3B). The thick layer 3 model shows a noticeable shallowing of all events relative to our preferred model, whereas the thick layer 2 model shifts the dip of the dipping feature to ~40°, which terminates at a depth similar to that of our preferred model. This behavior is consistent with a greater number of unrelocatable “air quake” events and events at <2 km depth using the thick layer 3 model, which suggests this model poorly matches the region’s actual velocity structure.

Our test of varying the thickness of oceanic layer 3 suggests that our results may be sensitive to variations in this layer’s thickness and velocity. We investigate the effect of varying thickness in the section below by varying the depth to the Moho and focus here on two models with varying layer 3 velocities. We assign layer 3 a velocity of 6.5 km/s in one model and 7.5 km/s in the second, with both falling within the range of oceanic layer 3 velocities reported by Mauffret and Leroy (1997) for oceanic crust in the Caribbean. The slow layer 3 model successfully relocates ~9% more events than our preferred model with somewhat lower average uncertainty, while the fast layer 3 model successfully relocates ~40% fewer events with a much greater average uncertainty. Comparison of the plotted relocations (Fig. S4, see footnote 1) shows that the slow layer 3 model results in a shallower pattern of events overall, a shift of the intersection point between the near vertical and dipping features to ~10 km depth, a decrease in the dip angle of the dipping feature to ~40°, and a termination of the dipping feature at ~35 km depth. In contrast, the fast layer 3 model results in a poorly resolved pair of features with a possible change in dip at ~20 km depth. The large uncertainties, poor ability to relocate events, and inability to resolve structure led us to consider the fast layer 3 model as inaccurate for the region; however, the slow layer 3 model may be a reasonable representation for the region. We retain our preferred 7.0 km/s model, as overthickened oceanic crust in the Caribbean tends to have a greater portion of higher velocity lower crust in its oceanic layer 3 than the overall regional average (Mauffret and Leroy, 1997), which is more consistent with the 7.0 km/s model.

The relative success of the slow layer 3 model suggests that a model based on Carvajal-Arenas and Mann’s (2018) interpretation of the Lower Nicaragua Rise as extended continental crust may be an appropriate model for this region. To test this, we constructed a 3-layer model with a 6.0 km/s layer and a 7.0 km/s layer, respectively, corresponding to the upper and lower continental crustal materials of Carvajal-Arenas and Mann (2018). This model successfully located two more events than the slow layer 3 model with slightly higher uncertainties. Fig. S5 (see footnote 1) compares the relocations obtained using the extended continental crust model and our preferred model. The configuration of events below ~17–18 km depth is largely similar in both models, with a significant number of events reaching ~40 km depth. Shallower events show greater differences, and the vertical feature observed with our preferred model is eliminated. Events above ~13 km depth relocated with the extended continental crust model are also displaced somewhat to the southeast. Those nearest the surface do not locate near the mapped fault trace. This geometry is difficult to reconcile with the strike-slip motion observed within the focal mechanisms of the mainshock, largest foreshocks, and largest aftershocks (see Fig. 1C), which suggests that this model may be inaccurate for the region. These results, along with those obtained from varying the velocity of layer 3, also suggest that the region is more complicated than our simple 3-layer model. However, without additional constraints, we cannot further divide the layers of our model to account for this complexity.

Finally, we test both the low 1.71 Vp/Vs and high 1.73 Vp/Vs values permitted by the 2σ uncertainty from our Vp/Vs calculation described above. Comparison of results using these models (Fig. S6, see footnote 1) shows few differences between our preferred model and the high Vp/Vs model, with many events plotting in essentially the same location. The high Vp/Vs model relocates ~1% more events than our preferred model and has uncertainties that are nearly identical to those of our preferred model. In contrast, the low Vp/Vs model results are noticeably different from those of our preferred model, with ~1% fewer events than our preferred model, and lateral uncertainty is lower and vertical uncertainty is higher. The differences between the low Vp/Vs model and our preferred model do not significantly change our interpretation of a near-vertical feature intersecting a dipping feature. These features are shifted slightly to the southeast, while the dipping feature dips slightly steeper at ~60° and extends slightly deeper to >40 km depth.

Effects of Moho Depth

One of the more striking features observed in our results regardless of the velocity model examined in the preceding section is the intersection of a near-vertical feature and a feature dipping at a moderate to high angle. This feature could be an artifact related to the position of our velocity model Moho relative to the true Moho in the region. In addition, we have shown in the preceding section that events consistently relocate to form a feature that cross-cuts the Moho. The layer boundary positions of velocity models are known to some-times produce artifacts in event locations using other methods (e.g., Lomax, 2020); however, exploration of this problem in the existing literature using the HypoDD method has been limited. To explore these two issues, we present relocations calculated using four different Moho depths (Fig. 4).

We restricted Moho depths to values between 17.5 km and 23.5 km, bracketing the 20–21 km Moho depths reported beneath Isla Providencia (Blanco et al., 2017) and more generally within the Lower Nicaragua Rise (Mauffret and Leroy, 1997). Within this range, the depth to the Moho has little effect on the position of the intersection of the vertical and dipping features (Fig. 4). The intersection consistently remains above the modeled Moho, shifting slightly deeper with deeper Moho values. This trend is slightly obscured in the 17.5 km Moho model (Fig. 4A) due to an apparent gap in events around the Moho discussed in the latter part of this section and an increase in the dip of the dipping feature. The persistence of this feature in all models indicates that the shift from a near-vertical to dipping feature is not an artifact and that it occurs within oceanic layer 3 of the Lower Nicaragua Rise.

The cross-sections in Fig. 4 also highlight a clear artifact related to the Moho depth in the velocity models. A nearly event-free patch is apparent around the Moho in the 17.5 km Moho model (Fig. 4A). This gap is less pronounced but still apparent in the 19.5 km Moho case (Fig. 4B) and not present in the 21.5 km or 23.5 km Moho cases (Figs. 4C and 4D). This behavior is clearly more pronounced and spatially limited than the approximately ± 0.5 km variations induced by mismatched true and assumed velocity models examined by Michelini and Lomax (2004), and we have included a first-order examination of this problem in the Supplemental Material. We find that double-difference approaches tend to affect events occurring at or near a velocity boundary within the Earth in opposite ways if the assumed velocity boundary is an underestimate or overestimate due to differences between true and calculated ray paths (see Fig. S7 and Tables S1 and S2, footnote 1). If the velocity model underestimates the depth of a velocity boundary, the differences in ray path will deflect the relocated position of events occurring near the true boundary away from the assumed boundary. In contrast, if the velocity model overestimates the depth of a velocity boundary, the differences in ray path will attract the relocated position of events occurring near the true boundary toward the assumed boundary.

The deflecting effect of an underestimated boundary depth is much more visually apparent than the attracting effect of an overestimated boundary. To avoid using an overestimated boundary, we selected the shallowest Moho model with no apparent deflection of events: the 21.5 km model. Further attempts to refine the model and estimate the true Moho depth would be excessive given the vertical uncertainties in the relocated events. We also note that similar behavior should be expected at other prominent velocity boundaries in the structure of the Lower Nicaragua Rise, and our model boundary between oceanic layer 2 and oceanic layer 3 may have deflected some events in the 21.5 km Moho case (Fig. 4C). However, a similar deflection is not observed in the other cases (Figs. 4A, 4B, and 4D), nor is it consistently observed in models that vary the depth of the upper layer contacts (Figs. S2–S5). With no stations closer to these shallow events, we cannot determine whether the decrease in the number of events located around shallow-layer boundaries in some of the relocation results is a real feature, an artifact of our velocity model, or a result of poor constraints on depth for these shallow events. Therefore, we do not interpret this gap.

Effects of Input Event Locations, Origin Times, and Arrival Times

Having refined our velocity model, we now investigate the effects of input data from our 450-km-radius data set, which contains all nine stations in our study area (Fig. 1). As a relocation technique, results from the HypoDD method are dependent on input locations for events as well as observations of P- and S-wave arrival times. As an extreme example, Figure 5 compares a map view of relocations run based on events in the RSNC catalog shifted 10 km to the north, east, south, and west. These changes move the centroid point for the Providencia sequence—the reference point that the HypoDD method uses for its relative relocation calculation. As seen here, events cannot be relocated back to the RSNC catalog locations and remain centered on points ~10 km from these catalog locations. This particular pathology is not likely relevant to our study as the RSNC catalog locations match the location and orientation of a known fault from bathymetric data (see Fig. 1); however, less extreme cases could be relevant to our results.

To further test the effect of varying initial event locations, we prepared three sets of perturbed event locations to use as starting locations. These events have been moved a random distance in three dimensions by up to 5 km, 10 km, and 20 km with the condition that the perturbed locations do not have a depth of less than 2 km. This condition introduces a small downward bias to the events’ centroids, but it keeps the input data sets the same size as our unperturbed input data set. Figure 6 shows a comparison of results from all three perturbed data sets (Figs. 6A6C) and the unperturbed data set (Fig. 6D). The data set with a 5 km perturbation relocates ~1% fewer events than the unperturbed data set (see Table 3) with essentially identical average uncertainties, while the 10 km perturbation yields ~3% fewer events than the unperturbed data set with comparable lateral uncertainties but notably higher depth uncertainty. The 20 km perturbation continues this trend, yielding ~6% fewer events than the unperturbed data set and higher uncertainties. The relocations are variable at the level of individual events, as was also shown by the 10 km perturbation. They are potentially as variable at the scale of the entire sequence as the results obtained by varying the velocity model discussed above. However, the overall configuration of the sequence retains a vertical feature that intersects a dipping feature that continues to ~40 km depth. As the maximum lateral location uncertainty from the RSNC initial locations is 10.9 km and the maximum vertical location uncertainty is 15.4 km, the results of our randomly perturbed initial location data sets indicate that a possible mislocation of the input events should not significantly affect our results.

We tested the effect of errors in catalog event origin time by adding random values of up to ± 0.1 s, ± 0.5 s, ± 1.0 s, and ± 2.0 s to all arrivals from each event to create four perturbed data sets. The two smaller perturbations may also be interpreted as simulating alternative criteria for selecting arrival picks (e.g., selection of a zero-crossing or of an initial up/down motion). Figure 7 presents the relocation results of these data sets. The ± 0.1 s, ± 0.5 s, and ± 1.0 s data sets (Figs. 7A7C) show no significant change in the overall form of relocations and increasing progression in average uncertainty (see Table 3). This progression continues with the ± 2.0 s perturbed data set, but the events shift slightly southeastward (Fig. 7D). These minor changes again do not significantly affect our interpretation.

Finally, we tested the effect of errors in individual phase arrival times, simulating possible errors in picking introduced by high noise levels, by adding random values of up to ± 0.1 s, ± 0.5 s, ± 1.0 s, and ± 2.0 s to each arrival, again creating four perturbed data sets. Figure 8 presents the relocated results from these data sets. There is a clear decline in quality for the ± 1.0 s (Fig. 8C) and ± 2.0 s (Fig. 8D) perturbed data sets; in the ± 2.0 s case, no coherent features can be identified. The ± 0.1 s (Fig. 8A) and ± 0.5 s (Fig. 8B) perturbed data sets relocate a similar number of events with comparable average uncertainties. The ± 1.0 s and ± 2.0 s perturbed data sets’ relocation quality declines rapidly with the former relocating only ~84% of events with much higher average uncertainties than the unperturbed data set and the latter relocating only 50% of events with extreme average uncertainties. Investigation of average arrivaltime pick errors for P-wave arrivals by Zeiler and Velasco (2009) has shown that experienced network analysts on average have picking errors of <0.5 s for events at local distances with low signal-to-noise (as well as at regional and telesesimic distances generally) and <0.1 s for events at local distances with high signal-to-noise. Combined with our above analysis, these error values indicate that picking errors do not significantly affect our results.

These tests illustrate that by far the largest control on relocation quality are errors in individual arrivaltime picks. While errors in initial location, event origin time, or differences in event arrival selection criteria can be quite large before significantly affecting the HypoDD method’s results, errors in individual arrival times in excess of ~0.5 s significantly degrade the method’s ability to relocate events. We note, however, that errors in excess of ~0.5 s are unlikely to affect a significant fraction of our data set. If this type of error were to affect a significant part of our data set, it would not generate alternate discrete structures in the relocated events. Instead, it would result both in an inability to clearly identify discrete structures and in extreme location uncertainty for each event. Both of these are easily identifiable pathologies that are not evident in the unperturbed data set.

Improvements Compared to Catalog Locations and the Providencia Sequence’s Progression in Time and Space

After establishing a velocity model for the Providencia region and assessing the potential effects of uncertainty in hypocenter locations and arrival times, we ran our 450-km-radius data set with our preferred velocity model (see Table S3 for relocated events’ locations and uncertainties, footnote 1). Internal quality checks within the HypoDD program eliminated all picks from stations JTS and BCIP located in mainland Costa Rica and Panama, respectively. Their removal may reflect these two stations’ locations within a different terrane and above an active subduction zone, which differs significantly from our preferred velocity model. The addition of the three RSNN stations west of our study area allowed for the relocation of 19 more events than was possible using only the four stations of the RSNC. Sixteen of these occur in 2019, when data gaps for the RSNC stations would otherwise prevent event relocation. The addition of the RSNN stations has little effect on the average uncertainty of our event locations (compare Table 2 “Preferred” and Table 3 “Unperturbed, All Stations”) or on the location of events (compare Fig. 6D and Fig. 9B).

The relocation we obtained using the 450-km-radius data set and our preferred velocity model shows a large improvement in location uncertainty compared to the RSNC catalog’s initial event locations (Figs. 9B, 9D, 9A, and 9C, respectively), which allows us to clearly resolve structures in the region. Our results show that the Providencia sequence is significantly less clumped than the RSNC catalog locations indicate. Events in the sequence are less concentrated near the RSNC-reported hypocenter depth for the mainshock (14.31 km ± 3.80 km) and extend more clearly from near the seabed to ~40 km depth. While the mainshock’s relocation is too shallow (3.21 ± 0.31 km depth) given its distance to the nearest station to be certain its depth is accurate, this overall pattern is likely robust. In addition, the combination of a vertical feature transitioning to a dipping feature hinted at in the reported initial event locations is much clearer in our results, and event location uncertainty has been decreased by >3×. These improvements allow for further interpretation of spatial and temporal patterns within the Providencia sequence.

Foreshock activity appears to be largely limited to the deeper crust and upper mantle (Fig. 10A) and may have been limited to the dipping feature. However, depths reported in the RSNC catalog for foreshocks (see Fig. 1) prior to the deployment of additional stations in the region suggest that the events unrelocatable with our method may have extended to much shallower depths. Aftershock activity (Figs. 10B10F) extends from near the surface to >30 km depth within the first few hours after the mainshock and continues across this depth range to at least seven months after the mainshock’s occurrence. Aftershock activity appears to be somewhat more concentrated around the intersection of the subvertical and dipping features near 15 km depth. Taken together, the foreshocks and aftershocks indicate that rock fracture occurs concurrently from the near surface to ~20 km below the Moho.

Viewed along strike (Fig. 11), more interesting temporal patterns are clear. Both foreshocks (Fig. 11A) and initial aftershocks (Fig. 11B) appear to be limited to the middle of the sequence’s along-strike extent. Over the day following the mainshock (Figs. 11C and 11D), the primary center of aftershock activity moves away from the mainshock toward the northeast before becoming more widely distributed across the sequence’s extent over the following weeks (Fig. 11E) to months (Fig. 11F). Further analysis of these patterns would require a slip model of the mainshock.

An additional feature of the Providencia sequence is clear in Figure 11: its large spatial extent. Aftershocks for most earthquakes show a correlation with the internal and external boundaries of the mainshocks’ slip distributions (Das and Henry, 2003; Woessner et al., 2006) and may be used as a first-order approximation of a mainshock’s rupture area (Neo et al., 2021).

Following this relationship, the Providencia mainshock slip area would be on the order of 2000 km2. This area is significantly larger than the 100 km2 predicted slip area for a 6.0 MW event and comparable to that of a 7 MW event (Wells and Coppersmith, 1994). One possible explanation for this discrepancy may lie in the scattered distribution of aftershocks. Aftershocks are often interpreted as occurring at the boundaries of barriers that do not rupture during the mainshock (e.g., Aki, 1979; Das and Henry, 2003), and the scattered distribution of aftershocks in the Providencia sequence may indicate a very patchy or rough fault surface that allows only a fraction of the total area of the sequence to rupture during the mainshock. Alternately, this distribution may potentially indicate a complex triggering mechanism for the aftershocks.

Our relocation of the Providencia sequence earthquakes presented above indicates that the fault they occurred on runs from the sea floor to at least ~10 km into the mantle lithosphere. However, these relocations on their own cannot determine whether the fault acts as a channel for fluids to enter the mantle. RFs’ variable response to seismic velocity discontinuity characteristics provides us with a way to test for the effects of fluids on the mantle immediately below the Moho in the vicinity of the fault.

RFs are highly sensitive to discontinuities in the seismic velocity structure, like the Moho, below a seismic station. RFs calculated using the radial component of a teleseismic P-wave and associated P-to-SV conversions may be used to detect these discontinuities (e.g., Langston, 1979; Cassidy, 1992), while those calculated using the tangential component of the P-wave and associated P-to-SH conversions can be used to detect anisotropy at these discontinuities (e.g., Porter et al., 2011; Eckhardt and Rabbel, 2011). RFs have a frequency-dependent sensitivity to these discontinuities in both vertical and horizontal directions.

We calculated RFs for station PRV on Isla Providencia (Fig. 12A) using a 2.5 Gaussian value. This Gaussian value has a corner frequency of 1.2 Hz, which permits discrete layers of >~1 km thickness to be resolved. Complex discontinuities also have a significant effect on the vertical sensitivity of RFs, as shown in Fig. 12B by a set of synthetic RF responses to varying velocity structures. An abrupt discontinuity yields a single, relatively narrow peak in the RF’s response. A more gradational discontinuity yields a broader, lower amplitude peak. Sets of equivalent shouldered peaks result from complex boundaries composed of either a pair of abrupt discontinuities or an abrupt discontinuity overlying a more gradational discontinuity. We also note that both shouldered peak cases require an increase in velocity with depth. These shouldered peak cases emphasize both the non-uniqueness of the RF approach and how the widths of RF peaks in time may interfere with one another while still providing some constraints on the structure of a velocity discontinuity at depth.

Horizontal sensitivity of an RF is controlled by the first Fresnel zone of the converted S-wave, which is dependent on the overlying layers’ average S-wave velocity, the RF’s frequency content, and the depth to the boundary (Ryberg and Weber, 2000). For our preferred velocity model, the first Fresnel zone corresponds to a circle with a diameter of ~10 km at the Moho depth. Fig. 12C shows the first Fresnel zone for each RF in Fig. 12A calculated at 21.5 km depth. Since Isla Providencia lies in the footwall of the westward-dipping normal faults to the island’s west, RFs coming from back azimuths to the northwest sample an area largely undisturbed by faulting. In contrast, RFs coming from back azimuths to the southeast directly sample the area affected by the fault ruptured during the Providencia sequence (Fig. 12D).

The back azimuthal plot in Fig. 12A shows a clear difference between RFs arriving from the northwest and those arriving from the southeast. The radial components for the RFs arriving from the northwest show that the Moho there is likely a sharp discontinuity, while the tangential components show little to no evidence of anisotropy across the Moho. This observation contrasts significantly with the RFs sampling the Moho near the fault ruptured by the Providencia sequence. These radial RFs show reduced, shouldered arrivals to nearly non-existent arrivals, which suggests a more complex Moho discontinuity, while the tangential RFs show significant evidence for anisotropy across the Moho. The three RFs arriving from west–southwest back azimuths are too few to reliably interpret and may represent interference between features to the northwest and southeast of the island given the overlap of their first Fresnel zones (Fig. 12C). The variability of these RFs makes further quantitative analysis of the velocity structure beneath station PRV difficult (see Fig. S8 and associated discussion, footnote 1).

Taken together, these RFs indicate two very distinct velocity structures at the Moho to the northwest and southeast of Isla Providencia. To the northwest of Isla Providencia, the Lower Nicaragua Rise’s Moho is relatively pristine, representing a rapid transition from crustal to mantle velocities. To the southeast, the rise’s Moho has been affected by the Southern Roncador fault and requires a stepped or gradational transition from crustal to mantle velocities. This, in turn, requires a decrease in the average velocity of mantle materials just below the Moho, which indicates either the presence of a fluid at these depths or hydration and serpentinization of the mantle. Both possibilities are consistent with the evidence of anisotropy on the tangential RFs sampling to the southeast as both fluid-filled fractures (Crampin, 1981) and serpentinized mantle (Horn et al., 2020) produce S-wave anisotropy.

Given the correspondence between the unusual Moho sampled by RFs arriving from the southeast and the location of the Southern Roncador fault at depth inferred from our earthquake relocations, it is evident that the fault has acted as a conduit for fluids to at least reach the mantle material at >20 km depth. Fluid penetration along the fault to greater depths is possible but cannot be demonstrated with the present data set.

Understanding the potential extent of faulting in the region around Isla Providencia in both time and space is necessary for comparing our results with those of other regions of oceanic intraplate deformation. These comparisons would allow us to develop a new model for lithospheric hydration to depths of 35–40 km, which is necessary to explain instances of widely separated planes of Wadati-Benioff zone seismicity. Taken together, our observations of the Isla Providencia region and prior observations from other sites of oceanic intraplate deformation emphasize that the full tectonic history of an oceanic plate plays a key role in its hydration and behavior during subduction.

Regional Significance of Faults Revealed by the Providencia Sequence

Focal mechanism solutions for the Providencia sequence mainshock, as well as those calculated for foreshocks and aftershocks, correspond to a left-lateral, strike-slip fault with a minor component of extension. However, our results clearly show a more complex structure at depth, with a near-vertical fault intersecting a deeper fault plane dipping at ~50° to the northwest. This steep dip is most consistent with a normal fault; however, active normal faulting in the San Andres Rift immediately adjacent to this fault system follows a north–south trend (Carvajal-Arenas and Mann, 2018). This discrepancy in fault orientation suggests that the deeper fault is an inherited extensional structure, either a deeply penetrating normal fault or an extensional mantle shear zone (e.g., Reston, 1990) that is presently reactivated to accommodate strike-slip displacement.

Constraining the origin of this inherited deep fault allows for initial integration of our earthquake relocation results into the broader tectonic context of the Lower Nicaragua Rise. As the fault cross-cuts the overthickened oceanic layer 3 portion of the Lower Nicaragua Rise’s crust, it is unlikely that this inherited structure predates the emplacement of the Caribbean Large Igneous Province in the region at ca. 80–90 Ma (Nerlich et al., 2014; Dürkefälden et al., 2019). Subsequent to Caribbean Large Igneous Province emplacement, three extensional episodes were observed in the vicinity of Isla Providencia (Carvajal-Arenas and Mann, 2018). The most recent of these is associated with an increase in activity of the San Andres Rift and the establishment of significant left-lateral motion along the Pedro Bank Fault Zone in the mid-Miocene (Carvajal-Arenas and Mann, 2018; Torrado et al., 2019). This Miocene–present episode represents the initiation of the present kinematic regime and as such, the inherited fault must predate it. The earliest episode is characterized by early Eocene rifting and was centered to the west of the present San Andres Rift along a north–south orientation (Torrado et al., 2019). Both the early Eocene rifting’s distance from the Providencia region and its different orientation suggest that the kinematic regime in the early Eocene is unlikely to be responsible for the deep, reactivated fault. With the first and third extensional episodes ruled out, the late Eocene to Oligocene episode (Carvajal-Arenas and Mann, 2018) is the most likely origin of the reactivated fault. Normal faults interpreted as active during this episode based on seismic reflection survey data (Carvajal-Arenas and Mann, 2018) have a more complex set of orientations than those active in the other two intervals: the majority show a north–south trend, but a large subset show a northeast–southwest trend similar to the strike of the deep, reactivated fault south of Isla Providencia.

A number of similar Eocene–Oligocene faults in the region show evidence of reactivation. Faults bounding the eastern and southern sides of the Roncador Bank, including the eastward continuation of the fault ruptured by the Providencia sequence, cut slope deposits formed along the edges of the bank and indicate ongoing activity (Idárraga-García and León, 2019). The fault bounding the southern edge of the Serrana Bank likewise shows evidence of reactivation following the establishment of the Miocene–present kinematic regime (Carvajal-Arenas and Mann, 2018; Idárraga-García et al., 2021). Taken together, these faults define an area at least ~130 km × ~200 km. If similar, though possibly no longer active, features extend farther eastward into the more poorly studied sections of the Lower Nicaragua Rise, the region of deeply penetrating deformation may cover much of the ~900 km × ~250 km rise.

Comparison to Other Oceanic Intraplate Settings with Significant Deformation

Before comparing the Providencia sequence to sequences in other oceanic intraplate settings, the relationship between the age and thermal structure of the lithosphere within the Caribbean Large Igneous Province must be considered. This relationship may depart significantly from a simple plate-cooling model. If the lithosphere here is ~110 m.y. old, as is called for in some models of Caribbean Large Igneous Province formation (e.g., Nerlich et al., 2014; Whattam and Stern, 2015), the extensive volcanism associated with the emplacement of the Caribbean Large Igneous Province would have largely reset the thermal structure of the lithosphere. After the end of Caribbean Large Igneous Province and Lower Nicaragua Rise emplacement at ca. 75–80 Ma (Nerlich et al., 2014; Dürkefälden et al., 2019), cooling would resume as in a plate-cooling model. If the Caribbean Large Igneous Province formed contemporaneously with the lithosphere as other models suggest (e.g., García-Reyes and Dyment, 2021), the components of the lithosphere of the Lower Nicaragua Rise would have an age of ca. 73 Ma to <92 Ma and would have then cooled as in a plate-cooling model. The present observed depth to the lithosphere–asthenosphere boundary of ~80 km (Blanco et al., 2017) is approximately within measurement uncertainty for the predicted plate thickness of a 70–80-m.y.-old cooling plate (McKenzie et al., 2005). Estimates for the Curie Point (575 °C to 585 °C) depth for the region lie between ~45 km and ~55 km depth (Salazar et al., 2017), which again are comparable to estimates for a 70–80-m.y.-old cooling plate (McKenzie et al., 2005).

This thermal structure should be approximately comparable to that of the Central Indian Basin, where evidence for oceanic intraplate compressional deformation has been documented (Chamot-Rooke et al., 1993) within ~60–90 m.y. old oceanic lithosphere (Müller et al., 2008). Within the basin, a combination of low- and high-angle thrusts has accommodated 2.5% to 4.3% shortening over the last 7–8 m.y. (Chamot-Rooke et al., 1993; Krishna et al., 2001). They fall into two populations: small (<100 m), vertical throw faults spaced, on average, 3 km apart, and large, vertical throw faults spaced ~20 km apart (Delescluse and Chamot-Rooke, 2008; Delescluse et al., 2008). These faults extend 8–15 km below the surface of the oceanic crust (Delescluse and Chamot-Rooke, 2008), ~2–9 km into the mantle, and are associated with pervasive serpentinization to ~20 km beneath the surface of the oceanic crust (Louden, 1995; Delescluse and Chamot-Rooke, 2008), or to ~14 km into the mantle.

In many ways, these observations are more similar to those of fault throw, spacing, and extent of serpentinization found in outer-rise normal faulting (see Boneh et al., 2019, and references therein) than to our observations from the Lower Nicaragua Rise. Fault throw, taken as a proxy for full fault displacement, can be used as an estimate for the width of a fault’s damage zone (Savage and Brodsky, 2011; Faulkner et al., 2011)—an area in which repeated rock fracture allows enhanced fluid flow (e.g., Wibberley and Shimamoto, 2003). Piercing points have not been determined for the strike-slip systems in the Isla Providencia area; however, an approximate estimate can be made based on the relationship between fault length and displacement of strike-slip faults (Kolyukhin and Torabi, 2012). This relationship indicates that strike-slip systems with lengths of ~10+ km correspond to >1 km of long-term displacements, which indicates that the ~40 km fault ruptured during the Providencia sequence is likely associated with >1 km of displacement. This feature is comparable to the population of large offset thrust faults observed in the Central Indian Basin and is likely to produce comparably wide damage zones. Faults in the Providencia region are likely to be equally effective fluid conduits into the oceanic lithosphere, but their comparatively larger spacing and deeper extent likely leads to more focused and deeper serpentinization than is observed in the Central Indian Basin.

Observations of intraplate deformation from the Wharton Basin may also help to contextualize our results. Six large strike-slip earthquakes ranging from 6.8 MW to 8.6 MW have occurred over the past century within the Wharton Basin, extending from the outer rise associated with the Sunda Trench to >1000 km from the trench (see compilation in Heidarzadeh et al., 2017). These events appear to reactivate fracture zone faults and, in some cases, rupture multiple, near orthogonal faults (e.g., Wei et al., 2013; Hill et al., 2015; Kwong et al., 2019) within lithosphere ranging in age from ~35 m.y. to >80 m.y. old (Jacob et al., 2014).

Modeled rupture for the most recent of these events may show an interesting trend related to lithospheric age. The 7.8 MW 3 March 2016 event occurring in ~35–40 m.y. old lithosphere (Jacob et al., 2014) and modeled by Heidarzadeh et al. (2017) shows rupture to only 25–30 km depth. In contrast, within ~50–65 m.y. old lithosphere, models of rupture for the pair of >8.0 MW April 2012 events show rupture to 50–60 km depth (Wei et al., 2013; Hill et al., 2015).

A possible link between these latter events and the deep hydration of the oceanic lithosphere has been suggested based on interpretations of an active-source seismic survey line obtained along the outer rise of the Sunda Trench, ~110 km trenchward of the larger of the two April 2012 events (Carton et al., 2014; Qin and Singh, 2015). In this location, dipping reflections arriving at times corresponding to 45 km depth, ~25–30 km below the oceanic plate’s Moho, were interpreted as partially serpentinized faults lying below a more extensively serpentinized layer (Carton et al., 2014; Qin and Singh, 2015). These faults may have been reactivated at the same time deformation began within the neighboring Central Indian Basin (Carton et al., 2014). Both studies examining the survey failed to find a clear relationship between faults imaged above the Moho and these dipping features, and subsequent reinterpretation (Sibuet et al., 2019) has argued that these features are in fact artifacts from out-of-plane reflections of the serpentinized, shallow (~25 km depth, ~10 km below the Moho) portion of a vertical fracture zone that may be unimageable with available active source data (Sibuet et al., 2019).

Our observations from the Providencia sequence show that events of much smaller magnitude than those observed in the Wharton Basin are capable of rupturing through the brittle lithosphere to depths well below the oceanic Moho, and that these smaller events may rupture across intersecting faults. This observation, along with the observations of the 3 March 2016 Wharton Basin event, suggests that the age and thermal structure of the plate is a key factor that influences the depth of rupture and, by extension, the depth that a fault damage zone may reach. Our observations of the Providencia sequence also provide clearer support for the interaction of near vertical faults cutting through the oceanic crust and dipping structures extending to >40 km depth within the mantle, which are directly comparable to proposed structures in the Wharton Basin. In addition, our RF observations indicate that these faults can act as fluid pathways to allow serpentinization.

Taken together, comparison of the Providencia sequence and our RF results with observations from the Central Indian and Wharton Basins suggests that intraplate deformation within oceanic plates may be an effective way to hydrate oceanic lithosphere to depths >30 km below the oceanic crust. This deeply penetrating deformation appears to be limited to lithosphere with a thermal structure equivalent to a plate at least 50 m.y. old at the time of deformation. This age corresponds to the difference between the emplacement of the Caribbean Large Igneous Province in the Providencia region and the creation of late Eocene–Oligocene faults in the region (ca. 70–80 Ma and ca. 37–24 Ma, respectively), the difference between the youngest lithosphere in the Central Indian Basin and the onset of deformation (ca. 60 Ma and 7–8 Ma, respectively), and the 20–30 km rupture depth observed for younger (ca. 35–40 Ma) versus 50–60 km rupture depth observed for older (ca. 50–65 Ma) lithosphere in the Wharton Basin.

Implications for the Origin of the Lower Plane of Wadati-Benioff Zones

The upper plane of Wadati-Benioff zone seismicity is nearly a universal feature of subducting oceanic lithosphere, and its association with hydration due to outer-rise normal faulting is well established (e.g., Boneh et al., 2019, and references therein). Essentially all subduction zones for which the outer rise has been imaged using seismic methods have an identified low-velocity layer that extends to at most ~10–25 km below the oceanic Moho; it does not extend seaward beyond the outer rise (e.g., Middle America Trench: Van Avendonk et al., 2011; Japan Trench: Obana et al., 2019; Japan-Kuril Trenches: Fujie et al., 2018; Aleutian Trench: Shillington et al., 2015; Cascadia Trench: Horning et al., 2016; Marianas Trench: Cai et al., 2018; Zhu et al., 2021; Chile Trench: Contreras-Reyes et al., 2008).

In contrast, the separation observed between the upper and lower planes of the Wadati-Benioff zone is highly variable. Subducted material younger than ~70–80 m.y. in age shows a separation between these planes of ~10–20 km (Brudzinski et al., 2007; Florez and Prieto, 2019), which is consistent with an origin related to the pervasive uppermost mantle hydration observed between the outer rise and trench. Older subducted material shows a separation of ~25–40 km (Brudzinski et al., 2007; Florez and Prieto, 2019), which cannot be explained by observations of hydration at trenches. In addition, the lower plane is much less laterally continuous than the upper plane beneath the Japan and Kuril trenches (see plots in Kita et al., 2010) and non-existent beneath the Lesser Antilles and Hellenic trenches (Bie et al., 2019, and Halpaap et al., 2019, respectively), where material much older than 80 m.y. subducts. This difference strongly suggests that some other process not necessarily linked to a near-trench setting is responsible for deeper hydration in locations where separation of the upper and lower Wadati-Benioff zone planes exceeds 25 km.

We suggest that this deep hydration occurs in regions with oceanic lith-osphere older than 50 m.y. undergoing extensive intraplate deformation. As discussed in the preceding section, our observations from the Lower Nicaragua Rise indicate that ongoing faulting in this intraplate region has resulted in the creation of fault damage zones that act as fluid conduits and reach >20 km depth and possibly >40 km depth. While this extent is only ~20 km deeper than the Caribbean Large Igneous Province’s Moho, similar processes in a location with a more typical oceanic crustal thickness would result in localized hydration >30 km below the oceanic Moho. As noted in the preceding section, similar processes are ongoing in the Central Indian Basin, where serpentinization extends over a wide area, and in the Wharton Basin, where individual faults may allow localized serpentinization 25–30 km below the Moho.

A key feature of this process may be the relatively prolonged and repeated periods of activity that faults in these regions undergo. Studies of serpentinization in olivinepyroxene and natural peridotite compositions at high pressures and temperatures have shown that rates of serpentinization are 1–2 orders of magnitude higher than was determined by earlier studies of serpentinization in pure olivine and may extend to temperatures of up to 500–580 °C (Schwartz et al., 2013; Nakatani and Nakamura, 2016; Huang et al., 2020). These laboratory results strongly suggest that serpentinization is controlled by the rate at which water is supplied to the mantle. Water supply may be largely controlled by the formation of serpentine minerals themselves, as the increase in volume associated with their formation may clog pathways for water to move through the mantle until tectonic processes refracture the rock (Farough et al., 2016). Alternately, overpressure driven by the formation of serpentine minerals within rock may form additional cracks that promote flow (Kelemen and Hirth, 2012). This latter, self-sustaining process requires an increasing amount of over-pressure with depth and likely occurs at depths <10 km at low temperatures (Kelemen and Hirth, 2012) and shallower at higher temperatures (Skarbek et al., 2018). As a result, deeper hydration is likely reliant on earthquake activity to reopen fractures and sustain fluid flow deeper into the oceanic lithosphere.

Faulting in the intraplate settings discussed in the preceding two sections has remained active on the order of 10 m.y. or more. In contrast, faulting initiated around an outer rise, located ~50 km to 100 km from the trench (Lay et al., 2009), occurs over shorter time periods. Given that the convergence rates (DeMets et al., 2010) at subduction zones with widely separated Wadati-Benioff zone planes range from ~9 cm/yr (Japan) to ~8.2 cm/yr (Kuril) to ~1.9 cm/yr (Marianas), these outerrise faults may remain active for at most 1.1–5.3 m.y. Slip rates may provide some sense of how much refracturing of clogged fluid pathways could occur at these settings. Given an average outer-rise fault throw between ~300 m (Marianas) and ~450 m (Japan) for these subduction zones (Boneh et al., 2019), average slip rates for these faults range from ~0.6 mm/yr to ~4.1 mm/yr. These average slip rates are comparable to those estimated for the Providencia area of 1.19 mm/yr (Carvajal-Arenas and Mann, 2018). Thus, intraplate faults may be as active as outer-rise faulting over much longer periods of time, which would allow for more extensive hydration of deeper parts of the intraplate faults.

With this in mind, an interesting possible explanation for the 30–40 km separation of the Wadati-Benioff zone planes in the Northwest Pacific emerges. A model investigating possible past intraplate deformation of the Pacific plate associated with the subduction of the Izanagi-Pacific spreading center between 50 Ma and 60 Ma predicts extensive deformation within the region recently subducted by the Kuril, Japan, and Marianas trenches (Butterworth et al., 2014). During this period of deformation, this now subducted region would have been >60 m.y. old (Seton et al., 2012), which suggests that this deformation may have formed structures comparable to those discussed in the preceding two sections in relation to the Lower Nicaragua Rise, Central Indian Basin, and Wharton Basin.

To summarize, we propose a model that combines both the thermal structure and deformation history of an oceanic plate (Fig. 13) to explain the possible origin of deep lithospheric hydration necessary for the wide separation observed between the upper and lower planes of some Wadati-Benioff zones. Outside of slow to ultraslow spreading ridges with a spreading rate of <5 cm/yr and the immediate vicinity of fracture zones (e.g., Mével, 2003; Boschi et al., 2013; Rüpke and Hasenclever, 2017), serpentinization of the uppermost mantle is possible away from a spreading ridge but is likely limited (Mével, 2003). In the context of our model (Fig. 13A), the depth extent of serpentinization can be viewed as a result of the crust and uppermost mantle reaching a temperature below 600 °C, which thus triggers a switch from ductile to brittle behavior (Kohlstedt et al., 1995) in this portion of the lithosphere and then allows fluid access to the uppermost mantle. At the same time, this cooling allows the uppermost mantle to reach temperatures below 500–580 °C, at which serpentine minerals are stable (Schwartz et al., 2013; Nakatani and Nakamura, 2016; Huang et al., 2020). The slow retreat of the 600 °C isotherm deeper into the lithosphere prevents deep fracture and serpentinization >30 km below the oceanic Moho until the lithosphere is ~50 m.y. old and >40 km until it is ~60 m.y. old. Deformation prior to this age would result in a hydration depth nearly indistinguishable from that observed at an outer rise. For older lithosphere (Fig. 13B), periods of intraplate deformation result in deeply penetrating faults, and hydration around the damage zone of these faults becomes possible. Hydration of the uppermost mantle may be more extensive than in the fault damage zone due to both the greater permeability of serpentinized peridotite at depths of <12 km (Hatakeyama et al., 2017) and the potential for serpentinization to induce rock fracture at similarly shallow depths (Kelemen and Hirth, 2012). The deeper portion of these faults may be difficult to image using passivesource seismic techniques without optimally placed seismic stations as the hydrated areas may be narrow and widely spaced. As in other models of outer-rise hydration (e.g., Ranero et al., 2005; Iyer et al., 2012), as cold, old lithosphere approaches a trench (Fig. 13C), tightly spaced outer-rise faults provide a more extensive network of fluid pathways to hydrate the upper ~20–25 km of the oceanic plate, which allows for imaging of a shallow low-velocity zone using passive-source seismic methods (e.g., Fujie et al., 2018; Cai et al., 2018; Obana et al., 2019; Zhu et al., 2021). Once subducted, our proposed model converges with existing models of dehydration-induced (e.g., Hacker et al., 2003) or dehydration-driven, stress transfer-induced (e.g., Ferrand et al., 2017), intermediate-depth seismicity. Dehydration of the oceanic crust and of mantle hydrated by relatively uniform outer-rise faulting would form the upper plane of Wadati-Benioff zone seismicity (e.g., Hacker et al., 2003; van Keken et al., 2012), while dehydration of the sporadic, deeper intraplate deformation-induced hydration would form the lower plane of Wadati-Benioff zone seismicity (Hacker et al., 2003; Kita and Ferrand, 2018).

The sequence of seismicity associated with the left-lateral, strike-slip 25 November 2018 6.0 MW Providencia mainshock ruptured to a depth >40 km, cutting through the overthickened oceanic crust of the Lower Nicaragua Rise and the upper ~20 km of the mantle lithosphere. This rupture occurred along a vertical fault plane from the surface to ~15 km depth and an intersecting, ~50° dipping fault plane from 15 km to 40 km depth. This dipping plane is likely a reactivated Eocene–Oligocene normal fault. Additional mapped faults with the same orientation in the Lower Nicaragua Rise may also have the potential to be reactivated in this manner. RF results sampling the Moho where it intersects the fault on which the Providencia sequence occurred indicate that the fault has allowed fluids to reach depths >20 km, which potentially has serpentinized the mantle lithosphere.

This sequence provides the best-constrained observations of surface to ~40 km depth rupture in an intraplate oceanic setting. Similar deeply penetrating intraplate deformation has been proposed for the Central Indian and Wharton Basins. In all three settings, this deformation appears to be limited to oceanic lithosphere with a thermal profile equivalent to a >50 m.y. old plate. This style of deformation is capable of providing the deepreaching fluid path-ways necessary to drive hydration and serpentinization at a depth that may be associated with the lower plane of Wadati-Benioff zone seismicity within a subducted oceanic plate. This deep lower plane cannot be explained solely by observations of outer-rise–related hydration at trenches.

The long-term deformation history of oceanic lithosphere far from subduction zones plays a key role in its hydration, and therefore in its seismogenic behavior during subduction.

1Supplemental Material. PDF contains a map with labeled seismic station locations, plots of results obtained using alternate velocity models, an illustration of ray path differences induced by misestimated Mohorovičić Discontinuity depths, and an Hκ-plot for seismic station PRV. Excel file (Table S3) contains a catalog of relocated hypocenters obtained using the preferred velocity model. Please visit to access the supplemental material, and contact with any questions.
Science Editor: Andrea Hampel
Associate Editor: Giampiero Iaffaldano

This work was supported by National Science Foundation grant EAR-1760802. We also thank the staff of the Servicio Geológico Colombiano’s Red Sismológica Nacional de Colombia (RSNC) for assistance using the RSNC earthquake catalog; the hypocenter and arrival time data referenced here are openly available through the RSN’s online seismic catalog at We thank Giampiero Iaffaldano, Josip Stipčević, and an anonymous reviewer for providing comments that greatly strengthened this paper.

Gold Open Access: This paper is published under the terms of the CC-BY-NC license.