Continental block motion in the Northern Andes from GPS measurements

SUMMARY

Northwestern South America is a plate boundary zone where the Nazca, Caribbean and South American plates interact to produce a wide area of active continental deformation from the Gulf of Guayaquil (latitude 3°S) to Venezuela. Previous studies have identified a ∼2000-km-long continental sliver, referred as the North Andean Sliver (NAS), squeezed between the Nazca, Caribbean and South American plates and escaping at ∼1 cm yr⁻¹ northeastward with respect to South America. Subduction of the Nazca Plate beneath the NAS has produced a sequence of large and great earthquakes during the 20th century along the coast of Ecuador and Colombia. Large crustal earthquakes up to magnitude 7.7 have been documented along the proposed eastern boundary of the NAS. However, active tectonics data, historical and recent earthquakes all indicate active fault systems within the NAS, possibly resulting from the interaction of several tectonic blocks. Here, we derive an extensive horizontal velocity field using continuous and episodic GNSS data from 1994 to 2019.9, covering northern Peru, Ecuador, Colombia, Panama and Venezuela. We model the GNSS velocity field using a kinematic elastic block approach that simultaneously solves for rigid tectonic block rotations and interseismic coupling along the subduction interfaces and along major crustal faults. In contrast to previous results that considered a single rigid NAS, our dense GNSS velocity field demonstrates that the NAS undergoes significant internal deformation and cannot be modelled as single rigid block. We find that block kinematics in the northern Andes are well described by the rotation of 6 tectonic blocks, showing increasing eastward motion from south to north. The Eastern boundary of the sliver is defined by a right-lateral transpressive fault system accommodating 5.6–17 mm yr⁻¹ of motion. Fragmentation of the NAS occurs through several fault systems with slip rates of 2–4 mm yr⁻¹. Slow reverse motion is found across the sub-Andean domain in Ecuador and northern Peru at 2–4 mm yr⁻¹, marking a transitional area between the NAS and stable South America. In contrast, such a transitional sub-Andean domain does not exist in Colombia and western Venezuela. At the northwestern corner of Colombia, fast (∼15 mm yr⁻¹) eastward motion of the Panama block with respect to the NAS results in arc-continent collision. We propose that the Uramita fault and Eastern Panama Deformed Zone define the current Panama/NAS boundary, accommodating 6 and 15 mm yr⁻¹ of relative motion, respectively. A fraction of the Panama motion appears to transfer northeastward throughout the San Jacinto fold belt and as far east as longitude ∼75°W. Along the Caribbean coast, our model confirms, slow active subduction at ∼4.5 mm yr⁻¹ along the South Caribbean Deformed Belt offshore northern Colombia and a relatively uniform rate of ∼1–2 mm yr⁻¹ offshore northern Venezuela. Along the Nazca/NAS subduction interface, interseismic coupling shows a first-order correlations between highly locked patches and large past earthquake ruptures. These patches are separated by narrow zones of low/partial coupling where aseismic transients are observed. Compared to previous studies, our interseismic coupling model highlights the presence of deep coupling down to 70 km in Ecuador.

1 INTRODUCTION

Along the Ecuador–Colombia Pacific margin, the convergence of the oceanic Nazca Plate beneath the South American continental plate (SOAM) occurs at a rate ranging from 52 to 56 mm yr⁻¹ towards ∼N83°E, a direction which is 30–40° oblique to the trench perpendicular direction (Nocquet et al. 2014; Yepes et al. 2016; Jarrin et al. 2022). In many ocean–continent subduction settings, oblique convergence results in the development of an upper plate continental sliver and partitioning of the plate convergence vector into reverse slip at the subduction interface and trench parallel motion of the sliver, accommodated along strike-slip crustal faults within the overriding plate (e.g. McCaffrey 1992). Such oblique convergence partitioning is observed along the whole ∼6000 km length of the Nazca/SOAM Plate boundary (Nocquet et al. 2014). Indeed, in Chile south of latitude ∼38°S, the Liquiñe-Ofqui fault, a right-lateral strike-slip fault zone, delimits the Chiloe sliver moving northward at ∼6 mm yr⁻¹ with respect to SOAM (Wang et al. 2007; Moreno et al. 2008; Melnick et al. 2009). In central and northern Chile, several studies concurrently identified microblocks with motion up to ∼1 cm yr⁻¹ (Brooks et al. 2003; Chlieh et al. 2011; Métois et al. 2013). 5–6 mm yr⁻¹ of southeastward motion is observed throughout Peru and southernmost Ecuador for the Inca sliver (Nocquet et al. 2014), possibly fragmented into two separate domains (Villegas-Lanza et al. 2016b).

In the northern Andes, Pennington (1981) recognized the existence of a North Andean Block or Sliver (hereinafter as NAS) lying at the northwestern edge of the South American Plate, wedged between the Nazca, Caribbean and South American plates. The NAS is ∼2200 km long and 300–1000 km wide, comprising the Andean Cordillera and its margin in Ecuador north of the Gulf of Guayaquil, Colombia and western Venezuela (Figs. 1 and 2). Large megathrust earthquakes occurred along its western boundary in 1906 (M_w ∼ 8.6), 1942 (M_w 7.8), 1958 (M_w 7.7), 1979 (M_w 8.2) and 2016 (M_w 7.8; Kanamori & McNally 1982; Beck & Ruff 1984; Sewnson & Beck 1996; Nocquet et al. 2016; Ye et al. 2016, Fig. 1), with several smaller events (M_w ∼ 7.2) at the southern and northern edges of the 1906 earthquake rupture area in 1998, 1991 and 2004 (Fig. 1). In addition to megathrust events, the NAS hosts numerous crustal earthquakes (M_w ∼6.2–7.7; Fig. 1), among the largest documented in the Andes (Dimaté et al. 2005; Audemard et al. 2008; Choy et al. 2010; Beauval et al. 2010, 2013; Colón et al. 2015; SGC 2021).

GPS survey mode measurements performed during the 1980s and 1990s provided first-order constraints on the continental deformation in Ecuador and Colombia (Freymueller et al. 1993; Trenkamp et al. 2002; White et al. 2003). Trenkamp et al. (2002) quantified the northeastward escape of the NAS as a single rigid block moving at ∼6 ± 2 mm yr⁻¹ to the northeast. Analysis of 4 yr of a GPS network extending from central Peru to northern Ecuador provided an updated horizontal velocity field at the regional scale (Nocquet et al. 2014). In order to determine the NAS motion, Nocquet et al. (2014) noted that elastic strain induced by interseismic locking at the megathrust is negligible in Ecuador at latitude 1.5°S–3°S, south of the Carnegie ridge and north of the Gulf of Guayaquil (Fig. 2). Then, they found that these velocities are consistent at the 1 mm yr⁻¹ level with velocities in southern and central Colombia ∼300 km away from the trench, which they used to estimate an Euler pole at longitude −83.4°E, latitude 15.21°N, 0.29 deg Myr⁻¹, predicting northeastward (N70° ± 10E) 8.5 ± 1 mm yr⁻¹ motion with respect to South America. Their kinematic model for the NAS implies that ∼20 per cent of the trench-parallel motion is taken up by the sliver motion. More recently, Mora-Páez et al. (2019) used additional GPS data in northcentral Colombia north of latitude 4°N. Their proposed Euler pole (long. −185.2°E, lat. 58.6°N, 0.07 deg Myr⁻¹) predicts 8.6 mm yr⁻¹ towards N60°E relative to South America. Despite the large difference in location, both Euler poles predict similar velocities at latitudes 2–4°N. However, differences reaching 5 mm yr⁻¹ in northern Colombia show better predictions from Mora-Páez et al. (2019). Symithe et al. (2015) and Pérez et al. (2018) estimated kinematic models at the northernmost part of the sliver in northern Colombia and Venezuela, showing large differences with the Nocquet et al. (2014) or Mora-Páez et al. (2019) predictions for the northern part of the NAS. Such results suggest that the NAS kinematics are not well described using the assumption of a single idealized rigid block and calls for a study providing a model consistent at the regional scale.

Analysis of focal mechanism solutions, seismicity distribution, and tectonic studies suggest that the translational motion of the NAS is accommodated throughout several fault systems running along the Eastern Andean Cordillera and merging with the subduction zone offshore the Gulf of Guayaquil in Ecuador (Pennington 1981; Ego et al. 1995; Velandia et al. 2005; Tibaldi et al. 2007; Audemard et al. 2008; Audemard 2009; Egbue & Kellogg 2010; Alvarado et al. 2016; Fig. 2). This major boundary hosts most of the large crustal earthquakes documented in the northern Andes (Fig. 1). However, active crustal deformation within the NAS is indicated by tectonic studies (Audemard & Audemard 2002; Audemard 2014a; Alvarado et al. 2016) and both recent and historical seismicity (for example the 1868 M_w ∼7.3 Ibarra earthquake, the 1911 M_w ∼6.4 Yarumal earthquake and the 1885 M_w ∼ 6.4 El Tambo earthquake, Fig. 1; Dimaté et al. 2005; Beauval et al. 2010; SGC 2021). Focal mechanism solutions of moderate-sized earthquakes also attest to regular M_w >5 events occurring on active faults within the NAS (Ekström et al. 2012; Yepes et al. 2016; Vaca et al. 2019; Fig. 2). However, how much of the deformation is taken up by these faults and whether they are significantly less active than the faults delimiting the eastern boundary of the NAS remains largely unknown.

Understanding and modelling the NAS’s current motion and internal deformation requires first knowing the kinematic conditions along its boundary, provided by the motion of the major tectonic plates. Aside from the oblique convergence of the Nazca Plate described above, there are two additional boundaries to account for: (1) in northwestern Colombia, crustal earthquakes, neotectonic studies and previous geodetic results point out the collision of the Panama block against the NAS at ∼1 cm yr⁻¹ (Pennington 1981; Kellogg & Vega 1995; Kobayashi et al. 2014; Mora-Páez et al. 2019). The ongoing collision that started ∼15 Ma induces broad-scale deformation in northwestern Colombia (Taboada et al. 2000; Trenkamp et al. 2002; Kobayashi et al. 2014; Mora-Páez et al. 2019; Kellogg et al. 2019). (2) Northeast of the Panama block, the kinematic boundary acting along the northern edge of the NAS is provided by the eastward motion of the Caribbean Plate, at 17–18 mm yr⁻¹ with respect to South America (Symithe et al. 2015). The boundary between the NAS and the Caribbean Plate consists of an accretionary prism offshore western Venezuela and northern Colombia, referred to as the ‘South Caribbean Deformed Belt’ (SCDB: Fig. 2, Kroehler et al. 2011), overriding a low angle subduction interface (Kellogg & Bonini 1982; Hilst & Mann 1994; Taboada et al. 2000; Bernal-Olaya et al. 2015; Syracuse et al. 2016; Mora et al. 2017; Kellogg et al. 2019; Mora-Páez et al. 2019; Lizarazo et al. 2021). This subduction appears atypical with the absence of a developed magmatic arc and the lack of documented large megathrust earthquakes. Both the SCDB and the subduction interface are expected to accommodate shortening between the NAS and the Caribbean Plate. Recent elastic block models from Symithe et al. (2015) and Pérez et al. (2018) estimate Caribbean/NAS convergence at 5–9 mm yr⁻¹ offshore northern Colombia. However, these models neither include the Panama block motion nor consider the possible contribution of the Nazca subduction in the GPS velocity field, possibly resulting in a biased estimate.

To date, a model that accurately describes the present-day kinematics of the NAS based on geodetic observations is not available. This study presents an updated and extensive interseismic velocity field derived from continuous and episodic GPS measurements. The velocity field encompasses the North Andean Sliver and neighbouring regions, benefiting from (1) a longer time-series than available for previous studies and (2) new GPS sites in Ecuador, Colombia, Panama and Peru. Using the velocity field as input, we then build a kinematic elastic block model that simultaneously estimates rigid block rotations and spatially variable coupling at the subduction interfaces, providing crustal fault slip rates consistent with the derived kinematics. Our model is not constrained by using prior information on geological slip rates or focal mechanisms, but we use this information to validate our model. We finally discuss our results in terms of how tectonic plate motions are accommodated along active tectonic structures and draw some conclusions about the geodynamics governing the region.

2 ACTIVE TECTONIC SETTING

Here, we summarize the active seismo-tectonic information that is then used to build the geometry of our models described in Section 5. Fig. 1 shows the location of the major earthquakes mentioned in the following paragraphs, while Fig. 2 outlines the major active plate boundaries and tectonic features along the North Andean Sliver compiled from Audemard (1996), Audemard et al. (2000), Taboada et al. (2000), Alvarado et al. (2016), and the South American Risk Assessment (SARA) fault database (Alvarado et al. 2017), together with focal mechanisms solutions from Ekström et al. (2012), Vaca et al. (2019) and SGC (2021).

2.1 North Andean Sliver (NAS) eastern boundary

The eastern boundary of the NAS marks the transition from actively deforming areas to undeformed stable South America. Starting from the northeastern part of the NAS, the east–west trending El Pilar and San Sebastian strike-slip faults accommodate most of the relative motion between the Caribbean and South American plates (Pérez et al. 2001; Audemard 2009; Jouanne et al. 2011; Symithe et al. 2015). The San Sebastian fault marks the eastern limit between the northeasternmost part of the NAS and South America (Fig. 2). South of the San Sebastian fault, second-order subparallel structures such as the La Victoria fault are thought to accommodate very slow motion (Audemard et al. 2000; Pérez et al. 2018), which is negligible in our context. The westernmost extension of the San Sebastian fault merges with the NE–SW trending Boconó fault (BF), a predominantly right-lateral ∼500-km-long fault system, that also shows significant shortening. Such a shortening seems compatible with subparallel reverse faults flanking the Venezuelan Andes (Audemard et al. 2000; Audemard & Audemard 2002). The Boconó fault (BF, Fig. 2) appears to be the most active structure in the Mérida or Venezuelan Andes. Two large historical earthquakes with magnitudes estimated as large as M ∼7.5 occurred in 1610 and 1894 (Fig. 1) impacted several cities between Mérida and San Cristobal. They probably ruptured the southern portion of the Boconó fault (Schubert 1982; Audemard 1997, 2014b, 2016; Audemard et al. 2008).

South of the Boconó fault, the 1785 Cundinamarca and 1827 Altamira earthquakes (M ∼7; Dimaté et al. 2005; SGC 2021, Fig. 1) highlight the Guaycaramo (GCF), Algeciras (AF) and Afiladores-Sibundoy faults (ASF) as major active faults in the eastern Andes in Colombia (Chorowicz et al. 1996; Paris et al. 2000; Velandia et al. 2005; Acosta et al. 2007; Tibaldi et al. 2007; Diederix & Romero 2009; Diederix et al. 2021; Audemard 2021; Fig. 2). Several focal mechanism solutions (Fig. 2) support the view of a transition from transpression across the Guaycaramo fault to pure right-lateral motion along the Afiladores-Sibundoy fault.

Further south, the Afiladores-Sibundoy fault connects to the Chingual, Cosanga, Pallatanga and Puna faults in Ecuador (Fig. 2). These faults have been defined as the CCPP fault system by Alvarado et al. (2016) and Yepes et al. (2016). The Pallatanga and Puna faults delimit the boundary between the NAS and the Inca sliver, which encompasses the Andean Cordillera and its margin throughout Peru (Nocquet et al. 2014; Villegas-Lanza et al. 2016b). Although the Chingual, Cosanga and Pallatanga faults have been recognized as a transpressive right-lateral fault system, normal faulting occurs on secondary faults in the actively opening Gulf of Guayaquil (Eguez et al. 2003; Witt et al. 2006; Tibaldi et al. 2007; Alvarado et al. 2016). From historical records dating back 5 centuries, the ∼800-km-long CCPP fault has hosted large devastating crustal earthquakes as the 1987 M_w ∼7.1 Reventador earthquake that caused 1000 deaths and the 1797 M ∼7.6 earthquake, which destroyed the city of Riobamba and affected many towns in Chimborazo, Tungurahua and Cotopaxi provinces (Beauval et al. 2010, 2018; Alvarado et al. 2016; Fig. 1).

In most regions along the Andean orogen, fold-and-thrust belts have developed over tens or even hundreds of kilometres east of the Cordillera, referred as the sub-Andean domain (Baby et al. 1997; Beck & Ruff 2005; Brooks et al. 2011; Weiss et al. 2016). However, no sub-Andean domain is documented north of latitude 4°N. In Ecuador and Peru, active seismic exploration profiles reveal the presence of active thrust overlying a shallow (<10 km) continuous flat decollement as far as 600 km from the trench and 300 km from the eastern flank of Andean Cordillera (Baby et al. 1997). Active deformation within the sub-Andean domain is further witnessed by M_w ∼6 earthquakes in the Amazonia basin and for well developed morphology in the Eastern sub-Andean Belt (ESB), like the Napo-Cutucu fault and Macas fault. The latter fault hosted the M_w 7.0 1995 earthquake (Figs 1 and 2; Legrand et al. 2005).

2.2 Colombia

Additional active faults are documented within the area encompassed by the NAS. In northern Colombia, the main active structures are the Oca-Ancon (OAF) and the Santa Marta-Bucaramaga (SMB) faults (Fig. 2). The east–west trending Oca-Ancon prolongates to the west from the San Sebastian fault. The Oca-Ancon fault undergoes right-lateral strike-slip motion with reverse motion documented along the Oca segment (Audemard 1996). Crustal shallow earthquakes with M_w$<\sim$5 define this ∼600-km-long structure (Engdahl et al. 2020; SGC 2021; Fig. 4). However, evidence of Holocene activity with M ∼7 earthquakes is suggested by palaeoseismological studies (Audemard 1996). The OAF intersects the SSE-NNW Santa Marta-Bucaramanga fault at longitude 74°W. Although seismicity has remained low for the last decades along the Santa Marta-Bucaramanga fault (Vargas et al. 2019; Engdahl et al. 2020), palaeoseismological studies suggest that several significant earthquakes have occurred on the Bucaramanga segment since the Quaternary (Diederix et al. 2020). Both the Oca-Ancon and Santa Marta-Bucaramanga faults delimit the proposed Bonaire and Maracaibo blocks (Fig. 2), accommodating the differential motion of the northernmost part of the NAS with the Caribbean Plate and probably the Panama block (Audemard 1993, 2009, 2014a; Symithe et al. 2015).

In southcentral Colombia (latitudes 1°N to ∼7°N), the Romeral fault system (RFS, Fig. 2) is the main structure running along the western slope of the central Cordillera. Geometrical complexities arise southward of latitude ∼4°N where the fault changes from left-lateral strike-slip motion to right-lateral with a reverse component (Ego et al. 1995; Paris et al. 2000; Taboada et al. 2000), in agreement with the focal mechanism solutions displayed in Figs 1 and 2 at latitudes ∼3°S and ∼4.5°N (Ekström et al. 2012). Active faulting for this structure is evidenced by several M_w > 6 earthquakes that affected several departments in the Cauca river valley (Figs 1 and 2), like the 1994 M_w 6.8 Paéz earthquake and the 1999 M_w 6.1 Armenia earthquake which caused 921 deaths (SGC 2021).

2.3 Ecuador

Aside from the Chingual, Cosanga, Pallatanga, and Puna faults (CCPP) defining the eastern boundary of the NAS, additional active faults have been recognized within the Interandean Depression in central Ecuador, like the Latacunga fault (LF) and the Quito fault (QF; inset Fig. 2). Active deformation along the Quito fault is evidenced through several M_w > 5 shallow damaging earthquakes during the last decades and local InSAR/GNSS results showing a sharp velocity gradient across it (Mariniere et al. 2019). The Quito fault is a blind reverse fault, dipping west that includes minor dextral strike-slip faulting and extending along ∼60 km (Alvarado et al. 2014, 2016). Further north, the Otavalo-El Angel fault is considered as being the northernmost extension of the Quito fault. It hosted the 1868 M ∼7.2 earthquake that destroyed the city of Ibarra (Beauval et al. 2010). The Quito-Latacunga and the Cosanga faults are the eastern and western boundaries of an independent Quito-Latacunga microblock proposed by Alvarado et al. (2016).

In southern Ecuador (southward of latitude 4°S), the Macas fault follows the eastern flank of the Andean Cordillera, extending to northern Peru and merging with the transpressive Marañon fault system. The Marañon fault is a transpressional left-lateral strike-slip structure that follows the boundary between the Western and Eastern cordilleras. It probably extends further south in central Peru along the Chontal fault (Villegas-Lanza et al. 2016b). Although no large historical earthquake has been reported for the Marañon fault, M_w 4–5 shallow earthquakes witness of its current seismic activity (Engdahl et al. 2020; Ekström et al. 2012).

3 GPS DATA ANALYSIS

We use data from several regional continuous (cGPS) and episodic GPS data from a collaborative project between the French National Research Institute for Sustainable Development (IRD), the Geophysical Institute (IG-EPN; Mothes et al. 2013; Alvarado et al. 2018), the Military Geographical Institute of Ecuador (IGM 2021) and the Geophysical Institute of Peru (IGP). We integrate data of 21 cGPS sites from the GeoRED Project (Geodesia: Red de Estudios de Deformación) in Colombia, which is run by the Space Geodesy Research Group from the Colombian Geological Survey (SGC; Mora-Páez et al. 2018). We also include 180 cGPS sites from the National Geographical Institute of Peru (IGN), the Agustin Codazzi Geographical Institute of Colombia (IGAC 2021), the Tommy Guardia Geographical Institute of Panama (IGNTG 2021), the Argentine Network for Continuous Satellite Monitoring (Piñón et al. 2018), the Low-Latitude Ionosphere Sensor Network (LISN; Valladares & Chau 2012), the COCONet Project (UNAVCO Community 2008) and regional stations from the global network of the International GNSS Service for Geodynamics (IGS; Dow et al. 2009).

In addition, we integrate the GPS survey-mode data already used by Nocquet et al. (2014) but now benefiting from one to three additional epochs of measurement in Ecuador and Peru. The final data set includes 416 sites covering the 1994–2019.9 period, among which 123 are located within the NAS.

Continuous and episodic GPS data were homogeneously analysed using the GAMIT/GLOBK software release 10.71 (Herring et al. 2015, 2018), applying the methodology described in Jarrin et al. (2022), to derive time-series expressed in the international reference frame ITRF2014 (Altamimi et al. 2016) through its weekly updated realization from the IGS (Rebischung et al. 2016). We select cGPS benefiting from at least 3 yr of observations to mitigate the impact of seasonal variations (Blewitt & Lavallée 2002). For the final velocity field, we select 83 campaign sites with at least four periods of measurements spanning at least 4 yr. Daily position time-series are visually inspected in order to identify offsets induced by earthquakes or antenna changes. Most sites in Ecuador show significant post-seismic deformation following the 2016 M_w 7.8 Pedernales earthquake (Mothes et al. 2018; Rolandone et al. 2018). We therefore excluded post-Pedernales data for sites located within a radius of 750 km surrounding the Pedernales epicentre. Periods of slow-slip events occurring at the Ecuador-northern Peru subduction interface (Vallée et al. 2013; Segovia et al. 2015; Villegas-Lanza et al. 2016a; Vaca et al. 2018) were also removed and offsets estimated for them. For cGPS time-series, realistic velocity uncertainties are derived by estimating parameters for white and power-law noise, using the maximum likelihood estimator implemented in the CATS software (Williams 2008). In order to account for correlated noise in the velocity estimates for campaign sites, we add (quadratically) a coloured noise contribution taken from the median value of surrounding cGPS sites to the uncertainty obtained using classical least-square (white noise). Our obtained uncertainties are of the order of ∼0.3 mm yr⁻¹ at the 1σ confidence level for most cGNSS sites. For episodic GPS sites, velocity uncertainties are between 0.3 and 1.1 mm yr⁻¹ at the 1σ confidence level for 10-15 yr of observation.

In order to derive the most comprehensive velocity field at the scale of the North Andean Sliver, we combine our velocity field with the most recent published velocity field solutions from Mora-Páez et al. (2019), Pérez et al. (2018) and Mencin (2018) in a single consistent reference frame using common sites between solutions and applying the methodology proposed by Nocquet (2012). It is worth noting that the Pérez et al. (2018) solution is mostly derived from episodic GPS measurements in Venezuela, resulting in a significant number of velocities with uncertainty exceeding 3 mm yr⁻¹. In some cases, several nearby sites do not show coherent sense of motion with respect to the overall tectonic motion. Therefore, velocities with uncertainties greater than 2.0 mm yr⁻¹ and velocities showing opposite motion to nearby cGPS sites were rejected. In the case of significant velocity discrepancies at common sites between solutions, we keep the most recent estimate. Applying these criteria, our combination shows an average weighted root mean square (wrms) of 1.0 mm yr⁻¹, indicating a good agreement between our solution with the Mora-Páez et al. (2019), Pérez et al. (2018) and Mencin (2018) solutions. We also removed obvious outliers for which the velocity departs a few  mm yr⁻¹ from nearby sites.

Finally, the resulting velocity field is expressed with respect to the stable part of the South American Plate, by applying the SOAM/ITRF2014 Euler pole from Jarrin et al. (2022) (−133.28°E, 18.4°S, 0.121 deg Myr⁻¹), which shows a wrms of 0.15 mm yr⁻¹. Velocities are listed in the Supporting information (Table S1).

Compared to previous studies, our velocity field provides a better spatial sampling of the NAS, with an average intersite distance of 100–200 km in Colombia and ∼35 km in Ecuador. It also includes more sites in the Amazon basin, allowing us to test how far deformation spreads inside the South American Plate east of the Andes. Many Slow Slip Events (SSEs) have been documented at specific locations on the Ecuadorian coast (i.e. the La Plata island and the Punta Galera-Mompiche zone) and northern Peru (Vallée et al. 2013; Segovia et al. 2015; Villegas-Lanza et al. 2016a; Rolandone et al. 2018; Vaca et al. 2019), with recurrence times of a few years at certain locations (2–3 yr the La Plata island). Compared with the Nocquet et al. (2014) solution, eventually used in Chlieh et al. (2014), Villegas-Lanza et al. (2016b), and Mora-Páez et al. (2019), our updated solution is more robust to potential biases induced by undetected SSEs because: (1) several survey sites from Nocquet et al. (2014) have been progressively equipped with cGPS from 2011 onwards, allowing better detections of SSEs, (2) the survey sites used in our study have additional epochs of measurements, since only sGPS with at least 4 epochs of measurements with linear behaviour are kept and (3) cGPS have longer time-series against the ∼3 yr available in the Nocquet et al. (2014) solution for most sites. A possible influence of transient motion is further suggested by comparing our solution with respect to Nocquet et al. (2014) velocity field. The median bias (median of the differences) is less than 0.5 mm yr⁻¹ for sites located in the Amazonia and the Andean Cordillera but is 1.3 mm yr⁻¹ westward for sites located in the coastal region of Ecuador. The westward bias even exceeds 3 mm yr⁻¹ in southern coastal Ecuador. Search for potential transients during the 2009–2013 period is left for future studies. However, we note that the area where significant bias is found corresponds to the one impacted by the 2015 deep SSE described in Rolandone et al. (2018). Using the data selection criteria described above, our new velocity field for Ecuador can be considered as reflecting the steady interseismic regime between SSEs.

4 MAIN PATTERNS OF THE VELOCITY FIELD

The obtained regional velocity field with respect to Stable South America (SOAM; Fig. 3) shows several large-scale patterns we describe below. Across the Gulf of Guayaquil (latitude ∼3°S, Puna and Pallatanga faults), a clear velocity direction change occurs from N70°E to N130°E, highlighting the boundary between the NAS and the Inca Sliver (Nocquet et al. 2014; Villegas-Lanza et al. 2016b). With respect to the Nazca/SOAM convergence direction (N83°E), the velocity field is rotated by 10° counter-clockwise rotation in the southern part of the NAS and 40° clockwise in northern Peru. Although this rotation is coherent with the change of convergence obliquity induced by the change of the trench strike from the northern Andes to Peru, it is interesting to note that the direction change of inland GPS velocities also occurs in southern Ecuador at latitude ∼4°S, that is ∼200 km north of the apex of south America where the convergence obliquity changes. As noted in Nocquet et al. (2014) and Villegas-Lanza et al. (2016b), little if any elastic strain rate induced by locking at the megathrust occurs in that area.

Within the Peruvian Andes, velocities at latitude 6°S (LMAS, MOYB and YRMG, Fig. 3) depart from the overall 4–5 mm yr⁻¹ southeastward constant motion within the Inca sliver to 2–3 mm yr⁻¹ eastward dominant motion. These velocities therefore witness internal deformation of the Inca sliver. The relative motion of LMAS and MOYB with respect to YRMG imply ∼2 mm yr⁻¹ of shortening accommodated through folding and thrusting in the sub-Andean domain (Villegas-Lanza et al. 2016b).

In central and northern Ecuador (latitudes 2°S–1°N), velocities rapidly decrease from 20 to 25 mm yr⁻¹ at the coastline down to 6–7 mm yr⁻¹ at ∼180 km inland in the eastern Cordillera. This systematic decrease defines a velocity gradient that reflects significant locking on the megathrust interface (Chlieh et al. 2014; Nocquet et al. 2014, 2016; Jarrin et al. 2016; Gombert et al. 2018; Chlieh et al. 2021). East of the Chingual and Cosanga faults, the direction of velocities within the subandean domain is similar to those observed in the sub-Andean domain of northern Peru (Fig. S1), suggesting ∼2 mm yr⁻¹ of shortening accommodated across the Eastern subandean Belt.

Along the Pacific coast of Colombia, the few available cGPS sites (TUCO, GUAP, BUGT, latitude 2–4°N) are located at almost the same distance from the trench. Previous kinematic models for the NAS motion show that the rigid block contribution is the same within 0.1 mm yr⁻¹ for those sites (Mora-Páez et al. 2019), so that any change in their velocity magnitude can be used as a proxy for lateral changes of interseismic coupling at the megathrust. For example, GUAP and BUGT show velocity magnitudes 35 and 45 per cent lower than TUCO, indicating that interseismic coupling decreases from south to north. To illustrate this point, insets A and B in Fig. 3 show two profiles perpendicular to the trench and a best-fitting 2-D back-slip model prediction (Savage 1983) for a ∼23° constant dipping with homogeneous locking. Both profiles show low-velocity gradients (3 mm yr⁻¹), consistent with ∼50 per cent partial interseismic coupling at the subduction interface, down to ∼40 and ∼35 km depth for profile A and B, respectively.

The overall velocity field within the NAS shows a progressive counter-clockwise rotation from latitude 0° in Ecuador to latitude 5.5°N in Colombia and then a clockwise rotation from latitude 5.5°N to latitude 12°N. This pattern demonstrates that neither the Euler pole from Nocquet et al. (2014) predicting a progressive counter-clockwise rotation to the north nor the Euler pole from Mora-Páez et al. (2019) predicting a progressive clockwise rotation to the north can correctly describe the kinematics of the NAS (see Fig. 7 for a comparison). Indeed, fitting such a pattern requires either at least two blocks and/or internal deformation within the sliver.

Where available, GPS sites located east of the Andean Cordillera indicate different levels of residual velocity with respect to stable SOAM. In Peru at latitude 6°S, site YRMG shows negligible velocity (∼1.5 mm yr⁻¹) in a SOAM fixed frame, indicating that all shortening is accommodated across the westernmost part of the sub-Andean domain there. On the contrary in Ecuador, velocities at latitude 0.5°S show systematic 1–2 mm yr⁻¹ eastward motion, as far as NORE located ∼250 km east of the Andean Cordillera foothills in Ecuador. In Colombia, sites VIVI, OCEL, VSJG surrounding latitude 4°N appear to be part of stable SOAM. Together with the very small velocity at CAPI (latitude ∼5.3°N), all of them indicate very little deformation east of the Eastern Cordillera. A similar observation can be made east of the Boconó fault in Venezuela, where GPS sites (ARCA, BARI, CN41 and GUAC) appear to be part of the stable SOAM (∼0.5 mm yr⁻¹). Thus, the existence of an actively deforming sub-Andean domain appears restricted to Peru and Ecuador.

Nine cGPS velocities in Panama show a 22–24 mm yr⁻¹ eastward motion for the Panama block with respect to SOAM. The magnitude and direction of VORA and APTO velocities (located in the Atrato region, in Colombia, Fig. 3) suggest that both sites also belong to the Panama block. Close to the Uramita fault and San Jacinto fold belt, 5 GPS sites (SINC, CORO, BACO, BARU and CART) are moving at 17–18 mm yr⁻¹ with the same orientation as the observed velocities within Panama. On the contrary, velocities northward of latitude 10.5°N experience a slight counter-clockwise rotation with respect to the direction of the Caribbean motion. Therefore, the velocity field in northernmost Colombia and western Venezuela appears to undergo the effects of the Panama block collision and the Caribbean Plate subduction, in addition to the NAS translation (Trenkamp et al. 2002; Mora-Páez et al. 2019).

5 MODEL SETUP

5.1 Modelling approach

We model the velocity field expressed in a South America fixed reference frame by assuming that the horizontal velocity at any GPS site results from a rotation of a rigid block hosting the site and an elastic contribution induced by interseismically locked patches at faults. Such an approach called ‘elastic block modelling’ has been implemented by McCaffrey (2002), Meade & Hager (2005) and Meade & Loveless (2009). Here, we use our own implementation (Nocquet, in preparation), which follows the linear formulation described in Meade & Loveless (2009) with two minor differences: (1) the matrix relating the unit slip at a triangular dislocation to the surface displacement is built using the artefact-free formulation of Nikkhoo & Walter (2015) instead of one from Meade (2007). (2) While Meade & Loveless (2009) use Laplacian-like regularization constraints for solving spatially variable slip-deficit at subduction interface, we use a regularization constraint through a decreasing exponential model variance–covariance matrix that simultaneously imposes some level of smoothness and minimizes the departure from a prior model m₀, as in Radiguet et al. (2011) and Nocquet et al. (2014). This approach allows further exploration of the range of models allowed by the data by varying the prior model from null coupling (fully creeping) to the fully coupled interface (Nocquet et al. 2014). Furthermore, this approach is useful to test whether deep coupling induces a trade-off with the estimated block rotations by varying regularization parameters in the inversion. In our block model, the mentioned approach for estimating the spatial distribution of coupling is performed for (1) the Nazca subduction interface in northern Peru, Ecuador and southcentral Colombia and (2) the Caribbean subduction interface in northern Colombia.

5.2 Subduction interface geometry

The Slab2.0 Nazca/SOAM subduction interface (Hayes et al. 2018) is discretized into 1021 triangular dislocations of 30-km-long edges down to 80 km depth. For the La Plata area in central Ecuador (latitude 0.2S°–3.5°S), we modified the slab2.0 contours by including contours derived from marine seismic profiles and hypocentral solutions of microseismicity (Collot et al. 2017; Segovia et al. 2018; Font et al. 2019). Our triangular grid, therefore, extends southward as far as latitude ∼13°S in order to re-evaluate the motion of the Inca sliver in northern and central Peru. For the Caribbean/SOAM subduction interface westward of longitude 73°W, the geometry is derived from hypocentral solutions from global (Engdahl et al. 2020) and local (Colombian Geological Survey, SGC 2021) earthquake catalogues between 1970 and 2019. Revisiting the shallow seismicity distribution along AA’ and BB’ cross-sections shown in Fig. 4, we define a subduction interface with an average dip of 11° down to ∼40 km depth. We discretize this subduction interface into 25 km × 20 km rectangular dislocations with a constant dip of 11°. Eastward of longitude 73°W, the Caribbean/SOAM subduction interface is discretized using 100 km × 60 km rectangular dislocations dipping at 13° down to 30 km (Schmitz et al. 2008; Pérez et al. 2018). We use dipping fault planes at 20° and 15° for the western and eastern portions of the North Panama Deformed Belt (NPDB; Kobayashi et al. 2014). The South Panama Deformed Belt (SPDB) is constrained as single vertical faults locked down to 13 km based on the average hypocentres provided by the ISC reviewed earthquake catalogue (Engdahl et al. 2020).

5.3 Crustal fault geometry

Crustal block boundaries are modelled using rectangular fault planes with constant dip, assumed to be locked from the surface down to a prescribed depth. Fault geometry parameters are assigned from published information derived from geological studies where available. In the case of no published information, we use vertical faults locked down to 15 km depth as a starting value. When GPS data are sparse, we impose the locking depth to be equal for adjacent fault segments (Elliot et al. 2010). Iterative inversions were run to adjust fault parameters that provide the lowest velocity residuals (observed minus modelled). Similarly to Symithe et al. (2015), the best overall locking depth is found to be 14 ± 1 km for crustal faults bounding the eastern NAS boundary (Fig. 5). However, crustal faults within and outside the NAS are best modelled by prescribing significantly lower locking depths ranging from 3 down to 14 km, which suggest the presence of creeping or partial creeping on some active structures (e.g. the Cosanga and Quito faults). Crustal fault parameters are available in the Supporting Information (Text S1 and Table S2) for the best-fitting model (model D) discussed in Section 6.1.

5.4 Neighbouring plate and block motion

In all inversions, we fix the relative motion of the Caribbean/SOAM to the values published by Symithe et al. (2015) and use the values of Jarrin et al. (2022) for the NAZCA/SOAM relative motion. Both Euler poles were carefully determined from both plate scale studies and our data set is not expected to provide any additional contribution. For the PANAMA/SOAM and INCA/SOAM relative motions, we use values from Kobayashi et al. (2014) and Villegas-Lanza et al. (2016b) as prior, but allow them to be estimated in the inversion. We make this choice because our data set includes sites at those two block boundaries that might require adjustment of previously published values.

5.5 Block model geometry and model selection

Defining the geometry of a block model remains, at least partly, a subjective choice. We choose to move from the simplest model consisting of a single NAS block and progressively increase the complexity by introducing step-by-step faults delimiting new blocks. For that, we use the available information about mapped fault traces from active tectonic studies and seismicity distribution described in Section 2, together with the velocity gradients observed in the GPS data. In order to assess whether additional complexity is required by the data, we use the classical F_ratio test (eq. A3). The F_ratio test quantifies whether the decrease of chi-square (χ²) associated with adding more adjustable parameters, defined as the weighted quadratic sum of model residuals obtained when adding one or several blocks, is statistically significant (Stein & Gordon 1984; Nocquet et al. 2001). Available Holocene slip rates for faults and slip vectors from focal mechanism solutions are not included in our inversions. They are used for a posteriori validation of our best model described in Section 6. Fig. 6 shows our selection of the most relevant geometries.

5.6 Results for increasing complex models

Model A (Fig. 6a) considers the North Andean Sliver (NAS) as a single block outlined by large-scale tectonic plate boundaries with the Nazca, South American and Caribbean plates, and the Panama block (Audemard et al. 2000; Taboada et al. 2000; Bird 2003; Machare et al. 2003; Alvarado et al. 2016). North of latitude ∼7.0°N, the Panama block collides with the NAS producing a relatively broad deformation zone (Pennington 1981; Audemard 1993; Kellogg & Vega 1995; Trenkamp et al. 2002; Mora-Páez et al. 2019). As a first attempt, we define the Panama/NAS boundary following the trend of shallow seismicity clusters (see Fig. 4) as far as the southernmost extension of the Uramita fault southward of VORA and APTO GPS sites (see Fig. 3). At the southern tip of the NAS, we follow the Peruvian Forearc (PF) boundary proposed by Villegas-Lanza et al. (2016b) along the Marañon river in Peru, but choose to extend it to the sub-Andean domain along the Macas reverse fault, merging to the Cosanga fault, instead of cutting the Andean Cordillera in southern Ecuador (Fig. 2). This choice is grounded in the well-documented compressive Macas fault (Legrand et al. 2005) and the absence of known large active faults or crustal seismicity in the western Cordillera in southern Ecuador. As expected from the description of the general pattern of the velocity field, Model A cannot satisfactorily fit the data. Despite a relatively low average residuals (wrms = 1.90 mm yr⁻¹), larger velocity residuals are observed in several areas of the sliver, especially in northwestern Colombia and Venezuela, reaching up to ∼7 mm yr⁻¹ (velocity residuals available in Fig. S3).

In northern Colombia, clusters of shallow seismicity are observed along the Oca-Ancon and Santa Marta-Bucaramaga faults (Fig. 4), and geological studies concluded that these two fault systems delimit an independent Maracaibo block (Audemard 1993, 2009, 2014a). Adding the Oca-Ancon fault allows us to delimit an independent Bonaire block between the Maracaibo block and the Caribbean Plate (Fig. 2). A second area of obvious misfit is the Interandean Depression in Ecuador, where ∼3 mm yr⁻¹ of very localized east–west shortening (Mariniere et al. 2019) suggests the existence of an additional block outlined by the Latacunga, Quito and El Angel faults (Alvarado et al. 2016, Fig. 2). Our model B therefore includes these 3 additional blocks to model A. Model B (Fig. 6b) shows significant fit improvement (χ² = 6116 and variance reduction of 98.4 per cent). The F_ratio test with respect to model A is positive at the 99 per cent confidence level (Table S3). However, 3–7 mm yr⁻¹ velocity misfits are observed south of the Santa Marta-Bucaramanga fault (SMB) at several sites surrounding the Panama/Colombia boundary and the Guaycaramo fault (Fig. S3). Similarly, 2–3 mm yr⁻¹ of residuals persist within the sub-Andean region of Ecuador and northern Peru. In addition, model B predicts ∼6 mm yr⁻¹ of opening along the Puna and the southernmost extension of Pallatanga fault, at odds with the dominant strike-slip faulting inferred by geological studies (Winter et al. 1993; Baize et al. 2015, 2020; Alvarado et al. 2016).

In an attempt to solve these problems, we note that GPS velocities suggest ∼2 mm yr⁻¹ of shortening accommodated throughout the Romeral fault system in central Colombia. Furthermore, active faulting along this ∼1000-km-long structure is confirmed by the 1994 and 1999 shallow earthquakes (Section 2; Taboada et al. 2000; Ekström et al. 2012). To define a new block, we merge the Romeral fault to the north with the southernmost extension of the Uramita fault and southward with the El Angel fault in northern Ecuador. For the sub-Andean region of Ecuador and northern Peru, the velocity field suggests 2–3 mm yr⁻¹ of shortening accommodated by the Eastern Subandean Belt (Fig. 3). Focal mechanisms also confirms active reverse faulting in that area (Fig. 2). Therefore, we define the EPSUB block in order to model the sub-Andean zone following the Eastern Subandean Belt boundary from southernmost part of the Andean Cordillera in Colombia to northern Peru. Including these new blocks in model C (Fig. 6c), velocity residuals improve at ∼1 mm yr⁻¹ within the sub-Andean and at several sites close to the Guaycaramo fault (Fig. S3). Model C variance reduction is 99.1 per cent, and the improvement of chi-square decrease is significant, well above the 99 per cent confidence level (χ² = 2780.4). The unrealistic opening along the Puna and Pallatanga faults is now reduced to ∼1.5 mm yr⁻¹, making this fault system predominantly right-lateral strike-slip in agreement with tectonic studies. However, we find two remaining issues for model C: (1) large velocity residuals (4–5 mm yr⁻¹) are still observed at several GPS sites (BACO, BARU, CART, CORO, SINC and URRA: Fig. 3) close to the Panama/Colombia border (Fig. S3), requiring additional complexity in that area and (2) model C predicts 6–7 mm yr⁻¹ of shortening along the southernmost extension of the Romeral fault between latitudes ∼1°N and ∼4°N, at odd with the geological study from Paris et al. (2000) which proposes slip rate of the order ∼1 mm yr⁻¹.

In order to improve the fit to GPS data close to the Panama/Colombia boundary, we follow the proposition from several studies in which the San Jacinto Fold Belt (SJFB) is the northward prolongation of the Romeral fault merging with the westernmost extension of the Oca-Ancon fault surrounding the Santa Marta city (latitude ∼11°N). However, the absence of surface fault traces and low seismicity make it difficult to image the northernmost segments (Taboada et al. 2000; Cedriel et al. 2003; Mora et al. 2017; Kellogg et al. 2019). By contrast, GPS sites (BQLA, VPOL and CN37, Fig. S3) at latitudes 10.6°N–11°N show velocities consistent with sites located in central Colombia east of the Romeral fault. We therefore propose a boundary merging the northernmost extension of the San Jacinto Fold Belt (SJFB) with the Caribbean Plate boundary at latitude ∼10.5°N, hence defining a new San Jacinto (SANJ) block (Fig. 6d). On the other hand, the low slip rate proposed by Paris et al. (2000) for the southern segments of Romeral fault suggests that it accommodates only a fraction of the deformation within the Andean Cordillera in southern Colombia. Furthermore, two sites POPA and POVA show velocities consistent with the kinematics of the westernmost block in Ecuador and Colombia (Figs 3 and S3). Both observations suggest that some deformation should also take place east of the southern Romeral fault, but the scarce GPS data in that part of Colombia prevent us to determine an even first-order location for that boundary. In our model D, this is technically achieved by attributing sites west and east of the Romeral fault to two different blocks without accounting for elastic strain for that part of the boundary separating them. Finally, the geometry for our final model is shown in Fig. 6(d).

6 BEST-FITTING MODEL DESCRIPTION AND UNCERTAINTY

6.1 Best-fitting model geometry summary

Our best-fitting model includes 11 blocks that rotate with respect to South America. The SOAM, NAZCA and Caribbean plates define kinematic boundary conditions at the model edges. The Panama block and Inca sliver are treated as intermediate boundary conditions, taking previously published Euler poles from Kobayashi et al. (2014) and Villegas-Lanza et al. (2016b), but allowing a small adjustment to their published values. For the Panama block, our pole is now at (lon: −82.69°E, lat: 43.65°N, 0.346 deg Myr⁻¹), almost identical to Kobayashi et al. (2014) (lon: −88.358°E, lat: 43.447°N, 0.364 deg Myr⁻¹) but predicting ∼ 1 mm yr⁻¹ slower velocity. Another difference with Kobayashi et al. (2014) is that we do not find the need for an additional Chocó block south of the Panama/Colombia border. Instead, an intermediate San Jacinto (SANJ) block is required to explain GPS velocities north of the Uramita fault (Figs 8 and S3). Although our new estimate for the Inca sliver is close to that of Villegas-Lanza et al. (2016b), sites in the southern Ecuador Cordillera indicate a different boundary delimiting the proposed Peruvian forearc sliver.

The kinematics within the NAS can be modelled by the rotation of six blocks, accommodating their relative motion across major faults or active structures identified by previous active tectonics studies. Two block boundaries appear to be required by the geodetic data, despite a lack of seismological or geological observational evidence for their activity: a ∼150 km segment at the northern tip of the San Jacinto Fold Belt in northern Colombia and a ∼200 km segment south of latitude 4°N in southern Colombia. Both are indicated by dotted lines from Figs 8 to 17.

An additional complexity to previous models is also observed south of latitude 2°N, where contractional deformation occurs within the sub-Andean domain, possibly penetrating the SOAM Plate over a few hundred kilometres in Ecuador and northern Peru.

6.2 Euler poles and comparison with previous results

Table 1 summarizes Euler poles with respect to South America and the associated uncertainties for our preferred model. Figs 7 and S4(d) show velocities at selected sites predicted from rigid rotation contribution only (no elastic strain contribution). The two smallest blocks (UIOL and SANJ) within the NAS have the largest uncertainties, which are directly related to their size and the reduced coverage of GPS sites within them.

Predicted velocities within the North Andean Sliver (NAS) are the fastest within the BONA block. For this block, the Euler pole predicts velocities at ∼17.5 mm yr⁻¹ towards N83.5°±1.5E, similar to the 17–18 mm yr⁻¹ towards N87.5°±2.5 predicted by the Caribbean/SOAM relative motion along the northernmost boundary of the NAS (the South Caribbean Deformed Belt). Such agreement in velocity predictions suggests a strong influence of the Caribbean subduction in that part of the NAS. By contrast, the pole of the UIOL block predicts the lowest velocities throughout the NAS at ∼5.5 mm yr⁻¹ towards N59±7°E. The MARA block undergoes fast clockwise rotations. Velocity azimuths increase from west to east by 15° (N70°–85°E) and velocity magnitudes increase from 10.6 mm yr⁻¹ at latitude 8°N to 13 mm yr⁻¹ at latitude ∼10°N (Figs 7 and S4d). Inside the SANJ block, the Euler pole predicts velocities at 17 ± 0.2 mm yr⁻¹ towards N87°E, almost the same direction as the Panama block (22 mm yr⁻¹ towards N86°E), hence suggesting a predominant control of the Panama collision in that part of the NAS. Along the NAW block, our model prediction is 8.4 ± 0.1 mm yr⁻¹ towards N68°±1E as rigid motion, consistent with the average predicted rigid motion from simple models reported by Mora-Páez et al. (2019) (8.0 mm yr⁻¹ towards N59°E) and Nocquet et al. (2014) (8.5 ± 1 mm yr⁻¹ towards N70°±10E). Despite this agreement, we find discrepancies in the predicted rigid motion from Nocquet et al. (2014) in southcentral Ecuador (latitudes 3°S–0.5°N) at 1.5–2.0 mm yr⁻¹ in the east and north velocity components, and only ∼1 and 0.4 mm yr⁻¹ in the east and north velocity components with respect to the prediction from the model of Mora-Páez et al. (2019).

Regarding the Inca sliver, our predicted velocities for the Peruvian Forearc West block (PFW) in southern Ecuador and northern Peru are very similar both in magnitude and direction (4.2 ± 0.3 mm yr⁻¹ towards N125°±2E) to those proposed by Villegas-Lanza et al. (2016b) (4.5 ± 0.5 mm yr⁻¹ towards N125°±5E).

6.3 Model fit and uncertainties

Our preferred model (Fig. 6d) shows post-fit residual statistics significantly better (Table S3) than previous simpler models that included a smaller number of blocks (model A–C in Fig. 6). For our preferred model, the data variance reduction is 99.6 per cent, the wrms is 0.99 mm yr⁻¹, and the reduced χ² = 2.07, indicating a good agreement between residual magnitude and data uncertainties. The spatial distribution of velocity residuals does not indicate any systematic pattern or area suffering obvious mismodelling (Fig. 8). Indeed, the histograms of velocity residuals follow a normal distribution in which most are lower than 1.0 mm yr⁻¹ (see Fig. 8b). Compared to simpler models, the F_ratio test shows significance above the 99 per cent confidence level, suggesting an appropriate division of blocks (Table S3).

Although our model provides an average fit of ∼1 mm yr⁻¹ to the observed GNSS velocities, statistically significant residuals, exceeding 2 mm yr⁻¹, are still observed at several GPS sites (BASO, BACO, DAR2, VORA and CASI; Fig. S3). All these sites (except BASO) are located within a broad deformation zone caused by the collision between the Panama block and the NAS. There, deformation might be more diffuse, involving additional secondary faults not accounted for in our model (Duque-Caro 1990; Paris et al. 2000; Taboada et al. 2000).

Assessing the uncertainty and validating a kinematic model is difficult because the choice of the boundaries is subjective. However, the uncertainty can be provided assuming that our choice of block boundaries and locking depths is correct, and we detail it below. Since our approach uses a linear formulation, the posterior variance–covariance matrix provides formal estimates of the uncertainties. The full variance–covariance matrix is provided as Supplementary Data (Table S4). Various elements of the variance–covariance matrix are useful to look at. First, the full subvariance–covariance matrix includes information about the precision of each pole. Using the propagation law of covariance, an error ellipse can be provided for the predicted motion at the centroid of each block. Fig. 9 shows that each block has its kinematics formally determined at the ∼0.5 mm yr⁻¹ to 95 per cent confidence level.

We further address potential trade-off between the interseismic locking at the Carribean and Nazca subduction interfaces and the derived continental kinematics. Indeed, high coupling at depth might induce surface motion with low-velocity gradient, similar to a translation that could bias the estimates of block rotations. To assess this potential impact, we predict the motion at several sites within the North Andean West (NAW) block (Fig. 6d) using Euler pole estimates from a range of interseismic coupling (ISC) models allowed by the GPS data.

This exploration of ISC models is performed by varying the weight given to damping parameter (σ_m) with respect to a priori model m₀ (regularization parameters and inversion approach described in the Supporting Information). Fig. 10(a) shows the L-curve (χ² values versus σ_m constraints) for the dip-slip component with a priori model m₀ = 0 (null coupling). The best ISC model if found for σ_m = 20 mm yr⁻¹. However, ISC models with a σ_m between 8 and 35 mm yr⁻¹ are also possible (Text S4 and Fig. S7). Predicted velocity magnitudes (Fig. 10b) show negligible discrepancies for ISC models with σ_m values of 8, 20 and 35 mm yr⁻¹ (for example 0.2 mm yr⁻¹ southward of latitude 2°N for the ISC model with σ_m = 8). Slight variations in predicted azimuths that do not exceed 2° (Fig. 10c) are observed for the ISC model with a σ_m = 8 mm yr⁻¹ between latitude 2°N and 2°S, where the coupling distribution is high (∼80 per cent). The same test of velocity predictions was performed for the UOIL and ROME blocks, providing similar results. In conclusion, we find negligible biases in the block rotations coming from the trade-off between interseismic coupling and block rotations, attesting to Euler poles’ stability for the different blocks forming the NAS.

6.4 Comparison with earthquake slip vectors

Where available, earthquake slip vectors provide information, independent of geodesy, about the direction of the relative motion of blocks surrounding a fault. Earthquake slip vectors can therefore be used to validate or point out problems in our model. We use focal mechanism solutions from global (Ekström et al. 2012) and local (Vaca et al. 2019; SGC 2021) catalogues for the NAS and neighbouring regions. We select earthquakes with M_w > 5, shallower than 30 km for crustal faults, and shallower than 50 km for megathrust earthquakes. We exclude aftershock solutions following the 2016 Pedernales earthquake (M_w 7.8) to avoid an over representation of central Ecuador in the statistics. Applying these criteria leaves 91 focal mechanism solutions, listed in Table S5. The surface projection of earthquake slip vectors is then compared to the direction of relative motion predicted by our model for pairs of blocks calculated at the location of the epicentre. Fig. 11 shows relative block directions inferred from focal mechanism solutions together with the prediction from our best-fitting model (model D).

In general, we note a good consistency between earthquake slip vectors and the relative motions predicted by our model. Statistical analysis of focal mechanism solutions worldwide suggests that uncertainties of ∼15° in the axial orientation of tensor solutions must be considered (Frohlich & Davis et al. 1999). Assessing the relative motion at block boundaries, we note a very good agreement between the geodetic model prediction and focal mechanisms for crustal faults (inset b in Fig. 11 indicated by a mean bias 1° and a standard deviation of 12°). At the Ecuador–Colombia and northern Peru subduction zones, earthquake slip vector directions agree within ∼2.5° (bias = −3° and standard deviation of 3.2°) to the predicted NAZCA/NAW and NAZCA/PFW relative motions (inset a within Fig. 11), validating our model estimates for the trench parallel motion. The agreement between the earthquake slip vectors and our model prediction is excellent along the eastern boundary of the NAS and in the sub-Andean domain. It is also surprisingly good along our proposed Panama/NAS boundary (PANA, SANJ, NAW and CARI blocks) supporting our choice for the block geometry in that area. Areas with lower consistency (20–30°) are found along the Oca-Ancon fault for the BONA/MARA relative motion and the central Romeral fault.

7 COMPARISON AGAINST ACTIVE TECTONIC STUDIES AND GEOLOGICAL SLIP RATES

In this section, we compare the prediction of our best kinematic block model with available results from active tectonic studies in terms of Quaternary or Holocene fault slip rates and style of faulting. Fig. 12 summarizes the slip rates predicted from our preferred model together with geological slip rates available from the literature. In the Supporting Information, geodetic and geological (strike-slip and dip-slip) slip rates are available in Table S2.

7.1 Eastern boundary of the NAS in Ecuador, Colombia and Venezuela

At the southernmost boundary of the NAS, our best model estimates right-lateral slip rates at 7.5 ± 0.1 mm yr⁻¹ along the Puna fault and the southern section of Pallatanga fault, resulting from the relative motion between the North Andean West (NAW) and Peruvian Forearc West (PFW) blocks. This estimate agrees with the range of 5.8–8 mm yr⁻¹ proposed by Dumont et al. (2005) from river channel offsets at the Puna island (at lon. 80.23°W, lat. 3°S) since the mid-Holocene (5–6 kyr). In addition to the predominantly strike-slip motion, our model also predicts fault-normal opening rates across the Puna and southern Pallatanga faults at 1.5 ± 0.3 mm yr⁻¹. From the thickness of sediment deposited in the Gulf of Guayaquil since 5.3 Myr, Lavenu et al. (1995) estimated an opening rate normal to the Puna fault of 2.4 ± 1.1 mm yr⁻¹, consistent with our model prediction. Available seismic reflection profiles offshore the Gulf of Guayaquil have also identified normal faulting delimiting subsiding basins (Witt et al. 2006).

The northern segment of the Pallatanga fault delimits the boundary between the Peruvian Forearc West and the Quito-Latacunga block, whose motion is slower than the western NAS. Slip rates along that section are therefore expected to be slower than the southern Pallatanga and Puna faults. Our estimate is right-lateral slip at 5.6 ± 0.2 mm yr⁻¹. Baize et al. (2020) summarized slip rates along a 90-km-long segment of the Pallatanga fault from the Rumipamba section (lon. 78.9°E, lat. 1.8°S) to the Pisayambo section (lon. 78.3°E, lat. 1.1°S) from dated offsets of morphological markers and lava flow. Their proposed slip rates range from 2 to 6 mm yr⁻¹, compatible with our model prediction. However, in detail, their new Holocene slip rate estimates for the Rumipamba section are 2 mm yr⁻¹ slower than previous estimates of 2.9–4.6 mm yr⁻¹ proposed by Winter et al. (1993) and Baize et al. (2015). At this location, our estimate is ∼7 mm yr⁻¹ if we consider that this fault segment is part of the NAW/PFW boundary and ∼5 mm yr⁻¹ if it is part of the UIOL/PFW (Quito-Latacunga/Peruvian Forearc west) boundary. Igualata volcano Holocene lava flow offsets (lon. −78.3, lat. −1.5) indicate slip rates of 4–6 mm yr⁻¹, in agreement with our estimates. In addition, Baize et al. (2020) also propose slip rates of 3–4 mm yr⁻¹ for the northernmost Pisayambo segment that would be 1–2 mm yr⁻¹ slower than our model prediction. In summary, our model predictions agree with the upper range of Holocene estimates in that area.

The Cosanga and Chingual faults define the boundary between the sub-Andean domain and the Quito-Latacunga block from latitude 1°S to 1°N. Our model predicts right-lateral strike-slip rates at 6 ± 0.2 mm yr⁻¹ and 6.3 ± 0.1 mm yr⁻¹ along the Cosanga and Chingual faults respectively, with a minor reverse component of 2.2 ± 0.2 mm yr⁻¹ for the Cosanga fault vanishing along the Chingual fault. Geological estimates for the Cosanga fault are not available. Further north, average slip rates (7 ± 3 mm yr⁻¹) estimated from volcanic deposits dated at 37 and 8.6 kyr from the Soche volcano at the southwestern section of the Chingual fault (Ego et al. 1995) are consistent with our model prediction.

Dominant right-lateral motion at 8.7 mm yr⁻¹ along the Algeciras-Sibundoy faults agree with late Pleistocene–Holocene rates (7.7–11.9 mm yr⁻¹) based on geomorphological studies (Tibaldi et al. 2007). Further north, our model finds a transition from dominant right-lateral strike-slip motion to transpression along the Guaycaramo fault. We estimate right-lateral strike-slip rates decreasing from 6.2 mm yr⁻¹ at latitude ∼3°N to 5.8 mm yr⁻¹ at latitude ∼6°N, while reverse slip rates increase from 2.1 to 4.7 mm yr⁻¹. These values are compatible with the upper bounds of average Quaternary right-lateral motion (3 ± 2 mm yr⁻¹) from landforms and offsets of alluvial terraces reported by Paris et al. (2000). Despite the scarce distribution of GNSS sites in Venezuela, our model predicts 10.5 ± 0.5 mm yr⁻¹ of average right-lateral motion along the Boconó fault separating the Maracaibo block from stable SOAM. Several studies have proposed geological slip rate estimates along the Boconó fault for its different sections. For the northern section (the Yaracuy valley), Pousse-Beltran et al. (2017) use absolute dating of two offset alluvial fans to propose a Pleistocene right-lateral slip rates of 5–11.2 mm yr⁻¹. Along this section, our model predicts 11 ± 0.1 mm yr⁻¹. On the central section (Apartaderos at Mesa del Caballo), Audemard et al. (2008) estimate Holocene slip rates of 9–10 mm yr⁻¹ using palaeoseismological studies that agree with our estimate of 10.2 ± 0.3 mm yr⁻¹. On the southern section (between Mérida and San Cristobal), our model predicts 10.8 ± 0.1 mm yr⁻¹, which is slightly higher than the upper bound of the 4–8 mm yr⁻¹ estimated from palaeoseismological studies (Audemard 1997; Pousse-Beltran et al. 2017). Our model also finds sections (southern and central) of the Boconó fault partially locked down to 14 km depth (Coupling coefficient fixed at 60 per cent for our model: Table S2 and Text S1), suggesting a possible strain accumulation in secondary structures on either side of the main fault. In addition to strike-slip motion, our model predicts shortening across the Boconó fault ranging from 2.5 to 5 mm yr⁻¹ at the southern and central sections, increasing up to ∼7.0 mm yr⁻¹ at the northernmost section.

Fast and dominant right-lateral strike-slip motion occurs along the San Sebastian fault (17.2 mm yr⁻¹). Although no geological slip rate estimates are available for the San Sebastian fault, the local geodetic study from Pérez et al. (2018) proposed that the San Sebastian and the La Victoria faults (Fig. 2) accommodate ∼17 and ∼ 2 mm yr⁻¹ of the Caribbean–South America relative motion. Our model finds velocity residuals lower than 1 mm yr⁻¹ for all available GPS sites south of the San Sebastian fault (CRCS, USB0, USB1 and CAS2: Fig. S3), which are located on either side of the La Victoria fault (Fig. 3). Such a result prevents us from estimating the La Victoria fault slip rate, but our velocity residuals agree with the available geological slip rate (0.4–1.1 mm yr⁻¹) of Audemard et al. (2000). Outside the North Andean Sliver, sparse velocities do not allow solving shallow locking found for the El Pilar fault (Jouanne et al. 2011; Reinoza et al. 2015; Pousse-Beltran et al. 2016). However, our model is compatible with variable locking at their western and eastern fault sections. Regarding slip rates, we find right-lateral strike-slip motion at 18 mm yr⁻¹ with 1.7 mm yr⁻¹ of normal opening. Overall, our results support that both the San Sebastian and the El Pilar faults accommodate almost all of the Caribbean/SOAM relative motion.

In summary, we find good agreement between present-day fault slip rates and Pleistocene–Holocene slip rates for almost all faults along the eastern boundary of the NAS. Such an agreement suggests that a limited number of major faults accommodate the deformation, with little diffuse strain. At the large scale, our model finds that the eastern boundary of the North Andean Sliver is a dominant right lateral transpressive fault system with increasing slip rates from south to north. The relative magnitude of the strike-slip and shortening components changes with the changes of the strike of the fault with respect to the block motion (Fig. S4c).

7.2 Deformation within the Interandean Depression in Ecuador

In central Ecuador, the existence of the Quito-Latacunga (UIOL) block outlined by the foothills of the Western Cordillera, the Cosanga and the Chingual faults is required to explain decreasing geological slip rates observed from the Puna to the Pallatanga faults. Together, the surface traces of the El Angel fault and the westward dipping thrust faults across the Latacunga and Quito faults define the western boundary of the UIOL block (Lavenu et al. 1995; Alvarado et al. 2014, 2016). Seismic activity recorded by instrumental networks since the 1990s shows recurrent M_w ∼5 earthquakes (e.g. 1990 M_w 5.3, 2014 M_w 5.1) along the Latacunga and Quito faults (Ekström et al. 2012; Alvarado et al. 2014; Vaca et al. 2019), but greater earthquakes are documented along the whole western boundary of the UIOL block in the past centuries [e.g. the 1868 M_w ∼7.3 Ibarra, 1587 M_w ∼6.4 Guayallabamba and 1757 M_w ∼6.2 Latacunga earthquakes (Beauval et al. 2010, 2013)].

Regarding the Latacunga fault (LF), the relative motion between the NAW and UIOL blocks predicts shortening at 2.5 ± 0.1 mm yr⁻¹ at the latitude of Latacunga city (∼0.9°S), which can be compared to several geological estimates. Lavenu et al. (1995) proposed a shortening rate of 1.4 ± 0.3 mm yr⁻¹ constrained either by outcrops surrounding the Latacunga fault since the Late Plio-Quaternary or 2 ± 1 mm yr⁻¹ for the shortening integrated across the Interandean Depression based on the difference of Quaternary slip rates among the Pallatanga and Chingual faults (Ego et al. 1995). Lavenu et al. (1995) further indicated that these values could be considered as a lower bound. More recently, Baize et al. (2020) used InSAR data in the Interandean Depression and proposed horizontal shortening rates of 2–2.4 mm yr⁻¹. These values agree with the prediction from our model.

Although there is no geological slip rate estimate available for the Quito fault, detailed GNSS studies complemented by InSAR results proposed a range of 3–5 mm yr⁻¹ for east–west shortening rate (Alvarado et al. 2014; Mariniere et al. 2019). Mariniere et al. (2019) further found that very shallow locking is required to explain the sharp velocity gradient across the Quito fault observed in both GNSS and InSAR data. Although the GPS data used as input for our block model has a density around Quito scarcer than the dataset used by Mariniere et al. (2019), our model offers the advantage of integrating the GNSS measurements at a larger scale, enforcing the consistency between slip rates along the Latacunga and Quito faults. Our model predicts 3.4 ± 0.2 mm yr⁻¹ of E–W shortening at across the southern section of the Quito fault but with 1 ± 0.2 mm yr⁻¹ of right-lateral strike-slip motion, changing to 2.3 ± 0.2 and 2.2 ± 0.1 mm yr⁻¹, respectively on its northern section. Our model finds an average coupling coefficient not higher than 50 per cent for the Quito fault. Thus, our block model is in good agreement with the result of Mariniere et al. (2019), and provides a narrower range of relative motion across the fault (2.3–3.4 mm yr⁻¹).

Further north, 2.4 ± 0.2 mm yr⁻¹ of shortening are found across the El Angel fault with 2.9 ± 0.1 mm yr⁻¹ of right-lateral motion, but no geological data are available there for a comparison.

7.3 Southernmost Romeral fault system in Colombia

The Romeral fault system is a complex and continuous structure, inherited from sequences of oceanic accretions since the Cretaceous (Taboada et al. 2000; Cedriel et al. 2003; Suter et al. 2008). It consists of several parallel faults running along the western foothills of the Central Cordillera in Colombia. The southernmost extension of this fault is thought to be connected to the northernmost extension of the El Angel fault in Ecuador (Paris et al. 2000; Taboada et al. 2000; Suter et al. 2008). By considering this geometry in model C (Fig. 6c), we find unrealistic 6.5 ± 1.5 mm yr⁻¹ reverse slip rates between latitudes 4°N and ∼0.7°N. Such a model neither agrees with slip rate predictions for the northern section of the Romeral fault (reverse slip of 1 ± 0.5 mm yr⁻¹ between latitudes 4°N and 6°N) nor its southernmost extension at the El Angel fault (∼2.0 mm yr⁻¹ of reverse slip). In addition, model C also predicts 4–6 mm yr⁻¹ of shortening across the Algeciras-Sibundoy faults, which disagrees with a predominantly strike-slip motion inferred from geological mapping, landforms and aerial images (Velandia et al. 2005; Tibaldi et al. 2007). Thus, we conclude that the lack of geodetic observations together with the block geometry of model C is not able to constrain the motion in southern Colombia adequately.

An alternative block geometry is provided by our best model (Fig. 6d). Within this configuration, we find left-lateral motion at 2 ± 0.1 mm yr⁻¹ with a shortening component at 1 ± 0.6 mm yr⁻¹ across the Romeral fault between latitudes 4°N and 6°N (Figs 12 and S4). Long-term slip rates remain largely unknown for the Romeral fault. Only Paris et al. (2000) suggest Quaternary rates of ∼1 mm yr⁻¹ for the Paraiso and Piendamó sections (latitude ∼3°N–5°N) based on geomorphic analysis, but lacking precise dating. Therefore, the southern section of the Romeral fault possibly accommodates only a fraction of the North Andean West/Romeral (NAW/ROME) relative motion. The scarce geodetic data available in that part of Colombia prevents us from defining the boundary, which might consist of several faults. Understanding the deformation in southern Colombia is left for future studies.

7.4 Sub-Andean deformation in Ecuador and northern Peru

The sub-Andean domain is characterized by low-angle thrust ramps rising as fold and thrust belts of sedimentary rocks, east of the Eastern Cordillera in Ecuador and Peru (Suárez et al. 1983; Beck & Ruff 2005; Calderón et al. 2013; Eude et al. 2015; Alvarado et al. 2016; Villegas-Lanza et al. 2016b; Baby et al. 2018). In our model, the kinematics of this domain is described by the motion of an idealized and certainly oversimplified rigid Subandean block (EPSUB). Across the Eastern Subandean Belt (ESB), E–W shortening occurs at 1.5 ± 0.4 mm yr⁻¹ with 1.6 ± 0.2 mm yr⁻¹ of left-lateral motion in northern Ecuador surrounding the Napo region, east of the Chingual fault (effective slip rate = 1.8 ± 0.1 mm yr⁻¹, see Figs 12 and S4). Shortening across the ESB increases to 2.2 ± 0.2 mm yr⁻¹ with left-lateral strike-slip motion to 1.7 ± 0.1 mm yr⁻¹ (effective slip rate = 3 ± 0.2 mm yr⁻¹, Fig. 12) south of the Puyo and Pastaza regions (latitude 2°S–4°S). Although not statistically significant, this increase in slip rates (∼1 mm yr⁻¹) is consistent with a faster and more active deformation in the southern Ecuador sub-Andean domain compared to that observed at latitude 0°. Within the sub-Andean domain, the western boundary of the EPSUB is defined by the Macas fault in Southern Ecuador. Our best model estimates dominant E–W shortening across this structure at 1.1 ± 0.1 mm yr⁻¹, confirming that fault as the eastern limit of the Peruvian Forearc West block (PFW) in southern Ecuador.

In northern Peru, the Marañon fault accommodates 1 ± 0.2 mm yr⁻¹ of transpressive motion (reverse slip at 1 ± 0.1 with left-lateral strike-slip at 0.6 ± 0.2 mm yr⁻¹). In addition to the observed shortening component across the ESB, model D predicts 3.1 ± 0.1 mm yr⁻¹ of left-lateral shear surrounding the Yurimaguas region (latitudes 5°S–7°S, Figs 12 and S4). A left-lateral component of motion across the Eastern SubAndean Belt in northern Peru is kinematically consistent with a normal (opening) component of slip along the the Gulf of Guayaquil across the Puna and southern Pallatanga faults. It is also consistent with the partitioning of the NASCA/SOAM oblique convergence in Peru.

8 RELATIVE MOTION AT THE NAS BOUNDARIES

8.1 Nazca subduction

Our model combines the improved determination of the Nazca Plate kinematics presented in Jarrin et al. (2022) together with a finer description of the NAS fragmentation, both contributing to a more precise determination of the long-term slip rate at the Nazca subduction interface. Our model predicts convergence rates from 49.0 mm yr⁻¹ at latitude ∼4°S to 46.0 mm yr⁻¹ at latitude ∼1.5°N in Ecuador (average rate of 47.5 mm yr⁻¹). In southcentral Colombia, rates decrease from 45.5 mm yr⁻¹ at latitude ∼2°N to 43.1 mm yr⁻¹ at latitude ∼6.0°N (average rate of 44.3 mm yr⁻¹; Fig. 12).

The predicted Nazca/NAW (NAW: North Andean West, Fig. 6) trench-parallel motion increases from 4.2 mm yr⁻¹ at latitude ∼3.5°S to 6.3 mm yr⁻¹ at latitude ∼6°N, while the trench-normal motion decreases from 48.3 to 42.6 mm yr⁻¹ at the same latitudes. Compared to Nocquet et al. (2014), we find that their predicted velocities are 2–3 mm yr⁻¹ faster for the trench-normal motion and 0.5–1.5 mm yr⁻¹ faster for the trench-parallel motion. These discrepancies arise from: (1) the Euler pole used by Nocquet et al. (2014), for the Nazca Plate, is ∼1 mm yr⁻¹ faster than the pole used by our model (Jarrin et al. 2022) and (2) the NAS kinematics from Nocquet et al. (2014) predicts constant trench-parallel motion at 6 mm yr⁻¹.

In northern Peru, we estimate constant Nazca/Peruvian Forearc West (Nazca/PFW) convergence rates of ∼55 mm yr⁻¹ between latitudes 5°S and 8°S. These predicted convergence rates are consistent within ∼0.6 mm yr⁻¹ with the values reported by Villegas-Lanza et al. (2016b). Our additional data and analysis therefore confirms the view of an INCA sliver made of two separated blocks, namely the Peruvian Forearc West (PFW) and the Subandean (EPSUB), separated by the Marañon fault (Figs 2 and 12).

Our model further confirms that trench-parallel motion reverses at the latitude ∼5°S, with right-lateral strike-slip motion along the North Andean Sliver and left-lateral strike-slip along the Peruvian sliver, as would be expected from partitioning of the convergence vector for oblique subduction.

8.2 Interseismic coupling along the Ecuador–Colombia Subduction Zone

Fig. 13 shows our preferred ISC model along the Ecuador–Colombia subduction zone from latitude 3°S to 5°N. Alternative ISC models are available in the Supporting Information. We discuss the main patterns of the coupling distribution in terms of its relationship with historical and past seismic ruptures at the plate interface. The ISC map images three major segments (La Plata island, Galera-Tumaco and Buenaventura) of high coupling, separated by narrow zones of low or partial coupling.

The La Plata Island segment extends from the Salinas Peninsula (latitude ∼2°S) to the Manta Peninsula (latitude ∼1°S). It is characterized by a shallow (4–10 km) and fully locked patch located below La Plata Island. This isolated patch lies within a ∼1400 km long-creeping segment extending from central Ecuador to central Peru (latitude ∼10°S; Nocquet et al. 2014, 2016; Villegas-Lanza et al. 2016b). In the vicinity of La Plata island, three episodic shallow slow slip events were geodetically observed in the last decade (Vallée et al. 2013; Segovia et al. 2015; Rolandone et al. 2018). All of them suggest that the La Plata segment periodically releases the stress accumulated in a dominant aseismic mode (Rolandone et al. 2018), which is consistent with the absence of documented large historical seismic ruptures. Compared to the Nocquet et al. (2014, 2016), Chlieh et al. (2014), Gombert et al. (2018), Collot et al. (2017), Sagiya & Mora-Páez (2020) models, we found a second and deep patch of low coupling (∼20 per cent) between 50 and 80 km depth. In addition, along strike, a narrow strip (∼20 km long) of low coupling (10–20 per cent) is observed at the northernmost extension of the La Plata segment, extending along the down-dip direction from the trench down to 80 km. This transition zone likely acts as a barrier to seismic ruptures, as observed during the 2016 earthquake seismic rupture (Nocquet et al. 2016; Gombert et al. 2018).

From northern Ecuador to southern Colombia (from lat. 0.6°S to lat. 4°N), the second segment is fragmented into three sub-segments that extend along the great 1906 M_w ∼8.6 earthquake rupture area (Kanamori & McNally 1982; Ye et al. 2016). The first one, the Bahia-Pedernales subsegment, spreads along the 1942 and 2016 earthquake rupture areas (M_w 7.8). Here, this sub-segment is defined by a highly locked patch (ISC > 70 per cent) confined between 15 and 30 km depth. On this patch, a locked asperity (ISC: 80–90 per cent) overlaps most of the region of the maximum coseismic slip (>5 m) from the 2016 M_w 7.8 earthquake (Nocquet et al. 2016; Gombert et al. 2018). The second subsegment (Galera-Tumaco) extends from latitude ∼0.7°N to ∼2.5°N. This 200-km-long subsegment shows a highly coupled zone (70–90 per cent) between 10 and 40 km depth, in which three locked asperities (ISC > 80 per cent) are well identified (Fig. 13). Highly coupling areas on the Galera-Tumaco subsegment encompass major moment release areas from the 1958 M_w 7.7 and 1979 M_w 8.2 earthquakes (Beck & Ruff 1984; Sewnson & Beck 1996). The northernmost extension of the Galera-Tumaco subsegment seems bounded by a narrow strip (∼25 km long) of partial coupling (∼30 per cent) at the Sanquianga peninsula in southern Colombia that agrees with the end of the major moment release area from the 1979 earthquake (Beck & Ruff 1984; Sewnson & Beck 1996). Similarly, its southernmost extension is defined by a transition zone of partial coupling (ISC: 30–50 per cent). Therefore, this transition zone would act as a barrier separating two highly coupling subsegments (Bahia-Pedernales and Galera-Tumaco) that hosted large historical earthquakes.

The density of GPS sites decreases north of Tumaco, preventing us from imaging the possible continuity of locked ISC patches along central Colombia. The model shown in Fig. 13 uses a null coupling model as prior, meaning that ISC tends to be null when patches are not resolved by data. In central Colombia, it indicates a minimum coupling required by the GPS data at a few areas. For instance, the single GPS velocity at Guapi requires some intermediate (30–50 per cent) coupling to be explained, certainly less than in the Tumaco segment. Interseismic coupling at Buenaventura (lat. 4°N) is constrained by a profile including several GPS sites (Fig. 3). The model shows intermediate but significant coupling (50–60 per cent) from the trench down to 50 km depth. How this coupling extends along strike is not resolved. North of Buenaventura, no large historical earthquakes have been reported in the Bucaramanga segment, but this segment experienced two M_w 7.2 earthquakes in 1991 and 2004 (Fig. 1), leaving open the potential for large earthquakes in that area.

Overall, our ISC model is similar to previous models (Chlieh et al. 2014; Nocquet et al. 2016; Sagiya & Mora-Páez 2020; Chlieh et al. 2021). Aside from minor differences in the size of locked patches and the amount of shallow coupling close to the trench, our model finds moderate (20–40  per cent) 40–70 km deep coupling in Ecuador. We tested whether this deeper coupling is required by the data by running successive inversions limiting the subduction interface from 40 to 100 km depth (Fig S8) and we found the best fit for 80 km depth allowing partial deep coupling. Deep intermediate coupling is also consistent with the location of the 2015 deep slow slip event found by Rolandone et al. (2018).

8.3 Caribbean subduction

Tomography results and plate reconstructions support active subduction of the Caribbean Plate beneath the northern Andes since ∼75 Myr (Boschman et al. 2014; Spikings et al. 2015; Kellogg et al. 2019). Seismic reflection profiles image underthrusting of the Caribbean crust below a deformed accretionary prism in the northermost part of the NAS (Kellogg & Bonini 1982; Ladd et al. 1984; Hilst & Mann 1994; Kroehler et al. 2011; Boschman et al. 2014; Bernal-Olaya et al. 2015; Kellogg et al. 2019). Whether subduction is still active nowadays is crucial for both understanding the geodynamics of the Northern Andes and for seismic and tsunami hazard assessment. Previous GPS-based models support active subduction beneath northern Colombia but differ about the convergence rate and the part of the relative motion taken up by the deformation within the South Caribbean Deformed Belt (SCDB; Symithe et al. 2015; Pérez et al. 2018; Lizarazo et al. 2021). Symithe et al. (2015) modelled the northern part of the NAS as 3 blocks similar to our model C (Fig. 6c). Their model predicts a 9 mm yr⁻¹ convergence rate decreasing eastward to zero at the longitude 68°W and an increasing extension up to 5 mm yr⁻¹ at the junction with the El Pilar fault. Oppositely, Lizarazo et al. (2021) propose that the northern NAS is made of a single Macondo block with an ill-defined southern boundary. In this model, the convergence rate occurs at 7 mm yr⁻¹ and the velocity gradient observed at few available GPS data requires significant coupling down to 20 km depth.

For our model, the Caribbean/NAS relative motion in northwestern Colombia and Venezuela is similar to Symithe et al. (2015), showing a transition from active subduction (Fig. 14) in the western part of the SCDB to slow opening in the easternmost part (Fig. 12). In the details, our model shows a few differences compared to previous studies: the convergence in northwestern Colombia is accommodated partly along the subduction interface (∼3 mm yr⁻¹) but also by shortening across the San Jacinto Fold Belt (4–6 mm yr⁻¹). This pattern is required to improve the fit on GPS velocities (CART, BARU, BACO, URRA, SINC and CORO: Figs 14 and S3). It is further consistent with the view of a mechanically weak deforming accretionary prism. Our subduction rate (∼3 mm yr⁻¹) is therefore slower than Symithe et al. (2015) (9 mm yr⁻¹) or Lizarazo et al. (2021) (7 mm yr⁻¹). Oppositely to Lizarazo et al. (2021), we find no significant interseismic coupling at the Caribbean subduction interface, that was left as a free parameter in our model (Fig. 14). Relative motion is also found to be slower than previous estimates offshore Venezuela with a maximum convergence rate of 2 mm yr⁻¹ against 4 mm yr⁻¹ for Symithe et al. (2015) at longitude 72°W and 2 mm yr⁻¹ of oblique opening against 5 mm yr⁻¹ for Symithe et al. (2015) at the eastern end of the Bonaire block.

To summarize, the northernmost part of the NAS shows a complex deformation pattern that accommodates the 17 mm yr⁻¹ of the Caribbean Plate motion with respect to South America. Slip along the subduction interface accounts for at most 30 per cent of the relative motion in the western part to 0 per cent in the eastern part. Another 30 per cent is taken up along the Guaycaramo and Boconó faults. The remaining part occurs through a combination of shortening along the San Jacinto Fold Belt and strike-slip motion along the boundaries of the Maracaibo block.

8.4 Panama block collision

The on-going collision of the Panama block against the NAS produces an area of active and distributed deformation in northwestern Colombia (Pennington 1981; Audemard 1993, 2009; Adamek et al. 1988; Mann & Kolarsky 1995; Audemard & Audemard 2002; Trenkamp et al. 2002; Mora-Páez et al. 2019). Several shallow earthquakes are reported by the ISC-GEM catalogue [time window 1901–2019: Storchak et al. (2013, 2015), Di Giacomo et al. (2018)]. Available focal mechanisms are consistent with east–west shortening, also observed for the two largest earthquakes recorded during the instrumental period (1992 M_w 6.6 and 7.1 Mutata earthquakes, Fig. 15).

In our model D, the Panama/NAS boundary is defined from the location of clusters of shallow seismicity, sharp gradients found in GPS velocities and available neotectonic information. We name this block boundary the East Panama Deforming Zone (EPDZ; Figs 14 and 15a). Along the EPDZ, right-lateral strike-slip motion decreases from 15.5 mm yr⁻¹ at longitude 77°W down to 11 mm yr⁻¹ at longitude 76.4°W, whereas compression increases from ∼1 to 10 mm yr⁻¹ at the same longitudes. The Uramita fault [the easternmost Panama/NAS boundary: Duque-Caro (1990), Mann and Mann & Kolarsky (1995)] accommodates dominant compression at 5.6 ± 0.4 mm yr⁻¹ with 0.3 ± 0.2 mm yr⁻¹ of left-lateral motion (Figs 12, 15, 17 and S4). Together with the shortening at the San Jacinto fold belt described in the previous paragraph, the EPDZ defines a broad deforming zone accommodating the Panama and Caribbean eastward motion. That part of the NAS appears to be among the most active both in terms of velocity gradient with 1 cm yr⁻¹ accommodated over 100 km, slip rate as large as 1.5 cm yr⁻¹ and occurrence of large crustal earthquakes.

We also evaluate alternative boundaries for the EPDZ. We first test in Fig. 16(b) a boundary following the reverse Pierre and Unguia faults suggested in Page (1986), Mann & Corrigan (1990) and Mann & Kolarsky (1995). This model results in an increase of residuals by 3–5 mm yr⁻¹ for sites close to the collision area compared to the residuals from model D (Fig. 16a). We then test the proposition of a single Panama-Chocó block extending from Panama to the northern Colombia Pacific coast, as suggested from the long-term evolution of the Panama Arc (Duque-Caro 1990; Taboada et al. 2000; Montes et al. 2012; Barat et al. 2014; León et al. 2018; Kellogg et al. 2019). This block boundary follows the Uramita fault trace and its southward prolongation along the Garrapatas fault (Duque-Caro 1990; Paris et al. 2000; Taboada et al. 2000, Fig. 16c). In this case, we find a model unable to explain most of the observed velocities in the collision area, producing large systematic mismodelling with opposite residuals between Panama and the Colombia coast (Fig. 16c). As a final test, we consider two independent Panama and Chocó blocks separated by the Panama–Colombia border (Fig. 16d). Results indicate residuals up to 4 mm yr⁻¹ at sites within the Chocó block. Thus, GPS results do not support the existence of a Chocó block. Our exploration of several geometries for the eastern Panama supports that model D is the most suitable Panama–NAS boundary able to explain the observed GPS velocity field.

9 CONCLUSIONS

Our study provides the first kinematic model consistent at the scale of the northern Andes, simultaneously solving for the effects of rigid block rotations and elastic strain accumulation at the subduction interface and major crustal faults. Our preferred model includes 12 blocks that can explain more than 80 per cent of the observed interseismic velocity field. It is proposed that the North Andean Sliver is actually fragmented into 6 blocks.

The fast and oblique subduction of the Nazca Plate beneath South America controls the overall deformation regime. Trench-normal motion is mostly accommodated at the megathrust, but also by shortening taking place within the sub-Andean domain in northern Peru and Ecuador. The trench-parallel motion is accommodated by right-lateral slip producing the lateral translation-like motion of the sliver from 8.5 ± 0.1 mm yr⁻¹ towards N68° ± 1E (NAW block) to 17.5 ± 0.2 mm yr⁻¹ towards N84° ± 2.5E (BONA block). In addition, the Caribbean/South American convergence is accommodated offshore across the Southern Caribbean Deformed Belt and along onshore faults including the Oca-Ancon, Santa Marta-Bucaramanga faults and San Jacinto Fold Belt. That last fault system also accommodates some fraction of the Panama collision. Our results confirm that the Panama block is colliding against NW Colombia as a single rigid block based on the available GPS observations. However, future studies are required to evaluate possible motion variations in western Panama. The kinematics of the Inca sliver in northcentral Peru (as far south as latitude 14°S) is consistent with the motion of two rigid blocks accommodating deformation across known Quaternary fault systems.

In northern Venezuela, rapid slip rates take place along the San Sebastian-El Pilar faults. Both fault systems are moving at almost the same rate as the Caribbean Plate. However, a fraction of the Caribbean/SOAM relative motion is also transferred to the easternmost extensions of the Oca-Ancon and Boconó faults. The Nazca subduction and the Carnegie ridge entering the subduction in Ecuador trench often have been interpreted as the responsible mechanisms for the N–NE escape of the NAS since ∼2 Myr (Pennington 1981; Gutscher et al. 1999; Kellogg & Mohriak 2001; Trenkamp et al. 2002; Egbue & Kellogg 2010). However, the observed deformation pattern from GPS north of latitude ∼7°N suggests that the Panama collision and/or the Caribbean motion also contributes significantly to the ‘increase’ of eastward motion.

Within the NAS, two deformation patterns are observed: (1) localized deformation along faults showing slip rates of 2-4 mm yr⁻¹ (the Oca-Ancon, Santa Marta-Bucaramanga, Romeral, Latacunga, Quito, El Angel faults) and (2) distributed deformation broadening as far as hundreds of kilometres surrounding the Panama-northeastern Colombia suture zone (Uramita fault and Eastern Panama deformed zone). Outside the NAS, the Eastern Subandean Belt accommodates 2–3 mm yr⁻¹ of the Nazca/SOAM convergence. However, the convergence also appears to control a pattern of distributed deformation as far as ∼220 km east of the Andean Cordillera within the Ecuadorian Amazonia. The velocity field in that region shows velocities at ∼2 mm yr⁻¹, directed east–southeastward (Figs 3 and 17). Future coverage of continuous geodetic observations is required to better constrain the sub-Andean domain’s width and to delimit where its easternmost front reaches the stable part of South America in this region.

Slight motion along the Southern Caribbean Deformed belt at ∼4.5 mm yr⁻¹ offshore northern Colombia, decreasing eastward to 2 mm yr⁻¹ offshore northern Venezuela, is consistent with a low angle subduction interface without significant interseismic coupling. The Nazca Plate subducts beneath the NAS at ∼47.5 mm yr⁻¹ in Ecuador and ∼44.5 mm yr⁻¹ in southcentral Colombia. Our interseismic coupling model confirms a segmentation of the Nazca Plate interface consistent with the segmentation observed from the history of large historical and recent seismic ruptures. Detailed investigations beyond the scope of our study are required to investigate the link between the coupling distribution, seismicity and the numerous aseismic processes found in Ecuador.

As a main outcome, our model determines the slip rates at many already-identified active crustal faults. It provides new and precise inputs for future probabilistic seismic hazard assessment (PSHA) updates. The good agreement found between present-day slip rate estimates from our model and long-term geological estimates suggests that the modelled GPS velocity field accurately describes the NAS motion and deformation for the Quaternary period. Finally, our study points out a few areas where scarce geodetic observations were insufficient to propose a boundary, to estimate slip rates and level of interseismic coupling. Some of these regions have experienced major historical earthquakes in the past like the Panama/NAS collision zone or the coast of southern-central Colombia. Our study therefore helps to define where future GPS densifications are most useful.

SUPPORTING INFORMATION

suppl_data

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.

ACKNOWLEDGMENTS

This research was supported by the Secretaría Nacional de Educación Superior, Ciencia y Tecnología (SENESCYT-Ecuador) in the frame of P. Jarrin’s doctoral fellowship (grant number: IFTH-DFN-2018-0096/092-2017). We acknowledge support from the French National Research Institute for Sustainable Development (IRD) and Agence Nationale de la Recherche (ANR) in the frame of ANR-19-CE31-0003-01 S5 project and the ANR-15-CE04-0004 REMAKE project. We thank to UNAVCO, the Geophysical Institute of Ecuador (IG-EPN), the Military Geographical Institute of Ecuador (IGM), the Geophysical Institute of Peru (IGP), the Space Geodesy team (GeoRED) from the Colombian Geological Survey, the Agustin Codazzi Geographical Institute (IGAC) of Colombia, the Argentine Network for Continuous Satellite Monitoring and the Tommy Guardia Geographical Institute (IGNTG) of Panama for providing valuable geodetic data. We are grateful to Jack Loveless for his initial comments and suggestions in the block model. We thank B. Mercier de Lepinay for fruitful discussions about the Panama and Caribbean Plate. We thank the Editor and the two reviewers, Jack Loveless and Jeffrey Freymueller, for thoughtful reviews and comments that improved the paper. We finally thank Christophe Vigny and François Jouanne for their constructive comments on an earlier version of the manuscript.

DATA AVAILABILITY

Site information, GNSS velocities and block model parametrization are provided in the Supporting Information. Additional data related to this paper may be requested from the corresponding author. GNSS data are available online thought: (1) the GAGE Facility (UNAVCO), operated by EarthScope Consortium, with support from the National Science Foundation, the National Aeronautics and Space Administration and the U.S. Geological Survey under NSF Cooperative Agreement EAR-1724794, (2) Instituto Geográfrico Militar of Ecuador (IGM 2021), (3) Instituto Geográfrico Agustin Codazzi of Colombia (IGAC 2021), (4) Instituto Geográfrico Tommy Guardia of Panama (IGNTG 2021) and (5) Red Argentina de Monitoreo Satelital Continuo (Piñón et al. 2018). Earthquake catalogue and focal mechanism solutions were retrieved from the International Seismological Centre (Storchak et al. 2013; 2015; Di Giacomo et al. 2018), the Global Centroid Moment Tensor (Ekström et al. 2012) and the Servicio Geológico Colombiano (SGC 2021). The South American Risk Assessment (SARA) fault database is available in Alvarado et al. (2017).

References

APPENDIX

Eq. (A1) is used to quantify the goodness of fit from residual values in terms of χ². V_obs and V_mod are the observed and modelled velocities, respectively, σ is the standard error assigned to the observed velocities and n is the number of observations.

We use the variance reduction R to assess the accuracy increase of the model based on the block geometry evolution (Bougrine et al. 2019):

The F_ratio test (Stein & Gordon 1984; Nocquet et al. 2001) is used to compare how well two models derived from least squares estimation fit the data according to their degree of freedom. ${\chi }_{p_1}^2$ and ${\chi }_{p_2}^2$ are the chi-square statistics from the two models with p₁ and p₂ degrees of freedom, respectively

The empirical value of F is compared to the expected value from a Fischer–Snedecor distribution with p₁, p₂ degrees of freedom at a given confidence level [1 − P(per cent)].