Bridging the gap between North and Central Chile: insight from new GPS data on coupling complexities and the Andean sliver motion

GPS surveys have been extensively used over the past 20 yr to quantify crustal deformation associated with the Andean subduction zone in Chile. Such measurements revealed the coupling variations associated with the seismic segmentation of the subduction. However, because of data gaps mostly due to access difﬁculties, the Atacama–Antofagasta regions of North Chile remain poorly known. We present here an upgraded interseismic velocity ﬁeld aggregating new data acquired between 2012 and 2016 in the region of Taltal (24 ◦ S–26 ◦ S), over a small-scale network of 20 benchmarks. This denser data set reveals a new complexity regarding the modelling methodology commonly used. We ﬁrst show that a large-scale rigid Andean sliver, running from central to North Chile, does not allow to explain the velocities measured near the cordillera in the region of Taltal. This region exhibits an additional coherent block motion of almost 8 mm yr − 1 with respect to the inland motion generated by the rotation of the sliver proposed by, for example, Brooks et al. 2003; M´etois et al. 2013, 2014, which works well everywhere else. Second, once this local block motion is taken into account, the coupling in the Taltal area is reﬁned, which brings new insights about the subduction segmentation there. The Taltal area shows as a relative low in coupling (although coupling values are still high), potentially cutting a long section of the subduction into two independent highly coupled segments: the Paranal segment—north of Taltal, between 23 ◦ S and 25 ◦ S—and the Cha˜naral segment—south of Taltal, between 26 ◦ S and 28 ◦ S. These segments may rupture individually with magnitude ∼ 8 earthquakes or simultaneously which would produce a larger earthquake, especially if a third segment (Atacama—more to the south—between 28 ◦ S and 29 ◦ S) is also involved.


I N T RO D U C T I O N
Twenty years of GPS campaign carried out along the Chilean subduction interface allowed to estimate a very dense and quasicontinuous interseismic velocity field from the Peninsula of Arauco (38 • S) up to Arica (18 • S) (Klotz et al. 2001;Khazaradze et al. 2002;Brooks et al. 2003;Ruegg et al. 2009;Vigny et al. 2009;Chlieh et al. 2011;Métois et al. 2012Métois et al. , 2014Métois et al. , 2016. Interseismic velocities are commonly modelled as a combination of elastic deformation of the upper plate due to coupling on the subduction interface and the * Now at: Institute of Geophysics and Planetary Physics, Scripps Institution of Oceanography, University of California, San Diego, La Jolla, CA 92093-0225, USA. rigid rotation of an Andean sliver block (Brooks et al. 2003;Chlieh et al. 2011;Métois et al. 2016). The coupling variations reveal a clear segmentation of the interface: long highly coupled segments separated by shorter low coupled zones (LCZs) (Moreno et al. 2011;Métois et al. 2012;Béjar-Pizarro et al. 2013;Métois et al. 2016). LCZs seem to behave as barriers to the propagation of megathrust ruptures and highly coupled segments correlate very well with the rupture zones of the recent megathrust earthquakes: Maule (M w 8.8 2010, Moreno et al. 2010;Vigny et al. 2011;Métois et al. 2012), Iquique (M w 8.2 2014, Schurr et al. 2014;Ruiz et al. 2014) and Illapel (M w 8.3 2015, Ruiz et al. 2016;Métois et al. 2016;Klein et al. 2017). Moving forward, a challenging possibility opened by recent geodetic measurements in poorly known areas is to identify the potential location and extent of future large earthquakes. The last three megathrust earthquakes (2010,2014,2015)   Background seismicity  from USGS, represented in colour scale as function of depth. Topography and bathymetry from ETOPO5, white lines represent the main bathymetric features of the Nazca plate (Álvarez et al. 2015). Black lines depict the Atacama Fault Zone (AFZ) and the Taltal Fault (TF) (Lemke et al. 1968;Armijo & Thiele 1990). (b) Segmentation (Métois et al. 2016) and historical seismicity as function of time, M w ∼ 7 earthquake epicentres are depicted by grey stars, seismic swarms by blue pentagrams. Black lines represent rupture zones of historical (dashed) and instrumental (continuous) of major events, their epicentres are depicted by black crosses. apparent seismic gaps: (i) the Metropolitan region offshore Valparaiso (from 32 • S to 34 • S, last ruptured in 1906 and 1985); (ii) North Chile (between 20.5 • S and 23 • S unruptured since 1877); (iii) the Atacama region (between 30 • S and 25 • S). There, a seismic cycle of about a century seems to emerge with the occurrence of two M+ 8.5 earthquake in 1819 and 1922 (Lomnitz 2004;Engdahl & Villaseñor 2002;Fig. 1). A range of magnitudes have been proposed for the last 1922 earthquake: M s 8.3 (Abe 1981), M w 8.4 (Lomnitz 2004), recently re-evaluated to M w 8.5-8.6 by (Carvajal et al. 2017).
This variety may also relate to a real complexity of the seismic sources. The 1819 earthquake is in fact a triple event occurring several days apart (3, 4 and 11 April) and the 1922 earthquake is also described as two great shocks following each other within 8 min (Lomnitz 1970). Possibly in adequacy with this complexity, GPS measurements in this region revealed high coupling values between 30 • S and 23 • S, but with a complex coupling pattern and at least two identified segments between 30 • S and 25 • S: Atacama and Chañaral (see Fig. 1b; Métois et al. 2016). Going north, another highly coupled segment has been imaged by the GPS data south of 23 • S, the Paranal segment, where the 1995 M w 8 Antofagasta earthquake could have broken only the deepest portion of the interface (Ruegg et al. 1996;Chlieh et al. 2004;Métois et al. 2013). Because of the late human settlement in this desertic region, whether and when the last megathrust rupture occurred there remains unknown.
The boundary between Atacama and Chañaral segments is marked by a low coupling portion, the Baranquilla LCZ (see Fig. 1b; Métois et al. 2014Métois et al. , 2016. By contrast, a boundary between Chañaral and Paranal segments may exist somewhere in the region of the city of Taltal (25.4 • S) but is not clearly identified since this area remained one of the last data gap with only four geodetic markers installed and surveyed over a 200 km long area until now (Klotz et al. 2001). Therefore, the northern limit of the Chañaral segment and the southern limit of the Paranal segment remain unclear and the behaviour of the portion of interface in front of Taltal unknown, opening ways for two distinct scenarios. First, an LCZ separates the Chañaral and Paranal segments, behaving as a barrier to the propagation of past earthquake ruptures (1819 and 1922 along the Atacama segment and 1995 along the Paranal segment). This LCZ could act similarly for future large earthquakes, occurring on either sides. This hypothesis could be supported by the presence of a known geological structure, the Sala y Gomez Ridge, which is subducted beneath the South American plate offshore the Taltal bay (see Fig. 1). Such a correlation between bathymetric structures and LCZs was already observed in Chile but remains debated (i.e. the Copiapó Ridge faces the LCZ of Baranquilla and the Challenger Fracture Zone faces the LCZ of La Serena, Métois et al. 2016;Collot et al. 2017). The second scenario implies a single and continuous 500-700 km long segment, highly coupled, capable of producing a giant earthquake of M w +9, and only partially broken by the last known earthquakes. Such an earthquake, breaking the two segments of Chañaral and Paranal (and possibly a third one, Atacama, if going across the Baranquilla LCZ also), could trigger a major tsunami with devastating impact along a very large section of the Chilean coastline.
Here we present new interseismic velocities acquired in the area since 2012 and that help arbitrate between both hypothesis.

S E I S M O T E C T O N I C C O N T E X T
In this study, we focus on the Taltal area (24 • S-26 • S). At these latitudes, the Nazca plate subducts beneath the South American plate with a convergence rate of 67 mm yr −1 (Angermann et al. 1999;Vigny et al. 2009). This area is characterized by the presence of one of the main crustal fault system in Chile, the Atacama Fault Zone (AFZ; Fig. 1a). Extending over more than 1000 km, it separates the uplifting Coastal Cordillera from the compressional fore-arc, but is characterized by very little seismic activity to date and most probably presently inactive (e.g. Arabasz 1968;Okada 1971;Armijo & Thiele 1990). It is divided in two segments cut at the latitude of the Taltal bay by the eponymous NW-SE fault, corresponding also to the latitude where the Sala y Gomez Ridge enters the subduction (Fig. 1). In the region of Taltal, maybe due to the limited completeness and historical continuity of the catalogue, only one major event is known, a M w 7.7 earthquake in 1966 but several events of M w close to 7 have occurred in the last century (Deschamps et al. 1980, Fig. 1). In the whole area considered here, the shallowest part of the interface shows little background seismicity, at the exception of the repeated occurrence of seismic swarms offshore Caldera (in 1973Caldera (in , 1979Caldera (in and 2006Holtkamp et al. 2011, Fig. 1b). On the contrary, the deep seismic activity appears very heterogeneous: we observe a clear gap of deep seismicity (>100 km depth) between 27 • S and 25 • S, contrasting with the northernmost part where intense deep seismicity is observed.

GPS data acquisition and processing
We installed and surveyed the first markers in the Atacama region immediately after the Maule earthquake of 2010 (Métois et al. 2014). Starting in 2012, 20 additional benchmarks were installed in the region of Taltal and measured annually. We also remeasured three CAP (Central Andes GPS Project, see Brooks et al. 2003) and eight SAGA (South America Geodynamic Activities, see Klotz et al. 2001) markers in the area, building a data set of almost 70 points. Overall, all included sites were measured at least three times and up to seven times (see Supporting Information Table S3).
We also include observations from the recently installed regional continuous stations operated by the Chilean National Seismolog-ical Center (CSN) and from a selection of well-distributed permanent stations across the South American continent [IGS (Beutler et al. 1999), Red Argentina de Monitoreo Satelital Continuo (RAM-SAC) and Rede Brasileira de Monitoramento Contínuo (RBMC) networks]. This data set is processed using the GAMIT/GLOBK software following the classical MIT methodology (Herring et al. 2010a,b). Daily solutions are combined using a global stabilization approach because of the post-seismic deformation following the Maule earthquake which affects many reference stations on the South American continent (see the Supporting Information and Klein et al. 2016). Plotted relative to the South American Plate defined by the NNR-Nuvel-1A model (25.4 • S, 124.6 • W, 0.11 • Myr −1 , DeMets et al. 1994), estimated velocities of stations located on the stable part of the South American continent present a systematic northward residual of about 2 mm yr −1 (see Supporting Information Fig. S1a). Therefore, we map our solution in a new Stable South America reference frame. The Euler pole location and angular velocity of this plate (20.9 • S, 128.6 • W, 0.122 • Myr −1 ) were inverted using a subset of seven sites that exhibit rigid block motion with respect to the ITRF2008 (Altamimi et al. 2012, Supporting Information Fig. S1b).
Interseismic velocities from Métois et al. (2013) north of the Peninsula of Mejillones are combined with this data set by applying the rotation vector estimated in this study. Note that data from previous and older studies (Klotz et al. 2001;Brooks et al. 2003, represented in Fig. 1) were not included in the solution for consistency purpose. We also excluded 2016 position measurements at several sites between Copiapó (27.3 • S) and Chañaral (26.4 • S) because they seemed to have been affected by a transient deformation of unknown origin occurring in this area between the two measurement campaigns of 2015 and 2016. Finally, we excluded 2016 measurements south of Copiapó to avoid any contamination due to the post-seismic deformation following the 2015 Illapel earthquake and the time-series of sites in the region of interest show that they are not contaminated further north (see Supporting Information Fig. S3). Because the 2014 Iquique earthquake that occurred at more than 400 km from the north of our network was only M w ∼ 8.1, we are confident that no post-seismic motion associated with this distant quake could contaminate our data set.
As mentioned previously, post-seismic deformation at continental scale has been highlighted following the 2010 Maule earthquake, not only directly in front of the rupture zone up to the east coast of Argentina, but also north of the rupture zone (Klein et al. 2016). This effect was indeed measurable in the region of the 2015 Illapel earthquake, where an increase of north-eastward velocity of some 2 mm yr −1 was reported between 2010 and 2015 (Ruiz et al. 2016). Yet, north of La Serena (30 • S), the residual velocities between pre-2010 and post-2010 velocities no longer show any consistent pattern expected from such processes, especially from campaign measurements (cf. fig. 4 in Ruiz et al. 2016). The studied region is located even further north, that is, some 1000 km from the Maule rupture zone. Time-series of both continuous stations and campaign sites are provided to confirm their linearity and therefore the quality of velocities estimation, also regarding seasonal signals (see Supporting Information Figs S2 and S3). A prediction of the viscoelastic model of post-seismic deformation following the Maule earthquake (Klein et al. 2016) also shows that this effect is not any more affecting significantly the surface deformation at these latitude (cf. Supporting Information Fig. S5).
Therefore, our data set encompasses only the linear trend determined at all sites over a relatively short time window (6 yr) but constrained with yearly measurements (up to 7 positions), hopefully representative of long term inter-seismic deformation. The vertical component of GPS data being less precise and more strongly affected by any hydrological signal than the horizontal component, a longer time span or the knowledge of the long-term hydrological signal derived from neighbouring cGPS stations is necessary to constrain a robust vertical velocity field. Therefore, since very few and recent cGPS stations are in place in the study area, we do not use vertical velocities in this study.

Velocity field description
In general, velocities show a typical interseismic loading pattern, with a north-eastern orientation and coastal vectors roughly parallel to the plate convergence direction (Fig. 1a). We observe the usual clockwise rotation and decrease of amplitude going inland. However, the velocity gradient varies along strike. A very strong gradient of almost 10 mm yr −1 over the first 30 km off the coast is detected between Caldera (27 • S) and Chañaral (26.3 • S) (cf. Supporting Information Fig. S6-profile E) when a much weaker gradient is observed north of (26.3 • S). Moreover, coastal velocities increase significantly from 25 mm yr −1 in the south of the network  Fig. S6-profile C), compared to less than 20 mm yr −1 at the same distance from the trench but few hundreds of kilometres to the south (Supporting Information Fig. S6-profile F). This atypical pattern (weak gradient associated with high velocities) does not correspond to the standard elastic accumulation model and has often been interpreted as due to the motion of a uniformly rigid block in the Andes. However, other unidentified deformational processes could also be invoked such as distributed deformation on crustal faults over the mountain belt.

I N T E R S E I S M I C C O U P L I N G A L O N G T H E TA LTA L A R E A ( 2 • S -2 6 • S )
Using the TDEFNODE code developed by McCaffrey (2015) (version 2015.04.20), the SLAB1.0 geometry (Hayes et al. 2012) and the modelling settings described in the supplementary material and in Métois et al. (2016), we first calculate the residuals produced by a simple bimodal coupling model (zero or 100 per cent coupling) derived from the previously published large-scale models by (Métois et al. 2016;Fig. 2). In this model, the upper-crust deformation is mimicked by the combination of elastic deformation due to partial locking of the Nazca-South America megathrust together with the rigid motion of an Andean sliver relative to the stable interior of South America. The sliver motion is described by the Euler pole determined by Métois et al. (2016) and located in the South Atlantic ocean (56.37 • S, 41.27 • W) with an angular velocity of −0.12 • Myr −1 . It generates an average 8 mm yr −1 northeastward motion in our study area and nearly 10 mm yr −1 northeastward motion over the Bolivian fold-and-thrust belt front, consistently with the few GPS velocities available in the far-field (Brooks et al. 2011;Weiss et al. 2016).
This Andean sliver is a widely used modelling trick to remove from the observed velocity field the a-typical translation-like motion observed in the middle to far field, that is, in the main cordillera, and that is roughly parallel to the convergence velocity between the Nazca and South American plates (Chlieh et al. 2011;Nocquet et al. 2014;Métois et al. 2016). Since the eastern boundaries of such a sliver remain vague and that it is rather likely that the Chilean upper plate deforms in a more complex way than a purely rigid body, we call for caution in interpreting such sliver motion. This coarse firstorder model presented in Fig 2 highlights the defining feature of the Taltal area, while the residuals are small north of the Mejillones Peninsula (23.1 • S) and south of Caldera (27 • S), large (>5 mm yr −1 ) and systematic eastward pointing residuals are observed over the Taltal area (24.5 • S-25.5 • S, see Fig. 2). The transition between the residuals on the edges and those in the centre of the zone draw two counter-clockwise rotation cells. Such residuals demonstrate that the 8 mm yr −1 translation-like motion attributed to the large scale Andean sliver motion is not sufficient to explain the large eastward velocities observed inland in the Taltal area. Whatever the coupling distribution used on the seismogenic zone (0-80 km depth), no acceptable fit could be obtained for the inland velocities at the latitude of Taltal if keeping a 8 mm yr −1 eastward moving sliver (see Supporting Information Fig. S8). Instead, a locally faster rotation of the Andean sliver relative to stable South America seems to be needed, thus challenging the sliver rigidity assumption used in previous studies (Chlieh et al. 2011;Nocquet et al. 2014;Métois et al. 2016).
Such a pattern of residuals could also suggest a viscous effect related to post-seismic deformation due to distant earthquakes. The viscoelastic relaxation following the Maule earthquake indeed predicts a whirlpool effect on both sides of the rupture zone (Klein et al. 2016): a clockwise rotation in the north of the rupture zone which ends with eastward velocities roughly aligned with interseismic velocities in the Illapel region (Ruiz et al. 2016), and a counter-clockwise rotation in the south of the Maule rupture zone. In the present case, two events would be necessary to explain the observed residuals: one in the north responsible for the counterclockwise part of the rotation, and a second in the south generating the clockwise rotation (Fig. 2b). The amplitude of the residuals would imply recent and large events (no more than 50 yr ago and M w ∼ 8.5) for which we do not have any historical record. The 1922 Atacama earthquake occurred about a century ago, which is too long ago to produce such an effect. The 1995 Antogafasta earthquake that occurred 22 yr ago could only explain the residuals in the northern part of our study area (between Antofagasta and Taltal) but not south of Taltal. Therefore, the unusually large middle to far field eastward velocities observed in the Taltal area could not be simply explained by post-seismic rebound induced by any of the known seismic events in the area (1922 Atacama or 1995 Antofagasta).
To date, too few GPS velocities spanning the inner cordillera deformation are available to build a consistent alternative model of the Andes deformation (Brooks et al. 2003;McFarland et al. 2017). Therefore, the rigid-sliver modelling trick remains the simplest assumption in the frame of a purely elastic modelling strategy, to account for the non-elastic deformation of the Chilean upper plate and correctly invert for the coupling distribution on the megathrust (see Supporting Information Fig. S10). Since this cannot be done here with a single large scale rigid sliver over the entire margin (38 to 18 • S), we chose to model the velocity field observed in the Taltal area by inverting simultaneously for the coupling distribution on the megathrust and the Euler pole of a local sliver. Fig. 3 presents three models obtained with varying constrains on the coupling distribution: (i) imposed down-dip decrease, (ii) Gaussian distribution of coupling with depth inverting for amplitude, central depth and width of the Gaussian, and (iii) independent functions of coupling with depth. All three distributions lead to reasonable fit to the data with nRMS < 1.7 and non-systematic residuals.
All three models (Fig. 3) highlight significant along-strike and along-dip changes in the amount and distribution of coupling on the subduction plane north and south of the Taltal bay (25.5 • S). First, when coupling is forced to decrease with depth, the best model exhibits a highly coupled zone offshore the Paranal segment, while coupling is lower and deepens in the Chañaral segment. When coupling varies as a Gaussian function of depth, the fit to the data is slightly improved by reducing the sliver motion to 16.6 mm yr −1 and by allowing full coupling on the deep part of the seismogenic zone, from 20 to 60 km, in the Chañaral segment. When only a 2-D smoothing is applied on nodes following the shape spread smoothing technique proposed by McCaffrey (2015), the nRMS is reduced (1.57 versus 1.69 and 1.65 for previous models), the best sliver motion is ∼16 mm yr −1 , and the same shallowing of the highly coupled zone is observed in the Paranal segment compared to the Chañaral segment. Increasing step by step the complexity of the coupling distribution brings out the common and robust features that are really required to fit the data and helps defining a reliable segmentation of the margin: two highly coupled segments named after Métois et al. (2016), the Paranal segment from 23 • S to 25 • S where the highly coupled zone is shallow, and the Chañaral segment from 25 • S to 27 • S where the highly coupled zone is deeper, are separated by an anomalous low in coupling in the Taltal bay.

D I S C U S S I O N
Several aspects of the best model (thus corresponding to the third model presented on the previous section) showed in Fig. 5(c) have to be considered carefully: first, coupling is not resolved over the first 10 km depth of the slab (see sensitivity map in Fig. 4 and checker board tests in Supporting Information Fig. S7) preventing us from interpreting the value of the very shallow coupling in terms of tsunami hazard for instance. Second, the hypothesis of a local rigid sliver distinct from a large scale Andean micro-plate is questionable (see Section 5.1). However, all models agree on two major points: (i) the observed velocities in the Taltal area require the local Andean sliver to rotate around a Euler pole still located in the South Atlantic ocean but producing a 16-17 mm yr −1 northeastward translation towards South America, that is, more than 8 mm yr −1 higher than the rigid motion expected by Métois et al. (2016); (ii) the Taltal bay (25.5 • S) is a transition zone between two highly coupled segments: the Paranal segment to the north and the Chañaral segment to the south.

The Andean sliver: one rigid block?
There is a trade-off between the amount of coupling, its distribution with depth on one side and the amplitude of the sliver rotation on the other side: inverting for coupling variations only while using the predicted 8 mm yr −1 Andean sliver motion in the area would fail to retrieve the observed velocities and would require a deep and probably unrealistic highly coupled zone (down to ∼70 km depth, nRMS 2.39; see Supporting Information Fig. S8). Opposite, a larger than 16 mm yr −1 Andean sliver motion would decrease the average coupling and reduce the highly coupled zones to the shallowest portion of the slab. This trade-off has been explored in details by several authors (for instance in Chlieh et al. 2011;Métois et al. 2013Métois et al. , 2016. Therefore, our main finding showing a shallowing of the highly coupled zone in the Paranal segment compared to the deepening of the southernmost Chañaral segment when using a local homogeneous sliver rotation could also result from a progressive increase in the back-arc rotation amount from north to south, between 25 • S and 27 • S. As stated before, too few GPS velocities are still available in the medium to far field away from the subduction zone to untangle this issue. It is to note however that the published velocity fields spanning the eastern front of the Andes so far (Brooks et al. 2011;McFarland et al. 2017) are consistent with the average motion predicted locally by our rigid-Andes approach.
To the best of our knowledge, our study raises a new concern: a simple large scale rigid Andean sliver scheme does not provide a satisfactory physical model for describing the large scale straining of the Andes. Indeed, no common Euler motion could retrieve  the surface velocities observed on the western cordillera front in North Chile and metropolitan area together with the ones observed in the Taltal region. The motions for each of these zones inverted by (Métois et al. 2013(Métois et al. , 2014 and this study imply a 10 mm yr −1 north-eastward motion of the North Chile area (from 18 • S to 24 • S, confirmed by the recent velocities acquired by Weiss et al. (2016) in Bolivia), ∼16 mm yr −1 in the Taltal area (24 • S to 27 • S, confirmed by the recent paper by McFarland et al. (2017) and less than 5 mm yr −1 south of La Serena (30 • S). No large or small scale shear zone, parallel to the convergence velocity and that could accommodate such gradients on the long-term has been detected to date. A possible explanation for this discrepancy is that the so-called 'Andean sliver' which north and south boundaries remain furthermore unclear, is rather straining in a more distributed way that remains to be understood. This study therefore challenges the commonly used rigid Andean sliver hypothesis.
Another possibility would be that part of the motion interpreted so far as sliver-related would results from the improper elastic assumption used in our modelling. Previous studies showed that viscoelastic assumption of the asthenosphere properties would predict significantly larger horizontal velocities in the medium field at the end of the cycle (e.g. in the Indonesian subduction zone or in North Chile, Trubienko et al. 2013;Li et al. 2015). Therefore, at least part of the deformation associated with rigid blocks in previous elastic models (eg Simons et al. 2007;Métois et al. 2016) could actually be related to the seismic cycle on the subduction interface, assuming viscoelastic properties of the asthenosphere. This hypothesis, going beyond the scope of the present work, will be further investigated in future studies. It is to note that assuming a purely elastic behaviour rather that a viscoelastic behaviour impacts only very slightly the along-strike coupling variations that depends on the strain rate in the very near field of the fault (<150 km from the trench, see for instance Li et al. 2015). We are therefore confident that our analysis in terms of kinematic segmentation of the margin remains correct.

The Taltal transition zone
The Taltal bay marks in all models a depth transition between the highly coupled segments of Chañaral and Paranal (Fig. 5b). Indeed, the average depth of the highly coupled zone is close to 20 km depth for the Paranal segment and 30 km depth (close to being under the coastline) for the Chañaral segment, as highlighted in Fig. 5(a). The shallower locking depth in the Paranal segment compared to the Chañaral segment could result from the 1995 Antofagasta event that may have increased the stress on the unruptured shallowest part of the interface in the case the rupture did not reach the trench, which remains debated (Ruegg et al. 1996;Delouis et al. 1997). The Taltal area presents also a local minimum in the average coupling calculated on the convergence direction on 20 km along-strike wide sliding windows, whatever the depth interval chosen to integrate the coupling coefficient (see Fig. 5b). However, it is to note that the average coupling in the Taltal area remains high (∼ 50 per cent) compared to many low coupling zones (LCZ) previously identified (e.g. the Mejillones LCZ (23.1 • S) or the Baranquilla LCZ (27.5 • S); Métois et al. 2016, see Fig. 5c). As stated in Section 1, the Taltal bay is the place where (i) the Sala y Gomez volcanic ridge, roughly parallel to the convergence direction, subducts (Fisher & Norris 1960;Maksymowicz 2015), (ii) the Antofagasta segment of the Atacama Fault Zone ends (Arabasz 1968;Armijo & Thiele 1990), (iii) several seismic events of M w <7 (Deschamps et al. 1980) occurred in the past, and (iv) the 1995 Antofagasta and the 1922 Copiapó earthquakes stopped (see Figs 1b and 5d). The whole Taltal region correlates also with a deficit in very deep seismicity compared to the La Serena and North Chile area (Fig. 1a).
Together, these observations confirm that the Taltal bay (25.5 • S) plays an important role in the segmentation of the megathrust interface and opens the question of the physical control of the coupling coefficient. Indeed, the Taltal low coupling anomaly could either be controlled by the subducting sea mounts of the Sala y Gomez ridge, which is supported by the apparent parallelism between the low coupling area and the ridge orientation, or controlled by the upperplate properties and complexities like the Atacama Fault Zone for instance. A detailed analysis of the background seismicity and indepth imaging of the margin would be required to untangle these controlling factors (interface roughness, fluid amount, connection between crustal faults and the subduction plane etc.).

Plausible scenarios for future ruptures
Better understanding the coupling distribution from the Mejillones LCZ (23.1 • S) to the Baranquilla LCZ (27.5 • S) has a direct impact on rupture scenarios that could be assessed for the region. Previous studies had detected the onset of the Paranal and Chañaral segments by imaging high coupling in front of Antofagasta (Métois et al. 2013) and north of Copiapó (Métois et al. 2014). However, whether these two segments were continuous and therefore able to rupture with a very large earthquake or separated by an LCZ capable of blocking rupture propagation was unknown. Our results show that the Taltal area has a unique behaviour in terms of kinematics, seismicity and morphology and that a narrow low in coupling is observed in the Taltal bay. This Taltal anomaly could result from a more velocity-strengthening behaviour of the interface under the bay, therefore able to creep during the interseismic period and to slow-down or stop a rupture coming from the neighbouring segments depending on its size. However, the average coupling under the Taltal bay is ∼50 per cent, that is, higher than other LCZ detected along the margin, suggesting that this zone may not be fully free to creep during the interseismic period. Rather, we suggest that this intermediate coupling area is apparently partially locked due to the stress shadows associated with the neighbouring velocityweakening and highly locked asperities (see Métois et al. 2012;Hetland & Simons 2010, for discussion).
These highly coupled segments north and south of the Taltal bay could therefore rupture alone with an M w 8+ earthquake that would be stopped by the Taltal velocity-strengthening area, alike the 1819, 1922, 1987 and 1995 earthquakes; or jointly, if the rupture would propagate through the Taltal LCZ.
Such a velocity-strengthening narrow zone surrounded by velocity-weakening asperities is also a good candidate for nucleation and initiation of a rupture since stress is prone to accumulate more in these zones of mechanical transition. At least two of the three recent megathrust earthquakes that stroke Chile since 2010 nucleated in a zone of intermediate coupling, that is, at the transition between velocity-strengthening and velocity-weakening behaviour: the 2010 Maule earthquake nucleated in an area that experienced very high rates of aftershocks and large post-seismic slip (Vigny et al. 2011;Métois et al. 2012) and the 2014 Iquique earthquake nucleation sequence developed in the Iquique LCZ and around the Camarones segment (Ruiz et al. 2014;Schurr et al. 2014).
These observations suggest that the Taltal bay could very probably host the nucleation point of a future megathrust earthquake that could therefore propagate north and south, rupturing the Paranal and Chañaral segments jointly. Such megaearthquake rupturing the interface from Mejillones to Baranquilla (i.e. 525 km) would reach at least a magnitude close to M w 8.5, and possibly even larger if it also cuts across the Baranquilla LCZ and involves the Atacama segment.
The two previous megathrust earthquakes of 1819 and 1922 likely ruptured the two southern segments of Atacama and Chañaral. Their ruptures were complex, involving two or three main shocks (Willis 1929), which is compatible with simultaneous or rapidly subsequent ruptures of several segments, first Atacama, then Chañaral. Since no damages were reported north of the city of Taltal, we conclude that the rupture was probably stopped by the Taltal LCZ and did not propagate north of it, into the Paranal segment. A repeat of these earthquakes is plausible to occur soon, since both segments are fully coupled and have accumulated enough deformation since 1922 to produce M w 8+ earthquakes. However, a propagation across the Taltal LCZ into the northern most segment of Paranal is a distinct possibility given its relative weakness. Another scenario of a northern rupture involving only Paranal and Chañaral segments but not Atacama is also plausible even though not supported by any known historical earthquake, because of the presence of the stronger Baranquilla LCZ to the south. Any of these two scenarios involving rupture of the Paranal segment (on top of other segments) would not only increase the magnitude of the expected earthquake but more importantly the magnitude of the subsequent tsunami, since the coupling in this northernmost segment seems to extend at shallow depths. Recently, (Carvajal et al. 2017) have revised the magnitude of the 1922 earthquake to M w 8.6 based on far-field tsunami data. A large tsunami produced by the rupture of the Atacama and Chañaral segments could result from coupling extending to shallower depth before 1922 than presently in the ruptured segments (Atacama and Chañaral) implying either (i) temporal variations of the coupling coefficient during the seismic cycle or (ii) underestimated coupling in these shallow area due to low sensitivity of our data, or could result from dynamic propagation of the rupture up to the trench.
In all cases, because of the poor knowledge of historical megathrusts in the region, the Taltal area offers a unique opportunity to test the capability of coupling models alone to sustain correct rupture scenarios.

A C K N O W L E D G E M E N T S
This work was performed in the frame of the French-Chilean LiA 'Montessus de Ballore' with financial support of the French Agence Nationale de la Recherche (ANR-2012-BS06-004) and the European Marie Curie Initial Training Network Zooming in between plates (ITN FP7 Grant Agreement no.: 604713). We are also thankful to the French Institut National des Sciences de l' Univers (INSU-CNRS) and the Réseau Sismologique & Géodésique Français (RE-SIF) for providing the geodetic instruments for all campaigns. We would like to warmly thank E. Saldaño, M.C. Valderas-Bermejo and all the CSN personnel, for their precious help for the field work. We thank M. Moreno for access to previous data sets on GFZ benchmarks that were remeasured for this study. We finally thank our two anonymous reviewers for their constructive comments, which helped to improve this manuscript.

S U P P O RT I N G I N F O R M AT I O N
Supplementary data are available at GJI online. Table S1: Stations used to define our reference frame. The first 7 are located in South-America, ISPA and GLPS are located on the Nazca plate. Velocities are indicated with respect to ITRF08. Others stations are located worldwide. Table S2: Nazca-South America relative angular velocities and respective predicted velocities on the Chilean trench at 28 • S.  Figure S1: Horizontal velocities (mm/yr) with respect to stable South-America a) defined in NNR-Nuvel-1A (DeMets et al. 1994); b) estimated in this study (20.9 • S; 128.6 • W; 0.122 • Myr −1 ). Red numbers indicate the velocity in mm/yr, red dots highlight the stations used to estimate the rotation vector (KOUR, BRFT/FORT, SAVO, TOPL, BRAZ, CHPI, CUIB), ellipses depict the region of 99% confidence Figure S2: Horizontal time series in ITRF of stations located on the brazilian craton that were used to estimate the rotation vector (Fig. S1). Figure S3: Horizontal time series in ITRF of A) cGPS stations COPO (Copiapó), UCNF (Antofagasta), TUCU (Argentina); B) campaign sites. All stations and sites are located on Fig. S4. Figure S4: Location map of permanent stations and campaign sites for which time series are presented on Fig. S3. Figure S5: Prediction of Klein et al. (2016) viscoelastic model of post-seismic deformation following the 2010 Maule earthquake: amplitude of horizontal velocity in the area of interest during the same period (2010)(2011)(2012)(2013)(2014)(2015)(2016). The dashed line highlights the limit north of which less than 1 mm/yr of postseismic signal is predicted. Figure S6: Left: Interseismic velocity field (mm/yr) relative to stable South America (same as in Fig. 1 of the main text). Right: Profiles of horizontal velocities (mm/yr) as function of distance from the trench (km) along 7 profiles. The corresponding profiles are depicted on the map. Figure S7: Checker board resolution tests. Left: synthetic checker board coupling distributions showing large scale (top, 80 × 30 km) or small scale (bottom 40 × 30 km) checkers. Centre: coupling distribution inverted using the synthetic velocities. Right: coupling distribution inverted using the synthetic velocities plus white noise Figure S8: Left: Best coupling distributions obtained when imposing the Andean sliver motion to the one proposed by Métois et al. (2016), with (bottom) or without (upper) allowing deeper than 60 km depth coupling. Right: associated residuals and nRMS. Figure S9: Left: Best coupling distributions obtained when fixing the Nazca-South America relative pole to the new value calculated in this study (53.9 • N, 86.1 • W, 0.605 • Myr −1 ). Right: associated residuals and nRMS. The coupling distribution is very similar to the best "model' presented in Figure 5 of the main text Figure S10: Decomposition of the modelled velocity field: adding the translation like motion associated to the Andean sliver rotation (left) to the deformation produced by elastic loading only (centre) gives the overall modelled velocity field (right). Figure S11: Sketch presenting the Andean rigid rotation inverted for three distinct regions of the margin: North Chile only Métois et al. (2013), Metropolitan and Atacama (Métois et al. 2014(Métois et al. , 2016, and Taltal area (this study). The inverted Euler poles are all located in South Atlantic, and the rotation produced over the margin is close to be parallel to the South America-Nazca convergence velocity. However, the amount of this sliver-induced motion significantly varies from North to South without any active structure able to accommodate these along-strike variations being identified. Therefore, the current study presenting new GPS data in the Taltal area challenges the hypothesis of an Andean belt rigidly deforming.
Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.