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

The Central Andes subduction has been the theater of numerous large earthquakes since the beginning of the 21st Century, notably the 2001 Mw = 8.4 Arequipa, 2007 Mw = 8.0 Pisco and 2014 Mw = 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 Mw = 8.4–9.0 earthquake potential depending principally on the considered seismic catalog and the seismic/aseismic slip ratio.


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 w ≥ 8.5) earthquakes that occurred in 1746 in the Central and in 1604 and 1868 in the Southern segments.Large (M w ≥ 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).In the North of Chile where the great 1877 M w = 8.8 earthquake occurred, the 2014 M w = 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 w = 8.2 and the 2001 M w = 8.4 Arequipa earthquakes are the two largest events since the great 1868 M w = 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 w ∼ 8.0 occurred in 1940in , 1966in , 1974in , and 2007 (Pisco) (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  (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 w = 8.1 2014 Iquique and M w = 8.4 2001 Arequipa earthquakes have been modeled in this study.M w = 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.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).

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.All the GNSS time series with their trajectory model are available in Supporting Information S1.

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 −1 (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).

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: (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.
• 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 τ Ra and a displacement a i at time t i .A value of τ 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)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019)(2020)(2021)(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.

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 λ = 10 km and a damping parameter σ 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.

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.

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 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 μ 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: 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: and the modeled slip rates m on each source element: 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: (5) d(k,l) being the distance between source elements k and l, λ 0 a scaling factor (characteristic size of fault patch), λ the spatial correlation length, and σ m the damping parameter controlling the relative weight of the distance to initial model versus data misfit.In practice, we adjust the regularization length λ and σ 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: d being the data, m being the model, σ 2 the variance, and ν = 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: 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 Id, 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.

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 χ 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 λ of 15, 30, or 200 km that all provide an acceptable GNSS misfit with χ 2 ranging from 0.28 to 0.47. Figure 5 reports the coupling distributions for various smoothing factors λ = 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: 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 × 10 20 Nm/yr respectively for models with λ = 15 and 200 km.

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 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 (χ 2 = 0.24 with the sliver, while χ 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.(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).

Comparison With Previously Published Interseismic Coupling Models
Compared to Villegas-Lanza et al. ( 2016) and Chlieh et al. (2011), our results confirm some features while refining 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) 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).

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)(1998)(1999)(2000)(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.  .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 (λ 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.

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 χ 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 χ 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).

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 × 10 20 Nm/yr (∼M w = 7.4-7.6),respectively for λ varying between 15 and 200 km (Table 1 and Table S4 in Supporting Information S1).The MDR increases with λ, however this increase is relatively low in the South Peru segment where it evolves from 1.4 to 1.8 × 10 20 Nm/yr, compared to the Central Peru segment where the MDR dispersion between models is much higher (0.6-1.1 × 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.

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 zone stalled the 2001 M w = 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 vs = ΔT vs × L vs ) is the product to its length L vs (measured along the trench strike) and to the stress increase ΔT vs required for the velocity strengthening patch to sustain the seismic slip to propagate through it.The efficiency B of a barrier is defined as: where ΔT vw is the stress drop transferred from the ruptured velocity weakening segment of length L vw 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 vs = 60 km, L vw = 180 km and β = 1/6 (from 3D models of Kaneko et al. ( 2010)), we find ΔT vs > ΔT vw /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 vs tends to zero.
Because ΔT vs 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 w ≥ 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 Note.Each line corresponds to a specific smoothing (λ and σ m ), combined with a given reference frame (SSA or PS) and slip prior (creeping or locked).

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 10.1029/2023JB027114 16 of 25 events, including small magnitude events, the IGP catalog being complete at M w ≥ 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).

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 w = 7.9 Lima and 2007 M w = 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 w = 8.2 Nazca (200 km), 2001 M w = 8.4 Arequipa (300 km) and 1833 M w = 7.8 Ilo-Arica (150 km) earthquakes.The largest events of 1604 M w = 8.7 and 1868 M w = 8.8 appear to have ruptured simultaneously the Arequipa and Ilo-Arica sub-segments.
In the 1942 rupture area, the 1996 M w = 7.7 event as well as the 2013 and 2018 M w = 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 ∼ 8.0 ± 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 ∼ 8.7-8.8 megathrust earthquakes, one way to estimate the return period T(Mw char ) of such characteristic earthquakes can be defined as: where dMo 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 Tr ∼ 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.

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., 2013Storchak et al., , 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 max ) is defined as (Avouac, 2015;Molnar, 1979): where α is the fraction of transient slip that is seismic during the cycle and is written as: with  Mo seismic the seismic moment,  Mo aseismic the aseismic moment including postseismic afterslip and slow slip episodes that occur during the cycle (Matsu'ura, 2015).When α = 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 α that we still need to bound.The moment deficit rates of the family of acceptable coupling models range from 1.15 to 2.81 × 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) × 10 20 Nm/yr (models with λ = 200 km), we varied α 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 α 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 α = 0.21/0.27with the IGP catalog and α = 0.37/0.47 with the ISC-GEM catalog.An upper bound for α can be estimated considering that the afterslip moment released during the postseismic period corresponds to at least 25% of the seismic moment (i.e., α ≤ 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 α = 0.8 we find an earthquake potential of M w = 9.6/9.4(Tr ∼5,300/3,400 years) with the IGP catalog and M w = 8.8/8.7 (Tr ∼190/170 years) with the ISC-GEM catalog.
If now we suppose that during a full cycle, Mo seismic = Mo aseismic , then α = 0.5 and we find an earthquake potential of M w = 9.2/9.0(Tr ∼2,100/1,350 years) with the IGP catalog and M w = 8.6/8.4 (Tr ∼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 ∼125 ± 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 ≥ 9.0 with a return period of Tr ≥ 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).• 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 λ = 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 μ = 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).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 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.

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 α.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).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.

Figure 1 .
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 w = 8.1 2014 Iquique and M w = 8.4 2001 Arequipa earthquakes have been modeled in this study.M w = 8.0 2007 Pisco rupture is reported fromChlieh 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 fromVeloza 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.

Figure 2 .
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.

Figure 3 .
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.

Figure 4 .
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 λ = 30 km) being depicted by the gray line.The topography is plotted above.

Figure 5 .
Figure5.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 w = 8.4 2001 Arequipa has been estimated in this study, while the M w = 8.0 2007 Pisco earthquake rupture is reported fromChlieh 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.
Norabuena et al. (1998) andBevis 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, whileBevis 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 fromChlieh et al. (2011) andVillegas-Lanza et al. (2016) were able to present more detailed interseismic coupling distributions thanks to new developments.

Figure 6 .
Figure6.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 λ = 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).

Figure 7 .
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.

Figure 8 .
Figure 8. Interseismic coupling distribution from our inversion of interseismic velocities in the Peruvian Sliver (PS) frame, with three different smoothing (λ 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 fromChlieh 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 (λ 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.

Figure 9 .
Figure 9. (left) Approximate rupture lengths of large subduction earthquakes (M w > 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.

Figure 11 .
Figure 11.Estimation of maximum expected magnitude with α = 0.7 and μ = 40 GPa, in the Peruvian Sliver frame and with a smoothing λ = 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.

Figure 10 .
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).

Figure 12 .
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 λ = 200 km, a rigidity μ = 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.

Figure 13 .
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 (λ = 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 α (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.

Figure 14 .
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 (λ = 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 α (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.