Seismic and Aseismic Cycle of the Ecuador–Colombia Subduction Zone

The Colombia–Ecuador subduction zone is an exceptional natural laboratory to study the seismic cycle associated with large and great subduction earthquakes. Since the great 1906 Mw = 8.6 Colombia–Ecuador earthquake, four large Mw > 7.5 megathrust earthquakes occurred within the 1906 rupture area, releasing altogether a cumulative seismic moment of ∼35% of the 1906 seismic moment. We take advantage of newly released seismic catalogs and global positioning system (GPS) data at the scale of the Colombia–Ecuador subduction zone to balance the moment deficit that is building up on the megathrust interface during the interseismic period with the seismic and aseismic moments released by transient slip episodes. Assuming a steady-state interseismic loading, we found that the seismic moment released by the 2016 Mw = 7.8 Pedernales earthquake is about half of the moment deficit buildup since 1942, suggesting that the Pedernales segment was mature to host that seismic event and its postseismic afterslip. In the aftermath of the 2016 event, the asperities that broke in 1958 and 1979 both appears to be mature to host a large Mw > 7.5 earthquakes if they break in two individual seismic events, or an Mw∼7.8–8.0 earthquake if they break simultaneously. The analysis of our interseismic-coupling map suggests that the great 1906 Colombia–Ecuador earthquake could have ruptured a segment of 400 km-long bounded by two 80 km wide creeping segments that coincide with the entrance into the subduction of the Carnegie ridge in Ecuador and the Yaquina Graben in Colombia. These creeping segments share similar frictional properties and may both behave as strong seismic barriers able to stop ruptures associated with great events like in 1906. Smaller creeping segments are imaged within the 1906 rupture area and are located at the extremities of the large 1942, 1958, 1979, and 2016 seismic ruptures. Finally, assuming that the frequency–magnitude distribution of megathrust seismicity follows the Gutenberg–Richter law and considering that 50% of the transient slip on the megathrust is aseismic, we found that the maximum magnitude subduction earthquake that can affect this subduction zone has a moment magnitude equivalent to Mw ∼8.8 with a recurrence time of 1,400 years. No similar magnitude event has yet been observed in that region.


INTRODUCTION
The Colombia-Ecuador subduction zone extends about 900 km, measured along the trench axis, from the Gulf of Guayaquil (3°S) in South Ecuador to Buenavantura (4°N) in Colombia (Figure 1). This subduction segment has experienced the great (Mw > 8.5) 1906 Colombia-Ecuador megathrust earthquake; however, many uncertainties remain on its exact moment magnitude that fluctuates from 8.5 to 9.0 in the scientific literature and on its rupture length that is proposed between 300 and 500 km long (Kanamori and McNally, 1982;Kelleher, 1972;Mendoza and Dewey, 1984). From the analysis of long-period surface wave records, it was suggested that the 1906 earthquake did not exceed Mw 8.5 (Okal, 1992), a magnitude similar from its associated tsunami's height modeling with a rupture length of ∼300 km (Tsuzuki et al., 2017;Yoshimoto et al., 2017). Reassessments of the 1906 seismic moment suggest a moment magnitude Mw 8.6 (Ye et al., 2016) that will be our reference moment magnitude for this study (Table 1). Recent cores collected offshore Esmeraldas at 3,000 m below sea level that samples the last 800 years indicate that the comparison of the main features of the 1906 turbidite with older earthquake-triggered turbidites suggests that a great earthquake similar to the 1906 event might have occurred ∼600 years ago (Migeon et al., 2017).
During the 20th century, the 1906 rupture segment hosted three juxtaposed large (Mw > 7.5) seismic ruptures of about 100-150 km long (Figure 1 and Supplementary Figure S1). This sequence started off the Ecuadorian coast in 1942 with the Mw∼7.8 event (Mendoza and Dewey, 1984;Swenson and Beck, 1996), followed northward by the 1958 Mw 7.6 event and finished off the coast of Colombia in 1979 with the Mw 8.1 event (Beck and Ruff, 1984;Kanamori and McNally, 1982) ( Table 1). The same segment that ruptured in 1942 failed again during the 2016 Mw 7.8 Pedernales earthquake which initiated ∼50 km north of the equatorial line and propagated 120 km southward toward the Pedernales city which was severely impacted ( Figure 2). The Harvard focal mechanism (from the Global Centroid Moment Tensor catalog, gCMT, https://www. globalcmt.org/CMTsearch.html) associated with the 2016 Mw 7.8 Pedernales earthquake is consistent with a subduction earthquake with an average rupture strike 27°, dip 21°, slip rake 124°, and with a moment magnitude Mw 7.8 ( Figure 2). The occurrence of the 2016 event naturally raises the question about the beginning of a new seismic sequence in that region.
Geodetic measurements collected in the last decade in Ecuador have shown that the friction in the seismogenic zone along the plates interface is highly heterogeneous and confined between the trench axis and the first 40 km depth Nocquet et al., 2014). This heterogeneity reflects variations in the steadystate frictional properties of the interplate contact with velocityweakening behavior in areas that appear highly locked and velocity-strengthening behavior in areas that are creeping. Creeping is the most prominent process at the scale of the whole slab interface and the locking occurs only in very specific places of the seismogenic zone. It has been shown in many subduction zones and continental faults that there is always a high spatial correlation between the location of large earthquake sources and the location of locked asperities during the interseismic period (Chlieh et al., 2008;Chlieh et al., 2011;Konca et al., 2008;Loveless and Meade, 2011;Metois et al., 2016;Moreno et al., 2010). Here, we take advantage of newly released global positioning system (GPS) data acquired at the scale of the Colombia-Ecuador subduction zone to map the interseismic locking over the great 1906 rupture and neighboring subduction segments (Figures 1, 2). This map will provide us a better perception of the seismogenic zone FIGURE 1 | Seismotectonic setting of the Ecuador-Colombia subduction zone where the Nazca plate subducts below the North Andean Sliver (NAS) and the South America Craton (SOAM). The Carnegie Ridge and the Yaquina Graben are the two main geomorphological features on the oceanic seafloor. The Yaquina graben is an abandoned transform fault that behaves as a fracture zone that contributes to the reactivation of normal faults associated with paleo-rifts as the Malpelo rift and the Buenaventura rift that intersects the trench at 3.5°N (Lonsdale, 2005;MacMillan et al., 2004;Yepes et al., 2016) (see also Supplementary Figure S1). At the Colombia-Ecuador border, the relative Nazca/SOAM plate convergence rate is about 55 mm/yr (Kendrick et al., 2003). This oblique convergence is partitioned between quasi-normal motion accommodated on the subduction megathrust interface at the trench axis and right-lateral transpressional strike-slip motion on the Chingual-Cosanga-Pallatanga-Puná (CCPP) fault system Yepes et al., 2016). The ellipses indicate the rupture and aftershocks areas with associated epicenters (stars) of the great 1906 Mw 8.6 (black dashed lines), of the 1942 Mw 7.8 (light blue), of the 1958 Mw 7.6 (light green), and of the 1979 Mw 8.1 (green) megathrust earthquakes (Beck and Ruff, 1984;Kanamori and McNally, 1982;Mendoza and Dewey, 1984;Swenson and Beck, 1996). segmentation which is fundamental to understand the physical parameters that control the genesis, extent, and arrest of large megathrust earthquakes.
Aseismic slip has an important role in the seismic cycle and it is only with the advent of modern space-based geodesy that we have been able to quantify its associated aseismic moment. Aseismic slip appears in the form of afterslip during the postseismic period and of spontaneous slow slip events (SSEs) during the interseismic and pre-seismic periods. It is often accompanied by a specific seismicity signature: postseismic aftershocks due to the afterslip and tremors and very low frequency or repeaters seismicity during SSEs episodes. Systematically, the cumulative seismic moment associated with this seismic activity is found to be much lower than the afterslip and SSEs moments derived from geodesy, suggesting that these transient slip processes are dominated by aseismic slip. Because permanent geodetic observations are slowly developing in many countries during the last 2 decades, our actual knowledge of the aseismic slip process is limited to sparse reports documented in various subduction zones. When the spatial resolution of the data is high, we often found that the postseismic afterslip moment is about 10-30% of the coseismic moment (Bürgmann et al., 2000;Chlieh et al., 2007;Hsu et al., 2006;Lin et al., 2013;Perfettini et al., 2010;Remy et al., 2016;Shrivastava et al., 2016). For SSE, their recurrence is not the same from one region to another and their moment magnitude varies from one SSE to another in the same area (Ide et al., 2007;Michel et al., 2019;Radiguet et al., 2012;Vallée et al., 2013), making it difficult to have a precise long-term quantification of the aseismic moment associated with SSEs. This estimation has a direct impact on the determination of the moment magnitude and recurrence time of large seismic events, which makes it a fundamental parameter that is unfortunately often neglected in many seismic hazard studies.
In active tectonics, the seismic moment conservation principle, which arises from the elastic rebound theory, requires that the cumulative interseismic moment buildup on a fault plane balances the sum of all the transient seismic and aseismic moments released by pre-seismic slip, coseismic slip (earthquake), postseismic afterslip, and spontaneous slow slip event. Deterministic approach to estimate the recurrence time of large earthquakes relies on the use of the historical and instrumental seismic records. This assumption is reasonable in places where the earthquake catalog is long enough to represent the average seismicity over the long term which is unfortunately rarely the case. Indeed, the largest events are often not well represented as their recurrence time exceeds the typical duration of earthquake catalogs. To overcome such limitation, statistical methods are applied. When the magnitude-frequency distribution of earthquakes obeys the Gutenberg-Richter (G-R) law over the moderate magnitude range, rates of earthquakes in the upper magnitude range can be obtained by extrapolation (Aki, 1965;Gutenberg and Richter, 1954). The G-R distribution must be bounded, otherwise the moment released by earthquakes would exceed the moment deficit that builds up due to locking during the interseismic period, which is unphysical. The earthquake recurrence model should, therefore, be constrained to satisfy the principle of seismic moment conservation (Avouac, 2015). We adopt this approach here to balance the moment buildup during the interseismic period with the moment released through earthquakes and aseismic transient episodes along the Ecuador-Colombia subduction zone.  Figure S2) that complement the Ecuadorian velocity field published by Nocquet et al., 2014. This unified GPS  Newton-meter (N-m) and corresponding moment magnitude Mw (Kanamori, 1977) for transient seismic and aseismic events observed along the Ecuador-Colombia subduction zone. Note that the cumulative seismic moment of the 1942, 1958, 1979, and 2016  Frontiers in Earth Science | www.frontiersin.org September 2021 | Volume 9 | Article 701720 field covers the whole Ecuadorian subduction zone and the rupture area of the great 1906 earthquake that extended along the Colombian subduction zone ( Figure 2). Based on that geodetic data set, a reestimation of the location and the angular velocity of the North Andean Sliver (NAS) Euler pole was proposed to be located at 58.6°N, 174.8°W and rotates counter-clockwise at a rate of 0.072°/Ma relatively to a fixed South American (SOAM) reference frame (Mora-Páez et al., 2019). This pole predicts a NAS rigid motion of about 8.6 mm/yr toward N 60°E, a result relatively similar with the 8.5 ± 1.0 mm/yr proposed by Chlieh et al. (2014). The resulting Nazca/NAS convergence rate shows that along the southern Ecuador Trench in the Gulf of the Guayaquil area, the Nazca plate subducts beneath the North Andean Sliver at a rate of about 50 mm/yr in a direction nearly orthogonal to the trench. The direction of convergence becomes increasingly oblique toward the northwest and the relative plate convergence is partitioned into thrusting at the trench and right-lateral transpressional strike-slip motion on the Chingual-Cosanga-Pallatanga-Puná (CCPP) fault system Yepes et al., 2016). Offshore the Colombia-Ecuador border, the Nazca/NAS convergence rate is about 48 mm/yr and decreases northward up to 46 mm/yr offshore Buenaventura in Colombia ( Figure 2). To explore models of interseismic coupling on the subduction megathrust interface, we selected 43 horizontal GPS velocities of Mora-Paez et al. (2019) with their respective 1-sigma error ellipses collected in Ecuador and Colombia between latitudes 3°S and 4°N (Supplementary Table S1 and Figure 2 and Supplementary Figure S2). In the NAS reference frame, GPS velocities highlight significant variations in the elastic strain accumulation along the Nazca subduction interface. Near the Colombia-Ecuador border, coastal sites Esmeraldas and Tumaco moved east-northeastward at 22.2 ± 0.5 mm/yr and 18.3 ± 0.4 mm/yr, respectively. Velocities decrease northward with the Cali and Buenaventura velocities showing 7.6 ± 1.3 mm/yr and 9.7 ± 0.3 mm/yr, respectively. The spatial resolution of the slip distribution on the megathrust interface depends on the spatial density of observations at the surface. To evaluate the power of slip resolution on the megathrust interface of the Ecuadorian and Colombian GPS networks, we run several checkerboards resolution tests for characteristic rectangular elements of 120 × 60 km 2 (Supplementary Figure S3A) or square elements of 80 and 60 km 2 , respectively (Supplementary Figure S3B). Each of the elements are either fully locked or fully creeping at the plate convergence rate. The predicted velocities of the initial checkerboard models are computed at all the GPS sites and then inverted to evaluate the power of restitution of the GPS networks. We found that the lateral variations of the slip patches are relatively well reproduced for all characteristic size patches between slab depths of 20 and 100 km. Patches as small as 60 km 2 are also well reproduced for the shallower 20 km depths of the slab interface offshore Ecuador where the trench-coast distances are the shortest, that is, in Isla La Plata and in Punta Galera, and Manta and Punta Salinas Peninsulas (see C and D in Supplementary Figure S3B). Offshore Colombia, where the trench-coast distances are the highest, even large 120 × 60 km 2 patches are not resolved for the shallower 20 km depths of the megathrust interface located between 2°N and 4°N (Supplementary Figure S3A). The sparse distribution of coastal GPS sites in Colombia (only three sites, Tumaco, Guapi, and Buenaventura) together with five continental GPS sites allow, however, to restitute relatively well all patches of 80 km 2 located between 20 and 100 km depths of the slab interface as well as patches of 60 km 2 when GPS sites are lying right on top or nearby (Supplementary Figure S3A). The moment restitution between the initial and inverted checkerboard models is always higher than 97% suggesting that even if the spatial location of locked asperities can be locally lost near the trench, the Ecuador-Colombia GPS networks will remain quite sensitive to the degree of locking during the interseismic period.

INTERSEISMIC MODELING
Assuming that the surface deformation measured by geodesy is the result of slip on the megathrust interface, we perform FIGURE 2 | GPS velocities (black arrows) in the North Andean Sliver reference frame from Mora-Paez et al. (2019) with their respective error ellipses. The relative Nazca/NAS convergence rate taken up by the subduction interface is about 48 mm/yr at the Colombia-Ecuador border. Interseismic GPS velocities pointing towards the east show variable decreasing gradients which attest for variations on the frictional coupling along the megathrust interface. The epicenter location (red star) and focal mechanism (from the Harvard gCMT catalog) associated with the 2016 Mw 7.8 Pedernales earthquake are consistent with a subduction earthquake that occurred on the megathrust interface. The megathrust geometry is reported from Slab 2 model geometry (Hayes et al., 2018) with iso-depth contours every 20 km.
Frontiers in Earth Science | www.frontiersin.org September 2021 | Volume 9 | Article 701720 nonlinear inversions of the GPS data based on a stochasticsimulated annealing algorithm to determine coupling/slip on the megathrust, as previously done in Chlieh et al. (2014). First, we describe the geometrical and kinematics parameters before introducing the inversion procedure and the results.

Megathrust Geometry and Medium Rheology
The geometry of the slab was built considering the position of the trench axis and the depth contours of the Slab 2 geometry model from Hayes et al. (2018). This new slab geometry shows significant differences with the geometry used in previous studies Font et al., 2013;Hayes et al., 2012). The geometry of the Slab 2 model is much steeper and about 10-15 km deeper than the Font et al. (2013) and Slab 1 geometries between latitudes −1°S and 2°N. Slab 2 geometry becomes shallower than the Slab 1 model along the Colombian coast between 2°N and 4°N (Supplementary Figure S4). Our modeled megathrust interface here is composed of 3 juxtaposed rectangles extending along strike from southern Ecuador northward 240 km, 240, and 480 km, respectively, and about 300 km along the megathrust dip (Supplementary Table S2). The average dip angle measured from the trench axis located at 7 km below sea level to slab iso-depth contour of 50 km increases northward from 18°, to 20°and 22°, respectively. We discretize the megathrust interface into 700 square elements of 20 km × 20 km with a source point fixed at the center of each element. The megathrust geometry is then embedded in a Poisson solid-layered half-space (Supplementary Table S3) with a Poisson coefficient of 0.25 and an average shear modulus µ of 36 GPa computed between 0 and 50 km depth of the subduction wedge. With that formalism, static displacements are computed following the analytical expressions of (Ji et al., 2002;Xie and Yao, 1989). The slip rake is authorized to vary ±10°from the average rake direction given by local gCMT slip vectors. When all the parameters described above are set, the full representation of the fault response relies principally on the dislocation amplitude parameter and on the rake angle that is partially controlled. We can invert these parameters by fitting the modeled static displacements to the observations considering a reduced chisquare χ r 2 criterion defined as: where n is the number of data, obs and pred are the observation and the prediction, respectively, and σ is the standard deviation error associated with each GPS component.

Inversion Procedure
The general forward problem is written as follows: where G represents the Green functions linking the observables to the model parameters, d is the data vector and m is the model parameter vector (here, the back-slip amplitude and rake direction at each of the 700 source points). In each inversion, we weighted the GPS measurements by their respective 1-sigma uncertainties (Supplementary Table S1). The inverted back-slip distribution follows two types of regularized constraints: one that minimizes the difference between the slip on adjacent sub-faults and a second that minimizes the final moment (or moment deficit rate in case of interseismic). The cost function to minimize is: where Dc represents the differences in the slip (or back-slip for interseismic) of adjacent cells, Mo apriori is an a priori moment deficit rate, and Mo is the final moment. The coefficient λ 1 controls the smoothing coefficient through an L1+L2 norm, and the coefficient λ 2 modulates the weight assigned to minimize the final moment deficit rate. We found that for λ 1 0.1 the fit to the GPS data is optimal without over-estimating the roughness of the final slip solution (Supplementary Figure S5).

Back-Slip Inversions
For the interseismic modeling, we applied a back-slip approach (Savage, 1983) following the inversion procedure described above. We bound the back-slip rate (V back ) to the relative Nazca/NAS plate convergence rate (V pl ) of 50 mm/yr in the southern Ecuador, of 48 mm/yr in the central-north Ecuador, and of 46 mm/yr in Colombia. The local interseismic coupling (IC) is then defined as: resulting in pure creeping on the megathrust interface, when V back V pl (IC 0), and in full locking of the plate contact, when V back 0 (IC 1). Using that definition, the moment deficit rate integrated over a specific fault surface S embedded in rigid medium µ is: where V pl -V back is the back-slip offset rate across the fault. Inversions of the GPS velocities provide a highly heterogeneous interseismic coupling map at the scale of the Ecuador-Colombia subduction (Figures 3-5). Considering the 1-sigma error uncertainties of the GPS components and solutions with a reduced Chi-square χ r 2 <1, we found that models with an average moment deficit rate of 8.5 ± 1.5 × 10 19 N-m/yr (Supplementary Table S4) for the whole Colombia-Ecuador subduction zone provide the best fit to the GPS data ( Figure 4, Supplementary Figures S6,  S7). The rate of moment deficit corresponds to the equivalent of Mw∼7.2 ± 1 earthquake that is accumulating annually on the megathrust interface from Guayaquil to Buenaventura. Models with a lower moment deficit rate would degrade the fit to the data and reduce spatially the locking at the centers of highly locked patches (Supplementary Figure S6). Models with a higher moment deficit rate would also increase the data misfit (Supplementary Figure S7) and saturate the coupling near the trench (Supplementary Figure S6). Because the shallow portion of the plate interface lies in the stress shadow of deeper-locked asperities (Bürgmann et al., 2005;Chlieh et al., 2008), it has been for long time impossible to assess the level of locking in the shallower portion of the slab interface without near-trench seafloor observations (Chlieh et al., 2011;Gagnon et al., 2005). A promising approach suggests adding some physical constraints to Frontiers in Earth Science | www.frontiersin.org September 2021 | Volume 9 | Article 701720 evaluate more precisely the level of locking that remains in the shallow stress shadow of deeper asperities (Almeida et al., 2018;Herman and Govers, 2020;Lindsey et al., 2021). Our average and maximum solutions, shown in the inset C and D of Figure 4, share some similarities with this new approach and indicates that the shallow coupling can be higher than 0.4 nearly everywhere and even reach 1 near the trench.
Source Inversion of the 2016 Mw = 7.8 Pedernales Earthquake To analyze the seismic source of the 2016 Mw 7.8 Pedernales earthquake, we inverted the horizontal and the vertical coseismic displacements of permanent and campaign GPS stations of the Ecuadorian GPS network (listed in Supplementary Table S5) published by Nocquet et al. (2016) using the same procedure described above. In the rupture area of the Pedernales earthquake, there are significant differences between slab geometries that we explored by running series of dislocations geometries in which the strike, dip, and rake angles are varied (Supplementary Figure S8). We found tan optimum fit to the coseismic GPS data for rupture geometries with an average strike angle of 25°± 5°, dip angle of 18°± 3°, and rake angle of 120°± 5°. These angles are relatively well consistent with the 27°strike, 21°dip, and 124°rake angles determined by the gCMT focal mechanism and also closer to the average dip angle of Slab 2 geometry. Modeled focal mechanism using the Slab 2 geometry and observed gCMT are relatively well consistent and shown in Figure 6A. Then, using the same Slab 2 megathrust geometry, layered medium and inversion procedure for the interseismic models, we fixed the seismic moment to magnitude Mw 7.8 (with λ 2 > 100) and varied the smoothing factor λ 1 from 0.01 to 10 (see Supplementary Figure S7). The fit to the GPS data becomes optimal for a smoothing coefficient λ 1 lower than 1. Best-GPS fitting slip distributions exhibit either one slip patch that peak up to 3.5 m or two slip patches for rougher solutions, typically when λ 1 ≤ 0.01 (see Supplementary Figure S7). Many publications of the seismic source of the 2016 Mw 7.8 Pedernales event suggest either one slip patch when the inversions is derived only from geodetic (InSAR or GPS) data (He et al., 2017) or two slip patches when the inversions include seismic waveforms and/ or tsunami data (Gombert et al., 2018;Heidarzadeh et al., 2017;Nocquet et al., 2016;Ye et al., 2016;Yi et al., 2018). Here, for consistency with the interseismic modeling and with the goal to balance properly each phase of the seismic cycle using the same inversion procedure, slab geometry and boundary conditions, we decide to retain the coseismic slip distribution for the smoothing coefficient λ 1 0.1 ( Figure 6A) that fit very well the coseismic GPS data ( Figure 6B) and that is the same than for interseismic models.

On the Ecuador-Colombia Interseismic Coupling Map
The pattern of interseismic coupling allows us to distinguish with confidence several persistent asperities of about 50-100 km width and with a relatively high level of coupling typically above 0.6. From south to north, these asperities appear on the megathrust interface offshore Punta Salinas, below La Plata island, between the Bahia de Caraquez and the Atacames promontory (Pedernales asperities), offshore Esmeraldas, and offshore the Tumaco Peninsula ( Figure 5). Further north in Colombia, between the Sanquianga Peninsula and Buenaventura, it is less clear whether there are two 50 km-wide asperities located between the slab depths of 20 and 50 km or if the coupling is much higher on a large asperity that extends diagonally from about 10 km depth offshore the Sanquianga Peninsula to 40 km depth offshore Buenaventura. The lack of resolution of the shallow megathrust in this area leaves a high uncertainty in the spatial location of the Cali-Buenaventura asperities; however, it appears quite clear that the Cali-Buenaventura asperities accumulate a significant amount of moment deficit ( Figure 5) and that the locking never extends below a slab depth of 20 km offshore the Sanquianga Peninsula. These asperities are disconnected by creeping segments with widths varying along strike from 20 to 80 km that are either fully creeping (level of coupling 0) or partially creeping with a level of coupling typically lower than 0.4 (Figures 4, 5). On the Coseismic Source Figure 6A shows one of the best GPS fitting solutions with a seismic source that reaches a maximum slip of 3.5 m at about 25 km depth. The maximum slip amplitude is consistent with the 3.55 m predicted by the finite plate convergence considering a relative Nazca/NAS convergence rate of 48 mm/yr between 1942 and 2016, suggesting no overshoot. Inverted slip rake directions are collinear with the gCMT slip rake directions of background seismicity and aftershocks that have preceded and followed the Mw 7.8 Pedernales earthquake (Figure 6 and Supplementary Figure S1). These rake directions are also relatively well consistent with the relative Nazca/NAS plate convergence vector, suggesting no significant slip partitioning during that earthquake. One of the most striking features of the gCMT aftershocks are two normal-trench alignments south of the Atacames promontory and offshore the Bahia de Caraquez ( Figure 6 and Supplementary Figure S1). These two aftershock ribbons bound laterally the 2016 Mw 7.8 seismic source, supporting the presence of two creeping transition zones there. Many recent studies document very well that aftershock sequence including large aftershocks, seismic swarms, and slow slip events activated by the afterslip on the megathrust interface that occurred during the postseismic period, but also crustal seismic events associated with the reactivation of shallow crustal faults (  (B-D) shows the minimum, the average, and the maximum interseismic coupling models, respectively, with the fully locked asperities (in red) and the fully creeping segments in white. Overall, the interseismic coupling remains confined within the first 50 km depth of the slab interface with some robust characteristics that remain stable in all models (i.e., large asperities location and creeping segments) and significant variations in areas of low spatial resolution (shaded gray rectangle) as deduced from the checkerboard tests shown in Supplementary Figures S3A,B.
Frontiers in Earth Science | www.frontiersin.org September 2021 | Volume 9 | Article 701720 geomorphological features as ridges or seamounts Marcaillou et al., 2016). The spatial distribution of our 2016 seismic source correlates relatively well with the location of a highly-coupled and persistent asperity well identified in our interseismic models ( Figure 5) and located right below the coastline between 20 and 40 km depth of the megathrust interface. The rupture extremities of the 2016 Pedernales earthquake are laterally well bounded by two creeping segments: the Carnegie creeping segment to the south and the Mompiche creeping segment to the north (Figures 5, 7) that may have played the role of seismic barriers to the rupture propagation. The reactivation of these areas during the postseismic period is an additional evidence that support the presence of these creeping segments. Because our knowledge of the 1942 seismic source is quite limited compared to the recent 2016 seismic source analysis, it is hard to say whether these two large seismic events unambiguously broke the same locked asperity, however, some seismic and geodetic evidences, let suppose that, it is likely the case. Between the Atacames promontory and Bahia de Caraquez located at 0.5°N FIGURE 5 | Seismic segmentation of the Ecuador-Colombia subduction zone from the interseismic coupling map. (A) High moment deficit rates (in N.m/yr) indicate the location of highly locked asperities since low moment deficit rates highlight the presence of creeping segments. (B) Average interseismic coupling map resulting from geodetic inversion with a 68% confidence level. Creeping segments of ∼80 km-width and with an average level of coupling lower than 0.2 coincide with the entrance into the subduction of the Carnegie Ridge in Ecuador and the young oceanic crust east of the Yaquina Graben in Colombia. These two wide creeping segments have the characteristics of relatively strong seismic barriers able to stop great rupture like the great 1906 Colombia-Ecuador earthquake. Thinner creeping segments of about 20-40 km large appear between Tumaco and Esmeraldas and offshore Mompiche.
Frontiers in Earth Science | www.frontiersin.org September 2021 | Volume 9 | Article 701720 and 0.5°S, respectively, the interseismic coupling map indicates two distinctively locked asperities on the megathrust (Figure 7). One asperity is shallow and extents from the trench to about 20 km depth, and the second asperity is deeper and extends from about 20 to 40 km depths. The location of this deeper asperity overlaps relatively well with the coseismic source distribution of the 2016 Pedernales event ( Figure 7) and overlap with the rupture and aftershocks area of the 1942 Mw 7.8 earthquake Swenson and Beck, 1996). It is unlikely that one of these seismic events did rupture the shallow asperity without producing an important local tsunami as it was observed from well identified tsunami-earthquakes with similar or lower moment magnitude during the 1996 Mw 7.5 Chimbote event in Peru (Ihmlé et al., 1998) or the 1994 Mw 7.6 and 2006 Mw 7.8 Java events in Indonesia (Abercrombie et al., 2001;Ammon et al., 2006). The relatively low tsunami generated in 1942 and 2016 rather support that both ruptures are associated with the failure of the deeper mapped asperity. Assuming a steady-state interseismic moment deficit rate, we found that over the 2016 rupture area, the moment deficit buildup since 1942 is equivalent to 1.3 ± 0.1 10 21 N-m which is nearly twice the seismic moment released during the 2016 event (Figure 7). This result is in contradiction with recent studies that suggest an overshoot of slip and then excess of seismic moment released through earthquakes since 1906 in that region (Gombert et al., 2018;Nocquet et al., 2016).
Postseismic afterslip associated with the 2016 Pedernales event occurred up-dip of the coseismic rupture area and has released a postseismic moment of 0.16 × 10 21 N-m which represents about 23% of the 2016 seismic moment (Rolandone et al., 2018). From a high-rate GPS analysis, it has been shown that the early postseismic afterslip that occurred in the 12 to 72 first hours following the initiation time of the 2016 rupture would increase the postseismic/coseismic moments ratio to about ∼30% (Tsang et al., 2019;Twardzik et al., 2019). Then, the cumulative 2016 coseismic and postseismic moment would represent about 60-70% of the moment deficit buildup between 1942 and 2016. This suggests that the Pedernales asperities were well mature to host the 2016 Mw 7.8 earthquake and its associated afterslip. The remaining moment deficit can be released by the background seismicity and/or aseismically through slow slip episodes as reported in many places of the Ecuadorian subduction zone Rolandone et al., 2018;Vaca et al., 2018;Vallée et al., 2013). It is also possible that a fraction of that remaining moment deficit is simply stored for future large subduction earthquakes and long-term topography building (Jolivet et al., 2020;Melnick, 2016;Saillard et al., 2017).
With similar steady-state interseismic loading rate assumptions, we found that the Esmeraldas asperity that failed during the 1958 Mw 7.6 event has built up a cumulative subduction earthquake derived from the joint inversion of horizontal and vertical GPS data using Slab 2 megathrust geometry. The color indicates the slip amplitude and the gray arrows the inverted slip rake directions at each source points that mesh the megathrust interface. Inverted slip rake directions are collinear with the gCMT slip rake directions of background seismicity (black bars) and aftershocks (green bars) that have preceded and followed the Mw 7.8 Pedernales earthquake (see the gCMT catalog in Supplementary Figure S1). These rake directions are also collinear with the relative Nazca/NAS plate convergence vector, suggesting no significant slip partitioning during that earthquake. The modeled focal mechanism (green beach-ball) is well consistent with the observed gCMT (red beach-ball). (B) Observed horizontal (black arrows) and vertical (red arrows) coseismic GPS displacements compared to the modeled displacements (green arrows) associated with the source model shown in A. Slip contours of the seismic sources are reported every 0.5 m.
Frontiers in Earth Science | www.frontiersin.org September 2021 | Volume 9 | Article 701720 1958-2020 moment deficit of 0.9 ± 0.1 10 21 N-m (∼Mw 7.9) and that the Tumaco asperity that failed during the 1979 Mw 8.1 event has built up a cumulative 1979-2020 moment deficit of 0.7 ± 0.2 10 21 N-m (∼Mw 7.8). Considering that only half of these moment deficits is released seismically through large earthquakes as observed in the Pedernales segment, then the Esmeraldas asperity would be mature to host a Mw∼7.7 event and the Tumaco asperity would be mature for a Mw∼7.6 event (Figure 7). If we consider that the two asperities fail simultaneously, a Mw∼7.9 seismic event could be generated. Alternative scenarios can be built considering the moment rate deficit associated with each asperity listed in Supplementary  Table S4 with an uncertainty determined at 1-sigma (68% confidence level) and considering different fractions of aseismic/seismic moments.

Segmentation of Large and Great Megathrust Earthquakes and Barriers Efficiency
The GPS-derived average interseismic coupling model indicates the locations of large apparent asperities associated with a relatively high rate of moment deficit. These asperities are bounded by creeping segments of few tens of kilometers where the rate of moment deficit is relatively low. These creeping segments promote velocity Frontiers in Earth Science | www.frontiersin.org September 2021 | Volume 9 | Article 701720 strengthening on the megathrust interface which in the longterm play the role of seismic barrier able to slow down and inhibit seismic rupture propagation. The barrier efficiency of a velocity strengthening patch has been shown to be proportional to its along-strike length and an average level of interseismic coupling (Kaneko et al., 2010). Then, large creeping barriers with low interseismic coupling will be much more efficient than smaller barriers and would stop systematically all seismic ruptures, even those associated with great earthquakes like in 1906. Smaller barriers have a dual frictional behavior and will work efficiently in case of a large Mw∼7.5-8.0 seismic rupture but will fail during great Mw > 8.5 rupture that break, simultaneously, many large juxtaposed asperities. These smaller barriers will slow down the rupture speed but will not be able to totally inhibit it as it was observed in Peru during the 2001 Mw 8.4 Arequipa and the 2007 Mw 8.0 Pisco earthquakes (Perfettini et al., 2010;Robinson et al., 2006). In Ecuador, the Carnegie Ridge creeping segment is characterized by an along-strike length of about 50-60 km and an average interseismic coupling lower than 0.2 and has been considered as a strong and efficient seismic barrier Yepes et al., 2016). This is supported by the fact that the 1906, 1942, 1998, and 2016 ruptures all abut at the northern extremity of the Carnegie creeping segment. In Colombia, a similar creeping segment that shares the same length and frictional properties with the Carnegie Ridge appears offshore the Sanquianga peninsula ( Figure 5). Offshore Colombia, from the Tumaco Peninsula to Buenaventura, the gCMT catalog indicates that the outer rise seismicity located between the Yaquina Graben and the trench is dominated by trench-parallel normal faulting (see Supplementary Figure S1). This seismicity is due to reactivation of normal faults associated with buried extinct rift that intersects the Colombian accretionary complex at latitude 3.5°N (Hardy, 1991;Lonsdale, 2005). Normal faulting crosses the trench and follows the slab geometry suggesting a highly fractured subducting slab interface (Collot et al., 2002;Yepes et al., 2016). The roughness of the megathrust interface induced by the subduction of the Carnegie ridge in Ecuador and the highly fractured oceanic crust east of the Yaquina Graben in Colombia, coupled with poorly consolidated sediment at the plate contact may promote aseismic slip, a process governed by rate-strengthening friction . This suggests that the Sanquianga creeping segment can also be considered as a strong and efficient seismic barrier and might probably have limited the northern extent of the 1906 rupture.
If the Sanquianga creeping segment bounded the northern extremity of the great 1906 Colombia-Ecuador earthquake, then the 1906 great rupture would have not exceeded an along-strike length of 400 km distance measured along the trench axis from the northern edge of the Carnegie Ridge creeping segment to the southern edge of the Sanquianga creeping segment ( Figure 5). Comparatively, the 2005 Mw 8.6 Nias earthquake in Indonesia ruptured a segment of ∼350-400 km of the Sunda megathrust (Chlieh et al., 2008;Hsu et al., 2006). Along the South America subduction zone, the 2001 Mw 8.4 Arequipa earthquake in Peru ruptured a segment of ∼300-350 km (Chlieh et al., 2011;Pritchard et al., 2007), and the 2010 Mw 8.8 Maule earthquake in Chile ruptured a segment of ∼500-550 km (Delouis et al., 2010). By contrast, "weaker" barriers as the Mompiche and the Tumaco-Esmeraldas creeping segments appear no longer than about 20-40 km and with an average interseismic coupling lower than 0.4 ( Figure 5). The proximity between the locked asperities that have respectively failed in 1942/2016, 1958, and 1979 may produce a lateral stress shadow that is translated as saddle-like shape in the lateral variation of the moment deficit rate ( Figure 5A). It appears that these creeping segments have a dual behavior since they have worked efficiently during the large earthquakes of 1942, 1958, 1979, and 2016 but have both failed during the great 1906 rupture. Aftershocks that have followed the 2016 Mw 7.8 Pedernales earthquake cluster along trench-normal alignments that fit quite well with the Mompiche creeping segment supporting the presence of some afterslip and then aseismic slip there (Agurto-Detzel et al., 2019;Ye et al., 2016). In addition, a 2-year recurrent slow slip events that trigger seismic swarms on the megathrust interface have been reported in that same segment (Vaca et al., 2018), suggesting that the apparent creeping behavior revealed by the coupling models is not a continuous process but is likely the result of higher coupling during inter-SSE periods that alternate with recurrent transient creeping episodes that contribute to release locally the accumulated moment deficit.

Seismic Catalog and Balance Between the Seismic Rate and the Moment Buildup Rate
To model the magnitude-frequency distribution along the Colombia-Ecuador subduction zone, we used a unified and homogeneous earthquake catalog ( Figure 8A) recently updated in Beauval et al. (2018), which is based on instrumental and early instrumental earthquake catalogs (Beauval et al., 2018;Beauval et al., 2013;Yepes et al., 2016). This catalog is built from four global earthquake catalogs: the ISC event catalog, the ISC-Global Earthquake Model (GEM) catalog covering the time window 1900-2013 (Storchak et al., 2015), the Global Centroid Moment Tensor (gCMT) catalog covering the time window 1976-2014, and finally the NEIC catalog. For earthquakes with mb and Ms teleseismic magnitudes, Mw proxies are estimated applying conversion equations of (Lolli et al., 2014). The final catalog contains 2,531 events with Mw ≥ 4.2. For the first half of the century , only magnitudes down to Mw 5.5 are reported in the catalog. More details on the building of the catalog can be found in Beauval et al. (2018). For our specific analysis of the slip budget in the seismogenic zone of the subduction interface, we excluded outer rise seismicity and selected the earthquakes located in the 0-60 km depth of the Colombia-Ecuador megathrust interface.
The cumulative seismic moment rate is also calculated by adding the moment rates of earthquakes estimated over their completeness period. Assuming that the earthquake catalog is Frontiers in Earth Science | www.frontiersin.org September 2021 | Volume 9 | Article 701720 complete with relatively well-constrained magnitudes and that the interseismic model derived from decadal geodetic measurements reflects the average long-term interseismic buildup, we balance the interseismic moment buildup rate and the seismic moment released rate. In Figure 8B, the red line show the cumulative seismic moment release rate in the period 1900-2016 and the blue lines show the moment buildup rates calculated from our best interseismic coupling models shown in Figure 4. We found that, during the time period 1900-2016, seismicity does not balance locking in the interseismic period. One great or many large megathrust earthquakes are needed to close the budget. Alternative possibilities would invoke that aseismic slip is more prominent or that interseismic moment buildup may vary with time.

Seismic Moment Balance at the Scale of the 1906 Colombia-Ecuador Subduction Zone
Slip budget closure is reached when the long-term slip at any point on the fault equals the summation of seismic and aseismic slips (Avouac, 2015). This principle of moment balance is verified when the cumulative interseismic moment buildup balances the sum of transient seismic and aseismic moments released either through earthquakes, afterslip, or spontaneous SSEs. The moment balance can be then written as: where Mo Deficit is the cumulative moment deficit associated with the interseismic loading, Mo Seismic the cumulative seismic moment, and Mo aseismic the cumulative aseismic moment due to postseismic afterslip (and various slow slip events that occur during the interseismic and pre-seismic periods). Assuming that the frequency-magnitude distribution of megathrust earthquakes follows the Gutenberg-Richter (G-R) law on the long-term (Gutenberg and Richter, 1944): where N is the number of earthquakes larger than Mw, Mw 2/3 log10 (Mo)-6.07 with Mo in N.m (Kanamori, 1977), and M max the maximum magnitude on the megathrust. The slope of the distribution is the b-value and a describes the productivity. Different forms can be assumed for predicting the rates in the upper magnitude range close to M max (Mariniere et al., 2021); however, in this study only this form is investigated.
To derive a long-term earthquake recurrence model of the Colombia-Ecuador subduction zone, we selected subduction earthquakes in a box that includes the 0-60 km depth of the megathrust interface (purple dots in Figure 8A). Earthquake rates are obtained by normalizing the catalog by its timespan. The frequency-magnitude distribution indicates a b-value of 0.93 ( Figure 8C), computed over the magnitude interval 5.2-7.5 considering the whole Colombia-Ecuador subduction zone from latitudes 3°S to 4°N. FIGURE 9 | Moment budget closure for the Ecuador-Colombia subduction zone for various seismic/aseismic ratio. The blue line represents the frequency of the maximum magnitude earthquake to fulfill the moment budget closure condition. The red lines follow the GR law with a b value of 0.93 (found in the instrumental catalog, Figure 8C) and show where the seismicity should lie given an earthquake with maximum magnitude of 6, 7, 8, 9, or 10 from top to bottom, respectively. We assume different amount of seismic vs. aseismic slip: (A) α 1, transient slip is fully seismic, (B) α 0.77 (30% of the transient slip is aseismic), (C) α 0.5 (transient slip is half seismic and half aseismic), and (D) α 0.37, lower α limit to be consistent with an Mw 8.6 event like in 1906.

Maximum Magnitude Earthquake and Recurrence Time
In a region where the seismicity follows a Gutenberg-Richter law, the maximum magnitude earthquake (M max ) can be determined looking for the upper bound of the recurrence model that implies a seismic budget closure (Stevens and Avouac, 2017). Knowing the a-and b-values, the seismic budget closure determined by the interseismic moment buildup rate gives the necessary constraint to find M max . Assuming that the frequency-magnitude distribution follows Equation 7, the mean recurrence time of the maximum magnitude earthquake M max is defined as (Avouac, 2015;Molnar, 1979): where α is the fraction of transient slip that is seismic during the seismic cycle and is written as: α Mo seismic /(Mo seismic +Mo aseismic ), Assuming a b-value of 0.93 found from the instrumental catalog and extrapolated for larger magnitudes and an average moment deficit rate (dMo inter /dt) of 8.5 ± 1.5 × 10 19 N m/yr, we explore what would be the maximum magnitude earthquake for various α between 0.1 and 1 (Supplementary Figure S10 and Supplementary Table  S6). As the minimum bound for M max is the magnitude of the largest earthquake observed in that subduction zone, that is, the 1906 Mw 8.6 event, then we found that α cannot be lower than 0.37. When α < 0.37, the frequency-magnitude distribution of the seismicity overshoots the frequency of the largest earthquake and the maximum magnitude earthquake would be lower than the 1906 event. For α 0.1, 02, or 0.3 the model would predict a M max of 7.6, 8.1 and 8.4, respectively, which is not consistent with the M max of the historical seismic catalog. This indicates that aseismic slip cannot be more than 1.7 larger than seismic slip over the long-term. Alternatively, if aseismic slip is considered as negligible compare to seismic slip, then α would tend to 1. In the specific case where α 1, we found a maximum magnitude event Mw max of 9.36. This case is certainly not valid since it does not take into account any aseismic slip. We know from the 2016 Pedernales postseismic afterslip and many other postseismic afterslip examples in different subduction zones that the postseismic afterslip would release an aseismic moment of about 10-30% of the coseismic moment Hsu et al., 2006;Perfettini et al., 2010;Remy et al., 2016). Considering that Mo afterslip corresponds to about 30% of Mo seismic , like in the case of the 2016 Pedernales event, then α that cannot be higher than 0.77. In that specific case, we found a maximum magnitude event M max of 9.16.
Moment budget closure are conducted for specific alpha value. In Figure 9, the red lines represent the seismicity rates required for the maximum magnitude earthquake to be equal to 6, 7, 8, 9 and 10 respectively from top to bottom. The blue line represents the frequency of the largest earthquake derived from Equation 8 with the geodetically-derived moment deficit rate. The intersection between the blue line and the red line fitting the moderate-to-large earthquakes indicates the largest magnitude earthquake that could occur on the whole Ecuador-Colombia subduction zone. If we consider that in the long-term, the aseismic moment released equals the seismic moment released, then α would equal to 0.5 and M max 8.8 with a recurrence time of about 1,400 years. In that specific case, the model would predict that the recurrence time of an event equivalent to the 1906 Mw 8.6 would be of 635 ± 110 years, a recurrence time that is quite consistent with paleo-seismic records from a core collected in the trench offshore Esmeraldas that indicate that a 1906-type event have occurred ∼600 years ago (Migeon et al., 2017).

CONCLUSION
We Found 1. The 1906 rupture length probably did not exceed ∼400 kmlong and appears to be bounded by two 80 km-large creeping segments that coincide with the subduction of the Carnegie Ridge in Ecuador and of the young and fractured oceanic crust east of the Yaquina Graben in Colombia (Lonsdale, 2005;Yepes et al., 2016). 2. The 2016 Mw 7.8 Pedernales earthquake appears to be the successor of the 1942 Mw 7.8 event since they both broke the same asperity of the Pedernales segment 3. The 2016 Mw 7.8 Pedernales earthquake and its associated postseismic afterslip released a cumulative moment of about 60-70% of the moment build-up since the last rupture of that segment in 1942 suggesting that this area was well mature to host the 2016 seismic event and its afterslip. 4. In the aftermath of the 2016 Pedernales event, the 1958 rupture area is mature to host an Mw 7.7 to 7.9 event and the 1979 rupture area is mature for a Mw 7.6 to 7.8 event. The simultaneous rupture of these two segments could host a Mw∼7.8-8.0 earthquake. 5. The cumulative seismic moment of the 1942, 1958, 1979 and 2016 earthquakes represents ∼35% of the 1906 seismic moment. 6. In the long-term, if 30% of the slip is aseismic, the moment balance suggests that the maximum magnitude event in this subduction zone could be around ∼9.2 with a mean recurrence time around 2,800 years. 7. In the long-term, if 50% of the slip is aseismic, the moment balance suggests that the maximum magnitude event in this subduction zone could be around ∼8.8 with a mean recurrence time around 1,400 years. In that case, the recurrence time of Mw 8.6 like in 1906 would be 635 ± 110 years.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material; further inquiries can be directed to the corresponding author.