Crustal motion and deformation in Ecuador from cGNSS time series

Accepted Manuscript Crustal motion and deformation in Ecuador from cGNSS time series Alejandra Staller, José Antonio Álvarez-Gómez, Marco P. Luna, Marta Béjar-Pizarro, Jorge M. Gaspar-Escribano, Sandra Martínez-Cuevas PII:

S0895-9811(17)30357-7

DOI:

10.1016/j.jsames.2018.05.014

Reference:

SAMES 1936

To appear in:

Journal of South American Earth Sciences

Received Date: 27 August 2017 Revised Date:

23 May 2018

Accepted Date: 28 May 2018

Please cite this article as: Staller, A., Álvarez-Gómez, José.Antonio., Luna, M.P., Béjar-Pizarro, M., Gaspar-Escribano, J.M., Martínez-Cuevas, S., Crustal motion and deformation in Ecuador from cGNSS time series, Journal of South American Earth Sciences (2018), doi: 10.1016/j.jsames.2018.05.014. 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

1

CRUSTAL MOTION AND DEFORMATION IN ECUADOR

2

FROM cGNSS TIME SERIES

3

Abstract

5

In this paper, we present the first velocity field from cGNSS (Continuous GNSS) stations in the

6

Continuous Monitoring GNSS Network (REGME) in Ecuador. We have analyzed data from 33

7

cGNSS REGME stations for the 2008-2014 period in order to characterize horizontal crustal motion

8

and deformation in Ecuador. Prior to this, we analyzed the time series for the 33 REGME stations

9

in order to determine their seasonality and the type of spectral noise. For most stations, we found

10

a predominance of uncorrelated white noise with annual and semi-annual variations as the

11

predominant first and second periods. Velocity was estimated by introducing the trend,

12

seasonality and noise in each series in the general model, thus allowing us to improve accuracy as

13

well as magnitude. The velocity and strain distribution correspond to the transpressive right-

14

lateral slip of the westward-dipping faults of the Major Dextral System and the NNE movement of

15

the North Andean Block relative to the South American plate. The distributions of our deformation

16

rate and velocity field indicate a differentiated tectonic behavior between northern, central and

17

southern Ecuador. In northern Ecuador, there is an estimated right-lateral motion of 7.6 ± 0.5

18

mm/yr, consistent with the NNE movement of the NAB relative to the South American plate. In

19

central Ecuador, the right-lateral motion decreases to 5.3 ± 0.4 mm/yr. In the southern region of

20

Ecuador (from the Guayaquil Gulf to Peru) there is no strain accumulation, GNSS velocities

21

decrease and turn to south, this zone belongs to the so-called Inca or Peru sliver.

22

These results are consistent with the distinct behavior of subduction in Ecuador, with no coupling

23

south of Ecuador, and increased coupling towards the north in the zone where megathrust

24

earthquakes have occurred over the last century. The southern part of the Carnegie Range marks

25

the limit between the two zones.

26

We suggest that the main driving force responsible for ongoing crustal deformation in Ecuador is

27

the convergence between the Nazca and South American plates with the variable coupling pattern

28

and the collision of the Carnegie Range. This results in different velocity patterns for the north and

29

the south.

AC C

EP

TE D

M AN U

SC

RI PT

4

1

ACCEPTED MANUSCRIPT

30 31

Key Words: GNSS time series; velocity field; crustal deformation; strain rates; Ecuador; North

32

Andean Block.

RI PT

33

1. INTRODUCTION

35

Ecuador is located in a complex tectonic setting, controlled by the interaction between the Nazca,

36

South American plates and North Andean block (NAB) (Figure 1). The main tectonic mechanism is

37

the subduction of the Nazca plate under the South American plate at 55-58 mm/yr in a N83°E

38

direction along the Ecuadorian margin. Of particular significance for the study area is the oceanic

39

Carnegie Ridge, located within the Nazca plate and formed by its motion over the Galapagos

40

hotspot (Hey, 1977). This ridge contributed to the uplift of the Coastal Ranges (Gutscher et al.,

41

1999) at rates of 0.1 – 0.5 mm/yr (Pedoja et al., 2009), to the generation of an erosive subduction

42

margin (Sage et al., 2006), and to the northeastward “escape” of the North Andean crustal block

43

(Trenkamp et al. 2002). The motion of the North Andean crustal block produces a complex system

44

of active faults that generate shallow-focus earthquakes on the eastern front of the Andes and in

45

the Sub-Andean region (Nocquet et al., 2014).

46

The rapid oblique convergence between the Nazca and South American plates causes two distinct

47

effects that contribute to the observed GNSS velocity field in Ecuador: 1) locking along the

48

subduction interface, which causes elastic stress accumulation along the subduction interface and

49

induces shortening on the overriding plate. This elastic stress will be released in future

50

earthquakes, 2) long-term motion of the NAB towards NNE with respect to the South American

51

plate.

52

The Ecuador-Colombian margin encompasses two seismically and tectonically contrasted

53

segments (Bethoux et al., 2011): a southern segment (latitude: 0.5°S-3°S) with lower seismicity

54

activity and a northern segment (latitude: 3.5°N-0.5°S) with higher seismicity activity, including

55

several magnitude 7 or greater earthquakes since 1900, such as the 1906 M 8.7 Esmeraldas or the

56

2016 Mw 7.8 Pedernales earthquakes. These events caused more than 7,500 casualties, damaged

57

thousands of buildings and resulted in billions of dollars of economic losses.

AC C

EP

TE D

M AN U

SC

34

2

ACCEPTED MANUSCRIPT

In order to understand the mechanisms controlling seismicity patterns and seismic hazard, it is

59

important to study the deformation field in the region and discuss it the seismotectonic context. In

60

this regard, the availability of cGNSS constitutes a timely and precise resource for the analysis.

61

Continuous and campaign-type GPS measurements have been used to determine crustal motions

62

and deformation in Ecuador and adjacent areas over the last 14 years, (Trenkamp et al., 2002.

63

Nocquet et al., 2009; Nocquet et al., 2014; Drewes and Heidbach., 2012; Sánchez and Drewes,

64

2016). Focusing on Ecuador, Cisneros and Nocquet (2011) presented a velocity field for the

65

country based on measurements at RENAGE (the National GPS Network of Ecuador) passive-

66

landmarks This field is based on 3-hour observation periods carried out in years 1994, 1996 and

67

1998. More recently, the availability of GNSS recording stations operating in Ecuador has

68

increased significantly in the last years. This includes the REGME network. It is a continuous

69

monitoring GNSS network managed by the Geographic Military Institute of Ecuador (IGM)with 44

70

stations in operation to date.

71

In this study, we present the first velocity field from cGNSS stations in the REGME network in

72

Ecuador. We have analyzed data from the available cGNSS REGME stations for the 2008-2014

73

period (only 33 stations) in order to characterize horizontal crustal motion and deformation in

74

Ecuador. We document these new results and discuss them in terms of time series, velocities and

75

strain rates. We propose an analysis of the GPS measurements in Ecuador to explain the state of

76

continental deformation, including velocity partitioning and locking in the subduction interface.

77

These results are interpreted in the light of the Mw 7.8 2016 Pedernales earthquake.

SC

M AN U

TE D

EP

78

RI PT

58

2. SEISMOTECTONIC SETTING

80

Ecuador is located in northwestern South America, in a tectonic setting characterized by the

81

subduction of the oceanic Nazca plate under the continental South American plate (Figure 1). The

82

location and geometry of the subducting slab exert primary control over: 1) active volcanism in the

83

North Central Andes Mountains; 2) transpressional deformation in the upper plate; and 3)

84

seismicity at various depths (Hayes et al., 2012; Guillier et al., 2001).

AC C

79

3

ACCEPTED MANUSCRIPT

The trench is located 40 – 100 km offshore, where the Nazca plate subducts with a N80°–83°E

86

direction beneath South American plate at a convergence rate of 58 mm/yr (Trenkamp et al.,

87

2002).

88

The Carnegie Ridge, which extends from the Galapagos Islands to the subduction trench with a EW

89

direction, subducts beneath South America at latitude 1°N-2°S. This ridge plays a relevant role in

90

the subduction, the degree of interplate coupling and continental deformation (Gutscher et al.,

91

1999; Pedoja et al., 2009; Sage et al., 2006; Michaud et al. 2009; Goyes, 2009; Trenkamp et al.

92

2002). There have been five major earthquakes (1906, 1942, 1958, 1979, 2016) in the northern

93

part of the ridge collision and one in the southern part. Gutscher et al. (1999) indicated that none

94

of these events seems to have broken through the cordillera itself, which seems to exert a prime

95

control in propagating the ruptures of large subduction earthquakes.

96

The Grijalva Fracture Zone is located south of this ridge, with an approximate orientation of N60°E,

97

and its subduction starts at the latitude of the Gulf of Guayaquil (3°S). This fracture zone marks the

98

boundary between the Neogene Nazca crust to the north and the Oligocene oceanic crust to the

99

south (e.g. Pedoja et al., 2009; Collot et al., 2009; Gutscher et al., 1999).

M AN U

SC

RI PT

85

One consequence of the oblique subduction of the Nazca plate and the coupling of the subduction

101

interface, and the subduction of the Carnegie Ridge (Witt and Bourgois, 2009) is the motion of the

102

North Andean block (NAB) towards the NE (Figure 1). This movement results in the opening of the

103

Gulf of Guayaquil and also in the formation of Quaternary NW-SE normal faulting in the region

104

(e.g. Dumont et al., 2005), which delimits the NAB and the South American plate. The structural

105

lineaments between these tectonic segments begin in the Gulf of Guayaquil, crossing the coastal

106

forearc basin and entering through the inter-Andean depressions towards the North through the

107

Andes Mountain Range (e.g. Chunga et al., 2016; Nocquet et al., 2014; Trenkamp et al., 2002)

108

(Figure 1). It is important to indicate that these lineaments re are no neatly defined and

109

correspond to an intraplate collision band that some authors refer to as the Major Dextral System

110

(e.g. Pennington et al., 1981; Kellog and Bonini, 1982; Mann and Burke; 1984; Chunga et al., 2016)

111

and others call the Guayaquil-Caracas Megashear (e.g. Dumont et al., 2005). The boundaries of

112

this active intraplate continental margin is itself subject to numerous controversies (Eguez et al.,

113

2003; Espinoza, 1992; Mendoza & Dewey, 1984). Trenkamp et al. (2002) delimit the NAB by the

AC C

EP

TE D

100

4

ACCEPTED MANUSCRIPT

Boconó Fault, the East Andean Fault System and the Guayaquil-Caracas Megashear to the east, the

115

South Caribbean Deformed Belt to the north, and the Colombia-Ecuador Trench and the micro-

116

block of Panama to the west. The Guayaquil-Caracas Megashear is a system of dextral strike-slip

117

faults with a northeast orientation and north direction reverse faults and constitutes the eastern

118

boundary along which the NAB moves. We will use the Major Dextral System (MDS) denomination

119

to refer to this fault system that delimits the east of NAB.

120

Several authors (e.g. Nocquet et al., 2009; Chunga et al., 2009; Trenkamp et al., 2002; Witt and

121

Bourgois, 2009; Egbue & Kellogg, 2010) suggest that this block is moving N35°E at a rate of 6-10

122

mm/yr. Nocquet et al., 2014 estimates that this block is moving NE at a rate of 7.5-9.5 mm/yr and

123

describes its kinematic motion as starting from the center of Ecuador and heading towards

124

Colombia.

125

The subduction zone has produced damaging events in the past. The strongest event registered to

126

date took place on January 31, 1906 with an estimated moment magnitude of Mw 8.8 (Kanamori

127

and McNally, 1982). Others relevant events in the area include: The May, 14 1942 (Mw 7.8); the

128

January 19, 1958 (Mw 7.8); the December 12, 1979 (Mw 8.1) (Mendoza and Dewey, 1984) and the

129

April 16, 2016 (Mw 7.8) (Nocquet et al., 2017) (see red stars in Figure 1). However, events of lesser

130

magnitudes occurring in the Ecuadorian highlands have caused even greater damage (e.g. Schuster

131

et al., 1996; Beauval et al., 2010) because of their shallow foci and its proximity to populated

132

areas. These events are mainly associated with the MDS that accommodates the motion between

133

the NAB and the South American plate (red circles in Figure 1).

SC

M AN U

TE D

EP

134

RI PT

114

3. GNSS DATA AND PROCESSING

136

At present, the Continuous Monitoring GNSS Network (REGME) is made up of 44 stations

137

distributed throughout the continental and insular territory (Figure 2). The network is managed by

138

the IGM and is part of the SIRGAS Continental Network. It was created, installed and set up

139

according to the international and national standards established by IGM. Most of the stations

140

have been installed on concrete columns measuring 2.0 m high by 0.30 m wide and located on the

141

roofs of public buildings (Figure 3). The first REGME Network station was established in 2008. For

142

this study, we selected 33 stations (Table 1) (active and operating stations in the period covered by

AC C

135

5

ACCEPTED MANUSCRIPT

this study) that had continuous GNSS data available for the period 2008-2014. In addition, we

144

analyzed 17 continuously recording GPS stations (Table 1) from the International GNSS Service

145

(IGS) (Dow et al., 2009) network to establish a link with the global frame of reference ITRF2008

146

(IGb08). These stations are distributed throughout the South American (6), North American (5),

147

Caribbean (4) and Nazca (2) plates (Figure S1 Supplementary Material). To ensure a similar

148

timespan for IGS and REGME data, we selected IGS cGNSS (continuous GNSS) sites with at least 5-

149

years of data available. These IGS cGNSS stations were analyzed for the entire timespan covered

150

by the REGME data: from 2008 to the end of 2014.

151

The data were processed on a daily basis using the Processing Engine (BPE) of the Bernese GPS

152

software 5.0 (Dach et al., 2007), which applies a double-difference processing strategy. We applied

153

IGS precise orbits and earth orientation parameters, with an absolute antenna phase center and a

154

FES2004 ocean tide-loading model (Lyard et al, 2006). Only observations made above an elevation

155

cut-off angle of 5° were used to estimate parameters, and an elevation-dependent weighting was

156

applied. The tropospheric effect was modeled on a prior dry-Niell model and completed by

157

estimating zenith delay corrections for each site at 1-hour intervals using the wet-Niell mapping

158

function (Niell, 1996). The ambiguity resolution is based on the Quasi-Ionosphere-Free (QIF)

159

baseline-wise analysis. Baselines were constructed based on a strategy of maximum observations.

160

Variance-covariance information was taken into account when processing the whole network. As a

161

result, we obtained daily solutions with daily coordinates that were only weakly constrained to the

162

reference frame. In order to ensure adequate internal accuracy within the network, constraints on

163

implementing the chosen reference frame were only imposed a posteriori. The coordinates of the

164

IGS stations were constrained (NNT- No Net Translation option) to their ITRF2008 values. As a

165

result of this processing, we obtained daily SINEX files (Kouba, 1996) with geocentric Cartesian

166

coordinates and variance-covariance matrixes for all REGME and IGS stations for the 2008-2014

167

period. This is the starting point for analyzing the time series.

168

AC C

EP

TE D

M AN U

SC

RI PT

143

169

4. TIME-SERIES ANALYSIS

170

We use the daily SINEX files from Bernese, which contain coordinates and covariances. In order to

171

obtain topocentric time series, the geocentric coordinates (X, Y, Z) were first transformed into

6

ACCEPTED MANUSCRIPT

topocentric coordinates (E, N, Up). Then, we eliminated the remaining outliers, estimating and

173

correcting the offsets defined by changes in the antenna, dome and receiver, and visually

174

inspecting the series. Finally, classified daily position time series by trend, seasonality and type of

175

noise present using spectral analysis, thus improving the magnitude and accuracy of estimated

176

velocities.

177

4.1.

178

The remaining outliers were removed from the time series. Following Wdowinski et al. (1997),

179

outliers were defined as: (a) daily coordinates with formal error (from GNSS evaluation) greater

180

than three times the mean formal error for the site, and/or (b) coordinates that differ by more

181

than three times the weighted root mean square (WRMS) from a regression line. For time series

182

with an offset, separate linear trends were calculated for each of the two parts. Outliers were

183

removed in all three components for every station analyzed, with an average 3% of values

184

eliminated for each station. This process further improved uncertainty by about 3% (WRMS and σ

185

of velocities).

186

Any discontinuous change in position that appeared within the time series that could be attributed

187

to unknown variations of phase center offsets due to a change of antenna or other actions were

188

estimated and corrected. Offset was detected in 7 REGME stations (Figure S2 Supplementary

189

Material), with values between 6.36 mm and 10.76 mm. Only two stations (CUEC and PDEC) had

190

more than two offsets. The corresponding log-files were checked in order to determine the

191

possible causes of these offsets, but since no changes of antenna or other important actions were

192

reported for these dates, it is possible that stations may have been manipulated without due

193

notification. On average, eliminating time series offsets improved the accuracy (RMS) by about

194

22% in the horizontal component and about 4% in the vertical component, as indicated by the

195

RMS values before and after correcting time series offsets Figure S2 Supplementary Material).

196

4.2.

197

It is widely acknowledged that all periodic signals are made up of both the principal frequency and

198

its higher harmonics. Because of this, we assumed that the time series contain both determinant

199

parts and thus applied a functional model that includes trend and seasonality as well as

200

background noise (Amiri-Simkooei et al., 2007). Consequently, after removing outliers and

AC C

EP

TE D

M AN U

SC

Eliminating Remaining Outliers and Correcting Offsets.

RI PT

172

Time-Series Model.

7

ACCEPTED MANUSCRIPT

201

correcting offsets, we modeled the series analyzing the trends, seasonal variations and types of

202

noise, so that the sum of all these components represents the real time series. The theoretical

203

model is: =

204

(1)

is seasonality and

noise, that is to say

206

irregular components and random fluctuations of innovations not explained by other components.

207

A weighted regression line was used to determine an initial value for the trend. Each daily

208

component in the series is associated with an uncertainty derived from the data processing

209

variance-covariance matrix. On average, these values are less than one millimeter in all three

210

components, which is unrealistic and incongruent with the repeatability of the series (see below).

211

Nevertheless, this trend is taken as an initial velocity estimate with its corresponding uncertainty,

212

with average values of 2.4 mm for the horizontal component and 6.5 mm for the vertical

213

component. As it is well known, uncertainties in GNSS velocities are a function of the position

214

repeatability, the length of the time series, the number of data points, the presence of steps in the

215

time series and the noise model (Agnew, 1992).

216

Having determined the first value for the trend, the seasonal signal is modeled, as this can affect

217

relevant parameters that are estimated based on the time series, especially station velocity

218

(Blewitt and Lavallée, 2002; Bos et al., 2008). By observing the time series for each station (Figure

219

S3 Supplementary Material) it is possible to detect any possible deterministic cycles for each

220

component of the series. In order to determine the periodogram and calculate the harmonics, the

221

trend is eliminated from the series, working with the residuals instead. We followed the procedure

222

proposed by Lomb (1976), determining the power spectrum using a periodogram for unequally

223

spaced non-discreet series, as is our case. Peaks in the power spectrum ( ) represent

224

predominant frequencies and periods. Seasonal variations are expressed as:

AC C

EP

TE D

M AN U

SC

RI PT

where

226

is the trend,

+

205

225

is the observed value,

+

( )=∑

=

∑!

( !

. sin(2

∗ sin(2

)+

! )

= ∑!

)) + (2) !

∗ cos(2

! ) (3)

227

where

228

component may have a different fundamental period, it is possible to detect any possible

#

= $/ ,

is the fundamental period,

. cos (2

!

is the residual and $ = 1,2, … /2. As each

8

ACCEPTED MANUSCRIPT

deterministic cycle for a time series. To avoid overestimating the magnitude and accuracy of

230

velocities, we used only a single harmonic for each season () = 1) when analyzing the series.

231

A spectral analysis was carried out for each time series component of all stations in order to

232

estimate the spectral index and amplitude of colored noise. The temporal correlation of

233

coordinates has been detected and studied in recent years, as determining this correlation entails

234

taking various numerous geophysical and/or meteorological effects dependent on the weather

235

into account (García, 2016). Some authors (e.g. Johnson and Agnew, 1995; Zhang et al. 1997; Mao

236

et al., 1999; Williams et al., 2004; Hackl et al., 2011) have already established a correlation

237

between the weather and the presence of noise as well as its effect on estimating uncertainty in

238

the series, indicating that uncertainty is underestimated when only the white noise in the series is

239

taken into account. For many geophysical phenomena, the power spectrum (P) may be fairly

240

accurately calculated by a power-law dependence on frequency formulated as follows (Agnew,

241

1992)

M AN U

SC

RI PT

229

( )=

242

⁄

+ * (4)

243

where f is temporal frequency,

244

Usually, the spectral index is not an integer, and takes values between -3 and 0, for most

245

geophysical phenomena. If −1 < 0 < 0 it is called fractional Gaussian noise and it is generally

246

considered stationary, i. e., its statistical properties do not vary over time. There are three specific

247

types of noise: white noise, flicker and random walk, that result in spectral index values of 0, -1

248

and -2, respectively. It is important to determine the spectral index because this allows us to

249

identify the type of noise present in the series so that we can model it and take it into account

250

when estimating velocities and improving accuracy. The method used to obtain the spectral index

251

value 0 consists of adjusting a straight line log (P (f)) – log (f), whose slope determines its value

252

(Figure S4 and Table S1 Supplementary Material).

253

Once the trend, seasonality and noise type have been obtained, they are introduced in the general

254

model for the series in order to calculate the model parameters as well as periodic variations,

255

series trends and noise. This was done using a minimum quadratic adjustment, estimating the

256

intercept

257

*

are normalization constants and α is the spectral index.

AC C

EP

TE D

-

and

*

*,

slope 2, sinusoidal terms , !

=

*

and discontinuities.

+ 2 ! + . sin(2

!)

+ . cos(2

!)

+ ! (5)

9

ACCEPTED MANUSCRIPT

The values estimated initially using the simple regression line, without taking into account

259

seasonality and noise, vary in absolute value up to 2.6 mm, 1.0 mm and 13.4 mm in the east, north

260

and vertical components, respectively. The improvement in the vertical component of stations

261

MHEC and GYEC (Figure 5) is especially noteworthy. These stations are located on the coast, and

262

their values vary by -3.11 mm and 2.44 mm, respectively. In general, the model described reduces

263

the uncertainties (WRMS) obtained when estimating velocity for all three components, particularly

264

in the vertical component, with improvements of up to 2.8 mm. Figure S3 Supplementary Material

265

shows the time series improvement of daily coordinates for all sites in Ecuador and IGS stations,

266

along with the trend, seasonality, adjusted values and uncertainties for the three components. As

267

may be observed, the adjustment is better suited to series with more dispersed data.

SC

RI PT

258

M AN U

268 269

5. RESULTS

270

5.1.

271

Figure S5 Supplementary Material shows the predominant periods for each component. A peak

272

during the annual seasonal period (~365 days) is observed for most stations: for the east and north

273

components, the annual seasonal period is presented as the first or second predominant period,

274

with some stations containing semi-annual seasonal periods (~180 days) or quarterly seasonal

275

periods (~90 days). Annual periods are more evident in the vertical component. Station QUEM

276

(Figures S3 and S4 Supplementary Material) is worth noting, as it does not present any of these

277

seasonal periods. This may be due to the abundant gaps contained in the observation period (73%

278

of the values are missing in a discontinuous observation period of 2.4 years). No significant

279

differences in cyclical variations between the stations located in the different geographic regions

280

(coast, Andean Range and Amazon) are observed.

281

The periodicity of the stations is more evident in the vertical component, within which we have

282

classified the sites into three categories according to periodicity:

283

(1) Unmodeled periodic ground movement: This category includes series with annual and

284

semiannual seasonality, and refers to the fact that the site moves periodically (Amiri-Simkooei

285

et al., 2007). The stations with annual seasonality are: ALEC (1.03 years), BAHI (1.13 years),

286

CHEC (1.06 years), COEC (1.05 years), CUEC (1.03 years), CXEC (1.07 years), ECEC (1.05 years),

AC C

EP

TE D

Seasonality and Spectral Analysis.

10

ACCEPTED MANUSCRIPT

287

EPEC (1.11 years), ESMR (1.13 years), GZEC (1.03 years), IBEC (0.87 years), LJEC (0.99 years),

288

MAEC (1.11 years), PDEC (1.03 years), PJEC (1.05 years), RIOP (1.03 years), SEEC (1.03 years),

289

STEC (0.92 years), TNEC (1.08 years). The stations with semiannual seasonality were: CLEC

290

(0.73 years), LREC (0.65 years), QUEM (0.59 years), QUI1 (0.41 years). (2) Periodic variation of estimated time series: These correspond to stations that “appear” to

292

move periodically, but only as a result of errors in the system due to the tidal effect produced

293

by time series with spurious systematic effects (Penna and Stewart, 2003; Stewart et al., 2005;

294

Amiri-Simkooei et al., 2007). This category includes the following series: AUCA (1.54 years),

295

EREC (1.55 years), MHEC (1.97 years), MTEC (1.20 years), NJEC (0.21 years), PREC (1.56 years),

296

SNLR (0.06 years).

SC

RI PT

291

(3) Periodic variations correlated with years of observation: This category includes stations whose

298

first fundamental period falls halfway through the observation period or after the full

299

observation period. The category includes the following stations: GYEC (period and years of

300

observation =5.75 years), PTEC (period = 2.27 years; years of observation = 4.53 years); QVEC

301

(period and years of observation = 3.95 years). The second fundamental period for these

302

stations corresponds to shorter time periods; we thus have GYEC (1.92 years), PTEC (0.91

303

years); QVEC (1.32 years).

TE D

M AN U

297

The time series periodicity analysis for some stations is inconclusive because the observation

305

periods are very low (less than 1 year, stations LREC and QUI1) or the time series presents many

306

gaps (station BAHI, Figure S3 in Supplementary Material).

307

Estimated spectral index values are within the range −1 < 0 < 0 (Table S1 and Figure S4 in

308

Supplementary Material), corresponding to fractional Gaussian noise, which is considered

309

stationary. Several authors have suggested that white noise and flicker dominate the noise

310

spectrum for GNSS coordinate time series, as does, to a lesser degree, random walk noise.

311

Similarly, it has been demonstrated that for shorter observation times, white noise is the

312

predominant source of noise, while for longer observations there is a preponderance of random

313

walk (Hackl et al., 2011). Thus, we may expect to find flicker noise in our study, since many of the

314

series have observation periods between 2 and 7 years. Furthermore, 15 stations are located on

315

the rooftops of public buildings and 5 are on the ground (Figure 3). Nevertheless, we did not find

AC C

EP

304

11

ACCEPTED MANUSCRIPT

any correlation between the observation period, the type of noise and station location. Therefore,

317

we can assume that the noise value obtained is not a reflection of the type of noise that was

318

actually present. This may be mostly due to the great number of gaps in most of the series

319

analyzed, given that, despite observation periods of over 5 years, between 22% and 37% of the

320

data are missing (e.g. CUEC, ESMR, GYEC, LJEC, PTEC and RIOP began their operations in 2008 and

321

2009) (Table 1). On average, the percentage of data that could be recovered during the entire

322

observation period was 68.4% (Table 1), because, there is a great number of gaps or missing data,

323

which has an especially strong influence on this analysis. Figure S4 and Table S1 in Supplementary

324

Material show the spectral index values for all components and stations analyzed.

325

5.2.

326

The ITRF2008 velocities for all stations (33 REGME stations and 17 IGS stations) were determined

327

using the model described in Section 4.2. The ITRF2008 velocities (Table S2 in Supplementary

328

Material) were then transformed into a South American fixed reference frame by rotating them

329

about the SOAM/ITRF2008 pole, which was calculated by minimizing the rigid motion of five IGS

330

stations located in South America (Goudarzi et al., 2014), namely BRAZ, BRFT, LPGS, PARC and

331

UNSA (Figure S1 in Supplementary Material). In this frame, the five core stations show residuals

332

below 1.5 mm/yr. This pole differs from that obtained by Sánchez and Drewes, 2016, mainly due

333

to the difference in the number and distribution of stations across South American plate used.

334

However, the difference in the velocity values obtained using either pole is lower than 0.3 mm/yr

335

in the horizontal component. This difference is negligible given the uncertainties resulting for the

336

velocities (see Figure S6 in Supplementary Material).

337

Figure 2 shows our present-day velocity field of Ecuador in a South American reference frame, the

338

95% confidence ellipse is plotted at the tip of each velocity vector. The results are derived using

339

the assumption of constant velocities over the seven-year span (2008–2014). The velocities and

340

their uncertainties are provided in Table 1.

341

Given the short time interval of the time series for the stations LREC, BAHI, QUI1 and EREC (Figure

342

S3 in Supplementary Material), the velocities for these stations have been excluded from analyses

343

and subsequent discussion.

SC

RI PT

316

AC C

EP

TE D

M AN U

GNSS Velocities.

12

ACCEPTED MANUSCRIPT

The predominant characteristic of our velocity field with respect to the South American fixed

345

reference frame is a prevailing eastward movement and a decreased of velocities towards the

346

interior. The orientation turns from north to south towards NE-ENE and SE. This will be analyzed in

347

greater detail in the Discussion section.

348

Our velocity field is consistent with previous studies (Trenkamp et al., 2002; Nocquet et al., 2014;

349

Sánchez and Drewes, 2016). Our velocities have been compared to those obtained by Nocquet et

350

al., 2014 in the 7 stations common to both studies, obtaining average differences of less than 1.4

351

mm/yr for the horizontal component. The greatest difference was found in station GYEC with -2.9

352

mm/yr and 3.9 mm/yr, respectively in the east and north components. Moreover, the comparison

353

of our velocities with the velocity model for South America, VEMOS2015 (Velocity Model of South

354

America 2015, Sánchez and Drewes, 2016), obtained by applying the multiannual solution

355

SIR15P01 by SIRGAS (Geocentric Reference System for the Americas, Sánchez and Drewes, 2016),

356

shows average differences lower than 2 mm/yr in the horizontal component. The greatest

357

differences are found in stations located in the MDS deformation zone. This suggests that the

358

VEMOS2015 model is not capable to register the behavior of stations located in the MDS

359

adequately. There is a noteworthy maximum difference in station ESMR, located on the north

360

coast, with -9,1 mm/yr and -0.6 mm/yr, respectively in the east and north components. This

361

record may be indicating a significant deformation prior to the April 2016 earthquake (see Section

362

7), which cannot be reflected by the model. Nevertheless, it is important to keep in mind that,

363

although the VEMOS2015 model is based on solution SIR15P01, which includes a great number of

364

REGME stations, the observation times for the time series used to obtain said solution are shorter

365

than those used in the time series for the present study.

366

We identified three main groups of stations with more or less homogeneous behavior (Figure 2):

367

(1) The eastern group (east of the MDS) is characterized by low velocities between 2.7 and 3.4

368

mm/yr, which reflect the wavelength of interseismic strain accumulation on the MDS. (2) The NAB

369

group with an average velocity of 16.8 mm/yr and an average orientation of 80 ± 12° (E-ENE).

370

Stations on the coast showed a highly significant eastward component, which suggests that part of

371

the velocity of the Nazca plate is transferred directly to the upper continental plate. This group

372

clearly evinces the NE “escape” of the NAB (Trenkamp et al., 2002; Nocquet et al., 2014), as we

AC C

EP

TE D

M AN U

SC

RI PT

344

13

ACCEPTED MANUSCRIPT

will discuss later. (3) The central group, over the Andes Mountains and over the MDS, show

374

intermediate velocities, with values ranging from 11 to 4 mm/yr and orientations of 75° (E-ENE) to

375

120° (ESE-SE) for stations located to the north and to the south, respectively. These stations are

376

located along a wide deformation zone that marks the transition between the motion of the NAB

377

and the SOAM plate (Trenkamp et al., 2002).

378

To facilitate the interpretation of the velocity filed and how the deformation is accommodated

379

within the study area, we analyze the velocities relative to the South American plate projected

380

along three profiles perpendicular to the trench (Figure 2) in the north (profile A-A’), center

381

(profile B-B’) and south (profile C-C’) of Ecuador. Figure 5a–c depicts profile-parallel (black dots)

382

and profile-normal (triangles) velocity components at all stations plotted against the distance

383

along the profile.

384

In profile A–A′ (Figure 5a) the profile-parallel components show a progressive velocity decrease

385

towards the interior of the country, with rates from 17 to 2 mm/yr, and an average compression

386

rate of -0.04 (mm/yr)/km. The profile-normal component shows a similar decreasing pattern

387

although, in general, somewhat lower in the stations located in the DMS. The normal component

388

of the SNLR and EPEC stations is slightly lower than expected. The parallel and normal components

389

of this profile show the “escape” of NAB towards the NE, with a clear partition of the velocity

390

vectors.

391

In Profile B–B′ (Figure 5b), profile-parallel and profile- velocities also decrease inwards, although

392

with a lower rate, mainly in the normal component (from 7 mm/yr to -2 mm/yr). The parallel

393

component shows a compressional rate similar to the north profile (-0,04 mm/yr/km). The normal

394

component shows a velocity decrease to the south (profile orientation 105ºN), from the stations

395

located in the MDS with negative normal components. The distribution of velocities in this profile

396

differs from the other profiles, where the parallel component decreases and the normal velocities

397

have a lower rate, with respect to the north profile.

398

In C-C 'profile (Figure 5c), the profile-normal and -parallel components are clearly different from

399

those obtained in the north (A-A') and central (B-B') profiles. The parallel-component has an

400

average value of 6 mm/yr (with a maximum in the MHEC station located on the coast in the Gulf of

401

Guayaquil), and there is no decrease in this component along the profile. However, the profile-

AC C

EP

TE D

M AN U

SC

RI PT

373

14

ACCEPTED MANUSCRIPT

normal component is negative (except for the MHEC station), which reflects a higher southern

403

component at these velocities compared to the previous two profiles, and similar values with an

404

average of -2 mm/yr (with the exception of CLEC with a maximum value of -5 mm/yr). Given their

405

location, magnitude and orientation, these stations may well be part of to the so-called INCA Sliver

406

(Nocquet et al., 2014) or the PERU Block (Sánchez and Drewes, 2016), as we discuss below.

407

Along with the profiles of the velocities, it is show a profile of the topography in each profile. In

408

these, we can appreciate the height difference of the Andes mountain range, which decreases

409

from 6000 m in the north profile to almost 3000 m in the south profile. This also highlights the

410

difference in tectonic behavior between these zones.

411

5.3.

412

In contrast to the displacement data, the strain rate tensor is independent of the reference frame.

413

The strain rate field reveals local strain accumulation rates and their possible connection to

414

seismic hazard potential (Ward, 1994). Therefore, we calculated the velocity gradient tensor based

415

on the horizontal velocity field of REGME station, based on the formulation proposed by Feigl et

416

al. (1990). To this end, a Delaunay triangulation method was used and displacements at triangle

417

vertices were assumed homogeneous. Station EPEC was not taken into account in the

418

triangulation due to its proximity to station QUEM. Subsequently, the velocity gradient was used

419

to calculate strain and rotation rates within the network. Based on these parameters, we also

420

calculated dilatation ( 3) and maximum shear strain rates (43 ) and their directions.

421

The horizontal principal strain rate axes ( 3 and 3 ) and dilatation ( 3)̅ are shown in Figure 6. A

422

convention of positive strain rates indicating extension is used. Henceforth, we will refer to

423

positive dilatation as “extension” and negative dilatation as “shortening”. There is an overall

424

predominance of negative dilatation that we believe must be related to a compressional regimen

425

in the zone. The mean shortening rates are -0.05 ± 0.01 μstrain/yr. The main orientation of the

426

extension and shortening axes are mostly NNE and ESE, respectively, which are congruent with the

427

orientation of the direction of subduction and the main traces of the MDS. Within the NAB, the

428

maximum extension is located on the coast, between latitudes -2° and 0°, with a maximum

429

extension rate of 0.34 ± 0.01 μstrain/yr and a direction that is practically perpendicular to the

430

trench 115 ± 1° (ESE). Nevertheless, from latitude 0° to the north, the shortening rate

AC C

EP

TE D

M AN U

Strain Rate Calculation.

SC

RI PT

402

15

ACCEPTED MANUSCRIPT

predominates with an average value of -0.05 ± 0.01 μstrain/yr and an average orientation of the

432

compression axis of 54 ± 9°. In the Gulf of Guayaquil, the estimated extension rate is 0.08 ± 0.01

433

μstrain/yr and the shortening rate is estimated at -0.03 ± 0.01 μstrain/yr with a 115 ± 3° (ESE)

434

direction for this axis, consistent with the direction of the main tectonic structures in the zone. On

435

the MDS (over the Andes Mountains), the deformation zone between the NAB and the SOAM

436

plate, the shortening rate predominates with an average value of -0.04 ± 0.01 μstrain/yr and an

437

average orientation for this axis of 80 ± 5° (ENE-E). Nevertheless, it may be observed that (Figure

438

6) the orientation of this axis changes from north to south within this zone, from azimuths of ~70°

439

(ENE) to ~115° (ESE). As may be expected, deformation rates are decreasing towards the interior

440

of the continent.

441

The dilatation rate ( 3) is calculated by adding the shortening and extensional rates. Negative

442

dilatation rates predominate in the study area (Figure 6), with a maximum of -0.13 ± 0.01

443

μstrain/yr. The maximum dilation rate is located on the coast, between latitudes -2° and 0°, with a

444

positive dilation rate of 0.34 ± 0.01 μstrain/yr. In the south, dilation rates are positive (~0.05

445

μstrain/yr), although the absence of data in Peru, makes it difficult to confirm and analyze these

446

values.

447

Figure 7 shows maximum shear strain rates (43 ) and their directions. Maximum values of 43 are

448

obtained around the MDS, with the exception of the value obtained in the coastal zone (between

449

latitudes -2° and 0°), with a maximum value of 0.34 ± 0.01 μstrain/yr and a NNW-SSE orientation.

450

With the exception of this value, the predominant orientation of the maximum shear formation is

451

~50° (NE-ENE), which is consistent with the orientation of the main structures of the MDS and with

452

the movement of the North Andean Block.

453

In terms of the rigid body deformation, clockwise rotation predominates (Figure S7 Supplementary

454

Material) along the study area, with a mean value of ~1°/Myr.

455

5.4.

456

To jointly analyze the velocity and strain fields obtained with the seismicity of the study area, we

457

have resort to the focal mechanisms catalog of Global CMT (Ekström et al., 2012). In Figure 8a we

458

have plotted the focal mechanisms for the region classified by rupture mechanism (Álvarez-

459

Gómez, 2014), which allows to highlight the different tectonic processes.

AC C

EP

TE D

M AN U

SC

RI PT

431

Focal Mechanisms slip vectors.

16

ACCEPTED MANUSCRIPT

We have computed also the horizontal slip vector for each focal mechanism (Figure 8b). Both

461

nodal planes of a focal mechanism are physically equivalent, meaning that without further

462

information or assumptions we cannot know which fault plane is responsible for the earthquake.

463

Working with a populated database, the simplest way to decide which fault plane to choose is by

464

applying an Andersonian mechanical assumption. This assumption implies that the responsible

465

planes for reverse ruptures are those with the lowest dip (vertical minimum stress axis), while

466

steeper planes (vertical maximum stress axis) are responsible for normal earthquakes. In the case

467

of strike-slip earthquakes, both planes should be sub-vertical (vertical intermediate stress axis) and

468

further information is needed to determine which plane is responsible. We have decided to

469

represent both sub-orthogonal slip vectors for strike-slip events.

470

The horizontal slip can be obtained by breaking the slip vector over the plane down into horizontal

471

and vertical components. The net slip is calculated based on the seismic moment:

M AN U

SC

RI PT

460

7

6̅ = 9 :8 (6)

472

μ is the rigidity (assumed as 30 GPa for continental crust) and A is the rupture area, obtained from

474

the Blaser et al. (2010) scaling relations. The slip thus obtained must be considered as an average

475

over the plane, although, for the sake of simplicity, its representation on Figure 8b is as points.

476

It seems from Figure 8b that the slip vectors of trench-related seismicity are parallel to the

477

subduction direction of the Nazca Plate under the South American Plate (see Figure 2) and

478

consequently parallel to the vectors of the GNSS network near the trench; suggesting the absence

479

of slip partitioning in front of Ecuador (in contrast with the proposition made by Nocquet et al.,

480

2014). This is congruent with the small obliquity angle (γ) of the subduction in most of the area

481

(McCaffrey, 1994) and the oblique rake of the subduction interface thrust ruptures (~100°). The

482

remarkable amount of thrusting seismicity offshore from the Coastal Range of Ecuador is

483

associated with the subduction of the Carnegie Ridge and the higher subduction interface coupling

484

(Nocquet et al., 2014).

485

The trailing edge of the Dolores-Guayaquil Megashear is located in the Gulf of Guayaquil, where

486

an extensional fault splay is thinning the crust forming the gulf and producing the normal rupture

487

earthquakes shown in Figure 8. The slip vectors of strike-slip earthquakes located in the Gulf of

488

Guayaquil show a northeastward right-lateral motion, coinciding with the northeastern motion of

AC C

EP

TE D

473

17

ACCEPTED MANUSCRIPT

the NAB through the Dolores-Guayaquil Megashear, with the orientation of the subducting

490

Grijalva, Alvarado and Sarmiento ridges (Lonsdale, 2005), and with the shear strain orientation

491

obtained (Figure 7).

492

The seismicity of the MDS can be divided into three groups, as seen in Figure 8. A southern group

493

is located in the Peruvian Andean, thrusting over the stable South American plate on the Amazon

494

Basin. It is composed mostly by reverse fault earthquakes, showing a E-NE convergence. The

495

northern group, in Colombia, shows mainly right-lateral strike-slip events on NE-striking nodal

496

planes. One such event was the MW 6.9 1994 Paez earthquake, associated with a set of NE striking

497

structures in the axis of the Cordillera Central (Irlanda Fault, Paris et al., 2000). Finally, there is a

498

group of reverse and strike-slip events along the Andes Mountains in Ecuador. Most of the

499

thrusting events evince a shallow dipping plane towards the west, with a slip vector towards the

500

east-southeast, but also towards the east-northeast. Towards the south of the group, the slip

501

vectors seem to point more towards the southeast. On the opposite side of the group, in the

502

north, events seem to point towards the northeast (see Figure 8). This geometric configuration

503

seems to indicate a “escaping” disposition of the vectors, which is confirmed by our GNSS results

504

(Figure 2).

SC

M AN U

TE D

505

RI PT

489

6. MODELLING OF THE VELOCITY FIELD.

507

To make a more robust interpretation of the crustal deformation obtained from our GNSS

508

velocities, we inverted the horizontal components of the GNSS velocity vectors using the software

509

code DEFNODE developed by McCaffrey (2002). GNSS velocities are considered from the result of

510

the combination of relative block rotations and elastic deformation due to coupling at the block

511

boundaries. The relative block motions are defined by spherical Earth-angular velocity vectors

512

(Euler rotation poles and rates) while the interseismic deformation is modeled as backs slip on the

513

faults that separate blocks (Okada, 1985; Savage, 1983). The faults bounding the finite blocks are

514

defined in 3D by a series of nodes along the fault planes. Fault locking is parametrized at each

515

node by a coupling factor (φ) that represents the fractional part of the relative block motion,

516

which is not accommodated by steady, aseismic slip. Values for φ ranges between 0 (no coupling)

517

and 1 (full coupling). Block-angular velocities and coupling factors φ can be inverted by minimizing

AC C

EP

506

18

ACCEPTED MANUSCRIPT

the misfit between observations (GNSS velocities) and predicted data using a simulated annealing

519

method.

520

We constrained our model by applying the horizontal GNSS velocities and their associated

521

uncertainties listed in Table 1. We defined a geometry for our model in four blocks: SOAM (South

522

American Plate), NZ (Nazca Plate), NAB (North Andean Block) and INCA (Inca Sliver) (Figure 9). The

523

limit between SOAM and NAB was simplified by a 70° west-dipping fault plane (Figure 9) that

524

represents the MDS, with 35 km depth and fully coupled up to a 20 km depth. We did not model

525

the limit between INCA and SOAM. The geometry and location of the first 100 km depth of the

526

Colombian-Ecuador subduction zone is constrained using the geometry of Hayes et al. (2012).

527

The angular velocity vectors relative to the SOAM plate are from Kendrick et al. (2003) for NZ and

528

from Nocquet et al. (2014) for NAB and INCA.

529

We inverted only for the coupling along the Ecuadorian Subduction Zone. We have assumed that

530

the coupling can occur between the trench and 40 km depth (e.g. Chlieh et al., 2014), with no

531

coupling near the trench.

532

A series of resolution tests (Figure S8 Supplementary Material) indicate that our data distribution

533

allows us to constrain variations of coupling on the subduction interface between latitudes 4°S and

534

1.5°N and with patches measuring 50 km along-strike and 20 km along-dip.

535

The best-fitting model (Figure 9) shows a main coupling asperity between Muisne (0.5°N) and

536

Bahía de Caráquez (0.5°S), north of Manta, between 10 and 40 km in depth, coinciding with the

537

epicenter and rupture area of the 2016 Pedernales earthquake.

538

To the north of this main asperity, the model shows two other small asperities; east of Esmeraldas

539

between 30 and 40 km in depth, and north of Esmeraldas at a depth of 20 km, which coincide with

540

the epicenters of the 1979, and 1906 and 1958 earthquakes, respectively (Figures 1 and 9).

541

However, as indicated in the resolution tests, our distribution of data does not allow us to

542

constrain coupling variations beyond 1.5°N.

543

In line with the results obtained by Chlieh et al., 2014, our model shows another isolated, small

544

and highly coupled patch between Manta and the peninsula of Santa Elena (between 1.5° and 2°S)

545

at a depth of 20 km. However, we do not have GNSS data in the La Plata Island, which could

546

constrain better our model.

AC C

EP

TE D

M AN U

SC

RI PT

518

19

ACCEPTED MANUSCRIPT

North of this asperity, there is a large creeping corridor that lies immediately south of the trail

548

along the shallow axis of the Carnegie Ridge Track (see fig. 4 in Chlieh et al., 2014).

549

Our model shows no interseismic coupling at the interface to the south of Ecuador, from the

550

peninsula of Santa Elena (south of latitude 2°S) to the south of Ecuador.

551

Our results are in agreement with the results obtained by Chlieh et al. (2014) despite the different

552

time span and spatial distribution of the GNSS stations used in both studies.

RI PT

547

553

7. DISCUSSION

555

The interseismic velocity field in most subduction plate boundaries is usually interpreted in terms

556

of two major processes: strain accumulation due to coupling on the shallow (less than 40 km

557

depth) subduction zone, and trench-parallel translation of the forearc (e.g. McCaffrey et al., 2000).

558

In Ecuador, there are three main models that explain the north-northeast directed forearc motion

559

of the NAB: (1) the strain partitioning caused by oblique convergence and mechanical coupling on

560

the Colombia-Ecuador subduction zone, (2) collision of Carnegie Ridge and resulting tectonic

561

escape, or (3) variable coupling along the Colombia-Ecuador trench. (e.g. Lonsdale, 1978;

562

Pennington, 1981; Ego et al., 1996; Trenkamp et al., 2002; Witt et al., 2006; Witt and Bourgois,

563

2009; Nocquet et al., 2014; Chliet et al., 2014).

564

In this section, we describe the main characteristics of our GNSS velocity field and strain rate field

565

and discuss the significance and interpretation of these results in terms of the tectonics of the

566

region.

567

7.1.

568

The dominant (eastward) direction of motion in the convergence of the Nazca and South American

569

plates (Figure 1) is the main driving force responsible for ongoing crustal deformation in Ecuador

570

and is observed as the most prominent feature in our velocity field.

571

One of the stations analyzed is located on the Nazca plate: station IGS on the Galapagos Islands

572

(GLPS) (Figure 2), with a horizontal velocity of 55.8 ± 0.1 mm/yr and a practically eastward

573

direction (87.4°) (Table 1 and Figure 2), congruent with the eastern movement of the Nazca plate

574

relative to fixed South America (Trenkamp et al., 2002; Sánchez and Drewes, 2016).

EP

TE D

M AN U

SC

554

AC C

Strain and velocity Fields.

20

ACCEPTED MANUSCRIPT

The rest of the observed velocities, located on the overriding plate, reach their highest values

576

along the coast, mainly in the north, with a maximum rate of ~23 mm/yr (which is almost a 40% of

577

the rate motion between Nazca and South American plates) at ESMR station. We observed a

578

notable velocity decrease eastward, reaching values close to zero (statistically insignificant, at a

579

95% confidence level) inland, east of the MDS. The change in the velocity orientations is especially

580

noteworthy, with an observable rotation from NE-ENE to SE from north to south, which may

581

reflect the coupling difference in the subduction interface and the escape towards the SE away

582

from the Carnegie Ridge subduction of the Inca Sliver.

583

The movement registered by station ESMR, located in a highly coupled zone near the epicenter of

584

the Mw 7.8 earthquake that occurred in April 2016 (Figure 9), stands out for its clearly higher

585

easterly component (22.3 ± 0.1 mm / yr.). This is related to a significant accumulation of

586

deformation prior to the April 2016 earthquake. It would be important to unravel the present-day

587

strain distribution, after the 2016 earthquake, to confirm this conclusion and to draw further

588

implications for the seismic hazard in the area.

589

The stations located in the deformation zone of the boundary of the NAB display a clockwise

590

rotation in their velocities from north to south (Figures 2 and 5). The movements recorded by

591

these stations clearly indicate a right-lateral motion distributed along the fault zone, which is

592

congruent with the directions obtained for maximum shear deformation (Figure 7). These stations,

593

located in the transition zone, show a greater eastern component in their movement, which could

594

be due to the ongoing shortening deformation associated with the MDS.

595

In order to determine the interseismic deformation accumulated in the study area, we have

596

removed of our velocity field the rotational part of each block used in section 6 (Figure 10a).

597

Figures 10b-d show these velocities projected on three profiles (Figure 2) perpendicular to the

598

trench. The maximum variation of the velocity component normal to the profile (parallel to the

599

main MDS faults) in northern Ecuador (red triangles in Figure 10b) indicates a right-lateral motion

600

of 7.6 ± 0.5 mm/yr (between station ESMR, on the coast, and station AUCA, in the interior). The

601

maximum variation of the velocity component parallel to the profile (perpendicular to the main

602

MDS faults) indicates a shortening of 9.6 mm/yr ± 0.5 mm/yr (black circles in Figure 10b), which is

603

consistent with the results of the strain analysis.

AC C

EP

TE D

M AN U

SC

RI PT

575

21

ACCEPTED MANUSCRIPT

In the Central Ecuador, the variation of the normal component decreases (Figure 10c), with an

605

estimated right-lateral motion of 5.3 ± 0.4 mm/yr. The parallel component also decreases with a

606

shortening of 6.9 mm/yr. These results are consistent with the distinct behavior of subduction in

607

Ecuador, with no interseismic coupling in southern Ecuador, and increased coupling towards the

608

north in the zone where megathrust earthquakes have occurred over the last century (stars in

609

Figures 1 and 6). The southern part of the Carnegie Range marks the limit between the two zones.

610

In the south of Ecuador, the velocities projected in the profile (Figure 10d) display a different

611

behavior, with values of the velocity components parallel and normal to the profile practically

612

zero, with the exception of MHEC and NJEC stations, which are very close to the main trace of

613

MDS.

614

Although this study does not address vertical deformation rates, the results obtained for vertical

615

velocities (see Figure S9 Supplementary Material) indicate the uplift of stations located in the

616

transition zone between the NAB and the SOAM (stations located over the Andes Range), which is

617

consistent with the existence of transpression, which would cause elevation in this zone. Stations

618

located in the northern half of the coast of Ecuador, closer to the trench (ESMR, PTEC) register

619

uplift of ~3 ± 0.2 mm/yr, which is consistent with the coupling of the subduction zone in this area

620

(Figure 9); however, station SEEC registers a slight subsidence of -1.2 ± 0.6 mm/yr, consistent with

621

a weaker coupling in this zone (as calculated in the previous section) and agrees with the results of

622

Chlieh et al., 2014.

623

The orientation and magnitude of the principal strain rate axes (Figures 6 and 7) and rotation rates

624

(Figure S7 Supplementary Material) obtained by inverting GNSS data are in agreement with other

625

regional tectonic and seismic studies (e.g. Trenkamp et al., 2002; White et al., 2003) and

626

consistent with the active regime in the area.

627

Principal strain and dilatation rates (Figure 6) clearly indicate a difference in behavior between the

628

half north and south of the coast, the central zone corresponding to the Andes mountain range

629

and the southern part of Ecuador, in which the deformation values decrease and the N-S

630

extension predominates.

631

In the northern part (from Bahía de Caráquez towards Colombia), the values registered were

632

│ 3 │<│ 3 │ with a negative dilation rate and a predominance of the shortening component, with a

AC C

EP

TE D

M AN U

SC

RI PT

604

22

ACCEPTED MANUSCRIPT

maximum of -0.09 ± 0.01 μstrain/yr. This is consistent with the coupling values obtained for this

634

zone.

635

In the southern part (between Bahía de Caráquez and the Gulf of Guayaquil), the values registered

636

were │ 3 │>│ 3 │ with a positive dilation rate, reaching its peak value of 0.34 ± 0.01 μstrain/yr near

637

Isla de la Plata. These values could be related to the 2010 SSE of La Plata Island (Chlieh et al., 2014)

638

although, as shown in the coupling model (Figure 9), this zone is currently coupled.

639

In the MDS deformation zone that forms the boundary between the NAB and SOAM, we

640

obtain│ 3 │<│ 3 │, with a predominance of negative dilation rates and shortening behavior.

641

In the Gulf of Guayaquil, the extensional component predominates, │ 3 │>│ 3 │, with a NNE (~25°)

642

direction, which can be explained as trailing edge deformation of the NAB on its motion towards

643

the NE.

644

7.2.

645

Our results show a difference in the subduction dynamics in Ecuador due to the heterogeneity of

646

the coupling pattern at depths of up to 40 km.

647

Our model suggests a main asperity of 200 km along-strike (between 0.5°N and 0.5°S) and

648

between 10 and 40 km depth, that coincides with the epicenter and spatial slip distribution (black

649

star and discontinues lines in Figure 9) of the 2016 Pedernales earthquake (Nocquet et al., 2016).

650

Our GNSS data cover the 2008-2014 period and are clearly registering the interseismic strain

651

accumulation before the 2016 earthquake.

652

Our study indicates the absence of coupling in southern Ecuador (southwards the southern border

653

of the Carnegie Ridge) and suggests that coupling increases towards the north, in the zone where

654

megathrust earthquakes have occurred over the last century (stars in Figure 6). Our results show

655

that the epicenter of the 1998 earthquake marks the boundary between the two zones that

656

exhibit different behavior, practically coinciding with the southern border of the Carnegie Ridge.

657

In line with the results obtained by Chlieh et al. (2014), our model suggests a small asperity south

658

of Manta at a depth of 20 km. This is related to the La Plata area (latitude, 1.3°S; longitude,

659

81.1°W, Figure 1), which is the best-documented region along the Ecuador subduction zone for

660

frequent SSEs, seismic swarms, and/or repeating earthquakes, which have been recorded in 1977,

661

1998, 2002, 2005, 2010, 2013 and 2016 (Holtkamp et al., 2011; Vallée et al., 2013; Segovia et al.,

M AN U

SC

RI PT

633

AC C

EP

TE D

Coupling Along the Ecuadorian Seismogenic Subduction Zone.

23

ACCEPTED MANUSCRIPT

2015; Rolandone et al., 2018). Nevertheless, analyzing the time series of the stations closest (SEEC,

663

PJEC and PTEC) (see Figure S2 in Material Supplementary) to this zone, no significant displacement

664

is found that could be related to the 2010 and 2013 SSE. However, this is an important zone in

665

understanding the seismic hazard in that region.

666

7.3.

667

The distribution of strain and velocity is consistent with the transpressive right-lateral slip of the

668

westward-dipping faults of the MDS and the NNE movement of the NAB relative to the SOAM

669

plate, as suggested by previous studies (e.g. Pennington, 1981, Trenkamp et al., 2002).

670

Our results show a different behavior between the north, central and south of Ecuador, which

671

seems to be controlled by the subduction of the Carnegie Ridge and the degree of interseismic

672

coupling of the subduction interface.

673

In order to check the possible slip partitioning at the trench we have developed an obliquity

674

analysis. In Figure 11a the different parameter azimuths are shown. We have followed the

675

McCaffrey (1993) notation shown in Figure 11b. The blue shaded area, and the thick blue line,

676

represents the azimuth of the subduction interface normal on the upper 10 km from the SLAB 1.0

677

model (Hayes et al., 2012). The plate motion vector azimuth with respect to a fixed SOAM plate on

678

each one of the previous subduction interface points are represented as a red line. Finally, the slip

679

direction from the reverse earthquakes of the Global CMT focal mechanisms on the first 40 km of

680

the subduction zone are shown as orange circles. A trend line using a gaussian filter has been

681

obtained to be used in the following calculations.

682

As is clearly shown on Figure 11a the azimuth of the seismicity slip vectors presents a variability

683

range of ~40°, overlapping the plate motion azimuths, explaining the previously done visual

684

interpretation of the coincidence of the earthquake slip vectors and the plate motion azimuth.

685

Nevertheless, when we compute the general trend of the earthquake slip vectors they form as

686

average an angle of 5 – 10° with the plate motion vector. In Figure 11b the angles between the

687

earthquake slip vectors and the trench normal (ψ), between the plate motion vector and the

688

trench normal (γ) and the difference between both (δ) are shown. This angular difference is the

689

responsible of the forearc sliver motion in a strain partitioning context. If δ is 0 it means that the

RI PT

662

AC C

EP

TE D

M AN U

SC

Tectonic Interpretation

24

ACCEPTED MANUSCRIPT

convergence is absorbed entirely as oblique thrusting in the subduction interface. The higher the

691

difference the higher the degree of partitioning.

692

The degree of partitioning can be defined by kinematic (pk) or seismic (ps) criteria (Chemenda et

693

al., 2000). The kinematic paritioning is defined as pk = vs/(V sin γ), where vs is the trench parallel

694

velocity of forearc sliver translation, V is the plate motion velocity and γ the plate motion vector

695

obliquity. As is shown in this work the maximum strike-slip rate of the MDS faults is 7.6 mm/yr (vs)

696

and using 55 mm/yr as plate convergence rate (V) the kinematic partitioning (pk) varies from 0.40

697

in southern Ecuador (γ = 20º), 0.28 in front of the Carnegie Ridge subduction (γ = 30º), to 0.18 in

698

northern Ecuador (γ = 50º). The seismic partitioning is defined as ps = 1 – (γ – δ)/(γ). The ps value

699

has been computed and plotted in Figure 11c. The maximum values of ps are obtained in southern

700

Ecuador (ps ~ 0.5 – 0.6) while in central and northern Ecuador the ps values are between 0.3 and

701

0.15. Both partitioning criteria values, pk and ps, are consistent and show a low degree of

702

partitioning.

703

The capacity of a subduction to generate slip partitioning is directly related to the interplate

704

coupling, the obliquity and the presence of a weak zone in the upper plate, usually a volcanic arc

705

(Jarrard, 1986; Chemenda et al., 2000; Philippon and Corti, 2016), although this point is discussed

706

by McCaffrey et al. (2000). In Ecuador the volcanic arc acts as the weak zone to promote the

707

formation of the trench parallel structures allowing the NE motion of the NAB. The obliquity in

708

general is lower than 30°, and only around the latitude 1°N the obliquity reach values of 40° – 50°.

709

The value of 30° (or even higher) is sometimes considered as a critical limit and below it the

710

subduction is unable to produce a partitioning by itself (Fitch, 1972; Philippon and Corti, 2016) or

711

maybe the degree of obliquity is directly related to the amount of partitioning, the higher the

712

obliquity the higher the partitioning (Jarrard, 1986). From McCaffrey et al. (2000) we can also infer

713

that when the obliquity and subduction interface coupling are high there is no need for a weak

714

zone to develop the strike-slip faults that allow the forearc sliver motion.

715

The coupling must be sufficiently high to enhance the efficient stress transfer from lower to upper

716

plates (Chemenda et al., 2000; Philippon and Corti, 2016). In Figure 11c the degree of coupling

717

obtained in this work is represented. The dots show the discrete values on the interface and the

718

line is the trend line showing the maximum coupling. The main patches described on the coupling

AC C

EP

TE D

M AN U

SC

RI PT

690

25

ACCEPTED MANUSCRIPT

can be clearly seen. Taking into account this coupling variation, and using the geometrical relations

720

shown in Figure 11b, we have estimated the velocities on the upper plate from the plate motion

721

vectors: vn is the trench-normal component of the plate subduction vector, vp is the trench-

722

parallel component, and vs is the predicted forearc sliver (NAB) velocity as described by McCaffrey

723

(1993): vs = V (sin γ – cos γ tan ψ).

724

The maximum contribution to the forearc sliver motion is obtained between latitudes -1 and 1;

725

and also in a patch around latitude -1.5 (Figure 11d). These areas correspond to the Carnegie

726

Ridge subducting on a N80°E direction, although the partitioning in this area is low, its impact on

727

the upper plate seems to be directly related to the NAB motion towards the NE. The maximum

728

NAB velocity predicted by this geometrical method, assuming a rigid forearc sliver, is 14 mm/yr in

729

northern Ecuador. This value is several mm/yr higher than the NAB motion obtained with GNSS in

730

this work, which is reasonable as the block is internally deformed as shown in section 5.3 and

731

Figure 6.

732

Earthquake slip vectors suggest that the NAB is being pushed by the Carnegie Ridge subduction.

733

This leads to a higher subduction interface coupling in this area. The higher amount of seismicity in

734

the trench offshore the Coastal Range and the radial disposition of slip vectors in the Andean

735

Range supports this model.

736

Furthermore, stations located in southern Ecuador present a negative north velocity component,

737

indicating that these stations present a clearly distinct behavior, outside the NAB, and could

738

perhaps be located within the Inca Block in Peru, as suggested by Nocquet et al., 2014.

739

Based on our results, we suggest that the main driving force responsible for ongoing crustal

740

deformation in Ecuador is the convergence between the Nazca and South American plates with

741

the variable coupling pattern and the collision of the Carnegie Ridge. This produces a drastic

742

change in the velocity pattern obtained from north to south. This pattern is similar to that

743

obtained in the Pacific coast of Central America with the collision of the Cocos Ridge in Costa Rica

744

(Kobayashi et al., 2014).

AC C

EP

TE D

M AN U

SC

RI PT

719

745 746

8. CONCLUSIONS

26

ACCEPTED MANUSCRIPT

We have presented a new GNSS velocity fields for Ecuador based on the time series analysis for 33

748

REGME stations for the period between 2008 and 2014. Our analysis of these data confirms and

749

quantifies current tectonic activity in Ecuador.

750

The velocity and strain distribution clearly correspond to the transpressive right-lateral slip of the

751

westward-dipping faults of the MDS and the NNE movement of the NAB relative to the SOAM

752

plate.

753

The distributions of our deformation rate and velocity field indicate a differentiated tectonic

754

behavior between northern, central and southern Ecuador. In northern Ecuador, there is an

755

estimated right-lateral motion of 7.6 ± 0.5 mm/yr, consistent with the NNE movement of the NAB

756

relative to the South American plate. In central Ecuador, the right-lateral motion decreases to 5.3

757

± 0.4 mm/yr. In southern Ecuador (from the Guayaquil Gulf to Peru) there is no strain

758

accumulation, the GNSS velocities decrease and turn to the south. This zone belongs to the so-

759

called Inca or Peru sliver.

760

These results are consistent with the distinct behavior of subduction in Ecuador, with an absence

761

of coupling in southern Ecuador, and increased coupling towards the north in the zone where

762

megathrust earthquakes have occurred over the last century. The southern part of the Carnegie

763

Ridge marks the limit between the two zones.

764

We suggest that the main driving force responsible for ongoing crustal deformation in Ecuador is

765

the convergence between the Nazca and South American plates with the variable coupling pattern

766

and the collision of the Carnegie Ridge. This produces the significant change in the velocity pattern

767

obtained from north to south.

EP

TE D

M AN U

SC

RI PT

747

AC C

768 769

ACKNOWLEDGEMENT

770

We would like to thank the Military Geographic Institute of Ecuador (IGM), for providing access to

771

the Ecuador -REGME GNSS Continuous Monitoring Network data. The figures were produces using

772

the Generic Mapping Tools (GMT) software.

773 774

REFERENCIAS

27

ACCEPTED MANUSCRIPT

775 776 777 778

Agnew, D. C. (1992). The time-domain behavior of power-law noises. Geophysical research letters, 19(4), 333-336. Álvarez-Gómez, J. A. (2014). FMC: a one-liner Python program to manage, classify and plot focal mechanisms. In EGU General Assembly Conference Abstracts (Vol. 16, p. 10887). Amiri-Simkooei, A. R., Tiberius, C. C. J. M. and Teunissen, S. P. (2007). Assessment of noise in GPS

780

coordinate time series: methodology and results. Journal of Geophysical Research: Solid

781

Earth, 112(B7). doi: 10.1029/2006JB004913.

RI PT

779

Beauval C., Yepes H., Bakun W., Egred J., Alvarado A., Singaucho, J. (2010) Locations and

783

magnitudes of historical earthquakes in the Sierra of Ecuador (1587–1996). Geophysical

784

Journal International 181: 1613-1633. doi: 10.1111/j.1365-246X.2010.04569.x.

SC

782

Bethoux, N., Segovia, M., Alvarez, V., Collot, J. Y., Charvis, P., Gailler, A., & Monfret, T. (2011).

786

Seismological study of the central Ecuadorian margin: evidence of upper plate

787

deformation. Journal

788

10.1016/j.jsames.2010.08.001.

790

South

American

Earth

Sciences, 31(1),

139-152.doi:

Bird, P. (2003). An updated digital model of plate boundaries. Geochemistry, Geophysics, Geosystems, 4(3).

TE D

789

of

M AN U

785

791

Blaser, L., Krüger, F., Ohrnberger, M. and Scherbaum, F. (2010). Scaling Relations of Earthquake

792

Source Parameter Estimates with Special Focus on Subduction Environment. Bulletin of the

793

Seismological Society of America, 100(6): 2914–2926.

795

Geophysical Research: Solid Earth, 107(B7). doi: 10.1029/2001JB000570. Bos, M. S., Fernandes, R. M. S., Williams, S. D. P. and Bastos, L. (2008). Fast error analysis of

799

Chemenda, A., Lallemand, S., & Bokun, A. (2000). Strain partitioning and interplate friction in

800

AC C

796

Blewitt, G. and Lavallée, D. (2002). Effect of annual signals on geodetic velocity. Journal of

EP

794

oblique subduction zones: Constraints provided by experimental modeling. Journal of

801

Geophysical Research: Solid Earth, 105(B3), 5567-5581.

797 798

continuous GPS observations. Journal of Geodesy, 82(3), 157-166. doi: 10.1007/s00190-0070165-x.

28

ACCEPTED MANUSCRIPT

802

Chlieh, M., Mothes, P. A., Nocquet, J. M., Jarrin, P., Charvis, P., Cisneros, D., ... & Vallée, M. (2014).

803

Distribution of discrete seismic asperities and aseismic slip along the Ecuadorian megathrust.

804

Earth and Planetary Science Letters, 400, 292-301. doi:10.1016/j.epsl.2014.05.027. Chunga, K., Michetti, A., Pazmiño, N., Martillo, C., Romero, A., Quiñonez, M., Gruppo di Geologia

806

Ambientale (2009). Estimación de máximos niveles de sismicidad para el litoral ecuatoriano a

807

través de la integración de datos geológicos, sismológicos y sismotectónicos. International

808

Journal: Oro y Petróleo 19: 46-57.

RI PT

805

Chunga, K., M.F. Quiñónez, Huaman, F., Besenzon, D., Mulas, M., Garcés, D., Larreta, E., Gorshkov,

810

A., Michetti, A.M. (2016). Geología de Terremotos y Tsunamis. Instituto Panamericano de

811

Geografía e Historia. Sección Nacional de Ecuador.

SC

809

Cisneros, D., and Nocquet, J. M. (2011). Campo de velocidades del Ecuador, obtenido a través de

813

mediciones de campañas GPS de los últimos 15 años y medidas de une red GPS

814

permanente. Revista Geoespacial, 9, 30-49.

M AN U

812

Collot J. Y., Michaud F., Alvarado A., Marcaillou B., Sosson M., Ratzov G., Pazmino A. (2009). Visión

816

general de la morfología submarina del margen convergente de Ecuador-Sur de Colombia:

817

Implicaciones sobre la transferencia de masa y la edad de la subducción de la Cordillera de

818

Carnegie. En: Geología y Geofísica Marina y Terrestre del Ecuador desde la Costa Continental

819

hasta las Islas Galápagos. Editores: Collot J. Y., Sallares V., Pazmiño N. Impreso: Argudo &

820

Asociados, Guayaquil-Ecuador. pp. 47-74. ISBN-978-9978-92-737-3.

823 824

EP

822

Dach, R., Hugentobler, U., Fridez, P. & Meindl, M., (2007). Bernese GPS Software Version 5.0. Astronomical Institute. University of Berne. Dow, J.M., Neilan, R.E. y Rizos, C. (2009) The International GNSS Service in a changing landscape of

AC C

821

TE D

815

Global Navigation Satellite Systems. Journal of Geodesy, 83:191–198.

825

Drewes, H., & Heidbach, O. (2012). The 2009 horizontal velocity field for South America and the

826

Caribbean. In Geodesy for Planet Earth (pp. 657-664). Springer Berlin Heidelberg. doi

827 828 829

10.1007/978-3-642-20338-1_81.

Dumont, J. F., Santana, E., & Vilema, W. (2005). Morphologic evidence of active motion of the Zambapala Fault, Gulf of Guayaquil (Ecuador). Geomorphology, 65(3-4), 223-239.

29

ACCEPTED MANUSCRIPT

830 831

Egbue,

O.,

and

Kellogg,

J.

(2010).

Pleistocene

to

present

North

Andean

“escape”. Tectonophysics, 489(1), 248-257. doi: 10.1016/j.tecto.2010.04.021. Ego, F., Sébrier, M., Lavenu, A., Yepes, H., & Egues, A. (1996). Quaternary state of stress in the

833

Northern Andes and the restraining bend model for the Ecuadorian Andes. Tectonophysics,

834

259(1-3), 101-116.

RI PT

832

835

Egüez, A., Alvarado, A., & Yepes, H. (2003). Mapa de fallas y pliegues cuaternarias de Ecuador y

836

regiones oceánicas adyacentes. US Geological Survey–Escuela Politécnica Nacional,

837

Programa Internacional de la Litosfera, Grupo de Trabajo II-2. OFR, 03-289.

Ekström, G., Nettles, M., and Dziewoński, A.M. (2012). The global CMT project 2004-2010:

839

Centroid-moment tensors for 13,017 earthquakes. Physics of the Earth and Planetary

840

Interiors, 200-201: 1-9.

842

M AN U

841

SC

838

Espinoza, J. (1992). Terremotos tsunamigénicos en el Ecuador. Instituto Oceanográfico de la Armada INOCAR, Guayaquil, Ecuador.

Feigl, K. L., King, R. W., & Jordan, T. H. (1990). Geodetic measurement of tectonic deformation in

844

the Santa Maria fold and thrust belt, California. Journal of Geophysical Research: Solid Earth,

845

95(B3), 2679-2699. doi: 10.1029/JB095iB03p02679.

848 849 850 851

southeast Asia and the western Pacific. Journal of Geophysical research, 77(23), 4432-4460. García, L., (2016). Análisis de series temporales en estaciones permanentes GPS. Tesis Doctoral. Universidad Complutense de Madrid. Facultad de Ciencias Matemáticas. Madrid, España.

EP

847

Fitch, T. J. (1972). Plate convergence, transcurrent faults, and internal deformation adjacent to

Goudarzi, M.A., Cocard, M., Santerre, R., 2014. EPC: Matlab software to estimate Euler pole parameters. GPS Solution, 18, 153–162.

AC C

846

TE D

843

852

Goyes P. (2009). Fondos Marinos de Soberanía y Jurisdicción del Ecuador. En: Geología y Geofísica

853

Marina y Terrestre del Ecuador desde la Costa Continental hasta las Islas Galápagos.

854 855

Editores: Collot J. Y., Sallares V., Pazmiño N. Impreso: Argudo & Asociados, GuayaquilEcuador. pp. 131-149. ISBN: 978-9978-92-737-3.

856

Guillier B., Chatelain J. L., Jaillard E., Yepes H., Poupinet G., Fels J. F. (2001). Seismological evidence

857

on the geometry of the Orogenic System in central-northern Ecuador (South America).

858

Geophysical Research Letters. 28 (19): 3749-3752.

30

ACCEPTED MANUSCRIPT

859

Gutscher M A, Malavielli J, Lallemand S, Collot J (1999) Tectonic segmentation of the North

860

Andean margin: impact of the Carnegie Ridge collision. Earth and Planetary Science Letters

861

168: 255-270. doi: 10.1016/S0012-821X(99)00060-6. Hackl, M., Malservisi, R., Hugentobler, U., and Wonnacott, R. (2011). Estimation of velocity

863

uncertainties from GPS time series: Examples from the analysis of the South African TrigNet

864

network. Journal

865

doi: 10.1029/2010JB008142.

of

Geophysical

Research:

RI PT

862

Solid

Earth, 116(B11).

Hayes G. P., Wald D. J., Johnson R. L. (2012). Slab1. 0: A three-dimensional model of global

867

subduction zone geometries. Journal of Geophysical Research. 117: 1-15, doi:

868

10.1029/2011JB008524.

873 874 875 876 877 878 879 880

M AN U

872

Holtkamp, S. G., M. E. Pritchard, R. B. Lohman (2011). Earthquake swarms in South America. Geophys. J. Int. 187, 128–146.

Jarrard, R. D. (1986). Relations among subduction parameters. Reviews of Geophysics, 24(2), 217– 284.

TE D

871

1404–1420, 1977.

Johnson, H. O., and Agnew, D. C. (1995). Monument motion and measurements of crustal velocities. Geophysical Research Letters, 22(21), 2905-2908. Kanamori, H. and McNally, K. C. (1982). Variable rupture mode of the subduction zone along the Ecuador-Colombia coast. Bulletin of the Seismological Society of America, 72(4), 1241-1253.

EP

870

Hey, R.: Tectonic evolution of the Cocos-Nazca spreading center, Geol. Soc. Am. Bull., 88(10),

Kellogg, J. N., & Bonini, W. E. (1982). Subduction of the Caribbean plate and basement uplifts in the overriding South American plate. Tectonics, 1(3), 251-276.

AC C

869

SC

866

881

Kendrick, E., Bevis, M., Smalley Jr, R., Brooks, B., Vargas, R. B., Laurıa, E., & Fortes, L. P. S. (2003).

882

The Nazca–South America Euler vector and its rate of change. Journal of South American

883

Earth Sciences, 16(2), 125-131.

884

Kobayashi, D., LaFemina, P., Geirsson, H., Chichaco, E., Abrego, A. A., Mora, H., & Camacho, E.

885

(2014). Kinematics of the western Caribbean: Collision of the Cocos Ridge and upper plate

886

deformation. Geochemistry, Geophysics, Geosystems, 15(5), 1671-1683.

31

ACCEPTED MANUSCRIPT

887

Kouba, J. (1996). SINEX–Solution-Independent Exchange Format Version 1.00. In Proceedings of

888

the IGS Analysis Center Workshop, Silver Spring, Maryland, USA, IGS Central Bureau, JPL,

889

Pasadena, California, USA.

893 894 895 896 897 898 899 900 901

RI PT

892

space science, 39(2), 447-462.

Lonsdale, P. (2005). Creation of the Cocos and Nazca plates by fission of the Farallon plate. Tectonophysics, 404(3–4), 237–264.

Lyard, F., Lefevre, F., Letellier, T., & Francis, O. (2006). Modelling the global ocean tides: modern insights from FES2004. Ocean Dynamics, 56(5-6), 394-415. doi: 10.1007/s10236-006-0086-x.

SC

891

Lomb, N. R. (1976). Least-squares frequency analysis of unequally spaced data. Astrophysics and

Mann, P., & Burke, K. (1984). Neotectonics of the Caribbean. Reviews of Geophysics, 22(4), 309362.

M AN U

890

Mao, A., Harrison, C. G., & Dixon, T. H. (1999). Noise in GPS coordinate time series. Journal of Geophysical Research: Solid Earth, 104(B2), 2797-2816. doi: 10.1029/1998JB900033. McCaffrey, R. (1994). Global variability in subduction thrust zone-forearc systems. Pure and Applied Geophysics, 142(1), 173–224.

McCaffrey, R. (2002). Crustal block rotations and plate coupling. Plate boundary zones, 101-122.

903

McCaffrey, R. (1993). On the role of the upper plate in great subduction zone earthquakes. Journal

904

TE D

902

of Geophysical Research, 98(B7), 11953–11966. McCaffrey, R., Zwick, P. C., Bock, Y., Prawirodirdjo, L., Genrich, J. F., Stevens, C. W., ... & Subarya, C.

906

(2000). Strain partitioning during oblique plate convergence in northern Sumatra: Geodetic

907

and seismologic constraints and numerical modeling. Journal of Geophysical Research: Solid

908

Earth, 105(B12), 28363-28376.

AC C

EP

905

909

Mendoza, C. and Dewey J.W. (1984). Seismicity associated with the great Colombia-Ecuador

910

earthquakes of 1942, 1958, and 1979: Implications for barrier models of earthquake rupture.

911

Bulletin of the Seismological Society of America. 74 (2): 577–593. Retrieved 2010-03-07.

912

Michaud F., Pazmiño N., Collot J. Y. (2009). El karst submarino de mega depresiones circulares de

913

la Cordillera de Carnegie (Ecuador): Posible Origen por disolución submarina. En: Geología y

914

Geofísica Marina y Terrestre Del Ecuador desde la Costa Continental hasta las Islas

32

ACCEPTED MANUSCRIPT

915

Galápagos. Editores: Collot J. Y., Sallares V., Pazmiño N. Impreso: Argudo & Asociados,

916

Guayaquil-Ecuador. pp. 29-45. ISBN-978-9978-92-737-3. Niell, A. E. (1996). Global mapping functions for the atmosphere delay at radio

918

wavelengths. Journal

of

919

doi: 10.1029/95JB03048.

Geophysical

Research:

Solid

Earth, 101(B2),

3227-3246.

RI PT

917

920

Nocquet, J., Mothes, P., Alvarado, A. (2009) Geodesia, geodinámica y ciclo sísmico en Ecuador. In:

921

Collot J, Sallares V, Pazmiño N (eds) Geología y Geofísica Marina y Terrestre del Ecuador, 1st

922

edn. Argudo & Asociados, Guayaquil-Ecuador, pp. 83-94

Nocquet, J. M., Villegas-Lanza, J. C., Chlieh, M., Mothes, P. A., Rolandone, F., Jarrin, P., ... & Martin,

924

X. (2014). Motion of continental slivers and creeping subduction in the northern Andes. Nat

925

Geosci, 7(4), 287-291.

M AN U

SC

923

926

Nocquet, J. M., Jarrin, P., Vallée, M., Mothes, P. A., Grandin, R., Rolandone, F., ... & Régnier, M.

927

(2017). Supercycle at the Ecuadorian subduction zone revealed after the 2016 Pedernales

928

earthquake. Nature Geoscience, 10(2), 145.

930 931

Okada, Y. (1985). Surface deformation due to shear and tensile faults in a half-space. Bulletin of the seismological society of America, 75(4), 1135-1154.

TE D

929

Paris, G., Machette, M. N., Dart, R. L., & Haller, K. M. (2000). Map and database of Quaternary

932

faults and folds in Colombia and its offshore regions, 61. U.S. Geological Survey Open-File

933

Report 00-0284

Pedoja K., Dumont J. F., Ortlieb L. (2009). Levantamiento Cuaternario costero del Arco de Talara

935

(Ecuador y norte del Perú): cuantificaciones con las secuencias de terrazas marinas. En:

936

Geología y Geofísica Marina y Terrestre del Ecuador desde la Costa Continental hasta las

938 939 940

AC C

937

EP

934

Islas Galápagos. Editores: Collot J. Y., Sallares V., Pazmiño N. Impreso: Argudo & Asociados, Guayaquil-Ecuador. pp. 107-129. ISBN-978-9978-92-737-3.

Penna, N. T. and Stewart, M. P. (2003). Aliased tidal signatures in continuous GPS height time series. Geophysical Research Letters, 30(23). 2184, doi:10.1029/2003GL018828.

941

Pennington, W. D. (1981). Subduction of the eastern Panama Basin and seismotectonics of

942

northwestern South America. Journal of Geophysical Research: Solid Earth, 86(B11), 10753-

943

10770. doi: 10.1029/JB086iB11p10753.

33

ACCEPTED MANUSCRIPT

944

Philippon, M., & Corti, G. (2016). Obliquity along plate boundaries. Tectonophysics, 693, 171-182.

945

Rolandone, F., Nocquet, J. M., Mothes, P. A., Jarrin, P., Vallée, M., Cubas, N., ... & Font, Y. (2018).

946

Areas prone to slow slip events impede earthquake rupture propagation and promote

947

afterslip. Science advances, 4(1), eaao6596. Sage, F., Collot, J. Y., & Ranero, C. R. (2006). Interplate patchiness and subduction-erosion

949

mechanisms: Evidence from depth-migrated seismic images at the central Ecuador

950

convergent margin. Geology, 34(12), 997-1000.

RI PT

948

Sánchez, L., and Drewes, H. (2016). Crustal deformation and surface kinematics after the 2010

952

earthquakes in Latin America. Journal of Geodynamics. DOI: 10.1016/j.jog.2016.06.005.

953

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

M AN U

954

SC

951

955

Schuster, R. L., NietoThomas, A. S., O'Rourke, T. D., Crespo, E., & Plaza-Nieto, G. (1996). Mass

956

wasting triggered by the 5 March 1987 Ecuador earthquakes. Engineering geology, 42(1), 1-

957

23.

Segovia, M., Y. Font, M. M. Regnier, P. Charvis, J.-M. Nocquet, A. Galve, Y. Hello, A. Ogé, P. Jarrin,

959

M. C. Ruiz (2015). Intense microseismicity associated with a SSE at La Plata Island in the

960

central subduction zone of Ecuador. Am. Geophys. Union Fall Meet. 2015, S31A-2736 (2015).

961

Stewart, M. P., N. T. Penna, and D. D. Lichti (2005), Investigating the propagation mechanism of

962

unmodelled systematic errors on coordinate time series estimated using least squares, J.

963

Geod., 79, 479–489. doi:10.1007/s00190-005-0478-6.

EP

TE D

958

Trenkamp R, Kellogg J, Freymuller J, Mora H (2002) Wide plate margin deformation, South Central

965

America and Northwestern South America, CASA GPS observations. Journal of South

969

AC C

964

970

Geophys. Res. Solid Earth 118, 2965–2981.

966

American Earth Sciences 15: 157–171. doi: 10.1016/S0895-9811(02)00018-4.

967

Vallée, M., J.M. Nocquet, J. Battaglia, Y. Font, M. Segovia, M. Régnier, P. Mothes, P. Jarrin, D.

968

Cisneros, S. Vaca, H. Yepes, X. Martin, N. Béthoux, M. Chlieh (2013). Intense interface seismicity triggered by a shallow slow slip event in the Central Ecuador subduction zone. J.

971

Wdowinski, S., Bock, Y., Zhang, J., Fang, P., and Genrich, J. (1997). Southern California permanent

972

GPS geodetic array: Spatial filtering of daily positions for estimating coseismic and

34

ACCEPTED MANUSCRIPT

973

postseismic displacements induced by the 1992 Landers earthquake. Journal of Geophysical

974

Research: Solid Earth, 102(B8), 18057-18070. White, S. M., Trenkamp, R., & Kellogg, J. N. (2003). Recent crustal deformation and the earthquake

976

cycle along the Ecuador–Colombia subduction zone. Earth and Planetary Science Letters,

977

216(3), 231-242. doi: 10.1016/S0012-821X(03)00535-1.

RI PT

975

978

Williams, S. D., Bock, Y., Fang, P., Jamason, P., Nikolaidis, R. M., Prawirodirdjo, L., Miller, M. &

979

Johnson, D. J. (2004). Error analysis of continuous GPS position time series. Journal of

980

Geophysical Research: Solid Earth, 109(B3). doi:10.1029/2003JB002741.

Witt C., Bourgois J., Michaud F., Ordoñez M., Jiménez N., Sasson M. (2006). Development of the

982

Gulf of Guayaquil (Ecuador) during the Quaternary as an effect of the North Andean block

983

tectonic escape. Tectonics. 25 (TC3017): 1-22. doi: 10.1029/2004TC001723.

M AN U

SC

981

984

Witt, C. and Bourgois, J. (2009). Relaciones entre la evolución de la cuenca del Golfo de Guayaquil-

985

Tumbes y el escape del Bloque Nor-Andino. In: Collot J, Sallares V, Pazmiño N (eds) Geología

986

y Geofísica Marina y Terrestre del Ecuador, 1st edn. Argudo & Asociados, Guayaquil-Ecuador,

987

pp 95-106.

Zhang, J., Bock, Y., Johnson, H., Fang, P., Williams, S., Genrich, J., ... & Behr, J. (1997). Southern

989

California Permanent GPS Geodetic Array: Error analysis of daily position estimates and site

990

velocities. Journal

991

doi: 10.1029/97JB01380.

TE D

988

Geophysical

Research:

Solid

Earth, 102(B8),

18035-18055.

AC C

EP

of

35

ACCEPTED MANUSCRIPT

TABLES Table 1. GNSS-velocities of REGME/IGS stations relative to South American plate and related information. East, north and vertical GNSS velocity components (VE, VN, VUP) and uncertainties at

RI PT

95% confidence (svE, svN, svUP) are given in mm/yr. wrms E, wrms N and wrms Up: long-term repeatability for the east, north and vertical components in mm. Obs span: observation period for the time series in decimal years. #obs: number of observations, after eliminating outliers, used to estimate velocity. Alt. (m)

VE

VN

VUp

svE

svN

svUp

wrms E

wrms N

wrms Up

Obs span

#obs

-78.847

-2.202

2383.7

9.90

-3.70

2.00

0.34

0.28

AUCA 42017M001(2)

-76.883

-0.641

313.9

2.70

0.50

-3.42

0.51

0.41

0.79

2.37

1.96

5.56

2.06

596

1.32

2.12

1.71

5.44

1.54

356

BAHI 42018M001(3)

-80.398

-0.659

53.5

10.70

7.50

-3.24

0.70

0.90

BOGT 41901M001(4)

-74.081

4.64

2576.8

5.80

3.30

-33.04

0.06

0.04

BRAZ 41606M001(4)(*) BRFT 41602M002(4)(*)

-47.878 -38.426

-15.948 -3.877

1106.0 21.7

-1.50 0.70

0.10 -0.40

2.30 0.62

0.09 0.09

0.06 0.06

0.23 0.20

BRMU 42501S004(4)

-64.415

32.221

-11.6

-2.80

-3.10

-2.22

0.08

0.08

CHEC 42030M001(2)

-77.814

-0.339

1643.7

9.70

1.10

-1.11

0.61

0.63

CLEC 42031M001(1)

-79.956

-4.103

2024.4

3.80

-6.10

5.07

0.45

0.41

COEC 42023M001(1)

-77.787

0.716

3657.1

10.70

0.90

3.10

0.37

0.26

0.64

2.53

1.80

4.50

2.10

451

CRO1 43201M001(4) CUEC 42009M001(1)

-64.35 -79.003

17.453 -2.883

-31.5 2631.2

16.80 4.30

0.70 -2.10

-0.60 1.65

0.06 0.06

0.05 0.05

0.18 0.13

2.29 1.88

1.81 1.62

6.95 4.12

6.24 6.18

1,849 1,472 264

1.76

1.46

1.89

3.74

1.13

125

0.18

2.23

1.68

6.62

6.24

1,793

3.33 3.30

2.34 2.35

7.54 10.21

6.24 6.24

1,742 1,618

0.15

3.20

3.11

5.48

6.24

1,807

1.59

1.70

1.81

4.51

1.06

345

1.18

2.00

1.84

5.39

1.45

463

M AN U

ALEC 42029M001

(1)

Lat. (°N)

SC

Long. (°E)

Site name

-78.615

-0.935

2813.7

9.20

-0.70

-0.31

0.75

0.66

1.78

1.90

1.69

4.57

1.07

-79.452

-0.272

284.9

17.20

4.50

-8.95

0.29

0.30

0.72

2.02

2.08

5.12

2.10

610

EPEC 42039M001(1)

-78.446

-0.315

2519.1

9.90

-1.80

-3.16

0.71

0.67

1.92

1.90

1.81

5.22

1.11

279

EREC 42037M001(1)

-78.651

-1.671

2801.0

9.80

-5.00

7.47

0.46

0.37

1.03

2.00

1.71

4.69

1.55

426

ESMR 42011M001(1) GLPS 42005M002(4)

-79.724 -90.181

0.935 -0.444

252.9 1.8

22.30 55.70

6.00 2.50

3.03 -0.18

0.06 0.05

0.05 0.04

0.16 0.11

1.82 1.67

1.65 1.64

5.28 4.19

5.65 6.23

1,611 1,383

GOLD 40405S031(4)

-116.532

35.253

986.7

-4.60

-7.90

0.00

0.08

0.08

0.18

3.28

3.03

6.86

6.24

1,835

GUAT 40901S001(4)

-90.311

14.353

1519.9

11.50

-5.80

-1.68

0.07

0.06

0.16

2.46

2.00

5.63

6.22

1,642

GYEC 42007M001(3)

-79.892

-2.149

34.9

10.50

-0.50

-3.98

0.06

0.07

0.15

2.10

2.18

4.78

5.75

1,583

GZEC 42032M001(1)

-78.581

-3.401

884.0

4.20

-3.80

1.59

0.28

0.25

0.70

1.91

1.76

4.89

2.06

596

IBEC 42024M001(1) ISPA 41703M007(4)

-78.116 -109.204

0.35 -27.073

2246.3 112.5

10.00 65.20

2.90 -11.30

-1.69 -0.64

0.24 0.09

0.21 0.07

0.59 0.23

2.16 2.68

1.84 1.72

5.44 8.04

2.60 6.24

689 1,647

EP

TE D

CXEC 42038M001(1) ECEC 42027M001(1)

-79.199

-3.988

2143.5

3.30

-2.20

2.28

0.08

0.05

0.22

3.15

1.71

7.91

5.95

1,617

LPGS 41510M001(4)(*)

-57.932

-34.907

29.9

-1.80

0.40

-3.80

0.12

0.12

0.31

3.04

3.16

7.58

4.73

1,182

LREC 42014M001(2)

-75.986

-1.615

206.5

-6.50

-3.30

3.19

1.44

1.30

3.58

1.98

1.79

5.00

0.65

215

MAEC 42013M001(1)

-78.118

-2.305

1062.5

5.30

-1.70

1.60

0.10

0.09

0.27

2.11

1.82

5.54

4.45

957

MANA 41201S001(4) MDO1 40442M012(4)

-86.146 -104.005

12.086 30.405

71.0 2004.5

14.10 -1.50

-1.90 -15.30

-1.66 3.37

0.05 0.11

0.06 0.10

0.16 0.23

1.89 2.87

2.08 2.72

6.11 6.38

6.24 5.18

1,764 1,537 552

AC C

LJEC 42010M001(1)

MHEC 42022M001(2)

-79.957

-3.261

55.1

8.30

-0.10

0.70

0.29

0.28

0.77

1.97

1.95

5.26

1.97

MTEC 42015M001(2)

-76.982

-2.069

301.7

-3.10

-1.40

7.29

1.03

0.95

2.73

2.04

1.90

5.44

1.20

271

NJEC 42028M001(1)

-79.621

-2.675

51.4

10.10

-3.20

1.36

0.29

0.26

0.85

2.04

1.87

6.11

2.12

607

PARC 41716S001(4)(*)

-70.88

-53.137

22.3

1.30

0.00

-2.35

0.10

0.12

0.26

3.46

4.64

9.02

6.24

1,612

PDEC 42033M001(2) PJEC 42034M001(3)

-79.131 -80.425

-4.648 -1.552

1143.4 146.8

6.00 20.80

-3.70 -0.10

6.26 -4.28

0.34 0.28

0.33 0.34

1.01 0.75

1.99 1.81

1.90 2.18

5.80 4.86

2.06 2.10

466 434

PREC 42035M001(1)

-77.963

-1.708

904.4

3.70

-5.90

13.75

0.48

0.40

1.48

2.10

1.80

6.46

1.56

379

PTEC 42008M001(3)

-80.475

-1.058

61.3

13.00

4.00

3.41

0.11

0.09

0.22

2.41

1.93

4.96

4.53

1,045

QUEM 42020M001(1)

-78.497

-0.237

3054.7

14.40

4.20

-3.50

0.59

0.51

1.41

2.11

1.83

5.08

2.37

231

QUI1 42003S003(3)

-78.494

-0.215

2922.6

17.70

-3.20

-11.01

0.48

0.54

1.83

0.99

1.12

3.72

0.82

288

QVEC 42012M001(2) RIOP 42006M001(1)(4)

-79.47 -78.651

-1.012 -1.651

120.3 2793.0

13.60 6.30

-1.90 -3.20

1.00 0.69

0.16 0.07

0.13 0.05

0.34 0.18

2.58 2.16

2.13 1.63

5.66 5.73

3.95 6.20

866 1,688

SCUB 40701M001(4)

-75.454

20.004

21.9

1.60

-5.80

-1.36

0.08

0.09

0.28

1.81

2.30

6.65

4.70

1,333

SEEC 42036M001(1)

-80.904

-2.22

29.8

14.50

3.40

-1.18

0.38

0.25

0.62

2.71

1.71

4.31

2.06

481

36

ACCEPTED MANUSCRIPT

(2)

SNLR 42021M001 SSIA 41401S001(4)

-78.847 -89.066

1.293 13.415

STEC 42016M001(2)

-78.011

-3.051

321.7

5.00

-1.90

2.51

0.49

0.56

1.53

1.91

2.19

5.94

1.84

277

(2)

-77.816

-0.99

547.1

4.70

-2.20

0.12

0.17

0.17

0.48

2.09

2.03

5.90

3.22

821

UNSA 41514M001(4)(*)

-65.408

-24.728

1257.8

1.50

-0.20

-1.04

0.08

0.07

0.19

2.89

2.58

6.16

6.12

1,817

USNO 40451S003(4)

-77.036

38.551

48.9

-5.40

-8.00

0.49

0.06

0.07

0.23

1.88

2.42

6.89

6.24

1,638

TNEC 42026M001

23.1 664.4

15.50 15.90

-1.50 -4.80

-3.55 2.14

0.84 0.05

0.76 0.06

2.19 0.16

2.10 1.69

1.93 2.34

5.46 5.85

1.07 6.24

249 1,361

(1) REGME station presently in operation. (2) REGME station in maintenance at present. (3)

RI PT

REGME station presently not in operation. (4) IGS station. (*) IGS stations used to calculate the

AC C

EP

TE D

M AN U

SC

Euler pole of the South American plate.

37

ACCEPTED MANUSCRIPT

FIGURES

AC C

EP

TE D

M AN U

SC

RI PT

Figure 1

Figure 1. Tectonic Setting of Ecuador. Black vector indicates the relative rate and azimuth between the Nazca and South America plates (Trenkamp et al., 2002), white vector indicates the movement of the North Andean Block relative to the South American plate (Nocquet et al., 2014). 38

ACCEPTED MANUSCRIPT

Black lines outline main tectonic structures (from Trenkamp et al., 2002). Thin black lines outline the fault map (Egüez et al., 2003). The yellow triangles indicate the position of the volcanoes. Black squares show the location of the REGME stations used in this study. Red circles show instrumental seismicity (Mw> 5, depth <40 km, period 2008-2014) from the USGS (United States

RI PT

Geological Survey) catalog of the National Earthquake Information Center (NEIC). Red stars represent the location of megathrust occurred over the last century. The red dashed line shows the Carnegie Ridge Track (from Chlieh et al., 2014). GG – Guayaquil Gulf, DGM - DoloresGuayaquil Megashear, LPI – La Plata Island, CB – Bahía de Caráquez, SEP – Santa Elena Peninsula,

SC

NAB – North Andean Block, SOAM – South American Plate, CRT – Carnegie Ridge Track. Inset shows the regional tectonic setting, including the plates, micro-plates and blocks that interact in

M AN U

this region: Caribbean, South American, Nazca and Cocos plates, Panama Micro-plate and North Andean Block. The Major Dextral System (MDS) is located in the central part of the country, over

AC C

EP

TE D

the Andes Mountains. Plate boundaries are defined by Bird (2003).

39

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

Figure 2

Figure 2. GNSS horizontal velocity field (black vectors) of REGME stations relative to the South American plate reference frame with 95% confidence error ellipses. Upper inset shows the GNSS horizontal velocity of GLPS station located on the Nazca plate. NAB – North Andean Block, SOAM plate – South American plate. Dashed blue lines mark the position of the profiles of Figures 5 and

40

ACCEPTED MANUSCRIPT

10. The eastern and southern limits of the North Andean Block correspond to a simplification of

AC C

EP

TE D

M AN U

SC

RI PT

the Major Dextral System (MDS).

41

ACCEPTED MANUSCRIPT

Figure 3

b

SC

RI PT

a

M AN U

Figure 3. Photos of two types of monumentation of REGME stations. a) Rooftop monumentation

AC C

EP

TE D

with a 2-meter high concrete pillar, station CUEC. b) Monumentation on the ground, station RIOP.

42

ACCEPTED MANUSCRIPT

M AN U

SC

RI PT

Figure 4

TE D

Figure 4. Station GYEC time series. Black points show daily positions. East, north and vertical components with 95% confidence errors given in global ITRF2008 reference frame. Red line shows

AC C

EP

the trend and green line indicates seasonality.

43

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

Figure 5

Figure 5. Topography and GNSS velocities relative to SOAM projected along three profiles shown in Figure 2. Profile-parallel velocity components (black circles) and profile-normal velocity components (red open triangles) with 95% confidence level error bar are plotted versus distance along the profile. Vertical dashed lines indicate approximate location of the limit east and west of

44

ACCEPTED MANUSCRIPT

the Mayor Dextral System (MDS). Blue line shows the coast line in each profile. NAB – North

AC C

EP

TE D

M AN U

SC

RI PT

Andean Block. SOAM plate – South American plate, INCA S. – Inca Sliver.

45

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

Figure 6

Figure 6. Dilatation rate field (see color bar) and strain rates in the study area calculated based on inter-seismic velocities using Delanunay triangulation (in grey). Vectors show principal axes of the horizontal strain rate tensors. Inward pointing arrows represent compression (in red); outward pointing arrows depict extension (in blue). Negative dilatation rate values indicate compression 46

ACCEPTED MANUSCRIPT

(in red) and positive values show extension (in blue). Black lines outline main tectonic structures. Thin black lines outline active faults (Egüez, et al., 2003). Red stars represent the location of megathrust earthquakes occurred over the last century. Focal mechanisms (Harvard CMT) are

AC C

EP

TE D

M AN U

SC

RI PT

plotted for Mw >7 earthquakes.

47

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

Figure 7

Figure 7. Maximum shear strain rate field (see color bar) in the study area calculated from interseismic velocities using Delaunay triangulation (in grey). Bold red lines indicate orientation of the right-lateral shear strain rate within the triangle. Black lines outline main tectonic structures. Thin black lines outline active faults (Egüez, et al., 2003). 48

ACCEPTED MANUSCRIPT

TE D

M AN U

SC

RI PT

Figure 8

EP

Figure 8. a) Focal mechanisms of the Global CMT catalog (Ekström et al., 2012) classified by type of rupture: red for reverse, yellow for strike-slip and blue for normal earthquakes. b) Calculated

AC C

horizontal slip vector for each focal mechanism. The size and color of the arrow is proportional to the slip amount (see text for details). Date labels indicate major events.

49

ACCEPTED MANUSCRIPT

EP

TE D

M AN U

SC

RI PT

Figure 9

AC C

Figure 9. Preferred interseismic coupling model along the Ecuadorian subduction zone obtained from the inversion of our GNSS velocities. Blue indicate uncoupled areas, red represent highly coupled areas (asperities) and green and yellow represent partially coupled regions. Black and red vectors represent observed and modelled horizontal velocities respectively, with 95% confidence error ellipses. DGM - Dolores-Guayaquil Megashear, CRT – Carnegie Rigde Track, CB – Bahía de Caráquez, LPI – La Plata Island, SEP – Santa Elena Peninsula, GG – Guayaquil Gulf, NAB – North Andean Block.

50

ACCEPTED MANUSCRIPT

M AN U

SC

RI PT

Figure 10

TE D

Figure 10. Interseismic deformation accumulated in the study area obtained removing to the horizontal velocity field relative to SOAM, the rotational part of each block used in the modelling (Figure 9). a) GNSS horizontal velocities observed vectors minus rotational part with 95% confidence error ellipses. NAB – North Andean Block, SOAM plate – South American plate. Dashed

EP

blue lines mark the position of the profiles in Figures10b-d. The eastern and southern limits of the North Andean Block correspond to a simplification of the Major Dextral System (MDS). b-d)

AC C

Topography and Figure 10a velocities projected along three profiles shown in Figure 10a. Profileparallel velocity components (black circles) and profile-normal velocity components (red open triangles) are plotted versus distance along the profile. Vertical dashed lines indicate approximate location of the limit east and west of the Mayor Dextral System (MDS). Blue line shows the coast line in each profile. NAB – North Andean Block. SOAM plate – South American plate, INCA S. – Inca Sliver.

51

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

Figure 11

Figure 11. Obliquity and partitioning analysis of the Nazca – North Andean Block subduction. a) Azimuth of different slip vectors on the Nazca – NAB subduction interface. The blue shaded area corresponds to the dip direction of the subduction interface for depths up to10 km (equivalent to the usual approximation of trench normal orientation). The dark blue line represents the azimuth trend obtained with a gaussian convolution filter. The orange circles represent the azimuth plus 180° of the slip vector on the eastward dipping plane of the reverse focal mechanisms of the

52

ACCEPTED MANUSCRIPT

subduction interface. The size of the circle is proportional to the earthquake magnitude. The orange line represents the focal mechanisms slip vector trend. The red line shows the azimuth of the velocity vector on the subduction interface computed for the Nazca – South America. b) Obliquity defined as the angle between the trench normal and the motion vector. γ is the plate

RI PT

convergence obliquity computed for the NZ-SA pole. Ψ is the focal mechanisms slip vector obliquity (see inset for details), δ is the difference between both angles. c) Values of the seismic partitioning (ps) and the subduction interface coupling (φ). The dots show the coupling computed at different interface depths. The purple line shows the maximum computed interface coupling. d)

SC

Velocities. vp and vn are the margin parallel and margin normal components of the plate vector motion. vs is the forearc sliver slip rate relative to the upper plate. The dotted line show the rates

M AN U

obtained assuming a fully coupled subduction interface. The solid lines show the actual rates

AC C

EP

TE D

taking into account the subduction interface coupling shown in (c).

53

ACCEPTED MANUSCRIPT The highlights of our work “CRUSTAL MOTION AND DEFORMATION IN ECUADOR FROM cGNSS TIME SERIES” are: – First velocity field from time series analysis of cGNSS stations belonging to REGME network (Continuous Monitoring GNSS Network) in Ecuador.

RI PT

– In northern Ecuador, there is an estimated right-lateral motion of 7.6 ± 0.5 mm/yr, consistent with the NNE movement of the NAB relative to the South American plate. In central Ecuador, the right-lateral motion decreases to 5.3 ± 0.4 mm/yr. In southern Ecuador (from the Guayaquil Gulf to Peru) there is no strain accumulation, the GNSS velocities decrease and turn to the south. This zone belongs to the so-called Inca or Peru sliver.

AC C

EP

TE D

M AN U

SC

– The main driving force responsible for ongoing crustal deformation in Ecuador is the convergence between the Nazca and South American plates with the variable coupling pattern and the collision of the Carnegie Ridge. This produces the significant change in the velocity pattern obtained from north to south.