Repeating Earthquakes at the Edge of the Afterslip of the 2016 Ecuadorian M W 7.8 Pedernales Earthquake

Repeating earthquakes repeatedly rupture the same seismic asperity and are strongly linked to aseismic slip. Here, we study the repeating aftershocks of the April 16, 2016 M W 7.8 Pedernales earthquake in Ecuador, which generated a large amount of afterslip. Using temporary and permanent stations, we correlate waveforms from a one-year catalog of aftershocks. We sort events with a minimum correlation coefficient of 0.95 into preliminary families, which are then expanded using template-matching to include events from In total, 376 repeaters are classified into 62 families of 4–15 events. They are relocated, first using manual picks, and then using a double difference method. We find repeating earthquakes during the whole period, occurring primarily within large aftershock clusters on the edges of the areas of largest afterslip release. Their recurrence times, shortened by the mainshock, subsequently increase following an Omori-type law, providing a timeframe for the afterslip's deceleration. Although they are linked temporally to the afterslip, repeater-derived estimates of slip differ significantly from GPS-based models. Combined with the fact that repeaters appear more spatially correlated with the afterslip gradient than with the afterslip maxima, we suggest that stress accumulation at the edge of the afterslip may guide repeater behavior. locked, unmoving spots called asperities, that then break intermittently, causing repeating earthquakes. This means that by studying repeating earthquakes, we can measure the slow slip happening around them. In our study, we look for repeaters during the year before and after the Pedernales earthquake. We find repeaters primarily in large earthquake groups at the edge of the areas where afterslip happened. However, the measures of slow slip we calculated using repeaters and the estimates we calculated with GPS do not match. This, and the location of the repeaters around the slow slip areas rather than inside them, means that the repeaters might not be directly related to the slip amount, but to the stress accumulation at the edge of the slow-slip area.

Repeating earthquakes (or repeaters) are families of two or more events, usually of small magnitude, that represent repeated ruptures on a single seismic asperity through time (Ellsworth, 1995;Nadeau et al., 1994;Vidale et al., 1994). They are identified by their nearly identical waveforms and overlapping rupture areas. Over the last 25 years, they have been observed in a large variety of tectonic settings, from strike-slip faults (Nadeau & Johnson, 1998) and normal faults (Duverger et al., 2018), to reverse faults (K. H. Chen et al., 2008) and megathrust faults (Dominguez et al., 2016;Uchida et al., 2003). In all these settings, repeaters are associated with aseismic slip, be it creep, afterslip, or slow slip events (SSEs). The recurrence of these earthquakes suggests that a seismic asperity is being continually loaded, likely by aseismic slip, and breaking at regular intervals (T. Chen & Lapusta, 2009). This makes repeaters ideal tools to probe the properties of a fault and its slip history.
Repeaters are usually quasi-periodic during the interseismic period, with log(Tr) proportional to log(M 0 1/6 ) (where Tr is recurrence time and M 0 is seismic moment). This relationship seems to be universal across tectonic contexts, further suggesting that repeaters depend on tectonic loading rates (K. H. Chen et al., 2007). This makes them useful to estimate creep at depth (K. H. Chen et al., 2008;Nadeau & McEvilly, 1999). Since repeaters likely represent slip on an isolated asperity being loaded by surrounding creep, repeaters measure slip directly on the fault, independently from global positioning system (GPS) surface measurements. They can therefore complement geodetic models, which have their own uncertainties linked to instrument location, fault geometry, model smoothing, and other parameterization.
Afterslip is a transient aseismic process that occurs in response to a large earthquake, usually in conjunction with aftershocks. Both geodetic and repeating earthquake data can be used to estimate afterslip (Uchida, Yui, et al., 2009). They offer complementary constraints on the geographical extent and amplitude of afterslip. Repeaters are an especially powerful tool along subduction margins as they provide additional insights into afterslip occurring in offshore regions. Recurrence times of repeaters drop after a mainshock, then slowly increase following an Omori-type law (Marone et al., 1995;Peng et al., 2005;Schaff et al., 1998;Taira et al., 2009). This study aims to better characterize the relationship between repeaters and afterslip in the context of the Ecuador-Colombia subduction zone by looking at the post-seismic sequence of the 2016 M W 7.8 Pedernales megathrust earthquake in Ecuador.
The Ecuador-Colombia subduction zone is a complex, spatially heterogeneous region that exhibits a wide variety of slip behaviors, from aseismic slip to large megathrust earthquakes. Over the last century, five earthquakes with magnitudes above 7.5 have occurred in the area. The 1906 M W 8.4-8.8 earthquake (Kanamori & McNally, 1982;Yoshimoto et al., 2017) is the largest known earthquake in the region. It is thought to have ruptured a 200-500 km long segment from central Ecuador to southern Colombia (Kelleher, 1972). During the following 70 years, 3 earthquakes of magnitude 7.7-8.2 each ruptured a distinct ∼200 km portion of this segment, starting with the southern segment in 1942, then the middle segment in 1958 and the northern segment in 1979 (Beck & Ruff, 1984;Kanamori & McNally, 1982). These formed a northward sequence that may only have released as little as a fifth of the moment released by the 1906 earthquake, suggesting that the latter broke through not only the three smaller asperities together, but also some adjacent creeping portions of the megathrust (Beck & Ruff, 1984;Kanamori & McNally, 1982).
The April 16, 2016 M W 7.8 Pedernales earthquake was the latest in the along-strike megathrust sequence. It ruptured a highly coupled patch of the subduction interface  corresponding to the approximate rupture area of the 1942 earthquake . Some studies have argued that the moment released by the Pedernales earthquake represents the entirety of the strain accumulated since 1942, meaning each megathrust earthquake that occurs in this area completely resets the fault loading (Ye et al., 2016;Yoshimoto et al., 2017). Others have proposed that the moment released by the 2016 earthquake far exceeds the strain accumulated since 1942 , suggesting the existence of an earthquake supercycle in the Ecuador-Colombia subduction zone. This would explain the relative seismic quiescence before 1906 and the subsequent enhanced seismic hazard, implying that the 1906 earthquake did not fully reset the strain accumulation along the plate boundary .
In addition to megathrust earthquakes, there is a large variety of interacting slip behaviors along the Ecuadorian margin. North of the 2016 rupture zone (Figure 1), the Punta Galera region exhibits low to intermediate coupling that hosts frequent SSEs accompanied by repeating earthquakes and seismic swarms . This periodic unloading via aseismic slip may cause the area to impede rupture propagation . Other areas prone to slip aseismically include La Plata island, where SSEs have been observed concurrently with seismic swarms Segovia et al., 2018;Vallée et al., 2013), the region downdip of the rupture zone where an SSE was detected in 2015 , and upper plate faults near Esmeraldas (Hoskins et al., 2021). Aseismic behavior was also found during the postseismic period of the 2016 earthquake, with several regions experiencing afterslip, imaged by GPS data . During the first month following the Pedernales earthquake, two main patches of afterslip emerged updip of the coseismic rupture zone in the north and south, about 100 km apart   (Figure 1). Both patches remain static in space for the first month, growing in amplitude while their edges remain fixed . The northern patch is situated in a low interseismic coupling area (5%-40%) while the southern patch is centered on a high interseismic coupling region (50%-90%) (Figure 1), yet in both cases the amplitudes, dimensions and behaviors are similar, with large and rapid early slip .
In this unique context, we explore the link between afterslip and repeating seismicity. We first extract repeaters from an existing earthquake catalog (Agurto-Detzel et al., 2019), relocate them, and enhance our detection through template matching. Repeaters are then compared with an existing afterslip model, in order to discuss the spatiotemporal relationship between them. We note that repeating earthquakes seem related to afterslip gradient rather than to afterslip itself. We therefore hypothesize that repeaters are influenced by the stress transferred from the afterslip.

Materials and Methods
To extract repeating earthquakes occurring during the postseismic phase of the Pedernales earthquake, we use both the permanent seismic network in place in Ecuador (Alvarado et al., 2018), and the one-year temporary deployment of land and ocean bottom seismic (OBS) stations deployed in the aftermath of the Pedernales earthquake (Agurto-Detzel et al., 2019;Meltzer et al., 2019). We examine 14 months of postseismic data, as well as 12 months of interseismic data. We start with the one-year long catalog produced by Agurto-Detzel et al. (2019) using both automatic and manual picks. We calculate cross-correlations between all catalog earthquakes within the region, using them to obtain preliminary repeating earthquake families. These are later expanded through template-matching. We then relocate these earthquakes, first using manual picking to get robust preliminary locations, and then using the double-difference algorithm HypoDD to CHALUMEAU ET AL.    is shown in brown color scale. Slow slip events (SSEs) shown in pink: 2013 SSE offshore Punta Galera (PG) , 2013 SSEs around La Plata Island (LPI) , 2015 deep SSE . The presence of a 2016 SSE near Esmeraldas (E) is suspected but not yet modeled (Hoskins et al., 2021). White stars and white lines show the epicenters and approximate rupture areas of past megathrust earthquakes (Kanamori & McNally, 1982;Mendoza & Dewey, 1984). The yellow star and yellow line show the epicenter and the 1 m contour of the rupture zone of the 2016 Pedernales earthquake . The purple lines show the 20 cm edges of the 1-month Pedernales afterslip . Plate convergence between the Nazca plate and the North Andean Sliver is from Chlieh et al. (2014). get finer relocations (Waldhauser & Ellsworth, 2000) and to calculate magnitudes. We then obtain a final classification of repeaters. A flowchart summarizing the processing steps can be found in the supplementary materials ( Figure S1).
Between April 16, 2016 and April 30, 2017, 4762 earthquakes were cataloged in the study region. The first month of data was only recorded by the permanent seismic network, while the temporarily deployed land stations were active from mid-May 2016 through mid-May 2017, and the OBSs were active from the end of May to November 2016. This mismatch in coverage between permanent and temporary stations means our seismic data does not uniformly cover the month with the most earthquakes ( Figure S2). In an effort to get a homogeneous, unbiased seismic coverage, we use a subset of 7 permanent and 10 temporary stations for repeater classification, ensuring at least 3 stations are simultaneously active at all times ( Figure S2).
We perform time domain cross correlation on the vertical seismic component in a window starting 2 s before the theoretical P arrival and encompassing both the P and S waves, for all events and for all stations within a 140 km radius of the target earthquakes. The window length is fixed at 30 s for each station to ensure that all earthquakes' S waves are included in the correlation. Only one Global Seismographic Network station, OTAV, is used for earthquakes at all distances, with a window length of 57 s. While OTAV is farther away from the correlated earthquakes, it recorded high quality, low noise data during the entire period of interest. A different filter is used for each station based on the frequency band with the highest S/N ratio, in an effort to homogenize the various stations' detection capabilities. The parameters used to compute correlation coefficients are shown in Table S1.
We initially set a lower correlation threshold of 0.9 to find preliminary repeating families. This condition is required to be met at 2 stations if 5 or more stations are recording, or at 1 station if 4 or fewer stations are recording, to classify events as preliminary repeaters. This initial pass sorted 888 of the 4762 catalog earthquakes into 364 families. Waveforms from each family are then stacked to form a single representative template. Template-matching is then performed from the April 16, 2015 to the June 30, 2017 using the fastmatched filter (FMF) code (Beaucé et al., 2017), using a correlation threshold of 0.9 and adding 432 new repeaters to our families. The use of stacks and the low correlation coefficient ensure that the catalog of repeaters is as complete as possible. We also search for new earthquakes in the year before the mainshock by lowering the detection threshold to ten times the daily average correlation between templates and the continuous data. New earthquakes are correlated together to find interseismic families, resulting in 3 new doublets.
All preliminary families are then relocated. Manual picking is first used to improve family locations with the NonLinLoc algorithm (Lomax et al., 2000), using the same 1-D velocity model as Agurto-Detzel et al. (2019). At least one repeater in every family is relocated, and its coordinates are used to set family locations. If several repeaters are relocated within a family, the average location is used. Cross correlations of P and S waves are later used to perform a higher precision relative relocation of the whole catalog using double-difference travel times and the HypoDD software (Waldhauser & Ellsworth, 2000). The region is separated into 4 subregions for relocation, and 63 stations are used, with 12 stations used on average to relocate each earthquake pair. 2483 events are relocated out of 5617, including 1119 preliminary repeaters. We show examples of relocated repeaters within families in Figure S3. With these new locations, local magnitudes (M L ) are calculated for all earthquakes, from which seismic moments are derived (see details in supporting text S1).
Lastly, we calculate correlation coefficients again on the preliminary repeaters to enforce a stricter classification for our definitive families. The threshold correlation coefficient used to classify events as repeaters is 0.95, high enough to confirm the rupture of one single source and low enough to avoid missing events during periods with poor coverage (Uchida, 2019 and references therein). Earthquakes are sorted into families if they correlate above 0.95 at one third of available stations. We only retain families with 4 or more repeaters which span more than 15 days. These additional criteria are applied because very short-lived families (with under 15 days between the first and last repeater), and a significant portion of families with 2 or 3 events, appear to be near simultaneous events, which may be the result of nearby asperities rupturing separately, rather than one single asperity rupturing repeatedly (Lengliné & Marsan, 2009). We are confident in larger families because they are based upon more cross-correlation similarity measurements. Thus, excluding doublets, triplets, and very short-lived families, we find a total of 376 repeaters grouped into 62 families. In the interseismic period, we only identified 8 repeaters each belonging to different repeating families. Examples of family waveforms are shown in Figure 2 for the PDNS station.
Although other studies sometimes use a spatial parameter to classify repeaters, like source overlap or hypocentral distance (Uchida, 2019 and references within), we have elected to base our classification solely on correlation coefficients. This is in part because 61 out of 376 repeating earthquakes could not be relocated with HypoDD, due to the small number of stations available in the first month. It is also because the location uncertainty was too large to accurately determine whether two earthquakes shared the same source. We estimated the relative relocation average error attributable to noise in the data by introducing a random error in our original travel time data and repeating the inversion 100 times for all 2483 relocated events. That error was found to be around 750 m on average. To that we must add the errors that stem from the unevenness of the network, since most of the stations used do not cover the entire period, and since all the stations used for cross-correlation differential times are on land. We performed 50 iterations with all 2483 relocated events during which we randomly take out one station for each earthquake pair. Doing so yielded an average error of about 580 m. Since the relative location errors were an order of magnitude higher than the repeaters' source areas, it was not possible to use hypocentral distance as a classification criterion without missing a significant number of repeaters. CHALUMEAU ET AL.

Location and Evolution of Repeaters
The proportion of repeaters among aftershocks through time is shown in Figure 3. The catalog used as the basis for this study (Agurto-Detzel et al., 2019) stopped on day 381 after the mainshock (April 30, 2017). Any later earthquake was found exclusively with template matching, making the proportion of repeaters past day 381 biased. The proportion of repeaters shows that the increased number of stations past the first month did not lead to a sharp increase in detections, which ensures that our detection capacities remain relatively constant. Although some of the earlier aftershocks are clearly missing due to the decreased detection capabilities right after the mainshock, the proportion of repeaters seems to remain unaffected ( Figure 3). Large aftershocks affect the number of repeaters and other aftershocks, but not necessarily their relative proportions, as in the first ∼100 days following the Pedernales earthquake the proportion of repeaters remains relatively constant at around 0.1. This suggests that the processes driving both aftershocks and repeaters are strongly linked in the region during that time.
As expected after a large earthquake, the recurrence times of repeating earthquakes drop, then gradually increase after the mainshock ( Figure S4). At least 50 families of repeaters show a consistent gradual increase of recurrence time. Some families show perturbations and small decreases in recurrence time after large aftershocks with magnitudes above 6, like in May and July 2016. These perturbations tend to be accompanied by small increases in magnitude. Most families follow an Omori-type law, like what has been documented in other parts of the world (Schaff et al., 1998), with recurrence times decaying faster as we get closer to the afterslip center. This is especially apparent near the northern afterslip patch. Due to the lack of interseismic recurrence times, we cannot determine how recurrence times at the end of our study compare to interseismic ones. However, stress perturbations from the mainshock seem to still be present, because recurrence CHALUMEAU ET AL.
10.1029/2021JB021746 6 of 20 times at the end of our study are still increasing. This is consistent with geodetic studies that suggest that postseismic slip is still happening in Ecuador in 2020 (Rolandone et al., 2020). Figure 4. Absolute aftershock relocations have a lateral uncertainty of 1.7 ± 1.3 km and a depth uncertainty of 4.3 ± 3.7 km on average. There are more significant errors in depth within 25 km of the trench, especially toward the south, where many events appear to be aligned on the 20 km depth discontinuity of the 1-D velocity model. The lateral uncertainty in this region is 2.4 ± 1.1 km, while the depth uncertainty is 7.3 ± 4.2 km. A combination of poorer station coverage near the trench, especially in the first month, and a 1-D model ill-fitted for the complex velocity structure in the area likely accounts for these large errors. Earthquakes closer to the dense land network have better constrained locations, with a lateral error of about 1.5 ± 1.1 km and a depth error of about 3.3 ± 2.9 km. These events typically cluster close to the plate interface.

Relocated aftershocks and repeaters are shown in map and cross section view in
In map view, repeater distribution appears to be similar to that of aftershocks. Most families are located in larger clusters, with only a few families being isolated elsewhere. This is consistent with suggestions that the Pedernales aftershocks are driven in large part by afterslip (Agurto-Detzel et al., 2019), similarly to repeaters. In fact, repeaters primarily occur updip of the mainshock rupture, where most of the afterslip occurs . With a few exceptions north-east of the coseismic rupture, both repeaters and aftershocks are contained within three trench-perpendicular seismicity streaks that stretch between the coseismic rupture zone and the trench. These are permanent features in the area, visible in both the interseismic and postseismic periods (Agurto-Detzel et al., 2019;Font et al., 2013). This area between the rupture zone and the trench is subject to only small coseismic Coulomb stress changes from the mainshock ( Figure S5a). Repeaters in particular mostly experience small stress increases. About 10% of repeaters and 20% of non-repeaters in the first month experience more than 10 bars of coseismic stress increase from the mainshock ( Figure S5b). Overall, the median stress increase experienced by repeaters is 3.5 bars, enough for coseismic stress changes to induce short-term triggering for some repeaters (K. H. Chen et al., 2013). However, the pattern of repeaters' and earthquakes' distribution is very different from the distribution of coseismic stress changes ( Figure S5a   . Blue circles show the two largest M6.9 and M6.7 aftershocks that occurred on May 18, 2016. The 200 mm limit of the afterslip is shown in purple . (b) Seismicity in 30 km wide cross sections. The black line is the plate interface (Hayes et al., 2012). The blue highlight corresponds to the portion of the interface that experienced more than 1 m of coseismic slip during the 2016 Pedernales earthquake . The coastline is depicted using inverted light blue triangles.
Three main regions contain repeaters: the northern afterslip patch (A-A' and B-B', Figure 4), the southern afterslip patch (C-C'), and the area north-east of the main coseismic rupture zone (near blue circles, Figure 4). Repeaters are present around the northern patch of afterslip throughout the whole period ( Figure 5). This includes observations near the trench. This area experienced SSEs during the interseismic period, though their slip amplitude was about ten times lower than that of the 2016 afterslip . While the 2013 SSEs did cause possible repeaters , neither they nor the repeaters we detected prior to the mainshock occurred closer than 40 km from the trench ( Figure S6).
Similarly, most repeaters in the south occur around the slip maximum of the southern patch of afterslip. Unlike in the north, this patch sits on a highly coupled region (Figure 1). While the repeaters in the cluster closest to the center of the patch are more numerous, and their magnitudes are generally higher, overall repeaters behave similarly in the northern and southern afterslip patches. One difference is that families near CHALUMEAU ET AL.  the trench in the south mostly stop after ∼200 days, while trenchward families in the north remain active until the end of the study period. In the cluster closest to the southern slip maximum, however, families are numerous and repeaters occur for at least ∼300 days ( Figure 5). This cluster was also active before the mainshock, as two families were activated there during the year preceding the mainshock.
Finally, the region north-east of the rupture contains a few repeaters, along with many aftershocks. It is likely that these repeaters, along with the other earthquakes, are caused directly by the M6.9 and M6.7 aftershocks on May 18, 2016, as they are only activated after the first month. Repeating families in this region, in addition to being sparse, stop after less than 200 days. This may reflect the short duration of the aseismic slip that occurred in the area, although we cannot know as no afterslip model exists past the first month.
Several regions have a high density of aftershocks but no repeaters. There is a high density of aftershocks between the two coseismic slip maxima, as well as at the southern end of the rupture zone, but no repeaters occur there. Within the coseismic rupture zone, it is likely that aseismic slip is prevented by the nearly total stress release during the mainshock. Even aftershocks within the rupture area may not occur at the plate interface, although this cannot be confirmed given the depth uncertainties and the lack of focal mechanisms (Agurto-Detzel et al., 2019;Cordero et al., 2020). The other clusters without repeaters are located north of the rupture, at −79.90°W (Figure 4a). Aftershocks north of 0.8°N mostly occur in December of 2016, starting on December 2nd, in the overlying plate  and are not likely to be directly linked to afterslip. These events do not seem to directly relate to the aftermath of a large earthquake, although three earthquakes above magnitude 5 occur in that region between December 12th and December 20th, during the peak of seismic activity. They may be related to a swarm occurring in the Esmeraldas sequence 30 km to the east, itself happening on a crustal fault and probably in concurrence with aseismic slip, which starts on November 28th (Hoskins et al., 2021). Meanwhile, aftershocks in the cluster between 0.6°N and 0.8°N are primarily associated with the large M6.9 and M6.7 aftershocks of May 18, 2016 and the M6 and M6.2 earthquakes on July 11, 2016.

Repeaters and Afterslip
Repeaters identified in this study seem strongly associated with the two main afterslip patches. We aim at quantifying that relationship by using the GPS-derived 30-days afterslip model developed by . Since no geodetic afterslip model currently exists past the first month, we primarily focus on families that are active within that time period. To calculate slip from repeaters, we use the equation developed by Nadeau and Johnson (1998) Here the seismic moment M 0 is in dyn.cm and the slip d is in cm. Although the equation is an empirical relation based on geodetic data from Parkfield, California, it has been used successfully in several subduction settings, including Japan (Igarashi, 2020;Matsubara et al., 2005;Nomura et al., 2017;Uchida et al., 2016;Uchida & Matsuzawa, 2013;Uchida, Nakajima, et al., 2009;Uchida, Yui, et al., 2009) and Chile (Huang et al., 2017;Kato et al., 2016;Meng et al., 2015). The model has been criticized in different studies as being physically implausible (Beeler et al., 2001;Sammis & Rice, 2001), as it assumes that no aseismic slip occurs on the asperities and that stress drop is a function of magnitude, therefore predicting unrealistic stress drops for small earthquakes (Beeler et al., 2001;Sammis & Rice, 2001). A competing model developed by Beeler et al. (2001) would not materially alter our estimates, since it tends to give similar results to the NJ equation Mavrommatis et al., 2015). We therefore feel confident in using the NJ equation here.
The resulting slip estimates are shown in space and time in Figures 6a and 7. From the geodetic model of afterslip , we estimate the amount of slip that occurred between the first and last repeater of the month within a family. Some agreement is found in the spatial distribution of slip from repeaters and GPS (Figure 6a), as well as with its temporal evolution (Figure 7). is also a good temporal agreement between smaller-scale surges in GPS-derived slip and in repeater-derived slip, around days 5, 10 and 15, although these perturbations appear slightly later in the GPS data. However, upon closer inspection of individual families, the relationship between GPS and repeater-derived slip breaks down. First, there is a lot of spatial heterogeneity in the slip calculated ( Figure 6a). This is true even when considering long-term slip ( Figure S7). More explicitly, GPS and repeater-derived slip of individual families do not appear to have a strong relationship (Figure 8a), even when considering families only from one region. We confirm this by plotting the average moment of individual families against their average experienced slip as modeled from GPS data (Figure 8b). We find a large scatter and a poor correlation, despite a similar overall trend to the repeaters at Parkfield (Figure 8b).
The complexity of the slip/repeaters relationship is reinforced by the absence of repeaters in areas with maximum afterslip as well as in areas with no slip (Figure 6). The distribution of non-repeaters follows the slip distribution, at least for slip values larger than 100 mm, which are better resolved. Repeaters are more concentrated in areas of moderately large slip, between 100 and 450 mm. They are however absent from areas which experienced more than 500 mm of slip. In map view, it is apparent that repeaters mostly surround the two afterslip maxima (Figure 6a). This is true in the southern as well as in the northern slip regions, although the different coupling values may imply different asperity densities and productivity of seismicity. This reinforces the idea that something beyond simply the amount of slip is influencing the repeaters and their distribution. Since repeaters appear to be mostly on the edges of the main patches of afterslip, afterslip gradient is a likely candidate mechanism.

Repeaters and Afterslip Gradient
Repeaters seem to concentrate in areas with large afterslip gradients, and the gradient associated with repeaters is higher than for the aftershock population ( Figure 9). More generally, it seems non-repeaters are present at all gradient values, while repeaters occur only where gradient is moderate to high. To better CHALUMEAU ET AL.

10.1029/2021JB021746
10 of 20  . Outlined patches represent slip calculated from repeating families with 2 or more events in the first month, using the NJ equation. For easier comparison both among families and with the geodetic model, NJ slip from each family was extrapolated to one month, using the ratio of slip calculated using the NJ equation to slip modeled by GPS between the first and last repeater of the month. Each cell averages slip from families within a 2 km radius. (b) Percentage of first month repeaters and non-repeaters located in areas experiencing a given amount of afterslip (sampled over 5 km 2 ), as estimated from the GPS model . The total slip is the amount of afterslip experienced between the mainshock and the time of a given earthquake at that location, according to the GPS model. Non-repeaters are shown in red and repeaters are shown in blue. The gray bars show the percentage of the total study area experiencing a given amount of afterslip.  illustrate this, and since the afterslip patches are roughly concentric circles, we plot the distribution of earthquakes and repeaters against distance from the center of the slip patch at the end of the month (Figure 10). We see that the distribution of repeaters has a stronger correlation than non-repeaters with the afterslip gradient, in the south especially. In both the north and south, there are almost no earthquakes at the center of the patch where the slip is highest. Were the amount of slip to most strongly influence the location of repeaters or non-repeaters, their numbers would be highest at the center and decay with distance, but that is not the case. Instead, in the north and south the number of repeaters peaks at 25 and 15 km away from the center of the patch respectively, which is nearer to where the afterslip gradient peaks. Since the two main patches of afterslip remain relatively stationary throughout the month, their edges remain fixed in space ( Figure S8). As a result, afterslip gradient, which is largest at these edges (Figure 9a), grows consistently throughout the month, while its spatial extent remains near constant. This implies that the location of repeaters with regards to slip is constant over the whole month. Overall, this suggests that aseismic slip gradient may be a major control on the distribution of repeaters.

Discussion
We have found repeating families in the aftershock sequence of the Pedernales earthquake, which qualitatively appear to be associated with afterslip, both spatially and temporally. These repeaters are mainly located primarily at the edges of the two main afterslip patches, where afterslip gradient is highest, rather than at their center, where slip is highest. Slip gradient seems to be an important factor controlling repeaters. We will now explore this relationship in detail.

Influence of Uncertainties on the Results
Before interpreting our results, it is important to assess their robustness by evaluating the impact of uncertainties. While gaps in station coverage ( Figure S2b) have been shown to be in large part mitigated, errors in repeater detection and classification are a potential issue. We mitigate errors by using a high correlation CHALUMEAU ET AL.

10.1029/2021JB021746
12 of 20 threshold. This likely means that some repeaters are missing from the data set, but it ensures we have high confidence in the families that we do interpret. Our data set is additionally limited by the repeaters' magnitudes. While the total catalog local magnitudes are mostly between 1.5 and 3.5, repeaters' local magnitudes are concentrated between 2.2 and 3.4 ( Figure S9a), suggesting that repeating families with magnitudes lower than 2.2 are likely incomplete or undetected. This could affect the spatial distribution of repeaters in areas where the average magnitude is low, like northeast of the mainshock rupture and within the mainshock area. The latter 2 regions, however, have a magnitude of completeness of 2 and 1.5 respectively ( Figure S9b), so it is less likely that repeating families were missed there. A lack of completeness within a family could impact its slip estimate, and may explain low values when compared to GPS-derived slip. On the other hand, CHALUMEAU ET AL.

10.1029/2021JB021746
13 of 20 Figure 10. Average afterslip, afterslip gradient and seismicity distribution with regards to distance from the center of the northern (left) and southern (right) patch of afterslip at the end of the first month. The top panels (a) show the evolution of geodetic afterslip with distance. The panels below (b) show the average afterslip gradient as a function of distance. Half of all radial profiles fall within the shaded area. The third panels (c) show the distribution of all repeaters (blue) and non-repeaters (black) with regards to distance. The bottom panels (d) show only repeaters from families with median magnitudes above 3 and minimum magnitudes above 2.7 (blue), and non-repeaters with magnitudes above 3 (black). some non-repeaters may have been misclassified as repeaters due to an inadequate upper threshold of the filter used for correlation. This may be an issue for small-magnitude earthquakes with corner frequencies significantly above the filter's upper frequency boundary (Uchida, 2019 and references within). However, even when examining families with median magnitudes larger than 3 and minimum magnitudes larger than 2.7 ( Figure S10), there still is no visible relationship between slip estimations from GPS and repeaters. The placement of repeaters with regards to slip and gradient also remains similar ( Figure 10) Since we compare the distribution of repeaters to that of non-repeaters, the misclassification of repeaters as non-repeaters might also have an impact on our analysis. To avoid this issue, we create a buffer by excluding from non-repeaters every earthquake that has a correlation coefficient above 0.9 with another earthquake. The threshold of 0.9 is chosen as it is used by some studies to classify repeaters (e.g., Uchida, 2019 and references within), and thus the nature of these earthquakes as repeaters or non-repeaters is uncertain.
The different sensitivities of GPS and repeaters necessarily impact our comparisons with GPS models. The GPS slip model is smooth by necessity, while the repeaters' slip estimates are not. A highly heterogeneous slip at the interface would be smoothed when recorded at the surface and therefore be poorly resolvable by GPS models, but would explain the high degree of heterogeneity of the repeaters' slip estimates. However, if families did concentrate in small heterogeneous areas of higher-than-average slip unresolvable by the GPS, we would probably expect repeaters to overestimate the overall slip, which is seemingly not the case here.
GPS models themselves have uncertainties to take into account. We consider the center of the southern patch of afterslip to be well constrained, as a different afterslip model shows it to remain within 5 km (Figure S11) . In the north, the center of the afterslip deviates by 10 km, but the landward edge of the patch remains in place, leaving the trenchward side as the most uncertain ( Figure S11).
Another source of error is the uncertainty of earthquake locations. Although lateral earthquake location errors are small on average, they have to be combined with errors in the geodetic model of afterslip, especially near the trench, in order to ensure our observations regarding the relationship of repeaters to slip and to slip gradient is correct. For example, the underestimation of slip from an earthquake cluster close to the center of an afterslip patch compared to one given by the GPS model (Figure 6a), could be explained by the cluster being wrongly located too close to the afterslip patch center. In fact, most of the repeaters within the high slip areas tend to have low estimations of slip compared to the GPS, while the families further from the patch tend to have slip estimations that are higher than the GPS (Figure 6a). If we were to assume that our repeaters' slip estimations are correct, and this discrepancy is due entirely to location errors, then if anything, this confirms that repeaters should be more concentrated on the edges and less in regions of high slip. Alternatively, the regularization parameters used to derive the GPS model could have led to an overconcentration of slip at the center of the two slip patches . A more spread out slip may explain some of the regional discrepancies between the GPS and family slip, although it would still not explain why repeaters remain located on the edge of the afterslip. To account for location and geodetic model errors and smoothing, we give each family an error for GPS slip that corresponds to the standard deviation of the GPS slip within a 15 km radius of the family's location (Figure 8a). Even when taking the average of both NJ slip and GPS over a 30 km region, the agreement between the two is still poor ( Figure S12), suggesting that location errors are neither the reason for the lack of repeaters near the center of the afterslip, nor the cause of the discrepancy between repeaters and GPS slip.

What Controls Repeaters Occurrence?
We have shown a significant correlation between the location of repeaters and the gradient of the afterslip. Additionally, we have found that the relationship of geodetic slip estimates to repeaters is not linear in the region, even taking uncertainties into account. These findings hint at afterslip gradient being an important factor controlling repeaters. We will now discuss how feasible that statement is, in light of studies in other regions.
A few factors are widely recognized as influencing the locations of repeaters. The most important factor is the presence of aseismic slip (Waldhauser & Schaff, 2008). This means repeating earthquakes are typically outside, or around, historic large earthquake rupture zones and highly coupled regions (Chaussard et al., 2015;Igarashi et al., 2003;Ryder & Bürgmann, 2008;Templeton et al., 2009;Uchida et al., 2003;Uchida, Yui, et al., 2009), although on occasion they have been found inside past (Uchida & Matsuzawa, 2013) or future rupture zones (Uchida & Matsuzawa, 2011). In Ecuador, repeaters are indeed found in slipping areas outside of the rupture zone, although interseismic coupling appears to have only a limited bearing on the locations of repeaters, which seemingly occur at all values of coupling ( Figure S13). Another common observation that fits with our data is the location of repeaters within areas containing dense microseismicity (Kato & Igarashi, 2012;Ryder & Bürgmann, 2008). Repeaters can only occur on asperities that are seismogenic, meaning that frictional features always remain the primary factor controlling their distribution. As such, repeaters cannot occur in completely aseismic areas, any more than they can occur in locked areas.
To our knowledge, there are no studies examining specifically the locations of repeaters with regards to afterslip gradient. However, repeaters occurring primarily away from the aseismic slip center is not a unique observation. Along the San Andreas fault (California), earthquakes, including repeaters, tend to align along horizontal streaks (Lengliné & Marsan, 2009;Rubin et al., 1999). These streaks have been proposed to be boundaries between seismic and aseismic slip (Sammis & Rice, 2001), although they are more generally thought to have a geological origin (Rubin et al., 1999). Joint geodetic and seismological studies show that repeaters appear to be absent from areas that slip the most, preferring instead areas of intermediate slip (Chaussard et al., 2015;Templeton et al., 2009). As a result, studies using only repeaters can underestimate the total slip, since they do not sample the slip maximum (Templeton et al., 2009). While the NJ relationship was derived first in Parkfield, its application is not straightforward in the whole of California. There is always some small-scale heterogeneity, so that in some cases the M 0 -Tr relationship only weakly follows the NJ trend (Lengliné & Marsan, 2009). When comparing InSAR models to repeaters slip, slip profiles tend to agree, but are different in detail (Ryder & Bürgmann, 2008), with sometimes large scattering in the calculated versus InSAR rates (Chaussard et al., 2015). In some cases, particularly for short-term, faster creep, the parameters of the NJ equation need to be changed to better reflect the data (Khoshmanesh et al., 2015), suggesting that the relationship of M 0 to slip is not the same everywhere at every time.
The most studied subduction zone with regards to repeating earthquakes is northeastern Japan. During interseismic times, repeaters appear to be located in slow and fast creeping areas within the seismogenic zone , with a good agreement between GPS measurements and repeaters' slip (Nomura et al., 2017), and relatively small scatter . After the 2011 Tohoku earthquake, however, the density of repeaters appears to have been low within patches of significant afterslip, according to several afterslip models Silverii et al., 2014;Uchida & Matsuzawa, 2013).
It is perhaps best to compare our study to repeating earthquakes occurring after a large earthquake, to determine whether they do tend to organize around the afterslip maxima or not. On the one hand, the aftermath of the 2015 M W 8.4 Illapel earthquake clearly shows repeaters occurring largely within the area of maximum afterslip, and the slip estimates from repeaters agree well with the GPS models (Huang et al., 2017). In this case, however, the afterslip was low in amplitude, staying below 50 cm, and large in extent (Huang et al., 2017;Shrivastava et al., 2016), meaning there likely was not a large regional slip gradient anywhere. On the other hand, the 2012 M W 7.6 Nicoya earthquake offers a clear case of repeaters being located on the afterslip edge. Costa Rican repeaters exist largely between the two main patches of afterslip (Chaves et al., 2020;Hobbs et al., 2017;, despite the fact that some non-repeating earthquakes occur closer to the two slip centers . Like in Ecuador, the two patches of afterslip in Costa Rica both have a high amplitude, here up to 1 m, and a small spatial extent, less than 50 km, meaning the slip gradient would be very high ( Figure S14). There is also a poorer agreement between the slip derived from repeaters using the NJ equation and the slip from the GPS model, further highlighting the similarities between Costa Rica and Ecuador ( Figure S14).
It therefore seems that cases similar to ours exist elsewhere, and that the behavior of repeaters with regards to slip and gradient may be dictated by the amplitude of the spatial gradient. The question remains of the extent to which slip gradient influences repeater behavior. There are two main explanations that we can propose. The first is that gradient is simply the indicator of a change in the mechanical properties of the fault, which are themselves responsible for the occurrence of repeaters. The second is that gradient, or the stress that accompanies it, is directly controlling repeaters.
The first possibility is that the high gradient results from aseismic slip is being stopped by a large locked asperity, with repeaters nucleating at its edge. This idea was proposed by Sammis and Rice (2001) as a way to explain the relationship between repeating earthquakes' recurrence times and seismic moments as proposed by Nadeau and Johnson (1998) without requiring uncommonly high stress drops for repeating earthquakes. This model implies that slip on a repeating asperity is not directly related to overall aseismic slip. It is instead being controlled largely by the presence of the large asperity. Similarly, Anooshehpoor and Brune (2001) model repeaters as asperities within creeping patches inside of larger asperities, invoking slip velocity shielding to explain the scaling of repeaters moment to the overall measured aseismic slip. Both models, however, necessitate the concentration of repeaters along the boundary between large asperities and creeping regions, which is not true for most of our repeaters ( Figure S13).
A similarly high gradient could be observed if, rather than being stopped by a single large asperity, aseismic slip was slowed down by a multitude of smaller, concentrated asperities. This may still result in slip velocity shielding between asperities of different sizes, explaining heterogeneities between families and a complex relationship of repeaters to slip. In particular, it would explain why, in some regions, nearby repeating aftershocks have very different behaviors but still, when averaged, follow Omori's law and the NJ scaling relationship (Lengliné & Marsan, 2009). However, in our case, slip averaged over a region still does not match with slip derived from geodetic models, putting into question the direct relationship of slip to repeaters ( Figure S12).
More generally, the absence of repeaters close to the afterslip maxima could be explained by that region being completely aseismic and therefore unable to produce any earthquakes. In that case, the high spatial gradient could be a passive indicator of a rheology change that allows for the presence of earthquakes and repeaters. Like with the previous possibility, this offers no explanation as to why the fit of geodetically derived and repeater-derived is so poor. More importantly, it raises the question of why repeaters would be present in areas where the afterslip is moderate but the interseismic coupling is close to 0% ( Figure S13), yet absent from areas where the interseismic coupling and the afterslip are both higher. The transition of repeater-prone to completely aseismic regions would need to be more closely investigated for this to be answered conclusively.
It is also possible that stress transferred from the afterslip, which is directly related to the gradient, is a driving force behind repeaters in the region. While stress is released within the afterslip patches, it induces local stress concentrations at its edges (Andrews, 1976;Rice, 1993;Scholz, 2019;Wynants-Morel et al., 2020). Since the edges of the afterslip patch, along with the areas experiencing the largest slip gradient, remain static throughout the whole period, a recurrent loading of asperities through continuous stress increase at the edge is possible. We know that stress increase from SSEs can trigger seismic swarms at their edges, as has been shown in Hawaii (Segall et al., 2006) and in New Zealand (Bartlow et al., 2014). No study has currently looked at the applicability of that model for repeaters, but it may be worth investigating in detail. This link between gradient and repeaters has important implications for the estimation of slip using repeaters.
Although the presence of repeaters itself remains an indication of aseismic slip occurring in its vicinity, the fact that repeaters occur preferentially at the edges of the slip, and not at its peak, necessarily makes it difficult to accurately quantify the total slip. Additionally, if stress from the afterslip does indeed drive repeaters, then getting an accurate estimation of slip from the latter would likely be challenging.

Conclusion
We have conducted a systematic search for repeating earthquakes in 14 months of postseismic data of the M W 7.8 2016 Pedernales earthquake in Ecuador, as well as 12 months of interseismic data. Repeating earthquake families were found in the vicinity of the two main patches of afterslip in the 14 months following the mainshock. We show that there appears to be a spatial relationship between the location of repeaters and the spatial gradient of the afterslip. Repeating earthquakes seem to concentrate primarily on the edges of the afterslip, where slip gradient is high, rather than at its center, where slip is high. While structural controls undoubtedly remain the most important factor leading to the presence of repeaters, our results suggest that stress accumulation on the edges of the afterslip may be a driving force behind repeater activity. working together with the French Institut de Recherche pour le Développement (IRD) through the Joint International Laboratory "Earthquakes and Volcanoes in the Northern Andes." (LMI SVAN), for providing instrumentation, logistics and field support. Thank you as well to the PASSCAL facility of the Incorporated Research Institutions for Seismology (IRIS) for supporting instrumentation during the aftershock deployment. The authors also want to thank the Institut National des Sciences de l'Univers (INSU of CNRS) and the University of Liverpool UK for their support. Finally, the authors thank everyone who manually picked events for the elaboration of the one-year catalog, especially Cristina Gonzalez, Monica Segovia , Guillermo Viracocha, Javier Santo, Sandro Vaca, Andrea Cordova, Yahyeh Osman, Cedric Bulois, Andrea Peuzin, Rim Husseini, and Alexandre Boughanemi.