# Seismological evidence for a multifault network at the subduction interface

Caroline Chalumeau, Hans Agurto-Detzel, Andreas Rietbrock, Michael Frietsch, Onno Oncken, Monica Segovia, Audrey Galve

► To cite this version:

Caroline Chalumeau, Hans Agurto-Detzel, Andreas Rietbrock, Michael Frietsch, Onno Oncken, et al. Seismological evidence for a multifault network at the subduction interface. Nature, 2024, 628 (8008), pp.558-562. 10.1038/s41586-024-07245-y. hal-04717945

HAL Id: hal-04717945

https://hal.science/hal-04717945v1

Submitted on 2 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Seismological evidence for a multifault network at the subduction interface

https://doi.org/10.1038/s41586-024-07245-y

Received: 7 December 2023

Accepted: 27 February 2024

Published online: 17 April 2024

Open access

<figure><img src="image_2.png" /></figure>

Check for updates

Caroline Chalumeau<fnref n="1" />, Hans Agurto-Detzel<fnref n="1" />, Andreas Rietbrock<fnref n="1" />,²,³, Michael Frietsch<fnref n="1" />, Onno Oncken², Monica Segovia³ & Audrey Galve⁴

Subduction zones generate the largest earthquakes on Earth, yet their detailed structure, and its influence on seismic and aseismic slip, remains poorly understood. Geological studies of fossil subduction zones characterize the seismogenic interface as a 100 m–1 km thick zone<fnref n="1" />⁻³ in which deformation occurs mostly on metres-thick faults<fnref n="1" />,³⁻⁶. Conversely, seismological studies, with their larger spatial coverage and temporal resolution but lower spatial resolution, often image the seismogenic interface as a kilometres-wide band of seismicity⁷. Thus, how and when these metre-scale structures are active at the seismic-cycle timescale, and what influence they have on deformation is not known. Here we detect these metres-thick faults with seismicity and show their influence on afterslip propagation. Using a local three-dimensional velocity model and dense observations of more than 1,500 double-difference relocated earthquakes in Ecuador, we obtain an exceptionally detailed image of seismicity, showing that earthquakes occur sometimes on a single plane and sometimes on several metres-thick simultaneously active subparallel planes within the plate interface zone. This geometrical complexity affects afterslip propagation, demonstrating the influence of fault continuity and structure on slip at the seismogenic interface. Our findings can therefore help to create more realistic models of earthquake rupture, aseismic slip and earthquake hazard in subduction zones.

On 27 March 2022, a 5.8 moment magnitude ($M_w$) earthquake occurred at a depth of 19 km in the erosive Ecuadorian subduction margin, near the coastal town of Esmeraldas. The earthquake triggered a short-lived but highly productive aftershock sequence. At the time, a dense array of 100 short-period seismometers (Fig. 1 and Extended Data Fig. 1) was in place for the continuing offshore seismic experiment 'high-resolution imaging of the subduction fault in the Pedernales Earthquake Rupture zone' (also known as HIPER).

We detect and locate more than 1,500 events in the epicentral area between 12 March and 26 May 2022 (Fig. 1) using a combination of machine learning techniques⁸,⁹. Earthquakes are subsequently relocated in a three-dimensional (3D) velocity model¹⁰ with a double-difference algorithm¹¹, providing an image of the subduction interface seismicity at an unprecedented resolution. The seismicity clearly fell primarily along the megathrust at depths between 16 and 23 km in a 20 km wide along-strike region. Focal mechanisms show mostly oblique thrust faulting with a small right-lateral component. The earthquake locations broadly match the plate interface modelled by Font et al.,¹², and are consistent with previous offshore seismic studies¹³. A secondary cluster of seismicity associated with shallow normal faulting occurred in the upper plate, to the south-west of the mainshock, and was active intermittently throughout the observation period but at different times from the main earthquake sequence at the plate interface.

Fitting a plane to earthquakes with a location uncertainty below 75 m (Methods), we find that 90% of all well-located seismicity is contained

within an 870 m thick band with a strike and dip of 18° and 23°, respectively. However, seismicity is not purely planar (Fig. 2 and Extended Data Fig. 2) as a geostatistical analysis shows that topographic features have a characteristic width of 1.8 km and a characteristic height of 210 m, whereas the thickness of the seismicity is about 240 m (Methods). This is smaller than what has been previously observed seismologically in subduction zones⁷, probably due to our better resolution. It is, however, comparable to geological estimates of the thickness of fossil plate interface zones at seismogenic depths<fnref n="1" />,³,⁴,¹⁴. It also matches geophysical estimates of active plate interface thickness¹⁵, although studies imaging the deeper part of the seismogenic zone at a high resolution are rare. We infer that the observed band of seismicity occurred primarily within the plate interface region, although some seismic activity within damaged portions of the upper and lower plate cannot be excluded.

## Detachment fault versus broad deformation

The subduction interface can be viewed as a region in which the distribution of strain is heterogeneous, containing both strongly deformed zones and domains experiencing little deformation⁴,¹⁶. Within this region, it remains debated whether seismic slip concentrates on a single detachment fault, or whether it is distributed within the plate interface region¹⁷. Some studies have found evidence of seismic deformation concentrated along the roof fault zone, whereas aseismic slip is distributed within the subduction channel mélange⁴,¹⁸. Others have found

<footnote n="1">¹Karlsruhe Institute of Technology, Karlsruhe, Germany. ²GeoForschungsZentrum (GFZ), Potsdam, Germany. ³Institute of Geophysics, Escuela Politecnica Nacional, Quito, Ecuador. ⁴Université Côte d'Azur, CNRS, Observatoire de la Côte d'Azur, IRD, Géoazur, Valbonne, France. ⁵e-mail: rietbrock@kit.edu</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 1** | Map and cross-sections of the Esmeraldas sequence. **a–e**, Cross-sections (a–d) and map (e) of earthquakes with a location error less than or equal to 200 m, between 12 March and 26 May 2022. Cross-sections shown in **a,b,c** and **d** are marked on the map in **e** as green lines A–A', B–B', C–C' and D–D'. The stars show the position of the mainshock (yellow) and the three largest aftershocks (green, pink and red). A subset of 11 well-defined focal mechanisms is shown, including the United States Geological Survey focal mechanism of the mainshock (in red) and the focal mechanism of the largest direct foreshock (in blue).</figcaption></figure>

The locations for these focal mechanisms are also shown as labels in the cross-section. The slab model from ref. 12 is plotted as a black line in each cross-section and as dashed contours in map view. Nearby seismic stations are shown as yellow squares. **f** The inset shows the convergence velocity between the Nazca plate and the North Andean Sliver. Cross-sections **a–d** each have a width of 5 km and a bearing of 108°, corresponding to the dip direction of the main plane of seismicity.

evidence for several anastomosing fault strands hosting seismic slip within a structural width of hundreds of metres¹,¹⁹. It is, however, difficult to determine whether the activity on these different fault strands is concurrent¹,¹⁴, or if the different fault strands represent the gradual stepping down of the plate interface over many seismic cycles⁵,⁶.

Our data show seismicity occurring on several distinct, often sub-parallel fault planes (Fig. 2). These planes are identified by using the three-point method on families of similar earthquakes (Methods). They are mostly between 0 and 40 m thick with a median of 21 m, which is close to our resolution limit but similar to the thickness of individual fault zones found within subduction shear zones¹,³,⁶,²⁰. We can confirm that in some places several parallel planes are active at once, whereas in others seismicity is concentrated on only one plane. The superposition occurs in two main areas: the northern area of seismicity (Fig. 2b), and downdip of the mainshock rupture (Fig. 2d). Meanwhile, the area south of the epicentre is characterized by a single, thin plane of seismicity (Fig. 2c). Variations in the thickness of the reflective band interpreted as the plate interface have been previously observed at depth by seismic imaging¹⁵,²¹. For example, in Alaska, seismic imaging shows a marked increase in thickness associated with the change from a singular localized shear zone in the seismogenic zone (100–250 m thick low-velocity zone) to a wide deformation zone downdip (roughly 4 km thick zone with several low-velocity zones)¹⁵. Here in Ecuador, we show that seismic slip in the seismogenic zone itself is not necessarily concentrated on a single main fault, and that the distinction between single-fault

and superposed-faults regions is not linked to depth at this scale. This complexity of the plate interface geometry probably has repercussions on the propagation of seismic and aseismic slip in the region and on the earthquake rupture process.

## Spatio-temporal evolution of seismicity

Aftershocks at the plate interface occurred around the epicentre, propagating primarily in the north-north-west direction that is neither aligned with the strike nor the dip of the megathrust (Fig. 3a). The timing of these earthquakes follows the modified Omori law, with a *P* value of 0.99 (Extended Data Fig. 4). However, the area they covered increases with log of time (Fig. 3b), as does their cumulative number (Extended Data Fig. 5), indicating afterslip as the likely mechanism behind their occurrence²²,²³. Fluid diffusion cannot be the driver for the aftershock expansion, as the area would have to increase linearly with time²⁴, which we did not observe. Figure 3b also shows that the aftershock area expansion with time does not fit the model by Perfettini et al.²⁵ well. Instead, it shows far more complexity, slowing and speeding up and in one case perhaps overtaking the afterslip front²⁶. The aftershocks first occurred near the mainshock epicentre, initially spreading mostly south but gradually expanding in other directions. After about 12–20 min, the area covered by the aftershocks reaches the expected rupture size for a *M*<sub>w</sub> = 5.8 earthquake²⁷,²⁸. Downdip of the epicentre, a *M*<sub>w</sub> = 5.2 aftershock occurred 85 min after the mainshock and lead to
<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

**Fig. 2 | Map and cross-sections of earthquakes with location errors less than 75 m, showing the detailed structure of the seismicity.** a, Map of all cross-sections. Cross-sections B–D all have a bearing of 108°, corresponding to the dip direction of the main plane of seismicity. Cross-section E has a bearing of 18°, corresponding to the strike direction of the main plane of seismicity. Blue and red circles show families of similar earthquakes (CC > 0.75), aligned with the main plane of seismicity with an angle lower or higher than 30°, respectively. Each family shown here is assigned a lowercase letter for identification. b–d, Small cross-sections of regions B–D marked in yellow in the map in a, all at the same scale. The size of the crosses represents the earthquakes’ relative relocation error. Blue and red crosses show families aligned with the main plane of seismicity with an angle lower or higher than 30°, respectively. More cross-sections are shown in Extended Data Fig. 3. e, Cross-section of earthquakes along the strike of the seismicity, for the whole region. 95% of earthquakes are found between the yellow lines. Blue and red crosses show families aligned with the main plane of seismicity with an angle lower or higher than 30°, respectively. A histogram of earthquake depths relative to the plane is shown on the right, with the same vertical scale. f, Plots showing the thickness and orientation of the fault planes defined by families of similar earthquakes. The thickness plot only includes families with at least five events with relative location errors less than 20 m among them.

a small acceleration of the aftershock front, although seismicity did not expand far in that direction. North of the mainshock, a $M_w$ 5.4 and a $M_w$ 4.8 aftershock occurred about 16 and 16.5 h after the mainshock. Before these two earthquakes, the aftershock propagation slowed down near the location of the $M_w$ 4.8 aftershock (Fig. 3a,c). In fact, a

<figure><img src="image_7.png" /></figure>

**Fig. 3 | Evolution of aftershock area at the interface (for events with location errors less than or equal to 200 m).** a, Map of the aftershocks and aftershock front coloured by time. Each contour is the convex hull containing 95% of events at a given time. Stars represent the mainshock (white) and the largest aftershocks (coloured by time). b, Aftershock area as a function of log(time) (left-hand axis). Red stars show the times and magnitudes of the largest aftershocks (right-hand axis). Horizontal lines are estimations of the mainshock rupture area from scaling relations (dashed lines) or from the estimated size of the seismic gap (full black line). The full blue line shows the best-fit model assuming afterslip drives aftershock expansion (propagation velocity of 1.1 km per decade). The full grey line shows the best-fit model assuming fluid pressure diffusion drives aftershock expansion (diffusion coefficient of $1.42 \text{ m}^2 \text{ s}^{-1}$). Scaling relations: B2019, ref. 28; W1994, ref. 49; A2017, ref. 27. c, Aftershocks coloured by their time relative to large aftershocks. The stars show the mainshock (dark red) and the $M_w$ 5.2, $M_w$ 5.4 and $M_w$ 4.8 (light red, orange and white) earthquakes. Circles show the rupture areas of aftershocks assuming a stress drop of 3 MPa (ref. 41), occurring after the mainshock (dark red), $M_w$ 5.2 aftershock (light red), $M_w$ 5.4 aftershock (orange) and $M_w$ 4.8 aftershock (white). The light blue box represents a gap in seismicity, which serves as a first-order approximation of the mainshock rupture. d, Earthquakes coloured by time, projected onto a plane orthogonal to the streak indicated in c. The projected location of the seismicity streak is indicated as a dashed ellipse, and the projected main plane of seismicity is shown as a dark grey line. Timelapses are shown in Supplementary Videos 1 and 2.

north-east–south-west alignment of seismicity appeared during this time in this region (Fig. 3c,d). After these two aftershocks, there was a large increase in aftershock production and a notable acceleration of the aftershock expansion in the northern direction. The aftershock area eventually stopped expanding about 2.5 days after the mainshock, when it reached around $120 \text{ km}^2$. Subsequently, small ($M < 4$) earthquakes started to propagate towards the south-west of the mainshock after about 20 days, outside the main aftershock area.

The distribution of early aftershocks has often been used as a proxy to infer the extent and geometry of major earthquake ruptures in a mainshock-aftershock sequence²⁹. Here we show that the early increase in the area of the aftershock sequence is continuous, probably due to
the prevalence of early afterslip³⁰. Using early aftershocks to define the main rupture area seems therefore inherently arbitrary, as it is completely dependent on the cutoff time. However, there is a clear, 3.5 × 6 km gap in aftershock seismicity downdip of the mainshock epicentre, visible for the entire postseismic period (Fig. 3c). This gap probably represents the area where most of the accumulated strain was released by the mainshock³¹,³², and is therefore a conservative estimate of the rupture size, where slip most likely exceeds one-third of the maximum slip³³. North of the seismic gap, there is a topographic high in seismicity (Fig. 3 and Extended Data Fig. 2). This topographic high could have played a role in the rupture nucleation at its leading edge, whereas the stress shadow updip could have prevented the propagation of the rupture in that direction³⁴. Meanwhile, the seismicity is 'thick' at the downdip edge of the seismic gap, especially in the south where several fault planes are superposed over a thickness of more than 600 m (Fig. 2d and Extended Data Fig. 2). Here the geometrical discontinuity of the rupture fault plane could have contributed to stopping the earthquake rupture³⁵.

The spatial variability of the aftershock expansion is also worth exploring. The clear deceleration and acceleration that we see towards the north (Fig. 3a) around the time of the largest aftershock resembles other cases of step-like aftershock expansion thought to be related to fault segmentation³⁶, to pore fluid pulses³⁷ or to variations in frictional properties and/or the presence of subducted structure³⁸. In our case, the aftershock front was stalled at what is likely to be the edge of the M<sub>w</sub> 4.8 aftershock, where earthquakes form a straight lineament (Fig. 3c), hinting at the presence of a frictional or structural barrier. This boundary occurs in a region where the interface seismicity falls on a single thin fault, whereas both to the north and south of the boundary seismicity clearly occurs on several connected faults (Fig. 3d and Extended Data Fig. 2), suggesting that fault structure influenced the expansion of aftershocks. Mapping and studying such structures in detail could be crucial for our understanding of seismic versus aseismic slip in those regions.

## Afterslip deformation processes

Considering the subduction interface as a shear zone of finite thickness with potentially several fault planes has important implications for the way that we understand and model seismic and aseismic slip. With our thickness estimates we can place bounds on the shear strain rate associated with the afterslip, given as $γ = v_H/H$ in which v is the sliding velocity and H the thickness³. For this, we estimate the average afterslip displacement in two ways. We first assume that the afterslip moment is equal to 10–30% of the mainshock moment³⁹, and that it is reached by the end of our observation period. The total area of aftershocks after 1,227 h (51 days) is 170 ± 19 km², from which we subtract the estimated rupture area of 21–61 km² to obtain an afterslip area $A_{afterslip}$ of 109–149 km² (ref. 28) (Fig. 3b). The average displacement of the afterslip is given as $D = \frac{0.2 \times M_0^{\text{mainshock}}}{\mu A_{\text{afterslip}}}$, where the mainshock moment $M_0^{\text{mainshock}}$ is equal to $6.93 \times 10^{17}$ N m and the shear modulus μ is equal to 26 GPa (ref. 40). We obtain a displacement of 1.6–8.8 cm, whereas the average sliding velocity is equal to $3 \times 10^{-9}$ to $2 \times 10^{-8}$ m s⁻¹. Alternatively, we estimate the afterslip displacement by assuming that slip in the aseismic portions of the faults is equal to the average aftershock slip on the seismic portions of the fault. For this, we assume circular ruptures with an average stress drop Δσ of 2.5 MPa with a log error of 0.4 for all aftershocks⁴¹. Then, from ref. 42, the rupture area of a single aftershock is given as $A = \pi \times (\frac{7 M_0}{16 \Delta\sigma})^{\frac{2}{3}}$ whereas the displacement⁴³ is given as $D = \frac{M_0}{\mu A}$. We can therefore calculate the average seismic displacement on the megathrust as $D_{\text{average}} = \frac{\sum D_i A_i}{A_{\text{total}}}$ where D<sub>i</sub> and A<sub>i</sub> are the displacement and area of aftershocks within 800 m of the megathrust and A<sub>total</sub> is the total area where seismic ruptures occur, overlapping or not. Using this approach, the average displacement on the seismic portions of the fault is equal to $2.7 \pm 1.7$ cm, which overlaps with our first estimate

<figure><img src="image_8.png" /><figcaption>**Fig. 4 | Summary cartoon of the main features illuminated by the seismicity.** Aftershocks occur primarily around the mainshock rupture and expand along the plate interface as a function of the log of time, reflecting the likely influence of afterslip. We show that plate interface earthquakes occur on a network of metres-thick, subparallel, anastomosing faults at the plate interface, occurring over a thickness of a few hundred metres.</figcaption></figure>

and yields an average sliding velocity of $6 \pm 4 \times 10^{-9}$ m s⁻¹ with a maximum sliding velocity of $3 \times 10^{-5}$ m s⁻¹ in the first 5 min after the mainshock. If the aseismic deformation is distributed throughout the whole roughly 20 m width of the fault zones we observe with the seismicity, the strain rate is about $10^{-10}–10^{-9}$ s⁻¹ on average, reaching roughly $10^{-6}$ s⁻¹ immediately after the mainshock and providing a lower bound for the strain rate. An upper bound for the strain rate is provided by assuming a fault thickness-displacement scaling relationship as classically found for brittle faults⁴⁴. The observed afterslip displacements at $2.7 \pm 1.7$ cm would indicate a thickness of up to roughly 1 cm within the afterslip fault core¹,³,⁴⁴ and the corresponding strain rate would be around $10^{-7}–10^{-6}$ s⁻¹ on average and reach up to roughly $10^{-3}$ s⁻¹ after the mainshock.

The lower bound values of $10^{-10}$ to $10^{-6}$ s⁻¹ would be consistent with solution precipitation creep at low pore fluid pressures³, but tomographic data along the Ecuadorian subduction zone suggest that afterslip regions tend to correlate with anomalously high $v_p/v_s$ (P- and S-wave velocity, respectively) ratios, and hence probably high pore fluid pressures⁴⁵. Although foliated argillites may contribute to high $v_p/v_s$ (ref. 46), high pore pressures at convergent plate boundaries are also indicated by low megathrust earthquake stress drops⁴⁷, the force balance at convergent margins⁴⁸ and the presence of characteristic hydrofractures in exhumed plate interface zones³. In a region of high pore fluid pressure, the above range of strain rates cannot be accommodated by solution precipitation creep. Depending on the true thickness of the afterslip fault core—somewhere between the full fault width of roughly 20 m and that deduced from scaling relationships, roughly 1 cm—brittle creep appears the most likely mechanism for afterslip³.

## Megathrust structure and impact on slip

Our findings are summarized in Fig. 4. Thanks to our high precision microearthquake locations in the aftermath of a moderate-sized earthquake, we can illuminate the structure of the subduction interface at a depth of 15–20 km and examine in detail the spatio-temporal evolution of the seismicity. Looking at the timing of earthquakes, we find that fluid diffusion cannot explain the expansion of aftershocks in the region, meaning these are mostly controlled by afterslip. Meanwhile, the locations of earthquakes show that the seismic deformation is not constrained to a single continuous fault plane, but rather occurs on anastomosing, mostly subparallel metres-thick active faults that are sometimes superposed, in agreement with the geological record.
This fault zone complexity probably affects the propagation of after-
shocks and afterslip in the region, thus challenging our understand-
ing of the link between these two phenomena and the structures that
host them.

Our findings show the importance of conceptualizing the megathrust
as a fault network, particularly when discussing strain accumulation and
seismogenesis, rather than viewing it as a single plane. This is especially
relevant for creating more realistic models of earthquake rupture and
aseismic slip in subduction zones, and for a better assessment of the
related earthquake hazard.

Online content

Any methods, additional references, Nature Portfolio reporting summa-
ries, source data, extended data, supplementary information, acknowl-
edgements, peer review information; details of author contributions
and competing interests; and statements of data and code availability
are available at https://doi.org/10.1038/s41586-024-07245-y.

1. Rowe, C. D., Moore, J. C., Remitti, F. & the IODP Expedition 343/343T Scientists. The thickness of subduction plate boundary faults from the seafloor into the seismogenic zone. *Geology* **41**, 991–994 (2013).

2. Agard, P., Plunder, A., Angiboust, S., Bonnet, G. & Ruh, J. The subduction plate interface: rock record and mechanical coupling (from long to short timescales). *Lithos* **320–321**, 537–566 (2018).

3. Oncken, O., Angiboust, S. & Dresen, G. Slow slip in subduction zones: reconciling deformation fabrics with instrumental observations and laboratory results. *Geosphere* **18**, 104–129 (2021).

4. Bachmann, R. et al. Exposed plate interface in the European Alps reveals fabric styles and gradients related to an ancient seismogenic coupling zone. *J. Geophys. Res.* **114**, B05402 (2009).

5. Angiboust, S., Glodny, J., Oncken, O. & Chopin, C. In search of transient subduction interfaces in the Dent Blanche-Sesia Tectonic System (W. Alps). *Lithos* **205**, 298–321 (2014).

6. Wakabayashi, J. & Rowe, C. D. Whither the megathrust? Localization of large-scale subduction slip along the contact of a mélange. *Int. Geol. Rev.* **57**, 854–870 (2015).

7. Nippress, S. E. J. & Rietbrock, A. Seismogenic zone high permeability in the Central Andes inferred from relocations of micro-earthquakes. *Earth Planet. Sci. Lett.* **263**, 235–245 (2007).

8. Woollam, J., Rietbrock, A., Leitloff, J. & Hinz, S. HEX: hyperbolic event extractor, a seismic phase associator for highly active seismic regions. *Seismol. Res. Lett.* **91**, 2769–2778 (2020).

9. Woollam, J. et al. SeisBench—a toolbox for machine learning in seismology. *Seismol. Res. Lett.* **93**, 1695–1709 (2022).

10. León-Ríos, S. et al. 1D-velocity structure and seismotectonics of the Ecuadorian margin inferred from the 2016 M₉.78 Pedernales aftershock sequence. *Tectonophysics* **767**, 228165 (2019).

11. Zhang, H. & Thurber, C. H. Double-difference tomography: the method and its application to the Hayward Fault, California. *Bull. Seismol. Soc. Am.* **93**, 1875–1889 (2003).

12. Font, Y., Segovia, M., Vaca, S. & Theunissen, T. Seismic patterns along the Ecuadorian subduction zone: new constraints from earthquake location in a 3-D a priori velocity model. *Geophys. J. Int.* **193**, 263–286 (2013).

13. Garcia Cano, L. C., Galve, A., Charvis, P. & Marcaillou, B. Three-dimensional velocity structure of the outer fore arc of the Colombia-Ecuador subduction zone and implications for the 1958 megathrust earthquake rupture zone. *J. Geophys. Res. Solid Earth* **119**, 1041–1060 (2014).

14. Meneghini, F. & Moore, J. C. Deformation and hydrofracture in a subduction thrust at seismogenic depths: the Rodeo Cove thrust zone, Marin Headlands, California. *Geol. Soc. Am. Bull.* **119**, 174–183 (2007).

15. Li, J. et al. Downdip variations in seismic reflection character: Implications for fault structure and seismogenic behavior in the Alaska subduction zone. *J. Geophys. Res. Solid Earth* **120**, 7883–7904 (2015).

16. Bebout, G. E. & Penniston-Dorland, S. C. Fluid and mass transfer at subduction interfaces—the field metamorphic record. *Lithos* **240–243**, 228–258 (2016).

17. Fagereng, Å. Geology of the seismogenic subduction thrust interface. *Geol. Soc. Lond. Spec. Publ.* **359**, 55–76 (2011).

18. Kitamura, Y. et al. Mélange and its seismogenic roof décollement: a plate boundary fault rock in the subduction zone—an example from the Shimanto Belt, Japan. *Tectonics* **24**, TC5012 (2005).

19. Rowe, C. D. & Griffith, W. A. Do faults preserve a record of seismic slip: a second opinion. *J. Struct. Geol.* **78**, 1–26 (2015).

20. Meneghini, F. et al. Record of mega-earthquakes in subduction thrusts: the black fault rocks of Pasagshak Point (Kodiak Island, Alaska). *Geol. Soc. Am. Bull.* **122**, 1280–1297 (2010).

21. Nedimović, M. R., Hyndman, R. D., Ramachandran, K. & Spence, G. D. Reflection signature of seismic and aseismic slip on the northern Cascadia subduction interface. *Nature* **424**, 416–420 (2003).

22. Peng, Z. & Zhao, P. Migration of early aftershocks following the 2004 Parkfield earthquake. *Nat. Geosci.* **2**, 877–881 (2009).

23. Perfettini, H., Frank, W. B., Marsan, D. & Bouchon, M. Updip and along-strike aftershock migration model driven by afterslip: application to the 2011 Tohoku-Oki aftershock sequence. *J. Geophys. Res. Solid Earth* **124**, 2653–2669 (2019).

24. Shapiro, S. A., Patzig, R., Rothert, E. & Rindschwentner, J. Triggering of seismicity by pore-pressure perturbations: permeability-related signatures of the phenomenon. *Pure Appl. Geophys.* **160**, 1051–1066 (2003).

25. Perfettini, H., Frank, W. B., Marsan, D. & Bouchon, M. A model of aftershock migration driven by afterslip. *Geophys. Res. Lett.* **45**, 2283–2293 (2018).

26. Lange, D. et al. Comparison of postseismic afterslip models with aftershock seismicity for three subduction-zone earthquakes: Nias 2005, Maule 2010 and Tohoku 2011. *Geophys. J. Int.* **199**, 784–799 (2014).

27. Allen, T. I. & Hayes, G. P. Alternative rupture-scaling relationships for subduction interface and other offshore environments. *Bull. Seismol. Soc. Am.* **107**, 1240–1253 (2017).

28. Brengman, C. M. J., Barnhart, W. D., Mankin, E. H. & Miller, C. N. Earthquake-scaling relationships from geodetically derived slip distributions. *Bull. Seismol. Soc. Am.* **109**, 1701–1715 (2019).

29. Das, S. & Henry, C. Spatial relation between main earthquake slip and its aftershock distribution. *Rev. Geophys.* **41**, 1013 (2003).

30. Neo, J. C., Huang, Y., Yao, D. & Wei, S. Is the aftershock zone area a good proxy for the mainshock rupture area? *Bull. Seismol. Soc. Am.* **111**, 424–438 (2021).

31. Agurto, H., Rietbrock, A., Ryder, I. & Miller, M. Seismic-afterslip characterization of the 2010 Mw 8.8 Maule, Chile, earthquake based on moment tensor inversion. *Geophys. Res. Lett.* **39**, L20303 (2012).

32. Wetzler, N., Lay, T., Brodsky, E. E. & Kanamori, H. Systematic deficiency of aftershocks in areas of high coseismic slip for large subduction zone earthquakes. *Sci. Adv.* **4**, eaao3225 (2018).

33. Woessner, J., Schorlemmer, D., Wiemer, S. & Mai, P. M. Spatial correlation of aftershock locations and on-fault main shock properties. *J. Geophys. Res.* **111**, B08301 (2006).

34. Sun, T., Saffer, D. & Ellis, S. Mechanical and hydrological effects of seamount subduction on megathrust stress and slip. *Nat. Geosci.* **13**, 249–255 (2020).

35. Wnosky, S. G. Displacement and geometrical characteristics of earthquake surface ruptures: issues and implications for seismic-hazard analysis and the process of earthquake rupture. *Bull. Seismol. Soc. Am.* **98**, 1609–1632 (2008).

36. Kato, A. & Obara, K. Step-like migration of early aftershocks following the 2007 Mw 6.7 Noto-Hanto earthquake, Japan. *Geophys. Res. Lett.* **41**, 3864–3869 (2014).

37. Bedford, J. et al. A high-resolution, time-variable afterslip model for the 2010 Maule Mw=8.8, Chile megathrust earthquake. *Earth Planet. Sci. Lett.* **383**, 26–36 (2013).

38. Huang, H., Xu, W., Meng, L., Bürgmann, R. & Baez, J. C. Early aftershocks and afterslip surrounding the 2015 Mw 8.4 Illapel rupture. *Earth Planet. Sci. Lett.* **457**, 282–291 (2017).

39. Churchill, R. M., Werner, M. J., Biggs, J. & Fagereng, Å. Afterslip moment scaling and variability from a global compilation of estimates. *J. Geophys. Res. Solid Earth* **127**, e2021JB023897 (2022).

40. Geist, E. L. & Bilek, S. L. Effect of depth-dependent shear modulus on tsunami generation along subduction zones. *Geophys. Res. Lett.* **28**, 1315–1318 (2001).

41. Chalumeau, C., Agurto-Detzel, H., De Barros, L., Charvis, P. & the Rapid Response Team of the 2016 Pedernales Earthquake Spatio-temporal evolution of aftershock and repeater source properties after the 2016 Pedernales earthquake (Ecuador). *J. Geophys. Res. Solid Earth* **128**, e2022JB025353 (2023).

42. Eshelby, J. D. The determination of the elastic field of an ellipsoidal inclusion, and related problems. *Proc. R. Soc. Lond. Ser. Math. Phys. Sci.* **241**, 376–396 (1957).

43. Aki, K. Generation and propagation of G waves from the Niigata Earthquake of June 16, 1964; Part 2. Estimation of earthquake moment, released energy, and stress-strain drop from the G wave spectrum. *Bull. Earthq. Res. Inst. Univ. Tokyo* **44**, 73–88 (1966).

44. Torabi, A. & Berg, S. S. Scaling of fault attributes: a review. *Mar. Pet. Geol.* **28**, 1444–1460 (2011).

45. León-Ríos, S. et al. 3D local earthquake tomography of the Ecuadorian margin in the source area of the 2016 Mw 7.8 Pedernales earthquake. *J. Geophys. Res. Solid Earth* **126**, e2020JB020701 (2021).

46. Miller, P. K. et al. P- and S-wave velocities of exhumed metasediments from the Alaskan subduction zone: implications for the in situ conditions along the megathrust. *Geophys. Res. Lett.* **48**, e2021GL094511 (2021).

47. Wang, K. et al. Stable forearc stressed by a weak megathrust: mechanical and geodynamic implications of stress changes caused by the M=9 Tohoku-Oki earthquake. *J. Geophys. Res. Solid Earth* **124**, 6179–6194 (2019).

48. Dielforder, A., Hetzel, R. & Oncken, O. Megathrust shear force controls mountain height at convergent plate margins. *Nature* **582**, 225–229 (2020).

49. Wells, D. L. & Coppersmith, K. J. New empirical relationships among magnitude, rupture length, rupture width, rupture area, and surface displacement. *Bull. Seismol. Soc. Am.* **84**, 974–1002 (1994).

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_9.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2024
# Methods

## Event detection and location

Detection and association of phases was carried out using the machine learning tools available as part of the SeisBench package⁹. The continuous waveforms were scanned using PhaseNet⁵⁰, obtaining more than 6 million picks. Phase association was performed using the HEX method⁸, initially obtaining 112,336 events. This initial set of events was then relocated using NonLinLoc⁵¹,⁵² in a local one-dimensional velocity model¹⁰, and filtered in terms of number of phases and location accuracy, obtaining a high-quality subset of 1,775 events.

## Magnitude calculation

Local magnitudes were calculated using maximum peak-to-peak amplitudes and the original relation proposed by Bakun and Joyner⁵³ for California earthquakes (MCAL). These magnitudes were then benchmarked against the local magnitudes calculated by the Instituto Geofisico, Escuela Politecnica Nacional, Quito, Ecuador (MIG) for a subset of 148 common events. The obtained linear regression was MCAL = MIG + 0.49, which was then used to calculate the final local magnitude. We calculate the b value and magnitude of completeness in Extended Data Fig. 6.

## Focal mechanisms

Focal mechanisms for 24 selected events were computed using the grid-search algorithm HASH⁵⁴,⁵⁵, which takes the P-wave first motion polarities and S/P amplitude ratios as input data. For this, we used only stations within 50 km epicentral distance and clear P-wave onset.

## Relative relocation

We use the TomoDD software¹¹ with an existing 3D velocity model⁴⁵ to obtain accurate relative earthquake locations by using catalogue and correlation double-difference times. Cross-correlations are calculated on data from the vertical component filtered between 2 and 10 Hz, using a 2.1 and 3.3 s window for the P- and S-waves, respectively. If the correlation is above 0.75, the differential time is then used for relocation in TomoDD, with a weight equal to the square of the correlation coefficient. We require a minimum of eight phase links to define a pair and set the maximum separation between pairs to 6 km. Doing so still keeps all 1,571 successfully relocated events within a single cluster during the relocation. We also calculate errors through bootstrapping, because TomoDD does not give accurate error estimates. Two tests are used to determine errors. In the first jackknife test, we remove one station per event pair for 100 iterations to test the impact of the network on our results. In the second test, we add random noise to all our double-difference measurements for 100 iterations. The error is taken as the standard deviation of the locations obtained throughout the 100 iterations. The median location errors for the jackknife and bootstrap tests are 37 and 70 m, respectively, and the median errors in interevent distance are 22 and 33 m.

## Determining the main plane of seismicity

We define a plane representing the new local plate interface based on our data. Because most earthquakes between 15 and 25 km in depth with a relative location error below 75 m seem well-aligned on the plate interface, we simply fit a plane through them with a least-squares inversion. We then remove outliers, defined as earthquakes with depths within the lowest and highest 5% relative to the main plane. With this new set of data, we fit our final plane of seismicity, again with a least-squares inversion. The error is then calculated through a bootstrapping and jackknife test, in which we remove 10% of events and add a random error proportional to the event's location error to the locations for 500 iterations. Overall, our plane has a strike of 17.7 ± 0.6° and a dip of 23.3 ± 0.2°, with the seismicity forming an 872 ± 27 m thick band. However, these errors are probably underestimated due to the flexure

and topography of the plate interface, as well as to the possible inclusion of interplate seismicity in the calculation.

## Finding similar and repeating earthquakes

We perform cross-correlation on the vertical component, this time with a filter of 4–12 Hz and a window of twice the S and P time. Low correlations are discarded if the associated signal to noise ratio is below 15 to ensure that noise does not limit our classification.

Similar earthquakes are defined as earthquakes with correlations above 0.75 for at least 80% of usable stations (example in Extended Data Fig. 7). The 0.75 threshold was chosen on the basis of the histogram of cross-correlations found in Extended Data Fig. 8. Indeed, this value separates two main groups of events: non-similar events, whose correlation can be positive or negative but with a peak absolute value around 0.25, and similar events, which have high correlations with a peak around 0.95. This indicates that they probably are close in space and have similar focal mechanisms, and thus probably occur on the same or on close, parallel faults. We further define these individual fault planes with the three-point method⁵⁶, using earthquakes within the family that have location errors of less than 20 m relative to each other and less than one tenth of the interevent distance. Fisher statistics are used to find the best-fit plane and its associated uncertainties. We retain well-defined fault planes if the value of K is greater than five, the angle uncertainty is under 20° and the probability is greater than 80%.

Repeating earthquakes are defined as earthquakes with correlations above 0.92 for at least 80% of usable stations. In addition, their interevent distance must be smaller or equal, within errors, to the radius of the largest event's rupture. The latter is calculated assuming that moment magnitude $M_w = M_t - 0.5$ where $M_t$ is the local magnitude, with a stress drop of 3 Mpa (ref. 41). By using both criteria, we ensure that repeaters truly rupture the same area. We find ten pairs of repeating earthquakes, one of which occurs exclusively in the preseismic period and one of which has one preseismic and one postseismic event. The preseismic repeater occurs in the crustal cluster, whereas all other repeaters occur near the plate interface. However, it is likely that this repeating earthquake catalogue is incomplete due to the small time-window examined, as many asperities probably did not have time to reload.

## Geostatistical analysis and seismicity thickness calculations

The variations in elevation of earthquakes relative to the best-fit plane are analysed using an experimental half-variogram, defined as:

$$ \gamma(d) = \frac{1}{2N(d)} \sum_{i=1}^{N} [h(x_0 + d) - h(x_0)]^2 $$

where $h(x_0)$ is the elevation of a given earthquake relative to the best-fit plane, and $h(x_0 + d)$ is the elevation of another earthquake at a distance $d$. $N(d)$ is the number of pairs at a distance $d$ and $\gamma(d)$ is the semivariance of the elevation difference as a function of interevent distance. The half-variogram reaches a constant value, known as the sill, at a given distance, known as the correlation length or the effective range. At that distance, separated values of elevation are no longer correlated at all, meaning this corresponds to a characteristic width of topographic features or of patches of a given thickness. Meanwhile, the sill is related to the characteristic height of topographic features, as it indicates the overall variance of the elevation. The value of the half-variogram at zero is known as the nugget, and in our case is linked to the thickness of the seismicity at a given point. We define the thickness as the width containing 95% of the seismicity, approximated as four times the square root of the nugget.

We calculated the half-variogram every 50 m over 4 km, using more than 108,000 earthquake pairs occurring within 800 m of the main plane of seismicity with a relative elevation error of less than 20 m and a distance error of less than 25 m (Extended Data Fig. 9). All errors
cited here are calculated by removing 10% of the data and adding noise
before fitting either model for 100 iterations. Our data are best fit by the
exponential model, which yields a correlation length of 1.82 ± 0.03 km
with a nugget of 3,511 m² and a sill of 43,980 m², for which we estimate
the thickness of the seismicity at around 237 ± 18 m and the character-
istic height of topographic features at around 210 ± 1 m. This estimate
of thickness is almost equal to the 227 m calculated using only event
pairs less than 50 m apart. By comparison, with the less well-fitting
spherical model, we obtain a seismicity thickness of 377 ± 11 m and a
characteristic height and correlation length of topography of 194 ± 1 m
and 1.59 ± 0.04 km, respectively.

As the thickness of the seismicity is particularly model dependent, we calculate an extra upper bound. To do so, we divide the megathrust into smaller regions of 1.8 × 1.8 km, equal to the exponential model's correlation length, thereby removing the effect of plate flexure and wide topographic features. We remove squares containing fewer than ten events and only keep results if the standard deviation of the jack-knife test we perform is smaller than 75 m. Doing this, we obtain an average thickness of 455 m with 95% of values falling between 164 and 958 m. This is larger than our previous estimates, as it is further affected by narrow or steep topographic features. This approach also allows us to map thickness variation and large-scale topographic features (Extended Data Fig. 2), thus giving an overall view of the structure of the plate interface.

Data availability

Earthquake catalogues created during this study are available in a data
supplement to this paper at https://doi.org/10.35097/1921, along with
the absolute and differential traveltimes of all earthquakes and the
location errors calculated by bootstrapping.

50. Zhu, W. & Beroza, G. C. PhaseNet: a deep-neural-network-based seismic arrival-time picking method. Geophys. J. Int. 216, 261–273 (2019).

51. Lomax, A., Virieux, J., Volant, P. & Berge-Thierry, C. in *Advances in Seismic Event Location* (eds Thurber, C. H. & Rabinowitz, N.) 101–134 (Springer, 2000).

52. Lomax, A., Michelini, A. & Curtis, A. in *Encyclopedia of Complexity and Systems Science* (ed. Meyers, R. A.) 1–33 (Springer, 2014).

53. Bakun, W. H. & Joyner, W. B. The ML scale in central California. *Bull. Seismol. Soc. Am.* **74**, 1827–1843 (1984).

54. Hardebeck, J. L. & Shearer, P. M. A new method for determining first-motion focal mechanisms. *Bull. Seismol. Soc. Am.* **92**, 2264–2276 (2002).

55. Hardebeck, J. L. & Shearer, P. M. Using S/P amplitude ratios to constrain the focal mechanisms of small earthquakes. *Bull. Seismol. Soc. Am.* **93**, 2434–2444 (2003).

56. Fehler, M., House, L. & Kaieda, H. Determining planes along which earthquakes occur: method and application to earthquakes accompanying hydraulic fracturing. *J. Geophys. Res. Solid Earth* **92**, 9407–9414 (1987).

**Acknowledgements** This work was supported by the ANR FLUID2SLIP project, grant no. ANR-21-CE49-0012 of the French Agence Nationale de la Recherche. Funding from the Geophysical Institute at KIT for field work and C.C. is acknowledged. We thank everyone involved in the deployment of land stations during the HIPER experiment, including W. Acero, I. Tapa, D. García, C. Viracucha, Christian Espin, J. Santo, F. Mejía, V. Reyes-Wagner, S. León-Ríos, L. Mejía and B. Braszus.

**Author contributions** C.C. performed the relative relocation with TomoDD, the search for repeaters and most of the data analysis (assisted by the other authors of the manuscript). H.A.-D. computed initial absolute locations, magnitudes and focal mechanisms. A.R. developed the conceptual idea of this study and acquired funding. M.S., A.G. and A.R. coordinated data collection. All authors assisted in discussing the results and editing the manuscript.

**Funding** Open access funding provided by Karlsruher Institut für Technologie (KIT).

**Competing interests** The authors declare no competing interests.

**Additional information**

**Supplementary information** The online version contains supplementary material available at
https://doi.org/10.1038/s41586-024-07245-y.

**Correspondence and requests for materials** should be addressed to Andreas Rietbrock.

**Peer review information** *Nature* thanks Thorne Lay and Christie Rowe for their contribution to the peer review of this work. Peer reviewer reports are available.

**Reprints and permissions information** is available at http://www.nature.com/reprints.
<figure><img src="image_10.png" /></figure>

**Extended Data Fig. 1 | Map of seismic stations (blue squares).** The focal mechanism from the Global CMT catalogue is shown in light blue, at the relocated location. The light blue rectangle shows the region mapped in Fig. 1.

Grey shading is used to illuminate the topography and bathymetry from the Global Earth Relief Grids dataset, available from GMT (NASA, 2013).
<figure><img src="image_11.png" /></figure>

**Extended Data Fig. 2** | Maps of seismicity within 800 m of the main plane with a location error of 75 m or less. A: Thickness containing 95% of plate interface events within a 0.9 km radius. While apparent thickness can also be affected by steep topography, intraplate earthquakes and sparse data, it can serve as a first-order proxy for fault structure. In this map, areas with thicknesses above ~450 m have multiple parallel planes of seismicity. B: Median elevation of events relative to the main seismicity plane, calculated within a 1.6 km radius.
<figure><img src="image_12.png" /></figure>

**Extended Data Fig. 3** | Additional cross-sections of earthquakes with location errors <75 m. Blue and red circles show families of similar earthquakes (CC > 0.75), aligned with the main plane of seismicity with an angle lower or

<figure><img src="image_13.png" /></figure>

higher than 30 degrees respectively. In the cross-sections for the whole region, 95% of earthquakes are found between the yellow lines.
<figure><img src="image_14.png" /></figure>

**Extended Data Fig. 4 | Omori plot for earthquakes occurring at the plate interface.** Blue dots show the real aftershock rates, while the blue line is the best-fit line with a least-squares inversion.
<figure><img src="image_15.png" /><figcaption>**Extended Data Fig. 5** | Cumulative number of earthquakes at the plate interface. A: Cumulative number of earthquakes in linear time, with large magnitude aftershocks as red stars. B: Cumulative number of aftershocks in log time, with large magnitude aftershocks as red stars.</figcaption></figure>
<figure><img src="image_16.png" /></figure>

**Extended Data Fig. 6** | B-value plot for earthquakes occurring at the plate interface. Blue dots show the real number of earthquakes above a given magnitude, while the blue line shows the best-fit line obtained with a least-squares inversion assuming a magnitude of completeness M_c of 0.75.
<figure><img src="image_17.png" /><figcaption>**Extended Data Fig. 7** | Family 5 of similar earthquakes, recorded at station H018. Earthquake magnitudes are between 0.74 and 1.29. Data is filtered between 2 and 10 Hz.</figcaption></figure>
<figure><img src="image_18.png" /></figure>

**Extended Data Fig. 8 | Histogram of cross-correlations calculated over the whole signal for all event pairs within a rupture length plus error of each other.** The rupture length is calculated assuming M<sub>w</sub> = M<sub>L</sub> - 0.5 and a stress drop of 3 Mpa.
<figure><img src="image_19.png" /></figure>

**Extended Data Fig. 9** | Semi-variogram of the earthquake elevation relative to the main plane of seismicity. Only earthquakes located within 800 m of the main seismicity plane are considered. Additionally, earthquake pairs are used only if their relative location error is below 25 m and their relative elevation error is below 20 m. A total of 108224 pairs are used.