# Self-reactivated rupture during the 2019 Mw = 8 northern Peru intraslab earthquake

Martin Vallée, Yuqing Xie, Raphaël Grandin, Juan Carlos Villegas-Lanza, Jean-Mathieu Nocquet, Sandro Vaca, Lingsen Meng, Jean-Paul Ampuero, Patricia Mothes, Paul Jarrin, et al.

► To cite this version:

Martin Vallée, Yuqing Xie, Raphaël Grandin, Juan Carlos Villegas-Lanza, Jean-Mathieu Nocquet, et al.. Self-reactivated rupture during the 2019 Mw = 8 northern Peru intraslab earthquake. Earth and Planetary Science Letters, 2023, 601, 10.1016/j.epsl.2022.117886. insu-03947012

HAL Id: insu-03947012

https://insu.hal.science/insu-03947012v1

Submitted on 23 Jul 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Self-reactivated rupture during the 2019 $M_w$=8
northern Peru intraslab earthquake

Martin Valléeª, Yuging Xieᵇ, Raphaël Grandinª, Juan Carlos Villegas-Lanzaᶜ, Jean-Mathieu Nocquetd,ª, Sandro Vacaᵉ, Lingsen Mengᵇ, Jean Paul Ampueroᵈ, Patricia Mothesᵉ, Paul Jarrinf,ᵈ, Ciro Sierra Farfang, Frédérique Rolandonef

## Abstract

The 2019/05/26 Northern Peru earthquake ($M_w$=8) is a major intermediate-depth earthquake that occurred close to the eastern edge of the Nazca slab flat area. We analyze its rupture process using high-frequency back-projection and seismo-geodetic broadband inversion. The latter approach shows that the earthquake propagated with almost purely normal faulting along the 60° eastward dipping plane. Both imaging techniques provide a very consistent image of the peculiar space-time rupture process of this earthquake : its 60-second long rupture is characterized both by a main northward propagation (resulting in a rupture extent of almost 200km in this direction) and by a reactivation phase of the hypocentral area, particularly active 35s to 50s after

<footnote>ªUniversité Paris Cité, Institut de physique du globe de Paris, CNRS, 1 rue
Jussieu, Paris, 75005, France</footnote>

<footnote>ᵇEarth, Planetary and Space Sciences, University of California, Los Angeles, CA
90025, USA</footnote>

<footnote>ᶜInstituto Geofísico del Perú, Lima, Peru</footnote>

<footnote>ᵈGeoazur, Université Côte d'Azur, IRD, CNRS, Observatoire de la Côte d'Azur, 250 rue
Albert Einstein, Valbonne, 06560, France</footnote>

<footnote>ᵉInstituto Geofísico, Escuela Politécnica Nacional, Quito, Ecuador</footnote>

<footnote>fSorbonne Université, CNRS, Institut des Sciences de la Terre Paris, ISTE P UMR
7193, 4 Place Jussieu, Paris, 75005, France</footnote>

<footnote n="9">⁹Instituto Geográfico Nacional, Lima, Peru</footnote>
origin time. Given the depth of this earthquake (125-140km), the reactivation
time window coincides with the arrival time of the surface-reflected elastic
wavefield. Computed values of the dynamic Coulomb stresses associated
with this wavefield are of the order of ten to several tens of kPa, in a range of
values where dynamic triggering has already been observed. The reactivation
phase of the Peru earthquake may thus originate from fault areas that were
brought close to rupture by the initial rupture front before being triggered by
stress increments provided by the reflected wavefield. Source time function
complexity observed for other large intermediate-depth earthquakes further
suggests that such a mechanism is not an isolated case.

Keywords: rupture propagation, dynamic triggering, array analysis,
seismic source inversion, intermediate depth earthquake, Nazca slab

# 1. Introduction

Earthquake rupture propagation results from the evolving stress field gen-
erated by the fault areas which already slipped. This stress field, character-
ized by large stress concentrations ahead of the rupture front that govern
the rupture propagation regimes (e.g. Andrews (1976); Dunham (2007)), is
modulated by structural heterogeneities at the fault or within the surround-
ing medium. For example, fault damaged zones trap waves that alter the
shear and normal stresses on the fault (Harris and Day, 1997; Huang and
Ampuero, 2011). At a broader scale, Earth stratification and the free surface
play an important role in the stress evolution and thus the rupture devel-
opment (Kaneko and Lapusta, 2010; Kaneko et al., 2011). The effect of the
free surface has therefore been considered since early studies of earthquake
dynamic rupture (Burridge and Halliday, 1971; Archuleta and Frazier, 1978).

The most prominent effects occur for dip-slip mechanisms, where both normal and tangential stresses are perturbed by the free surface (Nielsen, 1998; Oglesby et al., 2000).

Free surface effects have been well explored theoretically and numerically, but identifying their signature on the rupture process of real earthquakes is observationally difficult. Indeed, for large shallow earthquakes, the surface-reflected dynamic stress field overlaps in time and space with the direct one. As a consequence, it is challenging to attribute specific features of the rupture process to free surface effects. Some expected characteristics, such as re-rupturing or rupture jumps ahead of the rupture front (Nielsen, 1998), occur at times very close to the main rupture activation, resulting in an ambiguous imaging. This consideration motivates the analysis of large earthquakes occurring at depth, as the reflected stress field, even if of smaller amplitude, is much better separated in time from the direct stress field. For intermediate-depth events (70-300 km), the surface-reflected stress field reaches again the fault several tens of seconds after rupture initiation, offering the opportunity to better decipher its role on the complexity of the rupture process.

The 2019/05/26 ($M_w = 8.0$) North Peru earthquake is a large intermediate-depth earthquake (Figure 1), the largest one recorded since 1976 according to the Global CMT catalog (Ekström et al., 2012). Based on routine point source analyses (Ekström et al., 2012; Vallée and Douet, 2016), its depth is in the range 125-140 km, consistent with a location inside the Nazca slab, and it activated a North-South striking normal fault, dipping about 60° to the
East or 30° to the West. Such mechanisms are usually understood as a result of the negative buoyancy of the slab (slab pull effect): for slabs that do not extend all the way down to the bottom of the upper mantle, extension in the dip direction is expected at intermediate depth (Isacks and Molnar, 1971). In the case of northern Peru, an additional interpretation of the earthquake origin is related to its location at the eastern edge of the flat section of the Nazca slab (Figure 1), where the slab starts to dip more steeply (Liu and Yao, 2020; Jiménez et al., 2021). Bending stresses acting there, in a similar way as in the outer rise, should favor normal faulting in the upper side of the slab (e.g. Sandiford et al. (2020) and references therein).

Besides its unusual magnitude, the 2019 Peru earthquake also has an un-
common magnitude-duration relationship. Its duration of about 60 s (Fig-
ure 1, Liu and Yao (2020); Tavera et al. (2021); Ye et al. (2020)) would
be usual for a shallow subduction interplate earthquake, but is long for an
intermediate-depth earthquake (Hu et al., 2021; Persh and Houston, 2004;
Vallée, 2013). This feature, which can be hypothesized as a clue of a possible
free surface effect, further motivates a detailed analysis of this earthquake.
In the following, we first conduct high-frequency and broadband studies of
the rupture process, based on extensive seismological and geodetic data sets
(Figure 1). Previous studies of the 2019 Peru earthquake were only based on
telesismic data (Liu and Yao, 2020; Ye et al., 2020; Hu et al., 2021; Jiménez
et al., 2021) with the exception of the study of Tavera et al. (2021) who
used local and regional seismic records to determine the first order charac-
teristics of the rupture process. Our obtained images of the rupture process,
both at high and low frequencies, corroborate the teleseismic results of Hu
et al. (2021): the 2019 Peru earthquake ruptured the steep eastward dipping plane with a dominant northward propagation, but also involved a late reactivation of the hypocentral area. We show that both the actual fault plane and the delayed rupture reactivation are well resolved by the combined use of GNSS, InSAR, and seismic records across a wide range of distances and frequencies. We then evaluate the possibility that surface-reflected dynamic stresses triggered the late rupture reactivation that led to the complex rupture propagation of this earthquake.

## 2. High-frequency imaging

### 2.1. Data and method

We performed back-projection (BP) using teleseismic recordings to image the rupture of the 2019 North Peru earthquake. We selected the broadband seismograms located in Alaska and Europe (Figure 2), as these two areas are densely instrumented and in the appropriate teleseismic range (30-90° from the epicenter). We applied to these two arrays the Multitaper-MUSIC technique, which can resolve closely spaced simultaneous sources, and used the “reference window” strategy to eliminate “swimming” artifacts (Meng et al., 2011b, 2012). The vertical component of the direct P-wave phases is filtered between 1 and 4 Hz and a sliding window of 6 s duration is selected. BP provides source locations relative to the hypocenter, for which we adopt the NEIC determination (latitude=-5.812°N, longitude=-75.27°E, depth=122.6km). Due to the limited depth sensitivity of the method, the BP is performed at the fixed depth of 122.6 km. Supplemental material (Figures S1 and S2) shows that the use of other assumptions (e.g. imposing that the
BP emissions occur on the fault plane and not at a fixed depth) does not affect the main BP features.

## 2.2. Space-time organization of high-frequency radiation

For the Alaska array, Figure 3a) shows that the rupture propagates along the strike of the slab bilaterally. As the along-dip extension is small, the rupture propagation features can be analyzed in a time-distance diagram (Figure 3b). This shows that the rupture initially propagates northwards, during the first 25 s, over a length of ~70 km. Then, the northward front continues its propagation, but at the same time, rupture is also active in the hypocentral area and South of it. This secondary rupture appears to emerge close to the main northward rupture front and to propagate mostly southward; it could thus be interpreted as a branching back-propagating front (Hu et al., 2021). But as this secondary rupture is both less coherent and not clearly connected with a specific point of the main rupture front, this episode is here more generally referred to as a delayed rupture reactivation. We further show in the Supplemental Material (Figure S3) that, when using a synthetic complex rupture, the BP method adequately images the reactivation of the hypocentral area.

The reactivated rupture lasts for ~25 s and involves a 100 km-long segment. The northward rupture lasts for 35 s more and ruptures an additional length of 125 km. The total rupture, imaged with a duration of 60 s, is 200 km to the north and 50 km to the south of the epicenter. Careful analysis of the time-space radiation observed after 60 s (Figure S4) shows that rupture is unlikely to last much longer. As a matter of fact, the late radiators are well explained by the expected rupture replication by the *pP* depth phase,
which offsets the locations and times according to the green arrow of Figure 112
3b. During its whole northward propagation, the speed of the main rupture 113
front is close to 3.5 km/s. Its ratio to shear wave speed (about 80%) is in 114
the usual range for subshear earthquakes. 115

The rupture pattern resolved from the Europe array is similar to the one of 116
the Alaska array. It also shows a main northward rupture front with rupture 117
velocity of 3.5 km/s and a secondary rupture (Figure S6). Nevertheless, the 118
BP using the Europe array shows a less focused image, because of a lower 119
number of stations and strong contamination by the *pP* depth phase. Owing 120
to the radiation pattern, the *pP* phase has indeed a stronger and more stable 121
amplitude than the direct *P* phase. 122

All the previous BP analyses of the Peru earthquake (Ye et al., 2020; 123
Liu and Yao, 2020; Hu et al., 2021) imaged the main northward propagation 124
of the rupture, over a 150-200km length consistent with our determination. 125
The study of Hu et al. (2021) also detected the rupture reactivation in the 126
hypocentral area, but could locate there only one strong radiator 40s af- 127
ter rupture initiation. Our results rather show that the reactivated rupture 128
generated continuous high-frequency emission, in a similar way as the main 129
northward rupture front. 130

### 3. Broadband rupture process

#### 3.1. *Data*

##### 3.1.1. *InSAR*

We processed Sentinel-1 InSAR data from five satellite tracks spanning 134
the earthquake date and covering a broad area around the epicenter (Fig- 135
ure 4). Data were processed using the NSBAS software (Doin et al., 2011; Grandin, 2015). Processing details are provided in the Supplementary Material. In spite of the challenging surface conditions, characterized by dense vegetation cover, coseismic deformation is successfully imaged in both ascending and descending geometries for tracks A120, D069, D171, spanning temporal baselines of 30 days, 6 days and 12 days, respectively. Both descending tracks D069 and D171 cover most of the deformation zone, whereas descending track D098 appears to capture long-wavelength fringe patterns that cannot be related unambiguously to coseismic deformation. We fail to measure meaningful deformation on the remaining ascending track A047 due to low interferometry coherence, which is likely caused by long temporal baseline (30 days) and sub-optimal perpendicular baseline (69 days). Therefore, subsequent analysis is restricted to the three tracks A120, D069, D171.

Deformation is characterized by a broad elliptical pattern, elongated 400 km in the north-south direction, and 200 km in the east-west direction, consistent with the great depth of the earthquake. A maximum line-of-sight (LOS) displacement of ~20 cm is measured near the epicenter, directed away from the satellite, i.e. consistent with subsidence induced above a normal faulting earthquake. Due to the broad footprint of the deformation pattern, interferograms are affected by coseismic displacement on a large fraction of the covered area. Referencing is achieved by projecting displacement vectors derived from GNSS data into the LOS, subtracting the InSAR-derived LOS displacement at the location of the GNSS sites, fitting a bilinear plane to the residual LOS displacement, and finally adding the fitted plane to the interferograms. The resulting interferograms of the 2019 Northern Peru earthquake
(to our knowledge the deepest earthquake characterized by InSAR) are shown
in Figure 4.

3.1.2. GNSS

The GNSS data set includes 137 continuous stations, among which 53
belong to the Instituto Geográfico Nacional of Peru, 21 to the Instituto Ge-
ográfico Militar of Ecuador, 56 to the Instituto Geofísico de la Escuela Po-
litecnica de Quito, Ecuador (Mothes et al., 2018), 1 to the Instituto Geofísico
del Perú, 3 to the Instituto Geográfico Agustín Codazzi of Colombia and 3
from the Geored network from the Servicio Geológico Colombiano (Mora-
Páez et al., 2018). This data set is complemented by 4 survey-mode mea-
surements in the Amazonia basin that improve the coverage in the epicentral
area (LMAS, YRMG, LGN1 and SRM1 sites, Figures S16 and S20).

Both survey-mode and continuous raw GNSS data were simultaneously
analyzed using the GAMIT/GLOBK software (Herring et al., 2010) to pro-
duce daily loosely constrained solutions. In a second step, we expressed the
daily solutions in the International Terrestrial Reference Frame (ITRF2014,
(Altamimi et al., 2016)), using a seven-parameter Helmert transformation
estimated from the sites common to our solution and to global sites from
the International Global Navigation Satellite Systems Service. The average
repeatability for the east, north, and vertical components of daily positions
are 1.5, 1.5, and 5.0 mm, respectively.

For the continuous GNSS sites, co-seismic offsets were computed from
the difference of position averaged one week after and one week before the
earthquake. No significant post-seismic deformation was observed in the time
series. The 4 campaign sites only had one measurements in July 2018 before
the earthquake and were re-occupied in mid-July 2019, two months after the earthquake. For these sites, we took into account the motion predicted for the South America plate in the ITRF2014 frame. However, for Amazonia east of the Andean cordillera, tectonic deformation is expected to be small (< 2 mm/yr) with respect to the stable part of the South America plate (Nocquet et al., 2014; Villegas-Lanza et al., 2016) and can be neglected. Co-seismic offsets were therefore derived from the difference between the average position in July 2019 and July 2018.

The three-dimensional coseismic displacements are shown in Figure 4 and are reported, with their associated uncertainty in Table S2. Given its magnitude and depth, the earthquake caused displacement exceeding measurement uncertainties (1-2 mm) over a broad area extending ~1500 km from southern Colombia to central Peru. The largest displacement occurs at SRM1 (longitude -74.93°E, latitude -4.71°N) with almost 7 cm of horizontal motion and a coseismic subsidence of 15 cm, in agreement with the InSAR measurement.

### 3.1.3. Seismic records

Seismic data from broadband and accelerometer networks well cover the epicentral area (Figure 1). We use the data from the EC network in Ecuador (Alvarado et al., 2018), more specifically selecting the records from the REN-SIG broadband network (when signals are not clipped) and from the RENAC accelerometric network. In Peru, we use the data from the IGP and CISMID-SENCICO strong motion networks. Data from the broadband seismic Brazilian network (BR), which would complement the eastern azimuthal coverage, are clipped for the main shock and cannot be used for waveform modeling. We add to this local data set body-wave teleseismic records (11 for P waves
and 8 for SH waves) recorded at stations from the GSN and Geoscope networks. The records of these global stations were selected to avoid nodal planes and to provide a balanced coverage of azimuths. Credit to the seismic networks used is provided in the first section of the Supplementary Material.

## 3.2. Absolute location of the potential fault planes

Due to the large earthquake depth and its small number of aftershocks (Ye et al., 2020), little a priori information can be used to determine the exact location of the causative fault. Before inverting the whole data set for the space-time rupture process, we thus conduct static inversions using InSAR and GNSS data to determine the absolute fault locations for the two possible planes. Strike, dip and rake of the two possible planes, taken from GCMT (Ekström et al., 2012), are equal to (351°, 57°, -87°) and (166°, 33°, -94°), respectively. The former geometry (steep eastward-dipping plane) is hereafter referred to as FP1 and the latter one as FP2.

To do so, we first perform a Markov Chain Monte Carlo exploration of the latitude, longitude, depth and coseismic slip (taken as a constant over a rectangular area) to assess the uncertainties on the fault location and size. In this approach, strike and dip of both fault planes are kept fixed. Uncertainty on fault location reaches ~5 km in all three directions for both fault planes, whereas fault size (length and width) are uncertain within ~10 km. As a complementary approach, following the technique described in Lauer et al. (2020), we enlarge the fault and invert for spatially-variable slip distribution on a discretized fault plane, systematically exploring the fault location and degree of spatial smoothing to investigate potential trade-offs between fault location and discretization assumptions. Details of the static inversion are
provided in the Supplementary Material. From these static inversions, we conclude that InSAR and GNSS geodetic data well constrain the FP1 and FP2 absolute locations, but the geodetic data cannot discriminate alone the actual fault plane, since both fault planes achieve similar performance in explaining the deformation pattern at the surface.

### 3.3. Forward and inverse problem

We use the NEIC/USGS epicenter (lat=-5.812°N, lon=-75.27°E), and adapt the NEIC/USGS hypocentral depth (122.6km) so that it fits with the fault geometries determined from the static inversion: for FP1 and FP2, the hypocenter is fixed at 129.4 km depth and 121.6 km depth, respectively. FP1 and FP2 are then discretized into 12 km wide (along dip) and 12 km long (along strike) rectangular subfaults. With respect to the hypocenter, FP1 and FP2 models extend 186 km along strike in the NNW direction, 114 km in the SSE direction, and 54 km in both the updip and downdip directions. These values have been selected based on previous studies showing a dominant Northward rupture direction (Liu and Yao, 2020; Ye et al., 2020; Hu et al., 2021; Tavera et al., 2021; Jiménez et al., 2021) and on initial inversion tests, indicating that exploration of a larger rupture extent is unnecessary.

At the local scale, static and dynamic Greens's functions for each of the subfaults are computed with the EDCMP/EDGRN simulation code (Wang et al., 2003) and the discrete wavenumber method of Bouchon (1981), respectively. Teleseismic *P* and *SH* waves computation uses the reciprocity technique of Bouchon (1976) coupled with the reflectivity method (Müller, 1985). All these Green's functions are computed in the one-dimensional Crust1 structure model (Laske et al., 2013), in which the thin superficial
sedimentary layer (500m thickness) has been removed (Table S4). This simplified propagation model imposes a high-frequency limit to the waveform modeling, leading to the following selected frequency ranges: *P* and *SH* teleseismic displacements are filtered between 0.005 Hz and 0.125 Hz; at the local scale, the closest displacement records (BAGU, CHCA, IQU, JUAJ, RIOJ, TOCA) are filtered between 0.02 Hz and 0.1 Hz, and all the other data between 0.02 Hz and 0.05 Hz.

Kinematic analysis of the earthquake follows the approaches described in Delouis et al. (2002) and Grandin et al. (2015). The forward rupture model considers a subfault source time function built with two overlapping triangles of duration equal to 4 s, allowing each subfault to slip for a maximum of 6 s. For each of the 225 subfaults, the parameters controlling the space-time rupture evolution are the amplitudes of the two triangles, the rake, and the rupture onset time. This latter parameter makes the inverse problem non-linear, and a simulated annealing algorithm is used to converge toward the optimal model. The misfit function to be minimized is the sum, equally weighted, of the normalized *rms* errors for the four data sets (InSAR, GNSS, local seismic data, teleseismic data). A minimization constraint on the seismic moment, as well as modest smoothing constraints on rake and final slip, are enforced to prevent spurious slip, in particular on the fault borders.

## 3.4. *Exploration of rupture scenarios*

We first explore how our seismo-geodetic data set is able to discriminate between the two possible fault planes, FP1 and FP2. To do so, inversions are made for the two geometries with the same configuration. Rake is allowed to vary at ±30° compared to the GCMT value, and the only constraint
on the triggering times is based on a maximum average rupture velocity (with respect to hypocenter) equal to the *S* wave velocity of 4.5 km/s. As a result, subfaults located close to the hypocenter can slip at any time during the rupture, while a subfault e.g. 72 km away from the hypocenter has a minimum triggering time of 16 s. This loose constraint accounts for the fact that the earthquake propagation is dominantly along strike, in mode III, thus with rupture velocity limited by the *S* wave velocity. The model does not allow for the reactivation of the same subfault at different times; the formal inclusion of this possibility in the inversion would require a too large number of time windows, resulting in poorly constrained results. However, given the fine grid spacing of 12 km considered in the fault model, re-rupture is a valid interpretation when neighboring subfaults slip at very different times.

Figure 5 summarizes the misfits obtained for the four data sets after global optimization for the FP1 and FP2 cases. As expected from section 3.2, the geodetic misfits (InSAR and GNSS) can be made very similar for FP1 and FP2, and do not alone discriminate the fault planes. However, FP2 misfits are significantly worse for both the teleseismic and local seismic data sets, as further illustrated by waveforms comparison at all stations and components (Figures S21, S22, S24 and S25). Some individual stations, such as the key local station IQU (the only one located East of the rupture), are also definitely better modeled by the FP1 scenario.

The rupture evolution for FP1 is shown in Figure 6. In agreement with high-frequency imaging (section 2) and with the study of Hu et al. (2021), late activation of the hypocentral region is observed. This reactivation is particularly clear 36 s to 48 s after rupture initiation (Figure 6d) and appears
to last up to 60 s (see red filled STF in Figure 6h). In terms of seismic moment, the hypocentral area is activated in a similar way by its initial phase (up to 20s) and by its reactivation: both episodes carry a seismic moment of ~1.6 10²⁰ Nm, equivalent to $M_w$=7.4. Besides this important feature, most other aspects of the rupture process are consistent with teleseismic studies enforcing outward rupture propagation (Liu and Yao, 2020; Ye et al., 2020). The initial bilateral propagation during the first 10 s was followed by an episode of relatively low slip for about 15 s. Rupture then progressed northward with higher slip (up to 1.6 m), breaking in 30 s a 100 km-long segment of the seismic fault. The northern termination is located 170 km away from the hypocenter and the total seismic moment is found equal to 1.4 10²¹ Nm/s ($M_w$=8.03).

We finally show that the late reactivation is required by the data by constraining the minimum rupture velocity to be 2.5 km/s. In this scenario, 30 s or later after origin time, rupture occurs at least 75 km away from the hypocenter, thus excluding a delayed rupture in the hypocentral area. As for the FP2 scenario, misfits are similar for geodetic data, but significantly worse for the seismic data set (blue misfit bars in Figure 5). In particular, stations South of the earthquake (TOCA and JUAJ), which are expected to capture rupture complexities during the northward propagation, are not well modelled (see Figures S21 and S23).
332 **4. Possible mechanisms for delayed reactivation**

333 **4.1. *Dynamic origins without free-surface interaction***

334 Both high frequency imaging and broadband process inversion show that
335 the rupture history of the 2019 Peru earthquake does not simply consist in a
336 steady Northward propagation. It also involves a delayed reactivation of the
337 hypocentral area, which is particularly clear 35s to 50s after origin time.

338 Delayed ruptures have been reported for some past earthquakes, including
339 the 1984 Morgan Hill earthquake (Beroza and Spudich, 1988), the 1987 Su-
340 perstition Hills earthquake (Wald et al., 1990), the 1999 Chi-Chi earthquake
341 (Lee et al., 2006), the 2001 and 2007 Peru earthquakes (Lay et al., 2010), the
342 2008 Wenchuan earthquake (Zhang et al., 2012), the 2010 El Mayor - Cuca-
343 pah earthquake (Meng et al., 2011a) and the 2011 Tohoku earthquake (Lee
344 et al., 2011). This rupture pattern occurs in dynamic rupture simulations due
345 to a variety of mechanisms. One mechanism is the effect of heterogeneities of
346 fault strength and stress (e.g. Goto et al. (2012)); for instance, bilateral sec-
347 ondary rupture fronts can emerge from the interaction between a rupture and
348 the stress concentrations left by a previous earthquake (Kame and Uchida,
349 2008). Galvez et al. (2016) show that slip reactivation and back-propagating
350 fronts can also be triggered if the fault strength undergoes a second weak-
351 ening phase at large slip, for instance due to thermally activated weakening
352 processes. Dunham (2005) points out how the existence of interface waves
353 can be at the origin of rupture complexities.

354 Another proposed mechanism involves heterogeneities of the materials
355 that surround the fault. Idini and Ampuero (2020) and Huang and Ampuero
356 (2011) show that the presence of a damaged zone, modeled as a low rigidity
layer around the fault gives birth to the coexistence of crack-like and pulse-
like ruptures that involve multiple back-propagating rupture fronts. Finally,
even under uniform initial stress and frictional properties, hypocentral area
reactivation can occur due to the stress growing behind an outward propa-
gating pulse (Nielsen and Madariaga, 2003; Gabriel et al., 2012).

4.2. *Free-surface reflected dynamic stresses*

362

363

364

365

366

367

368

369

370

371

372

373

374

375

376

377

378

379

380

381

Reactivation of the hypocentral area occurred at a time very close to the
arrival time of the surface-reflected pP waves, i.e. at the two-way P travel
time from the hypocenter to the free surface. Given the hypocentral depth
(129.4 km) and the selected structure model based on Crust1 (Laske et al.
(2013), Table S4), the pP arrival time at the hypocenter is 34.6s. Owing to
the focal mechanism with a P axis close to the vertical, we expect the pP
wave to have a larger amplitude than the sS wave, and even larger than the
sP and pS waves because of the near-vertical incidence angle at the surface.
More specifically, because the pP wave induces vertical compression when
reflected from the free surface, its effect is to increase the shear stress τ for
a normal fault. On the other hand, the induced normal stress σ (taken as
positive for dilatation) is negative, but its absolute value is small for a steeply
dipping plane. Simple calculation in an homogeneous half-space for a vertical
pP wave and a 60° dipping normal fault plane (similar to the FP1 geometry)
shows that the shear stress amplitude is √3 larger than the normal stress
amplitude. In terms of Coulomb stress range, reactivation is then expected
to be favored for the FP1 scenario, whatever the friction coefficient μ. We
note that for a 30° dipping normal fault plane (similar to the FP2 scenario),
such reactivation would be favored only for low values of μ, as the normal
stress amplitude is $\sqrt{3}$ larger than the shear stress amplitude.

In order to quantitatively evaluate the surface-reflected dynamic stresses associated with the rupture process of the 2019 Peru earthquake, we use the QSSP program (Wang et al., 2017). This method can provide the full dynamic stress field, as well as the dynamic stress field without the free surface contribution, generated in a spherical Earth model by a source process described by an arbitrary number of point sources. We use as input the optimal FP1 scenario (Figure 6) and the crustal structure used in the seismogeodetic inversion (Table S4). Simulations are then successively run with and without the free surface contribution, in order to compute the stress field on the fault plane in both cases. The reflected stress field is then obtained by subtraction and the shear ($\tau$) and normal ($\sigma$) reflected stresses are calculated based on the FP1 mechanism.

Figure 7 shows the reflected dynamic Coulomb stress ($= \tau + \mu\sigma$, with $\mu$ here taken equal to 0.4) at one location of the hypocentral area. This location, indicated by the white circle in the snapshots of Figure 6, is 18 km along strike and 6 km updip from the hypocenter. Due to the large distance to the free surface, the stress field varies spatially slowly, and the Coulomb stress evolution shown is representative of all points located at similar distances, or closer, from the hypocenter. Coulomb stress takes positive values, of the order of 10-20 kPa at the time of hypocentral reactivation. As expected, these values are slightly lower for a larger $\mu$ (and larger for a smaller $\mu$) because $\tau$ and $\sigma$ tend to be of opposite signs. However, even for $\mu = 0.8$, values of 10 kPa are reached. The computed peak stress values most likely provide lower bounds of the actual peak stress, because high frequencies (> 0.2 Hz) have
407 not been inverted for in the source model.

408 Slip reactivation by surface-reflected waves can be understood as a case
409 of remote, instantaneous dynamic triggering, a phenomenon that has been
410 observed on the wake of several past earthquakes. Instantaneous dynamic
411 triggering occurs most often during the passage of surface waves (see Freed
412 (2005), Brodsky and van der Elst (2014) and Hill and Prejean (2015) for
413 review papers), but it has been also reported for intermediate-depth earth-
414 quakes during the passage of body waves (Luo and Wiens, 2020). Here, the
415 triggering distance is of at least 260 km, twice the hypocenter depth, which
416 is well within the distance range of observed dynamic triggering (up to thou-
417 sands of km). The dynamic Coulomb stresses estimated here are of the same
418 order of magnitude as the stresses reported to have dynamically triggered
419 earthquakes in the Coso geothermal field (California), during the passage of
420 surface waves of the 2002 Denali, Alaska earthquake (Prejean et al., 2004).
421 The corresponding dynamic strains are in the low range of values observed
422 to induce substantial seismicity rate increases in California (see Figure 4 of
423 Miyazawa et al. (2021)).

424 These stresses are however very small when compared to typical earth-
425 quake stress drops. Thus, rupture triggering through a small Coulomb stress
426 increment requires some areas of the fault to have already been brought very
427 close to failure when the triggering waves arrive. This requirement is ex-
428 pected to be satisfied anywhere along the edges of a rupture, where stresses
429 tend to concentrate. A plausible scenario for the Peru earthquake is that
430 its hypocentral reactivation started from such a critically stressed location,
431 where the surface-reflected stress field induced a secondary rupture.
4.3. *Correlations between source complexities and surface-reflected waves*

The hypothesis of slip reactivation by the *pP* waves for the 2019 Peru earthquake motivates the search for a similar mechanism for other complex past earthquakes. Although this would ideally require to know the space-time source process for a number of events, first insights can be gained by the observation of source time functions (STFs). We here consider the STFs of large earthquakes (*M*<sub>w</sub> > 7) in the depth range between 30km and 150km. Most of these STFs have simple shapes, in which case the role of surface-reflected phases cannot be easily discriminated, even when the STFs have long durations. We therefore focus on complex earthquakes, based on the criterion that their STFs do not grow monotonically toward their peak (see legend of Figure 8). We found 20 earthquakes in the SCARDEC database (Vallée and Douet, 2016) satisfying this criterion. Figure 8 shows their STFs along with the arrival times of the *pP* and *sS* surface-reflected phases (computed as the vertical two-way travel times to the surface, using the SCARDEC earthquake depths of Table S5).

The non-monotonic criterion tends to select long-duration earthquakes, which increases the likelihood that surface-reflected phases fortuitously arrive during earthquake rupture. We indeed observe that 13 of the selected STFs have durations longer than 30s. However, when looking at the SCARDEC database for the same magnitude range and non-monotonic criterion, a first interesting element is that no earthquakes with depth larger than 150km reach such a duration. This indicates that a mechanism responsible for sustained source emissions is present for earthquakes shallower than 150km. Such a mechanism may be due to intrinsic differences between shallow,
intermediate-depth and deep earthquakes (e.g. Houston (2015)), but the fact that surface-reflected waves arrive during the rupture process for 16 of the 20 earthquakes of Figure 8 offers the appealing interpretation that these waves played a role in the source complexity. In particular, besides the 2019 North Peru earthquake (index 20 in Figure 8), the 1995/08/16 Solomon (index 3), the 2001/08/21 New Zealand (index 8) and the 2006/05/03 Tonga (index 9) earthquakes have a major peak in their STFs occurring just after the *pP* arrival time. The three latter earthquakes, that occur at depths of 45-60km, further have an inverse mechanism, in which case both the shear and the normal reflected *pP* stresses favor rupture. Other earthquakes show a less obvious time correlation between surface-reflected phases and rupture complexities. However, this does not rule out the role of surface-reflected stresses. Indeed, a *pP* (or *sS*) wave which would be at the origin of an updip rupture complexity can have a shorter travel time than the one indicated. On the contrary, a *pP* (or *sS*) wave which would be at the origin of a downdip or an along-strike rupture complexity has a longer travel time than the one indicated.

The detailed case of the 2019 Peru earthquake and these STFs observations suggest that for shallow earthquakes, surface-reflected stresses should play an important role as their amplitudes are much larger in this case. Although rupture complexities can be caused by a number of alternative mechanisms (see section 4.1), surface-reflected stresses effects are likely not restricted to the cases where the fault intersects the free surface (Nielsen, 1998; Oglesby et al., 2000; Zhang and Chen, 2006).
5. Conclusion

We analyzed the source process of the 2019 Peru intermediate-depth earthquake using high-frequency back-projection and seismo-geodetic broadband kinematic inversion. Both techniques provide very consistent pictures of the rupture propagation. The 60s long rupture is characterized by a main ∼200 km Northward propagation and by a reactivation phase of the hypocentral area, particularly active 35s to 50s after origin time. Broadband inversion also favors the activation (with an almost purely normal mechanism) of the 60° Eastward-dipping plane rather than the 30° Westward dipping plane.

Reactivation phases have been previously reported for a number of shallow earthquakes, but the 2019 Peru earthquake provides an example of this phenomenon for an intermediate depth earthquake. A striking observation is that the initiation of the reactivation phase occurs at a time very close to the first *pP* surface-reflected wave arrival time. Using our kinematic inversion results, a simulation of the surface-reflected wavefield predicts that dynamic Coulomb stresses in the hypocentral area are of the order of 10-20kPa, and are certainly even larger at frequencies higher than 0.2 Hz, that are not considered in our calculation. Stress perturbations of the order of 10 kPa were previously shown to be sufficient for dynamic triggering. Thus, the observed reactivation phase may have occurred at fault areas that were brought close to the rupture by the initial front and were then triggered by the arrival of the reflected wavefield.

Source time functions of other large intermediate depth earthquakes worldwide provide clues of a similar mechanism, with rupture complexities correlated in time with the arrival of the surface-reflected waves. Such a mech-
506

anism has consequences for shallow earthquakes, for which the larger am-
plitude surface-reflected stresses have a larger potential to trigger secondary
ruptures.

507

508

509

Acknowledgments

510

511

512

513

514

515

516

517

518

519

520

521

522

523

524

Constructive reviews provided by two anonymous reviewers helped us to
improve the original manuscript. We thank Rongjiang Wang for his advises
on the use of the QSSP program. We thank Mario Mendoza for discussions
about the GNSS data in Peru. We are grateful to the Instituto Geografico
Nacional and the Instituto Geografico Militar for providing the GNSS data
in Peru and in Ecuador respectively. We thank David Cisneros for helping
with the IGM GNSS data. We are in dept to the people who went to the field
for the GNSS measurements in Amazonia. Funding from BQR OCA Andes-
width project is acknowledged, as well as the ANR through the S5 project
(ANR-19-CE31-0003-01). P. Jarrin was supported by the Secretaría Nacional
de Educación Superior, Ciencia y Tecnología (SENESCYT-Ecuador) through
the doctoral fellowship IFTH-DFN-2018-0096/092-2017. We acknowledge
logistic support from the Institut de Recherche pour le Development (IRD).
Numerical computations were partly performed on the S-CAPAD/DANTE
platform, IPGP, France.

525

Data availability

526

527

Seismic data are publicly available and mostly belong to the Federation of
Digital Seismograph networks (FDSN). Credit to the corresponding seismic
528 networks and to the data centers distributing their data is provided in the
529 first section of the Supplementary Material.

530 References

531 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, 6109–6131.

535 Alvarado, A., Ruiz, M., Mothes, P., Yepes, H., Segovia, M., Vaca, M., Ramos, C., Enríquez, W., Ponce, G., Jarrín, P., et al., 2018. Seismic, volcanic, and geodetic networks in ecuador: Building capacity for monitoring and research. Seismological Research Letters 89, 432–439.

539 Andrews, D., 1976. Rupture velocity of plane strain shear cracks. Journal of Geophysical Research 81, 5679–5687.

541 Archuleta, R.J., Frazier, G.A., 1978. Three-dimensional numerical simulations of dynamic faulting in a half-space. Bulletin of the Seismological Society of America 68, 541–572.

544 Beroza, G.C., Spudich, P., 1988. Linearized inversion for fault rupture behavior: Application to the 1984 morgan hill, california, earthquake. Journal of Geophysical Research: Solid Earth 93, 6275–6296.

547 Bouchon, M., 1976. Teleseismic body wave radiation from a seismic source in a layered medium. Geophysical Journal International 47, 515–530.
549

Bouchon, M., 1981. A simple method to calculate green's functions for elastic layered media. Bulletin of the Seismological Society of America 71, 959-971.

550

551

552

Brodsky, E.E., van der Elst, N.J., 2014. The uses of dynamic earthquake triggering. Annual Review of Earth and Planetary Sciences 42, 317-339.

553

554

Burridge, R., Halliday, G., 1971. Dynamic shear cracks with friction as models for shallow focus earthquakes. Geophysical Journal International 25, 261-283.

555

556

557

Delouis, B., Giardini, D., Lundgren, P., Salichon, J., 2002. Joint inversion of insar, gps, teleseismic, and strong-motion data for the spatial and temporal distribution of earthquake slip: Application to the 1999 izmit mainshock. Bulletin of the Seismological Society of America 92, 278-299.

558

559

560

561

Doin, M.P., Guillaso, S., Jolivet, R., Lasserre, C., Lodge, F., Ducret, G., Grandin, R., 2011. Presentation of the small baseline nsbas processing chain on a case example: the etna deformation monitoring from 2003 to 2010 using envisat data, in: Proceedings of the Fringe symposium, ESA SP-697, Frascati, Italy. pp. 3434-3437.

562

563

564

565

566

Dunham, E.M., 2005. Dissipative interface waves and the transient response of a three-dimensional sliding interface with coulomb friction. Journal of the Mechanics and Physics of Solids 53, 327-357.

567

568

569

Dunham, E.M., 2007. Conditions governing the occurrence of supershear ruptures under slip-weakening friction. Journal of Geophysical Research: Solid Earth 112.

570

571
Ekström, G., Nettles, M., Dziewoński, A., 2012. The global cmt project 2004–2010: Centroid-moment tensors for 13,017 earthquakes. Physics of the Earth and Planetary Interiors 200, 1–9.

Freed, A.M., 2005. Earthquake triggering by static, dynamic, and postseismic stress transfer. Annu. Rev. Earth Planet. Sci. 33, 335–367.

Gabriel, A.A., Ampuero, J.P., Dalguer, L.A., Mai, P.M., 2012. The transition of dynamic rupture styles in elastic media under velocity-weakening friction. Journal of Geophysical Research: Solid Earth 117.

Galvez, P., Dalguer, L.A., Ampuero, J.P., Giardini, D., 2016. Rupture reactivation during the 2011 m w 9.0 tohoku earthquake: Dynamic rupture and ground-motion simulations. Bulletin of the Seismological Society of America 106, 819–831.

Goto, H., Yamamoto, Y., Kita, S., 2012. Dynamic rupture simulation of the 2011 off the pacific coast of tohoku earthquake: Multi-event generation within dozens of seconds. Earth, planets and space 64, 1167–1175.

Grandin, R., 2015. Interferometric processing of slc sentinel-1 tops data, in: FRINGE'15: Advances in the Science and Applications of SAR Interferometry and Sentinel-1 InSAR Workshop, Frascati, Italy, 23-27 March 2015.

Grandin, R., Vallée, M., Satriano, C., Lacassin, R., Klinger, Y., Simoes, M., Bollinger, L., 2015. Rupture process of the mw= 7.9 2015 gorkha earthquake (nepal): Insights into himalayan megathrust segmentation. Geophysical Research Letters 42, 8373–8382.
595

Harris, R.A., Day, S.M., 1997. Effects of a low-velocity zone on a dynamic rupture. Bulletin of the Seismological Society of America 87, 1267-1280.

596

597

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, 58-61.

598

599

600

Herring, T., King, R., McClusky, S., et al., 2010. Introduction to gamit/globk. Massachusetts Institute of Technology, Cambridge, Massachusetts.

601

602

603

Hill, D., Prejean, S., 2015. 4.11 - dynamic triggering, in: Schubert, G. (Ed.), Treatise on Geophysics (Second Edition). second edition ed. Elsevier, Oxford, pp. 273-304.

604

605

606

Houston, H., 2015. 4.13 - deep earthquakes, in: Schubert, G. (Ed.), Treatise on Geophysics (Second Edition). second edition ed. Elsevier, Oxford, pp. 329-354.

607

608

609

Hu, Y., Yagi, Y., Okuwaki, R., Shimizu, K., 2021. Back-propagating rupture evolution within a curved slab during the 2019 m w 8.0 peru intraslab earthquake. Geophysical Journal International 227, 1602-1611.

610

611

612

Huang, Y., Ampuero, J.P., 2011. Pulse-like ruptures induced by low-velocity fault zones. Journal of Geophysical Research: Solid Earth 116.

613

614

Idini, B., Ampuero, J.P., 2020. Fault-zone damage promotes pulse-like rupture and back-propagating fronts via quasi-static effects. Geophysical Research Letters 47, e2020GL090736.

615

616
617

Isacks, B., Molnar, P., 1971. Distribution of stresses in the descending lithosphere from a global survey of focal-mechanism solutions of mantle earthquakes. Reviews of Geophysics 9, 103-174.

618

619

620

Jiménez, C., Luna, N., Moreno, N., Saavedra, M., 2021. Seismic source characteristics of the intermediate-depth and intraslab 2019 northern peru earthquake (mw 8.0). Journal of Seismology 25, 863-874.

621

622

623

Kame, N., Uchida, K., 2008. Seismic radiation from dynamic coalescence, and the reconstruction of dynamic source parameters on a planar fault. Geophysical Journal International 174, 696-706.

624

625

626

Kaneko, Y., Ampuero, J.P., Lapusta, N., 2011. Spectral-element simulations of long-term fault slip: Effect of low-rigidity layers on earthquake-cycle dynamics. Journal of Geophysical Research: Solid Earth 116.

627

628

629

Kaneko, Y., Lapusta, N., 2010. Supershear transition due to a free surface in 3-d simulations of spontaneous dynamic rupture on vertical strike-slip faults. Tectonophysics 493, 272-284.

630

631

632

Laske, G., Masters, G., Ma, Z., Pasyanos, M., 2013. Update on crust1.0—a 1-degree global model of earth's crust, in: Geophys. res. abstr, p. 2658.

633

634

Lauer, B., Grandin, R., Klinger, Y., 2020. Fault geometry and slip distribution of the 2013 mw 7.7 balochistan earthquake from inversions of sar and optical data. Journal of Geophysical Research: Solid Earth 125, e2019JB018380.

635

636

637

638

Lay, T., Ammon, C.J., Hutko, A.R., Kanamori, H., 2010. Effects of kinematic constraints on teleseismic finite-source rupture inversions: Great

639
640

641

peruvian earthquakes of 23 june 2001 and 15 august 2007. Bulletin of the
Seismological Society of America 100, 969-994.

642

643

644

Lee, S.J., Huang, B.S., Ando, M., Chiu, H.C., Wang, J.H., 2011. Evidence of large scale repeating slip during the 2011 tohoku-oki earthquake. *Geophysical Research Letters* 38.

645

646

647

Lee, S.J., Ma, K.F., Chen, H.W., 2006. Three-dimensional dense strong motion waveform inversion for the rupture process of the 1999 chi-chi, taiwan, earthquake. *Journal of Geophysical Research: Solid Earth* 111.

648

649

650

Liu, W., Yao, H., 2020. Rupture process of the 26 may 2019 mw 8.0 north-
ern peru intermediate-depth earthquake and insights into its mechanism.
Geophysical Research Letters 47, e2020GL087167.

651

652

653

Luo, Y., Wiens, D.A., 2020. High rates of deep earthquake dynamic triggering in the thermal halos of subducting slabs. *Geophysical Research Letters* 47, e2019GL086125.

654

655

656

Meng, L., Ampuero, J., Page, M., Hudnut, K., 2011a. Seismological evidence
and dynamic model of reverse rupture propagation during the 2010 m7. 2
el mayor cucapah earthquake, in: Agu fall meeting abstracts, pp. S52B-04.

657

658

659

660

Meng, L., Ampuero, J.P., Sladen, A., Rendon, H., 2012. High-resolution backprojection at regional distance: Application to the haiti m7. 0 earthquake and comparisons with finite source studies. *Journal of Geophysical Research: Solid Earth* 117.

661

Meng, L., Inbal, A., Ampuero, J.P., 2011b. A window into the complexity of
the dynamic rupture of the 2011 mw 9 tohoku-oki earthquake. Geophysical Research Letters 38.

Miyazawa, M., Brodsky, E.E., Guo, H., 2021. Dynamic earthquake triggering in southern california in high resolution: Intensity, time decay, and regional variability. AGU Advances 2, e2020AV000309.

Mora-Páez, H., Peláez-Gaviria, J.R., Diederix, H., Bohórquez-Orozco, O., Cardona-Piedrahita, L., Corchuelo-Cuervo, Y., Ramírez-Cadena, J., Díaz-Mila, F., 2018. Space geodesy infrastructure in colombia for geodynamics research. Seismological Research Letters 89, 446-451.

Mothes, P.A., Rolandone, F., Nocquet, J.M., Jarrin, P.A., Alvarado, A.P., Ruiz, M.C., Cisneros, D., Páez, H.M., Segovia, M., 2018. Monitoring the earthquake cycle in the northern andes from the ecuadorian cgps network. Seismological Research Letters 89, 534-541.

Müller, G., 1985. The reflectivity method: a tutorial. Journal of Geophysics 58, 153-174.

Nielsen, S., Madariaga, R., 2003. On the self-healing fracture mode. Bulletin of the Seismological Society of America 93, 2375-2388.

Nielsen, S.B., 1998. Free surface effects on the propagation of dynamic rupture. Geophysical Research Letters 25, 125-128.

Nocquet, J.M., Villegas-Lanza, J.C., Chlieh, M., Mothes, P., Rolandone, F., Jarrin, P., Cisneros, D., Alvarado, A., Audin, L., Bondoux, F., et al., 2014. Motion of continental slivers and creeping subduction in the northern andes. Nature Geoscience 7, 287-291.
685

Oglesby, D.D., Archuleta, R.J., Nielsen, S.B., 2000. Dynamics of dip-slip faulting: Explorations in two dimensions. Journal of Geophysical Research: Solid Earth 105, 13643-13653.

686

687

688

Persh, S.E., Houston, H., 2004. Deep earthquake rupture histories determined by global stacking of broadband p waveforms. Journal of Geophysical Research: Solid Earth 109.

689

690

691

Prejean, S., Hill, D., Brodsky, E., Hough, S., Johnston, M., Malone, S., Oppenheimer, D., Pitt, A., Richards-Dinger, K., 2004. Remotely triggered seismicity on the united states west coast following the m w 7.9 denali fault earthquake. Bulletin of the Seismological Society of America 94, S348-S359.

692

693

694

695

696

Sandiford, D., Moresi, L., Sandiford, M., Farrington, R., Yang, T., 2020. The fingerprints of flexure in slab seismicity. Tectonics 39, e2019TC005894.

697

698

Tavera, H., Delouis, B., Mercado, A., Portugal, D., 2021. Loreto Intermediate Depth Earthquake of 26 May 2019 (Northeast Peru): Source Parameters by Inversion of Local to Regional Waveforms and Strong-Motion Observations. Seismological Research Letters doi:10.1785/0220200459.

699

700

701

702

Vallée, M., 2013. Source time function properties indicate a strain drop independent of earthquake depth and magnitude. Nature communications 4, 1-6.

703

704

705

Vallée, M., Douet, V., 2016. A new database of source time functions (stfs) extracted from the scardec method. Physics of the Earth and Planetary Interiors 257, 149-157.

706

707
708 Villegas-Lanza, J.C., Chlieh, M., Cavalié, O., Tavera, H., Baby, P., Chire-
709 Chira, J., Nocquet, J.M., 2016. Active tectonics of peru: Heterogeneous
710 interseismic coupling along the nazca megathrust, rigid motion of the pe-
711 ruvian sliver, and subandean shortening accommodation. Journal of Geo-
712 physical Research: Solid Earth 121, 7371-7394.

713 Wald, D.J., Helmberger, D.V., Hartzell, S.H., 1990. Rupture process of
714 the 1987 superstition hills earthquake from the inversion of strong-motion
715 data. Bulletin of the Seismological Society of America 80, 1079-1098.

716 Wang, R., Heimann, S., Zhang, Y., Wang, H., Dahm, T., 2017. Complete
717 synthetic seismograms based on a spherical self-gravitating earth model
718 with an atmosphere-ocean-mantle-core structure. Geophysical Journal
719 International 210, 1739-1764.

720 Wang, R., Martin, F.L., Roth, F., 2003. Computation of deformation in-
721 duced by earthquakes in a multi-layered elastic crust—fortran programs
722 edgrn/edcmp. Computers & Geosciences 29, 195-207.

723 Ye, L., Lay, T., Kanamori, H., 2020. Anomalously low aftershock productiv-
724 ity of the 2019 mw 8.0 energetic intermediate-depth faulting beneath peru.
725 Earth and Planetary Science Letters 549, 116528.

726 Zhang, G., Vallée, M., Shan, X., Delouis, B., 2012. Evidence of sudden
727 rupture of a large asperity during the 2008 mw7. 9 wenchuan earthquake
728 based on strong motion analysis. Geophysical Research Letters 39.

729 Zhang, H., Chen, X., 2006. Dynamic rupture on a planar fault in three-
730 dimensional half-space–ii. validations and numerical experiments. Geo-

731 physical Journal International 167, 917–932.
<figure><img src="image_2.png" /><figcaption>Figure 1: Main earthquake characteristics and local observation configuration. Earthquake epicenter from USGS is shown by the red star, and is associated with the GCMT focal mechanism. Local seismic stations (strong motion and broadband) are indicated by blue triangles and named. Squares show the locations of coseismic GNSS observations. The footprints of the three InSAR scenes used in the present study are shown as orange polygons. Slab geometry from Slab2.0 (Hayes et al., 2018) is shown by depth contours every 20 km. The average SCARDEC source time function (STF), using the method of Vallée and Douet (2016), is shown as an inset, in units of 10¹⁹ Nm/s.</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>Figure 2: Teleseismic arrays used for back-projecting the source emissions of the 2019 North Peru earthquake (red star). The Alaska and Europe arrays are made of the broadband stations shown by the yellow and green triangles, respectively. Credit to the seismic networks used is provided in the first section of the Supplementary Material.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Figure 3: (a) Back-projection imaging of the high-frequency radiation of the 2019 Peru earthquake using the Alaska array, color-coded by time. The white dashed lines are contours of the Slab2.0 model (Hayes et al., 2018). (b) Distance and time of the seismic radiators relative to the hypocenter. The vertical axis shows the distances from hypocenter projected along the 352° strike direction. The dashed line shows a reference rupture speed of 3.5 km/s. The green arrow shows the space-time offset of the *pP* phases estimated from the PREM model. The associated effect (*pP* rupture replication) is shown in Figure S4, where a longer source time window is selected.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 4: Geodetic data set covering the 26 May 2019 Peru earthquake. Vectors show the GNSS-derived static displacement (blue: horizontal; red: vertical). Unwrapped Sentinel-1 interferograms acquired on two descending tracks (D069 – 2019/05/23 → 2019/05/29; D171 – 2019/05/18 → 2019/05/30) and one ascending track (A120 – 2019/05/02 → 2019/06/01) are shown in the background. Negative displacement (blue color) corresponds to motion away from the satellite.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 5: Misfits between data and synthetics for the four types of data (InSAR, GNSS, teleseismic and local seismic records) and the three tested scenarios (see main text for their description).</figcaption></figure>
<figure><img src="image_7.png" /></figure>

Figure 6: Broadband rupture process of the 2019 Peru earthquake. Snapshots of the coseismic slip occurring every 12 s are shown in the panels a)-f), with increasing time from top left to bottom right. The final coseismic slip is shown in g). In each of these panels, black dots are the locations of the inverted source parameters, white triangle is the hypocenter, and white circle is the location where the surface-reflected dynamic Coulomb stress is computed in Figure 7. The top right panel h) shows the total source time function (STF, thick black curve) and its decomposition as a function of key areas of the fault: the red-filled domain is the moment rate coming from the hypocentral area (red box in a)), the green-filled domain is the moment rate coming from the main slip area (green box in a)), and the remaining white-filled domain below the STF is the moment rate coming from other areas of the fault. At 40 s, the moment rate originating from the reactivated hypocentral area is as large as the one originating from the main slip area.
<figure><img src="image_8.png" /><figcaption>Figure 7: Surface-reflected dynamic Coulomb stress at the hypocentral area. Coulomb stress is computed 18 km along strike and 6 km updip from the hypocenter, at the location shown by the white circle in Figures 6a-g. The time window of the main hypocentral reactivation (see Figure 6h) is filled in grey.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 8: Source time functions (STFs) of complex large earthquakes ($M_w > 7$) at depths between 30 km and 150 km. All average STFs from the SCARDEC catalog (Vallée and Douet, 2016) in the corresponding depth and magnitude ranges are considered. They are plotted here if they meet the following criterion characterizing non-monotonic growth: noting $F$ the time-sampled STF and $i_m$ the time index of the STF maximum $F_m$, there exists indexes $i$ and $j$ such as $i < j < i_m$ and $F(j) < F(i) - \alpha F_m$. The value $\alpha$ quantifying the significance of an early secondary maximum is taken equal to 0.15. The dashed blue and green lines are the $pP$ and $sS$ times, respectively, computed by the two-way vertical travel times. Each STF is referred by an index and information on the corresponding earthquake is given in Table S5.</figcaption></figure>
Highlights

Self-reactivated rupture during the 2019 $M_w=8$ northern Peru intraslab earthquake

Martin Vallée, Yuging Xie, Raphaël Grandin, Juan Carlos Villegas-Lanza, Jean-Mathieu Nocquet, Sandro Vaca, Lingsen Meng, Jean Paul Ampuero, Patricia Mothes, Paul Jarrin, Ciro Sierra Farfan, Frédérique Rolandone

* The 2019 Peru intraslab earthquake did not propagate steadily in the northward direction

* Hypocentral region was reactivated several tens of seconds after rupture origin time

* Dynamic stresses induced by surface-reflected waves favored such delayed reactivation