Strain budget of the Ecuador–Colombia subduction zone: A stochastic view

a Institut de Physique du Globe de Strasbourg, UMR7516, Université de Strasbourg, EOST/CNRS, France b Laboratoire de géologie, Département de Géosciences, École Normale Supérieure, PSL Research University, CNRS UMR 8538, Paris, France c Seismological Laboratory, Geological and Planetary Sciences, California Institute of Technology, Pasadena, CA, USA d Institute of Geophysics and Planetary Physics, University of California, San Diego, La Jolla, CA, USA e Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA, USA


Introduction
A long standing question is the existence of persistent fault segments remaining locked in the inter-seismic period and failing suddenly during earthquakes while the surrounding interface creeps continuously. This conceptual model predicts so-called "character-  1942, 1958, 1979, and 1998events (Kanamori and McNally, 1982Chlieh et al., 2014). The location of these previous ruptures is still debated, and some alternative plausible locations are shown in Fig. 10 and S9. The thick gray line shows the along-strike extent of the 1906 M W = 8.6 earthquake. The focal mechanism of the 2016 Pedernales earthquake is presented in blue. Thin black lines are isocontours of the slab depth. The line with the adjacent black triangles shows the location of the trench. The black arrow illustrates the convergence direction of the Nazca plate toward the North Andean Silver plate (NAS, Chlieh et al., 2014). (For interpretation of the colors in the figure(s), the reader is referred to the web version of this article.) long segment of the subduction interface (Gutenberg and Richter, 1949;Ye et al., 2016). Several decades later, the same area was re-ruptured by a series of smaller M W ≤ 8.2 events in 1942, 1958, 1979and 1998(Kanamori and McNally, 1982Beck and Ruff, 1984;Mendoza and Dewey, 1984;Chlieh et al., 2014). In April 2016, the region in the vicinity of the 1942 Ecuador event was again ruptured by the M W = 7.8 Pedernales earthquake (Ye et al., 2016;He et al., 2017;Nocquet et al., 2017;Yi et al., 2018). Such variability among successive ruptures is also observed in other regions like Japan and Sumatra where recent M W ∼ 9 megathrust earthquakes ruptured large fault segments that previously experienced a series of smaller events (Simons et al., 2011;Lay, 2015).
In addition to such spatial variability among successive ruptures, major earthquakes in the Colombia-Ecuador subduction zone seem to be clustered in time. Specifically, it has been recently suggested that the seismic moment of the 1942, 1958 and 1979 earthquakes exceeds the deficit accumulated since 1906 and that the 2016 Pedernales event may be associated with more fault slip than the deficit accumulated since the 1942 earthquake . Similar observations are reported in other regions, for example in 1797 and 1833 earthquakes in Sumatra (Sieh et al., 2008), 1812 and 1857 earthquakes in California (Jacoby et al., 1988;Heaton, 1990), and for the 2003 M W = 7.6 and 2013 M W = 7.8 Scotia sea earthquakes (Vallée and Satriano, 2014). Such spatial and temporal clustering can be caused by spatial variations of fault coupling associated with heterogeneous frictional properties (Kaneko et al., 2010). Moreover, there can be fluctuations in the patterns of inter-seismic fault coupling before large earthquakes (Perfettini and Avouac, 2004;Mavrommatis et al., 2014;Yokota and Koketsu, 2015) or during the post-seismic response of nearby large earthquakes (Heki and Mitsui, 2013;Melnick et al., 2017).
Although the existence of an anomalously large co-seismic slip associated with a supercycle behavior is plausible, other studies suggest that the seismic moment of the 2016 Pedernales earth-quake is actually consistent with the strain accumulated in the region since the 1942 and 1906 earthquakes (e.g., Ye et al., 2016;Yoshimoto et al., 2017;Yi et al., 2018). These contrasting statements partly results from the ill-posed nature of inter-and coseismic slip inversions used to evaluate the strain budget along the megathrust. Such inferences are affected by the lack of resolution near the trench during the inter-seismic period but also by non-physics-based smoothing constraints used to regularize slip inversions. In addition, inter-and co-seismic estimates usually do not incorporate rigorous uncertainties (or very often, no uncertainty at all), which complicates a quantitative assessment of the overall strain budget. Strain budget analyses also suffer from the lack of information about past earthquakes (Yi et al., 2018). Incorrect considerations on the size and position of historical events can strongly affect the conclusion on the strain state of the plate boundary.
We propose a probabilistic exploration of the Colombia-Ecuador earthquake sequence, fully accounting for uncertainties, including measurement errors, modeling errors, but also uncertainties in the location or magnitude of past events. Using a Bayesian framework, we explore both the inter-seismic geodetic coupling of the subduction interface and the co-seismic slip distribution of the M W = 7.8 Pedernales earthquake. These estimates do not rely on any spatial smoothing and provide full posterior probability distributions describing the ensemble of plausible models that fit the observations and are consistent with simple prior constraints (e.g., slip positivity in the direction of convergence).

Stochastic inter-seismic modeling
We first compute a stochastic model of geodetic coupling along the Ecuadorian subduction interface. We use inter-seismic GPS velocities computed by Chlieh et al. (2014) and Nocquet et al. (2014) from 29 stations installed in Ecuador and Colombia and measured from 1994 to 2012 Mora-Páez et al., 2018). Considering that more than 20 years separate the 1979 earthquake and the first GPS measurements, it is unlikely that they are significantly affected by post-seismic deformation. The 1998 earthquake having a smaller magnitude, its impact on the data is also probably minimal. The fault geometry is based on a 3D surface following the Slab1.0 interface and discretized in triangles (cf. Fig. S1 in the electronic supplements). Using a back-slip approach (Savage, 1983), we invert for the inter-seismic slip rate along the direction of convergence between Nazca and North Andean Sliver (NAS) plates at each of the triangle knots assuming a barycentric interpolation scheme within the triangles. This approach avoids unphysical slip discontinuities associated with traditional parameterizations based on sub-faults with piecewise constant slip (Ortega Culaciati, 2013).
In our Bayesian inversion framework, the solution is the posterior ensemble of all plausible inter-seismic slip rate models (m I ) that fit the GPS data (d I ) and that are consistent with our prior hypotheses. This solution does not rely on any smoothing regularization and is based on a simple uniform prior for the inter-seismic slip-rate that writes p(m I ) = U (−0.05 · V p , 1.05 · V p ) M where V p is the plate rate and M is the number of triangle knots (260 knots). We thus restrict our posterior PDF to models in which slip on the fault aligns with the direction of plate motion. Following Bayes' theorem, the posterior PDF is given by where G I is the Green's function matrix and C I is the misfit covariance matrix combining observational errors and prediction uncertainties. Green's functions are computed for a semi-infinite stratified elastic medium derived from regional velocity models shown in Fig. S2 (Béthoux et al., 2011;Vallee et al., 2013;Nocquet et al., 2017). We account for prediction uncertainties due to inaccuracies in this layered model using the approach of Duputel et al. (2012Duputel et al. ( , 2014. The uncertainty on the elastic structure, presented as grey histograms in Fig. S2, is estimated by comparing previously published models in the region. We sample the posterior PDF p(m I |d I ) using AlTar, a parallel Markov Chain Monte Carlo (MCMC) algorithm following the CATMIP algorithm (Minson et al., 2013). More details on the application of AlTar to investigate inter-seismic deformations can be found in Jolivet et al. (2015b) and Klein et al. (2017). The resulting posterior ensemble of slip-rate models in eq. (1) is then converted into stochastic coupling maps (m C ) using m C = 1 − m I /V p .

Geodetic coupling results
Using our Bayesian framework, we generate 160 000 models corresponding to the posterior information on geodetic coupling given measured inter-seismic velocities. We find that this number is large enough to converge toward the posterior probability density. Representing the ensemble of posterior models is challenging for multidimensional problems such as those addressed in this study. To represent an ensemble solution, a common choice is to compute the posterior mean (i.e., the average of all model samples). The posterior mean coupling model is shown in Fig. 1 and Several features in our solution can be observed in previously published geodetic coupling models (e.g., Nocquet et al., 2014;Chlieh et al., 2014). In the South, there is a very clear highcoupling area offshore the Manta peninsula. This inter-seismically highly coupled region has been previously associated with transient slow-slip events (Vallee et al., 2013;Nocquet et al., 2014). As shown in Fig. 2a and Fig. 2e, this area is associated with small model uncertainties probably because a GPS station is located on La Plata Island, right above the coupled asperity. This coupled patch is bounded to the north by a low-coupling corridor that might have acted as a creeping barrier for the 1906,1942,1998 and 2016 earthquakes (cf. Fig. 1; Chlieh et al., 2014).
North of Bahía de Caráquez, we infer multiple patches of high geodetic coupling. Other coupled patches can be identified offshore of Bahía de Caráquez, North and South of Pedernales, and far offshore Esmeraldas. To first order, such heterogeneity is consistent with the "unsmoothed" solution of Chlieh et al. (2014). This is unsurprising since our modeling approach is not affected by any prior-induced spatial smoothing. The high coupling asperity directly offshore of Bahía de Caráquez probably ruptured individually during the 1998 M W = 7.2 earthquake while the coupled Reprint Fig. 3. GPS and tsunami observations used in this study. a) GPS data and model predictions. Black and red arrows show observed and predicted GPS horizontal displacements along with their 95%-confidence ellipses (representing observational and prediction uncertainties, respectively). For the permanent and High-rate GPS, the symbol color represents the vertical displacement. The outer symbol is the observation while the inner symbol is the mean model prediction. b) Observed and predicted tsunami waveforms. The red star defines the event epicenter while black diamonds are the locations of the two DART buoys that recorded the tsunami. For each of them, the amplitude of the first arrival is plotted as a thick black line. The surrounding shaded area marks the 2-σ confidence interval. Stochastic forward model predictions are plotted in red. areas closer to Pedernales could have failed during the 1942 and 2016 earthquake (Fig. 1). On the other hand, the large region of high coupling between Esmeraldas and Cap Manglares could be involved in the 1958 and 1979 ruptures (cf. Fig. 1).
However, we observe larger model uncertainties in this northern part due the lack of offshore measurements (Fig. 2b). This is quite clear in Fig. 2e showing that marginal PDFs close to the trench are nearly uniform. Coupling uncertainties are also illustrated in the supplementary movie M1 showing that large variations in our model ensemble can fit the GPS observations equally well. To quantify the robustness of our coupling map, we calculate the information gain from prior to posterior marginal PDFs using the Kullback-Leibler divergence, defined as: where m C i is the coupling sampled in i-th knot of the triangular mesh. The resulting map shown in Fig. 2c, indicates how much information is gained from the data in different regions of the model. It illustrates the difficulty to infer coupling properties close to the trench using land-based geodetic data. Still, the information gain remains significant within 30-40 km of the coast, and even sometimes almost up to the trench (e.g., offshore of the Manta peninsula and between Esmeraldas and Cap Manglares). This suggests that aforementioned asperities are reliable features of our solution. To evaluate the impact of the mesh size on our solution, we have conducted an inversion using a coarser fault discretization (cf. Fig. S4). Because the heterogeneities visible in Fig. 2 are also present for that coarse parameterization, we choose the finer one to get a better spatial resolution.

Data overview
We use several geodetic datasets covering both near-field and far-field static displacements (cf. Fig. 3a). We gather GPS data from 12 campaign stations and 14 permanent stations with daily solutions (CGPS; Mothes et al., 2018), and 8 high-rate stations (HRGPS; Alvarado et al., 2018). Static offsets from campaign and permanent stations are provided by Nocquet et al. (2017). We estimate our own static displacements from HRGPS by measuring co-seismic offsets from the position before and after the event. We use 1-σ errors provided by Nocquet et al. (2017) for the campaign and CGPS and estimate uncertainties for HRGPS offsets from the standard deviation measured in 20 seconds pre-and post-event time windows. Vertical components of campaign GPS are not used in the inversion as they show large uncertainties. In addition, we use three interferograms derived from ALOS-2 wide-swath descending acquisitions, from ALOS-2 strip-map ascending acquisitions and from Sentinel 1 descending acquisitions (cf. Fig. 4). Unwrapped interferograms are downsampled using a quad-tree algorithm (cf. Fig. S5; Lohman and Simons, 2005). We estimate uncertainties related to atmospheric noise by estimating empirical covariance functions for each interferogram (Jolivet et al., 2012(Jolivet et al., , 2015a. Estimated parameters are summarized in Table S1 and covariance functions are available in Fig. S6. Three nearby DART buoys (Deep-ocean Assessment and Reporting of Tsunamis) recorded the tsunami generated by this event. Unfortunately, the waveform recorded by the closest station (D32067) is unusable for modeling because of multiple data gaps and contamination by seismic waves. We use tsunami waveforms recorded at DART stations D32413 and D32411 (cf. Fig. 3b), as they provide important constraints on the up-dip part of the rupture. To remove tidal signals and reduce high-frequency noise, we band-pass filter the waveforms between 8 min and 3 hours using a third order Butterworth filter. We derive observational uncertainties from standard deviations computed in 140 and 100 min windows before the first arrivals respectively for buoys D32413 and D32411.
We also include near-field seismic waveforms recorded by 10 strong-motion accelerometers and 8 HRGPS stations (cf. Figs. 5 and 6; Alvarado et al., 2018). We integrate the accelerometric data twice and downsample them to 1 sps to match the HRGPS sampling rate. Waveforms are bandpass filtered between 0.015 Hz and 0.08 Hz, except for a few noisy records for which we increased the lower corner frequency to 0.037 Hz (Table S2)

Stochastic co-seismic modeling
Our kinematic modeling of the 2016 Pedernales earthquake is based on a non-planar fault geometry in which the dip varies from 10 • to 27 • between 10 and 50 km depth, following the bending of the Slab1.0 model (cf. Fig. S7; Hayes et al., 2012). The fault   is discretized in 15 × 15 km patches in which we sample static (m S ) and kinematic (m K ) model parameters. The static model vector m S includes two components of static slip in each patch (i.e., the final integrated slip) and extra nuisance parameters to account for InSAR orbital errors (i.e., 3 parameters per interferogram to model a linear function of range and azimuth). The two components of static slip are U , aligned with the direction of convergence between Nazca and NAS plates, and U ⊥ , which is perpendicular to U . The vector of kinematic parameters m K includes rupture velocity and rise time in each patch, along with hypocenter coordinates (i.e., the point of rupture initiation). Each point on the fault is only allowed to rupture once during the earthquake and we prescribe a triangular slip velocity function.
Following the approach of Minson et al. (2013), we first solve the final static slip distribution (i.e., m S ) given available static observations (d S ), i.e., InSAR, GPS offsets and tsunami data. Using AlTar, we thus sample the posterior distribution: where G S is the matrix including Green's functions that are computed using the same layered elastic medium than the one used for the inter-seismic coupling model (cf. section 2). Tsunami waveforms are simulated using COMCOT (Liu et al., 1998) assuming a time step of 1 sec and a 30-arc second GEBCO (General Bathymetric Chart of the Oceans) bathymetry. (Weatherall et al., 2015). As in eq. (1), the misfit covariance C S describes observational errors and prediction uncertainties due to inaccuracies of the assumed elastic structure (Duputel et al., 2012(Duputel et al., , 2014. As we want to promote a dominant thrust motion while allowing local variations of the slip direction, the prior PDF p(m S ) includes uniform prior U (−1 m, 15 m) along the direction of convergence (U ) and Gaussian prior N (0, 0.5 m) in the perpendicular direction (U ⊥ ).
In a second step, we address the full joint inversion problem by incorporating kinematic observations d K . HRGPS and strong motion data provide information on kinematic parameters m K and bring additional constraints on m S . The posterior PDF is then given by Minson et al. (2013): where g K (m S , m K ) is the (non-linear) forward predictions for HRGPS and strong motion waveforms that are based on the Herrmann (2013) implementation of the discrete wave-number method (Bouchon and Aki, 1977). As in eq. (3), C K is the misfit covariance describing measurement errors and predictions uncertainties due to Earth model inaccuracies. The prior p(m k ) is a combination of uniform priors U (1 s, 12 s) and U (1 km/s, 4 km/s) for rise-time and rupture velocity and a Gaussian PDF N (x h , σ = 5 km) for the hypocenter coordinates (x h ).

Co-seismic modeling results
The Pedernales rupture is mainly unidirectional with a significant southward directivity (see posterior mean model in Fig. 7, cumulative slip snapshots in Fig. 8, and supplementary movie M2; Ye et al., 2016;Nocquet et al., 2017;Yi et al., 2018). The inverted hypocenter (0.31 • N, −80.15 • W, depth = 19.6 km; indicated by the red star in Fig. 7), is consistent with estimates from the Instituto Geofísico de la Escuela Politécnica Nacional (0.35 • N, −80.16 • W, depth = 17.0 km; http://www.igepn .edu .ec). Our solution depicts two large slip asperities separated by ∼50 km that coincides roughly with two high-coupling zones north and south of the equator in Fig. 1 and supplementary movie M1. The first asperity is located close to the epicenter and fails within 15 s after the origin time (Fig. 8). The second slip asperity ruptures about 10 s later and contributes to more than 60% of the total seismic moment. The rupture directivity and the location of the southernmost asperity, with slip up to 8 m below the coastline, can probably ex- plain the large damages that have been reported south of the city of Pedernales .
Posterior model uncertainties indicate that we have good constraints on slip amplitude through the fault plane ( Fig. 7b and supplementary movie M3). Moreover, stochastic rupture fronts presented in Fig. 7a show that rupture initiation times are well resolved in large slip areas. There is however a tradeoff between rupture initiation times and rise times as illustrated in Fig. 7c-d. This is because our seismic observations are mostly sensitive on subfault centroid times rather than on rupture times and rise times, resulting in a negative correlation with a −1 slope between the two later parameters.
The southward directivity is clearly visible on HRGPS and strong motion data that show large ground motion amplitudes south of the rupture. This is well captured by our stochastic model predictions (Figs. 5 and 6). Some discrepancies are visible in the late arrivals, which are probably due to unaccounted 3D heterogeneities. Geodetic measurements provide good constraints on the static slip pattern, with large static displacements observed above the large slip asperity in the south. Our solution is able to predict GPS measurements (Fig. 3a) and InSAR data, with small residuals for Sentinel and ALOS-2 data (Fig. 4). We notice larger misfits for the ALOS-2 descending track, probably due to atmospheric noise since this interferogram is associated with significant spatially-correlated observational noise (cf. Fig. S6). Our solution also provides satisfactory fit to tsunami waveforms despite their relatively small amplitude (<1 cm, Fig. 3b). These tsunami observations are important since they clearly show the absence of slip in the shallow portion of the fault (shallow slip would produce large amplitude waves arriving too early at DART stations). This is also reported by Ye et al. (2016) that conducted trial and error teleseismic inversions, progressively removing shallow rows of patches to match the onset of tsunami signals.

Strain budget along the Colombia-Ecuador subduction zone
The Colombia-Ecuador subduction zone provides an outstanding opportunity to study the behavior of a megathrust fault over multiple earthquake cycles. As mentioned above, before the 2016 Pedernales earthquake, the subduction interface experienced a sequence of megathrust ruptures that started with a large M W = 8.6 event in 1906 followed by a series of smaller earthquakes in 1942, 1958, 1979 and 1998. Because these events seem to cluster in time, it has been suggested that strain released by most recent earthquakes exceeds the deformation that accumulated inter-seismically since 1906 Yi et al., 2018).
The strain budget along the megathrust can be investigated by comparing the co-seismic moment generated by earthquakes with the moment deficit accumulated during previous inter-seismic periods. We define the moment deficit accumulated over an interseismic time-span T over an area A of a fault as: where V p is the long-term convergence rate, μ(x) is the shear modulus along the subduction interface and m I (x) is the coupling model introduced in section 2. Using such approach,  propose that the co-seismic moment of the 1942 and 2016 earthquakes are much larger than the deficit accumulated since the 1906 earthquake (by a factor of 3 to 5 times for the 1942 event and 1.3 to 1.6 times for 2016). This seems also true for northern segments and in particular for the 1958 earthquake that has a seismic moment 1.5 to 1.8 times larger than the moment deficit estimated from the modeling of geodetic coupling. As discussed by Nocquet et al. (2017) and Yi et al. (2018), these estimates remain questionable given uncertainties on co-seismic slip and inter-seismic coupling. Hereafter, we use our stochastic co-seismic and inter-seismic solutions to fully account for posterior uncertainties and address the strain budget probabilistically. We assume a magnitude of M W = 7.8 ±0.2 for the 1942 earthquake (Swenson and Beck, 1996;Ye et al., 2016). We compare the probability distributions of seismic moment generated by the 1942 and 2016 earthquakes with the moment deficit accumulated since 1906 (Fig. 9). Assuming the two events are co-located, maximum a posteriori models indicate that the seismic moment for the 1942 and 2016 events are larger than the accumulated deficit by a factor of 2.0 and 1.2, respectively. Taken together for the 1906-2016 period, the moment generated co-seismically is 1.3 times larger than the moment deficit accumulated inter-seismically. Those estimates are subject to considerable uncertainties reflected by the overlap between the PDFs (Fig. 9). Although this overlap is not negligible, there is a relatively small probability of about 5% to have a moment deficit larger or equal than the cumulative seismic moment of the 1942 and 2016 earthquakes. In this scenario, an excess of co-seismic moment since 1906 is likely given available observations. This conclusion only holds if the 2016 rupture largely overlaps with the 1942 earthquake, whose location is still debated. In particular, Yi et al. (2018) suggests that the 1942 earthquake occurred at shallower depth than the 2016 rupture from the comparison of macroisoseismic maps of 1942 and 1958 events (Swenson and Beck, 1996). Therefore we test the alternative hypothesis suggesting that the 1942 earthquake occurred between lat 0.5 • S-0.5 • N at a depth shallower than 40 km . Fig. 10a, b shows that the negative moment balance no longer holds. In this case, the probability of having a deficit equivalent or larger than the co-seismic moment is larger than 70%. Fig. 10c, d shows that this remains true if we further restrain the location of the 1942 event to be located updip of the 2016 earthquake (as proposed by Yi et al., 2018).
We conducted a similar analysis for the 1958 northern Ecuador earthquake, assuming a magnitude M W = 7.6 ± 0.2 (according to Ye et al., 2016). Maximum a posteriori models in Fig. S9b show that the seismic moment generated by the 1958 earthquake is quite similar to the accumulated deficit between 1906 and 1958. This contradicts with Nocquet et al. (2017) that estimated that the 1958 earthquake had a seismic moment exceeding by 50% to 180% the moment accumulated inter-seismically. In our case, we clearly see that the PDF of the moment deficit falls within uncertainties of the 1958 co-seismic moment. As shown in Fig. S9, this still holds if we assume different location for the 1958 earthquake, which discards the negative moment balance issue reported for 1942 and 2016 earthquakes.

Discussion and conclusion
We develop stochastic models of the inter-seismic slip-rate along the Colombia-Ecuador subduction and of the 2016 Pedernales earthquake, which provide new constraints on uncertainties of inter-and co-seismic slip processes. Our results are to first order consistent with some previously published models (e.g., Chlieh et al., 2014). In particular, our coupling model presented in Fig. 2 is similar to the "unsmoothed" model of  since it is not affected by smoothing regularization. Our solution clearly depicts a heterogeneous coupling of the subduction interface (cf. Fig. 2). The heterogeneity of fault coupling properties seems to be a common feature to many subduction zones (Avouac, 2015), but is often blurred because of poor spatial resolution and smoothing constraints used in the inversion. Despite large uncertainties due to the lack of geodetic observations far offshore, all Such heterogeneity roughly correlates with the spatial complexity of the 2016 earthquake revealed by our co-seismic solution.
Our results indicate a unidirectional rupture towards the South with two large slip zones that coincide with two high-coupling asperities in the inter-seismic solution (cf. Fig. 1). We evaluate the possibility that the seismic moment generated by the 1942 and 2016 earthquakes is larger than the moment deficit accumulated since the great 1906 earthquake (as suggested by Nocquet et al., 2017). Our analyses show that this conclusion only holds if we assume that there is a large overlap between the 1942 and 2016 ruptures. If this particular assumption is loosened, results indicate that such an unbalanced moment budget is no longer required by observations. North of the Pedernales rupture, we also show that the seismic moment of the 1958 earthquake is not necessarily larger than the deficit accumulated since 1906 given uncertainties in co-and inter-seismic processes. The question therefore entirely lies within the accuracy of the location and extent of historical earthquakes.
One of the previously mentioned argument favoring an overlap between 1942 and 2016 earthquakes comes from the analysis of teleseismic waveforms recorded at a similar location for both events (see details in supplementary text T1). Ye et al. (2016) showed that 1942 and 2016 waveforms at the station DBN (De Bilt, Netherlands) present significant dissimilarities. Fig. S10 shows that such discrepancies can be explained by differences in the hypocenter location with the same slip distribution for both events (as previously suggested by Nocquet et al., 2017). However, the shape of the observed teleseismic P-wave is mostly controlled by the relative location between the hypocenter and the main slip asperities (i.e., by the corresponding apparent moment-rate function). In fact, Fig. S10c shows that the 1942 DBN waveform could be explained equally well if we assume that the 1942 rupture occurred updip of the 2016 earthquake as suggested by focal depth and macroisoseismic maps of the 1942 event (Yi et al., 2018). In this scenario, there is a probability of ∼70% to have a balanced moment budget since 1906 (i.e., a moment deficit that is larger or equal to the seismic moment of 1942 and 2016 events). On the contrary, if there is a large overlap between both earthquakes, our results show that there is a 95% probability that the moment generated by 1942 and 2016 ruptures is larger than the moment deficit accumulated since 1906. In this case such an unbalanced moment budget can possibly be explained by temporal variations in strain accumulation, which have been observed for example before and after the 2011 M W = 9.0 Tohoku earthquake (e.g., Mavrommatis et al., 2014;Heki and Mitsui, 2013) and after the 2010 Maule earthquake (Melnick et al., 2017;Loveless, 2017). Alternatively, Nocquet et al. (2017) propose a "supercycle" model where the apparent excess of co-seismic moment results from the fact that the 1906 and 1942 earthquakes did not release all of the accumulated strain along the megathrust. This is consistent with the modeling of historic tsunami records suggesting that the 1906 earthquake mainly ruptured the shallow part of the subduction without involving much slip close to the 2016 Pedernales event (Yoshimoto et al., 2017). However, these estimates might be biased by the poor sensitivity of tsunami data to deep slip, which can explain the relatively low magnitude of their resulting model (M W = 8.4). The fact that the surface wave magnitude M s = 8.6 is otherwise consistent with M W also suggests that the 1906 earthquake is not a typical "tsunami" earthquake and is therefore probably not associated with a predominantly shallow rupture (Kanamori, 1972).
The complex behavior of the Colombia-Ecuador subduction can be related to the large heterogeneity revealed by our coupling solution, which suggests significant spatial variability of fault frictional properties (Fig. 1) Such frictional heterogeneities could result from spatial variations in rheology, fluid pore pressure (e.g., Avouac, 2015) or fault roughness associated with the subduction of topographic features such as ridges, fracture zones, and seamounts (Collot et al., 2017;Graindorge, 2004). As shown for example by Kaneko et al. (2010), such frictional heterogeneity can produce earthquakes of different sizes re-rupturing the same fault region at short time intervals. Complex earthquake sequences may also be promoted by partial stress drop of past events that produces significant stress heterogeneity along the fault (Cochard and Madariaga, 1996). The fact that large earthquakes (like the 1906 event) are rapidly followed by sequences of smaller ruptures (e.g., in 1942, 1958, 1979, 1998 and 2016) can then be understood if static stress drop of the smaller events is small compared to the increase of dynamic stresses at rupture fronts (Heaton, 1990;Melgar and Hayes, 2017).
As instrumental observations accumulate, there is a growing record of large earthquakes that break portions of faults that experienced previously documented large ruptures. These earthquakes continuously provide new observations suggesting complex earthquake sequences with substantial spatial and temporal variability among successive ruptures of the same fault system. As shown here, the study of long-term earthquake sequences and the associated strain budget still relies on many assumptions and are affected by large uncertainties. To address the seismogenic behavior of active faults, we need to quantify how large observational and modeling uncertainties are and how much information we have gained in comparison to our preconceptions. Inaccuracies on historical earthquakes size and position can be substantial and also need to be properly considered. Such quantitative analysis is essential to understand how strain accumulates inter-seismically and is released by earthquakes, thereby improving seismic hazard assessment along subduction zones.

Acknowledgement
The ALOS-2 original data are copyright JAXA and provided under JAXA RA4 PI Project P1372002. The Copernicus Sentinel-1 data were provided by the European Space Agency (ESA). Contains modified Copernicus data 2016, processed by ESA and NASA/JPL. We thank the Instituto Geofísico de la Escuela Politécnica Nacional (IG-EPN Ecuador) and the IRD for making available the static GPS, high-rate GPS and accelerometers data used in this study. Aside from the data from the IG-EPN and IRD networks, the data set also includes accelerometer data from OCP and GPS data from the Instituto Geográfico Militar (IGM) and Servicio Geológico Colombiano, Proyecto GeoRED from Colombia, that we also thank. The waveform of the 1942 earthquake used on this study was provided by Lingling Ye and Hiroo Kanamori. This project has received funding from the Agence Nationale de la Recherche (ANR-17-ERC3-0010), the European Union's Horizon 2020 research and innovation program (grant agreement 758210), and the CNRS international program for scientific co-operation (PICS). This research was also supported by the NASA Earth Surface and Interior focus area and performed at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration. We thank the Editor, Jean-Philippe Avouac, and two anonymous reviewers for their constructive comments, which helped improve this manuscript.