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

B. Gombert, Z. Duputel, R. Jolivet, M. Simons, J. Jiang, C. Liang, E.J. Fielding, L. Rivera

► To cite this version:

B. Gombert, Z. Duputel, R. Jolivet, M. Simons, J. Jiang, et al. Strain budget of the Ecuador–Colombia subduction zone: A stochastic view. Earth and Planetary Science Letters, 2018, 498, pp.288–299. 10.1016/j.epsl.2018.06.046. hal-01920268

HAL Id: hal-01920268

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

Submitted on 20 Aug 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
Strain budget of the Ecuador-Colombia subduction zone: a
stochastic view

B. Gombert<sup>a,*</sup>, Z. Duputel<sup>a</sup>, R. Jolivet<sup>b</sup>, M. Simons<sup>c</sup>, J. Jiang<sup>d</sup>, C. Liang<sup>e</sup>, E. Fielding<sup>e</sup>, L. Rivera<sup>a</sup>

<sup>a</sup>Institut de Physique du Globe de Strasbourg, UMR7516, Université de Strasbourg, EOST/CNRS

<sup>b</sup>Laboratoire de géologie, Département de Géosciences, École Normale Supérieure, PSL Research University, CNRS UMR 8538, Paris, France

<sup>c</sup>Seismological Laboratory, Geological and Planetary Sciences, California Institute of Technology, Pasadena, California, USA

<sup>d</sup>Institute of Geophysics and Planetary Physics, University of California, San Diego, La Jolla, California, USA

<sup>e</sup>Jet Propulsion Laboratory, California Institute of Technology, Pasadena, California, USA

Abstract

The 2016 Pedernales earthquake (M<sub>W</sub>=7.8) ruptured a portion of the Colombia-Ecuador subduction interface where several large historical earthquakes have been documented since the great 1906 earthquake (M=8.6). Considering all significant ruptures that occurred in the region, it has been suggested that the cumulative moment generated co-seismically along this part of the subduction over the last century exceeds the moment deficit accumulated inter-seismically since 1906. Such an excess challenges simple models with earthquakes resetting the elastic strain accumulated inter-seismically in locked asperities. These inferences are however associated with large uncertainties that are generally unknown. The impact of spatial smoothing constraints on co-seismic and inter-seismic models also prevents any robust assessment of the strain budget. We propose a Bayesian kinematic slip model of the 2016 Pedernales earthquake using the most comprehensive dataset to date including InSAR and

<footnote>*Corresponding author</footnote>

<footnote>Email address: gombert@unistra.fr (B. Gombert)</footnote>
GPS offsets, tsunami waveforms, and kinematic records from high-rate GPS and strong-
motions. In addition, we use inter-seismic geodetic velocities to produce a probabilistic
inter-seismic coupling model of the subduction interface. Our stochastic co-seismic and
inter-seismic solutions include the ensemble of all plausible models consistent with our prior
information and that fit the observations within uncertainties. The analysis of these model
ensembles indicates that an excess of co-seismic moment during the 1906 - 2016 period is
likely in Central Ecuador only if we assume that 1942 and 2016 earthquakes are colocated. If
this assumption is relaxed, we show that this conclusion no longer holds given uncertainties
in co- and inter-seismic processes. The comparison of 1942 and 2016 teleseismic records
reveals large uncertainties in the location of the 1942 event, hampering our ability to draw
strong conclusions on the unbalanced moment budget in the region. Our results also show a
heterogeneous coupling of the subduction interface that coincides with two slip asperities in
our co-seismic model for the 2016 Pedernales earthquake and with the location of historical
ruptures in 1958, 1979 and 1998. The spatial variability in coupling and complexity in
earthquake history suggest strong heterogeneities in frictional properties of the subduction
megathrust.

Keywords: Ecuador-Colombia subduction zone, Strain budget, Bayesian inversion,
Kinematic source model, Geodetic coupling model

1. Introduction

2 A long standing question is the existence of persistent fault segments remaining locked

3 in the inter-seismic period and failing suddenly during earthquakes while the surround-
ing interface creeps continuously. This conceptual model predicts so-called “characteristic”
earthquakes repeatedly rupturing the same locked fault segments with either periodic, time-
predictable or slip-predictable behaviours (Shimazaki and Nakata, 1980; Schwartz and Cop-
persmith, 1984). This paradigm is contradicted by an increasing number of observations
showing that the same fault area can break entirely in a single large earthquake ($M_W>8.5$)
but also in a series of smaller ruptures. A remarkable example of such behaviour is the
Colombia-Ecuador subduction zone that experienced a complex sequence of earthquakes
since the beginning of the 20th century (see Figure 1). In 1906, the great $M_W=8.6$ earth-
quake ruptured a ~500-km-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 \le 8.2$ events in 1942, 1958, 1979 and 1998 (Kanamori and McNally, 1982; Beck
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 vari-
ability among successive ruptures is also observed in other regions like Japan and Sumatra
where recent $M_W \sim 9$ megathrust earthquakes ruptured large fault segments that previously
experienced a serie 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
25 more fault slip than the deficit accumulated since the 1942 earthquake (Nocquet et al., 2017).

26 Similar observations are reported in other regions, for example in 1797 and 1833 earthquakes
27 in Sumatra (Sieh et al., 2008), 1812 and 1857 earthquakes in California (Jacoby et al., 1988;

28 Heaton, 1990), and for the 2003 M<sub>W</sub>=7.6 and 2013 M<sub>W</sub>=7.8 Scotia sea earthquakes (Vallée
29 and Satriano, 2014). Such spatial and temporal clustering can be caused by spatial variations
30 of fault coupling associated with heterogeneous frictional properties (Kaneko et al., 2010).

31 Moreover, there can be fluctuations in the patterns of inter-seismic fault coupling before large
32 earthquakes (Perfettini and Avouac, 2004; Mavrommatis et al., 2014; Yokota and Koketsu,
33 2015) or during the post-seismic response of nearby large earthquakes (Heki and Mitsui,
34 2013; Melnick et al., 2017).

35 Although the existence of an anomalously large co-seismic slip associated with a su-
36 percycle behaviour is plausible, other studies suggest that the seismic moment of the 2016

37 Pedernales earthquake is actually consistent with the strain accumulated in the region since
38 the 1942 and 1906 earthquakes (e.g., Ye et al., 2016; Yoshimoto et al., 2017; Yi et al., 2018).

39 These contrasting statements partly results from the ill-posed nature of inter- and co-seismic
40 slip inversions used to evaluate the strain budget along the megathrust. Such inferences
41 are affected by the lack of resolution near the trench during the inter-seismic period but
42 also by non-physics-based smoothing constraints used to regularize slip inversions. In ad-
43 dition, inter- and co-seismic estimates usually do not incorporate rigorous uncertainties (or
44 very often, no uncertainty at all), which complicates a quantitative assessment of the overall
45 strain budget. Strain budget analyses also suffer from the lack of information about past
46

47

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.

48

49

50

51

52

53

54

55

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<sub>W</sub>=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).

56

2. Geodetic coupling

57

2.1. *Stochastic inter-seismic modeling*

58

59

60

61

62

63

64

65

66

We first compute a stochastic model of geodetic coupling along the Ecuadorian subduc-
tion 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 (Mothes et al., 2018; Mora-Páez et al., 2018). Considering that more
than 20 years separate the 1979 earthquake and the first GPS measures, 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 (c.f.,
Fig. S1 in the electronic supplements). Using a back-slip approach (Savage, 1983), we in-
vert 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 ($\mathbf{m}_I$) that fit the GPS data ($\mathbf{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(\mathbf{m}_I) = U(-0.05 \cdot V_p, 1.05 \cdot 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

$$p(\mathbf{m}_I|\mathbf{d}_I) \propto p(\mathbf{m}_I) \exp \left[ -\frac{1}{2} (\mathbf{d}_I - \mathbf{G}_I \mathbf{m}_I)^T \mathbf{C}_I^{-1} (\mathbf{d}_I - \mathbf{G}_I \mathbf{m}_I) \right] \quad (1)$$

where $\mathbf{G}_I$ is the Green's function matrix and $\mathbf{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. (2012, 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(\mathbf{m}_I|\mathbf{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 ($\mathbf{m}_C$) using $\mathbf{m}_C = 1 - \mathbf{m}_I/V_p$.

## 2.2. 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 Fig. 2a along with the associated 2-$\sigma$ posterior uncertainties in Fig. 2b. The posterior median model available in Fig. S3 is very similar to the posterior mean, confirming that most marginal PDFs are nearly Gaussians. The variability of the model population composing the solution is shown in supplementary movie M1.

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 high-coupling 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 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:

$$D_{KL}^i = \int p(m_c^i|\mathbf{d}_c) \log_2 \frac{p(m_c^i)}{p(m_c^i|\mathbf{d}_c)} dm_c^i \quad (2)$$
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 parameterisation, we chose the finer one to get a better spatial resolution and avoid any bias due to a coarse fault discretization.

## 3. Rupture process of the 2016 Pedernales earthquake

### 3.1. 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-\sigma$ 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 win-
dows. 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 de-
scending 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 interfer-
ogram (Jolivet et al., 2012, 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 contamina-
tion 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 (c.f. 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). Wave-forms are inverted in a 150 s-long time window starting from the origin time of the mainshock (23:58:36 UTC).

## 3.2. 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 ($\mathbf{m}_S$) and kinematic ($\mathbf{m}_K$) model parameters. The static model vector $\mathbf{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_\parallel$, aligned with the direction of convergence between Nazca and NAS plates, and $U_\perp$, which is perpendicular to $U_\parallel$. The vector of kinematic parameters $\mathbf{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., $\mathbf{m}_S$) given available static observations ($\mathbf{d}_S$), i.e., InSAR, GPS offsets and tsunami data. Using AlTar, we thus sample the posterior distribution:

$$ 
\begin{aligned}
p(\mathbf{m}_S|\mathbf{d}_S) & \propto p(\mathbf{m}_S) p(\mathbf{d}_S|\mathbf{m}_S) \\
& \propto p(\mathbf{m}_S) \exp \left[ -\frac{1}{2}(\mathbf{d}_S - \mathbf{G}_S\mathbf{m}_S)^T \mathbf{C}_S^{-1} (\mathbf{d}_S - \mathbf{G}_S\mathbf{m}_S) \right]
\end{aligned} 
\quad (3) 
$$
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, 2014). As we want to promote a dominant thrust motion while allowing local variations of the slip direction, the prior PDF $p(\mathbf{m}_S)$ includes uniform prior $\mathcal{U}(-1 \text{ m}, 15 \text{ m})$ along the direction of convergence ($U_\parallel$) and Gaussian prior $\mathcal{N}(0, 0.5 \text{ m})$ in the perpendicular direction ($U_\perp$).

In a second step, we address the full joint inversion problem by incorporating kinematic observations $\mathbf{d}_K$. HRGPS and strong motion data provide information on kinematic parameters $\mathbf{m}_K$ and bring additional constraints on $\mathbf{m}_S$. The posterior PDF is then given by (Minson et al., 2013):

$$p(\mathbf{m}_S, \mathbf{m}_K | \mathbf{d}_S, \mathbf{d}_K) \propto p(\mathbf{m}_K) p(\mathbf{m}_S | \mathbf{d}_S) p(\mathbf{d}_K | \mathbf{m}_S, \mathbf{m}_K) \quad (4)$$

$$\propto p(\mathbf{m}_K) p(\mathbf{m}_S | \mathbf{d}_S) \exp \left[ -\frac{1}{2} (\mathbf{d}_K - \mathbf{g}_K(\mathbf{m}_S, \mathbf{m}_K))^T \mathbf{C}_K^{-1} (\mathbf{d}_K - \mathbf{g}_K(\mathbf{m}_S, \mathbf{m}_K)) \right]$$

where $\mathbf{g}_K(\mathbf{m}_S, \mathbf{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), $\mathbf{C}_K$ is the misfit covariance describing measurement errors and predictions uncertainties due to Earth model inaccuracies.

The prior $p(\mathbf{m}_K)$ is a combination of uniform priors $\mathcal{U}(1 \text{ s}, 12 \text{ s})$ and $\mathcal{U}(1 \text{ km/s}, 4 \text{ km/s})$ for
204 rise-time and rupture velocity and a Gaussian PDF $N(\mathbf{x}_h, \sigma = 5 \text{ km})$ for the hypocenter
205 coordinates $(\mathbf{x}_h)$.

### 3.3. Co-seismic modeling results

207 The Pedernales rupture is mainly unidirectional with a significant southward directivity
208 (see posterior mean model in Fig. 7, cumulative slip snapshots in Fig. 8, and supplementary
209 movie M2; Ye et al., 2016; Nocquet et al., 2017; Yi et al., 2018). The inverted hypocenter
210 (0.31° N, -80.15° W, depth=19.6 km; indicated by the red star in Fig. 7), is consistent with
211 estimates from the Instituto Geofísico de la Escuela Politécnica Nacional (0.35° N, -80.16° W,
212 depth=17.0 km; http://www.igepn.edu.ec). Our solution depicts two large slip asperities
213 separated by ~50 km that coincides roughly with two high-coupling zones north and south
214 of the equator in Fig. 1 and supplementary movie M1. The first asperity is located close to
215 the epicenter and fails within 15 s after the origin time (Fig. 8). The second slip asperity
216 ruptures about 10 s later and contributes to more than 60% of the total seismic moment.
217 The rupture directivity and the location of the southernmost asperity, with slip up to 8 m
218 below the coastline, can probably explain the large damages that have been reported south
219 of the city of Pedernales (Nocquet et al., 2017).

220 Posterior model uncertainties indicate that we have good constraints on slip amplitude
221 through the fault plane (Fig. 7b and supplementary movie M3). Moreover, stochastic rupture
222 fronts presented in Fig. 7a show that rupture initiation times are well resolved in large
223 slip areas. There is however a tradeoff between rupture initiation times and rise times as
224 illustrated in Fig. 7c-d. This is because our seismic observations are mostly sensitive on
225 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.
226

227 The southward directivity is clearly visible on HRGPS and strong motion data that show
228 large ground motion amplitudes south of the rupture. This is well captured by our stochastic
229 model predictions (Figs. 5 and 6). Some discrepancies are visible in the late arrivals, which
230 are probably due to unaccounted 3D heterogeneities. Geodetic measurements provide good
231 constraints on the static slip pattern, with large static displacements observed above the
232 large slip asperity in the south. Our solution is able to predict GPS measurements (Fig. 3a)
233 and InSAR data, with small residuals for Sentinel and ALOS-2 data (Fig. 4). We notice
234 larger misfits for the ALOS-2 descending track, probably due to atmospheric noise since
235 this interferogram is associated with significant spatially-correlated observational noise (cf.,
236 Fig. S6). Our solution also provides satisfactory fit to tsunami waveforms despite their
237 relatively small amplitude (<1 cm, Fig. 3b). These tsunami observations are important since
238 they clearly show the absence of slip in the shallow portion of the fault (shallow slip would
239 produce large amplitude waves arriving too early at DART stations). This is also reported by
240 Ye et al. (2016) that conducted trial and error teleseismic inversions, progressively removing
241 shallow rows of patches to match the onset of tsunami signals.
242

4. Strain budget along the Colombia-Ecuador subduction zone
---

243 The Colombia-Ecuador subduction zone provides an outstanding opportunity to study
244 the behaviour of a megathrust fault over multiple earthquake cycles. As mentioned above,
245 before the 2016 Pedernales earthquake, the subduction interface experienced a sequence of
246 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 (Nocquet et al., 2017; 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 inter-seismic time-span $T$ over an area $A$ of a fault as:

$$M_0^{\text{deficit}} = T V_p \iint_A \mu(\mathbf{x}) m_c(\mathbf{x}) d\mathbf{x} \qquad (5)$$

where $V_p$ is the long-term convergence rate, $\mu(\mathbf{x})$ is the shear modulus along the subduction interface and $m_I(\mathbf{x})$ is the coupling model introduced in section 2. Using such approach, Nocquet et al. (2017) 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 \pm 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 (Nocquet et al., 2017). 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\pm0.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.

## 5. 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., Nocquet et al., 2017; Chlieh et al., 2014). In particular, our coupling model presented in Fig. 2 is similar to the "unsmoothed" model of Chlieh et al. (2014) 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 models in our posterior ensemble show a large spatial heterogeneity (cf., supplementary movie M1).

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 favouring 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 macro-isoseismic 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<sub>W</sub>=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<sub>W</sub>=8.4). The fact that the surface wave magnitude M<sub>s</sub>=8.6 is otherwise consistent with M<sub>W</sub> 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 behaviour 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 behaviour of active faults, we need to quantify how large observational and modelling 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.

## 6. Acknowledgment

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
392 area and performed at the Jet Propulsion Laboratory, California Institute of Technology. We

393 thank the Editor, Jean-Philippe Avouac, and two anonymous reviewers for their constructive

394 comments, which helped improve this manuscript.
Alvarado, A., Ruiz, M., Mothes, P., Yepes, H., Segovia, M., Vaca, M., Ramos, C., Enríquez, W., Ponce, G., Jarrín, P., et al. (2018). Seismic, volcanic, and geodetic networks in Ecuador: Building capacity for monitoring and research. *Seismological Research Letters*, 395 396
397
398 89(2A):432–439.

Avouac, J.-P. (2015). From geodetic imaging of seismic and aseismic fault slip to dynamic modeling of the seismic cycle. *Annual Review of Earth and Planetary Sciences*, 400 401
402
403 43(1):233271.

Beck, S. L. and Ruff, L. J. (1984). The Rupture Process of the Great 1979 Colombia Earthquake - Evidence for the Asperity Model. *J. Geophys. Res.*, 89(NB11):9281–9291.

Béthoux, N., Segovia, M., Alvarez, V., Collot, J.-Y., Charvis, P., Gailler, A., and Monfret, T. (2011). Seismological study of the central Ecuadorian margin: Evidence of upper plate deformation. *J. South Amer. Earth Sci.*, 31(1):139–152.

Bouchon, M. and Aki, K. (1977). Discrete wave-number representation of seismic-source wave fields. *Bulletin of the Seismological Society of America*, 67(2):259–277.

Chlieh, M., Mothes, P., Nocquet, J.-M., Jarrin, P., Charvis, P., Cisneros, D., Font, Y., Collot, J.-Y., Villegas-Lanza, J.-C., Rolandone, F., et al. (2014). Distribution of discrete seismic asperities and aseismic slip along the Ecuadorian megathrust. *Earth and Planetary Science Letters*, 400:292–301.

Cochard, A. and Madariaga, R. (1996). Complexity of seismicity due to highly
rate-dependent friction. *Journal of Geophysical Research: Solid Earth* (1978-2012), 101(B11):25321-25336.

Collot, J.-Y., Sanclemente, E., Nocquet, J.-M., Leprtre, A., Ribodetti, A., Jarrin, P., Chlieh, M., Graindorge, D., and Charvis, P. (2017). Subducted oceanic relief locks the shallow megathrust in central ecuador. *Journal of Geophysical Research: Solid Earth*, page 2016JB013849.

Duputel, Z., Agram, P. S., Simons, M., Minson, S. E., and Beck, J. L. (2014). Accounting for prediction uncertainty when inferring subsurface fault slip. *Geophys. J. Int.*, 197(1):464-482.

Duputel, Z., Rivera, L., Fukahata, Y., and Kanamori, H. (2012). Uncertainty estimations for seismic source inversions. *Geophys. J. Int.*, 190(2):1243-1256.

Graindorge, D. (2004). Deep structures of the Ecuador convergent margin and the Carnegie Ridge, possible consequence on great earthquakes recurrence interval. *Geophys. Res. Lett.*, 31(4):L04603.

Gutenberg, B. and Richter, C. F. (1949). *Seismicity of the Earth and Associated Phenomena*. Princeton University Press, Princeton, New Jersey.

Hayes, G. P., Wald, D. J., and Johnson, R. L. (2012). Slab1.0: A three-dimensional model of global subduction zone geometries. *J. Geophys. Res.: Solid Earth*, 117(B1).

He, P., Hetland, E. A., Wang, Q., Ding, K., Wen, Y., and Zou, R. (2017). Coseismic Slip
433 in the 2016 Mw 7.8 Ecuador Earthquake Imaged from Sentinel-1A Radar Interferometry.

434 *Seismol. Res. Lett.*, 88(2A):277-286.

435 Heaton, T. H. (1990). Evidence for and implications of self-healing pulses of slip in earthquake rupture. *Physics of the Earth and Planetary Interiors*, 64(1):1-20.

437 Heki, K. and Mitsui, Y. (2013). Accelerated pacific plate subduction following interplate thrust earthquakes at the Japan trench. *Earth Planet. Sci. Lett.*, 363:44-49.

439 Herrmann, R. B. (2013). Computer Programs in Seismology: An Evolving Tool for Instruction and Research. *Seismol. Res. Lett.*, 84(6):1081-1088.

441 Jacoby, G. C., Sheppard, P. R., and Sieh, K. E. (1988). Irregular Recurrence of Large Earthquakes Along the San Andreas Fault: Evidence from Trees. *Science*, 241(4862):196-199.

444 Jolivet, R., Candela, T., Lasserre, C., Renard, F., Klinger, Y., and Doin, M.-P. (2015a). The Burst-Like Behavior of Aseismic Slip on a Rough Fault: The Creeping Section of the Haiyuan Fault, China. *Bull. Seism. Soc. Am.*, 105(1):480-488.

447 Jolivet, R., Lasserre, C., Doin, M.-P., Guillaso, S., Peltzer, G., Dailu, R., Sun, J., Shen, Z.-K., and Xu, X. (2012). Shallow creep on the Haiyuan fault (Gansu, China) revealed by SAR interferometry. *Journal of Geophysical Research: Solid Earth*, 117(B6).

450 Jolivet, R., Simons, M., Agram, P., Duputel, Z., and Shen, Z.-K. (2015b). Aseismic slip and seismogenic coupling along the central San Andreas Fault. *Geophysical Research Letters*, 42(2):297-306.
453 Kanamori, H. (1972). Mechanism of Tsunami Earthquakes. *Phys. Earth Planet. Inter.*, 454 6:356-359.

455 Kanamori, H. and McNally, K. C. (1982). Variable Rupture Mode of the Subduction Zone Along the Ecuador - Colombia Coast. *Bull. Seism. Soc. Am.*, 72(4):1241-1253.

457 Kaneko, Y., Avouac, J.-P., and Lapusta, N. (2010). Towards inferring earthquake patterns from geodetic observations of interseismic coupling. *Nature Geosci.*, 3(5):ngeo843-369.

459 Klein, E., Duputel, Z., Masson, F., Yavasoglu, H., and Agram, P. (2017). Aseismic slip and seismogenic coupling in the Marmara Sea: What can we learn from onland geodesy? *Geophys. Res. Lett.*, 44(7):3100-3108.

462 Lay, T. (2015). The surge of great earthquakes from 2004 to 2014. *Earth Planet. Sci. Lett.*, 463 409:133-146.

464 Liu, P. L. F., Woo, S.-B., and Cho, Y.-S. (1998). Computer Programs for Tsunami Propagation and Inundation. Cornell University, USA.

466 Lohman, R. B. and Simons, M. (2005). Some thoughts on the use of InSAR data to constrain models of surface deformation: Noise structure and data downsampling. *Geochem. Geophys. Geosyst.*, 6(1):Q01007.

469 Loveless, J. P. (2017). Super-interseismic periods: Redefining earthquake recurrence. *Geo-phys. Res. Lett.*, 44(3):1329-1332.

471 Mavrommatis, A. P., Segall, P., and Johnson, K. M. (2014). A decadal-scale deforma-
472 tion transient prior to the 2011 Mw 9.0 Tohoku-oki earthquake. *Geophys. Res. Lett.*, 41(13):4486-4494.

474 Melgar, D. and Hayes, G. P. (2017). Systematic observations of the slip pulse properties of large earthquake ruptures. *Geophysical Research Letters*, 44(19):9691-9698.

476 Melnick, D., Moreno, M., Quinteros, J., Baez, J. C., Deng, Z., Li, S., and Oncken, O. (2017). The super-interseismic phase of the megathrust earthquake cycle in Chile. *Geophys. Res. Lett.*, 44(2):784-791.

479 Mendoza, C. and Dewey, J. W. (1984). Seismicity Associated with the Great Colombia-Ecuador Earthquakes of 1942, 1958, and 1979 - Implications for Barrier Models of Earth-quake Rupture. *Bull. Seism. Soc. Am.*, 74(2):577-593.

482 Minson, S., Simons, M., and Beck, J. (2013). Bayesian inversion for finite fault earthquake source models ITheory and algorithm. *Geophysical Journal International*, 194(3):1701-1726.

485 Mora-Páez, H., Peláez-Gaviria, J.-R., Diederix, H., Bohórquez-Orozco, O., Cardona-Piedrahita, L., Corchuelo-Cuervo, Y., Ramírez-Cadena, J., and Díaz-Mila, F. (2018). Space geodesy infrastructure in Colombia for geodynamics research. *Seismological Research Letters*, 89(2A):446-451.

489 Mothes, P. A., Rolandone, F., Nocquet, J.-M., Jarrin, P. A., Alvarado, A. P., Ruiz, M. C., Cisneros, D., Páez, H. M., and Segovia, M. (2018). Monitoring the earthquake cycle in
491

the northern Andes from the Ecuadorian cGPS network. *Seismological Research Letters*, 89(2A):534-541.

492

493

Nocquet, J. M., Jarrin, P., Vallée, M., Mothes, P. A., Grandin, R., Rolandone, F., Delouis, B., Yepes, H., Font, Y., Fuentes, D., Regnier, M., Laurendeau, A., Cisneros, D., Hernandez, S., Sladen, A., Singaucho, J. C., Mora, H., Gomez, J., Montes, L., and Charvis, P. (2017). Supercycle at the Ecuadorian subduction zone revealed after the 2016 Pedernales earthquake. *Nature Geoscience*, 10(2):145-149.

494

495

496

497

498

Nocquet, J. M., Villegas-Lanza, J. C., Chlieh, M., Mothes, P. A., Rolandone, F., Jarrin, P., Cisneros, D., Alvarado, A., Audin, L., Bondoux, F., Martin, X., Font, Y., Regnier, M., Vallée, M., Tran, T., Beauval, C., Mendoza, J. M. M., Martinez, W., Tavera, H., and Yepes, H. (2014). Motion of continental slivers and creeping subduction in the northern Andes. *Nature Geosci.*, 7(4):287-291.

499

500

501

502

503

Ortega Culaciati, F. H. (2013). *Aseismic Deformation in Subduction Megathrusts: Central Andes and North-East Japan*. PhD thesis, California Institute of Technology, Pasadena, CA, USA.

504

505

506

Perfettini, H. and Avouac, J. P. (2004). Stress transfer and strain rate variations during the seismic cycle. *J. Geophys. Res.: Solid Earth*, 109(B6).

507

508

Savage, J. (1983). A dislocation model of strain accumulation and release at a subduction zone. *Journal of Geophysical Research: Solid Earth*, 88(B6):4984-4996.

509

510

Schwartz, D. P. and Coppersmith, K. J. (1984). Fault behavior and characteristic earth-
511 quakes: Examples from the Wasatch and San Andreas Fault Zones. *Journal of Geophysical Research: Solid Earth* (1978-2012), 89(B7):5681-5698.

513 Shimazaki, K. and Nakata, T. (1980). Time-predictable recurrence model for large earthquakes. *Geophys. Res. Lett.*, 7(4):279-282.

515 Sieh, K., Natawidjaja, D. H., Meltzner, A. J., Shen, C. C., Cheng, H., Li, K. S., Suwargadi, B. W., Galetzka, J., Philibosian, B., and Edwards, R. L. (2008). Earthquake Supercycles Inferred from Sea-Level Changes Recorded in the Corals of West Sumatra. *Science*, 322(5908):1674-1678.

519 Simons, M., Minson, S. E., Sladen, A., Ortega, F., Jiang, J., Owen, S. E., Meng, L., Ampuero, J.-P., Wei, S., Chu, R., Helmberger, D. V., Kanamori, H., Hetland, E., Moore, A. W., and Webb, F. H. (2011). The 2011 Magnitude 9.0 Tohoku-Oki Earthquake: Mosaicking the Megathrust from Seconds to Centuries. *Science*, 332(6036):1421-1425.

523 Swenson, J. L. and Beck, S. L. (1996). Historical 1942 Ecuador and 1942 Peru subduction earthquakes and earthquake cycles along Colombia-Ecuador and Peru subduction segments. *Pure appl. geophys.*, 146(1):67-101.

526 Vallee, M., Nocquet, J.-M., Battaglia, J., Font, Y., Segovia, M., Regnier, M., Mothes, P., Jarrin, P., Cisneros, D., Vaca, S., et al. (2013). Intense interface seismicity triggered by a shallow slow slip event in the Central Ecuador subduction zone. *Journal of Geophysical Research: Solid Earth*, 118(6):2965-2981.
Vallée, M. and Satriano, C. (2014). Ten year recurrence time between two major earthquakes affecting the same fault segment. *Geophysical Research Letters*, 41(7):2312–2318.

Weatherall, P., Marks, K., Jakobsson, M., Schmitt, T., Tani, S., Arndt, J. E., Rovere, M., Chayes, D., Ferrini, V., and Wigley, R. (2015). A new digital bathymetric model of the world's oceans. *Earth and Space Science*, 2(8):331–345.

Ye, L., Kanamori, H., Avouac, J.-P., Li, L., Cheung, K. F., and Lay, T. (2016). The 16 April 2016, MW7.8 (MS7.5) Ecuador earthquake: A quasi-repeat of the 1942 MS7.5 earthquake and partial re-rupture of the 1906 MS8.6 Colombia–Ecuador earthquake. *Earth Planet. Sci. Lett.*, 454(Supplement C):248–258.

Yi, L., Xu, C., Wen, Y., Zhang, X., and Jiang, G. (2018). Rupture process of the 2016 Mw 7.8 Ecuador earthquake from joint inversion of InSAR data and teleseismic P waveforms. *Tectonophysics*, 722:163–174.

Yokota, Y. and Koketsu, K. (2015). A very long-term transient event preceding the 2011 Tohoku earthquake. *Nature Communications*, 6:5934.

Yoshimoto, M., Kumagai, H., Acero, W., Ponce, G., Vásconez, F., Arrais, S., Ruiz, M., Alvarado, A., Pedraza García, P., Dionicio, V., Chamorro, O., Maeda, Y., and Nakano, M. (2017). Depth-dependent rupture mode along the Ecuador-Colombia subduction zone. *Geophys. Res. Lett.*, 84(5):1561–2210.
<figure><img src="image_2.png" /><figcaption>**Figure 1: Interseismic coupling and historical earthquakes.** The colour scale indicates the geodetic coupling of the subduction interface obtained from inter-seismic GPS velocities (cf. section 2). Blue line and blue star are respectively the 2 m isocontours of co-seismic slip and hypocenter obtained for the 2016 Pedernales earthquake (cf., section 3). Grey dashed lines show the approximate extent of the 1942, 1958, 1979, and 1998 events (Kanamori and McNally, 1982; Chlieh 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<sub>W</sub> = 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).</figcaption></figure>
<figure><img src="image_3.png" /></figure>

Figure 2: Interseismic coupling of the Ecuadorian subduction margin. a) Posterior mean coupling model. Thin black lines represent the fault parametrization. Coupling values are inverted at each triangle knot. Interseismic GPS displacement and model predictions are plotted as black and blue arrows, respectively. b) 2-σ uncertainties of the coupling model. c) Kullback-Leibler divergence between the posterior and prior PDFs of coupling. Higher values indicate regions where the gain of information of the posterior PDF is significant relative to the prior distribution. d) et e) Marginal probability densities for the two nodes pointed out in b) and c).
<figure><img src="image_4.png" /><figcaption>Figure 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 colour 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.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 4: Model performance for InSAR. (a, d, g) InSAR observations. (b, e, h) Predictions for the posterior mean model. (c, f, i) Residuals of the Sentinel (top row), descending ALOS-2 wide-swath (middle row), and ascending ALOS-2 strip-map (bottom row) interferograms. Decimated observations, predictions, and residuals are shown in Fig. S5</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 5:** **High-rate GPS observations and model predictions.** The white diamonds on the top-left map indicate the position of the stations. The red star marks the inverted epicenter location. White rectangles are the fault parametrization. The East, North, and vertical components of each station are plotted around the map. For each waveform, the bold number indicates its maximum amplitude. The station azimuth $\Phi$ and distance $d$ to the epicenter are also given. The black line is the recorded waveform. The gray lines are the stochastic predictions for our posterior model. The red line is the mean of the stochastic predictions.</figcaption></figure>
<figure><img src="image_7.png" /></figure>

Figure 6: **Strong-motion observations and model predictions.** Same as Fig. 5. We show only the stations where three components are available. The remaining 5 waveforms from 3 additional stations and the associated model predictions are shown in Fig. S8
<figure><img src="image_8.png" /><figcaption>Figure 7: **Final co-seismic slip distribution.** a) The colour and arrows on the fault plane indicate the amplitude and direction of slip, respectively. Gray-scale lines are stochastic rupture fronts inferred from our model population plotted at 10s, 20s, and 30s. The darker the lines, the larger the slip at that location. The red star marks the hypocenter location. b) Slip uncertainty. The colour on the fault represents the absolute slip uncertainties. Black contour lines show the co-seismic slip every 1 m, starting from 2 m. c) Marginal probability distribution of rise time and initial rupture time in the first slip asperity (located close to the hypocenter). d) Posterior ensemble of source time functions at the same location of the fault. The source time functions labeled s₁ and s₂ in (d) correspond to rupture initiation times and rise times that are indicated with red stars in (c).</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 8: **Temporal evolution of co-seismic slip.** a) Cumulative slip on the fault 10 s, 15 s, 20 s, 25 s, and 30 s after the origin time. The red colour-scale indicates slip amplitude. The red star marks the epicenter location. b) Evolution of slip rate on the fault. c) Source time function (STF) of the event. Grey lines are stochastic STFs inferred from our model population while the black curve represents the posterior mean STF. Vertical red lines indicate the temporal position of each one of the snapshots</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 9: Comparison of co-seismic moment and moment deficit. a) The background colour represents the coupling posterior mean model. The blue stars show the hypocentre location. Blue lines are the 2m, 3m, and 4m co-seismic slip isocontours. The black dashed line delimits the area where the co-seismic moment and moment deficit are computed. b) Probability densities of the co-seismic moment released by the 1942 earthquake and the moment deficit accumulated between 1906 and 1942 within the dashed ellipse shown in a). c) Probability densities of the co-seismic moment released by the 2016 earthquake and the moment deficit accumulated between 1942 and 2016. d) Probability densities of the co-seismic moment released by the sum of the 1942 and 2016 events, and of the moment deficit during the 1906 - 2016 period.</figcaption></figure>
<figure><img src="image_11.png" /></figure>

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

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

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

Figure 10: Comparison of co-seismic moment and moment deficit considering the 1942 event happened at a different location. a) Same as Figure 9a. The black dashed line delimits the area between 0.5°S and 0.5°N where the co-seismic moment and moment deficit are computed. b) Probability densities of the co-seismic moment and moment deficit in the 1906-2016 period. The co-seismic moment is the sum of the 1942 and 2016 events moment. c) Same as a), but the dashed black area shows where the 1942 earthquake could have been located. d) Probability densities of the co-seismic moment and moment deficit.

The co-seismic moment is the sum of the 1942 and 2016 events moment. The moment deficit is the sum of the moment deficits computed in the updip section (shown in c)) for the 1942 - 2016 period and in the downdip section (ellipse in 9a) for the 1906 - 2016 period.
Supporting information for the main manuscript:

# "Strain budget of the Ecuador-Colombia subduction zone: a stochastic view"

B. Gombert<fnref n="1" />, Z. Duputel<fnref n="1" />, R. Jolivet <fnref n="2" />, M. Simons<fnref n="3" />, J. Jiang<fnref n="4" />, C. Liang<fnref n="5" />, E. Fielding<fnref n="5" />, and L. Rivera<fnref n="1" />

*This electronic supplement is a collection of additional figures referenced in the main article. These figures were added to ensure the precision of the description of our method and results.*

<footnote n="1">¹ Institut de Physique du Globe de Strasbourg, UMR7516, Université de Strasbourg, EOST/CNRS</footnote>

<footnote n="2">² Laboratoire de géologie, Département de Géosciences, École Normale Supérieure,
PSL Research University, CNRS UMR, Paris, France</footnote>

<footnote n="3">³ Seismological Laboratory, Geological and Planetary Sciences, California Institute of Technology,
Pasadena, California, USA</footnote>

<footnote n="4">⁴ Institute of Geophysics and Planetary Physics, University of California, San Diego, La Jolla, California,
USA</footnote>

<footnote n="5">⁵ Jet Propulsion Laboratory, California Institute of Technology, Pasadena, California, USA</footnote>
Supplementary text T1

Following Ye et al. (2014) and Nocquet et al. (2017), we compare waveforms
of the 1942 earthquake recorded at the DBN station (De Bilt, Netherlands) with
stochastic waveform predictions at the same station for the 2016 Pedernales slip
distribution.

We compute displacement Green's functions for each subfault patch using the Kikuchi-Kanamori program (Kikuchi and Kanamori, 2003; Kikuchi, Masayuki and Kanamori, Hiroo, 1982). For comparison, we then convolve predicted stochastic waveforms with the instrumental response of the Galitzin seismometer that recorded the 1942 earthquake (pendulum and galvanometer periods Tp=Tx=25 s and gain factor Vm=310; Charlier and Van Gils, 1953).

In Fig. S10a, we first compare 1942 waveforms with predictions of the kinematic slip model (i.e., for the posterior distributions of slip, rise-times, rupture velocities) and hypocenter location obtained for the 2016 Pedernales earthquake. Model predictions show poor fit to the 1942 earthquake waveform. In Fig. S10b, we then compute predictions for the same kinematic slip distribution, but with a hypocenter location between the two slip asperities. With that hypocenter location, model predictions have a very good fit to the 1942 waveform. Finally, in Fig. S10c, we predict waveforms for a slip distribution on the megathrust interface, but updip of the actual 2016 rupture. Notice, that the dip is different due to the variation of the slab interface geometry with depth. We also correct the slip amplitude for the variation of shear modulus in our velocity model (cf., Fig. S2). Similarly to the previous case, the hypocenter is located between the two slip asperities. In this scenario, we are also able to explain the 1942 waveform. It illustrates that the teleseismic P-waveform is mostly sensitive to the relative location of the hypocenter and slip asperity rather than the absolute location of the earthquake.
**Supplementary movie M1: Variability in the Ecuador-Colombia geodetic coupling solution** The animation is made with 150 models randomly selected in the posterior population represented by the background colour. Grey lines are the 2 m contour intervals of 150 co-seismic models also randomly selected in the posterior population.

**Supplementary movie M2: Temporal evolution of co-seismic slip of the 2016 Pedernales earthquake. (left)** Posterior mean model of the cumulative slip. The bottom-right inset shows the stochastic source time function. **(right)** Incremental slip on the fault. The red star marks the inverted posterior mean hypocenter location.

**Supplementary movie M3: Variability in the 2016 Pedernales earthquake co-seismic slip distribution solution** The animation is made with 200 models randomly selected in the posterior population.

<table><caption>Table S1: InSAR observations used in this study.</caption><thead><tr><th>Satellite</th><th>Orbit</th><th>Acquisition dates</th><th>N° of data</th><th>Std.</th><th>Corr. length</th></tr></thead><tbody><tr><td>ALOS-2</td><td>ascending</td><td>07/02/16 - 01/05/16</td><td>130</td><td>5.3 mm</td><td>2.88 km</td></tr><tr><td>ALOS-2</td><td>descending</td><td>01/04/16 - 29/04/16</td><td>483</td><td>9.2 mm</td><td>11.90 km</td></tr><tr><td>Sentinel-1A</td><td>descending</td><td>12/04/16 - 24/04/16</td><td>380</td><td>5.0 mm</td><td>15.0 km</td></tr></tbody></table>
<table><caption>Table S2: Seismological data and filtering used in this study. We use a 4th order Butterworth bandpass filter.</caption><thead><tr><th rowspan="2">Station</th><th rowspan="2">Type</th><th colspan="3">Filter corner frequencies</th></tr><tr><th>East</th><th>North</th><th>Up</th></tr></thead><tbody><tr><td>bahi</td><td>HRGPS</td><td>0.015Hz - 0.08Hz</td><td>0.015Hz - 0.08Hz</td><td>0.015Hz - 0.08Hz</td></tr><tr><td>cabp</td><td>HRPGS</td><td>0.015Hz - 0.08Hz</td><td>0.015Hz - 0.08Hz</td><td>0.015Hz - 0.08Hz</td></tr><tr><td>ecec</td><td>HRPGS</td><td>0.015Hz - 0.08Hz</td><td>0.015Hz - 0.08Hz</td><td>0.015Hz - 0.08Hz</td></tr><tr><td>flfr</td><td>HRPGS</td><td>0.015Hz - 0.08Hz</td><td>0.015Hz - 0.08Hz</td><td>0.015Hz - 0.08Hz</td></tr><tr><td>mlec</td><td>HRPGS</td><td>0.015Hz - 0.08Hz</td><td>0.015Hz - 0.08Hz</td><td>0.015Hz - 0.08Hz</td></tr><tr><td>momp</td><td>HRPGS</td><td>0.015Hz - 0.08Hz</td><td>0.015Hz - 0.08Hz</td><td>0.015Hz - 0.08Hz</td></tr><tr><td>onec</td><td>HRPGS</td><td>0.015Hz - 0.08Hz</td><td>0.015Hz - 0.08Hz</td><td>0.015Hz - 0.08Hz</td></tr><tr><td>pdns</td><td>HRPGS</td><td>0.015Hz - 0.08Hz</td><td>0.015Hz - 0.08Hz</td><td>0.015Hz - 0.08Hz</td></tr><tr><td>ISPT</td><td>Strong motion</td><td>N/A</td><td>0.015Hz - 0.08Hz</td><td>0.015Hz - 0.08Hz</td></tr><tr><td>PDNS</td><td>Strong motion</td><td>0.037Hz - 0.08Hz</td><td>0.037Hz - 0.08Hz</td><td>0.015Hz - 0.08Hz</td></tr><tr><td>LGCB</td><td>Strong motion</td><td>0.015Hz - 0.08Hz</td><td>0.015Hz - 0.08Hz</td><td>0.015Hz - 0.08Hz</td></tr><tr><td>AATC</td><td>Strong motion</td><td>0.028Hz - 0.08Hz</td><td>0.028Hz - 0.08Hz</td><td>0.032Hz - 0.08Hz</td></tr><tr><td>AES1</td><td>Strong motion</td><td>N/A</td><td>N/A</td><td>0.015Hz - 0.08Hz</td></tr><tr><td>AMNT</td><td>Strong motion</td><td>N/A</td><td>0.015Hz - 0.08Hz</td><td>0.015Hz - 0.08Hz</td></tr><tr><td>APED</td><td>Strong motion</td><td>0.035Hz - 0.08Hz</td><td>0.035Hz - 0.08Hz</td><td>0.035Hz - 0.08Hz</td></tr><tr><td>ATON</td><td>Strong motion</td><td>0.015Hz - 0.08Hz</td><td>0.015Hz - 0.08Hz</td><td>0.015Hz - 0.08Hz</td></tr><tr><td>AV18</td><td>Strong motion</td><td>0.015Hz - 0.08Hz</td><td>0.015Hz - 0.08Hz</td><td>0.015Hz - 0.08Hz</td></tr><tr><td>AV21</td><td>Strong motion</td><td>0.032Hz - 0.08Hz</td><td>0.015Hz - 0.08Hz</td><td>0.015Hz - 0.08Hz</td></tr></tbody></table>

<figure><img src="image_15.png" /><figcaption>Figure S1: Parametrization of the megathrust interface used for the coupling inversion. Coupling value is inverted at each nodes</figcaption></figure>
<figure><img src="image_16.png" /></figure>

Figure S2: Different models variability of the P-wave, S-wave, and density as a function of depth in central Ecuador. A layered model used in this study for Green's function [GF] calculations is plotted as a solid black line. The blue line represents the CRUTST2.0 model in the area (http://igppweb.ucsd.edu/~gabi/rem.html). The other models are from (Vallee et al., 2013; Bethoux et al., 2011; Nocquet et al., 2017). Grey histograms are the probability density function representing our confidence level on the elastic properties, as used to build the model prediction error.
<figure><img src="image_17.png" /><figcaption>Figure S3: **Posterior Median coupling model.** Thin black lines represent the fault parametrization. Coupling values are inverted at each nodes. Interseismic GPS displacement and predictions for the median model are plotted as black and blue arrows, respectively.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>Figure S4: **Posterior Mean coupling model for a coarse parametrisation.** Same as Figure 2a. in the main text but obtained with a coarser fault parametrisation.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>Figure S5: **Decimated InSAR observations, predictions, and residuals.** (a, d, g) Decimated InSAR observations inverted in this study. (b, e, h) Predictions for the posterior mean model. (c, f, i) Residuals of the Sentinel (top row), descending ALOS-2 (middle row), and ascending ALOS-2 (bottom row) interferograms.</figcaption></figure>
<figure><img src="image_20.png" /></figure>

Figure S6: **Empirical covariance functions for the InSAR observations 1D** empirical covariance functions and the associated best-fit exponential function for each tracks. For each image, we compute the empirical covariance as a function of the distance between pixels and then fit an exponential function to these covariances (Jolivet et al., 2012). This exponential function is then used to build the data covariance matrix used in the inversion.
<figure><img src="image_21.png" /><figcaption>Figure S7: Parametrization of the megathrust interface used for the co-seismic inversion

The coloured plane represent the slab1.0 model (Hayes et al., 2012). Each subfault patch is a 15 km x 15 km square</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>Figure S8: **Strong-motion observations and model predictions not presented in Figure 6 in the main text.** The North (left) and vertical (right) components of each station are plotted around the map. For each waveform, the bold number indicates it's maximum amplitude. The station azimuth $\Phi$ and distance $d$ to the epicenter are also given. The black line is the recorded waveform. The gray lines are the stochastic predictions for our posterior model. The red line is the mean of the stochastic predictions.</figcaption></figure>
<figure><img src="image_23.png" /></figure>

Figure S9: Comparison of co-seismic moment and moment deficit in the 1958 earthquake region. a) The background colour represents the coupling posterior mean model. The black dashed lines delimit four different areas where the co-seismic moment of the 1958 event and moment deficit for the 1906 - 1958 period are computed. b-e) Probability densities of the co-seismic moment released by the 1958 earthquake and the moment deficit accumulated between 1906 and 1958 in the different dashed area shown in a).
<figure><img src="image_24.png" /></figure>

Figure S10: Comparison of model predictions and 1942 earthquake waveform recorded in the DBN station, Netherlands. (top) Slip model and hypocenter location (red star) used to compute the predictions shown in the bottom row. The model presented in a) results from the kinematic slip inversion of the 2016 earthquake. The models in b) and c) use a different hypocenter located between the two main slip asperities. The slip model in c) is the same as in a) and b), but located updip along the megathrust interface. Black lines in c) are slip contours of the original slip model. (bottom) East component waveform recorded at DBN for the 1942 earthquake (in black) and stochastic predictions (in grey) for the model shown on top. The red line is the posterior mean prediction. Predictions were convolved with the instrumental response of the Galitzin that recorded the event.
References

Bethoux, N., Segovia, M., Alvarez, V., Collot, J.-Y., Charvis, P., Gailler, A., and Monfret, T. (2011). Seismological study of the central Ecuadorian margin: Evidence of upper plate deformation. *Journal of South American Earth Sciences*, 31(1):139-152.

Charlier, C. and Van Gils, J. (1953). Liste des stations seismologiques mondiales. *Association internationale de seismologie, Observatoire Royal de Belgique*.

Hayes, G. P., Wald, D. J., and Johnson, R. L. (2012). Slab1. 0: A three-dimensional model of global subduction zone geometries. *Journal of Geophysical Research: Solid Earth*, 117(B1).

Jolivet, R., Lasserre, C., Doin, M.-P., Guillaso, S., Peltzer, G., Dailu, R., Sun, J., Shen, Z.-K., and Xu, X. (2012). Shallow creep on the Haiyuan fault (Gansu, China) revealed by SAR interferometry. *Journal of Geophysical Research: Solid Earth*, 117(B6).

Kikuchi, M. and Kanamori, H. (2003). *Note on Teleseismic Body-Wave Inversion Program*. http://www.eri.u-tokyo.ac.jp/ETAL/KIKUCHI/.

Kikuchi, Masayuki and Kanamori, Hiroo (1982). Inversion of complex body waves. *Bulletin of the Seismological Society of America*, 72(2):491-506.

Nocquet, J.-M., Jarrin, P., Vallée, M., Mothes, P., Grandin, R., Rolandone, F., Delouis, B., Yepes, H., Font, Y., Fuentes, D., et al. (2017). Supercycle at the Ecuadorian subduction zone revealed after the 2016 Pedernales earthquake. *Nature Geoscience*, 10(2):145-149.

Vallee, M., Nocquet, J.-M., Battaglia, J., Font, Y., Segovia, M., Regnier, M., Mothes, P., Jarrin, P., Cisneros, D., Vaca, S., et al. (2013). Intense interface seismicity triggered by a shallow slow slip event in the Central Ecuador subduction zone. *Journal of Geophysical Research: Solid Earth*, 118(6):2965-2981.

Ye, L., Lay, T., Koper, K. D., Smalley Jr., R., Rivera, L., Bevis, M. G., Zakrajsek, A. F., and Teferle, F. N. (2014). Complementary slip distributions of the August 4, 2003 Mw 7.6 and November 17, 2013 Mw 7.8 South Scotia Ridge earthquakes. *Earth Planet. Sci. Lett.*, 401:215-226.