# Structure of the Ecuadorian forearc from the joint inversion of receiver functions and ambient noise surface waves

Clinton Koch, Colton Lynner, Jonathan Delph, Susan Beck, Anne Meltzer, Yvonne Font, Lillian Soto-Cordero, Mariah Hoskins, Josh Stachnik, Mario Ruiz, et al.

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

## To cite this version:

Clinton Koch, Colton Lynner, Jonathan Delph, Susan Beck, Anne Meltzer, et al.. Structure of the Ecuadorian forearc from the joint inversion of receiver functions and ambient noise surface waves. Geophysical Journal International, 2020, 222 (3), pp.1671-1685. 10.1093/gji/ggaa237. hal-02893205

HAL Id: hal-02893205

https://hal.science/hal-02893205

Submitted on 8 Nov 2021

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Structure of the Ecuadorian forearc from the joint inversion of receiver functions and ambient noise surface waves

Clinton D. Koch<sup>1</sup>, Colton Lynner<sup>1,2</sup>, Jonathan Delph,<sup>3</sup> Susan L. Beck,<sup>1</sup> Anne Meltzer,<sup>4</sup> Yvonne Font,<sup>6</sup> Lillian Soto-Cordero,<sup>4</sup> Mariah Hoskins,<sup>4</sup> Josh C. Stachnik,<sup>4</sup> Mario Ruiz,<sup>5</sup> Alexandra Alvarado,<sup>5</sup> Hans Agurto-Detzel<sup>6</sup>, Philippe Charvis,<sup>6</sup> Marc Regnier<sup>6</sup> and Andreas Rietbrock<sup>7</sup>

Accepted 2020 May 6. Received 2020 March 17; in original form 2019 November 22

## SUMMARY

The Ecuadorian forearc is a complex region of accreted terranes with a history of large megathrust earthquakes. Most recently, a $M_w$ 7.8 megathrust earthquake ruptured the plate boundary offshore of Pedernales, Ecuador on 16 April 2016. Following this event, an international collaboration arranged by the Instituto Geofísico at the Escuela Politécnica Nacional mobilized a rapid deployment of 65 seismic instruments along the Ecuadorian forearc. We combine this new seismic data set with 14 permanent stations from the Ecuadorian national network to better understand how variations in crustal structure relate to regional seismic hazards along the margin. Here, we present receiver function adaptive common conversion point stacks and a shear velocity model derived from the joint inversion of receiver functions and surface wave dispersion data obtained through ambient noise cross-correlations for the upper 50 km of the forearc. Beneath the forearc crust, we observe an eastward dipping slow velocity anomaly we interpret as subducting oceanic crust, which shallows near the projected centre of the subducting Carnegie Ridge. We also observe a strong shallow positive conversion in the Ecuadorian forearc near the Borbon Basin indicating a major discontinuity at a depth of ~7 km. This conversion is not ubiquitous and may be the top of the accreted terranes. We also observe significant north-south changes in shear wave velocity. The velocity changes indicate variations in the accreted terranes and may indicate an increased amount of hydration beneath the Manabí Basin. This change in structure also correlates geographically with the southern rupture limit of multiple high magnitude megathrust earthquakes. The earthquake record along the Ecuadorian trench shows that no event with a $M_w$ >7.4 has ruptured south of ~0.5°S in southern Ecuador or northern Peru. Our observations, along with previous studies, suggest that variations in the forearc crustal structure and subducting oceanic crust may influence the occurrence and spatial distribution of high magnitude seismicity in the region.

**Key words:** South America; Joint Inversion; Crustal Imaging; Seismicity and tectonics; Crustal Structure.

## INTRODUCTION

The $M_w$ 7.8 Pedernales earthquake, which occurred on 16 April 2016, caused significant damage and brought increased attention to the controls on the seismogenic behaviour of the megathrust

along the Ecuadorian margin. The Ecuadorian margin has hosted several large megathrust earthquakes in the past century along the Nazca/South America plate boundary. The oldest significant earthquake recorded is the 1906 $M_w$ 8.3–8.8 event that ruptured a segment nearly 500 km along the trench from southern Colombia to Central

<footnote n="1"><sup>1</sup>Department of Geosciences, University of Arizona, Tucson, AZ 85721, USA. E-mail: clinton.d.koch@gmail.com</footnote>

<footnote n="2"><sup>2</sup>Department of Earth Sciences, University of Delaware, Newark, DE 19716, USA</footnote>

<footnote n="3"><sup>3</sup>Department of Earth Sciences, University of Oregon, Eugene, OR 97403, USA</footnote>

<footnote n="4"><sup>4</sup>Department of Earth and Environmental Sciences, Lehigh University, Bethlehem, PA 18015, USA</footnote>

<footnote n="5"><sup>5</sup>Instituto Geofísico, Escuela Politécnica Nacional, Quito 170517, Ecuador</footnote>

<footnote n="6"><sup>6</sup>Geoazur, Université de Nice Sophia-Antipolis, IRD-CNRS-OCA, 06100 Nice, France</footnote>

<footnote n="7"><sup>7</sup>Karlsruhe Institute of Technology, Geophysical Institute, 76187 Karlsruhe, Germany</footnote>
Ecuador (Kanamori & McNally 1982; Yamanaka et al. 2017). Earthquakes with $M_w$ 7.6 and greater have occurred in 1942, 1958, 1979 and 2016 (Beck & Ruff 1984; Swenson & Beck 1999; Ye et al. 2016; Nocquet et al. 2017; Yamanaka et al. 2017). Collectively, these earthquakes have ruptured the majority of 1906 earthquake rupture zone. South of the southern extent of the 1906, 1942 and 2016 earthquakes (near 0.5°S), the occurrence of $M_w \ge 7.4$ earthquakes are absent throughout our study area until at least 3.5°S. The lack of high magnitude earthquakes between 0.5°S and 3.5°S suggests segmentation in megathrust behaviour along the margin. This segmentation has been previously attributed to the presence of the Carnegie Ridge (CR) acting as a barrier towards southward rupture of large earthquakes (Fig. 1; Collot et al. 2004, Gailler et al. 2007). In the area surrounding the CR near the La Plata Island, several slow slip events (SSE) have occurred, including an event in 2016 following the Pedernales Earthquake (Vallee et al. 2013; Vaca et al. 2018). SSEs have also been observed to the north in the Punta Galera-Mompiche Zone (PMGZ; Fig. 1) along with a deep SSE in 2016 (Fig. 1; Rolandone et al. 2018; Vaca et al. 2018).

The lack of high magnitude megathrust earthquakes south of the subduction of the CR appears to extend to at least 3.5°S, where a $M_w$ 7.4 earthquake occurred in 1953 that was likely not on the plate interface and therefore would not have ruptured much of the plate boundary (International Seismological Centre 2016). Recent GPS studies have suggested that plate coupling is low south of 0.5°S except for a small segment of the plate boundary near La Plata Island where SSEs have been observed (Nocquet et al. 2014), and is perhaps not capable of producing a great earthquake. It is important to understand the along-strike differences in the seismogenic zone to better characterize the slip behaviour seen throughout the region. Many factors, including topography on the downgoing plate, trench sediments, fluids, and properties of the upper plate, contribute to the rheological properties of the seismogenic zone (Bilek & Lay 2018, and references therein). In this study, we investigate the discontinuity and shear wave velocity structure of the upper plate along the Ecuadorian forearc to better understand how along-strike variability might contribute to variable slip behaviour along the margin.

The Ecuadorian forearc consists of several accreted oceanic terranes that were assembled between the Late Cretaceous and Early Eocene (Kerr et al. 2002; Jaillard et al. 2009). Many of the details related to the age of formation and accretion, and distribution of these terranes, comes from surface geological and geochemical data with minimal subsurface data (<5 km) available from oil exploration wells (Evans & Whittaker 1982). Previous geophysical studies, including Bouguer gravity analysis and earthquake relocation studies, have been used to constrain the subsurface geology (Feininger & Seguin 1983; Font et al. 2013). Numerous marine seismic studies have illuminated the offshore forearc structure (Collot et al. 2004, 2017; Graindorge et al. 2004; Gailler et al. 2007). These studies have provided a good first order understanding of the geology of the forearc in the offshore region. They have not, however, been able to elucidate the onshore portion of the forearc, and as a result have provided few insights on the role of the overriding crust on seismic behaviour along the margin.

Following the $M_w$ 7.8 Pedernales earthquake, an international collaboration arranged by the Instituto Geofísico at the Escuela Politécnica Nacional (IG-EPN) mobilized a rapid deployment of 55 broad-band, intermediate, and short-period seismometers and 10 ocean-bottom seismometers in order to capture the aftershock sequence for ~1 yr (Meltzer et al. 2019). This data provides an unprecedented opportunity to image Ecuadorian forearc structure

and study how it relates to the seismogenic zone and the range of slip behaviours observed at the plate boundary. To constrain forearc structure, we perform receiver function (RF) analyses at all available forearc stations and then combine the RF results with Rayleigh-wave dispersion measurements derived from ambient noise cross correlations (Lynner et al. 2020) in a joint inversion. By jointly inverting RFs and surface wave dispersion data, we produce a detailed shear wave velocity model that is sensitive to both vertical and lateral velocity variations. We then relate structures present in our shear wave velocity model to the seismogenic properties along the margin.

## TECTONIC BACKGROUND OF THE FOREARC

From the Western Cordillera to the present-day trench, Ecuador is made up of accreted oceanic plateau and arc rocks. Details of the timing and origin of the accreted terranes have been the subject of much debate, but the most recent evidence suggests that the oceanic terranes share a similar origin and were part of the Colombian Caribbean Oceanic Plateau (CCOP, Luzieux et al. 2006; Jaillard et al. 2009). The CCOP collided with South America in the Late Cretaceous and failed to subduct due to having a higher buoyancy than typical oceanic crust (Cloos 1993). The subsequent accretion onto the South American plate may have already been comprised of the amalgam of tectonic blocks that make up the Ecuadorian forearc, or the remainder of the blocks may have accreted in subsequent events (Luzieux et al. 2006; Jaillard et al. 2009). Regardless of the exact timing, the forearc of Ecuador is thought to largely be accreted oceanic lithosphere material that is heavily obscured by sedimentary basins.

Prior to the accretion of the forearc units, the Pallatanga block accreted onto South America and is currently part of the Western Cordillera (Fig. 1; Kerr et al. 2002). The Pallatanga block is bounded to the east by the Pujili Fault and the Inter-Andean valley and to the west by the Chimbo-Touchi shear zone. West of the Chimbo-Touchi shear zone lies the Macuchi Terrane, an Eocene-aged arc that formed by eastward subduction after the accretion of a fragment of the CCOP on to South America (Hughes & Pilatasig 2002, Vallejo et al. 2009). In the forearc, sedimentary cover obscures the underlying oceanic terranes, with a few exceptions of interpreted oceanic arc and plateau rocks outcropping in the north and south of our study area. Topography of the forearc is generally low with elevations reaching up to 650 m in a few isolated locations along the coastal cordillera. In northern coastal Ecuador, near the Pedernales earthquake epicentre and further south near the Chongón-Colonche Fault, outcroppings of interpreted oceanic plateau material are found (Fig. 1). These units are often referred to as the Piñon Formation. The Piñon Formation is made up of ultramafic, mafic and intermediate magmatic rocks (basalts, andesites, gabbros, dolerites, wherlites and peridotites; Luzieux 2007). In central Ecuador along the coast, the San Lorenzo Formation is in contact with the Piñon Formation (Luzieux et al. 2006). The San Lorenzo Formation is composed of island arc affinity rocks of a late Cretaceous age. The age and composition of the San Lorenzo Formation is similar to that of arc rocks in the Pallatanga block suggesting its formation was the result of westward driven subduction of South America beneath the CCOP (Luzieux et al. 2006). Together the Piñon and San Lorenzo Formations, as well as several other minor formations, comprise the Piñon Terrane. The Piñon Terrane is often mapped as underlaying the majority of the Ecuadorian forearc (Kerr et al. 2002; Jaillard et al. 2009).
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Left-hand panel: topographic map with seismic stations (triangles) deployed during Ecuador RAMP temporary deployment (networks XE and 8 G; white and blue), and permanent seismic stations (network EC; green). Blue stations' results are highlighted in Fig. 2. Inset map in upper left corner. Dashed black ellipse shows approximate rupture area of the 1906 earthquake (Chlieh *et al. 2014*). Dashed yellow outlines show the approximate aftershock areas of the 1942 and the 1958 earthquakes (Mendoza & Dewey 1984; Swenson & Beck 1996). Red polygon is the high slip (>1 m) area for the 2016 Pedernales earthquake; red star demarks the epicentre (Nocquet *et al. 2017*). Note the lack of major earthquakes south of 0.5°S. Slow slip event regions are shown as purple ellipses including the Punta Galera-Mompiche Zone (PMGZ; Vallee *et al. 2013*; Vaca *et al. 2018*; Segovia *et al. 2018*; Rolandone *et al. 2018*). -2200 m bathymetric contour outlines the Carnegie Ridge. Slab contours are from Slab2 (Hayes *et al. 2018*). Relative Nazca to fixed South America Plate Motion (PM; 55 mm yr⁻¹) is from Chlieh *et al. (2014)*. Right-hand panel: map of the major tectonic blocks, accreted terranes and basins in Ecuador significant to this study (Luzieux *et al. 2006*; Jaillard *et al. 2009*). Yellow lines: major fault systems; CCF: Chongón-Colonche Fault; JPF: Jipijapa Fault; JF: Jama Fault system; PF: Pujili Fault, CTSZ: Chimbo-Touchi Shear Zone. Solid black lines show the locations of cross sections shown in Figs 3, 4 and 6.</figcaption></figure>

## DATA AND METHODS

Following the M<sub>w</sub> 7.8 16 April 2016, Pedernales Earthquake, an international response coordinated by the IG-EPN, deployed seismic instruments immediately following the main event (Meltzer *et al. 2019*; Fig. 1). These instruments recorded continuously for approximately one year following installation. Data from this deployment has already been used to study the aftershock seismicity (Agurto-Detzel *et al. 2019*; Meltzer *et al. 2019*; Soto-Cordero *et al. 2020*) and velocity structure across the rupture area and Ecuador (León-Ríos *et al. 2019*; Lynner *et al. 2020*). We incorporate data from the temporary deployment as well as permanent stations from the Ecuadorian national network (Alvarado *et al. 2018*) to construct P-wave receiver functions (RFs) along the margin. We then jointly invert our calculated RFs with Rayleigh wave phase velocity measurements from Lynner *et al. (2020*) to develop the best constrained shear wave velocity model of the Ecuadorian forearc to date.

### Receiver function quality control

RF analysis is a commonly used technique to isolate P-to-S conversions from discontinuities beneath 3-component seismic stations (Langston 1979; Ammon 1991). The vertical component is deconvolved from the radial and tangential components to produce the

radial and tangential receiver functions respectively. In this study, we focus only on radial receiver functions.

RFs are calculated using data from 64 (50 rapid deployment and 14 permanent network instruments) broadband and short-period stations in the Ecuadorian forearc (Fig. 1). Teleseismic earthquakes from distances of 30°–90° with M<sub>w</sub> 5.9 and greater are used. The resulting backazimuthal event distribution is heavily concentrated to the southeast and northwest, with the majority of earthquakes coming from the Chilean and Aleutian subduction zones. Data from short-period stations had their instrument responses removed prior to calculating RFs in order to broaden the frequency band (Niu *et al. 2005*). Instrument responses for broadband and intermediate instruments are not removed, as the response of these instruments is already flat within our frequency range of interest. For each station-event pair, we cut and filter the data using a bandpass filter between 0.04 and 4 Hz. Initial quality control is performed manually for each event-station pair by verifying the presence of a strong initial P-wave arrival. If no arrival was visible, the event is discarded from further analysis. RFs are then calculated using the time-domain iterative deconvolution technique (Ligorria & Ammon *et al. 1999*) with Gaussian values of 5.0 and 2.8 (equivalent to a low pass filter of ~2.5 and 1 Hz, respectively). A second round of quality control was conducted on the RFs, keeping only RFs that show a positive initial peak (associated with the incoming P wave), and
contained RF amplitudes < 1. Amplitudes higher than 1 for the frequency content used in our study result from significant vertical-to-radial energy transfer due to complex shallow structure, such as anisotropy or dipping layers. RFs with significant oscillatory behaviour are also removed. After the quality control steps, 609 RFs remained.

## Adaptive common conversion point stacking

Adaptive Common Conversion Point (ACCP) stacks (Delph et al. 2015, 2017) are constructed from the radial RFs that passed our quality controls using a modified version of the Funclab package (Eagar & Fouch 2012; Porritt & Miller 2018). To create our ACCP stacking volume, RFs are ray-traced back along their theoretical ray path using a 1-D velocity model (Fig. S1) averaged from an *a priori* 3-D $V_p$ model of the Ecuadorian forarc with a $V_p/V_s$ of 1.77 (Font et al. 2013). Timing corrections are calculated based on the full 3-D velocity model following Eagar & Fouch (2012). These rays are then mapped to a 3-D grid and converted to depth. Their associated amplitudes are normalized to the highest value for each RF and linearly stacked (Dueker & Sheehan 1997). Our ACCP stacks have laterally sized bins starting at 0.1° in diameter spaced every 0.05° laterally and a 1 km thickness to a depth of 100 km. We allow the bins to laterally expand to a maximum of 0.2° if less than 5 RFs are present within a bin.

## Joint inversion

In order to obtain a high-resolution shear wave velocity model, we jointly invert the complementary data sets of RFs and ambient noise-generated Rayleigh wave dispersion data. Rayleigh wave dispersion data is able to recover the absolute velocities of vertically polarized shear waves, however, the broad sensitivity kernels of Rayleigh waves at different frequencies result in a weak sensitivity to velocity discontinuities. RFs, conversely, provide strong constraints on discontinuity structure, but are not sensitive to the absolute velocities on either side of these boundaries. By combining these data sets, we are able to reduce the non-uniqueness of the inverse problem by mitigating the shortcomings of each individual method (Julia et al. 2000) and produce a high-resolution shear wave velocity model that is sensitive to both vertical and lateral shear wave velocity variations.

We use the Rayleigh wave dispersion measurements from Lynner et al. (2020), who used intermediate and broadband seismic data from Ecuador to measure phase velocity dispersion throughout the forearc and arc. They extracted dispersion curves from cross-correlation of day-long records between all contemporaneously recording stations in the period range of 8–40 s. The interstation phase velocity measurements were then inverted for phase velocity maps on a 0.1° x 0.1° grid. From these phase velocity maps, we extract dispersion curves at each grid point in the model (examples in Fig. S2). Further details of the processing and inversion of the ambient noise data can be found in Lynner et al. (2020, and references therein).

A second common conversion point (CCP) volume is computed with 0.1° spatial sampling that matches that of the dispersion data to 100 km in depth so that we can extract and pair RFs from a gridded data set with the dispersion data. We use a bin spacing of 0.1° with bins 0.2° in diameter and a depth spacing of 0.5 km. Unlike the ACCP stacks described above, the bins in this volume are fixed at 0.2°. We use the same 1-D average velocity model

(Fig. S1) for the forearc region derived from Font et al. (2013) to map receiver functions to depth. The CCP stacks are constructed using non-normalized RFs with a Gaussian value of 2.8 in order to capture discontinuities on the order of ~1 km. From this CCP volume, RFs are extracted on a 0.1° x 0.1° grid to match the grid spacing of the dispersion data, and converted back to time using the same 1-D velocity model used to map RFs to depth. Unlike in the previously described ACCP stacks, we do not perform timing corrections to account for 3-D heterogeneity, as this would convolute which 1-D velocity structure would be most accurate for migrating RFs back to time. Extracting RFs from a CCP volume rather than using individual station RFs has the advantage of more accurately accounting for variations in backazimuth and ray parameter, thereby minimizing artefacts associated with imperfect data coverage and lateral heterogeneity beneath a station. CCP stacks can, however, suffer from unrealistic sharp changes in amplitude that may result from: (1) the discretization of the model into constant thickness bins, which will incorporate different time windows of the receiver functions based on the velocity model or (2) introducing new data into a bin, as rays spread out and number of receiver functions in a bin generally increases as a function of depth. These artefacts cannot be predicted by the receiver function forward operator. In order to mitigate these artefacts, the CCP-derived RFs are filtered with a Gaussian filter following the method described by Delph et al. (2017). This results in RFs with a slightly lower frequency (2.5 Gaussian alpha parameter), leading to a small loss of vertical resolution, which can then be paired with dispersion data.

The dispersion curves and the effective 2.5 Gaussian CCP-derived RFs are jointly inverted to construct a high-resolution shear wave velocity model (following Delph et al. 2017). We use an initial velocity model consisting of a 4.5 km s⁻¹ half-space discretized into 1 km thick layers. By assuming a simple starting velocity structure, we ensure that the features in the resulting shear wave velocity model are driven by the RFs and dispersion data and are not biased by *a priori* assumptions of velocity structure in the starting model. We prefer using a uniform starting model so that discontinuities, albeit smoothed, are placed at the depth required by the RF data. Using a starting model containing discontinuities severely biases the resulting shear wave velocity model, because discontinuities in the starting model remain fixed in depth with this linearized inversion technique. The joint inversion algorithm starts with the initial velocity model, computes predicted RF and dispersion curves at each grid point, then compares predicted and observed data to calculate a misfit. The misfit is used to iteratively update the velocity model. This process is repeated for 40 iterations or until the model change between iterations becomes negligible (<0.05 per cent). The model extends to 100 km to ensure the dispersion sensitivities are near 0 at the base of our model to avoid biases from boundary effects at the bottom of the model. The joint inversion requires a weighting factor between the RFs and dispersion curves in order to define the relative importance of each data set. After testing several weighting factors, we use a weighting factor of 0.3, indicating 30 per cent weight on the dispersion data and 70 per cent on the RFs, for the final model reported in this study. Higher weightings of the dispersion data produce minor improvements in the fit of the dispersion data but significantly degrade the RF fits and associated sensitivities to discontinuity structure (Fig. S3). We also run this inversion with different uniform velocity starting models (Figs S4 and S5) and a weighting factor of 0.7 (Figs S4 and S5) to test the effects of these parameters on our final model.
RESULTS

We show results of the radial RFs for some representative stations, the Adaptive Common Conversion Point stacks, and the shear wave velocity model resulting from the joint inversion of RF and Rayleigh wave dispersion data. We identify four along-strike segments (from north to south segments 1, 2, 3 and 4) based on the changes in seismic characteristics determined from the RF and variations in shear wave velocity structure.

Receiver functions

Fig. 2 shows examples of the 5.0 Gaussian RFs at representative stations in our study (labelled on Fig. 1). Station EC16 is located within the 2016 earthquake rupture region (Nocquet et al. 2017). This station shows two primary positive conversions (Ps₁ and Ps₂), at 1.4 and 2.2 s followed by a negative conversion, at ~3 s (Ps₃). In RFs, a positive amplitude P-to-S conversion results from an increase in velocity with depth at a discontinuity and a negative amplitude conversion results from a decrease in velocity with depth at a discontinuity. At around 4 and 5 s potential multiples of Ps₁ and Ps₂ are present (PpPs₁ and PpPs₂). Station EC14 is from the coastal area north of the rupture zone of the 2016 event and within the Borbon Basin. A series of positive peaks (Ps₁, Ps₂ and Ps₃) can be seen between 1 and 3 s and broad negative peaks (Ps₄) at ~4 s. A final peak at ~6.5 s (PpPs₂) may be a multiple of the Ps₂. Assuming a Vp/Vs of ~1.77 and a Vs of ~3.6 km s⁻¹, a primary P-to-S conversion at 6.5 s would map to a discontinuity at a depth of ~50 km but if the conversion is a multiple (PpPs₁) it would be generated at a discontinuity at ~16 km.

Stations EC08 and EC09 are farther from the coast, located within the Manabí Basin, and further demonstrate the complexity of the RFs in the basin area. EC08 has two strongly positive arrivals, Ps₁ at ~2.5 and Ps₂ at ~5 s that are consistent across most ray parameters (Fig. 2). The Ps₂ conversion is too close in time to be a multiple of Ps₁, indicating that it is most likely a primary conversion. The likely multiple of Ps₁ is seen at a time of ~6.5 s (PpPs₁). The negative multiple may also be present slightly later at 7 s, indicated by PsPs₁. The strength and time variability, however, suggests that there may also be a primary negative conversion at this time. EC09 has a very weak initial P arrival followed by a series of primary positive arrivals (Ps₁, Ps₂ and Ps₃) with a large peak at ~3 s. Other studies have shown similar complications in RFs associated with basin structures (Ma & Clayton 2016). At around 5.8 s, a negative conversion (Ps₄) is present that corresponds to a depth of ~43–50 km depending on the velocity model. Beneath this negative, a positive peak (PpPs₂) may be the multiple of Ps₂.

Stations in the south were placed nearer to or on solid bedrock rather than basin deposits and generally produced higher quality RFs. For example, the RFs for station VMON show a positive peak (Ps₁) at ~1 s and another (PpPs₁) at ~4 s. This second peak is likely a multiple of the first. Assuming a Vp/Vs of ~1.77 and a starting Vs of 3.0 km s⁻¹, the resultant multiple (PpPs₁) would be generated at a discontinuity at ~8 km. Station CAYO is located on the southern side of the Manta peninsula and has a positive peak (Ps₁) consistent across all ray parameters at ~1.4 s. This is followed by a negative peak around 2.2 s (Ps₂) and a positive peak after Ps₂ at ~3 s (Ps₃). Finally, a weak multiple of Ps₁, PpPs₁ may be present around ~4.3 s.

If the multiples can be identified correctly, they can constrain the Vp/Vs above the conversion the causes the multiple. Using the primaries and multiples we identify for the stations in Fig. 2, we approximate the Vp/Vs of using equations from Zandt et al. (1995)

and a range of Vp between 5.75 and 6.75 km s⁻¹. We note that the complexity of the tectonic structure in the region makes identifying multiples difficult; thus, not all primaries have associated multiples to use in the Vp/Vs calculation. For station EC16, Ps₁ has a Vp/Vs of between 1.86 and 1.94. Ps₂ in EC14 indicates Vp/Vs values between 1.71 and 1.77. Station EC08 indicates a Vp/Vs range of 1.70–1.77 for Ps₁ and 1.83–1.89 for Ps₂. Finally, the multiple for Ps₁ in stations VMON and CAYO indicates values between 1.75 and 1.83. These findings suggest that the Vp/Vs is strongly variable throughout the upper crust, which is not unexpected given the presence of basins and accreted terranes in the upper crust.

Adaptive common conversion point stacks

Trench parallel and trench perpendicular cross sections through the Gaussian 5.0 ACCP stacks are shown in Fig. 3. Several distinct features can be seen throughout the cross sections. The trench parallel profiles are divided into four segments based on the location and character of the shallow (<25 km) P-to-S conversions (Fig. 3; A–A' and B–B'). In segment 1 (S1), we observe two prominent primary positive conversions are at depths of ~7 and 17 km (D₁ and D₂; Fig. 3). In cross section C–C', discontinuity D₁ is present at ~5 km in the first ~40 km of the profile and again at 5–7 km on the eastern edge of the profile. Here, D₁ may be related to an increase in velocity at the base of the Borbon Basin in the west, and the Manabí basin in the east. Multiples related to the shallow basin structure make deeper interpretation difficult in this area. Where A–A' intersects the profile C–C' (Fig. 3), D₁ and D₂ are present and continue for ~20 km to the east. However, further to the east (~75 km), D₁ has shallowed and is no longer apparent in the cross section. D₂ also continues to weaken eastward.

In the segment 2 (S2) the D₂ conversion deepens to ~20 km and is overlain by a weak positive conversion in A–A' and a weak negative conversion in B–B'. In D–D', D₂ is present throughout the section, splitting into two conversions in the eastern 30 km of the profile (Fig. 3).

Between S2 and segment 3 (S3) D₂ shallows, reaching a depth of ~11 km in S3. As D₂ shallows another positive conversion, D₃, emerges, underlaying and paralleling D₂ for ~80 km at depths a depth of ~22 km in A–A' and ~25 km in B–B' (Fig. 3).

In segment 4 (S4) D₂ continues to shallow to a depth of ~7 km and is overlain by two more positive conversions at depths of 10 and 15 km and a negative conversion (N₁) at ~28 km in A–A' and B–B'. In E–E' a D₂ is continuous across the section, broadening and deepening further east (Fig. 3).

Transecting several of the segments at depths between ~22 and 35 km, a broad negative conversion, N₁, is observed throughout S1, S3, and S4 in A–A' and S3 and S4 in B–B' (Fig. 3). In the trench perpendicular sections N₁ is present in C–C' and D–D' near the western ends of the profiles above and below Slab2 respectively, and E–E' as a more continuous discontinuity that dips eastward. This discontinuity (N₁) closely parallels Slab2 (Hayes et al. 2018) in S1 (A–A') from the northern edge of the profile to the S1–S2 boundary and in S4 (E–E') for ~65 km from the western edge of the profile eastward.

Receiver function multiples

Multiples are a common concern in receiver function studies in areas with strong, shallow discontinuities. The arrival times of multiples from primary conversions at shallow depths can interfere with those of primary conversions from deeper boundaries or become misinterpreted as primary conversions. In our results, the shallowest primary
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Moveout plots for the 5.0 Gaussian radial receiver functions for stations EC16, EC14, EC08, EC09, CAYO, VMON (Fig. 1; blue triangles). Receiver functions are stacked in epicentral distance bins of 5°. Black lines delineate primary *P*-to-*S* conversions and dashed green lines delineate likely multiples. Primary conversions denoted by *Ps*<sub>x</sub>, and multiples by *PpPs*<sub>x</sub> and *PsPs*<sub>x</sub>, with subscripts indicating different conversions.</figcaption></figure>

conversions are seen in segments S1 and S4. Multiples from these conversions may impact the interpretation of features around depths of ~25 km and deeper in the northern section of A–A’ and in the southern region of A–A’ and B–B’ (Fig. S6). Other multiples are generally at or below the top of the subducting slab, placing them deeper than any features we focus on in this study.

## Joint inversion

Cross sections through the upper 50 km of our shear wave velocity model from our joint inversion at the same locations as the ACCP

stacks are shown in Fig. 4. As in the ACCP stacks, we separate each profile into four along-strike segments based on the crustal shear wave velocities and RF discontinuities. Side by side comparisons between the ACCP stacks and shear wave velocity model can be seen in Figs S7 and S8. Fits of the predicted data compared against the observations are shown in Fig. 5. In general, the final model fits the dispersion curve data well across the study area, while RF fits are degraded in regions with complicated shallow structure, such as around the Manabi and Borbon basins (Fig. 5). The fit for individual gridpoints at each of the stations shown in Fig. 2 can be seen in Fig. S2. In general, we observe good fits for periods 25 s
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Trench parallel cross sections, A–A' and B–B', and trench perpendicular cross sections, C–C', D–D' and E–E', through the 5.0 Gaussian ACCP stack. Dashed black lines are interpreted primary conversions from receiver functions. $D_x$ indicates a positive conversion $N_x$ indicates a negative conversion. Pink line is the top of the subducting oceanic crust from Slab2 (Hayes *et al.* 2018). Black dots are earthquakes from the Pedernales aftershock catalog shown only for trench perpendicular sections projected 10 km perpendicular to the section (Agurto-Detzel 2019). Approximate location of the Carnegie Ridge and the 2016 Pedernales high slip (>1 m) rupture region (A–A' only) are projected in the direction of plate motion. Velocity model and other parameters used for ACCP stacks described in text. Location map (bottom right) shows locations of cross sections and our interpretation of segments (S1, S2, S3 and S4) separated by red lines.</figcaption></figure>

and shorter, with worse fits for longer periods between 30 and 40 s, corresponding to depths of ~40–50 km. Figs S9 and S10 show the receiver function and dispersion fits along the cross sections seen in Fig. 4 along with corresponding cross sections through the ambient noise tomography (ANT) results from Lynner et al. (2020).

At a depth of ~20 km, the results from north to south show a moderately fast velocity region (S1, ~4.0 km s⁻¹), a slow velocity region (S2, <3.5 km s⁻¹) and fast velocity regions (segments 3 and 4. >4.4 km s⁻¹; Fig. 4). S1 extends from the northern edge of each

profile to ~140 km on cross-section A–A' and ~55 km on cross-section B–B'. The observed moderately fast seismic velocities range from ~10 to 25 km depth with shear wave velocities of up to ~4.2 km s⁻¹. Section C–C' (Fig. 4) transects this segment and the moderate velocities appear to extend to at least the western edge of the section.

In each profile, S2 shows slow shear wave velocities are present in the upper ~30 km. These velocities range from 2.8 to 3.5 km s⁻¹ and generally increase with depth. The shallowest portion of this slow
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Trench parallel cross sections, A–A' and B–B', and trench perpendicular cross sections, C–C', D–D' and E–E', through the shear wave velocity model as determined from the joint inversion. Contour interval is 0.2 km s⁻¹ and bolded contours are at 2.4, 3.0 3.6 and 4.2 km s⁻¹. Pink line is the top of the subducting oceanic crust from Slab2 (Hayes *et al.* 2018). Black dots are earthquakes from the Pedernales aftershock catalogue shown only for trench perpendicular sections projected from 10 km perpendicular to the section (Agurto-Detzel 2019). Approximate location of the Carnegie Ridge and the 2016 Pedernales high slip (>1 m) rupture region (A–A' only) are projected in the direction of plate motion. Location map (bottom right) shows locations of cross sections and our interpretation of segments (S1, S2, S3 and S4) separated by red lines.</figcaption></figure>

anomaly may be, in part, related to the Manabí basin. However, the lateral extent and depth make it unlikely that the entire slow velocity anomaly can be attributed to the basin. Furthermore, the upper 5 km have poor velocity resolution due to the low sensitivity at shallow depths of the dispersion data at the periods used in this study (Lynner *et al.* 2020). Section D–D' (Fig. 4) shows a thick slow velocity zone that thins eastward in the crust.

In S3, the slow velocities shallow to the upper ~8 km and a section of moderate shear wave velocities between 3.6 and 4.0 km s⁻¹ is present at depths between 8 and 20 km in sections A–A' and B–B'.

Beneath this, a sharp boundary at ~20 km is observed, and the velocities increase to ~4.4 km s⁻¹. The shear wave velocities below this anomaly begin to decrease at a depth of 35 km in the vicinity of Slab2 (Hayes *et al.* 2018).

In S4, the moderate shallow velocities in S3 are not present and the fast anomaly present in S3 shallows. In S4 this fast anomaly reaches shear wave velocities up to 4.5 km s⁻¹ at depths as shallow as ~10 km in both cross-sections A–A' and B–B' (Fig. 4). The transition between the S2 slow anomaly and the S3 and S4 fast anomaly shows a rapid lateral increase in velocity over ~20 km.
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Percent fit of our final joint inversion shear velocity model to the Rayleigh wave dispersion (left-hand panel) and receiver function (right-hand panel) input data.</figcaption></figure>

Cross-section E–E' through S4 shows that the fast velocities do not extend all the way to the coast, but are pervasive above the subducting oceanic crust in the east.

In cross-section A–A' there is a ~10–20 km thick slow (3.6 km s⁻¹) velocity anomaly that varies in depth along cross section A–A' (Fig. 4). In S1, the centre of this anomaly is at a depth of ~28–32 km and the base of the anomaly closely aligns with Slab2 (Hayes et al. 2018). In S2 this anomaly deviates from Slab2, rising to a depth of ~18 km and in S3 and S4 deepens again realigning with the top of Slab2 (Fig. 4). In B–B', the slow shear wave anomaly, and Slab2, are only present in S3 and S4. This feature is not imposed by the starting velocity model and is not clearly observed in ambient noise tomography even when starting models that impose a slab are used (supplemental info in Lynner et al. 2020). Thus, the anomaly is the result of the inclusion of the RFs in the joint inversion. Multiples from shallower conversions (D₁) may be present at these depths, especially in the S1 and S4 (Fig. S6), which reduces the confidence of any interpretation. Figs S4 and S5 show the resulting shear wave velocity models when periods greater than 25 s are removed. This shows that the low velocity feature at ~40 km in S2 is not dependent on poor fits to the longer periods in the dispersion data. The continuity, velocity, geometry, and close alignment with Slab2 away from the CR suggest this may be a real feature associated with the slab.

## Comparison with ambient noise tomography

The inclusion of RFs into the inversion produce more detailed shear wave velocity images than the inversion of dispersion data from ambient noise tomography alone, as would be expected due to the lack of constraints on discontinuity structure without the receiver functions. Figs S9 and S10 show direct comparisons along the cross sections presented here between the joint inversion and ANT-derived shear wave velocity models (Lynner et al. 2020). In S1 in cross-section A–A' the fast velocities between 10 and 25 km are not nearly as prominent in the ANT-only model. In cross-section C–C' the ANT-only results show a deep slow velocity feature at a distance of ~60 km from the start of the profile. In S2 both models (joint inversion and ANT) indicate slow velocities, while the boundaries of these features are much more clearly defined in the joint inversion results. Similarly, in S3, the sharp velocity gradients observed at ~10, 20 and 35 km in the joint inversion are much smoother in

the ANT-only model as expected. Finally, the fast velocities in S3 and S4 are not as prominent in the ANT-only model along these sections, only reaching velocities of ~4.3 km s⁻¹.

## DISCUSSION

In this section, we will discuss our interpretations of the ACCP stack and shear wave velocity model determined from the joint inversion in the context of the tectonic units and distribution of megathrust seismicity (Fig. 6).

### Terrane distribution in the forearc crust

The complexity of the Ecuadorian forearc appears to be reflected in the strong changes in shear wave velocity between the segments shown in the cross-sections in Figs 3 and 4, which suggest significant spatial variation in the composition and/or fluid content of the accreted terranes from north to south.

### Northern segment 1

In the upper ~25 km of A–A' and upper ~30 km of B–B', the velocities in the S1 increase rapidly, reaching 3.8–4.2 km s⁻¹ by 10–15 km depth. This northern segment extends south roughly to the location where the Piñon Formation outcrops north of the Jama Fault System (Fig. 1). In the A-A' line, the ACCP stacks show that conversion D₁ extends to roughly the same area and is likely the top of the high velocity feature. Given the tectonic history, it is likely that the crust in the S1 consists primarily of accreted oceanic plateau material as mapped in the Piñon Terrane, despite the fact that the Piñon Formation does not outcrop at the surface in S1. The velocities are consistent with this interpretation when compared to accreted oceanic plateaus elsewhere. The Siletzia Terrane, an accreted oceanic plateau in the Cascadia forearc, similarly shows shear wave velocities in the range of 3.6–3.9 km s⁻¹ (Rathnayaka & Gao 2017; Delph et al. 2018). A compilation of oceanic plateau studies shows that, in the deeper layers, P-wave velocities range from 6.5 to 7.0 km s⁻¹ (3.65–3.97 km s⁻¹ assuming a Vp/Vs = 1.77) and some oceanic plateaus contain a basal unit with P-wave velocities of up to 7.9 km s⁻¹ (Ridley & Richards 2010). We are not able to separate these layers distinctly, likely due to the strong deformation that occurs during the accretionary process (Tetreault & Buiter 2014), and/or because of the broad shear wave sensitivity kernels the Rayleigh waves. Discontinuity D₂ in A–A', however, suggests a strong velocity change that appears to be internal to the fast velocity body. At the base of the Pinon Formation, the model rapidly transitions into slower seismic velocities that may represent the subducting oceanic crust (Fig. 6).

### Central segment 2

The Ecuadorian forearc crust in S2 has slow velocities, ranging from 2.8 to 3.6 km s⁻¹, setting it apart from the northern and the two southern segments. The slow velocities are thickest towards the west, extending to a depth of ~30 km (D–D'; Fig. 4). The slow velocities transition rapidly to faster velocities in the north and south at the same depths. This rapid lateral transition is seen particularly well in B–B' (Fig. 4) where the slow velocities increase by ~0.6 km s⁻¹ over ~20–25 km. In A–A' S2 lines up with much of the rupture of the 2016 event and the deepest part of the slow velocities appear to be bounded in the north by the Jama Fault
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Interpreted cross sections A–A' and B–B' through the joint inversion shear wave velocity model (upper) and the 5.0 Gaussian ACCP stack (lower). The thin dashed black line outlines the approximate location of the Manabí Basin. Approximate location of the Carnegie Ridge and the 2016 Pedernales high slip (>1 m) rupture region (A–A' only) are projected in the direction of plate motion. Approximate segment boundaries are shown as thick dashed lines on each cross section.</figcaption></figure>

System (Fig. 4). The base of the slow anomaly is in close proximity to the velocities interpreted as the subducting oceanic crust in S1 (Fig. 6), but is ~20 km shallower than predicted by Slab2 (Hayes et al. 2018). At the surface, the Manabí Basin corresponds with the slow velocities (2.8–3.5 km s⁻¹) and the basin may account for some of the slow velocities in the upper ~10 km. The depth of the Manabí Basin is not well constrained, but it has been suggested to be at most 10 km thick (Jaillard et al. 1995; Deniaud 2000). As such, it is likely that the velocities seen beneath the Manabí Basin are due to accreted material that is slower than the faster forearc material in S1, S3 and S4. The sharpness of the northern and southern boundaries suggests a more abrupt transition, similar to what we would expect at a boundary between accreted terranes.

The sediment cover and lack of a surface expression of this deeper, slow feature beneath the Manabí Basin make it difficult to constrain. However, there are several possibilities that may have resulted in the observed slow velocity structure. Four possibilities discussed here include (1) fluids resulting from dehydration of the subducting slab infiltrating and altering the crust, (2) a zone of lithosphere that was altered prior to accretion, (3) the accretion and preservation of the upper layers of the Piñon Formation and (4) tectonic erosion of the forearc crust as the buoyant Carnegie Ridge subducts.

This slow region occurs above the projected centre of the Carnegie Ridge (CR, Fig. 4). As the CR subducts and heats up,

fluids are likely released through dehydration reactions within hydrated mineral assemblages in the subducted oceanic crust (Peacock 1989). The subducted structure of the CR may allow for pressure-driven flow of expelled fluids towards the centre of the ridge, causing fluids to accumulate in the lithosphere beneath the Manabí Basin. Subducting fluids have the potential to alter the overriding material and produce slower seismic velocities (Peacock 1989; Bostock et al. 2002; Hyndman & Peacock 2003). If the material beneath the Manabí Basin was originally that of the Pinon Terrane, which is thought to be comprised of mafic and ultramafic rocks, any hydration would lead to serpentinization and thus slower velocities.

An alternative reason for the slower seismic material is that this was a region of increased alteration prior to accretion onto the forearc. The presence of a serpentine-rich body could have developed via hydration and alteration of an ultramafic sliver of oceanic crust prior to accretion onto the margin (Cluzel et al. 2001; Guillot et al. 2015). Similar accreted serpentinized slivers have been suggested in the southwest Pacific (Cluzel et al. 2001) and in Cascadia (Nikulin et al. 2009).

A third interpretation is that the slow velocities result from the accretion of an overthickened oceanic plateau that may have relatively slow seismic velocities in the upper layers (e.g. the Shatsky Rise; Korenaga & Sager 2012; Tetreault & Buiter 2014). For example, Sallares et al. (2005) found that oceanic layer 3 in the CR crust has $V_p$ velocities between 6.5 and 7.2 km s⁻¹ ($V_s$ velocities of
3.67–4.07 km s⁻¹ assuming *V*<sub>p</sub>/*V*<sub>s</sub> of 1.77) as it enters the trench, similar to the forearc seismic velocities seen beneath the Manabí Basin. Some of the forearc in this segment could be underplated oceanic plateau crust from previous plateau accretion. Given the tectonic setting of the Ecuadorian forearc, it seems plausible that a portion of the oceanic plateau that makes up the Piñon Formation was overthickened and the upper layers were accreted onto the forearc while deeper, high velocity, and denser layers may have delaminated, subducted, or have since been tectonically eroded (Sage et al. 2006).

The final possibility suggests that as the Nazca plate subducts beneath the forearc, the buoyancy of the CR inhibits subduction, increasing tectonic erosion of the forearc crust. Increased erosion may aid in the transportation of fluids, as the upper plate is physically broken up to accommodate the bathymetry of the ridge. As in case 1, the influx of fluids may lead to alteration of the forearc and therefore slower seismic velocities. The observed slow velocities may also be in part related to the CR itself. If the ridge is in direct contact with the overriding plate, the observed velocities may be reflective of ridge material not the upper plate material.

From shear wave velocity results alone, it is not possible to distinguish between the mechanisms described above. Therefore, we estimated bulk $V_p/V_s$ measurements of the slow region from RFs. $V_p/V_s$ can provide additional constraints on the character and composition of the slow region below the Manabí basin. The relatively high $V_p/V_s$ ratios (1.83–1.94) found at stations EC16 and EC08, which lie above the slow anomaly in S2, suggest that either fluids or serpentinization are likely responsible. An alteration-based mechanism is, therefore, the most plausible explanation of the slower seismic velocities.

## Southern segments 3 and 4

Segments 3 and 4 consist of the fastest velocities seen in the forearc. These velocities rapidly increase with depth, reaching up to 4.6 km s⁻¹ at ~25 km depth in S3 and as shallow as ~10 km in S4. Given the fast seismic velocities, it is likely that the accreted terrane in this region includes a large amount mafic and ultramafic material consistent with lithospheric mantle material. The Piñon Formation is exposed south of S4 and the San Lorenzo Formation outcrops along the coast on the west side of S3 and S4 (Luzieux et al. 2006, Jaillard et al. 2009, Fig. 1). We note that Reyes & Michaud (2012) map many of the San Lorenzo Formation outcrops along the coast as Pinon Formation. The higher velocities seen in S3 and S4 are consistent with those of mantle lithosphere at the base of accreted extinct arc terranes, such as Kohistan Arc and Talkeetna Arc, where P-wave velocities in the interpreted uppermost mantle reach P-wave velocities as high as 8.5 km s⁻¹ (e.g. Miller & Christensen 1994; Christensen & Mooney 1995; Jagoutz & Behn 2013; $V_s$ velocity of 4.8 assuming $V_p/V_s$ of 1.77). As such, these velocities may represent a portion of accreted lithospheric mantle, suggesting that the San Lorenzo Formation may extend at depth further eastward beneath the forearc (labelled as San Lorenzo Mantle in Fig. 6). Additionally, a strong positive Bouguer gravity anomaly in S4 has previously been interpreted as shallow mantle material (Feininger & Seguin 1983). This region of high Bouguer gravity anomaly extends from ~1°S to ~2.5°S and closely corresponds with the fastest velocities in S4 (Fig. S11; Bonvalot et al. 2012). Given the terranes' position in the forearc and depth of the subducting oceanic crust (Slab2 from Hayes et al. 2018), we suggest that the forearc accreted material likely includes some mantle lithosphere

associated with the accreted terrane and is in contact with the subducting oceanic crust. In the ANT-only model, Lynner et al. (2020) saw a similar fast velocity anomaly in the forearc and interpreted a similar mantle source to explain the velocities.

In S3 and S4, the San Lorenzo Formation is 15–25 km thick (Fig. 6), comparable to the preserved thicknesses of the Bonanza arc (Canil et al. 2010) and Talkeetna arc (Greene et al. 2006). Similar to the S1, it is difficult to resolve layering in the accreted material. However, in S3, a small section of relatively slow velocity is present at 10–20 km depth (labeled as San Lorenzo crust in Fig. 6). The velocity of this section is ~3.8–4.0 km s⁻¹ and may be a crustal portion of the accreted San Lorenzo Formation.

## Subducting slab

Below the crustal structures discussed above, relatively slow velocities (~3.6 km s⁻¹) are observed at ~28–32 km depth in S1, ~18 km depth in S2 and ~45 km depth in S3 and S4 in cross section A–A'. The deeper widespread slow velocities are approximately 10–15 km thick and closely overlays Slab2 in the north and underlays Slab2 in the south of A–A' and also in the south of B–B' (Fig. 4). The shear wave velocity of this anomaly is consistent with offshore seismic refraction velocity estimates of the CR oceanic crust (P-wave average velocities between 6.5 and 7.0 km s⁻¹ in oceanic layer 3, corresponding to $V_s$ = 3.51–3.78 km s⁻¹ assuming a higher $V_p/V_s$ of 1.85 which may be more consistent for oceanic crust; Hyndman et al. 1979; Sallarès et al. 2005; Gailler et al. 2007). Sallarès et al. (2005) also found the thickness of the CR oceanic crust to be between 13 and 19 km, matching closely with our observations. Along E–E', this slow velocity anomaly is seen dipping eastward (Fig. 4). However, in D–D' the slow velocities are not observed dipping east. Negative multiples, which may be present in S1 and S4, may impart a low velocity feature in these regions. However, given the geometry, close alignment with Slab2 in the S1, S3 and S4, and velocities of this slow seismic anomaly, we interpret this as the subducted Nazca oceanic crust (Fig. 6). The inconsistencies with Slab2 are not entirely unexpected for this region, as Kwong et al. (2019) found NEIC catalogue epicentres to be off by up to 15–25 km in this area.

South of the Jama Fault system (Fig. 1), towards S2 in A–A', the slab slow velocity anomaly deviates from Slab2 (Hayes et al. 2018), appearing to shallow to a depth of ~22 km. The shallowest point is spatially correlated with the centre of the CR when projected in the direction of plate motion (Figs 3 and 4). If this is indeed the position of the subducting slab, the forearc crust seems significantly thinned in the centre of the profile, possibly indicating increased tectonic erosion as a result of the CR subduction. Further evidence for this being a CR-related structure is its correlation with coastal uplift seen along the margin (Gutscher et al. 1999).

Offshore deformation and tectonic erosion have been suggested for the CR in seismic reflection profiles along the Ecuadorian margin (e.g. Sage et al. 2006). Other studies, however, have suggested that the coastal uplift seen in Ecuador is not necessarily related to the CR and that it is unclear how far inland the CR actually subducts (Michaud et al. 2009). Furthermore, there is a wide range of proposed ages of trench-ridge collision from 1 to 15 Ma (Michaud et al. 2009, and references therein). In order to explain the slab anomaly described in S2 the CR would have had to have subducted a distance of at least 130 km from the trench. This would imply that the CR started subducting into the trench by at least ~2.4 Ma (assuming a rate of 55 mm yr⁻¹, Chlieh et al. 2014). As such, it is unclear if the slab anomaly is present at these shallow depths
(segment 2 in A–A' Fig. 4) or if the shallower slow anomaly in S2 discussed earlier extends to these depths and that the slab is obscured in this section.

## Forearc structure and megathrust earthquakes

Our shear wave velocity model shows structural forearc segmentation at ~0.5°S and ~1°S coincident with the change in megathrust behaviour. This change in seismic velocity is seen throughout the forearc crust and potentially into the subducting oceanic crust. A change is also observed in the Bouguer gravity data and in the ANT-only model (Figs S9 and S11; Lynner et al. 2020). The margin north of ~0.5°S has hosted several $M_w > 7.7$ earthquakes in the past century, whereas south of ~0.5°S no megathrust earthquakes with $M_w \ge 7.4$ are present until at least ~3.5°S. We propose that the upper plate intermediate shear wave velocity material (Piñon Formation; Fig. 6) in S1 and the slow shear wave velocity anomaly beneath the Manabí Basin in S2 are able to host high magnitude earthquakes (such as the 1906, 1942 and 2016 events). In contrast, the southern portion of our study area (S3 and S4), where we interpret accreted mantle lithosphere in the upper plate near the down dip edge of the megathrust plate interface, may be less conducive to producing high magnitude earthquakes. In Ecuador, previous studies (e.g. Mendoza & Dewey 1984; Collot et al. 2004) have suggested that the offshore extension of the Jama Fault system and heterogeneities along fault planes are responsible for the segmentation and large earthquake behaviour. As shown in A–A', the slow velocities beneath the Manabí Basin shallow near the Jama Fault, further supporting the idea that the offshore extension of the fault system may segment large earthquake behaviour. Other major transitions in our results are not as clearly aligned with known fault systems. However, the rapid lateral changes in seismic velocities and in the RF discontinuity structure that we observe may indicate previously unobserved faults at depth.

The role of the lower plate in megathrust earthquakes, such as the importance of asperities, trench sediments, and hydration, has been well studied (e.g. Bilek & Lay 2018, and references therein). This is also true in Ecuador, where Collot et al. (2017) found that irregularities in the seafloor bathymetry are spatially correlated with a highly coupled asperity beneath La Plata Island. Agurto-Detzel et al. (2019) also found that the slip mode in the Ecuadorian margin may be controlled by oceanic relief of the incoming slab. Variations in the hydration state may also alter the frictional properties of the megathrust contact, resulting in regions that are either more or less suited to large magnitude earthquakes (Audet & Schwartz 2013). Sediment thickness has also been correlated with high magnitude earthquakes along the margin. Regions with thicker sediments may lead to greater volumes of subducted sediments, creating a smoother interface with more uniform coupling along the megathrust (Ruff 1989; Scholl et al. 2015). To first order, this seems to hold true for the Ecuadorian trench. Using wide-angle seismic profiles Gailler et al. (2007) found that sediment thickness north of the CR is ~1 km greater than to the south. However, the thinnest sediments are observed where the CR enters the trench, which transects the proposed segmentation (Gailler et al. 2007). All of these factors need to be considered to better understand megathrust properties.

The role of the upper plate in the occurrence of megathrust earthquakes in Ecuador is less well understood. Wells et al. (2003) suggested that subsidence associated with forearc basins may indicate areas of subduction erosion and higher coupling. Additionally, variations in crustal geology has been associated with variations in the frictional properties of the megathrust (e.g. in Japan, Fujie et al.

2013; Bassett et al. 2016). It is possible that a similar process is happening in Ecuador, in which we see significant along-strike variations in seismic velocities (Fig. 6) as well as in gravity (Feininger & Seguin 1983; Tamay et al. 2018). These seem to correlate with variable megathrust behaviour. Similar to the effect of a large sediment influx into the trench, a more homogenous upper crust may result in smaller variations in friction along strike, allowing larger areas to slip in single events (Bassett et al. 2016). These frictional variations may cause favorable conditions for large earthquakes in the north, where we interpret accreted Pinon Formation. At the same time, they cause aseismic creep in the south, where accreted lithospheric mantle may be present. The base of the accreted lithospheric mantle may have a thin layer of serpentinized material, which would behave in a more ductile manner, inhibiting stress build up and brittle failure associated with large earthquakes. Antigorite-rich serpentinite may influence the slip behaviour of a subduction interface, widening the range at which slow slip behaviour can occur (Goswami & Barbot 2018). However, subduction channel rocks, such as calcareous rocks, quartzose rocks, and talc schist have been shown to localize deformation relevant to the depths at which SSEs occur (French & Condit 2019). Thus, while the properties of the subducting plate are undeniably important, variations in the upper plate also seem to be important in controlling the slip behaviour of the megathrust along the Ecuador forearc.

## CONCLUSIONS

We present adaptive common conversion point receiver function stacks and a 3-D shear wave velocity model derived from the joint inversion of receiver functions and surface wave dispersion for the upper 50 km of the Ecuadorian forearc. Our results present new evidence for along-strike segmentation in the composition and/or properties of the forearc crust. We identify four distinct segments along-strike that may play important roles in controlling the areal extent of rupture for large events along the convergent margin. In the northernmost portion of the forearc (Segment 1), we observe fast velocities consistent with the accreted oceanic lithosphere of the Piñon Formation. Slightly to the south (Segment 2), we observe a relatively slow velocity anomaly that may be related to subduction derived fluids or the accretion of a slow velocity material beneath the base of the Manabí Basin and the top of the subducting lithosphere. Further south, at ~0.5°S, above the southern half of the subducted Carnegie Ridge (Segments 3 and 4), the forearc shows the highest velocities seen in our results, and are interpreted as being associated with lithospheric mantle material likely associated with the San Lorenzo Formation. The transition from slow velocities beneath the Manabí Basin to faster velocities southward corresponds to the southern extents of the 1906, 1942 and 2016 earthquakes and may act as a barrier to rupture propagation. Furthermore, our results suggest that structure and composition of the upper-plate in the Ecuadorian forearc may play a role in controlling the behaviour of megathrust earthquakes along-strike. The presence of the Carnegie Ridge and other structure on the subducting crust are undoubtedly a critical component as well. Further study of megathrust behaviour is needed to better understand the importance the overriding plate.

## ACKNOWLEDGEMENTS

We used the Generic Mapping Tools developed by Wessel & Smith (1998) to make figures. Open and available seismic data [networks 8 G (Meltzer & Beck 2016) and EC (Alvarado et al. 2018)] were
accessed via the Data Management Center of the Incorporated Research Institutions for Seismology. We also use data from the XE seismic network (Regnier et al. 2016) and the Instituto Geofisico at the Escuela Politécnica Nacional (IG-EPN) seismic network. We are deeply grateful to the staff at the Instituto Geofisico at the Escuela Politécnica Nacional (IG-EPN) in Quito, Ecuador and thank all the field crews that helped in data collection. We acknowledge the thoughtful comments from the editor and two anonymous reviewers. This work was supported by (National Science Foundation) NSF EAR-1723065, NSF EAR-1723042 and NSF RAPID Program Award EAR-1642498.

## REFERENCES

Agurto-Detzel, H. et al., 2019. Ridge subduction and afterslip control after-shock distribution of the 2016 Mw 7.8 Ecuador earthquake, *Earth planet. Sci. Lett.*, **520**, 63–76.

Alvarado, A. et al., 2018. Seismic, volcanic, and geodetic networks in Ecuador: building capacity for monitoring and research, *Seismol. Res. Lett.*, **89**, 432–439.

Ammon, C.-J., 1991. The isolation of receiver effects from teleseismic P waveforms, *Bulletin-Seismological Society of America*, **81**(8):2504–2510.

Audet, P. & Schwartz, S.Y., 2013. Hydrologic control of forearc strength and seismicity in the Costa Rican subduction zone, *Nat. Geosci.*, **6**, 852.

Bassett, D., Sandwell, D.T., Fialko, Y. & Watts, A.B., 2016. Upper-plate controls on co-seismic slip in the 2011 magnitude 9.0 Tohoku-oki earthquake, *Nature*, **531**, 92.

Beck, S.L. & Ruff, L.J., 1984. The rupture process of the great 1979 Colombia earthquake: evidence for the asperity model, *J. geophys. Res.*, **89**, 9281–9291.

Bilek, S.L. & Lay, T., 2018. Subduction zone megathrust earthquakes, *Geosphere*, **14**, 1468–1500.

Bonvalot, S., Balmino, G., Briais, A., Kuhn, M., Peyrefitte, A., Vales, N. & Sarrailh, M., 2012. World gravity map, Bureau Gravimetrique International (BGI), Map, CGMW-BGI-CNES728, IRD, Paris.

Bostock, M., Hyndman, R., Rondenay, S. & Peacock, S., 2002. An inverted continental Moho and serpentinization of the forearc mantle, *Nature*, **417**, 536.

Canil, D., Styan, J., Larocque, J., Bonnet, E. & Kyba, J., 2010. Thickness and composition of the Bonanza arc crustal section, Vancouver Island, Canada, *Bulletin*, **122**, 1094–1105.

Chlieh, M. et al., 2014. Distribution of discrete seismic asperities and aseismic slip along the Ecuadorian megathrust, *Earth planet. Sci. Lett.*, **400**, 292–301.

Christensen, N.I. & Mooney, W.D., 1995. Seismic velocity structure and composition of the continental crust: a global view, *J. geophys. Res.*, **100**(B6), 9761–9788.

Cloos, M., 1993. Lithospheric buoyancy and collisional orogenesis: subduction of oceanic plateaus, continental margins, island arcs, spreading ridges, and seamounts, *Bull. geol. Soc. Am.*, **105**, 715–737.

Cluzel, D., Aitchison, J.C. & Picard, C., 2001. Tectonic accretion and under-plating of mafic terranes in the Late Eocene intraoceanic fore-arc of New Caledonia (Southwest Pacific): geodynamic implications, *Tectonophysics*, **340**, 23–59.

Collot, J.Y. et al., 2004. Are rupture zone limits of great subduction earthquakes controlled by upper plate structures? Evidence from multichannel seismic reflection data acquired across the northern Ecuador–southwest Colombia margin, *J. geophys. Res.*, **109**, doi:10.1029/2004JB003060.

Collot, J.Y. et al., 2017. Subducted oceanic relief locks the shallow megathrust in central Ecuador, *J. geophys. Res.*, **122**, 3286–3305.

Delph, J.R., Zandt, G. & Beck, S.L., 2015. A new approach to obtaining a 3D shear wave velocity model of the crust and upper mantle: An application to eastern Turkey, *Tectonophysics*, **665**, 92–100.

Delph, J.R., Levander, A. & Niu, F., 2018. Fluid controls on the heterogeneous seismic characteristics of the Cascadia margin, *Geophys. Res. Lett.*, **45**, 11 021–11 029.

Delph, J.R., Ward, K.M., Zandt, G., Ducea, M.N. & Beck, S.L., 2017. Imaging a magma plumbing system from MASH zone to magma reservoir, *Earth planet. Sci. Lett.*, **457**, 313–324.

Deniaud, Y., 2000. Enregistrements sédimentaire et structural de l' évolution géodynamique des Andes équatoriennes au cours du Néogène: étude des bassins d'avant-arc et bilans de masse.

Dueker, K.G. & Sheehan, A.F., 1997. Mantle discontinuity structure from midpoint stacks of converted P to S waves across the Yellowstone hotspot track, *J. geophys. Res.*, **102**, 8313–8327.

Eagar, K.C. & Fouch, M.J., 2012. FuncLab: A MATLAB interactive toolbox for handling receiver function datasets, *Seismol. Res. Lett.*, **83**, 596–603.

Evans, C. & Whittaker, J., 1982. The geology of the western part of the Borbón Basin, North-west Ecuador, *Geol. Soc. Lond., Spec. Publ.*, **10**, 191–198.

Feininger, T. & Seguin, M., 1983. Simple Bouguer gravity anomaly field and the inferred crustal structure of continental Ecuador, *Geology*, **11**, 40–44.

Font, Y., Segovia, M., Vaca, S. & Theunissen, T., 2013. Seismicity 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.

French, M.E. & Condit, C.B., 2019. Slip partitioning along an idealized subduction plate boundary at deep slow slip conditions, *Earth planet. Sci. Lett.*, **528**, 115828.

Fujie, G. et al., 2013. Along-trench structural variation and seismic coupling in the northern Japan subduction zone, *Earth Planets Space*, **65**, 75–83.

Gailler, A., Charvis, P. & Flueh, E.R., 2007. Segmentation of the Nazca and South American plates along the Ecuador subduction zone from wide angle seismic profiles, *Earth planet. Sci. Lett.*, **260**, 444–464.

Goswami, A. & Barbot, S., 2018. Slow-slip events in semi-brittle serpentinite fault zones, *Sci. Rep.*, **8**(1), 1–11.

Graindorge, D., Calahorrano, A., Charvis, P., Collot, J.Y. & Bethoux, N., 2004. Deep structures of the Ecuador convergent margin and the Carnegie Ridge, possible consequence on great earthquakes recurrence interval, *Geophys. Res. Lett.*, **31**, doi:10.1029/2003GL018803.

Greene, A., DeBari, S., Kelemen, P., Blusztajn, J. & Clift, P. 2006. A Detailed Geochemical Study of Island Arc Crust: the Talkeetna Arc Section, South–Central Alaska, *Journal of Petrology*, **47**(8):1051–1093.

Guillot, S., Schwartz, S., Reynard, B., Agard, P. & Prigent, C., 2015. Tectonic significance of serpentinites, *Tectonophysics*, **646**, 1–19.

Gutscher, M.-A., Malavieille, J., Lallemand, S. & Collot, J.-Y., 1999. Tectonic segmentation of the North Andean margin: impact of the Carnegie Ridge collision, *Earth planet. Sci. Lett.*, **168**, 255–270.

Hayes, G.P., Moore, G.L., Portner, D.E., Hearne, M., Flamme, H., Furtney, M. & Smoczyk, G.M., 2018. Slab2, a comprehensive subduction zone geometry model, *Science*, **362**, 58–61.

Hughes, R.A. & Pilatasig, L.F., 2002. Cretaceous and Tertiary terrane accretion in the Cordillera Occidental of the Andes of Ecuador, *Tectonophysics*, **345**, 29–48.

Hyndman, R., Davis, E. & Wright, J., 1979. The measurement of marine geothermal heat flow by a multipenetrating probe with digital acoustic telemetry and in situ thermal conductivity, *Mar. Geophys. Res.*, **4**, 181–205.

Hyndman, R.D. & Peacock, S.M., 2003. Serpentinization of the forearc mantle, *Earth planet. Sci. Lett.*, **212**, 417–432.

International Seismological Centre, 2016. On-line Bulletin, http://www.isc.ac.uk, Internatl. Seismol. Cent., Thatcham, United Kingdom.

Jagoutz, O. & Behn, M.D., 2013. Foundering of lower island-arc crust as an explanation for the origin of the continental Moho, *Nature*, **504**, 131.

Jaillard, E., Lapierre, H., Ordonez, M., Alava, J.T., Amortegui, A. & Vannelle, J., 2009. Accreted oceanic terranes in Ecuador: southern edge of the Caribbean Plate? *Geol. Soc. Lond., Spec. Publ.*, **328**, 469–485.

Jaillard, E., Ordoñez, M., Benitez, S., Berrones, G., Jiménez, N., Montengro, G. & Zambrano, I., 1995. Basin development in an accretionary, oceanic-floored fore-arc setting: southern coastal Ecuador during late Cretaceous-late Eocene time.

Julia, J., Ammon, C.J., Herrmann, R. & Correig, A.M., 2000. Joint inversion of receiver function and surface wave dispersion observations, *Geophys. J. Int.*, **143**, 99–112.
Kanamori, H. & McNally, K.C., 1982. Variable rupture mode of the subduction zone along the Ecuador-Colombia coast, *Bull. seism. Soc. Am.*, **72**, 1241–1253.

Kerr, A.C., Aspden, J.A., Tarney, J. & Pilatasig, L.F., 2002. The nature and provenance of accreted oceanic terranes in western Ecuador: geochemical and tectonic constraints, *J. Geol. Soc.*, **159**, 577–594.

Korenaga, J. & Sager, W., 2012. Seismic tomography of Shatsky Rise by adaptive importance sampling, *J. geophys. Res.*, **117**, doi:10.1029/2012JB009248.

Kwong, K.B., DeShon, H.R., Kim, J.W. & Lu, Z., 2019. Resolving teleseismic earthquake catalog and InSAR data discrepancies in absolute space to explore rupture complexity along the Ecuadorian Megathrust Fault, *J. geophys. Res.*, **124**(7), 6703–6719.

Langston, C., 1979. Structure Under Mount Rainier, Washington, Inferred From Teleseismic Body Waves, *Journal of Geophysical Research Solid Earth*, **84**(8):4749–4762.

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

Ligorria, J.P. & Ammon, C.J., 1999. Iterative deconvolution and receiver-function estimation, *Bull. seism. Soc. Am.*, **89**, 1395–1400.

Luzieux, L., 2007. *Origin and Late Cretaceous-Tertiary Evolution of the Ecuadorian Forearc*, ETH Zurich.

Luzieux, L., Heller, F., Spikings, R., Cruz, C.V. & Winkler, W., 2006. Origin and Cretaceous tectonic history of the coastal Ecuadorian forearc between 1 degrees N and 3 degrees S: Paleomagnetic, radiometric and fossil evidence, *Earth planet. Sci. Lett.*, **249**, 400–414.

Lynner, C. et al., 2020. Upper-plate structure in Ecuador coincident with the subduction of the Carnegie Ridge and the southern extent of large mega-thrust earthquakes, *Geophys. J. Int.*, **220**(3), 1965–1977.

Ma, Y. & Clayton, R.W., 2016. Structure of the Los Angeles basin from ambient noise and receiver functions, *Geophys. J. Int.*, **206**, 1645–1651.

Meltzer, A. & Beck, S.L., 2016. *Pedernales Earthquake Aftershock Deployment Ecuador*, International Federation of Digital Seismograph Networks. Other/Seismic Network. doi:10.7914/SN/8G.2016.

Meltzer, A. et al., 2019. The 2016 M w 7.8 Pedernales, Ecuador, earthquake: rapid response deployment, *Seismol. Res. Lett.*, **90**, 1346–1354.

Mendoza, C. & Dewey, J.W., 1984. Seismicity associated with the great Colombia-Ecuador earthquakes of 1942, 1958, and 1979: Implications for barrier models of earthquake rupture, *Bull. seism. Soc. Am.*, **74**, 577–593.

Michaud, F., Witt, C. & Royer, J.-Y., 2009. Influence of the subduction of the Carnegie volcanic ridge on Ecuadorian geology: Reality and fiction, *Backbone of the Americas: Shallow Subduction, Plateau Uplift, and Ridge and Terrane Collision*, **204**, 217.

Miller, D.J. & Christensen, N.L., 1994. Seismic signature and geochemistry of an island arc: a multidisciplinary study of the Kohistan accreted terrane, northern Pakistan, *J. geophys. Res.*, **99**, 11623–11642.

Nikulin, A., Levin, V. & Park, J., 2009. Receiver function study of the Cascadia megathrust: evidence for localized serpentinization, *Geochem. Geophys. Geosyst.*, **10**,

Niu, F., Levander, A., Ham, S. & Obayashi, M., 2005. Mapping the subducting Pacific slab beneath southwest Japan with Hi-net receiver functions, *Earth planet. Sci. Lett.*, **239**, 9–17.

Nocquet, J.-M. et al., 2014. Motion of continental slivers and creeping subduction in the northern Andes, *Nat. Geosci.*, **7**, 287.

Nocquet, J.-M. et al., 2017. Supercycle at the Ecuadorian subduction zone revealed after the 2016 Pedernales earthquake, *Nat. Geosci.*, **10**, 145.

Peacock, S.M., 1989. Numerical constraints on rates of metamorphism, fluid production, and fluid flux during regional metamorphism, *Bull. geol. Soc. Am.*, **101**, 476–485.

Porritt, R.W. & Miller, M.S., 2018, Updates to FuncLab, a Matlab based GUI for handling receiver functions, *Comp. Geosci.*, **111**, 260–271.

Rathnayaka, S. & Gao, H., 2017. Crustal-scale seismic structure from trench to forearc in the Cascadia subduction zone, *J. geophys. Res.*, **122**, 7311–7328.

Regnier, M., Font, Y., Charvis, P., Mercerat, D., Rietbrock, A., Ruiz, M. & Alvarado, A., 2016. Pedernales. International Federation of Digital Seismograph Networks, Dataset/Seismic Network, doi:10.7914/SN/XE.2016.

Ridley, V.A. & Richards, M.A., 2010. Deep crustal structure beneath large igneous provinces and the petrologic evolution of flood basalts, *Geochem., Geophys., Geosyst.*, **11**, doi:10.1029/2009GC002935

Rolandone, F. et al., 2018. Areas prone to slow slip events impede earthquake rupture propagation and promote afterslip, *Sci. Adv.*, **4**, eaao6596.

Ruff, L.J., 1989. Do trench sediments affect great earthquake occurrence in subduction zones?in *Subduction Zones Part II*, pp. 263–282, Springer.

Sage, F., Collot, J.-Y. & Ranero, C.R., 2006. Interplate patchiness and subduction-erosion mechanisms: Evidence from depth-migrated seismic images at the central Ecuador convergent margin, *Geology*, **34**, 997–1000.

Sallarès, V., Charvis, P., Flueh, E.R., Bialas, J. & Party, S.S., 2005. Seismic structure of the Carnegie ridge and the nature of the Galapagos hotspot, *Geophys. J. Int.*, **161**, 763–788.

Scholl, D.W., Kirby, S.H., von Huene, R., Ryan, H., Wells, R.E. & Geist, E.L., 2015. Great (≥ Mw8.0) megathrust earthquakes and the subduction of excess sediment and bathymetrically smooth seafloor, *Geosphere*, **11**, 236–265.

Segovia, M., Font, Y., Régnier, M., Charvis, P, Glave, Audrey, Nocquet, J.-M., Jarrin, P., Hello, Y., Ruiz, M. & Pazmino, A., 2018. Seismicity Distribution Near a Subducting Seamount in the Central Ecuadorian Subduction Zone, Space-Time Relation to a Slow-Slip Event, *Tectonics*, **37**(8):2106–2123.

Soto-Cordero, L. et al., 2020. Structural control on megathrust rupture and slip behavior: insights from the 2016 Mw 7.8 Pedernales Ecuador Earthquake, *J. geophys. Res.*, 125, doi:10.1029/2019JB018001, e2019JB018001.

Swenson, J.L. & Beck, S.L., 1999. Source characteristics of the 12 November 1996 M w 7.7 Peru subduction zone earthquake, in *Seismogenic and Tsunamigenic Processes in Shallow Subduction Zones*, pp. 731–751, Springer.

Tamay, J., Galindo-Zaldivar, J., Martos, Y. & Soto, J., 2018. Gravity and magnetic anomalies of ecuadorian margin: Implications in the deep structure of the subduction of Nazca Plate and Andes Cordillera, *J. South Am. Earth Sci.*, **85**, 68–80.

Tetreault, J. & Buiter, S., 2014. Future accreted terranes: a compilation of island arcs, oceanic plateaus, submarine ridges, seamounts, and continental fragments, *Solid Earth*, **5**, 1243–1275.

Vaca, S., Vallée, M., Nocquet, J.-M., Battaglia, J. & Régnier, M., 2018. Recurrent slow slip events as a barrier to the northward rupture propagation of the 2016 Pedernales earthquake (Central Ecuador), *Tectonophysics*, **724**, 80–92.

Vallejo, C., Winkler, W., Spikings, R.A., Luzieux, L., Heller, F. & Bussy, F., 2009. Mode and timing of terrane accretion in the forearc of the Andes in Ecuador, *Backbone of the Americas: Shallow Subduction, Plateau Uplift, and Ridge and Terrane Collision*, **204**, 197.

Vallée, M. et al., 2013. Intense interface seismicity triggered by a shallow slow slip event in the Central Ecuador subduction zone, *J. geophys. Res.*, **118**, 2965–2981.

Wells, R.E., Blakely, R.J., Sugiyama, Y., Scholl, D.W. & Dinterman, P.A., 2003. Basin-centered asperities in great subduction zone earthquakes: a link between slip, subsidence, and subduction erosion? *J. geophys. Res.*, 108.

Wessel, P. & Smith, W.H., 1998. New, improved version of Generic Mapping Tools released, *EOS, Trans. Am. Geophys. Un.*, **79**, 579–579.

Yamanaka, Y., Tanioka, Y. & Shiina, T., 2017. A long source area of the 1906 Colombia–Ecuador earthquake estimated from observed tsunami waveforms, *Earth, Planets Space*, **69**, 163.

Ye, L., Kanamori, H., Avouac, J.-P., Li, L., Cheung, K.F. & Lay, T., 2016. The 16 April 2016, MW 7.8 (MS 7.5) Ecuador earthquake: a quasi-repeat of the 1942 MS 7.5 earthquake and partial re-rupture of the 1906 MS 8.6 Colombia–Ecuador earthquake, *Earth planet. Sci. Lett.*, **454**, 248–258.

Zandt, G., Myers, S. & Terry, W., 1995. Crust and mantle structure across the Basin and Range- Colorado Plateau boundary at 370N latitude and implications for Cenozoic extensional mechanism, *Journal of Geophysical Research*, **100**:10529–10548.
SUPPORTING INFORMATION

Supplementary data are available at GJI online.

**Figure S1.** 1-D velocity model constructed by averaging a subset of the 3-D *a priori* P-wave velocity model from Font *et al.* (2013). The velocity model was averaged over a rectangular region in the Manabi basin area, defined by a box with latitude 0.04°N to 0.4°S and longitude 80.02°W to 79.8°W. This velocity model is utilized in both the ACCP stacks and CCP stacks.

**Figure S2.** Results of the joint inversion for individual gridpoints at locations nearest the stations shown in Fig. 2. (c) Predicted (red) and observed (black) receiver function, (b) predicted phase velocity (red) and observed (black) dispersion data and (c) initial (black) and final (red) shear wave velocity model.

**Figure S3.** Plot showing the average and standard deviation of the NRMS fit per cent for dispersion (purple) and receiver function (orange) data sets vs joint inversion weighting factor.

**Figure S4.** Cross sections A–A' and B–B' through different shear wave velocity models resulting from the joint inversion of receiver functions and dispersion data. Black line in each section is from the Slab2 model (Hayes *et al.* 2018). Unless otherwise specified, all other parameters are as described in the main text. (a) Cross sections using a uniform 4.8 km s⁻¹ starting velocity model. (b) Cross sections using a uniform 4.2 km s⁻¹ starting velocity model. (c) Cross sections using a weighting of 0.7 (70 per cent dispersion data, 30 per cent RF data). (d) Cross sections through the models with dispersion periods >25 s removed.

**Figure S5.** Cross sections C–C', D–D' and E–E' through different shear wave velocity models resulting from the joint inversion of receiver functions and dispersion data. Black line in each section is from the Slab2 model (Hayes *et al.* 2018). Unless otherwise specified all other parameters are as described in the main text. (a) Cross sections using a uniform 4.8 km s⁻¹ starting velocity model. (b) Cross sections using a uniform 4.2 km s⁻¹ starting velocity model. (c) Cross sections using a weighting of 0.7 (70 per cent dispersion data, 30 per cent RF data). (d) Cross sections through the models with dispersion periods >25 s removed.

**Figure S6.** Effective 2.5 Gaussian CCP stacks (data set used in the joint inversion) with the shallowest primary *P*-to-*S* conversion (black) and corresponding location of predicted multiples assuming

a $V_p/V_s$ of 1.77 (red and blue). The pink line is the top of the oceanic crust from Slab2 (Hayes *et al.* 2018).

**Figure S7.** Trench parallel cross sections A–A' (top panel) and B–B' (bottom panel) through the joint inversion shear velocity model and the 5.0 Gaussian ACCP stack. Pink line is the top of the subducting oceanic crust from Slab2 (Hayes *et al.* 2018). The velocity model used for ACCP migration and other parameters are described in the main text.

**Figure S8.** Trench perpendicular cross sections through the 5.0 Gaussian ACCP stack (left-hand panel) and joint inversion shear velocity model (right-hand panel). Black dots are earthquakes from the Pedernales aftershock catalog projected from 10 km perpendicular to the section (Agurto-Detzel *et al.* 2019). Red line shows the extent of the 2016 earthquake rupture (>1 m slip; Nocquet *et al.* 2017). The pink line is the top of the subducting oceanic crust from Slab2 (Hayes *et al.* 2018).

**Figure S9.** Cross sections A–A' and B–B' through shear wave velocity models derived from the joint inversion of receiver functions and dispersion data (left-hand panel) and the dispersion data only (right-hand panel; Lynner *et al.* 2020). Lynner *et al.* (2020) used the 3-D *a priori* P-wave velocity model from Font *et al.* (2013) as their starting model. Fit per cent for both RF (red) and surface wave dispersion data (blue) from the joint inversion are shown above the joint inversion profiles.

**Figure S10.** Cross sections C–C', D–D' and E–E' through shear wave velocity models derived from the joint inversion of receiver functions and dispersion data (left-hand panel) and the dispersion data only (right-hand panel; Lynner *et al.* 2020). Lynner *et al.* (2020) used the 3-D *a priori* P-wave velocity model from Font *et al.* (2013) as their starting model. Fit per cent for both RF (red) and surface wave dispersion data (blue) from the joint inversion are shown above the joint inversion profiles.

**Figure S11.** Bouguer gravity profiles through A–A' and B–B' cross sections along with the corresponding shear velocity profiles from this study (Bonvalot *et al.* 2012). Note the increase in the Bouguer gravity associated with the high shear wave velocity material in S3 and S4.

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.