# Heterogeneous Locking and Earthquake Potential on the South Peru Megathrust From Dense GNSS Network

B. Lovery, M. Chlieh, E. Norabuena, J C Villegas-lanza, M. Radiguet, N. Cotte, A Tsapong-tsague, W. Quiroz, C. Sierra Farfán, M. Simons, et al.

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

## To cite this version:

B. Lovery, M. Chlieh, E. Norabuena, J C Villegas-lanza, M. Radiguet, et al., Heterogeneous Locking and Earthquake Potential on the South Peru Megathrust From Dense GNSS Network. Journal of Geophysical Research: Solid Earth, 2024, 129 (2), pp.e2023JB027114. 10.1029/2023JB027114. hal-04976118

HAL Id: hal-04976118

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

Submitted on 4 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
JGR Solid Earth

RESEARCH ARTICLE

10.1029/2023JB027114

**Key Points:**

* We present a dense interseismic velocity field at the scale of the South Peruvian Andes, from new decadal global navigation satellite system data at 73 locations

* Low locking (~0.4) is estimated along the Nazca Ridge and the Nazca Fracture Zone, delimiting wide patches of high locking (~0.9)

* Moment budget analysis shows that the South Peru segment could host a $M_w = 8.4–9.0$ earthquake with a 100 to 1,000 years recurrence time

**Supporting Information:**

Supporting Information may be found in the online version of this article.

**Correspondence to:**

B. Lovery,
bertrand.lovery@univ-grenoble-alpes.fr

**Citation:**

Lovery, B., Chlieh, M., Norabuena, E., Villegas-Lanza, J. C., Radiguet, M., Cotte, N., et al. (2024). Heterogeneous locking and earthquake potential on the South Peru megathrust from dense GNSS network. *Journal of Geophysical Research: Solid Earth, 129*, e2023JB027114. https://doi.org/10.1029/2023JB027114

Received 17 MAY 2023
Accepted 12 JAN 2024

**Author Contributions:**

**Conceptualization:** B. Lovery, M.

Chlieh, M. Radiguet, A. Socquet

**Data curation:** B. Lovery, E. Norabuena, J. C. Villegas-Lanza, N. Cotte, A.

Tsapong-Tsague, C. Sierra Farfán, M.

Simons, J. M. Nocquet, H. Tavera, A.

Socquet

**Formal analysis:** B. Lovery, M. Chlieh

**Funding acquisition:** M. Chlieh, A.

Socquet

© 2024. The Authors.

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

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

# Heterogeneous Locking and Earthquake Potential on the South Peru Megathrust From Dense GNSS Network

B. Lovery¹, D., M. Chlieh¹, E. Norabuena², J. C. Villegas-Lanza² D, M. Radiguet¹ D, N. Cotte¹ D, A. Tsapong-Tsague¹ D, W. Quiroz² D, C. Sierra Farfán³, M. Simons⁴ D, J. M. Nocquet⁵ D, H. Tavera², and A. Socquet¹ D

¹CNRS, IRD, ISTERRE, Université Grenoble Alpes, Université Savoie Mont Blanc, Université Gustave Eiffel, Grenoble,

France, ²Instituto Geofísico del Perú, Lima, Peru, ³Instituto Geográfico Nacional, Lima, Peru, ⁴Caltech, Pasadena, CA, USA,

⁵Institut de Physique du Globe de Paris, IRD, Paris, France

**Abstract** The Central Andes subduction has been the theater of numerous large earthquakes since the beginning of the 21st Century, notably the 2001 $M_w = 8.4$ Arequipa, 2007 $M_w = 8.0$ Pisco and 2014 $M_w = 8.1$ Iquique earthquakes. We present an analysis of 47 permanent and 26 survey global navigation satellite system (GNSS) measurements acquired in Central-South Peru between 2007 and 2022 to better understand the frictional properties of the megathrust interface. Using a trajectory model that mimics the different phases of the cycle, we extract a coherent interseismic GNSS field at the scale of the Central Andes from Lima to Arica (12–18.5°S). Interseismic models on a 3D slab geometry indicate that the locking level is relatively high and concentrated between 20 and 40-km depth. Locking distributions indicate a high spatial variability of the coupling along the trench, with the presence of many locked patches that spatially correlate with the seismotectonic segmentation. Our study confirms the presence of a creeping segment where the Nazca Ridge is subducting; we also observe a lighter apparent decrease of coupling related to the Nazca Fracture Zone (NFZ). However, since the Nazca Ridge appears to behave as a strong barrier, the NFZ is less efficient to arrest seismic rupture propagation. Considering various uncertainty factors, we discuss the implication of our coupling estimates with size and timing of large megathrust earthquakes considering both deterministic and probabilistic approaches. We estimate that the South Peru segment could have a $M_w = 8.4–9.0$ earthquake potential depending principally on the considered seismic catalog and the seismic/aseismic slip ratio.

**Plain Language Summary** Using dense global navigation satellite system (GNSS) data collected in the South-Central Peru, we extracted a large scale interseismic velocity (surface velocity between two earthquakes) field at the scale of the Central Andes of Peru, where the oceanic Nazca plate goes under the continental South America plate at a velocity of about 6 cm/yr. This area has been the theater of several great subduction earthquakes and tsunamis, then estimating the stress build-up on the subduction interface is key to better anticipate future large earthquakes. Through a modeling of the GNSS velocities on a 3D slab geometry, we were able to obtain useful informations on the location, size, magnitude and return period of future great earthquakes in South Peru. Thereby, we obtained a very heterogeneous spatial distribution of interseismic coupling (degree of locking between the two tectonic plates), with low-coupled areas where the Nazca Ridge and the Nazca Fracture Zone are subducting, but highly-coupled areas close to the coasts of Lima and Arequipa. Finally, we estimate that the South Peru segment between the Nazca Ridge and the Arica band could have the potential to host a $M_w = 8.4$ to $M_w = 9.0$ earthquake, with a one century and one millennial recurrence time respectively.

## 1. Introduction

The size and timing of large earthquakes is a major scientific and societal issue. In subduction zones, the slab interface called megathrust is a gigantic thrust fault that can slip suddenly of many meters generating large $M_w > 7.5$, great $M_w > 8.5$ or giant $M_w > 9.0$ megathrust earthquakes. Among the most powerful earthquakes recorded, we can refer to the 1960 $M_w = 9.5$ Valdivia (Chile) earthquake (Cifuentes, 1989; Kanamori & Cipar, 1974), the 2004 $M_w = 9.2$ Sumatra (Indonesia) earthquake (Chlieh et al., 2007; Stein & Okal, 2005), or the 2011 $M_w = 9.1$ Tohoku (Japan) earthquake (Ide et al., 2011; Koketsu et al., 2011; Simons et al., 2011). The seafloor uplift and subsidence induced by these earthquakes trigger tsunamis that can be a threat to the population that live near coastal areas (Takabatake et al., 2018). It is now widely recognized that the seismic sources of major earthquakes
**Investigation:** B. Lovery, M. Chlieh, J. C. Villegas-Lanza, M. Radiguet, W. Quiroz, A. Socquet

**Methodology:** B. Lovery, M. Chlieh, M. Radiguet, J. M. Nocquet, A. Socquet

**Project Administration:** M. Chlieh, M. Radiguet, N. Cotte, M. Simons, H. Tavera, A. Socquet

**Resources:** M. Chlieh, E. Norabuena, J. C. Villegas-Lanza, N. Cotte, A. Tsapong-Tsague, W. Quiroz, C. Sierra Farfán, M. Simons, J. M. Nocquet, H. Tavera, A. Socquet

**Software:** B. Lovery, M. Radiguet, A. Tsapong-Tsague, J. M. Nocquet

**Supervision:** M. Chlieh, E. Norabuena, J. C. Villegas-Lanza, M. Radiguet, M. Simons, J. M. Nocquet, A. Socquet

**Validation:** B. Lovery, M. Chlieh, M. Radiguet, A. Socquet

**Visualization:** B. Lovery

**Writing – original draft:** B. Lovery, M. Chlieh, M. Radiguet, A. Socquet

**Writing – review & editing:** B. Lovery, M. Chlieh, M. Radiguet, J. M. Nocquet, A. Socquet

often overlap spatially well with locked patches issue from interseismic coupling models where the interseismic coupling is defined as the slip deficit over the long-term slip (Chlieh et al., 2008; Konca et al., 2008; Loveless & Meade, 2011; Métois et al., 2016; Moreno et al., 2010; Radiguet et al., 2016). It remains however challenging to deploy and maintain large and dense geodetic networks to obtain detailed mapping of the interseismic coupling, especially in subduction zone where a significant portion of the seismogenic zone lies below the seafloor and where seismic cycle deformation gradients affect distance of many hundreds of kilometers.

The heterogeneous pattern of the fault locking distributions indicate mechanical variations of the frictional properties of the fault interface. Highly locked fault patches with a coupling degree typically >0.4 are supposed to host the seismic sources of future large earthquakes and are characterized by velocity weakening behavior (Chlieh et al., 2011). These apparent locked zones are often called asperities by abuse of language and may rather reflect a high density of smaller size “real” asperities that lie in the stress shadow of each other making it difficult to dissociate them within an apparently locked patch (Bürgmann et al., 2005). Locked patches appear to be surrounded by creeping patches where the coupling degree is relatively lower. These creeping patches are dominated by velocity strengthening frictional behavior and supposed to act as partial or permanent barriers depending on their frictional and geometrical properties (Chlieh et al., 2008; Hetland & Simons, 2010; Kaneko et al., 2010; King, 1986; King & Nábělek, 1985; Saillard et al., 2017).

Here we present new global navigation satellite system (GNSS) measurements from the Central Andes of Peru issued from a long-term scientific collaboration between Peruvian, French and American institutions. The first GNSS stations were deployed in the aftermath of the 2007 M<sub>w</sub> = 8.0 Pisco earthquake with five continuous cGPS stations around the Paracas Peninsula (Perfettini et al., 2010). Between 2010 and 2012, 41 new continuous stations were installed along the coastline and in the Andes between Lima and Arica. In addition, 33 survey-mode sGPS points were deployed in 2012 and occupied every 2 years in average up to 2022, to densify the network along the coastline south from the Nazca ridge to Arica and in the Altiplano plateau up to the Titicaca Lake. This new data set is dense in the coastal areas and fills many gaps left by previous geodetic studies (Bevis et al., 2001; Chlieh et al., 2011; Nocquet et al., 2014; Norabuena et al., 1998; Villegas-Lanza et al., 2016). One primary objective of this paper is to process and analyze simultaneously these observations to extract a coherent interseismic signal at the scale of the Andes of Central-South Peru.

In this study, we provide a detailed mapping of creeping and locking patches and estimate the average moment deficit rate (MDR in Nm/yr) that is building up on a 3D megathrust interface. The uncertainties of the interseismic coupling models and associated MDR will be quantified including or not a Peruvian sliver (PS) (Villegas-Lanza et al., 2016) and discussed together with historical and instrumental seismic catalogs to better understand the mechanical properties of the slab interface. Finally, we use these results to conduct a moment budget analysis in South Peru to balance the present-day loading rate and the seismicity rate over the long-term (Avouac, 2015; Mariniere et al., 2021). By exploring the related uncertainties, this analysis provides a quantitative estimate of the earthquake potential (M<sub>max</sub>) and associated return period that could occur in this region.

## 2. Seismotectonic Context

Along the Western coast of the South America continent, the Nazca plate is subducting beneath the South American plate with a velocity of about 6 cm/yr (Altamimi et al., 2016; Kendrick et al., 2003). The Peruvian subduction zone is classically segmented into three main segments (Northern, Central, and Southern) of about 700 km each. The Northern segment extends from the Gulf of Guayaquil (3°S) to the Mendaña fracture zone (10°S). The Central segment is delimited by the Mendaña fracture zone (10°S) and the Nazca Ridge (~15°S), while the Southern segment extends from that ridge up to the Arica bend (18.5°S) (Dorbath et al., 1990; Villegas-Lanza et al., 2016). By contrast with the Northern segment that never experienced a great megathrust earthquake in the last 500 years (Nocquet et al., 2014; Villegas-Lanza et al., 2016), the Central and Southern segments are characterized by recurrent great (M<sub>w</sub> ≥ 8.5) earthquakes that occurred in 1746 in the Central and in 1604 and 1868 in the Southern segments. Large (M<sub>w</sub> ≥ 7.5) megathrust earthquakes are even more frequent and occur in average every 20 years during the last five centuries (Bilek, 2010; Chlieh et al., 2011; Dorbath et al., 1990; Nishenko, 1991; Perfettini et al., 2010; Pritchard et al., 2007; Remy et al., 2016; Silgado, 1978; Sladen et al., 2010; Tavera & Buforn, 1998; Tavera et al., 2002; Villegas-Lanza et al., 2016). The new GPS data presented here, recover the whole southern segment and the southern half of the Central segment (Figure 1).
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Seismotectonic background of the Central Andes subduction zone. Light gray ellipses indicate approximate rupture areas of major historical earthquakes since 1512 (Chlieh et al., 2011; Dorbath et al., 1990; Villegas-Lanza et al., 2016). Focal mechanisms from major recent earthquakes are from the gCMT catalog. The approximate rupture areas of the M<sub>w</sub> = 8.1 2014 Iquique and M<sub>w</sub> = 8.4 2001 Arequipa earthquakes have been modeled in this study. M<sub>w</sub> = 8.0 2007 Pisco rupture is reported from Chlieh et al. (2011). Slab isodepth contours are reported every 20 km from Slab2 model (Hayes et al., 2018). The Nazca/Stable South America convergence rate is depicted by the large white arrow (Kendrick et al., 2003). Major faults are from Veloza et al. (2011) and are depicted by blue (normal), green (reverse), and brown lines (strike-slip). Triangles show the location of global navigation satellite system stations used in this study, and are color coded according to their type (cGPS or sGPS) and their operator.</figcaption></figure>

In the North of Chile where the great 1877 M<sub>w</sub> = 8.8 earthquake occurred, the 2014 M<sub>w</sub> = 8.1 Iquique earthquake ruptured a segment between Iquique and Arica (Figure 1). This earthquake has significantly affected the GNSS network of Southern Peru. On the Southern Peru segment, the 1942 M<sub>w</sub> = 8.2 and the 2001 M<sub>w</sub> = 8.4 Arequipa earthquakes are the two largest events since the great 1868 M<sub>w</sub> = 8.8 event, and their cumulative ruptures have broken ~3/4 of the southern segment leaving an unbroken portion of about 150–200 km long between Ilo and Arica that can be considered as a seismic gap (Figure 1). Finally, in the Central segment offshore Lima, four large events with M<sub>w</sub> ~ 8.0 occurred in 1940, 1966, 1974, and 2007 (Pisco) with adjacent ruptures of about 150–200 km each that recover the whole Central segment (Beck & Nishenko, 1990; Beck & Ruff, 1989; Dorbath et al., 1990; Langer & Spence, 1995; Perfettini et al., 2010; Sladen et al., 2010; Villegas-Lanza et al., 2016).

Subandean shortening up to 10 mm/yr in South Peru has been reported by previous geodetic studies (Bevis et al., 2001; Métois et al., 2013; Norabuena et al., 1998; Villegas-Lanza et al., 2016), while paleomagnetic studies reported up to 30 mm/yr of shortening (Arriagada et al., 2008). In the last decade, rigid sliver motion has been introduced to explain the deformation pattern in Central Andes (Nocquet et al., 2014; Villegas-Lanza et al., 2016). Nocquet et al. (2014) defined a so-called Inca Sliver including southern Ecuador, Peru, and Bolivia extending from the trench to the Subandean belt. From an analysis of their GNSS velocities, they computed a Euler pole at 63.76°W, 22.47°N with a rate of 0.092°/Myr relative to the Stable South America (SSA) frame (SSA). While Villegas-Lanza et al. (2016) proposed a so-called PS associated to a Euler pole at 67.23°W, 8.36°N with a rate of 0.104°/Myr, from the residuals of their inversions of interseismic velocities. This PS encompass the whole GNSS network used in our study (12°–18°S), however it would be questionable to extend it further to the South, as the deformation pattern in Chile south of 19°S shows a different motion. Paleomagnetic studies
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Global navigation satellite system network used to extract the interseismic velocity field. Stations are color coded according to their type (cGPS or sGPS) and their operator.</figcaption></figure>

reported in Arriagada et al. (2008) indicate an anti-clockwise rotation in the long-term displacement vectors north of 19°S, and a clockwise rotation south of 19°S. This pattern is also reported in geodetic studies (Chlieh et al., 2011; Métois et al., 2016). Yáñez-Cuadra et al. (2022) proposed a different approach, with diffuse continental deformation and mantle wedge viscoelastic relaxation, leading to minimal rigid motion needed to explain the surface displacements in the Atacama region (23°–30°S). In this study, we remain on a rigid block motion approach, with purely elastic rheology, as the analysis of intracontinental deformation is beyond the scope of this paper, although we acknowledge that viscoelasticity is a physical process that plays a key role in the deformation of subduction zones. Our results are compared in Section 5.3 to the data-model misfits expected from a fully relaxed viscoelastic model (Li et al., 2015).

# 3. Data Analysis and Trajectory Models

## 3.1. Permanent and Campaign GNSS Networks

We analyze here the crustal motion recorded by 73 GNSS sites (Figure 2) deployed between 1993 and 2022 in the Andes of Central and Southern Peru between latitudes 12°–18°S and longitudes 60°–78°W. We also used continuous GNSS time series in South Peru and North Chile in order to assess the contribution of co and post-seismic deformation associated with recent earthquakes, notably the 2014 $M_w = 8.1$ Iquique earthquake in North Chile (~20°S) and the 2001 $M_w = 8.4$ Arequipa earthquake (Figure S20 in Supporting Information S1) in South Peru (~17°S) which have significantly affected the southern Peruvian stations. In Peru, this GNSS network is composed of:

* Twenty-one continuous GNSS (cGPS) stations deployed in the frame of a collaboration between the Instituto Geofísico del Perú (IGP), the Institut des Sciences de la Terre (ISTerre), and the California Institute of Technology (Caltech). They are maintained since 2015 by the IGP. Maximum data coverage is between 2007 and 2022, but several stations do not provide data for such a long time-range because of operating issues. Most of the stations are located between latitudes 16° and 18°S.

* Twenty-one continuous GNSS stations deployed by the Instituto Geográfico Nacional (IGN) of Peru covering the 2010–2018 time period. These stations are spread on a large area between latitudes 12° and 18°S.

* One continuous IGS (Johnston et al., 2017) GNSS station installed in 1993, the AREQ station located close to Arequipa city.
* Twenty-six survey GNSS points (sGPS) deployed by IGP and ISTerre and measured on average every 2 years, when possible, between 2012 and 2022. Campaigns have been performed in June 2012, June 2013, October 2013, April 2014, October 2015, November 2016, June 2018, November 2018, and March 2022.

* Four continuous GNSS time series from the Centro Sismológico Nacional (CSN) of Chile.

Finally, to correct properly the GNSS campaign measurements from the displacements generated by the 2014 $M_w = 8.1$ Iquique earthquake, we used 21 CSN cGPS stations located close to the rupture area (Figure S2 in Supporting Information S1).

All the GNSS time series with their trajectory model are available in Supporting Information S1.

## 3.2. Data Processing

Data were processed following a Precise Point Positioning approach at ISTerre using the GipsyX software v1.5 from the Jet Propulsion Laboratory (JPL) (Bertiger et al., 2020). To resolve the phase ambiguity, GipsyX imposes wide-lane phase bias constraints provided by JPL. Orbit and clock estimations are taken from JPL final products. NNR orbits are used. Ionospheric delays are canceled by using dual-frequency communications. No second order ionospheric correction is applied. The VMF1 (Boehm et al., 2006) is used to estimate the tropospheric zenith delay parameters. Tropospheric delays are inverted from the data: a delay is inverted every 5 min and two horizontal gradients are accounted for per session of 24 hr. We use the FES 2014b model to correct from ocean tide loading (Lyard et al., 2021). Phase center stability and multi-path protection are achieved using antenna calibration tables and by stacking GNSS measurements over the day to keep one position per day. In the case of GNSS campaigns, the positioning of the antenna precisely above the campaign point is usually an issue, that is minimized using a forced-centering system at most points. Scaling, rotation and translation factors provided by JPL are then applied in the Helmert transformation to map the solution in the IGS2014 (Altamimi et al., 2016). The processed GNSS time series indicate the daily positions of the station covering a time window of several years, as shown in Figure 3 and time-series in supplements. The displacements extracted from our GNSS data are given in the ITRF2014 reference frame (Altamimi et al., 2016). To express them relatively to the SSA reference frame, we used the ITRF2014/SSA Euler pole located at 18.68°S and 128.69°W with an angular velocity of 0.122°Myr⁻¹ (Altamimi et al., 2016). We double-checked the stability of our reference frame by processing various GNSS sites located in Brazil, and expressed in the SSA frame (Figure S14 in Supporting Information S1).

## 3.3. Trajectory Modeling of the Continuous GNSS Time Series

To model the GNSS time series, we applied a trajectory model that mimics the position of a GNSS station as a function of time using the ITSA software (Marill et al., 2021). The trajectory model used here is the sum of sub-models that sketch individually processes such as the secular motion of the station, instantaneous jumps either due to earthquakes or antenna changes, post-seismic relaxation and seasonal oscillations (Bevis & Brown, 2014; Marill et al., 2021). A schematic view of the process is given in Figure S1 of Supporting Information S1. Because no obvious signature of slow slip events has been detected in our processed GNSS time series, no term is assigned to that process here. Mathematically, the evolution of the station position $x(t)$ is expressed as the summation of linear, Heaviside, logarithmic and sinusoidal time functions that are described as follows:

$$x(t) = x_R + v(t - t_R) + \sum_{j=1}^{n_j} b_j H(t - t_j) + \sum_{l=1}^{n_j} a_l \log_{10} \left( 1 + \frac{t - t_l}{\tau_{Ra}} \right) + \sum_{k=1}^{2} [s_k \sin(2k\pi(t - t_R)) + c_k \cos(2k\pi(t - t_R))] + \sum_{a=1}^{n_A} b_a H(t - t_a) \quad (1)$$

* The first term $x_R$ being the station position at $t = 0$.

* The second term $v(t - t_R)$ is a linear term corresponding to the secular velocity, including a steady-state interseismic velocity with $t_R$ being the reference time.

* The third term corresponds to coseismic displacements modeled by a Heaviside function where the coefficient $b_j$ characterizes the coseismic jump at time $t_j$. Major seismic events ($M_w > 6.4$) reported in the ISC are considered, if they are in a radius of influence (in km) defined as $r(M_w) = 10^{(0.43Mw-0.7)/1.15}$, $M_w$ being the moment magnitude of the considered event.
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Time series on the norm of the horizontal components of continuous global navigation satellite system stations along the coast (a), relative to Stable South America, and corrected from antenna jumps and seasonal variations: LAGU, LOMI, SJUA, ATIC, CHRA, PMCA, and LYAR. Data points are displayed in blue, while the trajectory model is in red. The time series are arranged from North (top) to South (bottom). Coseismic and postseismic displacements from the $M_w = 8.1 2014$ Iquique earthquake at LYAR, PMCA and SJUA stations are clearly noticeable. At the LAGU station, we observe the postseismic relaxation of the $M_w = 8.0 2007$ Pisco earthquake. On the right panel (b), are displayed the station locations and the earthquakes accommodated in the time series.</figcaption></figure>

* The fourth term corresponds to the postseismic deformation detectable after major seismic events, that is classically modeled by a decaying logarithmic function with a characteristic time $\tau_{Ra}$ and a displacement $a_i$ at time $t_i$. A value of $\tau_{Ra} = 30$ days was found adequate for our data.

* The fifth term corresponds to annual and semi-annual seasonal oscillations, due to hydrological loading and temperature variations weighted by the coefficients $s_k$ and $c_k$.

* Finally, the last term is related to antenna changes leading to artificial jumps $b_a$ at time $t_a$ in the time series if the characteristics of the new antenna are different. These jumps are especially visible on the vertical component of the time series and are corrected from field work records.

To better appreciate the variability of the signals recorded by permanent stations, Figure 3 reports the norm of the horizontal components of four GNSS time series with their respective trajectory model, and the large earthquakes that occurred in that region during the time of observation (2007–2022). All the stations shown are located along the coastline between the Paracas Peninsula and Arica. These time series are all expressed in the SSA reference frame. Station LAGU, located near the Paracas Peninsula and installed right after the 2007 $M_w = 8.0$ Pisco earthquake shows the postseismic relaxation and reloading process that have followed. Station LYAR installed in South Peru has recorded the 2014 $M_w = 8.1$ Iquique earthquake and its postseismic relaxation. Finally, station SJUA located between LAGU and LYAR was sensitive to the 2013 and 2018 $M_w = 7.1$ Acari earthquakes.

It is worth noting that this trajectory model software is only relevant when the time series has enough data points (> 100) and an appropriate temporal distribution (at least 1 year of data). Although it is not a binding criterion for the data extracted from our cGPS stations, it prevents the use of the trajectory model for our sGPS time series. Consequently, we use another methodology to correct our sGPS time series.
### 3.4. Corrections of Survey-Mode GNSS Velocities

In the purpose of correcting properly the coseismic jumps and afterslip offsets generated by the 2014 Iquique seismic sequence at each of the survey-mode sGPS station, we estimate the displacements associated with the 2014 $M_w$ = 8.1 Iquique mainshock on 1 April, its main $M_w$ = 7.7 aftershock on 3 April, and the 30 first days of postseismic afterslip at 14 continuous cGPS stations in South Peru and 21 cGPS stations in north Chile (Figure S2 in Supporting Information S1). The coseismic displacements tables associated with each of these earthquakes as well as the afterslip motion after 30 days are listed in Table S2 of Supporting Information S1. Then, using the inversion procedure described below (Section 4.1), we performed static slip inversions of these displacements constrained by the seismic moment and rake vector from the gCMT. The derived slip sources are then used to predict through forward modeling the coseismic jumps to be corrected at each of the survey-mode GNSS sites in South Peru. Figure S2 in Supporting Information S1 shows the coseismic and afterslip recorded displacements by the permanent GNSS stations in South Peru used for the slip inversions, and the predicted displacements that have been used to correct survey GNSS points (Table S3 in Supporting Information S1). We retained a regularization length $\lambda$ = 10 km and a damping parameter $\sigma_m(k) = 10^{-1.2}$, as it provides a good balance between a low misfit and a reasonable consideration of inter-stations distances (see Section 4.1 for definition of the parameters). A 250° azimuth from North (clockwise) was selected, as it is consistent with the slip vector from gCMT and the plate convergence direction. The modeled slip distribution is in quite good agreement with Jara et al. (2018).

In addition to the modeling of the Iquique earthquake coseismic and postseismic displacements, we modeled the seasonal displacements at each sGPS point based on the approach of Hoffmann et al. (2018). We extract the seasonal signal modeled at each cGPS stations by the ITSA trajectory model, as depicted on Figure S15 in Supporting Information S1. Then, we spatially interpolate this seasonal signal at the sGPS locations, and we subtract it from the sGPS time series. Since the GNSS campaigns are often performed at the same season 1 year from another, the seasonal correction is not very strong. Nevertheless, we observed a reduction of 2% in the standard deviation when correcting the sGPS time series for the seasonal, as well as a velocity change of up to 5% at some sGPS points.

### 3.5. Extracting the Interseismic Velocity Field From cGPS and sGPS Analysis

Using our trajectory model for cGPS stations by keeping only the linear term, and performing a linear regression in our corrected sGPS position time series, we finally extracted the interseismic velocities on the horizontal components (East, North) for each of the 73 stations. We chose not to consider the vertical component for the inversions. Indeed, the errors on the vertical component is much larger, especially on sGPS sites. In addition, the vertical deformation is very sensitive to non-tectonic processes, such as hydrological processes, volcanic activity, or local terrain subsidence (Itoh et al., 2019). The observed vertical interseismic velocities at cGPS sites are reported on Figure S21 in Supporting Information S1. We observe that they are very heterogeneous, even at a relatively local scale, which would be an issue for the modeling. All these velocities are listed in a SSA reference frame in Table S1 of Supporting Information S1 and shown in Figure 4. Uncertainties on the velocities have been estimated through a statistical analysis of the time series using Median Interannual Difference Adjusted for Skewness (MIDAS) described in Blewitt et al. (2016). In a nutshell, MIDAS is a trend estimator method working with pairs of data, assuming a majority of the data have a Gaussian probability distribution function, adopting an unbiased Theil-Sen method. In SSA reference frame, all the interseismic GNSS velocities are oriented toward the Cordillera, coherently with the relative Nazca/South America convergence direction. Velocity gradients decrease with increasing distance from the trench axis consistently with gradients expected from interseismic locking models. Trench-lateral variations in the amplitudes of the velocities suggest lateral variations in the locking distribution (See cross section Figure 4).

It is worth noting that the AREQ station, located in Arequipa and recording since 1993, is still affected by the effects of the 2001 Arequipa earthquake. These effects are discussed in Section 5.2 and may affect other stations in the close area.

## 4. GNSS-Derived Interseismic Slip Inversion

### 4.1. Slip Inversion Procedure

In order to obtain an interseismic coupling model on the megathrust interface, we first build a 3D slab geometry based on the South American slab contours of the Slab2 model between latitudes 10° and 25°S (Hayes
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Interseismic velocity field relative to the Stable South America (SSA) frame. Velocities on continuous global navigation satellite system (GNSS) stations (cGPS) are displayed in red, while velocities on survey GNSS points (sGPS) are displayed in blue. Velocity uncertainties from Median Interannual Difference Adjusted for Skewness (Blewitt et al., 2016) are depicted by ellipses. Major faults are from Veloza et al. (2011). Four cross-sections (a–d) are plotted perpendicular to the trench, with 100-km wide tracks. The norm of the velocity at each station is displayed in red (cGPS) and blue (sGPS), the model shown in Figure 5 (SSA, with $\lambda$ = 30 km) being depicted by the gray line. The topography is plotted above.</figcaption></figure>

et al., 2018). Our modeled slab interface shown in Figure S3 of Supporting Information S1 is discretized into 6,700 15-km wide elementary triangular dislocations, with source points placed at the center.

The inversion procedure follows the formulation of Tarantola and Valette (Radiguet et al., 2011; Tarantola, 2005; Tarantola & Valette, 1982), implemented within the static inversion of the Principal Component Analysis Inversion Model software (Kositsky & Avouac, 2010). To model the interseismic coupling, we used a backslip approach (Savage, 1983) on the 3D slab geometry embedded in a purely elastic medium with a shear modulus $\mu$ of 40 GPa. We define interseismic coupling as the fraction between the slip rate on a patch and the plate convergence velocity, the interseismic coupling is therefore ranging from zero for a fully creeping patch to one for a fully locked patch. The shear modulus (a.k.a. rigidity) is difficult to estimate accurately, as it depends on the location, and more generally it increases with depth (Bilek & Lay, 1999). The retained value of 40 GPa is a compromise of the various values found in the literature (Bilek & Lay, 1999; Villegas-Lanza et al., 2016), considering the average depth of our coupling distribution. The slip amplitude and rake are fixed consistently with the horizontal projection of the relative plate convergence vector $V_{pl}$. The slip on the megathrust interface can be decomposed into a strike-slip and a dip-slip component for each sub-fault patch. The matrix $d$ of surface displacements can then be obtained by the following formula:

$$d = G \cdot m \qquad (2)$$

where $m$ is the matrix of slip rate on the fault at depth, and $G$ the associated Green's function. Green's functions are computed based on the analytical equations from Meade (2007). A linear least-squares inversion is used to minimize the cost function $S(m)$ defined as:

$$S(m) = \frac{1}{2} \left[ (Gm - d)^{\prime} C_{d}^{-1} (Gm - d) + (m - m_{0})^{\prime} C_{m}^{-1} (m - m_{0}) \right] \qquad (3)$$

and the modeled slip rates $m$ on each source element:

$$m = m_{0} + C_{m} G^{\prime} (GC_{m} G^{\prime} + C_{d})^{-1} (d - Gm_{0}) \qquad (4)$$
<figure><img src="image_8.png" /><figcaption>Figure 5. Interseismic coupling distribution from our inversion of interseismic velocities in the Stable South America frame, with various smoothing (λ = 15, 30, and 200 km). Slab contours (dashed lines) are reported every 20-km depth, from the Slab2 model (Hayes et al., 2018). The interseismic coupling is highly heterogeneous reflecting strong variations of the frictional properties on the slab interface. Thus, there are four highly coupled areas: a very-highly coupled area close to Lima, a less coupled area near Ica, a large band of high-coupling between the Nazca Ridge and the Nazca Fracture Zone (NFZ), and a highly coupled area between the NFZ and the Arica bend. GNSS sites are depicted by gray triangles. Discontinuities can be observed where the Nazca Ridge and the NFZ are subducting below the South America plate. The approximate rupture area of the M<sub>w</sub> = 8.4 2001 Arequipa has been estimated in this study, while the M<sub>w</sub> = 8.0 2007 Pisco earthquake rupture is reported from Chlieh et al. (2011). Patches with gray hachure pattern depict areas that are not well resolved by the model. Four cross-sections (a–d) are plotted perpendicular to the trench, with 100-km wide tracks. The norm of the velocity at each station is displayed in red (cGPS) and blue (sGPS), the three models (λ = 15, 30, and 200 km) are plotted with gray lines. The topography and slab profile are plotted above. Interseismic coupling from the model in panel 2 (λ = 30 km) is indicated by the color scale on the slab interface.</figcaption></figure>

$m_0$ being the prior model, $C_d$ and $C_m$ the co-variance matrices for data and model, $G'$ the transpose matrix of $G$. The data co-variance matrix $C_d$ is defined as the diagonal matrix of the associated uncertainties on the measurements. We neglect the co-variances, so all terms outside the diagonal become zero. The model co-variance matrix $C_m$ is used to regularize the slip distribution between neighboring elementary patches which is achieved by an exponential smoothing function expressed as:

$$C_m(k, l) = \left( \sigma_m \frac{\lambda_0}{\lambda} \right)^2 e^{-\frac{d(k,l)}{\lambda}} \quad (5)$$

$d(k,l)$ being the distance between source elements $k$ and $l$, $\lambda_0$ a scaling factor (characteristic size of fault patch), $\lambda$ the spatial correlation length, and $\sigma_m$ the damping parameter controlling the relative weight of the distance to initial model versus data misfit. In practice, we adjust the regularization length $\lambda$ and $\sigma_m$ in order to maintain an interseismic coupling between zero and one (cf., Section 4.2). Finally, we impose the coupling to taper at 60 km depth of the slab interface, as the thermodynamics conditions deeper should prevent fault locking (Gagnon et al., 2005).

To quantify the misfit between the observation and the model, we define the reduced chi-square criteria as:

$$\chi^2_{\text{red}} = \frac{\chi^2}{\nu} \quad \text{with} \quad \chi^2 = \sum_i \frac{(d_i - m_i)^2}{\sigma_i^2} \quad (6)$$

$d$ being the data, $m$ being the model, $\sigma^2$ the variance, and $\nu = n - p$ the degrees of freedom for $n$ observational data and $p$ fitted parameters.
In order to determine which areas are modeled with confidence, we can define the resolution matrix R:

$$R = C_m G^l (G C_m G^l + C_d)^{-1} G \quad (7)$$

where $C_m$ and $C_d$ are the co-variance matrices for model and data parameters and G the matrix of Green's functions. The closer the resolution matrix R is to the identity matrix $I_d$, the better the model is resolved (Radiguet et al., 2011). Then, we define the restitution as the sum of the rows of the resolution matrix R. Both the restitution and the diagonal elements of the resolution matrix are plotted on Figure S4 in Supporting Information S1. Whereas the diagonal elements of R inform us about how much the slip of a specific patch is accurately mapped, the restitution informs us if this slip has been correctly projected onto other patches.

Considering the 3D slab geometry and the location of the GNSS sites presented above, we found that the diagonal values of the resolution matrix R are relatively low due to the high number of small patches (6,700) for a relatively low amount of GNSS sites (73), in order not to impose too much constraint on the slip location. The diagonal values are higher where the density of station is high, indicating a better resolution there. The same goes for the restitution, except that the distribution is smoother. Overall, we can state that the spatial resolution of slip in the seismogenic zone above 60-km depth of the slab interface is relatively high between latitudes 12° and 18°S.

## 4.2. Results of Interseismic Coupling Models

Following the inversion procedure described above and considering a creeping slip prior model (i.e., $m_0 = 0$), we tested different values of regularization length and damping parameter in order to find the optimum GNSS fitting parameters. Figures S5 and S6 in Supporting Information S1 displayed L-curves that compare the $\chi^2$ for variable parameters. We searched for the parameter values which offer an adequate compromise between the misfit and model complexity (Florsch et al., 2014; Radiguet et al., 2011). In other words, we seek to minimize the misfit, while keeping a regularization length in the same order as the average distance between the GNSS sites and a maximal slip bounded by the plate the plate convergence rate.

Rather than showing the so-called best model, we propose a family of acceptable solutions for regularization lengths $\lambda$ of 15, 30, or 200 km that all provide an acceptable GNSS misfit with $\chi^2$ ranging from 0.28 to 0.47. Figure 5 reports the coupling distributions for various smoothing factors $\lambda = 15, 30, and 200$ km. The rough solutions may better reflect the coupling in the Southern segment where the slip resolution and restitution are much better. Low coupling is shown close to the trench, however this area is not very well resolved by the model (cf., Figure S4 in Supporting Information S1). Consequently, we also run inversions with a locked slip prior (Figure S16 in Supporting Information S1), showing high coupling close to the trench. In all these models, it appears that the coupling distributions tend to extend up to 60-km depth where it is supposed to creep due to the high temperatures there (>400°C according to Villegas-Lanza et al. (2016)). This may suggest that an element is missing in our model, like a sliver motion proposed by Nocquet et al. (2014) and Villegas-Lanza et al. (2016), that we are going to discuss in Section 5.1.

We define the MDR ($dMo_{inter}/dt$) of a coupling model as:

$$\frac{dMo_{inter}}{dt} = \int \mu(V_{pl} - V_{back}) dS \quad (8)$$

the integration being computed over the surface S of the slab interface, and where $(V_{pl} - V_{back})$ is the back-slip offset rate across the fault. In practice, we compute for each discretized patch k its local moment deficit that we sum up over a specific surface. For instance, if we integrate over the whole megathrust up to 60-km depth from Lima to Arica, we find moment deficit rates ranging from 2.76 to $3.58 \times 10^{20}$ Nm/yr respectively for models with $\lambda = 15$ and 200 km.

## 4.3. Extension of the Peruvian Sliver to South Peru?

Figure 6 and Figure S18 in Supporting Information S1 shows the GNSS residuals (data minus model) associated with the coupling distribution shown in Figure 5-2, inland GNSS residuals located in the Cordillera between 300 and 500 km from the trench indicate a southeastward motion of about 4 mm/yr (Figure 6a). This coherent motion is observable along ~700 km at all GNSS sites located in the Cordillera west of Lima up to Lake Titicaca. This
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Residuals from our inversions using global navigation satellite system (GNSS) data in the Stable South America (SSA) frame (a) and in the Peruvian Sliver (PS) frame (b), for a regularization length $\lambda = 30$ km. Most of the residuals are below 4 mm/yr, although some have higher values. Survey GNSS points are in blue, while continuous GNSS stations are in red. Uncertainties computed by Median Interannual Difference Adjusted for Skewness (Blewitt et al., 2016) are depicted by gray ellipses. Residuals at inland stations, within the black dashed rectangle, are highlighted in green on the roses. In a SSA frame, we see a southwestward direction that dominates in the residuals, especially well seen in the GNSS sites far away from the trench. In the PS reference frame, residuals are distributed randomly and might be due to internal deformation within the Andes or long-term viscous effect associated with recent large earthquakes that we do not account for here (Perfettini et al., 2005).</figcaption></figure>

tendency disappears when using the PS reference frame proposed by Villegas-Lanza et al. (2016) (Figure 6b), suggesting that there is room for a rigid sliver moving about 4 mm/yr southwestward, within which all our stations are included (Figure S13 in Supporting Information S1). Considering that PS as a reference frame, the relative Nazca/PS convergence rate would be reduced to ~55 mm/yr.

When the velocities rotated in the Peruvian reference frame are inverted to obtain the coupling, the orientation of the residual velocities is randomly distributed along the coastline.

If we consider the interseismic GNSS velocities in that PS reference and perform similar inversions to those described above, the data misfit is improved compared to previous models without sliver ($\chi^2 = 0.24$ with the sliver, while $\chi^2 = 0.30$ with the interseismic velocities in SSA reference). In addition, the new coupling distributions highlight much better the variations of the downdip limit of the locked subduction zone (compare Figures 5 and 8, or refer to Figure S19 in Supporting Information S1). The fact that these inversions with the sliver naturally contain the slip in the first 50–60 km of depth is a strong physical argument in favor of this sliver motion, compared to the model in the SSA frame which required a deeper coupling to fit the data. Increasing the smoothing factor tends to extend the coupling near the trench. The resulting global moment deficit rates are about 25% lower than for the previous models without sliver.

# 5. Discussion

## 5.1. Comparison With Previously Published Interseismic Coupling Models

Norabuena et al. (1998) and Bevis et al. (2001) have published preliminary insights on interseismic coupling in the area. Norabuena et al. (1998) obtained a good fit to the data by considering a locking of about 50% up to 20-km depth and 12 mm/yr of shortening, while Bevis et al. (2001) retained a locking close to 100% from 10 to 50-km depth and a 5–6 mm/yr shortening. However, both studies rely on only one continuous station (AREQ) and about 20 survey points between Lima and Arica with only two years of data (1994 and 1996 SNAPP measurements). In addition, no detailed slab geometry was available at the time, and they do not account for the seismic events or seasonal variations that could have affected their measurements. A decade later, new studies from Chlieh et al. (2011) and Villegas-Lanza et al. (2016) were able to present more detailed interseismic coupling distributions thanks to new developments.
Compared to Villegas-Lanza et al. (2016) and Chlieh et al. (2011), our results confirm some features while refin-
ing their restitution. The best fitting models from the three studies show a highly coupled area close to the coast
of Lima, one close to Pisco, one south of the Nazca Ridge, and one between the Nazca Fracture Zone (NFZ) and
the Arica bend. These four highly coupled areas are separated by three low coupled areas: one between Lima and
Pisco, and two where the Nazca Ridge and the NFZ are subducting.

In our study, the decoupling between Lima and Pisco is centered at about 13.5°S, while in Villegas-Lanza
et al. (2016) and Chlieh et al. (2011) this decoupling is observed at about 12.5°S. However, this area is not well
resolved, therefore the apparent decoupling may in fact be related to a lack of data and should be considered
with caution. Regarding the creeping area where the Nazca Ridge is subducting, it is a persistent feature in all
the models, with a width of at least 70 km and even 150–200 km in Chlieh et al. (2011). It is difficult to assess
accurately the width of this creeping area, as the GNSS sites are at least 50 km apart in this area.

It is on the Southern segment that our study provides the most clarification, with a significant increase in station density. There were seven stations between the Nazca Ridge and the Arica bend in Chlieh et al. (2011), eight in Villegas-Lanza et al. (2016), and 57 in our study. This improved resolution highlights a high coupling between the Nazca Ridge and the Arica bend, lowered where the NFZ is subducting but not as much as in Villegas-Lanza et al. (2016) or Chlieh et al. (2011). Having a lot of stations both on the coast and inland, up to Lake Titicaca, is a substantial improvement to constrain the depth of coupling compared to older studies. In addition, when accounting for a PS as defined by Villegas-Lanza et al. (2016), the maximum coupling depth is 10 km shallower. Apart from the number of stations, the type of station and observed time range also differ between the three studies. In South Peru, Chlieh et al. (2011) essentially based their analysis on survey points measured between 1993 and 2003, while Villegas-Lanza et al. (2016) used mostly continuous stations measured between 2007 and 2013. In our study, we analyzed a mix of continuous (47) and survey (26) points, measured in the 2007–2022 and 2012–2022 periods respectively. The extensive use of continuous measurements and longer time range reduce uncertainties compared to older studies. It is worth noting that the interseismic velocities from Chlieh et al. (2011) are representative of an earlier time period, mostly before the 2001 Arequipa earthquake. Our interseismic velocities at cGPS sites are extracted using a more comprehensive trajectory model, by accounting for various coseismic, postseismic, or seasonal displacements. While for sGPS sites we account for seasonal variations, as well as coseismic and postseismic displacements induced by the 2014 Iquique earthquake. Finally, our inversions are performed using a 3-D slab model (Hayes et al., 2018).

5.2. Long-Living Effects of Post-Seismic Viscous Relaxation

Some of the far-field residuals located on the Altiplano and Western Cordillera are pointing toward the trench (in the PS frame) and could be associated either with long-term postseismic viscous effects of the 2001 Arequipa and 2007 Pisco earthquakes (Hergert & Heidbach, 2006; Klein et al., 2016; Remy et al., 2016), or with internal deformation within the Andes (Kley & Monaldi, 1998). The Arequipa (AREQ) continuous station, installed in 1993 and operating until now, was the only one active during the 2001 earthquake in Southern Peru. We remark that the reversal time of motion of the AREQ station is about 7 years after the 2001 earthquake (Figure 7) suggesting that during that period the record of that station was dominated by postseismic (afterslip and viscous) relaxation. After 2008, the interseismic loading process became the most dominant process and the velocity increased gradually to nearly reach its pre-2001 velocity at the beginning of 2020. We quantify the velocity change at the AREQ station between the period prior to the 2001 event (1997–2001) and our period of observation on the GNSS network in South Peru (2012–2023 for most stations), by computing the average interseismic velocity (module of the North and East components) on the two periods at the station. We found that during our period of observations (2012–2023), the interseismic velocity measured at the Arequipa (AREQ) station is reduced by about 15% in the SSA frame, and 30% in the PS frame, compared to the velocity measured prior to 2001, suggesting that our estimates of MDR may be a lower bound in the region surrounding the 2001 rupture. However, the AREQ station is located close to the earthquake centroid (less than 200 km), and therefore probably more affected than most of the stations in the South Peru segment. These estimates of 15% (SSA) and 30% (PS) should therefore be seen as a maximum post-seismic effect on the MDR.
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Time series for the AREQ station, located in the city of Arequipa on the horizontal components (North, East), corrected from antenna jumps and seasonal variations. Data points are displayed in blue, while the trajectory model is displayed in red. North and East velocities are reported for two time ranges: 1997–2001 and 2012–2022, in order to compare the interseismic velocity before the 2001 event to the velocity in recent years.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 8.** Interseismic coupling distribution from our inversion of interseismic velocities in the Peruvian Sliver (PS) frame, with three different smoothing ($\lambda$ of 15, 30, and 200 km). Slab contours (dashed lines) are reported every 20-km depth, from the Slab2 model (Hayes et al., 2018). GNSS sites are depicted by gray triangles. The interseismic coupling is highly heterogeneous reflecting strong variations of the frictional properties of the slab interface. Thus, there are four highly-coupled areas: a very-highly coupled area close to Lima, a less highly-coupled area near Ica, a large band of high-coupling between the Nazca Ridge and the Nazca Fracture Zone (NFZ), and a highly-coupled area between the NFZ and the Arica bend. Discontinuities can be observed where the Nazca Ridge and the NFZ are subducting below the South America plate. The approximate rupture area of the $M_w = 8.4$ 2001 Arequipa has been estimated in this study, while the $M_w = 8.0$ 2007 Pisco earthquake rupture is reported from Chlieh et al. (2011). On the right panel are displayed the moment deficit rate by along-trench distance for models in the PS and Stable South America frame, and with various smoothing ($\lambda$ of 15, 30 and 200 km), computed over 10-km wide segments. Patches with gray hachure pattern depict areas that are not well resolved by the model. Rupture extents from major seismic events are reported as blue segments.</figcaption></figure>
### 5.3. Possible Effect of Viscoelastic Rheology on the Interseismic Deformation Field

In order to estimate the effect of viscoelastic rheology on the interseismic deformation field, Li et al. (2015) computed elastic and viscoelastic interseismic models on the North Chile subduction zone. They applied a uniform locking from the trench to a given maximum locking depth, varying from 30 to 80 km, and they quantified the misfit on the horizontal GNSS observations. From this approach, they observed an improved minimum of misfit when applying the viscoelastic model, especially when considering the backarc sites only. In addition, the locking depth corresponding to the minimum of misfit is lowered with the viscoelastic model: from ~55 to ~45 km with all the data, and from ~55 to ~50 km with backarc data only.

In a similar way, we tested models with uniform fully locked megathrust interface from the trench to various depths varying from 20 to 80 km. Then, we quantified the evolution of the reduced $\chi^2$ as function of the maximum locking depth. Results are reported on Figure S17 in Supporting Information S1, when considering all the GNSS sites (left) or only the backarc GNSS sites (right) indicated on Figure 6. We consider either SSA (in blue) or the PS (in orange) reference frames. We observe on both cases that the minimum value of $\chi^2$ is significantly reduced when considering the PS frame relative to the SSA frame: from 0.82 to 0.45 with all the GNSS sites, and from 0.32 to 0.18 when considering the backarc GNSS sites only. Furthermore, the maximum locking depth corresponding to this minimum is also reduced when considering the PS frame: from 45 to 40 km with all the GNSS data, and from 52 to 35 km with the backarc points only. This behavior is very similar to the one observed by Li et al. (2015) when comparing viscoelastic and purely elastic models. Consequently, considering a rigid PS block motion shows relatively similar geodetic fit.

Performing a full viscoelastic model of the loading on the subduction is beyond the scope of this paper. We therefore chose to model internal deformation by a sliver and estimate the variability in the models to assess the seismic potential and its uncertainties (cf., next section). A viscoelastic model that accounts for viscoelastic deformation, including the relaxation following the 2001 Arequipa and the 2007 Pisco earthquake, should ideally be performed to move toward a better mechanical interpretation of the deformation field in the central Andes that account for the whole seismic cycle (e.g., Li & Chen, 2023).

### 5.4. Along-Strike Variations of the MDR and Seismotectonic Parameters Controlling the Seismic Segmentation

We compute between Lima and Arica the along trench-strike variations of the MDR (Figure 8) by summing up the MDR of all nodes down to 60-km depth within 15-km wide strips normal to the trench. Although the final coupling distributions indicate noticeable differences, the MDR integrated at the scale of the studied area remains in a limited range of $2.0–2.9 \times 10^{20}$ Nm/yr (~$M_w$ = 7.4–7.6), respectively for $\lambda$ varying between 15 and 200 km (Table 1 and Table S4 in Supporting Information S1). The MDR increases with $\lambda$, however this increase is relatively low in the South Peru segment where it evolves from 1.4 to $1.8 \times 10^{20}$ Nm/yr, compared to the Central Peru segment where the MDR dispersion between models is much higher (0.6–1.1 $\times 10^{20}$ Nm/yr). On Figure S16 in Supporting Information S1 are displayed the models with a locked slip prior, the corresponding MDR values are reported in Table 1. Overall, a locked slip prior increases the MDR of ~80% when integrating from Lima to Arica, however this raise is limited to ~20% when integrating from Nazca to Arica, the later area being better constrained by the data. The moment budget analysis described in Section 5.5 will be carried out on the Nazca-Arica segment.

#### 5.4.1. Control of Lateral Variation of Coupling: Entering in Subduction of Seafloor Geomorphological Features

The dense GNSS network used in this study highlights creeping areas at the locations of the Nazca Ridge and the NFZ characterized by a rate-strengthening friction. These creeping areas are well resolved by the density of GNSS sites and appear to be persistent features of our coupling models. Because of its relatively high along-strike width of ~60 km and its historical behavior to stop systematically large and great earthquakes, the Nazca Ridge segment is considered a strong barrier. By contrast, the NFZ segment, which is ~20-km long is often crossed by the rupture of large events and rarely appears to have acted as a barrier in the last five centuries of seismic history. In addition, the decrease of coupling where the NFZ is subducting may be overestimated due to the visco-elastic relaxation following the 2001 Arequipa earthquake. Indeed, the location of this creeping intersegment falls in the middle of the 2001 rupture area, it is therefore difficult to discriminate between what may be related to the visco-elastic relaxation and what may be actually creeping. Robinson et al. (2006) have shown that this fracture
<table><caption>**Table 1**
Misfit and Moment Deficit Over the Whole Resolved Area (Lima to Arica) or Restrained to the South Peru Segment (Nazca to Arica) for Various Models</caption><thead><tr><th rowspan="2">Reference frame</th><th rowspan="2">Slip prior</th><th rowspan="2">λ (km)</th><th rowspan="2">log<sub>10</sub>(σ<sub>m</sub>)</th><th colspan="2">Moment (Nm/yr)</th><th rowspan="2">Displayed on:</th></tr><tr><th>Lima-Arica</th><th>Nazca-Arica</th></tr></thead><tbody><tr><td>SSA</td><td>Creeping</td><td>15</td><td>-1.53</td><td>2.76 × 10<sup>20</sup></td><td>2.04 × 10<sup>20</sup></td><td>Figure 5-1</td></tr><tr><td>SSA</td><td>Creeping</td><td>30</td><td>-1.32</td><td>3.16 × 10<sup>20</sup></td><td>2.24 × 10<sup>20</sup></td><td>Figure 5-2</td></tr><tr><td>SSA</td><td>Creeping</td><td>100</td><td>-1.00</td><td>3.47 × 10<sup>20</sup></td><td>2.39 × 10<sup>20</sup></td><td>/</td></tr><tr><td>SSA</td><td>Creeping</td><td>200</td><td>-0.70</td><td>3.58 × 10<sup>20</sup></td><td>2.44 × 10<sup>20</sup></td><td>Figure 5-3</td></tr><tr><td>PS</td><td>Creeping</td><td>15</td><td>-1.58</td><td>2.02 × 10<sup>20</sup></td><td>1.48 × 10<sup>20</sup></td><td>Figure 8-1</td></tr><tr><td>PS</td><td>Creeping</td><td>30</td><td>-1.35</td><td>2.40 × 10<sup>20</sup></td><td>1.67 × 10<sup>20</sup></td><td>Figure 8-2</td></tr><tr><td>PS</td><td>Creeping</td><td>100</td><td>-0.65</td><td>2.86 × 10<sup>20</sup></td><td>1.85 × 10<sup>20</sup></td><td>/</td></tr><tr><td>PS</td><td>Creeping</td><td>200</td><td>-0.35</td><td>2.92 × 10<sup>20</sup></td><td>1.87 × 10<sup>20</sup></td><td>Figure 8-3</td></tr><tr><td>SSA</td><td>Locked</td><td>15</td><td>-1.53</td><td>4.15 × 10<sup>20</sup></td><td>2.96 × 10<sup>20</sup></td><td>/</td></tr><tr><td>SSA</td><td>Locked</td><td>30</td><td>-1.32</td><td>3.93 × 10<sup>20</sup></td><td>2.74 × 10<sup>20</sup></td><td>/</td></tr><tr><td>SSA</td><td>Locked</td><td>100</td><td>-1.00</td><td>3.75 × 10<sup>20</sup></td><td>2.54 × 10<sup>20</sup></td><td>/</td></tr><tr><td>SSA</td><td>Locked</td><td>200</td><td>-0.70</td><td>3.67 × 10<sup>20</sup></td><td>2.49 × 10<sup>20</sup></td><td>/</td></tr><tr><td>PS</td><td>Locked</td><td>15</td><td>-1.58</td><td>3.37 × 10<sup>20</sup></td><td>2.04 × 10<sup>20</sup></td><td>Figure S16-1 in Supporting Information S1</td></tr><tr><td>PS</td><td>Locked</td><td>30</td><td>-1.35</td><td>3.22 × 10<sup>20</sup></td><td>1.99 × 10<sup>20</sup></td><td>Figure S16-2 in Supporting Information S1</td></tr><tr><td>PS</td><td>Locked</td><td>100</td><td>-0.65</td><td>3.15 × 10<sup>20</sup></td><td>1.92 × 10<sup>20</sup></td><td>/</td></tr><tr><td>PS</td><td>Locked</td><td>200</td><td>-0.35</td><td>3.12 × 10<sup>20</sup></td><td>1.89 × 10<sup>20</sup></td><td>Figure S16-3 in Supporting Information S1</td></tr></tbody></table>

Note. Each line corresponds to a specific smoothing (λ and σ<sub>m</sub>), combined with a given reference frame (SSA or PS) and slip prior (creeping or locked).

zone stalled the 2001 M<sub>w</sub> = 8.4 Arequipa seismic rupture but was not efficient enough to arrest its southward propagation. The probability that an earthquake ruptures through a creeping barrier depends on the resistance and the efficiency of this barrier to arrest the seismic rupture propagation (Kaneko et al., 2010). The resistance of a barrier (R<sub>vs</sub> = ΔT<sub>vs</sub> × L<sub>vs</sub>) is the product to its length L<sub>vs</sub> (measured along the trench strike) and to the stress increase ΔT<sub>vs</sub> required for the velocity strengthening patch to sustain the seismic slip to propagate through it. The efficiency B of a barrier is defined as:

$$B = (\Delta T_{vs} / \beta \Delta T_{vw}) \cdot (L_{vs} / L_{vw}), \qquad (9)$$

where ΔT<sub>wv</sub> is the stress drop transferred from the ruptured velocity weakening segment of length L<sub>wv</sub> to the velocity strengthening patch and β a geometrical factor taken equal to 1/6 in the 3D models of Kaneko et al. (2010). A barrier is supposed to be strong when the probability of an earthquake rupturing through the barrier is low, which is the case when B > 1. Applied to the Nazca Ridge segment and considering L<sub>vs</sub> = 60 km, L<sub>wv</sub> = 180 km and β = 1/6 (from 3D models of Kaneko et al. (2010)), we find ΔT<sub>vs</sub> > ΔT<sub>wv</sub>/2, meaning that the required stress increase on the velocity strengthening patch should be higher than half of the seismic stress drop to stop the seismic rupture propagation.

When the efficiency tends to zero, this means that the probability that an earthquake ruptures through a creeping segment is high. If we consider that the NFZ is an inefficient barrier this would mean that ΔT<sub>vs</sub> tends to zero. Because ΔT<sub>vs</sub> is proportional to (a - b), this would suggest that NFZ segment might tend to a velocity neutral frictional behavior, efficient enough to slow down the seismic rupture propagation (Robinson et al., 2006) but inefficient to arrest it as attested by historical records (Figure 9).

On Figure S7 in Supporting Information S1 are displayed the thrust events with M<sub>w</sub> ≥ 5.5 (from gCMT), on top of a map indicating the dominating friction law, velocity weakening or velocity strengthening. Plots of the seismicity for various catalogs are given in Figures S8–S10 of Supporting Information S1, with cross-sections. We consider areas with coupling below 0.4 to be in the velocity strengthening regime, while areas with coupling above 0.4 are considered in the velocity weakening regime. We observe a convincing overlap between the seismic data set and the velocity weakening domain, both in along-strike and along-dip variations, which support our interseismic coupling distribution. While on Figure S12 in Supporting Information S1 we displayed all seismic
<figure><img src="image_12.png" /><figcaption>**Figure 9.** (left) Approximate rupture lengths of large subduction earthquakes (M<sub>w</sub> > 7.5) since the sixteenth century, taken from Villegas-Lanza et al. (2016). Recent earthquakes are displayed in blue, while historical earthquakes are displayed in gray. Extension of the Nazca Ridge and the Nazca Fracture Zone are reported in yellow. (right) Moment deficit rate by latitude for models in the Peruvian Sliver and Stable South America frame, and with various smoothing (λ of 15, 30, and 200 km), computed over 5-km wide segments. Patches with gray hachure pattern depict areas that are not well resolved by the model. Rupture extents from major seismic events are reported as blue segments.</figcaption></figure>

events, including small magnitude events, the IGP catalog being complete at M<sub>w</sub> ≥ 4.5. In this case, we observe a high cumulative seismicity at the location of creeping areas, especially where the Nazca ridge and the NFZ are subducting. This result was expected, as creeping areas are known to host a high amount of low magnitude seismic events (Rubin et al., 1999; Scholz, 1990).

### 5.4.2. Seismic Segmentation, Characteristic Earthquakes, and Return Period

In the Central Peru segment, our coupling maps recover about half of that segment and specifically the rupture areas of the 1974 M<sub>w</sub> = 7.9 Lima and 2007 M<sub>w</sub> = 8.0 Pisco earthquakes. Both 1974 and 2007 seismic ruptures appear to correlate quite well with two highly locked patches. At the rates of moment deficit of our coupling models, we found that the 1974 rupture area can easily host a similar event if the locked patch off-shore Lima breaks individually.

In the Southern Peru segment, historical and instrumental seismic catalogs indicate that this segment can be subdivided into three sub-segments corresponding to the rupture areas of the 1942 M<sub>w</sub> = 8.2 Nazca (200 km), 2001 M<sub>w</sub> = 8.4 Arequipa (300 km) and 1833 M<sub>w</sub> = 7.8 Ilo-Arica (150 km) earthquakes. The largest events of 1604 M<sub>w</sub> = 8.7 and 1868 M<sub>w</sub> = 8.8 appear to have ruptured simultaneously the Arequipa and Ilo-Arica sub-segments.

In the 1942 rupture area, the 1996 M<sub>w</sub> = 7.7 event as well as the 2013 and 2018 M<sub>w</sub> = 7.1 Acari events which all ruptured within the 1942 rupture suggesting that the locked patch there is relatively fragmented. Considering the
actual MDR integrated over the 1942 Nazca rupture, this segment would be mature today to host a $M_w \sim 8.0 \pm 0.1$ seismic event.

In the wake of the 2001 Arequipa and 2014 Iquique earthquakes, the 150-km long Ilo-Arica portion is clearly identified as a seismic gap. In the historical catalog that segment appears to have broken individually in 1833 and together with the 2001 Arequipa rupture area during the great 1868 $M_w = 8.8$ mega-event. Considering the moment deficit rates of our models, we found that the Ilo-Arica segment is mature for a seismic event with a moment magnitude $M_w > 8.0$.

Now if we consider the 450-km long rupture area of the great 1604 and 1868 $M_w \sim 8.7–8.8$ megathrust earthquakes, one way to estimate the return period $T(M_{\text{char}})$ of such characteristic earthquakes can be defined as:

$$T(M_{\text{char}}) = M_{\text{0char}}/(dM_{\text{inter}}/dt), \qquad (10)$$

where $dM_{\text{inter}}/dt$ is the MDR integrated over the characteristic rupture area. At the actual MDR integrated over the 450-km long rupture, we find a return period of $T_r \sim 103–171$ years, a period much shorter than the 264 years that separates these two great events. This difference can be partially explained by the fact that the seismic moments released by other large events in that same rupture area, like the 1784 $M_w = 8.4$ and 1687 $M_w = 8.0$ events are not taken into account with that approach, neither aseismic slip as postseismic that follow such events.

This deterministic approach is informative, but it tends to underestimate the return period of characteristic earthquakes and does not provide insight on the maximum-magnitude earthquake. Also, this approach does not take into account the frequency-magnitude of the background seismicity that varies from one region to another, nor aseismic slip events which indicate the limitations of that deterministic approach.

## 5.5. Moment Budget Analysis and Earthquake Potential

To conclude our analysis, we perform a moment budget at the scale of the South Peru segment since this segment is well resolved by our data, and because it cannot be excluded that it could rupture at once. For that we integrated the MDR of our interseismic models within a polygon that goes from the trench down to 60-km depth, and from the Nazca Ridge to the Arica bend (see polygon contours on Figure S11 in Supporting Information S1, the three southern polygons are considered).

To balance the seismicity rate with the MDR that is accumulating annually on the Southern Peru megathrust interface, we considered four instrumental catalogs that are homogeneous in magnitude (Figure 10): ISC-GEM (Di Giacomo et al., 2018; Storchak et al., 2013, 2015), ISC (Willemann & Storchak, 2001), gCMT (Dziewonski et al., 1981; Ekström et al., 2012), and Instituto Geofísico del Perú (IGP), covering respectively the 1906–2022, 1978–2022, 1976–2022, and 1906–2022 periods. The two main differences are that the global catalogs (ISC-GEM, ISC and gCMT) have a similar magnitude of completeness of $M_w = 5.2$ and a b-value of 0.68, 0.73, and 0.72 respectively, while the local catalog from IGP has a magnitude of completeness of $M_w = 4.5$ and a b-value of 1.03. Compared to others subduction zone, the Peruvian subduction zone is characterized by a low b-value (Bilek & Lay, 2018).

Assuming the frequency–magnitude distribution for these catalogs, the mean recurrence time of the associated earthquake potential $T(M_{\text{max}})$ is defined as (Avouac, 2015; Molnar, 1979):

$$T(M_{\text{max}}) = \frac{1}{\left(1 - \frac{2b}{3}\right)\alpha} \frac{M_{\text{max}}}{M_0} \qquad (11)$$

where $\alpha$ is the fraction of transient slip that is seismic during the cycle and is written as:

$$\alpha = \frac{M_{\text{seismic}}}{M_{\text{seismic}} + M_{\text{aseismic}}} \qquad (12)$$

with $M_{\text{seismic}}$ the seismic moment, $M_{\text{aseismic}}$ the aseismic moment including postseismic afterslip and slow slip episodes that occur during the cycle (Matsu'ura, 2015). When $\alpha = 1$ and $b = 0$, one will find the formulation of the recurrence time for a characteristic earthquake. Then, the frequency of the earthquake potential depends on the MDR and the b-value of seismic catalogs (see Figure 11) for which we have some estimations and the parameter $\alpha$ that we still need to bound.
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Frequency-magnitude distribution for various seismic catalogs on the South Peru segment: (a) ISC-GEM, (b) ISC, (c) gCMT, (d) IGP. The *b*-value of the corresponding Gutenberg-Richter law is reported for each catalog in the magnitude range 5.2–7.4 (blue line). Red stars correspond to strong events above magnitude 7.5, including historical events from Villegas-Lanza et al. (2016).</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 11.** Estimation of maximum expected magnitude with $\alpha = 0.7$ and $\mu = 40$ GPa, in the Peruvian Sliver frame and with a smoothing $\lambda = 30$ km, for the ISC-GEM catalog (left) and the IGP catalog (right). The earthquake potential is $M_w = 8.48$ with a 240 years recurrence time for the ISC-GEM catalog, while the earthquake potential is $M_w = 9.04$ with a 4,000 years recurrence time for the IGP catalog. Instrumental earthquakes from the catalogs are depicted by cyan dots, historical earthquakes from Villegas-Lanza et al. (2016) by green dots. The red lines are Gutenberg-Richter laws for maximum magnitude of 6, 7, 8, 9, and 10 from top to bottom respectively. The GR law of our catalog is shown by the pink line. The blue line is the frequency of the maximum magnitude event.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Figure 12.** Variation of the estimated earthquake potential in the Peruvian sliver (left) and in the Stable South America (right) frames, for various values of alpha, from 0.1 to 1. The benchmarked model has a regularization length $\lambda = 200$ km, a rigidity $\mu = 45$ GPa, and a b-value of 0.68 (blue) or 1.03 (green) corresponding to the ISC-GEM and IGP catalogs respectively. The minimum boundary for alpha is determined from the maximum-magnitude event of instrumental catalogs (which is the 2001 $M_w = 8.4$ Arequipa earthquake). The maximum boundary is considering a minimum 20% aseismic slip.</figcaption></figure>

The moment deficit rates of the family of acceptable coupling models range from 1.15 to $2.81 \times 10^{20}$ Nm/yr. However, we discussed in Section 5.2 that the visco-elastic relaxation following the 2001 Arequipa earthquake may still affect the observed velocity at some of the stations, notably close to the epicenter. The rupture area of this event was in the middle of the region considered for moment budget analysis, between the Nazca Ridge and the Arica bend. It is unlikely that all the stations in the area have been affected in the same proportion than the AREQ station, located in Arequipa city, still we adopt here a conservative approach and consider for uncertainties that the total moment deficit on the South Peru segment may be underestimated by up to 15% in the SSA frame and up to 30% in the PS frame. Measuring a velocity trend on less than a decade is not supposed to be strictly representative of the long-term secular motion, subsequently the trend velocity measured before the Arequipa earthquake could also be accelerated and may not be representative either of the secular interseismic trend.

Using the maximum moment deficit rates found of 2.4 (SSA)/1.9 (PS) $\times 10^{20}$ Nm/yr (models with $\lambda = 200$ km), we varied $\alpha$ between 0.1 and 1 for both catalogs to explore its impact on the estimate of the earthquake potential (Figure 12). A lower bound for $\alpha$ can be assessed from the largest magnitude event recorded in the instrumental catalogs, in our case the 2001 $M_w = 8.4$ earthquake, for which magnitude is reached when $\alpha = 0.21/0.27$ with the IGP catalog and $\alpha = 0.37/0.47$ with the ISC-GEM catalog. An upper bound for $\alpha$ can be estimated considering that the afterslip moment released during the postseismic period corresponds to at least 25% of the seismic moment (i.e., $\alpha \le 0.8$). This value is an average from various afterslip studies in subduction zones (Jara et al., 2018; Perfettini et al., 2010; Pritchard et al., 2007; Remy et al., 2016), and from our own analysis. In the specific case where $\alpha = 0.8$ we find an earthquake potential of $M_w = 9.6/9.4$ (Tr $\sim$5,300/3,400 years) with the IGP catalog and $M_w = 8.8/8.7$ (Tr $\sim$190/170 years) with the ISC-GEM catalog.

If now we suppose that during a full cycle, $M_{\text{seismic}} = M_{\text{aseismic}}$, then $\alpha = 0.5$ and we find an earthquake potential of $M_w = 9.2/9.0$ (Tr $\sim$2,100/1,350 years) with the IGP catalog and $M_w = 8.6/8.4$ (Tr $\sim$150/100 years) with the ISC-GEM catalog. In this case, the ISC-GEM catalog has an earthquake potential lower or equal to the maximum-magnitude of historical seismic records which is possible considering the uncertainties on the seismic moment for the 1604 and 1868 events, with a recurrence time of about Tr $\sim 125 \pm 25$ years consistent with the four $M_w > 8.4$ of the 500-years historical catalog. Considering the IGP catalog, this would suggest a maximum-magnitude earthquake of $M_w \ge 9.0$ with a return period of Tr $\ge 1350$ years much longer than the period covered by the historical catalog.

Finally, Figure 13 reports the maximum-magnitude and its variability for the various parameters and uncertainties discussed:

* The influence of the sliver motion, by correcting or not from the PS. Applying the sliver correction tends to lower the moment deficit (see Table 1).
<figure><img src="image_16.png" /><figcaption>**Figure 13.** (left) Earthquake potential in function of various parameters, from left to right: sliver motion (with and without Peruvian sliver (PS)), slip prior (creeping or locked), smoothing ($\lambda$ = 15, 30, and 200 km), average medium rigidity (35, 40, and 45 GPa), neglecting or not viscous behavior following the 2001 Arequipa earthquake, rate of seismic events $\alpha$ (0.3, 0.55, and 0.8), *b*-value of seismic catalog (ISC-GEM: 0.68, gCMT: 0.72, ISC: 0.73, and IGP: 1.03), and the total range. Models in the PS and Stable South America frames are weighted by a factor 2/3 and 1/3 respectively, and models accounting or not for viscous behavior following the Arequipa earthquake are weighted by a factor 1/3 and 2/3 respectively. The dots represent the median $M_w$ value, while the vertical bars represent the variability between the 15th and 85th percentiles on the 864 models considered. (right) Distribution of the earthquake potential for the 864 computed models, accounting for all the discussed factors of uncertainty. Median, 15th and 85th percentiles are depicted by dashed lines.</figcaption></figure>

* The slip prior of the model, considering a creeping (no coupling) or a locked (full coupling) slip prior. This parameter has a strong influence on poorly resolved area, for instance close to the trench.

* The smoothing of the interseismic slip model, with regularization length $\lambda$ = 15, 30, and 200 km. All these values allow a good fit to the data and a realistic slip distribution, but increasing the smoothing tends to increase the moment deficit (see Table 1).

* The rigidity $\mu$ = 35, 40, and 45 GPa. We consider this range plausible considering the average depth of our slip distribution (Bilek & Lay, 1999). The moment deficit is proportional to the rigidity.

* The impact of a potential underestimation of the moment deficit, due to the postseismic viscous relaxation following the 2001 Arequipa earthquake. We consider as an upper boundary an underestimation of 15% on the moment deficit in the SSA frame, and 30% in the PS frame.

* The rate of seismic events $\alpha$ = 0.3, 0.55, and 0.8. The lower boundary and upper boundaries were selected based on the analysis described above (see Figure 12).

* The seismic catalog considered, ISC-GEM ($b$ = 0.68), gCMT ($b$ = 0.72), ISC ($b$ = 0.73), and IGP ($b$ = 1.03).

We computed 864 combinations of these parameters, following the logical tree on Figure 14, and we displayed the median (circles), as well as the 15th and 85th percentiles (vertical bars), for each value of these parameters. Models in the PS and SSA frames are weighted by a factor 2/3 and 1/3 respectively, as we consider these models more likely (best fit to the data). In the same way, models accounting or not for an eventual moment underestimation following the Arequipa earthquake are weighted by a factor 1/3 and 2/3 respectively, as the correction corresponds to an extreme upper boundary.

The *b*-value of instrumental catalogs as well as the amount of aseismic slip in the transient slips are the two parameters that show the largest variability in the determination of the Mmax, with a variability of ~0.6 in $M_w$. The presence or not of the sliver, as well as the rigidity, or the viscoelastic relaxation associated with the Arequipa earthquake can influence that determination of ~0.15 in $M_w$. The slip prior induced a variability of ~0.12 in $M_w$, and the smoothing factor of the interseismic models a variability of ~0.08 in $M_w$.

Overall, we found a median value of $M_w$ = 8.55, a 15th percentile at $M_w$ = 8.15, and a 85th percentile at $M_w$ = 9.0. Considering the largest earthquake recorded on the South Peru segment was the 2001 $M_w$ = 8.4 Arequipa earthquake, we can state that the earthquake potential on this segment is at least $M_w$ = 8.4. The recurrence time for a $M_w$ = 8.4 would be ~100 years, while it would be ~1,000 years for a $M_w$ = 9.0. A better characterization of the *b*-value and the amount of aseismic slip over the cycle are the two main parameters to improve the assessment
<figure><img src="image_17.png" /><figcaption>**Figure 14.** Logic tree used for exploration of uncertainties on Figure 13. We explore from left to right: sliver motion (with and without Peruvian sliver (PS)), specified slip prior (creeping or locked), smoothing ($\lambda$ = 15, 30, and 200 km), average medium rigidity (35, 40, and 45 GPa), neglecting or not viscous behavior following the 2001 Arequipa earthquake, rate of seismic events $\alpha$ (0.3, 0.55, and 0.8), and b-value of seismic catalog (ISC-GEM: 0.68, gCMT: 0.72, ISC: 0.73, and IGP: 1.03). Models in the PS and Stable South America frames are weighted by a factor 2/3 and 1/3 respectively, and models accounting or not for viscous behavior following the Arequipa earthquake are weighted by a factor 1/3 and 2/3 respectively.</figcaption></figure>

of the earthquake potential. It is also essential to note that our moment budget analysis being restricted to the South Peru segment, our computed magnitudes are bounded by the size of this segment. In fact, an earthquake rupturing both the South Peru segment and another segment, crossing the Nazca Ridge or the Arica bend, could exceed our estimates.

## 6. Conclusions

To summarize, we have extracted a dense interseismic velocity field at the scale of the South Peruvian Andes from time series analysis of 73 GNSS sites (47 cGPS and 26 sGPS), including a significant amount of unpublished data. From that, we obtained a map of interseismic coupling, with an unprecedented resolution, through an inversion process based on the Tarantola-Valette formulation. It highlights a highly heterogeneous coupling distribution, with strongly coupled patches offshore Lima (12°S), Atico (16°S) and Ilo (17.5°S) cities. Low coupling areas were observed where the Nazca ridge, and more tenuously where the NFZ, are subducting below the continent. From this distribution, we estimate that the South Peru segment which extends from the Nazca ridge to the Arica bend could host a $M_w$ = 8.4 (recurrence time of 100 years) to $M_w$ = 9.0 (recurrence time of 900 years). In addition, we observe that the region near the coast of Lima could have a high seismic potential due to its very high coupling, according to our model. The estimation of the earthquake potential depends on several factors, like the chosen seismic catalog, the rate of seismic events over all transient parameters, correcting or not for a sliver motion, or accounting for long-term visco-elastic relaxation following the 2001 Arequipa earthquake. Among these factors, the two main sources of uncertainty are the seismic catalog and the rate of seismic events.

Further developments may include the use of InSAR data, allowing for a large-scale mapping with a very high resolution in space compared to GNSS. It would help to better constrain the depth of the transition between brittle and ductile rheology, as well as the amount and extension of intracontinental deformation. In addition, it would help to estimate the extent of visco-elastic relaxation following megathrust events, like the 2001 Arequipa earthquake. The increased resolution would also be a key point to overcome the lack of constraint on the models we encounter in some areas, for example, South of Lima. We may also include a denser monitoring of creeping areas to assess if the aseismic creep is released continuously or through bursts of slow slip, in order to better constrain the frictional behavior of those barriers, and the actual value of the rate of seismic events $\alpha$. Finally, the use of more complex rheology, with visco-elastic or visco-elasto-plastic behaviors would be a significant progress in order to account for large-scale intracontinental deformation (Li et al., 2015). It would also be a step in the direction of a general interseismic coupling model at the scale of central Andes, extending south of the Arica bend where a rotation of the long-term slip direction was suggested by Arriagada et al. (2008) and Métois et al. (2016).
## Data Availability Statement

All GNSS time series used in this manuscript are available at (Socquet, Tsapong-Tsague, et al., 2023): https://doi.org/10.17178/GNSS.products.SouthAmerica_GIPSYX.daily. RINEX data from CSN network are available at http://gps.csn.uchile.cl. RINEX data from RBMC are available at https://www.ibge.gov.br/en/geosciences/geodetic-positioning/geodetic-networks/19213-brazilian-network-for-continuous-monitoring-of-the-gnss-systems.html?=&t=downloads. RINEX data from SNAPP campaigns are available from UNAVCO (Simons et al., 2010): 1994 campaign at https://www.unavco.org/data/gps-gnss/data-access-methods/dail/data_request.php?gid=1325&ds=1&parent_link=Campaign&pview=original, 1996 campaign at https://www.unavco.org/data/gps-gnss/data-access-methods/dail/data_request.php?gid=1493&ds=2&parent_link=Campaign&pview=original, and 2001 campaign at https://www.unavco.org/data/gps-gnss/data-access-methods/dail/data_request.php?gid=2128&ds=1&parent_link=Campaign&pview=original. RINEX data from campaigns in South Peru between 2012 and 2022 are available through the following links: 2012 campaign (Socquet, Cotte, et al., 2023) at https://doi.org/10.15148/12160e27-0951-41b7-be97-efc4fec7ff96, 2013 campaign (Socquet, Norabuena, & Cotte, 2023) at https://doi.org/10.15148/14ab8f86-3453-4a44-b921-ec3b0a133ef6, 2015 campaign (Socquet & Norabuena, 2023) at https://doi.org/10.15148/376f13e4-e333-4da1-aa10-76f4eae08ffc, 2016 campaign (Nocquet & Socquet, 2023) at https://doi.org/10.15148/f410d57f-2e67-4174-9e7e-17ae3fd55d99, 2018 campaign (Socquet, Nocquet, et al., 2023) at https://doi.org/10.15148/94cc6056-84d5-4c92-a043-6db8b769babf, and 2022 campaign (Socquet, Villegas, et al., 2023) at https://gpscope.dt.insu.cnrs.fr/campagnes/data/2022-041/rinex/. RINEX data from IGP/ISTerre/Caltech continuous stations in South Peru are available at https://doi.osug.fr/data/public/GNSS_products/Peru/RINEX/peru_igp/. RINEX data from IGN network are available at https://doi.osug.fr/data/public/GNSS_products/Peru/RINEX/peru_ign/. The trajectory model analysis has been performed using the ITSA software, hosted on GriCAD GitLab repository (https://gricad-gitlab.univ-grenoble-alpes.fr/isterre-cycle/itsa), and described in Marill et al. (2021). Geographical illustrations were made using the Generic Mapping Tools version 6 package (Wessel et al., 2019) licensed under LGPL version 3 or later, available at https://www-generic-mapping-tools.org/. Other illustrations were made with Matplotlib v3.7.1 (Hunter, 2007), available under the Matplotlib license at https://matplotlib.org.

## Acknowledgments

This work was supported by funding from European Research Council (ERC), CNES, CNRS-INSU, and LabEx Osug@2020 (Investissements d'avenir—ANR10LABX56). A. Socquet received funding from the ERC CoG 865963 DEEP-trigger. B. Lovery salary was covered by Grants from CNES and ERC DEEPtrigger. We thank all organizations and persons that are supporting and operating continuous GNSS stations and campaign measurements used in this study. The continuous GNSS stations in Chile used in this work are maintained and distributed by the Centro Simologico Nacional, Universidad de Chile. The IGP/ISTerre/Caltech GNSS stations in Peru were installed by ISTerre and Caltech in collaboration with IGP, and are now maintained by the Instituto Geofisico del Peru. The IGN network in Peru is maintained by the IGN. We thank the French national pool of portable geodetic instruments GPSMOB-RESIF (INSU-CNRS) for providing the geodetic instruments for the GNSS campaigns. Finally, we thank former engineers who worked on the operational implementation of the GipsyX GNSS processing at ISTerre, including Gaëlle Deschamps-Huygen and Rémi Molaro-Maqua.

## References

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

Arriagada, C., Roperch, P., Mpodozis, C., & Cobbold, P. R. (2008). Paleogene building of the Bolivian Orocline: Tectonic restoration of the central Andes in 2-D map view. *Tectonics*, **27**(6), TC6014. https://doi.org/10.1029/2008TC002269

Avouac, J.-P. (2015). 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. https://doi.org/10.1146/annurev-earth-060614-105302

Beck, S. L., & Nishenko, S. P. (1990). Variations in the mode of great earthquake rupture along the central Peru subduction zone. *Geophysical Research Letters*, **17**(11), 1969–1972. https://doi.org/10.1029/GL017i011p01969

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

Bertiger, W., Bar-Sever, Y., Dorsey, A., Haines, B., Harvey, N., Hemberger, D., et al. (2020). GipsyX/RTGx, a new tool set for space geodetic operations and research. *Advances in Space Research*, **66**(3), 469–489. https://doi.org/10.1016/j.asr.2020.04.015

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

Bevis, M., Kendrick, E., Smalley, R., Jr., Brooks, B., Allmendinger, R., & Isacks, B. (2001). On the strength of interplate coupling and the rate of back arc convergence in the central Andes: An analysis of the interseismic velocity field. *Geochemistry, Geophysics, Geosystems*, **2**(11), 1067. https://doi.org/10.1029/2001GC000198

Bilek, S. L. (2010). Invited review paper: Seismicity along the South American subduction zone: Review of large earthquakes, tsunamis, and subduction zone complexity. *Tectonophysics*, **495**(1–2), 2–14. https://doi.org/10.1016/j.tecto.2009.02.037

Bilek, S. L., & Lay, T. (1999). Rigidity variations with depth along interplate megathrust faults in subduction zones. *Nature*, **400**(6743), 443–446. https://doi.org/10.1038/22739

Bilek, S. L., & Lay, T. (2018). Subduction zone megathrust earthquakes. *Geosphere*, **14**(4), 1468–1500. https://doi.org/10.1130/GES01608.1

Blewitt, G., Kreemer, C., Hammond, W. C., & Gazeaux, J. (2016). MIDAS robust trend estimator for accurate GPS station velocities without step detection. *Journal of Geophysical Research: Solid Earth*, **121**(3), 2054–2068. https://doi.org/10.1002/2015JB012552

Boehm, J., Werl, B., & Schuh, H. (2006). 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*, **111**(B2), B02406. https://doi.org/10.1029/2005JB003629

Bürgmann, R., Kogan, M. G., Steblov, G. M., Hilley, G., Levin, V. E., & Apel, E. (2005). Interseismic coupling and asperity distribution along the Kamchatka subduction zone. *Journal of Geophysical Research*, **110**(B7), B07405. https://doi.org/10.1029/2005JB003648

Chlieh, M., Avouac, J.-P., Hjorleifsdottir, V., Song, T.-R. A., Ji, C., Sieh, K., et al. (2007). Coseismic slip and afterslip of the great Mw 9.15 Sumatra-Andaman earthquake of 2004. *Bulletin of the Seismological Society of America*, **97**(1A), S152–S173. https://doi.org/10.1785/0120050631

Chlieh, M., Avouac, J. P., Sieh, K., Natawidjaja, D. H., & Galetzka, J. (2008). Heterogeneous coupling of the Sumatran megathrust constrained by geodetic and paleogeodetic measurements. *Journal of Geophysical Research*, **113**(B5), B05305. https://doi.org/10.1029/2007JB004981
Chlieh, M., Perfettini, H., Tavera, H., Avouac, J.-P., Remy, D., Nocquet, J.-M., et al. (2011). Interseismic coupling and seismic potential along the Central Andes subduction zone. *Journal of Geophysical Research*, **116**(B12), B12405. https://doi.org/10.1029/2010JB008166

Cifuentes, I. L. (1989). The 1960 Chilean earthquakes. *Journal of Geophysical Research*, **94**(B1), 665–680. https://doi.org/10.1029/JB094iB01p00665

Di Giacomo, D., Engdahl, E. R., & Storchak, D. A. (2018). The ISC-GEM Earthquake Catalogue (1904–2014): Status after the Extension Project. *Earth System Science Data*, **10**(4), 1877–1899. https://doi.org/10.5194/essd-10-1877-2018

Dorbath, L., Cisternas, A., & Dorbath, C. (1990). Assessment of the size of large and great historical earthquakes in Peru. *Bulletin of the Seismological Society of America*, **80**, 551–576. https://doi.org/10.1785/BSSA0800030551

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

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

Florsch, N., Revil, A., & Camerlynck, C. (2014). Inversion of generalized relaxation time distributions with optimized damping parameter. *Journal of Applied Geophysics*, **109**, 119–132. https://doi.org/10.1016/j.jappgeo.2014.07.013

Gagnon, K., Chadwell, C. D., & Norabuena, E. (2005). Measuring the onset of locking in the Peru–Chile trench with GPS and acoustic measurements. *Nature*, **434**(7030), 205–208. https://doi.org/10.1038/nature03412

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

Hergert, T., & Heidbach, O. (2006). New insights into the mechanism of postseismic stress relaxation exemplified by the 23 June 2001 $M_w = 8.4$ earthquake in southern Peru. *Geophysical Research Letters*, **33**(2), L02307. https://doi.org/10.1029/2005GL024858

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

Hoffmann, F., Metzger, S., Moreno, M., Deng, Z., Sippl, C., Ortega-Culaciati, F., & Oncken, O. (2018). Characterizing afterslip and ground displacement rate increase following the 2014 Iquique-Pisagua $M_w = 8.1$ earthquake, Northern Chile. *Journal of Geophysical Research: Solid Earth*, **123**(5), 4171–4192. https://doi.org/10.1002/2017JB014970

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

Ide, S., Baltay, A., & Beroza, G. C. (2011). Shallow dynamic overshoot and energetic deep rupture in the 2011 Mw 9.0 Tohoku-Oki earthquake. *Science*, **332**(6036), 1426–1429. https://doi.org/10.1126/science.1207020

Itoh, Y., Wang, K., Nishimura, T., & He, J. (2019). Compliant volcanic arc and backarc crust in southern Kurile suggested by interseismic geodetic deformation. *Geophysical Research Letters*, **46**(21), 11790–11798. https://doi.org/10.1029/2019GL084656

Jara, J., Sánchez-Reyes, H., Socquet, A., Cotton, F., Virieux, J., Maksymowicz, A., et al. (2018). 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. https://doi.org/10.1016/j.epsl.2018.09.025

Johnston, G., Riddell, A., & Hausler, G. (2017). The International GNSS Service. In P. J. G. Teunissen & O. Montenbruck (Eds.), *Springer Handbook of Global Navigation Satellite Systems* (1st ed., pp. 967–982). Springer International Publishing. https://doi.org/10.1007/978-3-319-42928-1

Kanamori, H., & Cipar, J. J. (1974). Focal process of the great Chilean earthquake May 22, 1960. *Physics of the Earth and Planetary Interiors*, **9**(2), 128–136. https://doi.org/10.1016/0031-9201(74)90029-6

Kaneko, Y., Avouac, J.-P., & Lapusta, N. (2010). Towards inferring earthquake patterns from geodetic observations of interseismic coupling. *Nature Geoscience*, **3**(5), 363–369. https://doi.org/10.1038/ngeo843

Kendrick, E., Bevis, M., Smalley, R., Brooks, B., Vargas, R. B., Lauría, E., & Fortes, L. P. S. (2003). The Nazca–South America Euler vector and its rate of change. *Journal of South American Earth Sciences*, **16**(2), 125–131. https://doi.org/10.1016/S0895-9811(03)00028-2

King, G. (1986). Speculations on the geometry of the initiation and termination processes of earthquake rupture and its relation to morphology and geological structure. *Pure and Applied Geophysics*, **124**(3), 567–585. https://doi.org/10.1007/bf00877216

King, G., & Nábělek, J. (1985). Role of fault bends in the initiation and termination of earthquake rupture. *Science*, **228**(4702), 984–987. https://doi.org/10.1126/science.228.4702.984

Klein, E., Fleitout, L., Vigny, C., & Garaud, J. D. (2016). Afterslip and viscoelastic relaxation model inferred from the large-scale post-seismic deformation following the 2010 $M_w = 8.8$ Maule earthquake (Chile). *Geophysical Journal International*, **205**(3), 1455–1472. https://doi.org/10.1093/gji/ggw086

Kley, J., & Monaldi, C. R. (1998). Tectonic shortening and crustal thickness in the Central Andes: How good is the correlation? *Geology*, **26**(8), 723–726. https://doi.org/10.1130/0091-7613(1998)026<0723:TSACTI>2.3.CO;2

Koketsu, K., Yokota, Y., Nishimura, N., Yagi, Y., Miyazaki, S., Satake, K., et al. (2011). A unified source model for the 2011 Tohoku earthquake. *Earth and Planetary Science Letters*, **310**(3), 480–487. https://doi.org/10.1016/j.epsl.2011.09.009

Konca, A. O., Avouac, J.-P., Sladen, A., Meltzner, A. J., Sieh, K., Fang, P., et al. (2008). Partial rupture of a locked patch of the Sumatra megathrust during the 2007 earthquake sequence. *Nature*, **456**(7222), 631–635. https://doi.org/10.1038/nature07572

Kositsky, A. P., & Avouac, J.-P. (2010). Inverting geodetic time series with a principal component analysis-based inversion method. *Journal of Geophysical Research*, **115**(B3), B03401. https://doi.org/10.1029/2009JB006535

Langer, C. J., & Spence, W. (1995). The 1974 Peru earthquake series. *Bulletin of the Seismological Society of America*, **85**(3), 665–687.

Li, S., & Chen, L. (2023). How long can the postseismic and interseismic phases of great subduction earthquake sustain? Toward an integrated earthquake-cycle perspective. *Geophysical Research Letters*, **50**(11), e2023GL0103976. https://doi.org/10.1029/2023GL0103976

Li, S., Moreno, M., Bedford, J., Rosenau, M., & Oncken, O. (2015). Revisiting viscoelastic effects on interseismic deformation and locking degree: A case study of the Peru-North Chile subduction zone. *Journal of Geophysical Research: Solid Earth*, **120**(6), 4522–4538. https://doi.org/10.1002/2015JB011903

Loveless, J. P., & Meade, B. J. (2011). Spatial correlation of interseismic coupling and coseismic rupture extent of the 2011 $M_w = 9.0$ Tohoku-oki earthquake. *Geophysical Research Letters*, **38**(17), L17306. https://doi.org/10.1029/2011GL048561

Lyard, F. H., Allain, D. J., Cancet, M., Carrère, L., & Picot, N. (2021). FES2014 global ocean tide atlas: Design and performance. *Ocean Science*, **17**(3), 615–649. https://doi.org/10.5194/os-17-615-2021

Marill, L., Marsan, D., Socquet, A., Radiguet, M., Cotte, N., & Rousset, B. (2021). Fourteen-year acceleration along the Japan Trench. *Journal of Geophysical Research: Solid Earth*, **126**(11), e2020JB021226. https://doi.org/10.1029/2020JB021226

Mariniere, J., Beauval, C., Nocquet, J.-M., Chlieh, M., & Yepes, H. (2021). Earthquake recurrence model for the Colombia–Ecuador subduction zone constrained from seismic and geodetic data, implication for PSHA. *Bulletin of the Seismological Society of America*, **111**(3), 1508–1528. https://doi.org/10.1785/0120200338
Matsu'ura, T. (2015). Late Quaternary uplift rate inferred from marine terraces, Muroto Peninsula, southwest Japan: Forearc deformation in an oblique subduction zone. *Geomorphology*, **234**, 133–150. https://doi.org/10.1016/j geomorph.2015.01.012

Meade, B. J. (2007). Algorithms for the calculation of exact displacements, strains, and stresses for triangular dislocation elements in a uniform elastic half space. *Computers & Geosciences*, **33**(8), 1064–1075. https://doi.org/10.1016/j.cageo.2006.12.003

Métois, M., Socquet, A., Vigny, C., Carrizo, D., Peyrat, S., Delorme, A., et al. (2013). Revisiting the North Chile seismic gap segmentation using GPS-derived interseismic coupling. *Geophysical Journal International*, **194**(3), 1283–1294. https://doi.org/10.1093/gji/ggt183

Métois, M., Vigny, C., & Socquet, A. (2016). Interseismic coupling, megathrust earthquakes and seismic swarms along the Chilean subduction zone (38°–18°S). *Pure and Applied Geophysics*, **173**(5), 1431–1449. https://doi.org/10.1007/s00024-016-1280-5

Molnar, P. (1979). Earthquake recurrence intervals and plate tectonics. *Bulletin of the Seismological Society of America*, **69**(1), 115–133. https://doi.org/10.1785/BSSA0690010115

Moreno, M., Rosenau, M., & Oncken, O. (2010). 2010 Maule earthquake slip correlates with pre-seismic locking of Andean subduction zone. *Nature*, **467**(7312), 198–202. https://doi.org/10.1038/nature09349

Nishenko, S. P. (1991). *Circum-Pacific seismic potential*. U.S. Geological Survey, National Earthquake Information Center.

Nocquet, J.-M., & Socquet, A. (2023). EPOS-France - GPSMob data - Mission n° 16-153 - Perou (2016) - 2016-08-24/2016-12-06 - 36 points [Dataset]. OSU OREME. https://doi.org/10.15148/F410D57F-2E67-4174-9E7-17AE3FD55D99

Nocquet, J.-M., Villegas-Lanza, J. C., Chlieh, M., Mothes, P. A., Rolandone, F., Jarrin, P., et al. (2014). Motion of continental slivers and creeping subduction in the northern Andes. *Nature Geoscience*, **7**(4), 287–291. https://doi.org/10.1038/ngeo2099

Norabuena, E., Leffler-Griffin, L., Mao, A., Dixon, T., Stein, S., Sacks, I. S., et al. (1998). Space geodetic observations of Nazca-South America convergence across the central Andes. *Science*, **279**(5349), 358–362. https://doi.org/10.1126/science.279.5349.358

Perfettini, H., Avouac, J.-P., & Ruegg, J.-C. (2005). Geodetic displacements and aftershocks following the 2001 M<sub>w</sub> = 8.4 Peru earthquake: Implications for the mechanics of the earthquake cycle along subduction zones. *Journal of Geophysical Research*, **110**(B9), B09404. https://doi.org/10.1029/2004JB003522

Perfettini, H., Avouac, J.-P., Tavera, H., Kositsky, A., Nocquet, J.-M., Bondoux, F., et al. (2010). Seismic and aseismic slip on the Central Peru megathrust. *Nature*, **465**(7294), 78–81. https://doi.org/10.1038/nature09062

Pritchard, M. E., Norabuena, E. O., Ji, C., Boroschek, R., Comte, D., Simons, M., et al. (2007). Geodetic, teleseismic, and strong motion constraints on slip from recent southern Peru subduction zone earthquakes. *Journal of Geophysical Research*, **112**(B3), B03307. https://doi.org/10.1029/2006JB004294

Radiguet, M., Cotton, F., Vergnolle, M., Campillo, M., Valette, B., Kostoglakov, V., & Cotte, N. (2011). Spatial and temporal evolution of a long term slow slip event: The 2006 Guerrero Slow Slip Event. *Geophysical Journal International*, **184**(2), 816–828. https://doi.org/10.1111/j.1365-246X.2010.04866.x

Radiguet, M., Perfettini, H., Cotte, N., Gualandi, A., Valette, B., Kostoglakov, V., et al. (2016). Triggering of the 2014 M<sub>w</sub> 7.3 Papanoa earthquake by a slow slip event in Guerrero, Mexico. *Nature Geoscience*, **9**(11), 829–833. https://doi.org/10.1038/ngeo2817

Remy, D., Perfettini, H., Cotte, N., Avouac, J. P., Chlieh, M., Bondoux, F., et al. (2016). Postseismic relocking of the subduction megathrust following the 2007 Pisco, Peru, earthquake. *Journal of Geophysical Research: Solid Earth*, **121**(5), 3978–3995. https://doi.org/10.1002/2015JB012417

Robinson, D. P., Das, S., & Watts, A. B. (2006). Earthquake rupture stalled by a subducting fracture zone. *Science*, **312**(5777), 1203–1205. https://doi.org/10.1126/science.1125771

Rubin, A. M., Gillard, D., & Got, J.-L. (1999). Streaks of microearthquakes along creeping faults. *Nature*, **400**(6745), 635–641. https://doi.org/10.1038/23196

Saillard, M., Audin, L., Rousset, B., Avouac, J.-P., Chlieh, M., Hall, S. R., et al. (2017). From the seismic cycle to long-term deformation: Linking seismic coupling and Quaternary coastal geomorphology along the Andean megathrust. *Tectonics*, **36**(2), 241–256. https://doi.org/10.1002/2016TC004156

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

Scholz, C. H. (1990). *The mechanics of earthquakes and faulting* (p. 413). Cambridge University Press.

Silgado, E. (1978). Recurrence of tsunamis in the western coast of South America. *Marine Geodesy*, **1**(4), 347–354. https://doi.org/10.1080/01490417809387980

Simons, M., Galetzka, J., Genrich, J. F., Ortega, F., Comte, D., Glass, B., et al. (2010). Central Andean Tectonic Observatory Geodetic Array [Dataset]. UNAVCO, Inc. https://doi.org/10.7283/T50P0X037

Simons, M., Minson, S. E., Sladen, A., Ortega, F., Jiang, J., Owen, S. E., et al. (2011). The 2011 magnitude 9.0 Tohoku-Oki earthquake: Mosaicking the megathrust from seconds to centuries. *Science*, **332**(6036), 1421–1425. https://doi.org/10.1126/science.1206731

Sladen, A., Tavera, H., Simons, M., Avouac, J. P., Konca, A. O., Perfettini, H., et al. (2010). Source model of the 2007 M<sub>w</sub> 8.0 Pisco, Peru earthquake: Implications for seismogenic behavior of subduction megathrusts. *Journal of Geophysical Research*, **115**(B2), B02405. https://doi.org/10.1029/2009JB006429

Socquet, A., Cotte, N., Bondoux, F., & Norabuena, E. (2023). EPOS-France - GPSMob data - Mission n° 12-136b - Sud Pérou (2012) - 2012-06-23/2012-07-04 - 24 points [Dataset]. OSU OREME. https://doi.org/10.15148/12160E27-0951-41B7-BE97-EFC4FEC7FF96

Socquet, A., Nocquet, J.-M., Norabuena, E., Villegas, J. C., & Chlieh, M. (2023). EPOS-France - GPSMob data - Mission n° 18-135 - Perou (2018) - 2018-11-02/2018-11-21 - 26 points [Dataset]. OSU OREME. https://doi.org/10.15148/94CC6056-84D5-4C92-A043-6DB8B769BABF

Socquet, A., & Norabuena, E. (2023). EPOS-France - GPSMob data - Mission n° 15-265 - Pérou (2015) - 2015-10-22/2015-11-01 - 27 points [Dataset]. OSU OREME. https://doi.org/10.15148/376F13E4-E333-4DA1-AA10-76F4EAE08FFC

Socquet, A., Norabuena, E., & Cotte, N. (2023). EPOS-France - GPSMob data - Mission n° 13-136 - Sud Pérou (2013) - 2013-06-09/2013-06-15 - 2 points [Dataset]. OSU OREME. https://doi.org/10.15148/14AB8F86-3453-4A44-B921-EC3B0A133EF6

Socquet, A., Tsapong-Tsague, A., Lovery, B., Janex, G., & Radiguet, M. (2023). Metadata and GNSS daily position solutions for permanent GNSS stations in South America [Dataset]. These datasets include station metadata and GNSS daily solutions for over 700 GNSS stations in South America. CNRS, OSUG, ISTERRE. https://doi.org/10.17178/GNSS.PRODUCTS.SOUTHAMERICA_GIPSYX.DAILY

Socquet, A., Villegas, J. C., Langlais, M., Doubre, C., Donniol, B., Costantino, G., et al. (2023). EPOS-France - GPSMob data - Mission n° 22-041 - Perou (2022) - 2022-03-12/2022-08-07 - 69 points [Dataset]. OSU OREME. https://doi.org/10.15148/289f31dd-9109-49ba-8116-5217191226e5

Stein, S., & Okal, E. A. (2005). Speed and size of the Sumatra earthquake. *Nature*, **434**(7033), 581–582. https://doi.org/10.1038/434581a

Storchak, D. A., Di Giacomo, D., Bondár, I., Engdahl, E. R., Harris, J., Lee, W. H. K., et al. (2013). Public Release of the ISC-GEM Global Instrumental Earthquake Catalogue (1900-2009). *Seismological Research Letters*, **84**(5), 810–815. https://doi.org/10.1785/0220130034

Storchak, D. A., Di Giacomo, D., Engdahl, E. R., Harris, J., Bondár, I., Lee, W. H. K., et al. (2015). The ISC-GEM Global Instrumental Earthquake Catalogue (1900-2009): Introduction. *Physics of the Earth and Planetary Interiors*, **239**, 48–63. https://doi.org/10.1016/j.pepi.2014.06.009
Takabatake, T., Shibayama, T., Esteban, M., & Ishii, H. (2018). Advanced casualty estimation based on tsunami evacuation intended behavior: Case study at Yuigahama Beach, Kamakura, Japan. *Natural Hazards*, *92*(3), 1763–1788. https://doi.org/10.1007/s11069-018-3277-0

Tarantola, A. (2005). *Inverse problem theory and methods for model parameter estimation*. Society for Industrial and Applied Mathematics.

Tarantola, A., & Valette, B. (1982). Generalized nonlinear inverse problems solved using the least squares criterion. *Reviews of Geophysics*, *20*(2), 219–232. https://doi.org/10.1029/RG020i002p00219

Tavera, H., & Buforn, E. (1998). Sismicidad y sismotectónica de Perú. *Física de la Tierra*, *10*, 187–219.

Tavera, H., Buforn, E., Bernal, I., Antayhua, Y., & Vilacapoma, L. (2002). The Arequipa (Peru) earthquake of June 23, 2001. *Journal of Seismology*, *6*(2), 279–283. https://doi.org/10.1023/a:1015698621075

Veloza, G., Styron, R., Taylor, M., & Mora, A. (2011). Open source archive of active faults for northwest South America. *Geological Society of America Today*, *22*(10), 4–10. https://doi.org/10.1130/GSAT-G156A.1

Villegas-Lanza, J. C., Chlieh, M., Cavalié, O., Tavera, H., Baby, P., Chire-Chira, J., & Nocquet, J.-M. (2016). 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. https://doi.org/10.1002/2016JB013080

Wessel, P., Luis, J. F., Uieda, L., Scharroo, R., Wobbe, F., Smith, W. H. F., & Tian, D. (2019). The Generic Mapping Tools version 6. *Geochemistry, Geophysics, Geosystems*, *20*(11), 5556–5564. https://doi.org/10.1029/2019GC008515

Willemann, R. J., & Storchak, D. A. (2001). Data Collection at the International Seismological Centre. *Seismological Research Letters*, *72*(4), 440–453. https://doi.org/10.1785/gssrl.72.4.440

Yáñez-Cuadra, V., Ortega-Culaciati, F., Moreno, M., Tassara, A., Krumm-Nualart, N., Ruiz, J., et al. (2022). Interplate coupling and seismic potential in the Atacama Seismic Gap (Chile): Dismissing a rigid Andean sliver. *Geophysical Research Letters*, *49*(11), e2022GL098257. https://doi.org/10.1029/2022GL098257