# Detection of slow slip events along the southern Peru - northern Chile subduction zone

Jorge Jara, Romain Jolivet, Anne Socquet, Diana Comte, Edmundo

Norabuena

► To cite this version:

Jorge Jara, Romain Jolivet, Anne Socquet, Diana Comte, Edmundo Norabuena. Detection of slow slip events along the southern Peru - northern Chile subduction zone. Seismica, 2024, 3 (1), pp.17567-17579. 10.26443/seismica.v3i1.980 . hal-04979003

HAL Id: hal-04979003

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

Submitted on 20 Mar 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Detection of slow slip events along the southern Peru - northern Chile subduction zone

J. Jara ID 1, * R. Jolivet ID 1,2, A. Socquet ID 3, D. Comte ID 4,5, E. Norabuena ID 6

¹Laboratoire de Géologie, Départament de Géosciences, École Normale Supérieure, PSL Research University, CNRS UMR 8538, Paris, France,
²Institut Universitaire de France, Paris, France, ³Université Grenoble Alpes, Université Savoie Mont Blanc, CNRS, IRD, IFSTTAR, ITerre, Grenoble, France, ⁴Departamento de Geofísica, Facultad de Ciencias Físicas y Matemáticas, Universidad de Chile, Blanco Encalada 2002, Santiago, Chile, ⁵Advanced Mining Technology Center, Facultad de Ciencias Físicas y Matemáticas, Universidad de Chile, Av. Tupper 2007, Santiago, Chile, ⁶Instituto Geofísico del Perú, Lima, Perú

Author contributions: Conceptualization: A. S., J. J, R. J.. Methodology: J. J., R. J. Formal Analysis: J. J. Investigation: J. J. Resources: J. J., A. S., R. J., D. C., E. N. Writing - original draft: J.J. Writing - Review & Editing: J. J., R. J., A. S., D. C.

**Abstract** Detections of slow slip events (SSEs) are now common along most plate boundary fault systems globally. However, no such event has been described in the south Peru - north Chile subduction zone so far, except for the early preparatory phase of the 2014 Iquique earthquake. We use geodetic template matching on GNSS-derived time series of surface motion in Southern Peru - Northern Chile to extract SSEs hidden within geodetic noise. We detect 24 events with durations ranging from 17 to 36 days and magnitudes from $M_w$ 5.4 to 6.2. Our events, analyzed from a moment-duration scaling perspective, reveal values consistent with observations reported in other subduction zones. We compare the distribution of SSEs with the distribution of coupling along the megathrust derived using Bayesian inference on GNSS- and InSAR-derived interseismic velocities. From this comparison, we obtain that most SSEs occur in regions of intermediate coupling where the megathrust transitions from locked to creeping or where geometrical complexities of the interplate region have been proposed. We finally discuss the potential role of fluids as a triggering mechanism for SSEs in the area.

**Resumen** Hoy en día, las detecciones de eventos lentos (SSEs, por sus siglas en inglés) son comunes a lo largo de la mayoría de los sistemas de fallas activas a una escala global. Sin embargo, hasta ahora, no se han reportado eventos de este tipo en la zona de subducción del sur del Perú y norte de Chile (10°S-24°S), exceptuando aquellos ocurridos durante la fase de preparación del terremoto de Iquique de 2014. En el presente trabajo, nosotros utilizamos una técnica conocida como “Template Matching” en series temporales de desplazamiento medido por datos GNSS (Global Navigation Satellite System, GNSS por sus siglas en inglés) en el sur del Perú y el norte de Chile, para extraer la firma de eventos lentos asísmicos ocultos en el ruido geodésico. Nosotros detectamos 24 eventos asísmicos con duraciones de 17 a 36 días, y magnitudes de $M_w$ 5.4 a 6.2. El análisis de nuestros eventos utilizando leyes de escala momento-duración, revela valores consistentes con observaciones realizadas en otras zonas de subducción. El momento sísmico liberado por estos eventos es proporcional al cubo de su duración, lo que parece implicar una dinámica comparable con la de los terremotos clásicos. Los eventos detectados en este trabajo están principalmente localizados en zonas donde el acoplamiento intersísmico presenta valores en transición (0.3 - 0.8 de factor de acomplimiento), donde la zona de subducción transiciona de un estado bloqueado a uno de deslizamiento continuo. Finalmente, nosotros discutimos el rol potencial que podrían jugar los fluidos en el desencadenamiento de estos eventos lentos.

**Résumé** Depuis une vingtaine d'année, des événements de glissement asismiques ont été détectés le long de quasiment toutes les frontières de plaques au monde. Cependant, aucun n'a été décrit pour l'instant le long de la zone de Subduction allant du Perou au nord du Chili, si l'on omet le glissement mesuré lors de la période d'activité ayant mené au séisme d'Iquique en 2014. Nous utilisons une technique dite de Template matching sur des séries temporelles de déplacement mesuré par GNSS dans le nord du Chili pour extraire la signature d'événements de glissement asismiques cachés au sein du bruit géodésique. Nous détectons 24 événements asismiques avec des durées allant de 17 à 36 jours pour des magnitudes équivalentes allant de Mw 5.4 à 6.2. Nos événements ont des valeurs cohérentes avec les observations rapportées dans d'autres zones de subduction. Il apparaît que ces événements asismiques sont essentiellement localisés dans des zones de couplage intermédiaires où le megathrust est a mi-chemin entre un état bloqué et un état en glissement permanent. Nous discutons finalement de l'influence éventuelle de fluides profonds dans le déclenchement de ces événements asismiques.

Production Editor:

Gareth Funning

Handling Editor:

Wenbin Xu

Copy & Layout Editor:

Ethan Williams

Signed reviewer(s):

Piero Poli, Jeffrey

Freymueller

Received:

11 May 2023

Accepted:

3 May 2024

Published:

7 June 2024
# 1 Introduction

Overwhelming evidence suggest that the Elastic Re-bound Theory proposed by Reid (1910) after the 1906 California earthquake associated with the stick-slip behavior of frictional interface (Brace and Byerlee, 1966) is insufficient to explain the slip behavior along active faults. Geodetic measurements of surface motion have revealed the presence of aseismic, slow slip along all types of active faults. After the first descriptions in the mid-20th century from direct observations of damage to human-made structures crossing the San Andreas (Louderback, 1942; Steinbrugge et al., 1960) and North Anatolian (Ambraseys, 1970) faults, aseismic slip has been directly observed, or inferred, from geodetic measurements at different stages of the earthquake cycle. For instance, afterslip corresponds to the diffusion of slow slip during the post-seismic period accommodating a co-seismic stress perturbation (e.g., Heki et al., 1997; Bürgmann et al., 2001; Hsu et al., 2002, 2006). Creep, on the other hand, often refers to steady aseismic slip during the interseismic period (Steinbrugge et al., 1960; Ambraseys, 1970; Jolivet et al., 2015b). In addition, interseismic transients (i.e., slow slip events or SSEs) during this interseismic period were discovered in the 2000s along subduction zones. SSEs often locate in the deeper portion of the seismogenic zone (e.g., Hirose et al., 1999; Dragert et al., 2001), but some of these SSEs are associated with seismic signals that occur within the seismogenic zone, and may contribute to reducing geodetic coupling (Mazzotti et al., 2000; Bürgmann et al., 2005; Loveless and Meade, 2010; Radiguet et al., 2012; Béjar-Pizarro et al., 2013; Villegas-Lanza et al., 2016; Métois et al., 2016; Michel et al., 2019a; Jolivet et al., 2020; van Rijsingen et al., 2021; Lovery et al., 2024). This along-dip segmentation differs from one subduction zone to the other (Nishikawa et al., 2019) and we note more occurrences of SSEs along young, warm subduction zones (i.e., Nankai, Mexico, Cascadia), than old and cold ones. Finally, slow slip appears to be an important ingredient of the preparation phase of earthquakes (e.g., Ruegg et al., 2001; Ruiz et al., 2014; Radiguet et al., 2016; Socquet et al., 2017; Voss et al., 2018). More recently, it has been proposed that a significant fraction of observed geodetic displacement in seismically active regions results from the occurrence of slow slip events (Jolivet and Frank, 2020, and reference therein), suggesting a burst-like, episodic behavior of aseismic slip at all time scales from seconds to decades in places as varied as Mexico (Frank, 2016; Rousset et al., 2017; Frank and Brodsky, 2019), Cascadia (Michel et al., 2019a; Ducellier et al., 2022; Itoh et al., 2022), along the San Andreas Fault (Khoshmanesh and Shirzaei, 2018; Rousset et al., 2019; Michel et al., 2022), the Haiyuan fault in Tibet (Jolivet et al., 2015a; Li et al., 2021), on the Alto Tiberina and Pollino fault systems in Italy (Gualandi et al., 2017; Cheloni et al., 2017; Essing and Poli, 2022), or Japan (Nishimura et al., 2013; Takagi et al., 2019; Nishikawa et al., 2019; Uchida et al., 2020). All observations suggest the importance of accounting for

**Non-technical summary** Earthquakes correspond to a sudden release of elastic energy stored in the crust as a response to the relative motion of tectonic plates. However, this release of energy is not always sudden and accompanied by destructive seismic waves. It sometimes happens slowly during aseismic, slow slip events. It has been shown that SSEs can be associated with the nucleation, propagation, and termination of big earthquakes. SSEs have been detected along many subduction zones in the world but not in northern Chile, yet. Here, we use a template matching method to scan GNSS observations of ground motion to detect and characterize slow slip events along the southern Peru - northern Chile subduction zone. We find 24 aseismic events at depths comparable with that of SSEs in other subduction zones, as well as in regions that slip aseismically persistently. We discuss how our findings relate to past earthquake ruptures, the geometry of the subduction zone, and fluids circulating at depth. Our results show the importance of implementing methods to extract small aseismic signals in noisy data, key observations for a better understanding of fault mechanics.

aseismic slip in our understanding of earthquake cycle dynamics. However, the underlying physics controlling aseismic slip is still debated, mainly due to the lack of good, dense observational databases.

Nowadays, observations of aseismic slip in subduction zones are frequently documented over a wide range of slip amplitudes and at different stages of the earthquake cycle (Avouac, 2015; Obara and Kato, 2016; Bürgmann, 2018; Kato and Ben-Zion, 2021, and references therein). Regular slow slip events have been documented mainly along warm subduction zones such as Cascadia, Nankai (southwest Japan), Mexico, or New Zealand (e.g., Graham et al., 2016; Nishikawa et al., 2019; Wallace, 2020; Michel et al., 2022, and references therein). Instead, observations of slow slip events in cold subduction zones such as off-shore Japan or Chile are sparse or indirect, through seismic swarms, repeaters, or slow earthquakes (Kato et al., 2012; Kato and Nakagawa, 2014; Gardonio et al., 2018; Nishikawa et al., 2019), and rarely with geodetic observations (Hino et al., 2014; Ruiz et al., 2014; Socquet et al., 2017; Boudin et al., 2021). Geodetic displacement corresponding to such slow slip events are usually of mm to cm-scale amplitude and require the development of novel and systematic methods to extract SSEs from noisy time series of geodetic data (Frank, 2016; Rousset et al., 2017; Michel et al., 2019a; Uchida et al., 2020; Itoh et al., 2022).

We focus on the South Peru-North Chile subduction zone. The region is seismically active, with two historical earthquakes in 1868 (southern Peru), and 1877 (northern Chile), both tsunamigenic earthquakes of magnitude ~8.5 (Kausel, 1986; Comte and Pardo, 1991; Vigny and Klein, 2022) (Figure 1). Since these two events, the region has experienced several large earthquakes ($M_w > 7.5$) (Ruiz and Madariaga, 2018) accompanied by an important background seismic activity (Jara et al., 2017; Sippl et al., 2018, 2023) (Figure 1). In addition, coupling is highly variable along the subduction interface. Coupled regions overlap with the inferred rupture extent of the 2001 $M_w 8.1$ Arequipa and 2014 $M_w 8.1$ Iquique earthquakes (Schurr et al., 2014; Métois et al., 2016; Villegas-Lanza et al., 2016; Jolivet

<footnote>*Corresponding author: now at GFZ Potsdam, jorge@gfz-potsdam.de</footnote>
<figure><img src="image_2.png" /></figure>

**Figure 1** Seismotectonic map of the South Peru - North Chile subduction zone. White arrows show the extent of historical earthquakes (Comte and Pardo, 1991; Vigny and Klein, 2022). Gray contours are the rupture area of instrumental earthquakes with M>7.5, with corresponding epicenters (gray starts) and focal mechanisms (if available) (Dorbath et al., 1990; Beck and Ruff, 1989; Hartzell and Langer, 1993; Delouis et al., 1997; Chlieh et al., 2004; Pritchard et al., 2007; Dziewonski et al., 1981; Ekström et al., 2012; Peyrat and Favreau, 2010; Sladen et al., 2010; Béjar-Pizarro et al., 2010; Duputel et al., 2015; Jara et al., 2018). Yellow lines are the 0.1 m afterslip contours available in the region (Chlieh et al., 2004; Béjar-Pizarro et al., 2010; Remy et al., 2016; Hoffmann et al., 2018), whereas the green ones are the pre-seismic slip reported for Iquique earthquake by Socquet et al. (2017). Colored dots are earthquakes with M>4.0 from the International Seismological Centre (International Seismological Centre, 2016) over the period 1990 - 2016, color-coded by depth and scaled by magnitude. Large white arrow shows convergence direction and rate from Métois et al. (2016). SOAM: South America plate.

et al., 2020). A large coupled section is inferred where the 1877 earthquake is thought to have ruptured (Jolivet et al., 2020; Vigny and Klein, 2022). In addition, two low-coupling regions are observed. In southern Peru, low coupling coincides with the subduction of the Nazca ridge (~15°) (Villegas-Lanza et al., 2016; Lovery et al., 2024). In northern Chile, a reduction in coupling is inferred offshore Iquique and below the Mejillones peninsula (~21°) (Béjar-Pizarro et al., 2013; Métois et al., 2016; Jolivet et al., 2020).

In addition to low coupling, aseismic slip has been observed in South Peru and North Chile. Afterslip has been reported following large earthquakes, including the 1995 $M_w$ 8.1 Antofagasta (Chlieh et al., 2004; Pritchard and Simons, 2006), the 2001 $M_w$ 8.1 Arequipa (Ruegg et al., 2001; Melbourne, 2002), the 2007 $M_w$ 8.0 Pisco (Perfettini et al., 2010; Remy et al., 2016), the 2007 $M_w$ 7.7 Tocopilla (Béjar-Pizarro et al., 2010) and the 2014 $M_w$ 8.1 Iquique earthquakes (Hoffmann et al., 2018)

(Figure 1). Geodetic transients interpreted as the signature of aseismic slip occurred in the days to months preceding the $M_w$ 8.4 Arequipa earthquake in 2001, before one of its largest aftershock, and preceding the Iquique earthquake in 2014 (e.g., Ruegg et al., 2001; Melbourne, 2002; Ruiz et al., 2014; Schurr et al., 2014; Socquet et al., 2017). Aseismic slip is considered responsible for a significant fraction of such geodetic transients (Twardzik et al., 2022). There is therefore plenty of evidence of occurrences of aseismic slip in this broad region but, despite intense efforts to instrument the area, no obvious spontaneous slow slip events have been detected during the interseismic period.

A change in the interseismic surface velocity field was observed following the $M_w$ 7.5 intermediate-depth Tarapaca earthquake over a decade (Peyrat et al., 2006; Peyrat and Favreau, 2010) (Figure 1), an observation interpreted as the signature of a decoupling of the subduction interface (Ruiz et al., 2014; Jara et al., 2017). Com-
parable changes in surface velocity field, observed following the 2010 Maule earthquake, have also been observed in the regions affected by the 2015 Illapel (Ruiz et al., 2016) and 2016 Chiloé (Ruiz et al., 2017; Melnick et al., 2017) earthquakes. Such shifts in surface velocity may be linked to postseismic viscoelastic processes acting over long distances (Bouchon et al., 2018) in contrast to the localized behavior observed after the Tarapaca earthquake (Jara et al., 2017). Over the same period, we observed a significant increase in background seismicity (Jara et al., 2017), as well as an apparent synchronization of intermediate-depth and shallow seismic activities (Bouchon et al., 2016; Jara et al., 2017). Changes in background seismicity rates have been associated with the occurrence of aseismic slip events and fluid migration (Marsan et al., 2013; Reverso et al., 2016; Marsan et al., 2017). The synchronization of the seismicity is interpreted as related to aseismic slip events occurring along the subduction interface due to a broader slab deformation (Bouchon et al., 2016). These indirect observations suggest aseismic transients may occur in South Peru - North Chile during the interseismic period.

We aim to detect small, short-term aseismic slip events in this region and discuss their occurrence and location with respect to the interseismic coupling pattern and past seismic crises. We explore GNSS time series, searching for small transients, using a geodetic template matching approach (Rousset et al., 2017). We use GNSS and InSAR data to infer an updated distribution of interseismic coupling using a Bayesian framework following the approach of Jolivet et al. (2020), comparing the detected aseismic events with the coupling model, along with geophysical information available in the region (seismicity, Vp/Vs ratio, gravity models). We finally discuss potential mechanisms explaining the occurrence of aseismic events in the area.

# 2 Data, Methods and Results

## 2.1 GNSS processing and time series analysis

We process data from 119 continuous GNSS (cGNSS) sites in the central Andes region (Figure S1a) and worldwide (Figure S1b), using a double difference approach with the GAMIT/GLOBK software (Herring et al., 2015). 67 cGNSS sites are in the South Peru - North Chile region (Figure S1a and Figure 2, brown arrows), installed and maintained by the Integrate Plate boundary Observatory Chile (IPOC) (Klotz et al., 2017), the Laboratoire International Associé "Montessus de Ballore" (LIA-MB) (Klein et al., 2022), the Central Andean Tectonic Observatory (CAnTO, Caltech) (Simons et al., 2010), the Instituto Geofísico del Perú (IPG) (Jara et al., 2017; Socquet et al., 2017), the Institut des Sciences de la Terre (ISTerre) (Jara et al., 2017; Socquet et al., 2017), and the Centro Sismológico Nacional of Chile (CSN) (Báez et al., 2018). The remaining 52 stations are part of the International GNSS Service (IGS) (Teunissen and Montenbruck, 2017) global network. We separate these stations into three subnetworks (two locals and one global) with 33 overlapping stations, where the local separation depends on the station data span: one local network

with data from 2000-2014 and the other including data from 2007-2014. Global network processing includes 99 stations over the 2000 - 2014 period, with 22 stations in South America (Figure S1b). We use the GAMIT 10.6 software (Herring et al., 2015), choosing the ionosphere-free combinations and fixing the ambiguities to integer values. We use precise orbits from the IGS, precise earth-orientation parameters (EOPs) from the International Earth Rotation and Reference System Service (IERS) bulletin B, IGS tables to describe the phase centers of the antennas, FES2004 ocean-tidal loading corrections, and atmospheric loading corrections (tidal and non-tidal). We estimate one tropospheric zenith delay every two hours and one pair of horizontal tropospheric gradients per 24h session using the Vienna Mapping Function (VMF1) (Boehm et al., 2006). We use the GLOBK software to combine daily solutions and the PY-ACS software (Nocquet, 2018) to derive position time series in the ITRF 2008 reference frame (Altamimi et al., 2011). Finally, time series are referenced to fixed South America considering the Euler pole solution proposed by Nocquet et al. (2014).

We fit the time series with a parametric function of time for each component (N, E, and U) (Bevis and Brown, 2014). Each time series $x(t)$, function of time $t$, is modeled as

$$x(t) = x_R + v(t - t_R) + \sum_{j=1}^{n_j} b_j H(t - t_j) \\ + \sum_{k=1}^{n_F} [s_k \sin(\omega_k t) + c_k \cos(\omega_k t)] \\ + \sum_{i=i}^{n_T} a_i \log(1 + t_i / \Delta T), \quad (1)$$

where $x_R$ is a reference position at a time $t_R$ and $v$ is the interseismic velocity for each component. $H$ is a Heaviside function applied each time $t_j$ an earthquake (or antenna change) offsets the time series. The combination of $\sin$ and $\cos$ functions describes seasonal oscillations (with annual and semi-annual periods), while the logarithmic function models the transient, post-seismic signal following large earthquakes ($M_w \ge 7.5$) with a relaxation time $\Delta T$. For a given station, we consider a Heaviside function for all earthquakes of magnitude larger than 6 with an epicenter to station distance lower than $d(M) = 10^{\frac{M}{2}-0.8}$, as proposed by the Nevada Geodetic Laboratory (www.geodesy.unr.edu). We only include a post-seismic term for earthquakes of magnitude larger than 7.5. All inferred parameters for each component and each cGPS site are in Supplementary Information, Tables S1-S38. Figures S2-S17 compare the data and model at each station. We then estimate and remove a common-mode error by stacking all the time series (Bock and Melgar, 2016; Socquet et al., 2017; Jara et al., 2017). This procedure enables us to get residual time series (Figures S18-S19) as well as an interseismic velocity field (Table S1-S2). We use the obtained residual time series to search for geodetic transients compatible with slip on the megathrust and use the geodetic velocity field to update the last published coupling map (Jolivet et al., 2020).
<figure><img src="image_3.png" /><figcaption>**Figure 2** Geodetic data. (a) Colored dark green and pink arrows are the GNSS interseismic velocities from Métois et al. (2016) and Villegas-Lanza et al. (2016), respectively, while brown arrows are the continuous GNSS processed in this study. The inset shows the residual trench perpendicular displacement time series for GNSS station UAPE. (b) Line-of-sight (LOS) interseismic ground velocity from track 96 (Envisat data) from (Jolivet and Simons, 2018; Jolivet et al., 2020). Black arrows indicate the flight direction of the satellite and its line of sight (LOS).</figcaption></figure>

## 2.2 Fault Geometry and Green's Functions

Coupling map estimation and geodetic template matching methods need a fault geometry and Green's functions calculation, as described below. In both cases, we define the geometry of the megathrust using Slab 2.0 (Hayes et al., 2018) as a reference, but with different meshing strategies. For the coupling case, we use triangles with 10 km-long sides along the coast and 25 km-long sides, both at the trench and depth, between latitudes 17°S-25°S. In the northern part (10°S-17°S), we adapt the size to the GNSS station density, considering a constant 50 km-long triangle side. In contrast, in the geodetic template matching case, we use triangles with 10 km-long sides along the coast and 25 km-long sides in the entire region. Then, we consider slip on the fault as the linear interpolation of slip values at the mesh nodes. Finally, we compute the Green's functions assuming a stratified elastic medium derived from Husen et al. (1999) using the EDKS software (Zhu and Rivera, 2002).

## 2.3 Coupling map for Southern Peru - Northern Chile

We update the distribution of coupling from Jolivet et al. (2020) in order to compare short- (i.e., days to months) and long-term (i.e., years to decades) aseismic deformation in the region. We use the GNSS velocity fields

from Métois et al. (2016) (data span 1996 - 2013) and Villegas-Lanza et al. (2016) (data span 2008-2013), that we complement with our GNSS velocity field (Figure 2a, data span 2000-2016). Additionally, we use the line of sight (LOS) velocity map from Jolivet et al. (2020), derived from the processing of Envisat data covering the period 2003 - 2010 (Figure 2b).

We use the backslip approach to estimate the distribution of coupling (Savage, 1983). A coupling of 1 (resp. 0) corresponds to a fully locked megathrust (resp. a megathrust that slips at plate rate). We consider plate motion estimated by UNAVCO (www.unavco.org) under the ITRF 2014 model (Altamimi et al., 2016) to estimate the convergence rate, angle, and rake on each node of the fault mesh. The backslip rate is evaluated by subtracting the sliver movement proposed by Métois et al. (2016) in Chile (11 mm/yr) and by (Villegas-Lanza et al., 2016) in Peru (5.5 mm/yr) to the convergence rate. In the Arica bend (16°S - 18°S), at the boundary of the Chilean and Peruvian slivers, we build a gradient to make a smooth transition between the two slivers. We solve for the distribution of models that satisfy the geodetic data.

The forward problem is written as $d = Gm$, with $d$ the geodetic data (GNSS and InSAR velocities), $m$ the vector of parameters to solve for and $G$ the Green's functions (Section 2.2). Parameters include coupling at each mesh node and geometric transformations akin to those in Jolivet et al. (2020). We adopt a probabilistic approach to
estimate the parameters in order to evaluate the associated uncertainties. The *a posteriori* Probability Density Function (PDF) of a model **m** given a dataset **d**, p(**m**|**d**), writes as

$$p(\mathbf{m}|\mathbf{d}) \propto p(\mathbf{m})p(\mathbf{d}|\mathbf{m}), \qquad (2)$$

where p(**m**) is the *a priori* model PDF and p(**d**|**m**) is the data likelihood. The *a priori* PDF describes our knowledge of coupling along the megathrust before collecting geodetic data. We define the *a priori* PDF at each node for the coupling factor as follows:

$$\mathbf{X} \sim \begin{cases} \mathcal{N}(\mu_c, \sigma_c^2) & \text{if } -0.1 \le \mathbf{X} \le 1.1 \\ 0 & \text{otherwise} \end{cases} \qquad (3)$$

where $\mu_c$ and $\sigma_c$ are the mean and standard deviation of a normal distribution. We select the bounds of [-0.1, 1.1] to ensure an accurate sampling for the full range of coupling values between 0 and 1 (Dal Zilio et al., 2020a; Jolivet et al., 2020). We know the megathrust is decoupled below 60 km depth from geodetic (Chlieh et al., 2004; Béjar-Pizarro et al., 2013; Jolivet et al., 2020) and seismological evidence (Comte et al., 2016). Thus, we apply an *a priori* condition based on the depth of each node. If a node is deeper than 60 km, the *a priori* mean ($\mu_c$) is set to 0 and the standard deviation ($\sigma_c$) to 0.1. In cases where a node is shallower than 60 km, we assign an *a priori* mean ($\mu_c$) of 0.5 and a standard deviation ($\sigma_c$) of 0.5.

We adopt a Gaussian formulation for the data likelihood, p(**d**|**m**), which writes as

$$p(\mathbf{d}|\mathbf{m}) = \frac{1}{\sqrt{2C_{\chi}}} \exp\left\{\left[-\frac{1}{2}(\mathbf{G}\mathbf{m} - \mathbf{d})^T C_{\chi}^{-1} (\mathbf{G}\mathbf{m} - \mathbf{d})\right]^2\right\},$$

where $C_{\chi}$ is the misfit covariance matrix (Duputel et al., 2014) defined as $C_{\chi} = C_p + C_d$, where $C_d$ is the data covariance matrix (data uncertainties), while $C_p$ is the prediction error covariance matrix, representing uncertainties on the assumed elastic model (P and S wave velocities and density). We assume a 10% error on the elastic parameters following Jolivet et al. (2020).

We explore the model space using Altar (altar.readthedocs.io) to sample the *a posteriori* PDF of the coupling factor, generating 250000 models. Altar is based on the Cascading Adaptive Transitional Metropolis in Parallel (CATMIP) algorithm (Minson et al., 2013; Duputel et al., 2014; Jolivet et al., 2015b). These models enable us to perform statistics, derive the mean model for the interseismic coupling (Figure 3), and collect information about the model resolution (see Supporting Information for model GNSS and InSAR residuals, Figure S20-S23, as well as Standard Deviation, Mode, Skewness, and Kurtosis, Figure S24).

The mean coupling model (Figure 3a), is close to previously published models in the region (e.g., Chlieh et al., 2011; Béjar-Pizarro et al., 2013; Métois et al., 2016; Villegas-Lanza et al., 2016; Jolivet et al., 2020; Lovery et al., 2024), especially considering the along-strike segmentation. Our model differs from previously published models in the coupling intensity at locked patches, as well as the depth of these coupled patches.

In Peru, we observe three patches with interseismic coupling that varies between 0.5-0.75 (Figure 3a). Previous models report similar patches, although totally locked (coupling factor ~ 1) (Chlieh et al., 2011; Villegas-Lanza et al., 2016; Lovery et al., 2024). Unfortunately, the density of GNSS stations in this region is not anywhere near that in Chile, hence the large standard deviations in the Peruvian region (Figure S25). Analyzing the moments of the *a posteriori* PDF, including standard deviation, skewness and kurtosis confirms this (Figure S24). Similarly, these moments show that the resolution at the trench over the entire region is low. Additionally, our model varies from those constrained only by GPS data in Chile (e.g., Métois et al., 2016). The InSAR data helps constraining interseismic coupling at depth (Béjar-Pizarro et al., 2013; Jolivet et al., 2020) and the strong *a priori* coupling damps potential large variations at depth, which we consider not physical.

## 2.4 Detection of aseismic slip events with geodetic template matching

### 2.4.1 Methodology

We use a geodetic template matching approach to detect potential aseismic slip events on the residual GNSS time series (Section 2.1). We summarize here the method presented in detail by Rousset et al. (2019). We search for the spatio-temporal signature of slip events in surface displacement time series by cross-correlating synthetic templates with our GNSS residual time series, in velocity. These templates correspond to the surface displacement caused by slip on dislocations located on the subduction megathrust embedded in a stratified, semi-infinite elastic medium. We calculate such templates (w) by convolving the Green's functions (Section 2.2) with a time-dependent slip evolution s(t) defined as

$$\mathbf{s}(t) = \frac{1}{2} \left[ 1 - \cos \left( \frac{\pi t}{T} \right) \right], \qquad (5)$$

where T is the duration of a synthetic event. Following Rousset et al. (2019), we derive for each template the weighted correlation function for each fault node, defined as

$$\mathbf{C}_f(t) = \frac{\sum_{i=1}^{2N} |\mathbf{G}_i| \mathbf{C}_i(t)}{\sum_{i=1}^{2N} |\mathbf{G}_i|}, \qquad (6)$$

where **G** is the Green's functions and $C_i$ is the correlation between the time series and the synthetic template at a given fault node i given by

$$\mathbf{C}_i(t) = \frac{\sum_{k=1}^{T} \dot{\mathbf{w}}_i(t_k) \dot{\mathbf{d}}_i(t_k + \tau)}{\sqrt{\sum_{k=1}^{T} \dot{\mathbf{w}}_i^2(t_k) \sum_{k=1}^{T} \dot{\mathbf{d}}_i^2(t_k + \tau)}}, \qquad (7)$$

where $\dot{\mathbf{w}}$ and $\dot{\mathbf{d}}$ are the time derivatives of the template in terms of displacement (i.e., the template's velocity with duration T), and the time derivatives of the GNSS time
<figure><img src="image_4.png" /></figure>

**Figure 3** Location of detected aseismic slip events. Markers are color-coded by time of occurrence and scaled by magnitude. Four examples of weighted stacked correlations are shown with the event id number. Red line is the best fit model used to evaluate the event magnitude and duration, considering their estimated σ. Background color from white to dark through yellow and red is the mean coupling distribution. Black red areas (coupling factor ~1) are locked regions, while transparent areas (coupling factor ~0) are regions that slip aseismically at a rate equal to the plate convergence rate. Gray contours show instrumental ruptures. Yellow contours are afterslip regions, whereas green ones indicate slip inferred during the period preceding the Iquique earthquake. White arrows are the historical rupture extents.

series, respectively. τ denotes a moving time variable that enables the temporal matching search between templates and observations. We then search for peaks in C_f(t) corresponding to candidate slip events. As can be seen in the Supporting Information (see Fig. S31b, red and black lines), in the case of synthetic events, the correlation peaks in C_f arise from the geodetic noise using as many GNSS stations as possible.

For each candidate slip event, we stack the time series of displacement weighted by Green's functions around the time of detection (see Supporting Information Figure S31b, for an example of stacks on synthetic time series, purple and yellow lines). Such weighting accounts for displacement amplitude and direction, increasing the signal-to-noise ratio (Rousset et al., 2017). Stacks are computed over a period of 180 days, centered on each potential occurrence. On each stack, we estimate two linear trends, before and after the candidate occur-

rence, and the time dependent slip evolution of Eq. 5 to the weighted stack in order to determine the amplitude, the start and end date of each detected transient. We apply a non-linear regression to determine the posterior Probability Density Function of the model parameters given a stack of time series following Tarantola (2005). Effectively, we use an MCMC algorithm to derive 30,000 samples from the posterior PDF and evaluate the mean and standard deviation of the duration and magnitude of each candidate slow slip event.

In order to curate the potential detections from artefacts, we perform a sensitivity and resolution analysis, to determine the minimum magnitude of a slip event that can be detected for each fault node. Although the method above has been extensively described by Rousset et al. (2019), the novelty of our approach relies on the evaluation of uncertainties through a Bayesian exploration of all important parameters.
<figure><img src="image_5.png" /><figcaption>**Figure 4** Event magnitude as a function of the resolution magnitude of the node where the event is located. Red crosses are events that passed the resolution test. Dashed blue line is the 1:1 line that separates validated from excluded events.</figcaption></figure>

## 2.4.2 GNSS network sensitivity and resolution

We analyze the sensitivity of our approach by testing its ability to detect, locate, and estimate the source parameters (magnitude and duration) of synthetic aseismic slip events. We first evaluate the parameters characterizing the noise affecting each GNSS time series of displacement by building synthetic time series of noise on which we perform the tests. In order to generate synthetic noise, we model each component of the residual time series (Eq. 1) as a combination of white and colored noise (Williams, 2003), such as,

$$ \mathbf{P}(f) = P_0 (\mathbf{f}^{-\alpha} + f_0^{-\alpha}), \quad (8) $$

where P is the power spectrum as function of temporal frequency f, P₀ and f₀ are normalization constants, and α is the spectral index. We explore P₀, f₀, and α using Bayesian inference to estimate their mean and standard deviation at each station component (see the Supporting Information for further details and an example of the power spectrum and the probability density function (PDF) of parameters at the UAPE station in Figures S26 - S27, as well as Tables S39 - S42 for all the network noise parameters inferred). We use these inferred noise parameters to build 1000 synthetic time series of displacement at each GNSS station. We use these synthetic time series to estimate thresholds of detection for each fault node.

The number of GNSS stations in the study area has evolved during the observation period. We, therefore, must consider three periods independently depending on the number of active stations: 2000 - 2003 (four stations), 2004 - 2007 (20 stations), and 2008 - 2014 (55 stations). We first determine which stations are able to cap-

ture a slow slip event on a given node. For each period and fault node, we correlate the 1000 synthetic time series of noise with a template of a duration of 40 days and slip equivalent to a magnitude Mw 6.0. We evaluate the standard deviation of the resulting weighted correlation functions, σt, as a minimal threshold to be exceeded (i.e., when dealing with time series that might include slip events, a peak of correlation higher than 3σt is a positive detection).

Once this threshold has been defined, we compute the weighted correlation function for 1000 time series of noise to which we have added the signal of synthetic transients with different duration (10, 20, and 30 days) and magnitudes (5.0 - 7.0 Mw, every 0.1 of magnitude). In case of a detection, we stack the displacement time series around the detection time. We consider a synthetic event has been correctly detected and located if we can recover four quantities, including the slip event location, timing, duration, and magnitude. If the estimated location is within 150 km from the true location, if the estimated timing and duration are within five days of the actual ones, and if the estimated magnitude is within 0.25 of the actual one, we consider the detection to be valid. This procedure enables us to determine the minimum magnitude that can be detected over each of the three observation periods and build resolution maps for each period investigated (see Supporting Information, Figures S29-S30). For instance, in the Iquique region (~19°S - 71°W), the minimal magnitude Mw ranges from 6.6 to 6.8 from 2000 to 2003, decreases to 6.1-6.3 from 2004 to 2007 and again down to 5.9 to 6.1 from 2008 to 2014. Thus, as expected, we observe a significant improvement in detection sensitivity when the number of stations in a given region increases.

## 2.5 Application to GNSS time series

After exploring the network sensitivity to detect aseismic slip events, we search for transients in the residual time series obtained after subtracting the trajectory model described earlier. We fix the duration T of the template to 40 days and the slip to an event equivalent to Mw 6.0 (see Supporting Information, Figures S58-S59 for a test in the duration template sensitivity). By doing so, we detect 733 candidate slip events in the stacked correlation functions. Since some of these candidates may correspond to the same candidate slip event, we retain maximum occurrences within a radius of 150 km (i.e., if two maxima affect nodes separated by a distance higher than 150 km, they are considered as independent occurrences). After this selection step, we are left with 59 candidate slip events in the region. We evaluate their durations and magnitudes and compare these with our resolution maps. We keep candidates for which the obtained magnitude is higher than the minimum detectable magnitude for the corresponding node (Figure 4), leaving us with 24 validated slip events.

The duration of the slip events ranges from 17 to 36 days with magnitudes from Mw 5.4 to 6.2 and depths from 20 to 66 km. Figure 3 shows the location of the detected slip events along with four examples of weighted stacks. Figures 5 and 6 show two examples of stacks
<figure><img src="image_6.png" /></figure>

**Figure 5** Example of detected aseismic slip event #12 in the vicinity of the 2014 Iquique earthquake, its locations, and associated seismicity. Figure (a) features the weighted stack for the event #12, with the red line representing the preferred model used to estimate event duration and magnitude, as indicated at the top left. The dark green line denotes the correlation function where event detection is made. Figures (b) and (c) display the displacement time series for the North and East components, respectively. Displacement data from six stations contributing to the weighted stack are shown. The pink lines indicate the best-fitting model for each displacement time series, which incorporates a linear trend and a transient, in accordance with Eq. 5. Meanwhile, the green lines represent the displacements for the estimated magnitude of each event. Figure (d) illustrates the envet location (marked by white star), with dots indicating seismicity before and after the event (spanning half of the event's duration for each period), scaled by magnitude and color-coded by date. Inverted triangles mark the GNSS station locations. Pink arrows denote the GNSS-derived displacements from observations used to estimate the weighted stack during the detected slow slip event, whereas black arrows indicate displacements not used in the estimation. The green arrows show displacements resulting from dislocations for the estimated magnitudes at each event location (white star). Figure (e) displays the map view of the correlation peak within the correlation function (illustrated in dark green in Figures a) for the event, pinpointing the moment when the detection is made.
<figure><img src="image_7.png" /><figcaption>**Figure 6** Same caption as Figure 5, but for event #10.</figcaption></figure>

and correlation functions, along with the time series used to build the stacks and the map view of correlation peaks (see Supplementary information Tables S43 for the event parameters estimated with their uncertainties, and Figures S33 - S43 to see the data employed in the modeling, the data stack, and the model).

Following the methodology proposed by Nishimura et al. (2013), validated events are categorized into two types: probable and possible. This classification is achieved by comparing the displacement fields derived directly from observations with those generated by synthetics events of estimated magnitudes. Note that the magnitudes are estimated on the correlation stack and not directly on the measured displacements. A disagreement between the displacements corresponding to the detected magnitude on the detected node and the observed displacement would suggest our assumptions do not hold. Observed displacements are determined

directly on the GNSS time series by estimating a linear trend along with a time-dependent slip evolution (Eq. 5). To estimate the displacement field for a detected magnitude, the slip corresponding to that magnitude is applied at the inferred location of each event. Figures 5 (b) and (c) illustrate examples of these estimates, with the actual displacements shown in magenta, while the displacements predicted from the magnitudes of each event are shown in green for Events #10 and #12 (see Supplementary Information, Figures S33 - S43 for the rest of the events). Upon analysis, we find that the agreement between observed and modeled ground motion is acceptable for 10 of our events, leading us to classify these as probable (A events, Table S43). Meanwhile, we observe a weaker agreement for 14 events which we hence categorize as possible (B events, Table S43).

Since our template matching approach only considers GNSS observations, we must ensure that the de-
tected slip events (A and B) are mostly aseismic. We cross-check the 24 positive detections with the seismic catalog provided by the ISC (International Seismological Centre, 2016). We randomly generate 10000 synthetic locations for each slip event considering a normally distributed location uncertainty based on our resolution tests and estimate the sum of the seismic moment of all earthquakes occurring within at least a 2-σ radius of the detected slip event. We then compare this estimate of the seismic moment to the estimated aseismic one. All the detected slip events have an equivalent magnitude at least twice larger than the seismic magnitude (aseismic/seismic ratio for each event and further details on ratio estimation are in Supplementary Information, Table S43). Figures 5 and 6 (d) present the location of the two events detailed in Figures 5 and 6 (a) together with the seismicity that coincides with the occurrence of the slip event. These two events occur during the preparation phase of the 2014 Iquique earthquake (Event #12, Figure 1) and during the interseismic phase (Event #10). The combination of synthetic tests and the seismic vs. aseismic moment analysis confirms we detected 24 aseismic slip events (A and B) along southern Peru - northern Chile subduction zone over the period 2006 - 2014.

# 3 Discussion

## 3.1 Aseismic slip events and scaling laws

Aseismic slip events are now frequently observed along most subduction zones in the world, but the underlying physics is still debated. Among the points of debate, the comparison between slow slip and earthquakes should allow to point out whether comparable physics are involved. Ide et al. (2007) have proposed that, while the seismic moment of earthquakes is proportional to the cube of their duration, the moment of slow earthquakes, from tremors and low-frequency earthquakes to slow slip events, is proportional to the duration. Considering that simple considerations about size and stress drop led to the emergence of the observed scaling for earthquakes, the difference in moment-duration scaling should involve a fundamental difference between the mechanics of slow slip and that of earthquakes. Peng and Gomberg (2010) argued that the apparent moment duration scaling of slow earthquakes proposed by Ide et al. (2007) was only due to a lack of observations, suggesting both rapid and slow slip were driven by the same mechanism, namely a slip instability with variable speed and stress drop propagating along a weakened fault surface. In addition, Gomberg et al. (2016) proposed that seismic moment scales either with the duration or the cube of the duration depending on whether the rupture was elongated and pulse-like or mostly crack-like. Michel et al. (2019b) confirmed that the moment of slow slip events in Cascadia scales with the cube of their duration although being elongated and pulse-like. These observations agree with recent studies of aseismic slip and tremors in Japan (Takagi et al., 2019; Supino et al., 2020) and Mexico (Frank and Brodsky, 2019), as well as numerical modeling us-

ing dynamic simulations of frictional sliding (Dal Zilio et al., 2020b). Such numerical and observational evidence suggests that SSEs might exhibit comparable scaling as classical earthquakes, only with lower rupture speeds and stress drops.

We evaluate the scaling between moment and duration for the aseismic slip events we have detected. We estimate that the moment, *M*, is such as $M \propto T^{4.99 \pm 0.48}$, with *T* the duration for the 24 detected SSEs (refer to Figures 7, S45, and S46 in the Supporting Information for an in-depth explanation of the scaling estimation procedure). This scaling relationship remains consistent when analyzing events A ($M \propto T^{5.05 \pm 0.59}$, see Figures S47 and S48) and B ($M \propto T^{4.89 \pm 0.52}$, illustrated in Figures S49 and S50) independently. Our events seem to align with a moment-duration scaling $T^3$. However, as extensively discussed by Ide and Beroza (2023), uncertainties associated with the estimation of event duration might influence significantly our results. Consequently, it is challenging to definitively conclude that our findings adhere to the moment-duration $T^3$ scaling. That said, our detections are situated within the range of moment-duration observed in other subduction zones such as Cascadia, Japan, or Mexico (Ide and Beroza, 2023, and references therein). Building on this observation, we adopt the methodology outlined by Gomberg et al. (2016) to deduce the source properties of our events. We infer that the rupture velocities of our detections range between 0.5 and 10 km/day, accompanied by a stress drop of 0.1 MPa (see the Supporting Information for detailed information on the parameter estimation process). Although our method does not allow to detect events that would propagate, we observe our SSEs are more compatible with crack-like, unbounded ruptures than pulse-like, bounded ones. As a conclusion, our findings along southern Peru - northern Chile region align with SSEs observations from other subduction zones.

## 3.2 Aseismic slip and interseismic coupling distribution

Our coupling estimate corresponds to an average behavior over a decade, without accounting for potential slow slip events hidden within the noise. The slow slip events we detect hence correspond to fluctuations around this average. We compare the map of coupling to the location of our 24 aseismic events to explore how such fluctuations distribute with respect to locked and creeping asperities along the megathrust (Figure 3). We compare the distribution of coupling where our events are located to a distribution coupling at randomly picked locations (Figure 8, see the Supporting Information for a detailed explanation of the calculation of the PDF for coupling and detected events). The distributions differ but mostly when considering only events in northern Chile, where our estimate of coupling is much more robust. Detected slow slip events occur mostly in regions of intermediate coupling. This observation is not as clear for the Peruvian region, probably because of the sparsity of the data used here, although the same tendency is suggested on Figure 8. This result aligns with
<figure><img src="image_8.png" /><figcaption>**Figure 7** Seismic moment versus duration for our aseismic slip events following the scaling law proposed by Gomberg et al. (2016). Slow bounded/unbounded (SBG, SUG) and fast bounded/unbounded (FBG, FUG) regions are shown by light gray areas. Dashed lines are the theoretical relationship between moment and duration for a few selected stress-drop and rupture velocity values. The *M* ∝ *T* scaling is shown in green. The *M* ∝ *T*³ scaling is shown in red.</figcaption></figure>

Frank (2016) findings in the Mexico subduction zone, where a database of slow slip events seems to compensate the lack of slip deficit in transition zones with respect to coupled regions of the megathrust. Materna et al. (2019) describe a comparable behavior over longer periods where coupling variations seem to occur in regions of transitional coupling (Michel et al., 2019a). In addition, events offshore Peru tend to cluster spatially around locked asperities, areas that are generally of intermediate coupling (Figure 9). In general, slow slip events occur in transitional regions between seismic asperities and freely slipping areas. This is consistent with model predictions from rate-and-state friction in which slow slip events are expected to occur at the transition between seismic, rate-weakening and creeping, rate-strengthening asperities (e.g., Liu and Rice, 2005, 2007; Perfettini and Ampuero, 2008).

The average depth of the detected slow slip events is 33 km (Figure 8, see the Supporting Information for a detailed explanation of the PDF calculation). Separating the events, by region, yields an average depth of 37 km for Peru and 30 km for northern Chile with comparable standard deviations (19 and 10 km respectively, Figure 8). This result remains consistent when conducting separate analyses of events A and B (refer to Figures S56-S57 in the Supporting Information). Lay (2015) separates the subduction megathrust along depth into four domains (A, B, C, and D). Domain A, located between the trench and a depth of about 15 km, hosts either tsunami earthquakes or aseismic deformation. Domain B, between approximately 15 and 30 km depth, hosts large megathrust earthquakes. Domain C, between approximately 30 and 50 km depth, hosts in-

intermediate sized earthquakes. At greater depths, Domain D, between 50 and 70 km, hosts slow slip events, tremors, and very low-frequency earthquakes. Our slow slip events mainly occur in Domains C and D. It is understood that small, velocity weakening asperities in Domain C are embedded in conditionally stable regions of the megathrust, prone to host slow slip events. Domain D is dominated by aseismic sliding and potential slip rate variations could explain deeper detections. Therefore, the depth distribution of our events matches regions where slow slip events are expected in a subduction zone context.

Our resolution tests (Figures S24, S29-S30) suggest that it is impossible to capture aseismic slip near the trench, in domain A, with the current GNSS network. However, large, shallow slow slip events have been observed in Japan (Nishimura, 2014; Nishikawa et al., 2019) and New Zealand (Wallace, 2020). Seafloor geodesy might help to detect the occurrence of such large events and potentially for small, cm-scale ones comparable to our aseismic slip events (Araki et al., 2017). Additionally, stress-shadow induces apparent coupling in velocity-weakening regions, especially late in the interseismic period (Hetland and Simons, 2010; Lindsey et al., 2021). For this reason, we also cannot rule out the potential occurrence of aseismic slip event near the trench.

In addition to the depth-dependent segmentation, we observe an along-strike segmentation in the distribution of SSEs. In particular, we observe a lack of events within the rupture area of the 1877 earthquake, within the Arequipa rupture area and other detections gather around locked asperities, like in the doughnut model
<figure><img src="image_9.png" /></figure>

**Figure 8** Coupling, depth, and Vp/Vs ratio of the detected aseismic slip events. (a) Probability Density Functions (PDF) of 1000 coupling models for 24 random picks (gray) and PDF of coupling where 24 aseismic slip events are detected (green), with respective mean (μ) and standard deviations (σ). (b) and (c) are the same as (a) for the Peru region only (gray: random, blue: SSEs) and northern Chile only (gray: random, magenta: events), respectively. (d) PDF of the depths of 24 random events (gray) and aseismic slip events detected in the region (green). (e) and (f) Same as (d) but for Peru (gray: random, blue: events) and Chile (gray, magenta) regions. (g) PDF of the Vp/Vs ratio for the Chilean region (gray, 17 random events), and detected aseismic events in Chile (magenta).

for seismicity (Kanamori, 1981; Schurr et al., 2020). Such configuration is comparable to that of the Japan trench where the asperity that ruptured during the Tohoku earthquake in 2011 overwhelms the simple depth-dependent distribution of behavior from Lay (2015). In particular, Nishikawa et al. (2019) propose that, unlike the Nankai subduction interface which exhibits a depth-dependent segmentation due to a young, warm slab, the megathrust beneath Tohoku is not segmented at depth into four distinct domains. In our area of interest, the subducting slab is older than the Nankai slab and probably colder (Müller et al., 2008), which would explain why the behavior we unravel is not completely consistent with that of Lay (2015) and potentially closer to that of the Japan trench.

As an additional level of complexity, three events coincide with the subduction of the Nazca ridge (14°S, Figures 3 and 9a), six events are located beneath the Mejillones Peninsula (23°S, Figures 3 and 9d), and three events are within the Arica bend (17°S - 19°S, Figures 3 and 9b and c). These morphological structures are anomalies compared to the model proposed by Lay (2015) as they are considered as barriers to the propagation of large earthquakes (Armijo and Thiele, 1990; Comte and Pardo, 1991; Béjar-Pizarro et al., 2010;

Villegas-Lanza et al., 2016; Poli et al., 2017). In these regions, the depth of our detected slow slip events does not match the depth-dependency described by Lay (2015). We can speculate that local geometrical complexities may lead to the occurrence of slow slip events (Romanet et al., 2018) in the case of the subduction of the Nazca Ridge or that the apparent low coupling is the result of multiple slow slip events (Jolivet et al., 2020) in the case of the Arica Bend.

## 3.3 Aseismic slip events before and after large earthquakes

Among all the detected slow slip events, only events #7, and #12 (Figure 3, S36 and 5) do not occur during the steady interseismic period. Event #7 locates in the region struck by the Iquique earthquake in 2014 (Figure 9c, and S36) during the post-seismic relaxation that followed the mainshock (Meng et al., 2015; Hoffmann et al., 2018; Shrivastava et al., 2019) ($M_w$ 6.1 and duration of 28 days in June 2014). Such slow slip events embedded within a post-seismic sequence have already been observed following the Illapel earthquake (Tissandier et al., 2023) and in a completely different setting, following the 2004 Parkfield earthquake, along the San Andreas Fault (Michel et al., 2022).
<figure><img src="image_10.png" /><figcaption>**Figure 9** Zoom over a selection of regions of interest. Gray contours are instrumental ruptures. Yellow contours show reported afterslip. Our aseismic slip events are color-coded by time and scaled by magnitude. Background color shows our Bayesian inference of coupling. Inverted pink triangles are the GNSS stations used in this study. (a) Region struck by the Pisco (2007) and Nazca (1996) earthquakes. Our detections seem to cluster around asperities broken during earthquakes or afterslip regions. (b) Region struck by the Arequipa (2001) earthquake. (c) Region struck by the Iquique earthquake in 2014. Green contours show the preseismic slip reported by Socquet et al. (2017). Events occur around locked interseismic patches or low-coupled regions. (d) Region struck by the Antofagasta (1995) and Tocopilla (2007) earthquakes. Events surround broken asperities or locked interseismic patches, with a cluster beneath Mejillones Peninsula, potentially associated with earthquake afterslip. For citations of instrumental ruptures and afterslip, please refer to Figure 1</figcaption></figure>

Aseismic slip has been recognized as an important element of the earthquake preparation phase (Obara and Kato, 2016; McLaskey, 2019; Kato and Ben-Zion, 2021, and references therein). An 8-month-long slow slip event was reported before the Iquique earthquake in 2014 (Socquet et al., 2017), and event #12 coincides with one of the regions of the megathrust that slipped aseismically during that preparation phase (Figure 9 c). In addition, event #12 occurred where and when intermediate-depth and shallow seismicity synchronized before the Iquique earthquake (Bouchon et al., 2016; Jara et al., 2017) ($M_w$ 6.0 and duration of 30 days in January 2014). Such synchronization of seismicity began in January 2014, lasted for one month, and is interpreted as evidence of a slow, slab-wide deformation process prior to megathrust earthquakes (Bouchon et al., 2016). Furthermore, event #12 is coincident with the transient event reported by Boudin et al. (2021) using a long-base tiltmeter. Our epicentral location differs by ~50km from the one reported by (Boudin et al., 2021),

a difference that can be explained by different modeling strategies and/or uncertainties. We propose that event #12 is linked to the 8-month aseismic slip transient observed preceding the 2014 Iquique earthquake. Such detection suggests the growing instability preceding the Iquique earthquake exhibits a complex spatiotemporal behavior that hides within the noise of the data, in agreement with the hypothesis proposed by Jolivet and Frank (2020) and Twardzik et al. (2022).

## 3.4 Aseismic slip and fluids

Fluids may also play a role in the occurrence of aseismic slip events (Avouac, 2015; Harris, 2017; Jolivet and Frank, 2020, and references therein). Pore pressure affects fault normal stress, hence modify the probability of a slip instability as well as the nucleation size (Liu and Rice, 2007; Avouac, 2015; Bayart et al., 2016; Harris, 2017; Bürgmann, 2018; Jolivet and Frank, 2020; Behr and Bürgmann, 2021). An increase in pore pressure within the fault zone leads to a decrease in normal
<figure><img src="image_11.png" /><figcaption>**Figure 10** Map view of the depth of the continental Moho discontinuity from gravity-derived structural models by Tassara and Echaurren (2012). Magenta stars are the location of our 24 aseismic events. Black lines indicate the location of the profiles shown on the right. Colors indicate the structure at depth (upper and lower crusts, lithospheric mantle, asthenospheric wedge, and oceanic crust). White box indicates the id of events occurring along each profile.</figcaption></figure>

stress, which promotes slip but increases nucleation size, promoting slow slip. We compare our detections to the distribution of the Vp/Vs ratio and to gravity-inferred structural models in the region. We use the Vp/Vs ratio inferred by Comte et al. (2016) for the events located in Northern Chile. Statistically, the 17 aseismic events in northern Chile are not related to a specific Vp/Vs value (Figure 8, see the Supporting Information for a detailed explanation of the PDF calculation). In particular, no slow slip events are found to collocate with high Vp/Vs ratios (Vp/Vs > 1.8) (Comte et al., 2016) (Figure S44).

We also compare the location of our aseismic events to a 3-D density model in the region (Tassara and Echaurren, 2012). Figure 10 shows the location of aseismic events along ten different trench-perpendicular cross sections. The slow slip events are primarily located along the contact between the slab and the overriding lithospheric mantle (Figure 10, see Figure S51 for an analysis of depth uncertainties). This mantle corner is principally hydrated by the dehydration of the subducting slab due to water releasing metamorphic reactions (Peacock, 2001; Rüpke et al., 2004; Comte et al., 2016; Wang et al., 2019; Contreras-Reyes et al., 2021). The fact that our aseismic slip events tend to cluster at depths corresponding to the lithospheric mantle along the megathrust, and not deeper, might imply that fluids may be trapped and accumulate below the continental Moho, an hypothesis that would require further investi-

gations.

## 4 Conclusions

We have systematically analyzed GNSS time series in the region, searching for the occurrence of aseismic slip events with a template matching approach. We find 24 events in the period 2006 - 2014, with durations of 17 - 36 days, magnitudes of $M_w$ 5.4 - 6.2, and located at depths of 20-66 km. These events are mostly aseismic and are observed at all stages of the earthquake cycle, including during post-seismic periods (afterslip, one event), earthquake preparation phase (one event), and interseismic period (22 events). We compare those slow slip occurrence to a wide range of possible models of interseismic coupling based on GNSS and InSAR velocity fields and infer a distribution of coupling along the megathrust.

By conducting a moment-duration scaling analysis, we find that our observations are consistent with values reported in subduction zones globally. We do not find particular correlations with published seismic velocity structures but find that slow slip events cluster around past ruptures and locked asperities, where the megathrust transitions from sliding to locked. Additionally, our events are located in regions of intermediate coupling values and mean depths of 33 km, which match regions where slow slip events occur in the context of subduc-
tion zones.

Some of these events occur on the subduction interface deeper than than the continental MOHO, i.e. where the slab is in contact with the mantle wedge corner where fluids are supposedly trapped. This points toward the influence of fluids as it may explain both their spontaneous triggering and their long duration. However, as some events are found at shallower depth, the involvement of fluids might not be the only explanation. Other mechanisms such as geometrical complexities might be involved but more evidence are required.

The main outcome of this study is that we found numerous aseismic slip events in a place where none were found during the interseismic period before. As a consequence, aseismic slip events may be found elsewhere in subduction zone contexts where experts did not find any event, pending dedicated noise analysis methods. We provide here one piece of evidence supporting the hypothesis proposed by Jolivet and Frank (2020) which states that slow slip happens everywhere and at all times.

## Acknowledgements

The authors thank the Laboratoire International Associé “Montessus de Ballore” (LIA-MB), International Plate Boundary Observatory Chile (IPOC, www.ipoc-network.org), Central Andean Tectonic Observatory Geodetic Array (CAntO, www.tectonics.caltech.edu/resources/continuous_gps.html), and Instituto Geofísico del Perú (www.igp.gob.pe) for making the raw GNSS data available, as well as all researchers, technicians, and students involved in the installation and maintenance of such networks. J. Jara acknowledges a PhD scholarship granted by Chilean National Science Cooperation (CONICYT) through “Becas Chile” Program and MSCA Postdoctoral Fellowship (MSCA-101066069, project ERASMUS) A. Socquet and J. Jara benefited from a grant from ANR (ANR-17-CE31-0002-01, project AtypicSSE). J. Jara, R. Jolivet and A. Socquet acknowledge the funding of the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement 758210, project Geo4D and grant agreement 865963, project DEEPtrigger). R. Jolivet acknowledges funding from the Institut Universitaire de France. D. Comte acknowledges the funding provided by ANID via the projects PIA/ANID grant ACT-192169, and ANID Project AFB180004, AFB220002. The authors would like to thank J. Pina-Valdes, H. Bhat, H. Sanchez-Reyes, S. Michel, C. Vigny, R. Madariaga, and M. Bouchon for all the constructive discussions about this work. We thank the editors W. Xu, E. Williams, and G. Funning for handling the manuscript, and the reviewers P. Poli and J. Freymueller for their valuable comments, which greatly improved the present manuscript.

## Data and code availability

GNSS data are processed using GAMIT software (Herring et al., 2015) (http://geoweb.mit.edu/gg/), while the

reference frame is defined using PYACS (Nocquet, 2018) (https://github.com/JMNocquet/pyacs36). GNSS time series used in this work can be found at: https://doi.org/10.5281/zenodo.7898656. The modeling has been performed using the Classic Slip Inversion library (Jolivet et al., 2015b) (CSI, https://github.com/jolivetr/csi) and AlTar (Minson et al., 2013) (https://github.com/AlTarFramework/altar). All plots are made using Matplotlib (Hunter, 2007) and Cartopy (Office, 2010) Python packages.

## References

Altamimi, Z., Collilieux, X., and Métivier, L. ITRF2008: An improved solution of the international terrestrial reference frame. *Journal of Geodesy*, 85(8):457-473, 2011. doi: 10.1007/s00190-011-0444-4.

Altamimi, Z., Rebischung, P., Métivier, L., and Collilieux, X. ITRF2014: A new release of the International Terrestrial Reference Frame modeling nonlinear station motions. *Journal of Geophysical Research: Solid Earth*, 121(8):6109-6131, 2016. doi: 10.1002/2016JB013098.

Ambraseys, N. Some characteristic features of the Anatolian fault zone. *Tectonophysics*, 9(2-3):143-165, 3 1970. doi: 10.1016/0040-1951(70)90014-4.

Araki, E., Saffer, D. M., Kopf, A. J., Wallace, L. M., Kimura, T., Machida, Y., Ide, S., and Davis, E. Recurring and triggered slow-slip events near the trench at the Nankai Trough subduction megathrust. *Science*, 356(6343):1157-1160, 2017. doi: 10.1126/science.aan3120.

Armijo, R. and Thiele, R. Active faulting in northern Chile: ramp stacking and lateral decoupling along a subduction plate boundary? *Earth and Planetary Science Letters*, 98(1):40-61, 4 1990. doi: 10.1016/0012-821X(90)90087-E.

Avouac, J.-P. From Geodetic Imaging of Seismic and Aseismic Fault Slip to Dynamic Modeling of the Seismic Cycle. *Annual Review of Earth and Planetary Sciences*, 43(1):233-271, 2015. doi: 10.1146/annurev-earth-060614-105302.

Báez, J. C., Leyton, F., Troncoso, C., del Campo, F., Bevis, M., Vigny, C., Moreno, M., Simons, M., Kendrick, E., Parra, H., and Blume, F. The Chilean GNSS Network: Current Status and Progress toward Early Warning Applications. *Seismological Research Letters*, 89 (4):1546-1554, 7 2018. doi: 10.1785/0220180011.

Bayart, E., Svetlizky, I., and Fineberg, J. Slippery but Tough: The Rapid Fracture of Lubricated Frictional Interfaces. *Physical Review Letters*, 116(19):194301, 5 2016. doi: 10.1103/PhysRevLett.116.194301.

Beck, S. L. and Ruff, L. J. Great earthquakes and subduction along the Peru trench. *Physics of the Earth and Planetary Interiors*, 57 (3-4):199-224, 11 1989. doi: 10.1016/0031-9201(89)90112-X.

Behr, W. M. and Bürgmann, R. What's down there? The structures, materials and environment of deep-seated slow slip and tremor. *Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences*, 379(2193): 20200218, 3 2021. doi: 10.1098/rsta.2020.0218.

Béjar-Pizarro, M., Carrizo, D., Socquet, A., Armijo, R., Barrientos, S., Bondoux, F., Bonvalot, S., Campos, J., Comte, D., De Chabalier, J. B., and others. Asperities and barriers on the seismogenic zone in North Chile: state-of-the-art after the 2007 Mw 7.7 Tocopilla earthquake inferred by GPS and InSAR data. *Geophysical Journal International*, 183(1):390-406, 2010. doi: 10.1111/j.1365-246X.2010.04748.x.

Béjar-Pizarro, M., Socquet, A., Armijo, R., Carrizo, D., Genrich, J.,
and Simons, M. Andean structural control on interseismic coupling in the North Chile subduction zone. *Nature Geoscience*, 6 (6):462-467, 6 2013. doi: 10.1038/ngeo1802.

Bevis, M. and Brown, A. Trajectory models and reference frames for crustal motion geodesy. *Journal of Geodesy*, 88(3):283-311, 3 2014. doi: 10.1007/s00190-013-0685-5.

Bock, Y. and Melgar, D. Physical applications of GPS geodesy: A review. *Reports on Progress in Physics*, 79(10):106801, 2016. doi: 10.1088/0034-4885/79/10/106801.

Boehm, J., Werl, B., and Schuh, H. Troposphere mapping functions for GPS and very long baseline interferometry from European Centre for Medium-Range Weather Forecasts operational analysis data. *Journal of Geophysical Research: Solid Earth*, 111(2), 2006. doi: 10.1029/2005JB003629.

Bouchon, M., Marsan, D., Durand, V., Campillo, M., Perfettini, H., Madariaga, R., and Gardonio, B. Potential slab deformation and plunge prior to the Tohoku, Iquique and Maule earthquakes. *Nature Geoscience*, 9(5):380-383, 5 2016. doi: 10.1038/ngeo2701.

Bouchon, M., Marsan, D., Jara, J., Socquet, A., Campillo, M., and Perfettini, H. Suspected Deep Interaction and Triggering Between Giant Earthquakes in the Chilean Subduction Zone. *Geophysical Research Letters*, 45(11):5454-5460, 6 2018. doi: 10.1029/2018GL078350.

Boudin, F., Bernard, P., Meneses, G., Vigny, C., Olcay, M., Tassara, C., Boy, J. P., Aissaoui, E., Métois, M., Satriano, C., Esnoult, M.-F., Nercessian, A., Vallée, M., Vilotte, J.-P., and Brunet, C. Slow slip events precursory to the 2014 Iquique Earthquake, revisited with long-base tilt and GPS records. *Geophysical Journal International*, 228(3):2092-2121, 2021. doi: 10.1093/gji/ggab425.

Brace, W. F. and Byerlee, J. D. Stick-Slip as a Mechanism for Earthquakes. *Science*, 153(3739):990-992, 8 1966. doi: 10.1126/science.153.3739.990.

Bürgmann, R. The geophysics, geology and mechanics of slow fault slip. *Earth and Planetary Science Letters*, 495:112-134, 2018. doi: 10.1016/j.epsl.2018.04.062.

Bürgmann, R., Kogan, M. G., Levin, V. E., Scholz, C. H., King, R. W., and Steblov, G. M. Rapid aseismic moment release following the 5 December, 1997 Kronotsky, Kamchatka, earthquake. *Geophysical Research Letters*, 28(7):1331-1334, 4 2001. doi: 10.1029/2000GL012350.

Bürgmann, R., Kogan, M. G., Steblov, G. M., Hilley, G., Levin, V. E., and Apel, E. Interseismic coupling and asperity distribution along the Kamchatka subduction zone. *Journal of Geophysical Research: Solid Earth*, 110(7):1-17, 2005. doi: 10.1029/2005JB003648.

Cheloni, D., D'Agostino, N., Selvaggi, G., Avallone, A., Fornaro, G., Giuliani, R., Reale, D., Sansosti, E., and Tizzani, P. Aseismic transient during the 2010-2014 seismic swarm: evidence for longer recurrence of M > 6.5 earthquakes in the Pollino gap (Southern Italy)? *Scientific Reports*, 7(1):576, 4 2017. doi: 10.1038/s41598-017-00649-z.

Chlieh, M., De Chabalier, J. B., Ruegg, J. C., Armijo, R., Dmowska, R., Campos, J., and Feigl, K. L. Crustal deformation and fault slip during the seismic cycle in the North Chile subduction zone, from GPS and INSAR observations. *Geophysical Journal International*, 158(2):695-711, 2004. doi: 10.1111/j.1365-246X.2004.02326.x.

Chlieh, M., Perfettini, H., Tavera, H., Avouac, J. P., Remy, D., Nocquet, J. M., Rolandone, F., Bondoux, F., Gabalda, G., and Bonvalot, S. Interseismic coupling and seismic potential along the Central Andes subduction zone. *Journal of Geophysical Research: Solid Earth*, 116(12), 2011. doi: 10.1029/2010JB008166.

Comte, D. and Pardo, M. Reappraisal of great historical earthquakes in the northern Chile and southern Peru seismic gaps.

*Natural Hazards*, 4(1):23-44, 1991. doi: 10.1007/BF00126557.

Comte, D., Carrizo, D., Roecker, S., Ortega-Culaciati, F., and Peyrat, S. Three-dimensional elastic wave speeds in the northern Chile subduction zone: Variations in hydration in the supraslab mantle. *Geophysical Journal International*, 207(2):1080-1105, 2016. doi: 10.1093/gji/ggw318.

Contreras-Reyes, E., Díaz, D., Bello-González, J. P., Slezak, K., Potin, B., Comte, D., Maksymowicz, A., Ruiz, J. A., Osses, A., and Ruiz, S. Subduction zone fluids and arc magmas conducted by lithospheric deformed regions beneath the central Andes. *Scientific Reports*, 11(1):1-12, 2021. doi: 10.1038/s41598-021-02430-9.

Dal Zilio, L., Jolivet, R., and van Dinther, Y. Segmentation of the Main Himalayan Thrust Illuminated by Bayesian Inference of Interseismic Coupling. *Geophysical Research Letters*, 47(4):1-10, 2020a. doi: 10.1029/2019GL086424.

Dal Zilio, L., Lapusta, N., and Avouac, J. P. Unraveling Scaling Properties of Slow-Slip Events. *Geophysical Research Letters*, 47(10), 2020b. doi: 10.1029/2020GL087477.

Delouis, B., Monfret, T., Dorbath, L., Pardo, M., Rivera, L., Comte, D., Haessler, H., Caminade, J. P., Ponce, L., Kausel, E., and Cisternas, A. The Mw= 8.0 antofagasta (northern Chile) earthquake of 30 July 1995: A precursor to the end of the large 1877 gap. *Bulletin of the Seismological Society of America*, 87(2):427-445, 1997.

Dorbath, L., Cisternas, A., and Dorbath, C. Assessment of the size of large and great historical earthquakes in Peru. *Bulletin of the Seismological Society of America*, 80(3):551-576, 1990.

Dragert, H., Wang, K., and James, T. S. A Silent Slip Event on the Deeper Cascadia Subduction Interface. *Science*, 292(5521): 1525-1528, 5 2001. doi: 10.1126/science.1060152.

Ducellier, A., Creager, K. C., and Schmidt, D. A. Detection of Slow Slip Events Using Wavelet Analysis of GNSS Recordings. *Bulletin of the Seismological Society of America*, 112(5):2408-2424, 2022. doi: 10.1785/0120210289.

Duputel, Z., Agram, P. S., Simons, M., Minson, S. E., and Beck, J. L. Accounting for prediction uncertainty when inferring subsurface fault slip. *Geophysical Journal International*, 197(1): 464-482, 4 2014. doi: 10.1093/gji/ggt517.

Duputel, Z., Jiang, J., Jolivet, R., Simons, M., Rivera, L., Ampuero, J. P., Riel, B., Owen, S. E., Moore, A. W., Samsonov, S. V., Ortega Culaciati, F., and Minson, S. E. The Iquique earthquake sequence of April 2014: Bayesian modeling accounting for prediction uncertainty. *Geophysical Research Letters*, 42(19): 7949-7957, 2015. doi: 10.1002/2015GL065402.

Dziewonski, A. M., Chou, T.-A., and Woodhouse, J. H. Determination of earthquake source parameters from waveform data for studies of global and regional seismicity. *Journal of Geophysical Research: Solid Earth*, 86(B4):2825-2852, 4 1981. doi: 10.1029/JB086iB04p02825.

Ekström, G., Nettles, M., and Dziewoński, A. The global CMT project 2004-2010: Centroid-moment tensors for 13,017 earthquakes. *Physics of the Earth and Planetary Interiors*, 200-201:1-9, 6 2012. doi: 10.1016/j.pepi.2012.04.002.

Essing, D. and Poli, P. Spatiotemporal Evolution of the Seismicity in the Alto Tiberina Fault System Revealed by a High-Resolution Template Matching Catalog. *Journal of Geophysical Research: Solid Earth*, 127(10), 10 2022. doi: 10.1029/2022JB024845.

Frank, W. B. Slow slip hidden in the noise: The intermittence of tectonic release. *Geophysical Research Letters*, 43(19):125-10, 2016. doi: 10.1002/2016GL069537.

Frank, W. B. and Brodsky, E. E. Daily measurement of slow slip from low-frequency earthquakes is consistent with ordinary earthquake scaling. *Science Advances*, 5(10):eaaw9386, 10 2019. doi:
10.1126/sciadv.aaw9386.

Gardonio, B., Marsan, D., Socquet, A., Bouchon, M., Jara, J., Sun, Q., Cotte, N., and Campillo, M. Revisiting Slow Slip Events Occurrence in Boso Peninsula, Japan, Combining GPS Data and Repeating Earthquakes Analysis. *Journal of Geophysical Research: Solid Earth*, 123(2):1502–1515, 2 2018. doi: 10.1002/2017JB014469.

Gomberg, J., Wech, A., Creager, K., Obara, K., and Agnew, D. Re-
considering earthquake scaling. *Geophysical Research Letters*,
43(12):6243–6251, 6 2016. doi: 10.1002/2016GL069967.

Graham, S., DeMets, C., Cabral-Cano, E., Kostoglodov, V., Rousset,
B., Walpersdorf, A., Cotte, N., Lasserre, C., McCaffrey, R., and
Salazar-Tlaczani, L. Slow Slip History for the MEXICO Subduc-
tion Zone: 2005 Through 2011. *Pure and Applied Geophysics*,
173(10-11):3445–3465, 2016. doi: 10.1007/s00024-015-1211-x.

Gualandi, A., Michele, C., Serpelloni, E., Chiaraluce, L., Anderlini,
L., Latorre, D., Belardinelli, M. E., and Avouac, J.-P. Aseismic de-
formation associated with an earthquake swarm in the north-
ern Apennines (Italy). *Geophysical Research Letters*, 44(15):
7706–7714, 8 2017. doi: 10.1002/2017GL073687.

Harris, R. A. Large earthquakes and creeping faults. *Reviews of Geophysics*, 55(1):169–198, 2017. doi: 10.1002/2016RG000539.

Hartzell, S. and Langer, C. Importance of model parameterization
in finite fault inversions: application to the 1974 MW 8.0 Peru
earthquake. *Journal of Geophysical Research*, 98(B12), 1993.
doi: 10.1029/93jb02453.

Hayes, G. P., Moore, G. L., Portner, D. E., Hearne, M., Flamme, H., Furtney, M., and Smoczyk, G. M. Slab2, a comprehensive subduction zone geometry model. *Science*, 362(6410):58–61, 10 2018. doi: 10.1126/science.aat4723.

Heki, K., Miyazaki, S., and Tsuji, H. Silent fault slip following an
interplate thrust earthquake at the Japan Trench. *Nature*, 386
(6625):595–598, 4 1997. doi: 10.1038/386595a0.

Herring, T. A., King, R., Floyd, M., and McClusky, S. C. GAMIT Refer-
ence Manual. GPS Analysis at MIT GLOBK, Release 10.6, 2015.

Hetland, E. A. and Simons, M. Post-seismic and interseismic fault creep II: Transient creep and interseismic stress shadows on megathrusts. *Geophysical Journal International*, 181(1):99–112, 2010. doi: 10.1111/j.1365-246X.2009.04482.x.

Hino, R., Inazu, D., Ohta, Y., Ito, Y., Suzuki, S., linuma, T., Osa-
ada, Y., Kido, M., Fujimoto, H., and Kaneda, Y. Was the 2011
Tohoku-Oki earthquake preceded by aseismic preslip? Exam-
ination of seafloor vertical deformation data near the epicen-
ter. *Marine Geophysical Research*, 35(3):181–190, 2014. doi:
10.1007/s11001-013-9208-2.

Hirose, H., Hirahara, K., Kimata, F., Fujii, N., and Miyazaki, S. A slow thrust slip event following the two 1996 Hyuganada Earthquakes beneath the Bungo Channel, southwest Japan. *Geophysical Research Letters*, 26(21):3237–3240, 11 1999. doi: 10.1029/1999GL010999.

Hoffmann, F., Metzger, S., Moreno, M., Deng, Z., Sippl, C., Ortega-Culaciati, F., and Oncken, O. Characterizing Afterslip and Ground Displacement Rate Increase Following the 2014 Iquique-Pisagua Mw8.1 Earthquake, Northern Chile. *Journal of Geophysical Research: Solid Earth*, 123(5):4171–4192, 5 2018. doi: 10.1002/2017JB014970.

Hsu, Y.-J., Bechor, N., Segall, P., Yu, S.-B., Kuo, L.-C., and Ma, K.-F. Rapid afterslip following the 1999 Chi-Chi, Taiwan Earthquake. *Geophysical Research Letters*, 29(16):1–4, 8 2002. doi: 10.1029/2002GL014967.

Hsu, Y. J., Simons, M., Avouac, J. P., Galeteka, J., Sieh, K., Chlieh, M., Natawidjaja, D., Prawirodirdjo, L., and Bock, Y. Frictional afterslip following the 2005 Nias-Simeulue earth-

quake, Sumatra. *Science*, 312(5782):1921–1926, 2006. doi:
10.1126/science.1126960.

Hunter, J. D. Matplotlib: A 2D Graphics Environment. *Computing in Science & Engineering*, 9(3):90–95, 2007. doi: 10.1109/MCSE.2007.55.

Husen, S., Kissling, E., Flueh, E., and Asch, G. Accurate hypocentre determination in the seismogenic zone of the subducting Nazca Plate in northern Chile using a combined on-/offshore network. *Geophysical Journal International*, 138(3):687–701, 1999. doi: 10.1046/j.1365-246X.1999.00893.x.

Ide, S. and Beroza, G. C. Slow earthquake scaling reconsidered as a boundary between distinct modes of rupture propagation. *Proceedings of the National Academy of Sciences*, 120(32):2017, 8 2023. doi: 10.1073/pnas.2222102120.

Ide, S., Beroza, G. C., Shelly, D. R., and Uchide, T. A scaling law for slow earthquakes. *Nature*, 447(7140):76–79, 2007. doi: 10.1038/nature05780.

International Seismological Centre. On-line Bulletin, 2016.

Itoh, Y., Aoki, Y., and Fukuda, J. Imaging evolution of Cascadia slow-slip event using high-rate GPS. *Scientific Reports*, 12(1): 1–12, 2022. doi: 10.1038/s41598-022-10957-8.

Jara, J., Socquet, A., Marsan, D., and Bouchon, M. Long-Term Interactions Between Intermediate Depth and Shallow Seismicity in North Chile Subduction Zone. *Geophysical Research Letters*, 44(18):9283–9292, 9 2017. doi: 10.1002/2017GL075029.

Jara, J., Sánchez-Reyes, H., Socquet, A., Cotton, F., Virieux, J., Maksymowicz, A., Díaz-Mojica, J., Walpersdorf, A., Ruiz, J., Cotte, N., and Norabuena, E. Kinematic study of Iquique 2014 M w 8.1 earthquake: Understanding the segmentation of the seismogenic zone. *Earth and Planetary Science Letters*, 503: 131–143, 2018. doi: 10.1016/j.epsl.2018.09.025.

Jolivet, R. and Frank, W. B. The Transient and Intermittent Nature of Slow Slip. *AGU Advances*, 1(1), 2020. doi: 10.1029/2019av000126.

Jolivet, R. and Simons, M. A Multipixel Time Series Analysis Method Accounting for Ground Motion, Atmospheric Noise, and Orbital Errors. *Geophysical Research Letters*, 45(4):1814–1824, 2018. doi: 10.1002/2017GL076533.

Jolivet, R., Candela, T., Lasserre, C., Renard, F., Klinger, Y., and Doin, M. The Burst-Like Behavior of Aseismic Slip on a Rough Fault: The Creeping Section of the Haiyuan Fault, China. *Bulletin of the Seismological Society of America*, 105(1):480–488, 2 2015a. doi: 10.1785/0120140237.

Jolivet, R., Simons, M., Agram, P. S., Duputel, Z., and Shen, Z. K. Aseismic slip and seismogenic coupling along the central San Andreas Fault. *Geophysical Research Letters*, 42(2):297–306, 2015b. doi: 10.1002/2014GL062222.

Jolivet, R., Simons, M., Duputel, Z., Olive, J. A., Bhat, H. S., and Bletery, Q. Interseismic Loading of Subduction Megathrust Drives Long-Term Uplift in Northern Chile. *Geophysical Research Letters*, 47(8):1–11, 2020. doi: 10.1029/2019GL085377.

Kanamori, H. The Nature of Seismicity Patterns Before Large Earth-
quakes. In *Earthquake Prediction*, pages 1–19. Wiley Online Li-
brary, 3 1981. doi: 10.1029/ME004p0001.

Kato, A. and Ben-Zion, Y. The generation of large earthquakes. *Nature Reviews Earth & Environment*, 2(1):26–39, 1 2021. doi: 10.1038/s43017-020-00108-w.

Kato, A. and Nakagawa, S. Multiple slow-slip events during a foreshock sequence of the 2014 Iquique, Chile Mw 8.1 earthquake. *Geophysical Research Letters*, 41(15):5420–5427, 2014. doi: 10.1002/2014GL061138.

Kato, A., Obara, K., Igarashi, T., Tsuruoka, H., Nakagawa, S., and Hi-
rata, N. Propagation of slow slip leading up to the 2011 Mw9.0
Tohoku-Oki earthquake. *Science*, 335(6069):705–708, 2 2012.
doi: 10.1126/science.1215141.

Kausel, E. Los terremotos de agosto de 1868 y mayo de 1877 que afectaron el sur del Perú y norte de Chile. *Boletín de la Academia Chilena de Ciencias*, 3(1):8-13, 1986.

Khoshmanesh, M. and Shirzaei, M. Episodic creep events on the San Andreas Fault caused by pore pressure variations. *Nature Geoscience*, 11(8):610–614, 2018. doi: 10.1038/s41561-018-0160-2.

Klein, E., Vigny, C., Nocquet, J. M., and Boulze, H. A 20 year-long GNSS solution across South-America with focus in Chile. *BSGF- Earth Sciences Bulletin*, 193, 2022. doi: 10.1051/bsgf/2022005.

Klotz, J., Deng, Z., Moreno, M., Asch, G., Bartsch, M., and Ramaatschi, M. IPOC cGPS- Continuous Mode GPS data in the IPOC Region, Northern Chile. Technical report, GFZ Data Services, 2017.

Lay, T. The surge of great earthquakes from 2004 to 2014. *Earth and Planetary Science Letters*, 409(October 2016):133-146, 2015.
doi: 10.1016/j.epsl.2014.10.047.

Li, Y., Nocquet, J. M., Shan, X., and Song, X. Geodetic Observations of Shallow Creep on the Laohushan-Haiyuan Fault, Northeastern Tibet. *Journal of Geophysical Research: Solid Earth*, 126(6): 1–18, 2021. doi: 10.1029/2020JB021576.

Lindsey, E. O., Mallick, R., Hubbard, J. A., Bradley, K. E., Almeida, R. V., Moore, J. D. P., Bürgmann, R., and Hill, E. M. Slip rate deficit and earthquake potential on shallow megathrusts. *Nature Geoscience*, 14(5):321–326, 5 2021. doi: 10.1038/s41561-021-00736-x.

Liu, Y. and Rice, J. R. Aseismic slip transients emerge spontaneously in three-dimensional rate and state modeling of subduction earthquake sequences. *Journal of Geophysical Research: Solid Earth*, 110(8):1–14, 2005. doi: 10.1029/2004JB003424.

Liu, Y. and Rice, J. R. Spontaneous and triggered aseismic deformation transients in a subduction fault model. *Journal of Geophysical Research: Solid Earth*, 112(9), 2007. doi: 10.1029/2007JB004930.

Louderback, G. Faults and Eartquakes. *Bulletin of Seismological Society of America*, 32(4):305–330, 1942.

Loveless, J. P. and Meade, B. J. Geodetic imaging of plate motions, slip rates, and partitioning of deformation in Japan. *Journal of Geophysical Research*, 115(B2):B02410, 2010. doi: 10.1029/2008JB006248.

Lovery, B., Chlieh, M., Norabuena, E., Villegas-Lanza, J. C., Radiguet, M., Cotte, N., Tsapong-Tsague, A., Quiroz, W., Sierra Farfán, C., Simons, M., Nocquet, J. M., Tavera, H., and Socquet, A. Heterogeneous Locking and Earthquake Potential on the South Peru Megathrust From Dense GNSS Network. *Journal of Geophysical Research: Solid Earth*, 129(2), 2 2024. doi: 10.1029/2023JB027114.

Marsan, D., Reverso, T., Helmstetter, A., and Enescu, B. Slow slip and aseismic deformation episodes associated with the subducting Pacific plate offshore Japan, revealed by changes in seismicity. *Journal of Geophysical Research E: Planets*, 118(9): 4900–4909, 2013. doi: 10.1002/jgrb.50323.

Marsan, D., Bouchon, M., Gardonio, B., Perfettini, H., Socquet, A., and Enescu, B. Change in seismicity along the Japan trench, 1990-2011, and its relationship with seismic coupling. *Journal of Geophysical Research: Solid Earth*, 122(6):4645–4659, 6 2017. doi: 10.1002/2016JB013715.

Materna, K., Bartlow, N., Wech, A., Williams, C., and Bürgmann, R. Dynamically Triggered Changes of Plate Interface Coupling in Southern Cascadia. *Geophysical Research Letters*, 46(22):

12890–12899, 2019. doi: 10.1029/2019GL084395.

Mazzotti, S. S., Le Pichon, X., Henry, P., and Miyazaki, S.-I. Full interseismic locking of the Nankai and Japan-west Kurile subduction zones: An analysis of uniform elastic strain accumulation in Japan constrained by permanent GPS. *Journal of Geophysical Research: Solid Earth*, 105(B6):13159–13177, 2000. doi: 10.1029/2000jb900060.

McLaskey, G. C. Earthquake Initiation From Laboratory Observations and Implications for Foreshocks. *Journal of Geophysical Research: Solid Earth*, 124(12):12882–12904, 2019. doi: 10.1029/2019JB018363.

Melbourne, T. I. Precursory transient slip during the 2001 M w = 8.4 Peru earthquake sequence from continuous GPS. *Geophysical Research Letters*, 29(21):2032, 2002. doi: 10.1029/2002GL015533.

Melnick, D., Moreno, M., Quinteros, J., Baez, J. C., Deng, Z., Li, S., and Oncken, O. The super-interseismic phase of the megathrust earthquake cycle in Chile. *Geophysical Research Letters*, 44(2): 784–791, 2017. doi: 10.1002/2016gl071845.

Meng, L., Huang, H., Bürgmann, R., Ampuero, J. P., and Strader, A. Dual megathrust slip behaviors of the 2014 Iquique earthquake sequence. *Earth and Planetary Science Letters*, 411:177–187, 2015. doi: 10.1016/j.epsl.2014.11.041.

Métois, M., Vigny, C., and Socquet, A. Interseismic Coupling, Megathrust Earthquakes and Seismic Swarms Along the Chilean Subduction Zone (38°–18°S). *Pure and Applied Geophysics*, 173 (5):1431–1449, 2016. doi: 10.1007/s00024-016-1280-5.

Michel, S., Gualandi, A., and Avouac, J. P. Interseismic Coupling and Slow Slip Events on the Cascadia Megathrust. *Pure and Applied Geophysics*, 176(9):3867–3891, 2019a. doi: 10.1007/s00024-018-1991-x.

Michel, S., Gualandi, A., and Avouac, J.-P. Similar scaling laws for earthquakes and Cascadia slow-slip events. *Nature*, 574(7779): 522–526, 10 2019b. doi: 10.1038/s41586-019-1673-6.

Michel, S., Jolivet, R., Lengliné, O., Gualandi, A., Larochelle, S., and Gardonio, B. Searching for Transient Slow Slips Along the San Andreas Fault Near Parkfield Using Independent Component Analysis. *Journal of Geophysical Research: Solid Earth*, 127 (6):1–19, 2022. doi: 10.1029/2021JB023201.

Minson, S. E., Simons, M., and Beck, J. L. Bayesian inversion for finite fault earthquake source models I-theory and algorithm. *Geophysical Journal International*, 194(3):1701–1726, 2013. doi: 10.1093/gji/ggt180.

Müller, R. D., Sdrolias, M., Gaina, C., and Roest, W. R. Age, spreading rates, and spreading asymmetry of the world's ocean crust. *Geochemistry, Geophysics, Geosystems*, 9(4), 4 2008. doi: 10.1029/2007GC001743.

Nishikawa, T., Matsuzawa, T., Ohta, K., Uchida, N., Nishimura, T., and Ide, S. The slow earthquake spectrum in the Japan Trench illuminated by the S-net seafloor observatories. *Science*, 365 (6455):808–813, 2019. doi: 10.1126/science.aax5618.

Nishimura, T. Short-term slow slip events along the Ryukyu Trench, southwestern Japan, observed by continuous GNSS. *Progress in Earth and Planetary Science*, 1(1):1–13, 2014. doi: 10.1186/s40645-014-0022-5.

Nishimura, T., Matsuzawa, T., and Obara, K. Detection of short-term slow slip events along the Nankai Trough, southwest Japan, using GNSS data. *Journal of Geophysical Research: Solid Earth*, 118 (6):3112–3125, 2013. doi: 10.1002/jgrb.50222.

Nocquet, J. M. PYACS: A set of Python tools for GPS analysis and tectonic modelling. In *PYACS: A set of Python tools for GPS analysis and tectonic modelling*. 19th General Assembly of Wegener, 2018.
Nocquet, J. M., Villegas-Lanza, J. C., Chlieh, M., Mothes, P. A., Rolan-
done, F., Jarrin, P., Cisneros, D., Alvarado, A., Audin, L., Bondoux,
F., Martin, X., Font, Y., Régnier, M., Vallée, M., Tran, T., Beauval,
C., Maguiña Mendoza, J. M., Martinez, W., Tavera, H., and Yepes,
H. Motion of continental slivers and creeping subduction in the
northern Andes. *Nature Geoscience*, 7(4):287-291, 2014. doi:
10.1038/ngeo2099.

Obara, K. and Kato, A. Connecting slow earthquakes to huge earth-
quakes. *Science (New York, N.Y.)*, 353(6296):253–257, 2016. doi:
10.1126/science.aaf1512.

Office, M. Cartopy: a cartographic python library with a Matplotlib
interface, 2010.

Peacock, S. M. Are the lower planes of double seismic zones
caused by serpentine dehydration in subducting oceanic man-
tle? *Geology*, 29(4):299-302, 2001. doi: 10.1130/0091-
7613(2001)029<0299:ATLPOD>2.0.CO;2

Peng, Z. and Gomberg, J. An integrated perspective of the contin-
uum between earthquakes and slow-slip phenomena. *Nature*
*Geoscience*, 3(9):599-607, 2010. doi: 10.1038/ngeo940.

Perfettini, H. and Ampuero, J. P. Dynamics of a velocity strength-
ening fault region: Implications for slow earthquakes and post-
seismic slip. *Journal of Geophysical Research: Solid Earth*, 113
(9), 2008. doi: 10.1029/2007JB005398.

Perfettini, H., Avouac, J. P., Tavera, H., Kositsky, A., Nocquet, J. M., Bondoux, F., Chlieh, M., Sladen, A., Audin, L., Farber, D. L., and Soler, P. Seismic and aseismic slip on the Central Peru megathrust. *Nature*, 465(7294):78–81, 2010. doi: 10.1038/nature09062.

Peyrat, S. and Favreau, P. Kinematic and spontaneous rupture
models of the 2005 Tarapacá intermediate depth earthquake.
*Geophysical Journal International*, 181(1):369-381, 2010. doi:
10.1111/j.1365-246X.2009.04493.x.

Peyrat, S., Campos, J., de Chabalier, J. B., Perez, A., Bonvalot, S., Bouin, M. P., Legrand, D., Nercessian, A., Charade, O., Patau, G., Clévédæ, E., Kausel, E., Bernard, P., and Vilotte, J. P. Tarapacá intermediate-depth earthquake (Mw 7.7, 2005, northern Chile): A slab-pull event with horizontal fault plane constrained from seismologic and geodetic observations. *Geophysical Research Letters*, 33(22):1–6, 2006. doi: 10.1029/2006GL027710.

Poli, P., Jeria, A. M., and Ruiz, S. The M w 8.3 Illapel earthquake (Chile): Preseismic and postseismic activity associated with hydrated slab structures. *Geology*, 45(3):247–250, 3 2017. doi: 10.1130/G38522.1.

Pritchard, M. E. and Simons, M. An aseismic slip pulse in north-
ern Chile and along-strike variations in seismogenic behavior.
*Journal of Geophysical Research: Solid Earth*, 111(8), 2006. doi:
10.1029/2006JB004258.

Pritchard, M. E., Norabuena, E. O., Ji, C., Boroschek, R., Comte, D., Simons, M., Dixon, T. H., and Rosen, P. A. Geodetic, teleseismic, and strong motion constraints on slip from recent southern Peru subduction zone earthquakes. *Journal of Geophysical Research: Solid Earth*, 112(3), 2007. doi: 10.1029/2006JB004294.

Radiguet, M., Cotton, F., Vergnolle, M., Campillo, M., Walpers-
dorf, A., Cotte, N., and Kostogllov, V. Slow slip events
and strain accumulation in the Guerrero gap, Mexico. *Jour-
nal of Geophysical Research: Solid Earth*, 117(4), 2012. doi:
10.1029/2011JB008801.

Radiguet, M., Perfettini, H., Cotte, N., Gualandi, A., Valette, B., Kos-
togllov, V., Lhomme, T., Walpersdorf, A., Cabral Cano, E., and
Campillo, M. Triggering of the 2014 Mw7.3 Papanoa earthquake
by a slow slip event in Guerrero, Mexico. *Nature Geoscience*, 9
(11):829-833, 2016. doi: 10.1038/ngeo2817.

Reid, H. F. The Mechanism of the Earthquake. The California Earth-
quake of April 18, 1906: Rep. of the State Investigation Commiss.
Vol. 2. P. 1. Technical report, Carnigie Institution of Washington,

Washington, D. C., 1910.

Remy, D., Perfettini, H., Cotte, N., Avouac, J. P., Chlieh, M., Bon-
doux, F., Sladen, A., Tavera, H., and Socquet, A. Postseismic re-
locking of the subduction megathrust following the 2007 Pisco,
Peru, earthquake. *Journal of Geophysical Research: Solid Earth*,
121(5):3978-3995, 5 2016. doi: 10.1002/2015JB012417.

Reverso, T., Marsan, D., Helmstetter, A., and Enescu, B. Background seismicity in Boso Peninsula, Japan: Long-term acceleration, and relationship with slow slip events. *Geophysical Research Letters*, 43(11):5671-5679, 2016. doi: 10.1002/2016GL068524.

Romanet, P., Bhat, H. S., Jolivet, R., and Madariaga, R. Fast and Slow Slip Events Emerge Due to Fault Geometrical Complexity. *Geophysical Research Letters*, 45(10):4809-4819, 2018. doi: 10.1029/2018GL077579.

Rousset, B., Campillo, M., Lasserre, C., Frank, W. B., Cotte, N., Walpersdorf, A., Socquet, A., and Kostogllov, V. A geodetic matched filter search for slow slip with application to the Mexico subduction zone. *Journal of Geophysical Research: Solid Earth*, 122(12):498-10, 12 2017. doi: 10.1002/2017JB014448.

Rousset, B., Bürgmann, R., and Campillo, M. Slow slip events in the
roots of the San Andreas fault. *Science Advances*, 5(2):eaav3274,
2 2019. doi: 10.1126/sciadv.aav3274.

Ruegg, J. C., Olcay, M., and Lazo, D. Co-, Post- and Pre(?)-seismic Displacements Associated with the Mw 8.4 Southern Peru Earthquake of 23 June 2001 from Continuous GPS Measurements. *Seismological Research Letters*, 72(6):673-678, 11 2001. doi: 10.1785/gssrl.72.6.673.

Ruiz, S. and Madariaga, R. Historical and recent large megathrust earthquakes in Chile. *Tectonophysics*, 733(September 2017): 37-56, 2018. doi: 10.1016/j.tecto.2018.01.015.

Ruiz, S., Metois, M., Fuenzalida, A., Ruiz, J., Leyton, F., Grandin,
R., Vigny, C., Madariaga, R., and Campos, J. Intense fore-
shocks and a slow slip event preceded the 2014 Iquique
Mw8.1 earthquake. *Science*, 345(6201):1165-1169, 2014. doi:
10.1126/science.1256074.

Ruiz, S., Klein, E., del Campo, F., Rivera, E., Poli, P., Metois, M., Christophe, V., Baez, J. C., Vargas, G., Leyton, F., Madariaga, R., and Fleitout, L. The Seismic Sequence of the 16 September 2015 M w 8.3 Illapel, Chile, Earthquake. *Seismological Research Letters*, 87(4):789-799, 7 2016. doi: 10.1785/0220150281.

Ruiz, S., Moreno, M., Melnick, D., del Campo, F., Poli, P., Baez, J. C., Leyton, F., and Madariaga, R. Reawakening of large earthquakes in south central Chile: The 2016 Mw7.6 Chiloé event. *Geophysical Research Letters*, 44(13):6633-6640, 7 2017. doi: 10.1002/2017GL074133.

Rüpke, L. H., Morgan, J. P., Hort, M., and Connolly, J. A. Serpentine and the subduction zone water cycle. *Earth and Planetary Science Letters*, 223(1-2):17-34, 2004. doi: 10.1016/j.epsl.2004.04.018.

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

Schurr, B., Asch, G., Hainzl, S., Bedford, J., Hoechner, A., Palo,
M., Wang, R., Moreno, M., Bartsch, M., Zhang, Y., Oncken, O.,
Tilmann, F., Dahm, T., Victor, P., Barrientos, S., and Vilotte, J.-P.
Gradual unlocking of plate boundary controlled initiation of the
2014 Iquique earthquake. *Nature*, 512(7514):299-302, 8 2014.
doi: 10.1038/nature13681.

Schurr, B., Moreno, M., Tréhu, A. M., Bedford, J., Kummerow, J., Li, S., and Oncken, O. Forming a Mogi Doughnut in the Years Prior to and Immediately Before the 2014 M8.1 Iquique, Northern Chile, Earthquake. *Geophysical Research Letters*, 47(16), 2020. doi: 10.1029/2020GL088351.
Shrivastava, M. N., González, G., Moreno, M., Soto, H., Schurr, B., Salazar, P., and Báez, J. C. Earthquake segmentation in northern Chile correlates with curved plate geometry. *Scientific Reports*, 9(1):4403, 12 2019. doi: 10.1038/s41598-019-40282-6.

Simons, M., Galetzka, J. E., Genrich, J. F., Ortega, F., Comte, D., Glass, B., Gonzalez, G., and Norabuena, E. Central Andean Tectonic Observatory Geodetic Array - GPS/GNSS Observations. Technical report, Caltech, 2010.

Sippl, C., Schurr, B., Asch, G., and Kummerow, J. Seismicity Structure of the Northern Chile Forearc From >100,000 Double-Difference Relocated Hypocenters. *Journal of Geophysical Research: Solid Earth*, 123(5):4063-4087, 2018. doi: 10.1002/2017JB015384.

Sippl, C., Schurr, B., Münchmeyer, J., Barrientos, S., and Oncken, O. The Northern Chile forearc constrained by 15 years of permanent seismic monitoring. *Journal of South American Earth Sciences*, 126(December 2022):104326, 6 2023. doi: 10.1016/j.jsames.2023.104326.

Sladen, A., Tavera, H., Simons, M., Avouac, J. P., Konca, A. O., Perfettini, H., Audin, L., Fielding, E. J., Ortega, F., and Cavagnoud, R. Source model of the 2007 Mw8.0 Pisco, Peru earthquake: Implications for seismogenic behavior of subduction megathrusts. *Journal of Geophysical Research: Solid Earth*, 115(2), 2010. doi: 10.1029/2009JB006429.

Socquet, A., Valdes, J. P., Jara, J., Cotton, F., Walpersdorf, A., Cotte, N., Specht, S., Ortega-Culaciati, F., Carrizo, D., and Norabuena, E. An 8 month slow slip event triggers progressive nucleation of the 2014 Chile megathrust. *Geophysical Research Letters*, 44(9): 4046-4053, 5 2017. doi: 10.1002/2017GL073023.

Steinbrugge, K. V., Zacher, E. G., Tocher, D., Whitten, C. A., and Claire, C. N. Creep on the San Andreas fault. *Bulletin of the Seismological Society of America*, 50(3):389-415, 7 1960.

Supino, M., Poiata, N., Festa, G., Vilotte, J. P., Satriano, C., and Obara, K. Self-similarity of low-frequency earthquakes. *Scientific Reports*, 10(1):6523, 12 2020. doi: 10.1038/s41598-020-63584-6.

Takagi, R., Uchida, N., and Obara, K. Along-Strike Variation and Migration of Long-Term Slow Slip Events in the Western Nankai Subduction Zone, Japan. *Journal of Geophysical Research: Solid Earth*, 124(4):3853-3880, 2019. doi: 10.1029/2018JB016738.

Tarantola, A. *Inverse Problem Theory and Methods for Model Parameter Estimation*. SIAM, 2005. doi: 10.1137/1.9780898717921.

Tassara, A. and Echaurren, A. Anatomy of the Andean subduction zone: three-dimensional density model upgraded and compared against global-scale models. *Geophysical Journal International*, 189(1):161-168, 4 2012. doi: 10.1111/j.1365-246X.2012.05397.x.

Teunissen, P. J. and Montenbruck, O., editors. *Springer Handbook of Global Navigation Satellite Systems*. Springer International Publishing, Cham, 2017. doi: 10.1007/978-3-319-42928-1.

Tissandier, R., Nocquet, J., Klein, E., Vigny, C., Ojeda, J., and Ruiz, S. Afterslip of the M w 8.3 2015 Illapel Earthquake Imaged Through a Time-Dependent Inversion of Continuous and Survey GNSS Data. *Journal of Geophysical Research: Solid Earth*, 128(2):1-21, 2 2023. doi: 10.1029/2022JB024778.

Twardzik, C., Duputel, Z., Jolivet, R., Klein, E., and Rebischung, P. Bayesian inference on the initiation phase of the 2014 Iquique, Chile, earthquake. *Earth and Planetary Science Letters*, 600: 117835, 2022. doi: 10.1016/j.epsl.2022.117835.

Uchida, N., Takagi, R., Asano, Y., and Obara, K. Migration of shallow and deep slow earthquakes toward the locked segment of the Nankai megathrust. *Earth and Planetary Science Letters*, 531: 115986, 2020. doi: 10.1016/j.epsl.2019.115986.

van Rijsingen, E. M., Calais, E., Jolivet, R., de Chabalier, J., Jara, J., Symithe, S., Robertson, R., and Ryan, G. A. Inferring Interseismic Coupling Along the Lesser Antilles Arc: A Bayesian Approach. *Journal of Geophysical Research: Solid Earth*, 126(2): 1-21, 2 2021. doi: 10.1029/2020JB020677.

Vigny, C. and Klein, E. The 1877 megathrust earthquake of North Chile two times smaller than thought? A review of ancient articles. *Journal of South American Earth Sciences*, 117:103878, 8 2022. doi: 10.1016/j.jsames.2022.103878.

Villegas-Lanza, J. C., Chlieh, M., Cavalié, O., Tavera, H., Baby, P., Chire-Chira, J., and Nocquet, J.-M. Active tectonics of Peru: Heterogeneous interseismic coupling along the Nazca megathrust, rigid motion of the Peruvian Sliver, and Subandean shortening accommodation. *Journal of Geophysical Research: Solid Earth*, 121(10):7371-7394, 10 2016. doi: 10.1002/2016JB013080.

Voss, N., Dixon, T. H., Liu, Z., Malservisi, R., Protti, M., and Schwartz, S. Do slow slip events trigger large and great megathrust earthquakes? *Science advances*, 4(10):eaat8472, 10 2018. doi: 10.1126/sciadv.aat8472.

Wallace, L. M. Slow Slip Events in New Zealand. *Annual Review of Earth and Planetary Sciences*, 48(1):175-203, 5 2020. doi: 10.1146/annurev-earth-071719-055104.

Wang, H., Huismans, R. S., and Rondenay, S. Water Migration in the Subduction Mantle Wedge: A Two-Phase Flow Approach. *Journal of Geophysical Research: Solid Earth*, 124(8):9208-9225, 2019. doi: 10.1029/2018JB017097.

Williams, S. D. P. The effect of coloured noise on the uncertainties of rates estimated from geodetic time series. *Journal of Geodesy*, 76(9-10):483-494, 2003. doi: 10.1007/s00190-002-0283-4.

Zhu, L. and Rivera, L. A. A note on the dynamic and static displacements from a point source in multilayered media. *Geophysical Journal International*, 148(3):619-627, 3 2002. doi: 10.1046/j.1365-246X.2002.01610.x.

The article *Detection of slow slip events along the southern Peru - northern Chile subduction zone* © 2024 by J. Jara is licensed under CC BY 4.0.