Rupture process of the 2016 Mw 7.8 Ecuador earthquake from joint inversion of InSAR data and teleseismic P waveforms

Accepted Manuscript Rupture process of the 2016 Mw 7.8 Ecuador earthquake from joint inversion of InSAR data and teleseismic P waveforms

Lei Yi, Caijun Xu, Yangmao Wen, Xu Zhang, Guoyan Jiang PII: DOI: Reference:

S0040-1951(17)30453-5 doi:10.1016/j.tecto.2017.10.028 TECTO 127666

To appear in:

Tectonophysics

Received date: Revised date: Accepted date:

11 March 2017 24 October 2017 29 October 2017

Please cite this article as: Lei Yi, Caijun Xu, Yangmao Wen, Xu Zhang, Guoyan Jiang , Rupture process of the 2016 Mw 7.8 Ecuador earthquake from joint inversion of InSAR data and teleseismic P waveforms. The address for the corresponding author was captured as affiliation for all authors. Please check if appropriate. Tecto(2017), doi:10.1016/ j.tecto.2017.10.028

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Rupture Process of the 2016 Mw 7.8 Ecuador Earthquake from Joint Inversion of InSAR Data and Teleseismic P Waveforms

PT

Lei Yi1, 2, Caijun Xu1, 2, 3, *, Yangmao Wen1, 2, 3, Xu Zhang4, Guoyan Jiang3 1. School of Geodesy and Geomatics, Wuhan University, Wuhan, 430079, China.

RI

2. Key Laboratory of Geospace Environment and Geodesy of the Ministry of Education, Wuhan University, Wuhan, 430079, China.

SC

3. Collaborative Innovation Center of Geospatial Technology, Wuhan University, Wuhan, 430079,

NU

China.

4. Institute of Geophysics, China Earthquake Administration, Beijing, 100081, China.

Abstract: The 2016 Ecuador earthquake ruptured the Ecuador-Colombia subduction

MA

interface where several historic megathrust earthquakes had occurred. In order to determine a detailed rupture model, Interferometric Synthetic Aperture Radar (InSAR)

D

images and teleseismic data sets were objectively weighted by using a modified

PT E

Akaika’s Bayesian Information Criterion (ABIC) method to jointly invert for the rupture process of the earthquake. In modeling the rupture process, a constrained waveform length method, unlike the traditional subjective selected waveform length

CE

method, was used since the lengths of inverted waveforms were strictly constrained by the rupture velocity and rise time (the slip duration time). The optimal rupture

AC

velocity and rise time of the earthquake were estimated from grid search, which were determined to be 2.0 km/s and 20 s, respectively. The inverted model shows that the event is dominated by thrust movement and the released moment is 5.75 × 1020 Nm (Mw 7.77). The slip distribution extends southward along the Ecuador coast line in an elongated stripe at a depth between 10 and 25 km. The slip model is composed of two asperities and slipped over 4 m. The source time function is approximate 80 s that separated into two segments corresponding to the two asperities. The small magnitude of the slip occurred in the updip section of the fault plane resulted in small tsunami

ACCEPTED MANUSCRIPT waves that were verified by observations near the coast. We suggest a possible situation that the rupture zone of the 2016 earthquake is likely not overlapped with that of the 1942 earthquake. Keywords: 2016 Ecuador earthquake; Constrained waveform length method; Joint inversion; 1942 Ecuador earthquake.

PT

1 Introduction

RI

Recurrent megathrust earthquakes have produced widespread concern in the

SC

public and scientific communities. The North Andean convergent margin, i.e., the western margin of the South America Plate (SAP), experiences intense crustal

NU

deformation under a plate convergence rate of approximate 4.7 cm/yr (Nocquet et al., 2014; Gutscher et al., 1999; Trenkamp et al., 2002; Chlieh et al., 2014). Several large

MA

and great historic earthquakes have occurred along this subduction zone. Their temporal and spatial relationships have been extensively studied, but their relationships have not yet been clarified due to lack of direct seismic observations. In

to

the

results

PT E

According

D

2016, a large subduction earthquake happened on April 16 near the coast of Ecuador. of

the

U.S.

Geological

Survey

(USGS)

(http://earthquakes.usgs.gov; last accessed, 9 July, 2016), the moment magnitude of

CE

the 2016 Ecuador earthquake is Mw 7.8 and its hypocenter (0.382° N, 79.922° W; 21 km) is approximately 29 km south-southeast of Muisne, Ecuador (Fig. 1). Focal such

as

the

Global

Centroid

Moment

Tensor

(GCMT)

AC

mechanisms,

(http://www.globalcmt.org; last accessed, 9 July, 2016) and the W-phase moment tensor (USGS), are consistent with subduction of the Nazca Plate under the SAP. This earthquake provides an opportunity to obtain the first detailed slip model and analyze the relationships between the rupture ranges of large thrust earthquakes in this region. Before the 2016 Ecuador earthquake, this segment of the subduction zone had previously ruptured producing the 1906 (Mw 8.4) megathrust subduction earthquake (Yoshimoto et al., 2017) and the later 1942 (Mw 7.8) earthquake (Ye et al., 2016). Subsequently, the rupture zone extended to the north by the 1958 (Mw 7.6) and 1979

ACCEPTED MANUSCRIPT (Mw 8.1) earthquakes (Ye et al., 2016) (Fig. 1). The rupture zones of the latter three events indicate that they have re-ruptured the rupture zone of the 1906 event, having a 500 km rupture length along the Ecuador-Colombia subduction interface. But they released only one fifth of the total moment of the 1906 earthquake (Kelleher, 1972; Kanamori and McNally, 1982). The temporal and spatial relationships between the

PT

four earthquakes have been investigated in many studies and have been interpreted by two distinct models: the asperity and barrier rupture models (Kelleher, 1972;

RI

Kanamori and McNally, 1982; Beck and Ruff, 1984; Mendoza and Dewey, 1984). Although the source mechanisms of these events have been widely discussed, their

SC

detailed slip models are still not clear. After the 2016 earthquake, Ye et al. (2016) used teleseismic and tsunami data to invert for its rupture model. They claimed that the

NU

2016 earthquake is a quasi-repeat of the 1942 earthquake and partial re-rupture of the 1906 Colombia-Ecuador earthquake. Nocquet et al. (2016) utilized the teleseismic

MA

waveforms, InSAR images, nearfield high-rate and static GPS and acceleration records to invert for the rupture process. They also reviewed the evidence included in

D

Ye et al. (2016) and concluded that all information available concur to a large overlap

PT E

between the 1942 and 2016 earthquakes but cannot be demonstrated due to the scarce data available for the 1942 earthquake. And they proposed a supercycle model to emphasize large seismic hazard remains after great earthquakes in the subduction

CE

interface. Similar inference of the spatial correlation between the 1942 and 2016 earthquakes was also suggested in He et al. (2017) and Yoshimoto et al. (2017).

AC

Through simulating the tsunami records to infer the dominant slip area of the 1906 earthquake, however, Yoshimoto et al. (2017) suggested that the 1942 and 2016 earthquakes were not overlapped with the 1906 event. One goal of this study is to further investigate the correlations between these earthquakes by inverting the rupture model of the 2016 Ecuador earthquake and using a constrained waveform length method in our waveform inversions. Determining a reasonable rupture velocity, i.e., the propagation velocity of the first window, is essential in modeling the rupture process if we use the multi-time window inversion method (Hartzell and Heaton, 1983). When searching for an

ACCEPTED MANUSCRIPT optimal rupture velocity within a preset range, the rupture front delay time on the fault plane is significantly varied. The lengths of the inverted waveforms, however, are usually subjectively selected and are generally set to be equal, which means they are not adjusted to the various velocities. A bias underlies this method since the allowable ruptured space is larger when using a fast velocity rather than using a slow velocity,

PT

particularly when large multi-time windows are applied (Fig. 2). Whereas, in our waveform inversions, the first start time and last end time of the inverted waveforms

RI

are both objectively calculated, and they are constrained by the rupture velocity and rise time (duration time of slip) (Fig. 2). Our method introduces a modification to the

SC

traditional method in order to consider the various rupture velocities. In this study, three data sets were utilized, i.e., ascending and descending Interferometric Synthetic

NU

Aperture Radar (InSAR) images and 34 teleseismic P waveform records. We estimated the optimal rupture velocity and rise time of the earthquake from a grid

MA

search. Then, we jointly inverted for the rupture process and obtained a detailed rupture model for this subduction zone. Finally, based on the inverted slip model of

D

the 2016 Ecuador earthquake, we discuss the spatial correlations between the rupture

PT E

zones of large earthquakes in this region.

CE

2. Data and Method

Teleseismic P waves and InSAR images were both included in our inversions

AC

(Fig. 1). Using data downloaded from the Data Manager Center (DMC) of the Incorporated Research Institutions for Seismology (IRIS), we screened the teleseismic waves according to waveform quality and azimuth distribution and manually picked out their P wave arrival times. Then a 0.005 - 0.2 Hz Butterworth filter was applied to the obtained 34 vertical P waveforms. Instead of using these records all in displacement waveforms, we selected half of them to be used in velocity waveforms (Fig. S1). Thus displacement and velocity waveforms were inverted in an attempt to increase the resolution of the teleseismic body waveforms.

ACCEPTED MANUSCRIPT InSAR images, i.e., descending track T040D (acquired on 12 and 24 April, 2016) and ascending track T018A (acquired on 29 March and 22 April, 2016), acquired from the European Space Agency Sentinel-1A satellite were used to derive the line-of-sight (LOS) surface deformation caused by the earthquake. All SAR images were acquired in the novel Terrain Observation with Progressive Scans (TOPS) acquisition mode

PT

(Torres et al., 2012), which provided a large swath width of 250 km and resolutions of 5 m and 20 m along the range and azimuth directions, respectively. All of the

RI

interferograms were processed using the Switzerland GAMMA software (Werner et al., 2001). For TOPS interferometry, a method that considers the effects of the scene

SC

topography and a method that determines spectral diversity (Scheiber and Moreira, 2000) were used to ensure an accuracy of a few thousandths of a pixel co-registration

NU

in the azimuth direction. The effects of topography-correlated were removed from the interferograms using agency-provided orbits together with the 90-m resolution Shuttle

MA

Radar Topography Mission (SRTM) digital elevation model (DEM) (Farr et al., 2000). And the atmospheric noise in the interferograms is corrected using the

D

elevation-dependent empirical corrections (Scott and Lohman, 2016). Then, the

PT E

interferograms were filtered using a power spectrum filter (Goldstein and Werner., 1998) to reduce the effects of phase noise, after which they were unwrapped using the Minimum Cost Flow (MCF) algorithm. These interferograms were geocoded to the

CE

WGS84 geographic coordinates with a resolution of 90 m. To reduce the number and correlation of observations during the inversion, the InSAR data were downsampled

AC

using the quadtree sampling method (Jónsson et al., 2002). The location-dependent look vectors of the downsampled data points were calculated based on the differences between the incidence and azimuth angles according to the precise orbit data and local topography. Finally, we extracted 552 and 568 InSAR observations from the T040D and T018A interferograms, respectively (Fig. S2). A curved fault plane with varied dip angles was used to simplify the subduction interface in the Slab 1.0 model along the Ecuador coast (Hayes et al., 2012) (Fig. S3). Its strike angle (strike = 26°) was adopted from the W-phase moment tensor result (USGS). According to the IRIS, the fault plane cross the hypocenter location

ACCEPTED MANUSCRIPT (0.372°N, 79.940°W; 20 km) that occurred at the 10th and 6th grids along the strike and dip directions, respectively. Although the fault plane was discretized into subfaults with different intervals, the ground projected subfaults were square patches (15 km × 15 km). The fault plane parameter file is available in the supplementary material.

PT

In modelling the rupture process, after discretizing the assumed causative fault plane into subfaults, we used the traditional multi-time window method (Hartzell and

RI

Heaton, 1983) to parameterize the source time function in subfaults (Fig. 2a). Then we calculated the synthetic waveforms generated from a unit slip in each subfault.

SC

According to the superposition principle, all of the synthetic waveforms radiated from the fault plane are integrated into one waveform to fit the recorded waveform at each

NU

station. These synthetic waveforms are arranged on the rupture time axis starting from the initially ruptured hypocenter (Fig. 2). Furthermore, the start time (ST) and end

MA

time (ET) of these synthetic waveforms form the ST and ET of the inverted waveform (Fig. 2), which we focus on in this study.

D

To explain the constrained waveform length method, we first introduce the

PT E

convolution equation U  G * S for one subfault, where U , G , and S represent the synthetic waveform, Green function, and local source time function, respectively.

CE

Then, we use the equation LU  LG  TS  1 to denote the length relationship between these parameters, where LU is the waveform length, LG is the Green function

AC

length, and TS is the rupture duration time of the subfault. When assuming the initial time of the whole rupture process is starting from the hypocenter, the onset time of the synthetic waveform needs consider three delay times in the rupture process, i.e., the rupture front propagation delay time in the fault plane ( tdelay ), the onset delay time of the initial rupture time in the subfault ( tonset ) (the onset delay time is assumed to be included in tdelay ) and the propagation delay time in the earth ( t prop ) (see Fig. 2a). And for convenience, we assume the onset time of the Green function is its arrival

ACCEPTED MANUSCRIPT time without delay time, thus TS can be reformulated as TS  tdelay  t prop  trise , where trise is the rise time in the subfault. The ST of the synthetic waveform is constrained by the rupture front delay time tdelay and the propagation delay time t prop . The ET of the synthetic waveform is constrained by the ST, the rise time trise ,

PT

and the Green function length LG .

After setting the fault model, i.e., constructing the geometric position

RI

relationships between the stations and subfaults, the propagation delay time t prop is

SC

then calculated. After selecting the appropriate length for the Green function

NU

generated from an impulse point source, the Green function length LG is fixed. Thus, the ST and ET of the synthetic waveform are only restricted by tdelay and trise , i.e.,

MA

ST (Vrup) and ET (trise ,Vrup) , where Vrup is the rupture front expanding velocity

(i.e., rupture velocity) that determines the tdelay . Therefore, according to the described

D

relationships, all restricted synthetic waveforms are assembled together to fit the

PT E

recorded waveform. At the same time, the lengths of all inverted waveforms are objectively calculated that are constrained by the rupture velocity and rise time (Fig.

CE

2b).

Fig. 2 illustrates how the assembling process works in the waveform inversion as

AC

well as the difference between our method and the traditional one. In the traditional case, the ST of synthetic waveforms is calculated, but not the ET. The lengths of the inverted waveforms are always the same. The varied rupture velocity is not adjusted in the waveform length when searching for proper velocity and rise time. By enforcing equal contributions from all of the subfaults in our method, however, the ST and ET of the inverted waveforms are both strictly constrained and varied with the rupture velocity and the rupture velocity and rise time, respectively. By using constrained waveform lengths, we are able to use a two-dimensional grid search to determine the optimal rupture velocity and rise time. The rupture

ACCEPTED MANUSCRIPT velocity is searched in a range from 1.0 km/s to 3.5 km/s with an interval of 0.5 km/s. The local source time function in a subfault is discretized into mutually overlapping isosceles triangular functions with 2 s half widths. The number of triangular functions varies from 1 to 11 with an interval of 2. It is the same for all subfaults. Thus, the rise time is located in a range from 4 s to 24 s with an interval of 4 s. The slip directions of

PT

these triangular functions are restricted between 45° and 135°. Using the Crust 1.0 velocity model (Laske et al., 2013), the static and waveform Green functions were

RI

calculated using the code of Wang et al (2003) and Kikuchi and Kanamori (1991), respectively. The waveform Green functions were filtered the same way as the P

SC

waveforms.

After constructing the observation equations from different data sets, the

NU

covariance matrices were also calculated. Instead of determining the complex covariance components of the P waveforms, the covariance matrix of the P waveforms

MA

is simplified to a diagonal matrix and the variance of each waveform is taken as the variance calculated from their initial waveforms before the P wave arrival. The

D

uncertainty characteristics of the interferograms were generally described using a 1 D

PT E

covariance function in the non-deforming areas (Hanssen, 2001). The estimated standard deviations and decay lengths of the T040D and T018A images are (1.58 cm, 23.8 km) and (1.08 cm, 20.5 km), respectively. The covariance matrices of the T040D

CE

and T018A images were simplified the same way as the P waveforms. Then, the covariance matrices of the T040D and T018A images were determined using this

AC

function.

Spatial and temporal smoothing constraints were used to regularize the rupture model (Yagi et al., 2004). The weighted ratios between observations and constraints were objectively determined using a modified Akaika’s Bayesian Information Criterion (ABIC) expression (Yi et al., 2017). To apply this expression, we first formulate the inversion equation as

d  Ha+e ,

(1)

and 0=Sa+ ,

(2)

ACCEPTED MANUSCRIPT where d= dT 040 D

dT 018 A d P N1 is an N dimensional vector containing 3 data T

sets. The three subscripts, i.e., “T040D”, “T018A” and “P”, represent the three data

H=  HT 040 D

sets. S= ST

S S N

HT 018 A

H P N  M T

is

N M

an

coefficient

matrix.

is an N 2  M coefficient matrix containing 2 types of prior

T

2 M

constraints. The two subscripts, i.e., “T” and “S”, denote the time and space

varies with the time windows. e= eT 040 D

PT

smoothness constraints. a is an M dimensional model parameter vector that

eT 018 A eP N1 and   T

 S N

RI

T

T 2 1

are

In

e~Q( i 2 )

addition,

(i  T 040D, T 018A, P )

is

the

,

covariance

Q( i 2 )   T 040 D 2Q

where

matrix

NU

respectively.

SC

the vectors of the Gaussian distribution errors from observations and constraints,

N N

the

observations.

, where

QT 040 D ,

MA

Q=blkdiag QT 040 D  T 018 A2 /  T 040 D 2QT 018 A  P 2 /  T 040 D 2QP 

of

,

QT 018 A and Q P are the covariance matrices of the three data sets. At the same time,  ( i 2 )   T 2 

 S 2 / T 2I S  N N , where IT and I S are the covariance

PT E

constraints.  =blkdiag IT

(i  T , S ) is the covariance matrix of the

D

 ~ ( i 2 ) , where

2

2

matrices of the two constraints.  T 040 D 2 ,  T 018 A2 ,  P 2 ,  S 2 , and  T 2 are the

CE

hyperparameters to be determined.

Since the spatial and temporal smoothness constraints are utilized, the current

AC

ABIC function can be expressed as

ABIC  N log s(a* )  M log  2  log HT Q1H   2 (ST S) 1  log Q  log (ST S)  1  N  log(2 )  log N  1

,

(3)

where  2   T 040 D 2 /  T 2 and s(a)  (d  Ha)T Q1 (d  Ha)   2aT (ST S) 1a . When the hyperparameters are already known, the solution of as a* (Yabuki and Matsu’ura, 1992).

s(a) =0 is solved and noted a

ACCEPTED MANUSCRIPT The ABIC expression can be applied in our joint inversions when three observation data sets, i.e., two InSAR data sets and one waveform data set, are used simultaneously. By searching large enough ranges of hyperparameters, the optimal weight ratios are chosen when the ABIC value reaches the minimum, the similar to previous works (Yabuki and Matu’ura, 1992; Fukuhata et al., 2003, 2004; Funning et

PT

al., 2014). In all of our inversions, including the grid search and the final inversion, the weighted ratios between the observations and constraints are objectively

3. Results

NU

3.1 Rupture velocity and rise time

SC

RI

determined.

MA

In the grid search, we inverted 36 slip models by jointly applying the three data sets (Fig. 3a); meanwhile, another two groups were simultaneously carried out by jointly inverting the waveforms and T040D image (Fig. 3c) and individually inverting

D

the waveforms (Fig. 3d). The variance reductions (VR) (Cohee and Beroza, 1994) of

PT E

the waveforms, which contained the temporal rupture information, varied with respect to the rise time at different rupture velocities are denoted as curves in Fig. 3. These

CE

curves represent the influence of adopting different rise times and overall, the curves increase with the rise time. The VR curves in Fig. 3d are higher than the other two

AC

groups since they were not compromised by the added InSAR observations. The curves are very close when the rupture velocity is 2.5 km/s and 2.0 km/s, indicating that the rupture velocity is not well resolved when applying only the teleseismic waveforms. The pattern of these curves varying with the rupture velocity is consistent among the three cases. They first increase, and then decrease when the rupture velocity is smaller than 2.0 km/s. When the rupture velocity becomes too small (Fig. 4a), the length of inverted waveforms increases fast. Fitting these overly long waveforms might include other complex waves in waveform inversions which supposed to be the primary reason of the decrease of the VR curves. From a more

ACCEPTED MANUSCRIPT detailed analysis shown in Fig. 3a, it can be seen that when the rupture velocity is 2.0 km/s, the VR curve initially rises then reaches the highest turning point. The curve eventually drops, which is consistent with the rupture velocity being caused by fitting overly long waveforms when the rise time becomes too large (Fig. 4c). This feature also appears in another two groups (Figs. 3c and 3d). According to these results, the

PT

optimal rupture velocity and rise time are supposed to be 2.0 km/s and 20 s (Fig. 3b), respectively.

RI

As first stated, the optimal results are supported by the VR of the P waveforms, which is only one aspect of the grid search results. This indicates that it is not a

SC

critical reason for the determination of the optimal velocity and rise time. Using the back-projection method from teleseismic P waveforms, however, the optimal rupture

NU

velocity is also inverted to be within 2.0 km/s and 2.5 km/s (Ye et al., 2016). This supports our result independently, and is particularly similar to the results shown in

MA

Fig. 3d. Nocquet et al. (2016) also obtained a similar average rupture velocity of 2.3 km/s using the nearfield records. Therefore, our velocity is reliable considering our

D

results are supported from three inversion groups that the data sets are objectively

PT E

weighted and the constrained waveform length method is used in modelling the rupture process. The inverted slip models and the waveform fits of four selected waveforms are depicted in Fig. 4. The source time functions are shown in Fig. 5,

CE

which are separated into two 40 s long time windows.

AC

3.2 Source model for the Ecuador earthquake With the optimal search results, the detailed source model is inverted by narrowing down the search intervals in the weight ratio determination (Fig. 1). The result reveals that the seismic moment of the whole rupture is 5.75 × 1020 N m, which is equivalent to a moment magnitude of Mw 7.77. In addition, the earthquake is dominated by a thrust mechanism that is consistent with motion along the convergent boundary between the subducted Nazca Plate and the North Andes block. The inverted result is available in the supplementary material.

ACCEPTED MANUSCRIPT Jointly constrained by the displacement and velocity waveforms, the obtained slip model is composed of two distinct asperities spreading along the Ecuador coastline and its source time function is approximate 80 s that separated into two segments corresponding to the two asperities (Fig. 1a). The cumulated slip distributions at 10 s intervals show the expansion process of the rupture front (Fig. 1b).

PT

During the first 10 s, the rupture front of the first asperity surrounds the hypocenter in a small slip magnitude. Then the rupture of the first asperity is accomplished in the

RI

subsequent 30 s and the cumulated slip distribution mainly orients southward along the strike direction at a slightly shallower depth than the hypocenter. The slip pauses

SC

around the epicenter of the 1942 event located by Kelleher (1972) that around the coastline similar with the relocated epicenter from ISC-GEM (Storchak et al., 2015),

NU

while east of the epicenter given in Beauval et al. (2013). We do not distinguish between the epicenter of Beauval et al. (2013) and that of Mendoza and Dewey (1984)

MA

since both are obviously off the coast (Fig. 1). The second asperity released the largest portion of the event energy after 40 s and lasted in approximately 40 s. Its slip

D

continually propagates southward at a slightly deeper depth, which forms the second

PT E

large slip area. This is well recognized from the source time function that corresponded to the second pulse starting at around 40 s. The inverted slip distribution is well consistent with the grid search results given in Fig. 4. And it is well surrounded

CE

by the aftershocks (Fig. 1).

Most of the waveforms are well matched in both amplitude and phase (Fig. 6b),

AC

except for several waveforms from stations located in the Pacific Ocean. This is presumably because the seismic zone affected by the subducted Nazca Plate cannot be appropriately modelled using a 1-D velocity model. The LOS displacement fields are also well recovered (Fig. 6a), except for those to the east of the Coastal Ridge. This may be due to the local, highly variable topography and the adopted simple fault model. However, the amplitude of these misfits is quite small, so they may contribute little to the slip model.

ACCEPTED MANUSCRIPT 4. Discussion 4.1 Constrained waveform length The lengths of the inverted waveforms in our inversions vary since both the ST and ET of the inverted waveforms are calculated and constrained by the rupture

PT

velocity and rise time. Our waveform inversion method strictly obeys the

RI

superposition principle by using equal synthetic waveforms from all of the subfaults (Fig. 2b). To accomplish this goal, we use the equal length of the Green functions (our

SC

method) instead of the equal length of the inverted waveform records (the traditional method), which is an essential difference. After convoluting the Green functions with

NU

the local source time functions, the synthetic waveforms are all included in the inversion. On the same fault plane, the variable rupture velocities generate obviously

MA

variable rupture front delay times (Fig. 2c), particularly when the velocity is small, which led to notable variations in the waveform lengths (Fig. 4a) and the

D

corresponding slip models (Fig. 4b). At the same time, the variable rise times have

PT E

also introduced adjustments to the variable waveform lengths (Fig. 4c), and the slip models are also variable (Fig. 4d). Thus, estimating a proper rise time is also necessary.

CE

On the contrary, if the traditional method is used, the superposition principle will be distorted when the inverted waveform lengths are comparable or even shorter than

AC

the source time function length. For instance, according to the model of the 2016 Ecuador earthquake released by the USGS, the inverted waveform lengths are indeed close to the obtained source time function, which is approximate 70 s, similar to Ye et al. (2016). This is because the synthetic waveforms are not completely used since only the initial parts of the synthetic waveforms from the last ruptured subfaults are utilized in the inversion (Fig. 2). Hence, the traditional method is suggested to be inappropriate because the rupture signal is biased when short waveform records are inverted.

ACCEPTED MANUSCRIPT The problem of possibly losing rupture signals in the traditional case, however, can be reduced when adopting long enough lengths for all records. In this case, the actual significance of our method is probably small. While critically evaluating the influence of using these two methods is difficult due to many factors affect the obtained slip models, e.g., the subjectively selected lengths of inverted waveforms

PT

(Fig. 4), the constrained waveform length method is still an improvement of the traditional method (Nakahara, 2008) since it introduces the effects of variable rupture

RI

velocities and rise times. In addition, the essential advantage of our method is the complete and individual control of the synthetic waveforms generated from the

NU

improvement with respect to the traditional one.

SC

discretized source time functions. Therefore, our method potentially represents an

MA

4.2 Slip model for the 2016 Ecuador earthquake According to the inverted slip models in the grid search (Figs. 4b and 4d), we recognize that the dominant slipped area, i.e., the two asperities, is always to the south

D

of the epicenters of the 2016 earthquake. This reveals that the stabilized locations of

PT E

the two asperities are independent of the significantly varied rupture velocities and rise times. In the preferred rupture model for the 2016 earthquake when the rupture

CE

velocity is 2.0 km/s and the rise time is 20 s, its slip distribution is apparently separated into two distinct asperities in an oblique elongated stripe (Fig. 1). The two

AC

identified asperities are also observed in previous slip results (Ye et al., 2016; Nocquet et al., 2016; Yoshimoto et al., 2017). The large slipped areas in these models are both consistent with the large interseismic coupling (ISC) rate areas in this subduction zone (Nocquet et al., 2014; Chlieh et al., 2014). The inverted source time functions in the grid search are separated into two segments that both last approximately 40 s (Figs. 1 and 5), consistent with the successive rupture of the two separated asperities (Nocquet et al., 2016) (Fig. 1). Though two asperities are also observed in previous slip models (Nocquet et al., 2016; Ye et al., 2016; Yoshimoto et al., 2017), their ruptured duration times are less than 70

ACCEPTED MANUSCRIPT s and no obvious two segments is observed in their source time functions. The duration time is smaller than our results of approximate 80 s, including a rupture front delay time of at least 70 s propagating approximately 140 km with a 2 km/s rupture velocity and a 20 s rise time assumed for all subfaults. According to the slip snapshot in Fig. 1b, the slip magnitude and rupture area of the second asperity ruptured after 40

PT

s are both larger than these of the first one. This means the second segment of the source time function larger than the first one is reasonable. And the slip distribution

RI

after 70 s varies very slightly except small magnitude slips out of the main slip area (Fig. S4). This may be because the subfaults outside the main slip area are included in

SC

our inversions, so the duration time and the amplitude of the second segment (ruptured after 40 s) are probably overestimated. Thus, the proper size of the allowed

NU

ruptured fault plane is also an important factor, which needs to be considered. The reason we constructed a slightly larger fault plane is because we are

MA

concerned about the shallow slip near the Ecuador-Colombia trench axis, which is an important factor in tsunami warnings. In addition, the amplitudes of the shallow slips

D

that ruptured near the fault edge obtained from the preferred slip model are small,

PT E

which means the trench is not obviously ruptured (Fig. 1). This is consistent with the observed small tsunami waves off the coast (Ye et al., 2016; Nocquet et al., 2016). Using the preferred slip model, we forward synthesized the GPS displacements

CE

that were utilized in Nocquet et al. (2016) and compared with the recorded displacements as depicted in Fig. 7. The overall fitting of the GPS records is

AC

acceptable except for some stations’ records north of the epicenter of the 2016 earthquake, such as MOMP, MUIS, and ESMR. On the other side, the slip model inverted from the GPS records shows no slips north of the epicenter (Fig. S5a) while the slip result inverted from the GPS and InSAR data sets demonstrates some slips in that area (Fig. S5c). In addition, through comparing the synthetic records, we supposed that the slips north of the epicenter are the reason resulting discrepancies in fitting these records (Fig. S5). Thus the misfit records depicted in Fig. 7 are mainly due to the slips north of the epicenter. In the preferred slip model, the slips located

ACCEPTED MANUSCRIPT north of the epicenter are ruptured between 30 s and 50 s (Fig. 1b) and they vanish when fewer weights are assigned by the InSAR observations (Fig. S6). Therefore, we are uncertain about these slips considering the teleseismic P waveforms and inland InSAR LOS fields we used have low resolution for the slip away from the coast (Lay et al., 2011). The low resolution of inland observations was

PT

confirmed in Melgar et al. (2016) that the inland data sets have low resolution to recover the rupture process of the 2015 Illapel (Mw 8.3), Chile earthquake. In

RI

addition, some studies of the 2011 Tohoku-Oki earthquake indicate that though using plenty of inland geodetic observations, the accurate slip amplitude near the trench is

SC

uncertain unless constrained by the near trench observations (Yamazaki et al., 2011; Simons et al., 2011; Fujiwara et al., 2011). Thus, jointly analyzing or inverting the

NU

near trench records, such as GPS records and tsunami waveforms, is a more robust and reliable method of verifying the slips north of the epicenter and the slips near the

MA

trench (Ye et al., 2016). In addition to the observations’ resolution, a comprehensive investigation is made more certain by considering two factors that are simplified in

D

our modeling, i.e., the varied strike and dip angles of the subduction interface near the

PT E

equator (Hayes et al., 2012) and the possible presence of ruptured splay faults in the upper-plate (Collot et al., 2004, 2008).

CE

4.3 Spatial correlations between the 1942 and 2016

AC

earthquakes

Based on the obtained slip model of the 2016 earthquake, we attempt to analyze the uncertain range of large coseismic slip of the 1942 earthquake with the help of previous research results. By analyzing collected records of the 1942 earthquake, its source duration time is inverted and determined to be approximate 24 s (Swenson and Beck, 1996). If its rupture velocity is 2 km/s that approximates to the velocity in many subduction earthquakes, the ruptured asperity should occur in a hypothetical circular area approximately 50 km around the hypocenter of the 1942 earthquake (Swenson and Beck, 1996; Chlieh et al., 2014). The actual range is wider than this circle;

ACCEPTED MANUSCRIPT otherwise its coseismic slip that occurred in the circle is unreasonably large (Swenson and Beck, 1996). However, because of the unknown complexity of the source, we focus on this identified ruptured asperity around the hypocenter. Even so, since the epicenter position of the 1942 earthquake has two types of results, i.e., located near the coastline (Kelleher, 1972; Storchak et al., 2015) and off the coast (Mendoza and

PT

Dewey, 1984; Swenson and Beck, 1996; Beauval et al., 2013) (Fig. 1), this introduces discrepancies in identifying the location of the asperity (Fig. 8). It is important in our

RI

analyses since our obtained slip model of the 2016 earthquake has distinct patterns, i.e., overlapping and adjacent, to the epicenter of the 1942 earthquake (Fig. 8). More

SC

precisely, because of the two asperities in our slip model are located in deeper and smaller ranges than previous models that covering both the two epicenters (e.g., Ye et

NU

al., 2016; Nocquet et al., 2016; Yoshimoto et al., 2017), the indefinite positions of the epicenter of the 1942 earthquake generate particular concerns since showing a

MA

possible different situation discussed below.

The epicenter near the coastline (Kelleher, 1972; Storchak et al., 2015) was

D

supposed in Ye et al. (2016), which is surrounded by their inverted slip model.

PT E

Together with the similarities of the macro isoseismic maps and the recorded seismograms of the 1942 and the 2016 earthquakes, Ye et al. (2016) claimed that the 2016 Ecuador earthquake is a quasi-repeat of the 1942 earthquake and partial

CE

re-rupture of the 1906 earthquake. This indicates the three events largely re-rupture in the same subduction region. And the recurrence times of large earthquakes in this

AC

seismogenic region significantly vary between the 1906 and 1942 earthquakes (i.e., 36 years) and the 1942 and 2016 earthquakes (i.e., 74 years). Though Nocquet et al. (2016), Yoshimoto et al. (2017) and He et al. (2017) suppose the epicenter is off the coast, their slip models are still covering the two epicenter positions that are similar with Ye et al. (2016). Thus they also largely support the inference in Ye et al. (2016). But since the fault zone is not completely locked (Nocquet et al., 2014; Chlieh et al., 2014), the released seismic moments in the 1942 and the 2016 earthquakes are greater than the maximum accumulated moment in their previous interseismic periods (Nocquet et al., 2016). In addition to the varied recurrence times, this causes an

ACCEPTED MANUSCRIPT essential seismic moment deficit issue in this subduction interface. Nocquet et al. (2016) proposed a supercycle model to explain this essential issue, i.e., a super earthquake cycle rather than the general earthquake cycle between successively ruptured large earthquakes, e.g., the 1942, 1958 and 2016 earthquakes, to understand the seismic moment deficit problem. In other words, the released seismic moment in

PT

the 1942 event larger than the accumulated moment between 1906 and 1942 is coming from previously left moment in the great 1906 event, so too the 2016 event. It

RI

is due to the released seismic moments in the great 1906 earthquake is not all the accumulated seismic moment in a long interseismic period. Though we also enjoy the

SC

supercycle model, the varied recurrence times and the seismic moment deficit issue are both under the assumption that the rupture zones between these large and great

NU

earthquakes are largely overlapped, which we supposed a possible different situation. On the other hand, through simulating the recorded tsunami waveforms,

MA

Yoshimoto et al. (2017) inverted for the rough slip model of the 1906 earthquake. Their result indicates the dominant slip area of the 1906 event occurs in the shallow

D

subduction interface. This means that the 1906 event possibly has no spatial

PT E

correlations with the 1942 and 2016 events. This is partially supported from the macro isoseismic map of the 1906 event (Swenson and Beck, 1996) that the high intensity area not reaches to the equator where the 1942 mainly ruptured. If the same area of

CE

the 1942 earthquake was largely ruptured in the 1906 event, the isoseismic map of the 1906 event should have some corresponding responses similar to that of the 1942

AC

event which is apparently rejected. Therefore, the moment released in the 1942 earthquake is probably accumulated before 1906 and the released moment larger than the accumulated moment between 1906 and 1942 is not the unreleased strain from the great 1906 earthquake. Consequently, the moment deficit problem of the 1942 event is questionable and the previously mentioned moment deficit issue mainly reduces to the moment deficit problem of the 2016 earthquake since its seismic moment exceeds the accumulated moment since 1942. Reviewing the different epicenter positions of the 1942 earthquake previously stated, it has to be noted that though the relocated epicenter position of the 1942

ACCEPTED MANUSCRIPT earthquake in Storchak et al. (2015) is near the coastline, an accurate position is not reachable due to lack of direct observations. Similar with many studies (e.g., Swenson and Beck, 1996; Beauval et al., 2013; Chlieh et al., 2014; Nocquet et al., 2016; He et al., 2017; Yoshimoto et al., 2017), we also suggest the epicenter is possibly off the coast that supports from several less rigorous clues. Firstly, the source depth of the

PT

1942 earthquake is inverted to be 14 km (Swenson and Beck, 1996) indicating a relatively shallower hypocenter near the trench axis (Mendoza and Dewey, 1984).

RI

Secondly, the shallower hypocenter is also supported from the macro isoseismic maps of the 1942 and 1958 earthquakes. The seriously damaged area in the macro

SC

isoseismic map of the 1958 earthquake is obviously smaller than that of the 1942 earthquake (Swenson and Beck, 1996). The possible reason is the dominant slipped

NU

area of the 1958 earthquake is deeper than that of the 1942 earthquake. Different from the 1942 event, the rupture range of the 1958 event is less ambiguous and supposed

MA

around the hypocenter surrounded by aftershocks. And we are aware of the source depth of the 1958 earthquake that is similar to the near coastline epicenter position of

D

the 1942 earthquake while deeper than the shallower position off the Ecuador coast.

PT E

Thus it prefers a relative shallower ruptured asperity of the 1942 event. It is consistent with the macro isoseismic map of the 1942 earthquake showing wide open circles toward the ocean (Swenson and Beck, 1996; Ye et al., 2016). Then the reasonable

CE

epicenter position of the 1942 earthquake is shallow and supposed off the coast. Thirdly, the ruptured asperity in the 1942 earthquake occurred in shallower depth is

AC

partially supported from the ISC models in Nocquet et al. (2014) and Chlieh et al. (2014). These models show a consistent pattern that the ISC rates in the plate interface near the Ecuador coast gradually increase with depth. This suggests that in different regions of the interface, particularly different depths, different interseismic periods are needed to accumulate released energy in order to reach similar magnitudes, i.e., the 1942 and 2016 earthquakes. Then the largely varied interseismic periods between these two events indicate that they are ruptured in different areas. The slip model of the 2016 event in our work occurs in a slightly deeper position near the coastline than previous models (Ye et al., 2016; Nocquet et al., 2016; Yoshimoto et al.,

ACCEPTED MANUSCRIPT 2017; He et al., 2017), thus the mainly ruptured area in the 1942 event is suggested away from the coast, which is consistent with the macro isoseismic map of the 1942 earthquake. Therefore, the asperity around the epicenter of the 1942 earthquake is possibly off the coast. Besides these inferences, our obtained slip models in the grid search (Figs. 4b

PT

and 4d) and the preferred slip model (Fig. 1) explicitly do not cover the shallower epicenter. Therefore, a possible situation is suggested, i.e., the ruptured asperity of the

RI

1942 earthquake does not largely overlap with the two identified asperities of the 2016 earthquake (Fig. 8). Accordingly, the seismic moment deficit issue is defused

SC

since the released moment in the 2016 earthquake is probably accumulated before 1942. Thus, to understand the varied recurrence times and the seismic moment deficit

NU

issue, we need carry out a more credible investigation considering all possible situations between the ruptured zones of these successively ruptured large and great

MA

earthquakes in the subduction interface.

Leaving out the inaccurate epicenter position of the 1942 earthquake and the

D

assumed limited rupture area used in previous inference, a possible situation can be

PT E

that the ruptured asperity of the 1942 event is adjacent rather than overlapped to that of the 2016 event. Indeed, due to lack of direct observations and undoubted evidence, it has to be reminded that the rupture zone of the 1942 earthquake is undetermined

CE

and previous discussions are not rigorous and require further investigation and reanalysis based on more robust results. At the same time, the supercyle model

AC

proposed in Nocquet et al. (2016) is a well explanation to the seismic moment deficit problem and the high seismic hazard after great earthquake in the subduction interface they emphasized is well recognized and confirmed by the series large earthquakes following the great 1906 earthquake (Fig. 8). Meanwhile, there is a region near the trench (Chlieh et al., 2014) that is surrounded by the aftershocks of the 1942 earthquake, over a 200 km × 90 km area parallel to the coast (Swenson and Beck, 1996) (Fig. 8). This region has not been ruptured since at least 1942, indicating a high seismic risk. But the ISC rate along the shallow seismic zone is low (Chlieh et al., 2014), indicating less potential for the

ACCEPTED MANUSCRIPT accumulation of large amounts of strain in this region. Thus, whether this region is a potential site for the next large rupture depends on whether it has been ruptured and how fast tectonic strain is being accumulated. The recurrence time and source mechanism of large and great subduction earthquakes in this seismic zone can be further estimated if the correlations of this region and the ranges of the 1906 and 1942

PT

earthquakes are known.

RI

5. Conclusions

SC

The rupture model for the 2016 Ecuador earthquake was derived from the descending and ascending InSAR observations and the teleseismic P waveforms. We

NU

adopted a method that uses constrained lengths of inverted waveforms in the waveform inversion. Our results indicate that the slip model for the Ecuador

MA

earthquake is well represented by two distinct asperities and that its slip distribution mainly extends southward with increasing released seismic moments. The temporal and spatial relationships between the rupture zone of the 2016 earthquake and those of

D

previous large and great thrust earthquakes need to be investigated further. We suggest

PT E

a possible situation that the rupture zone of the 2016 earthquake is likely not overlapped with that of the 1942 earthquake. This detailed slip model could enhance

CE

our understanding of the heterogeneity of seismic risks in this subduction zone, which experiences recurrent megathrust earthquakes.

AC

Acknowledgement We thank K. Wang, J.M. Nocquet and another anonymous reviewer for their constructive comments and suggestions on the present work. This work is co-supported by the National Natural Science Foundation of China (key program) under grants No.41431069 and 4151101233, the National Key Basic Research Development Program (973 program) under grants No. 2013CB733303 and 2013CB733304. The teleseismic data are provided by the IRIS (DMS data center, http://www.iris.edu.hq). The Sentinel-1 SAR images are provided by ESA through

ACCEPTED MANUSCRIPT Sentinels Scientific Data Hub (SciHub, https://scihub.copernicus.eu). Aftershocks and focal mechanisms were obtained from the IRIS, GCMT and USGS. The numerical calculations have been done on the supercomputing system in the Supercomputing Center of Wuhan University. This work made use of GMT software (Wessel and Smith, 1998).

PT

Reference

RI

Beauval, C., Yepes, H., Palacios, P., Segovia, M., Alvarado, A., Font, Y., Aguilar, J.,

SC

Troncoso, L., Vaca, S., 2013. An earthquake catalog for seismic hazard assessment in Ecuador. Bulletin of the Seismological Society of America 103,

NU

773-786.

Beck, S.L., Ruff, L.J., 1984. The rupture process of the great 1979 Colombia

MA

earthquake: evidence for the asperity model. Journal of Geophysical Research 89, 9281-9291.

Chlieh, M., Mothes, P.A., Nocquet, J.M., Jarrin, P., Charvis, P., Cisneros, D., Font, Y.,

D

Collot, J.Y., Villegas-Lanza, J.C., Rolandone, F., Vallee, M., Regnier, M.,

PT E

Segovia, M., Martin, X., Yepes, H., 2014. Distribution of discrete seismic asperities and aseismic slip along the Ecuadorian megathrust. Earth and

CE

Planetary Science Letters 400, 292-301. Cohee, B.P., Beroza, G.C., 1994. Slip distribution of the 1992 Landers Earthquake

AC

and its implications for earthquake source mechanics. Bulletin of the Seismological Society of America 84, 692-712. Collot, J.Y., Marcaillou, B., Sage, F., Michaud, F., Agudelo, W., Charvis, P., Graindorge, D., Gutscher, M.A., Spence, G., 2004. Are rupture zone limits of great subduction earthquakes controlled by upper plate structures? Evidence from multichannel seismic reflection data acquired across the northern Ecuador-southwest Colombia margin. Journal of Geophysical Research 109(B11103), doi:10.1029/2004JB003060.

ACCEPTED MANUSCRIPT Collot, J.Y., Agudelo, W., Ribodetti, A., Marcaillou, B., 2008. Origin of a crustal splay fault and its relation to the seismogenic zone and underplating at the erosional north Ecuador-south Colombia oceanic margin. Journal of Geophysical Research 113(B12102), doi:10.1029/2008JB005691. Farr, T., Rosen, P., Caro, E., 2000. The shuttle radar topography mission. Reviews of

PT

Geophysics 45, 37-55. Fujiwara, T., Kodaira, S., No, T., Kaiho, Y., Takahashi, N., Kaneda, Y., 2011. The

334(6060), 1240, doi:10.1126/science.1211554.

RI

2011 Tohoku-Oki earthquake: displacement reaching the trench axis. Science

SC

Fukahata, Y., Yagi, Y., Matu’ura, M., 2003. Waveform inversion for ABIC for seismic source processes with two sorts of prior constraints: comparison between proper

NU

and improper formulations. Geophysical Research Letters 30(6): 2002GL16293. Fukahata, Y., Nishitani, A., Matu’ura, M., 2004. Geodetic data inversion using ABIC

MA

to estimate slip history during one earthquake cycle with viscoelastic slip-response functions. Geophysical Journal International 156, 140-153.

D

Funning, G.J., Fukahata, Y., Yagi, Y., Parsons, B., 2014. A method for the joint

PT E

inversion of geodetic and seismic waveform data using ABIC: application to the 1997 Manyi, Tibet, earthquake. Geophysical Journal International 196, 1564-1579.

CE

Goldstein, R., Werner, C., 1998. Radar interferogram filtering for geophysical applications. Geophysical Research Letters 25, 4035-4038.

AC

Gutscher, M.A., Malavieille, J., Lallemand, S., Collot, J.Y., 1999. Tectonic segmentation of the North Andean Margin: impact of the Carnegie Ridge collision. Earth and Planetary Science Letters 168, 255-270. Hanssen, R.F., 2001. Radar Interferometry: Data Interpretation and Error Analysis, Kluwer Academic Publishers: Dordrecht, The Netherlands. Hartzell, S.H., Heaton, T.H., 1983. Inversion of strong ground motion and teleseismic waveform data for the fault rupture history of the 1979 Imperial Valley, California, earthquake. Bulletin of the Seismological Society of America 73(6), 1553-15583.

ACCEPTED MANUSCRIPT Hayes, G.P., Wald, D.J., Johnson, R.L., 2012. Slab 1.0: a three-dimensional model of global subduction zone geometries. Journal of Geophysical Research 117, B01302. He P., Hetland, E.A., Wang, Q., Ding, K., Wen, Y., Zhou, R., 2017. Coseismic slip in the 2016 Mw 7.8 Ecuador Earthquake Imaged from Sentinel-1A Radar

PT

Interferometry. Seismological Research Letters 88, 277-286. Jónsson, S., Zebker, H., Segall, P., Amelung, F., 2002. Fault slip distribution of the

RI

1999 Mw 7.1 Hector Mine, California, earthquake, estimated from satellite radar and GPS measurements. Bulletin of the Seismological Society of America 92,

SC

1377-1389.

Kanamori, H., McNally, K.C., 1982. Variable rupture mode of the subduction zone

NU

along the Ecuador-Colombia coast. Bulletin of the Seismological Society of America 72, 1241-1253.

MA

Kelleher, J., 1972. Rupture zones of large South American earthquakes and some predictions. Journal of Geophysical Research 77, 2087-2103.

D

Kikuchi, M., Kanamori, H., 1991. Inversion of complex body waves-III. Bulletin of

PT E

the Seismological Society of America 81, 2335-2350. Laske, G., Masters, G., Ma, Z., Pasyanos, M., 2013. Update on Crust1.0 – A 1-degree Global Model of Earth’s Crust. Geophysical Research Abstracts 15, Abstract

CE

EGU2013-2658, 2013.

Lay, T., Ammon, C.J., Kanamori, H., Xue, L., Kim, M.J., 2011. Possible large

AC

neartrench slip during the 2011 Mw 9.0 off the Pacific coast of Tohoku earthquake. Earth Planets and Space 63, 687-692. Melgar, D., Fan, W., Riquelme, S., Geng, J., Liang, C., Fuentes, M., Vargas, G., Allen, R.M., Shearer, P.M., Fielding, E.J., 2016. Slip segmentation and slow rupture to the trench during the 2015, Mw8.3 Illapel, Chile earthquake. Geophysical Research Letters 43, 961-966. Mendoza, C., Dewey,

J.W., 1984. Seismicity associated with the great

Colombia-Ecuador earthquakes of 1942, 1958, and 1979: implications for barrier

ACCEPTED MANUSCRIPT models of earthquake rupture. Bulletin of the Seismological Society of America 74, 577-593. Nakahara, H., 2008. Seismogram envelope inversion for high-frequency seismic energy radiation from moderate-to-large earthquakes. Advance in Geophysics 50, 401-426.

PT

Nocquet, J.M., Jarrin, P., Vallee, M., Mothes, P.A., Grandin, R., Rolandone, F., Delouis, B., Yepes, H., Font, Y., Fuentes, D., Regnier, M., Laurendeau, A.,

RI

Cisneros, D., Hernandez, S., Sladen, A., Singaucho, J.C., Mora, H., Gomez, J., Montes., L., Charvis, P., 2016. Supercycle at the Ecuadorian subduction zone

SC

revealed after the 2016 Pedernales earthquake. Nature Geoscience 10, 145-149. Nocquet, J.M., Villegas-Lanza, J.C., Chlieh, M., Mothes, P.A., Rolandone, F., Jarrin,

NU

P., Cisneros, D., Alvarado, A., Audin, L., Bondoux, F., Martin, X., Font, Y., Régnier, M., Vallée, M., Tran, T., Beauval, C., Mendoza, J.M. Maguiña, Martinez,

MA

W., Tavera, H., Yepes, H., 2014. Motion of continental slivers and creeping subduction in the northern Andes. Nature Geoscience 7, 287-291.

D

Ratzov, G., Collot, J.Y., Sosson, M., Migeon, S., 2010. Mass-transport deposits in the

PT E

northern Ecuador subduction trench: result of frontal erosion over multiple seismic cycles. Earth and Planetary Science Letters 296, 89-102. Scheiber, R., Moreira, A., 2000. Coregistration of interferometric SAR images using

CE

spectral diversity. IEEE Trans. Geoscience Remote Sensing 38, 2179-2191. Scott, C.P., Lohman, R.B., 2016. Sensitivity of earthquake source inversions to

AC

atmospheric noise and corrections of InSAR data. Journal of Geophysical Research 121, 4031-4044. 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., Webb, F.H., 2011. The 2011 magnitude 9.0 Tohoku-Oki earthquake: Mosaicking the megathrust from seconds to centuries. Science 332, 1421-425. Storchak D A，Giacomo D Di，Engdahl E R，Harris J，Bondar I，Lee W H K，Bormann P，Villasenor A．2015．The ISC-GEM Global Instrumental Earthquake Catalogue

ACCEPTED MANUSCRIPT (1900-2009)：Introduction．Physics of the Earth and Planetary Interirors，239： 48-63． Swenson, J.L., Beck, S.L., 1996. Historical 1942 Ecuador and 1942 Peru Subduction earthquakes, and earthquake cycles along Colombia-Ecuador and Peru subduction segments. Pure and Applied Geophysics 146, 67-101.

PT

Torres, R., Snoeij, P., Geudtner, D., Bibby, D., Davidson, M., Attema, E., Potin, P., Rommen, B., Floury, N., Brown, M., Traver, I.N., Deghaye, P., Duesmann, B.,

RI

Rosich, B., Miranda, N., Bruno, C., L’Abbate, M., Croci, R., Pietropaolo, A., Hchler, M., Rostan, F., 2012. GMES Sentinel-1 mission. Remote Sensing of

SC

Environment 120, 9-24.

Trenkamp, R., Kellogg, J.N., Freymuellerm, J.T., Mora, H.P., 2002. Wide plate

NU

margin deformation, southern Central America and northwestern South America, CASA GPS observations. Journal of South America Earth Sciences 15, 157-171.

MA

Wang, R., Martin, F.L., Roth, F., 2003. Computation of deformation induced by earthquakes in a multi-layered elastic crust - Fortran programs EDGRN/EDCMP.

D

Computers & Geosciences 29, 195-207.

PT E

Werner, C., Wegmüller, U., Strozzi, T., Wiesmann, A., 2001. GAMMA SAR and interferometric processing software, In Proceedings of the ERS ENVISAT Symposium. Gothenburg, Sweden, 16-20 October 2001.

CE

Wessel, P., Smith, W.H.F., 1998. New, improved version of generic mapping tools released. EOS Trans, AGU, 79, 579.

AC

White, S.M., Trenkamp, R., Kellogg, J.N., 2003. Recent crustal deformation and the earthquake cycle along the Ecuador-Colombia subduction zone. Earth and Planetary Science Letters 216, 231-242. Yabuki T., Matu’ura, M., 1992. Geodetic data inversion using a Bayesian information criterion for spatial distribution of fault slip. Geophysical Journal International 109, 363-375. Yagi, Y., Takeshi, M., Javier, P., Gabriel, R., 2004. Source rupture process of the Tecoman, Colima, Mexico earthquake of 22 January 2003, Determined by joint

ACCEPTED MANUSCRIPT inversion of teleseismic body-wave and near-source data. Bulletin of the Seismological Society of America, 94, 1795-1807. Yamazaki, Y., Lay, T., Cheung, K.F., Yue, H., Kanamori, H., 2011. Modeling near-field tsunami observations to improve finite-fault slip models for the 11 March 2011 Tohoku earthquake. Geophysical Research Letters 38, 231-248.

PT

Yepes, H., Audin, L., Alvarado, A., Beauval, C., Aguilar, J., Font, Y., Cotton, F., 2016. A new view for the geodynamics of Ecuador: implication in seismogenic source

RI

definition and seismic hazard assessment. Tectonics 35, 1249-1279. Ye, L.L., Kanamori, H., Avouac, J.P., Li, L., Cheung, K.F., Lay, T., 2016. The 16 April

SC

2016, Mw 7.8 (Ms 7.5) Ecuador earthquake: A quasi-repeat of the 1942 Ms 7.5 earthquake and partial re-rupture of the 1906 Ms 8.6 Colombia-Ecuador

NU

earthquake. Earth and Planetary Science Letters 454, 248-258. Yi, L., Xu, C., Zhang, X., Wen, Y., Jiang, G., Li, M., Wang, Y., 2017. Joint inversion

MA

of GPS, InSAR and teleseismic data sets for the rupture process of the 2015 Gorkha, Nepal, earthquake using a generalized ABIC method. Journal of Asian

D

Earth Sciences 148, 121-130.

PT E

Yoshimoto, M., Kumagai, H., Acero, W., Ponce, G., Vasconez, F., Arrais, S., Ruiz, M., Alvarado, A., Garcia, P.P., Dionicio, V., Chamorro, O., Maeda, Y., Nakano, M., 2017. Depth-dependent rupture mode along the Ecuador-Colombia subduction

AC

CE

zone. Geophysical Research Letters 44, 2203-2210.

ACCEPTED MANUSCRIPT

PT E

D

MA

NU

SC

RI

PT

Figure caption

Fig. 1. (a) Surface map projection of the inverted slip model for the 2016 Ecuador earthquake. The blue triangles in the inset map are the teleseismic stations included in

CE

this study. The red boxes show the spatial coverage of the two SAR data. The red star and red dots denote the mainshock’s epicenter and its aftershocks in the first two

AC

months (USGS), respectively. The yellow and cyan stars represent epicenters of previous large earthquakes from Kelleher (1972) and digitized from Beauval et al. (2013), respectively. The thick purple line represents the rupture zones of these earthquakes modified from White et al. (2003). ECT = Ecuador-Colombia Trench, and DGM = Dolores-Guayaquil Megashear. (b) Cumulative slip snapshots of the rupture model. The circles and triangles are the aftershocks of the 1942 and 1958 earthquakes from Kelleher (1972), respectively.

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

D

Fig. 2. Illustration of the constrained waveform length method. (a) Three delay times

PT E

included in the rupture process, i.e., t-rupture delay ( tdelay ), t-onset delay ( tonset ) and t-propagation delay ( t prop ). Four selected subfaults in the hypothetic fault plane, i.e., 1,

CE

2, i, and n, denote the initial, second, arbitrary, and final ruptured subfaults, respectively. (b) The synthetic waveforms from all subfaults are arranged along the

AC

rupture time axis considering the delay times. The start time (ST) and end time (ET) of these synthetic waveforms, i.e., marked as Pi and Ei (i=1, n), respectively, are strictly calculated. In the constrained waveform length method, Pi and Ei are objectively calculated. But in the traditional method, Ei is subjectively selected and usually is fixed equally. (c) If the length of inverted waveforms is smaller than our method, i.e., incomplete synthetic waveforms are used in the inversion, the contributions from last ruptured subfaults are discounted. Thus the rupture space is distorted, i.e., the red line circled area exceed the preset rupture duration time.

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

Fig. 3. Variance reductions of the P waveforms from the inverted slip models in the

MA

grid search, using (a) two LOS fields and P waveforms, (c) T040D LOS data and P waveforms, and (d) P waveforms. Curves in (a), (c), and (d) with different colors represent different rupture velocities the values of which are labeled in the legend

D

(unit: km/s). The pentagram in (a) represents the optimal 2 km/s rupture velocity and

PT E

20 s rise time. (b) is the contour map of (a). The hexagram represents the optimal

AC

CE

result.

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

Fig. 4. (a, c) The observed (black) and synthetic (red) P waveforms of four stations

MA

derived from the slip models listed in (b, d). Waveforms in six horizontal groups derived from inversion results using (a) different rupture velocities marked at the left

D

side when the rise time is 20 s, and (c) different rise times when the rupture velocity is

PT E

2.0 km/s. The coefficient of the waveform similarity is labeled on the right side of each subfigure. The corresponding six slip models are listed in (b) and (d) with their velocities and rise times marked in top left corner of each subfigure. Large slips are

CE

always near the epicenter of the 1942 earthquake near the coastline (Kellher, 1972)

AC

rather than the epicenter off the coast (Beauval et al., 2013).

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

D

Fig. 5. Source time functions from slip models in Figs. 4b and 4d. They are inverted

PT E

results from using (a) different rupture velocities when the rise time is 20 s and (b) different rise times when the rupture velocity is 2 km/s. The moment magnitude of each result is labeled at middle right and the adopted (a) rupture velocity or (b) rise

CE

time is shown below. The blue dash line denotes the time around 40 s separating the

AC

two asperities in the rupture process.

PT E

D

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

Fig. 6. (a) Observed, synthetic and residual LOS fields for the T040D and T018A tracks from the preferred slip model. (b) The observed (black) and synthetic (red)

CE

velocity and displacement P waveforms. On the time axis under each column, zero represents the P wave arrival time. In each subfigure, the station name is located on

side.

AC

the left side and the coefficient of similarity (up) and distance (down) are on the right

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

Fig. 7. The observed (blue and green) and forward (red) horizontal (a) continue and (b)

AC

CE

PT E

D

campaign GPS displacements using the preferred slip model.

RI

PT

ACCEPTED MANUSCRIPT

Fig. 8. Schematic presentation of the approximate rupture zones (open circles) of the

SC

1942 (Mw 7.8) (red) and 1958 (Mw 7.6) (blue) earthquakes. The cyan open circle denotes areas un-ruptured since at least 1942. The closed green circle inside the

NU

rupture zone of the 1942 earthquake represents the ruptured asperity (Nocquet et al., 20016; Swenson and Beck, 1996; Chlieh et al., 2014). The red dots, circles, and

AC

CE

PT E

D

MA

triangles are aftershocks (the same as in Fig. 1).

ACCEPTED MANUSCRIPT

PT E

D

MA

NU

SC

RI

PT

Supplementary material

Fig. S1. Locations of teleseismic stations used in our inversions. The cyan circle dots

CE

and blue squares represent displacement and velocity records, respectively. Their

AC

station names are next to the dots and squares.

D

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

PT E

Figure S2. Retrieved LOS fields from the (a) T040D and (b) T018A images. The

AC

CE

amplitude decreases fast and the magnitude is small in remote area.

SC

RI

PT

ACCEPTED MANUSCRIPT

NU

Fig. S3. Side view of the curved fault plane. The red dots are the central of subfaults in the curved fault plane. The circles are the sampled points from the Slab 1.0 model

MA

(Hayes et al., 2012). The strike angle of fault plane is 26° that is the same as the

AC

CE

PT E

D

W-phase moment tensor result (USGS).

SC

RI

PT

ACCEPTED MANUSCRIPT

AC

CE

PT E

D

MA

NU

Fig. S4. Rupture animation of the preferred source model of the 2016 Ecuador earthquake.

PT E

D

MA

NU

SC

RI

PT

ACCEPTED MANUSCRIPT

CE

Fig. S5. Slip models from inverting the (a) continue and (b) campaign GPS records, and from jointly inverting the (c) continue and (d) campaign GPS records and the

AC

InSAR data used in the main text. The forward (red) and recorded (blue and green) horizontal displacements are compared in two different scales. Some records at stations north of the epicenter of the 2016 earthquake, i.e., MOMP, MUIS and ESMR (c, d), are far from well fitted which are supposed to be affected by the slips north of the epicenter after comparing the slip distributions and fitness in (a, b) and (c, d).

SC

RI

PT

ACCEPTED MANUSCRIPT

Fig. S6. Comparing slip distributions from (a) the preferred slip model (Fig. 1) that

NU

optimal weighted the three datasets and (b) one slip model contributed fewer from the InSAR observations (less weight). When sharing (b) fewer contribution from the

MA

InSAR observations, the slips north of the epicenter of the 2016 earthquake were vanishing and there were slips toward the fault edge near the epicenter of the 1998

AC

CE

PT E

D

earthquake.

PT

ACCEPTED MANUSCRIPT

AC

CE

PT E

D

MA

NU

SC

RI

Graphical abstract

ACCEPTED MANUSCRIPT Highlights: The rupture process of the 2016 Ecuador earthquake was weightedly inverted from two InSAR images and one teleseismic P waveform dataset. The lengths of inverted waveforms were constrained by the rupture velocity and rise time that optimal estimated from grid search.

PT

The two large asperities of the 2016 Ecuador earthquake is likely not overlapped with

AC

CE

PT E

D

MA

NU

SC

RI

that of the 1942 earthquake.