Probing the Northern Chile Megathrust With Seismicity: The 2014 M8.1 Iquique Earthquake Sequence

We used data from >100 permanent and temporary seismic stations to investigate seismicity patterns related to the 1 April 2014 M8.1 Iquique earthquake in northern Chile. Applying a multistage automatic event location procedure to the seismic data, we detected and located ~19,000 foreshocks, aftershocks, and background seismicity for 1 month preceding and 9 months following the mainshock. Foreshocks skirt around the updip limit of the mainshock asperity; aftershocks occur mainly in two belts updip and downdip of it. The updip seismicity primarily locates in a zone of transitional friction on the megathrust and can be explained by preseismic stress loading due to slow‐slip processes and afterslip driven by increased Coulomb failure stress due to the mainshock and its largest aftershock. Afterslip further south also triggered aftershocks and repeating earthquakes in several EW striking streaks. We interpret the streaks as markers of surrounding creep that could indicate a change in fault mechanics and may have structural origin, caused by fluid‐induced failure along presumed megathrust corrugations. Megathrust aftershocks terminate updip below the seaward frontal prism in the outer continental wedge that probably behaves aseismically under velocity‐strengthening conditions. The inner wedge locates further landward overlying the megathrust's seismogenic zone. Further downdip, aftershocks anticorrelate with the two major afterslip patches resolved geodetically and partially correlate with increased Coulomb failure stress, overall indicating heterogeneous frictional behavior. A region of sparse seismicity at ~40‐ to 50‐km depth is followed by the deepest plate interface aftershocks at ~55‐ to 65‐km depth, which occur in two clusters of significantly different dip.


Introduction
The frictional behavior of subduction megathrusts is heterogeneous, exhibiting regions slipping either in a stable, that is, aseismic, manner, or unstably during earthquakes (e.g., Lay et al., 2012;Perfettini et al., 2010). Mapping and understanding this variability is crucial because it ultimately determines the earthquake hazard potential. Measuring and inverting interseismic surface strain patterns due to variable plate locking based on GNSS measurements is one approach to achieve this (e.g., Bürgmann et al., 2005;Métois et al., 2013;Savage, 1983). However, sparse and one-sided data acquisition limits the resolution of such models. While the interseismic period is often seismically quiet, the coseismic and postseismic periods open a window in which the full spectrum of slip modes can be observed. Whereas coseismic slip inversions map the main asperities, aftershocks illuminate mainly the marginal areas (e.g., Asano et al., 2011;Beroza & Zoback, 1993;Das & Henry, 2003;Hsu et al., 2006;Kato et al., 2010;Mendoza & Hartzell, 1988;Perfettini et al., 2010;Wetzler et al., 2018;Woessner et al., 2006) for which slip behavior is varied and often not known.
Postseismic slip can be estimated from geodetic observations (e.g., Hoffmann et al., 2018, for the Iquique earthquake), albeit with the same limitations due to data coverage as for interseismic observations. Other indicators of aseismic slip may be repeating earthquakes (e.g., Igarashi et al., 2003), which are interpreted

Spatial Distribution of Events
The events define two main depth ranges, the first of which covers an offshore region of shallower seismicity (hypocentral depths of 0-65 km) concentrated between 19.5°S and 21°S and forms an eastward dipping triple seismic zone (Bloch et al., 2014;Sippl et al., 2018), which we subdivided into upper plate (here defined as events located more than 10 km above the slab surface and east of 70.7°W), intraslab (events more than 5 km below the slab surface), and presumed interface (remaining events) seismicity. The latter group of earthquakes belong mostly to the foreshock and aftershock sequences of the Iquique earthquake and will be analyzed in detail in this article (see Figure 2a). Since west of 70.7°W depth uncertainties are higher (>6 km; see Figure S2), we do not distinguish events located in this region and for simplicity color them as interface seismicity but will treat them with caution in our interpretation. The deeper onshore (east of~69.5°W) band of intermediate-depth seismicity (80-to 120-km depth) corresponds to intraslab background seismicity. These events have been previously examined in great detail Sippl et al., 2019), so we omit them from our analysis here. In the following, we will   are plotted in orange, yellow, and brown color, respectively. Gray shaded bathymetry is from . The limit between the lower and middle continental slopes (MLS) is taken from Maksymowicz et al. (2018). Magenta contours show the slab surface (Hayes et al., 2018) at 25-, 45-, and 65-km depth. Relocated earthquakes are presented in map view and as projections onto longitudinal and latitudinal planes, colored by days since the mainshock (events before the mainshock on top of events after). Lower right subplot shows a histogram of the event depths. Labeled regions outline an area of clearly decreased seismicity (SP) in the latitudinal section (dashed line in map view) and areas of low seismic activity (LN and LS) in the map view, which are discussed in the text. Red box outlines the region shown in Figure 2a. (a) Zoom into the study region, showing different seismicity clusters (labels C1-C4), streaks (labels S) and regions of high (labels R1-R3) and low (labels LS and SP) seismic activity discussed in the text. Cross Sections AA′ to FF′ of subfigure b are marked. Profile GG′ is shown in Figure 8. Slip contours of the mainshock (orange, 2, 6, and 10 m) and the M7.6 aftershock (brown-dashed, 0.56 and 0.85 m) are from Duputel et al. (2015). Epicenter of these events and the M6.6 foreshock are shown as in Figure 1. Bathymetry features are marked with green (embayments) and magenta (NW-SE antiforms and synforms; . Aftershocks are plotted as blue circles, on top of foreshocks in red. Presumed noninterface seismicity is plotted in green for upper-plate (circles) and intraslab (crosses) events. (b) Seismicity Cross Sections AA′ to FF′. Slab surface (Hayes et al., 2018) and swath bathymetry (upper curve; topography at depth > 0 km) are drawn in black. In Profile BB′, we also plot a histogram of event numbers, the coseismic slip (orange curve), and the location of the continental Moho (label M, Sodoudi et al., 2011). In Profile DD′, focal mechanisms are plotted as open colored circles (see Figure 3 for a definition of the event classes). provide a more detailed description of observed seismicity features (clusters, streaks, and regions of high and low seismic activity), moving from west to east, that is, in the downdip direction, through three different subregions. We will adopt the coseismic slip model from Duputel et al. (2015) for comparison with the seismicity patterns that we observe, as it was derived from the most complete set of observations to date (seismic strong motion and high-rate GNSS time series, InSAR-and GNSS-derived static surface displacements, tide gauge, and DART buoy measurements of tsunami height) using a Bayesian inversion scheme, which does not apply smoothing. We will further discuss some features of this model relevant to our interpretations in section 4.1. 3.1.1. Updip Seismicity (0-25-km Depth) Earthquake mechanisms in the region of the shallow plate interface just landward of the trench show a predominance of thrust faulting striking in a similar plane as the mainshock mechanism, but featuring a variety of dip angles. There are also smaller populations of differently oriented thrust events as well as some strikeslip and normal-faulting earthquakes (Figure 3).
In the northern part of the shallowest region of the megathrust, foreshocks and aftershocks occur in dense clusters that are located at the updip edge of the moderate-to-high (6-10 m) coseismic slip zone, with aftershocks slightly but systematically displaced updip with respect to the foreshocks (region labeled R1 in map view and in Cross Section BB′ in Figure 2).
Figures 4a and S5a show the location of earthquakes at the latitudes of the two main streaks (S1 and S2) in map view and cross sections. The continuous, eastward dipping geometry outlined by these events and most of the focal mechanisms in this area indicate thrust earthquakes presumably located along the plate interface.
Most of the seismicity updip of the mainshock main rupture appears to locate on the plate interface in the northern cross sections (Figure 2b). However, some events appear to branch away upward from the slab surface and intersect the seafloor landward of the trench in at least one of the southern sections, indicating possible activation of a splay fault (Profile DD′; Figure 2b).

Mainshock Main Rupture and Central Seismicity (25-to 45-km Depth)
Deeper on the megathrust, the highest coseismic slip area of the Iquique earthquake contains only a few foreshocks and aftershocks (Figures 2a and 2b, Cross Section BB′).
In the along-strike direction, this region is continuous into the area of highest slip of the M7.6 aftershock and its hypocenter to the south, which is surrounded by several minor~EW striking seismicity streaks, similar to the ones observed at shallower depths.
Downdip of the mainshock main rupture, an~SW-NE oriented region of increased aftershock seismicity that partially overlaps moderate-to-low (2-6 m) coseismic slip can be observed at depths of 25-45 km in the Cross Sections BB′ and CC′ and is labeled R2 in map view ( Figure 2). Interestingly, no foreshock activity is observed at such depths. A similar region of aftershocks (R3), slightly more oriented toward E, extends from north of the area of maximum coseismic slip to~19.4°S. Aftershocks in R3 connect to the concentration of foreshocks located updip of the mainshock's main slip and extend downdip to a similar depth as events in R2 (Cross Section AA′). In fact, between about 35 and 45 km, R2 and R3 merge, forming a belt of aftershocks just around the downdip edge of the deepest significant coseismic slip. Considerable upper plate seismicity is also noticeable above these two regions (Cross Sections AA′, BB′, and CC′).
To the south and north of Regions R2 and R3, respectively, two largely aseismic areas (labeled LS and LN in Figures 1 and 2) are observed that extend at least 50 km in both directions. Northward, low seismicity rates (LN) prevail at this depth until the limit of the region covered by our catalog, whereas a number of isolated clusters and small~EW seismicity streaks of aftershocks are observed to define the southern termination of the southern low seismicity area (LS) at~20.5°S (Figures 1 and 2).
At this depth range, the belt of intense aftershock activity formed by R2 and R3 downdip of the main rupture outlines a sharply defined inclined plane in cross section (the plate interface), which is shallower and shows significant topography relative to the slab surface model Slab2 (Hayes et al., 2018; profiles AA′ and BB′ in Figure 2b). Focal mechanisms in this region are mostly consistent with rupture along the east dipping plate . We classified the mechanisms by their rake angles into thrust with strike ±45°from north (magenta and green), thrust with significant rotation (orange; deviation >45°from north), normal (cyan), and strike-slip or oblique faulting (gray). We further subdivided the first group into events with dip differences ≤15°("interface," magenta) and >15°("steeper," green) between the mechanisms and the plate interface plunge (Slab2; Hayes et al., 2018). Large map covers the time period after the Iquique earthquake, smaller panel before the mainshock. Other features as in Figure 2a. The same colors and definitions are used in Figure 2b, profile DD′. . Detailed analysis of one seismicity streak (labeled S1 in Figure 2a). (a) Map view and depth profile of the events in the region of the streak (red: foreshocks; blue: aftershocks). Focal mechanism types are shown in cross section as bold colored circles (as in Figure 3). (b-d) Two-dimensional histograms for events in the streak (within cyan box shown in map view of subfigure a). Black curves show the results for the median of the data, which was calculated for a moving window with width of 1 km in subfigures b and d and 2°in subfigure c. The number of data used in each plot is shown in parentheses. (b) dt (S-P) versus interevent distance. Colored curves show the results for the horizontally (horiz), along-dip (slab), and vertically (vert) oriented artificial streaks. (c) P and S wave arrival time differences (dt) for all event pairs within S1, normalized by the interevent distance d divided by velocity V, versus azimuth. Red curve shows a theoretical cosine (cos) function expected for a streak oriented in EW direction. (d) CC coefficient versus interevent distance. Colored curves represent the noise level calculated from noise-noise (no-no) and noise-signal (no-sig) CC. Similar plots for the second major seismicity streak S2 are shown in Figure S5. interface ( Figure 3). There is, however, a cluster of presumed normal faulting events observed at~35-km depth in swath C, at the western edge of the low seismicity area LS (Figures 2a, 3, and 4a).

Deep Interface Seismicity (45-to 65-km Depth)
Further downdip, a depth interval with only sparse aftershock seismicity is observed at~40-to 50-km depth; it is labeled SP in Figure 1 (map view and latitudinal profile) and is also visible as a subtle minimum in the event depth histogram. Its depth extent is not constant along strike; it appears wider to the south (up to 25 km), where it connects to the southern postseismically aseismic patch LS described in section 3.1.2. This region of low activity separates the previously described offshore lower part of the central seismicity from the deepest interplate events, which are situated east of the coastline (Figures 1 and 2a). Here seismic activity concentrates in three clusters of aftershocks, marked C2, C3, and C4 in Figure 2. They each contain a large number of events (168, 320, and 539 events, respectively) produced by low-magnitude earthquakes (0.9 < M < 3.3, 0.7 < M < 3.9, and 0.5 < M < 4.5, respectively), which are well detected under the footprint of the network. The two deeper clusters form two fingers beneath the Coastal Cordillera to depths of 55-65 km. Even though we cannot be sure that events within C3 and C4 are located along the plate interface, they appear to connect further updip to the already described inclined planes of central-depth seismicity, but their on-profile plunges are observed to be slightly (C4) or clearly (C3) steeper than suggested by the Slab2 model ( Figure 2b, Cross Sections BB′ and DD′), implying a strong plate curvature. Alternatively, they could also indicate faulting in the lower plate, particularly the steeper dipping cluster C3. However, this alternative interpretation would imply a small angle between the plate surface and the presumed fault plane in the subducting oceanic plate, whereas known examples of elongated or sheet-like clusters in the lower plate tend to indicate faults at high angle to the plate interface (e.g., Fuenzalida et al., 2013), which also agrees with the focal mechanisms typically observed for intraplate events (Bloch et al., 2018). In addition, the preexisting fabric of the subducting plate is generally characterized by the high angle normal faults formed at the outer rise.

Temporal Evolution of Seismicity
The rich seismicity content of our catalog allows us to inspect the time evolution of the foreshock and aftershock series with a high level of detail.
We observe that the foreshocks (from 15 March to 1 April, the mainshock origin time) started updip and NW from where the mainshock rupture later nucleated. Then the region south of the future mainshock epicenter was rapidly activated, with the occurrence of the largest M6.6 foreshock ( Figure 5a and Movie S1). Events next propagated southward and northward for a couple of days, activating in the south one of the most prominent observed seismicity streaks (S1) on 17-18 March (events occurred 15.5-14 days before the mainshock; Figures 5a and 5b). From this point onward, seismicity migrated almost continuously northward to where the Iquique earthquake finally ruptured (Figure 5a and Movie S1). These patterns of foreshock seismicity migration have been associated with slow-slip events (Ruiz et al., 2014;Socquet et al., 2017), which would have increased the Coulomb stress in the future mainshock rupture area and thus contributed to its nucleation and propagation (Kato & Nakagawa, 2014;Yagi et al., 2014).
In the first 2 days after the mainshock, seismicity gradually progressed southward until the M7.6 aftershock nucleated on 3 April, while aftershocks also occurred in the updip and central parts of the main coseismic rupture, including a few events in the two main seismicity streaks (Figure 5c and Movie S2). After the M7.6 event, the region surrounding its epicenter, the minor streaks there as well as the cluster (C1) westward of it were rapidly activated. The southern and northern onshore cluster (C4 and C3), seismically inactive until then, started to form just few hours later (Figures 5d and 5e and Movie S3). During the following days and up to the 16th day after the mainshock (17 April), while the major streaks (S1 and S2) and clusters (C1, C3, and C4) remained active (Figure 5d), several small clusters of aftershocks progressively populated the region updip of the main rupture (R1) and formed the regions R2 and R3, previously described in section 3.1.2 (Movie S3). From 18 April until the end of the study period, all the major streaks and clusters, as well as regions updip and downdip of the mainshock asperity, remained active ( Figure 5f). Movie S4 shows a more detailed view of these aftershocks until 30 April. We note that earthquakes along the major streaks S1 and S2 did not propagate in any particular direction but occurred randomly over time (Movie S5).  Figure 6 shows an overlay of the aftershocks (since Day 17 after the mainshock) on a model of Coulomb Failure Stress change (ΔCFS) and the geodetically derived afterslip model (cumulative between Days 17 and 334 after the mainshock) of Hoffmann et al. (2018). The time period was selected for the best overlap of geodetic and seismic station coverage. Afterslip is concentrated within two aseismic regions (see section 3.1.2) NNE and SSE of the interseismically highly locked and coseismically ruptured asperity (Li et al., 2015;Schurr et al., 2014;Figures 6 and S6). Modeled afterslip is dominant in the downdip part of the rupture; however, synthetic tests (Hoffmann et al., 2018) indicate that farther offshore afterslip is not well resolved. The observed anticorrelation between aftershock seismicity and regions of high afterslip occurs mostly where afterslip is well resolved (Figure S15d of Hoffmann et al., 2018). For the ΔCFS model ( Figure 6), we projected the coseismic slip from the mainshock and the M7.6 aftershock (Duputel et al., 2015) onto a representation of the megathrust (Hayes et al., 2018) that consisted of 232 individual fault patches. We used the software Coulomb3.4 (Lin & Stein, 2004;Toda et al., 2005) and chose the following parameters: Young's modulus 8 × 10 4 MPa, Poisson's ratio 0.25 (Lin & Stein, 2004), and friction coefficient 0.1 (Lamb, 2006). The receiver fault orientation was prescribed to be identical to the plate interface model.

Relation Between Coseismic Slip, Afterslip, and Aftershock Seismicity
The regions that experienced coseismic slip collocate with decreased Coulomb stress. They are surrounded by two lobes of elevated Coulomb stress updip and downdip of the two coseismic ruptures. Both of these lobes extend primarily from~20°S southward, although the shallower branch also extends northward to the location of the mainshock epicenter ( Figure 6). South of 20°S, the updip distribution of seismicity correlates well with the region of elevated Coulomb stress, whereas parts of the shallow seismicity directly updip of the mainshock's highest coseismic slip appear to have occurred in areas that experienced negligible Coulomb stress change due to an apparent shallow slip patch in the Duputel et al.'s (2015) model (Region R1; Figure 6). However, the near-trench part of subduction earthquake slip models is notoriously underconstrained, and published slip models for the Iquique earthquake also show significant variability here. For example, the models of Liu et al. (2015) and Gusman et al. (2015) show a similar region of shallow slip as Duputel et al.'s (2015) Figure S12 for comparison of slip models). The most sensitive data to shallow slip are offshore DART buoy tsunami height measurements, which are incorporated in Duputel et al.'s (2015) and Gusman et al.'s (2015) models, but also in the model of An et al. (2014). ΔCFS on the megathrust is strongly linked to the assumed slip model, and hence, interpretations of details in the pattern should be made with due caution considering these uncertainties. What seems to be robust is a general pattern of increased ΔCFS updip and downdip of the mainshock rupture in regions of high aftershock seismicity.  Although it could not be resolved geodetically, this shallow region of the megathrust, which extends southward along the shallower lobe of positive ΔCFS, might have undergone significant afterslip, driven by the static stress increase due to the mainshock and its largest aftershock. As we move further south, the observation of decreasing interplate locking (Li et al., 2015; Figure S6) leads to the expectation of afterslip after a large event; under partial locking, creep must occur during the interseismic phase, making the region also susceptible to coseismic stress changes. This spatial correlation between regions of reduced interseismic coupling, significant afterslip, and increased Coulomb stress has been previously observed in the areas surrounding main ruptures associated with the 2005 M8.7 Nias-Simeulue (Chlieh et al., 2008;Hsu et al., 2006) and the 2015 M8.3 Illapel (Huang et al., 2017;Tilmann et al., 2016) earthquakes, which could point to a similar explanation. Consequently, we interpret our results as indicative of an upper portion of the megathrust that predominantly undergoes aseismic afterslip. The occurrence of repeating earthquakes there may be another indication that this part of the fault is dominated by stable slip except where the seismic clusters and streaks occur (see sections 4.2 and 4.3). To a certain degree, the two deepest seismicity clusters C3 and C4 may resemble the shallower seismicity streaks in fault mechanics and slip properties ( Figure S6; see also section 4.2). They occur close to the upper plate Moho (Figure 2b, Profiles BB′ and DD′) and thus where the downdip transition to aseismic behavior was hypothesized to occur, except for fast and young plates, due to the contact of the lower plate with a serpentinized mantle wedge (e.g., Oleskevich et al., 1999). In other subduction zones, for example, Sumatra, however, the seismogenic zone has been shown to extend into the forearc mantle, and the mantle wedge does not always seem to be serpentinized (Dessa et al., 2009;Simoes et al., 2004).
We observe more complicated seismicity patterns in the region downdip and south of the mainshock main rupture, where aftershocks located in Regions R2 and R3 form a belt of seismicity located between two nearly aseismic patches of high afterslip (> 0.4 m), which are found in the vicinity of the M7.6 epicenter in the south and NE of the mainshock in the northern part of the study region ( Figure 6). This belt of aftershocks locates in a region of increased Coulomb stress, significant locking (0.7-0.9) and moderate afterslip (0.3-0.4 m), possibly dominated by conditionally stable frictional behavior. However, considerable seismicity can also be observed in areas of negligible or negative ΔCFS. Near R2 high afterslip overlaps both with the coseismic rupture of the M7.6 aftershock and with high locking (>0.9), partly accompanied by aftershocks, partly not. This region is difficult to interpret in detail and maybe best explained by a velocitystrengthening behavior with small-scale frictional heterogeneity that is not well represented in the overly smooth coseismic, interseismic, and postseismic models as well as in the resulting ΔCFS calculations. High frictional heterogeneity seems to be a hallmark of the deepest part of the megathrust expressed, for example, in localized coherent short-period seismic radiation from the mainshock rupture (Lay et al., 2012) and aftershock accumulation (Schurr et al., 2012). In fact, the deepest interface seismicity was observed to closely coincide with the short-period emission points in the 2010 M8.8 Maule earthquake aftershock time series . The northern patch of high afterslip observed north of 19.5°S has slightly reduced CFS in our model due to some coseismic slip. This region may be conditionally stable, where coseismic slip may propagate to but that behaves as velocity-strengthening at interseismic low strain rates. This is corroborated by a moderate degree of interseismic locking observed in this patch (<0.7; Figure S6).
Downdip of the belt of aftershocks formed by R2 and R3, we see a region of very sparse aftershock seismicity at~40-to 50-km depth, approximately delineated by the coastline (Figures 1 and 2a). A similar along-dip seismicity gap, albeit more pronounced, was observed in the aftershock series of the 2010 M8.8 Maule earthquake  as well as before and after the 2015 M8.3 Illapel earthquake and has been suggested to be a general pattern of the megathrust seismicity in central Chile (Lange et al., 2016). The aseismic area separating the shallower and deeper clusters has been shown to approximately correlate with the continental Moho in south-central Chile (Lange et al., , 2016. Our relocated seismicity confirms the occurrence of an aseismic zone, suggesting that this particular depth segmentation of the megathrust seismicity may be generic along most of the Chilean subduction zone. However, in the northern Chile case, the continental Moho (Sodoudi et al., 2011) is imaged at a level below even the deepest interplate events, downdip of where the seismicity gap occurs (Figure 2b, Cross Section BB′).

Linear Seismicity Streaks on the Shallow Plate Interface as Markers of Aseismic Creep
A characteristic feature of our catalog is several~EW striking seismicity streaks (Figure 2). Similar alignments of microearthquakes have been observed in continental strike-slip systems like the San Andreas fault (e.g., Rubin et al., 1999), where they were interpreted as indicative of small regions of brittle failure in a predominantly creeping section of the fault. Since these streaks are located offshore, that is, outside the seismometer network, the event-station geometry is unfavorable, and location errors in EW direction are nonnegligible ( Figure S2). In order to check whether these alignments are not artifacts of the location process, we performed a number of tests using travel time differences of P and S phases from the catalog as well as cross-correlation (CC) time delays between event pairs of streak earthquakes. In particular, we analyzed streaks labeled S1 and S2 in Figure 2a  We calculated differential S-P travel times for event pairs (e1 and e2) detected at the same station, given by and plotted them against interevent distances (Figures 4b and S5b). dt (S-P) should be independent of a possible trade-off between event location and origin time, thus is a suitable proxy for the distance between event and station. For this calculation, we used travel times to stations within a small angle (10°) around the streak, so that the azimuth of the seismic rays connecting source and receiver is relatively uniform. We compared the observed linear trend to synthetically calculated dt (S-P) values for a perfectly horizontal chain of earthquakes, an alignment along the slab surface as well as a vertical alignment. While we can rule out the last of these geometries, we cannot distinguish between the two other scenarios based on our data. However, since most available focal mechanisms for intrastreak or neighboring earthquakes consistently show low-angle thrusting consistent with displacement along the plate interface (Figures 3 and 7c), we prefer to infer that these streaks are not purely horizontal but oriented along the megathrust.
Following Rubin et al. (1999), we also created plots of arrival time differences (dt) between event pairs in the streak registered at one station against the angle between interevent vector and recording station (azimuth), normalized by the interevent distance d divided by the seismic velocity V in the source region. In this calculation, we used P and S time delays for all the stations located between 19°S and 22°S that were calculated from CC, and we complemented those with corresponding delays obtained from catalog picks when a threshold CC coefficient (CC > 0.5) was not reached for the CC delays. In case of a true linear arrangement, the observed pattern of time delays should resemble a cosine function, with maximum positive time delay at an azimuth of 0°, zero time delay at 90°, and maximum negative delay at 180°. Figures 4c and S5c show such plots for the two main streaks in our data set, which suggest that the seismicity streaks we observe are most likely real features.
We additionally searched for repeating earthquakes in our data set, that is, sequences of very closely located seismic events with highly similar waveforms (CC > 0.95 at at least three stations out of a set of five: PB01, PB02, PB08, PB11, and PB12). The cross correlation utilizes a 35-s time window that includes both P and S waves and applies a 1-to 4-Hz bandpass filter). We found a total of 211 such sequences containing at least two events occurred since Day 17 after the mainshock. The large majority of these sequences is situated in the shallow part of the plate interface, updip of the highest coseismic slip (Region R1), and within the major EW oriented seismicity streaks, although some are also found in various of the minor lineaments of microearthquakes that surround the M7.6 aftershock epicenter as well as in the most distinct (C2, C3, and C4) and some minor clusters located at the deepest part of the plate interface ( Figure 6). Due to the fact that waveform similarity between two earthquakes increases the more closely located they are in space, CC coefficients for pairs of earthquakes along the main streaks show a tendency to decay with the distance between event epicenters (Figures 4d and S5d); thus, repeaters within these streaks lie in the uppermost left part of these plots.
The observation of numerous repeating earthquakes occurring inside seismicity streaks on the shallow portion of the megathrust most likely hints at the presence of afterslip (e.g., Huang et al., 2017;Igarashi et al., 2003;Meng et al., 2015), which loads a number of velocity-weakening patches on the interface, that then repeatedly rupture seismically (Nadeau & Johnson, 1998;Schaff et al., 1998). In this sense, seismicity streaks would be markers of surrounding creep along the plate interface. A number of repeaters were also found where the deepest onshore clusters C3 and C4 are located, in a region of moderate afterslip. This could suggest similar frictional properties in the deepest and shallowest parts of the plate interface. As previously seen by Waldhauser et al. (2004) in a strike-slip environment, the seismicity streaks observed in our catalog occur along the edges of the slip patches (both for the mainshock and the M7.6 aftershock), thus in regions with significant slip gradient. From this observation, and following Bennington et al. (2011), we speculate that such linear features may also be indicators of a change in the mechanical properties of the fault and consequently could have played an important role in limiting the updip and along-strike extents of both the Iquique mainshock and its largest aftershock ruptures.
If the streaks were caused by rock destruction induced, for example, by subducting seamounts, we would expect them to be subparallel to the plate convergence vector (76°, Angermann et al., 1999), coseismic slip vectors (253.7°direction of movement of the hanging wall, average between mainshock and M7.6 aftershock;  Angermann et al., 1999), the rake (gray; movement of the hanging wall) of the mainshock and M7.6 aftershock focal mechanisms , the afterslip (green; Hoffmann et al., 2018), and rupture direction estimates for repeating events that occurred before (red) and after (blue) the Iquique earthquake (Folesky et al., 2018). The red box outlines a region surrounding the major streaks, shown in subfigure c. Interface aftershocks are plotted as gray circles. (b) Rose diagram of directions plotted in the map of subfigure a. Peripheral numbers are direction azimuths measured in degrees clockwise from the north (0°), and radial numbers are counts of afterslip directions and augmented counts of rupture directions (x4, for enhancing visibility). Gray arrow points toward the average (253.7°) rake direction between the mainshock and the M7.6 aftershock. (c) Map view zoomed in the region of the major streaks (labeled S1 and S2). Focal mechanisms collected in this area are shown as colored beachballs (as in Figure 3). Interface foreshocks and aftershocks and focal mechanisms are scaled by magnitude as shown in inset on the left. (d) Rose diagram for the distribution of observed rupture directions in the region of subfigure c. The correlation between orientations of the major streaks and the preferential~EW rupture direction observed in subfigure b is maintained at this local scale. Hayes et al., 2014), or afterslip vectors (250-260°;Hoffmann et al., 2018). However, the~EW orientation of the seismicity streaks deviates slightly but consistently from all of these (Figures 7a and 7b). Recently published rupture direction estimates of Iquique earthquake aftershocks (Folesky et al., 2018) appear to strike also predominantly toward east, that is, in line with the streaks (Figures 7c and 7d). This could imply that streaks are caused by a local persistent mechanism, rather than related to coseismic or postseismic processes.

Role of Forearc and Incoming Plate Structures
The northern Chile convergent margin is sediment starved (<500 m of sediment thickness) due to the aridity of the Atacama desert and has been classified as erosive (von Huene & Scholl, 1991). Erosive-dominated margins generally have continental wedges composed of two segments: (1) an outer wedge at the most seaward part, composed mainly of framework rock of the upper plate and with a small frontal prism of slump debris derived from land and possibly some accreted seafloor sediments, and (2) an inner wedge further landward made of more consolidated framework rock (von Huene et al., 2004;Wang & Hu, 2006). Small frontal prisms formed of mostly sedimentary debris eroded from the continental slope have been suggested throughout northern Chile based on reflection seismics, bathymetry, and gravity data:~5 km wide at~23°S, the latitude of the Mejillones Peninsula (Sallarès & Ranero, 2005);~30 km wide at 22°S (Contreras-Reyes et al., 2012);~10-20 km wide at the latitudes of the Iquique earthquake rupture ; and on averagẽ 20-30 km wide, estimated from Arica to Mejillones Peninsula (Maksymowicz et al., 2018). We observe that the updip limit of seismicity (foreshocks and aftershocks) seems to correlate well with the approximate location of the limit between lower and middle continental slopes (MLS, Maksymowicz et al., 2018; Figure 1), which thus marks the boundary between a seaward aseismic (20-30 km wide) and a landward seismically activated portion of the plate interface (Figures 6 and 8a; Maksymowicz et al., 2018). The gradual increase in coseismic slip landward of this boundary (Figure 8b) also supports the idea of a small aseismic frontal prism. It has been suggested that this frontal prism is composed of highly fractured and hydrated rocks and acted as a barrier for the updip propagation of the coseismic deformation during the 2014 Iquique Maksymowicz et al., 2018), the 2010 Maule and the 1960 Valdivia earthquakes . A transitional zone may be present landward of the MLS, where seismicity rate and coseismic slip increase landward, and the degree of fracturing smoothly decreases in the same direction (Maksymowicz et al., 2018;Wang & Hu, 2006).
In the context of dynamic Coulomb wedge theory (Wang & Hu, 2006), we envisage that the frontal prism and the transitional region form the outer wedge, overlying an updip velocity-strengthening and a transitional (frictionally heterogeneous) segment of the plate interface, respectively. The inner wedge is located further landward, overlying the velocity-weakening (seismogenic) segment of the megathrust, where we observe the highest coseismic slip (Figure 8a). During the Iquique earthquake, the outer wedge was elastically compressed and brought to or beyond the critical state, with an increase in basal and internal stresses and pore fluid pressure. After the earthquake, once the seismogenic segment is again locked, the outer wedge is expected to creep in order to relax the coseismically induced stress, returning to a stable state.
The concentration of aftershocks located in the postulated transitional region updip of the main coseismic rupture (Region R1; Figure 6) may be mostly driven by afterslip due to postseismic stress relaxation, in a region of the plate interface with dominant conditionally stable behavior. Given the higher uncertainties of the shallow coseismic slip, it is also plausible that afterslip here occurs in small-scale velocity-

10.1029/2019JB017794
Journal of Geophysical Research: Solid Earth strengthening ("stable") areas hidden in the stress shadow of the nearby velocity-weakening coseismic asperity. In this case, the stable areas would not slip interseismically and thus would appear with high degree of interseismic locking (>0.9; Figure S6); however, they would be able to undergo considerable postseismic slip (Bürgmann et al., 2005;Hetland & Simons, 2010;Métois et al., 2012). The significant concentration of NS oriented thrust faulting aftershocks with dip angles considerably steeper than the plate interface in this region (Figure 3) may indicate additional activation of offshore upper plate faults, associated with the outer wedge deformation due to its postseismic stress relaxation.
Stress accumulation caused by slow-slip prior to the mainshock has been suggested to account for the concentration of foreshocks in the same region on the plate boundary (Kato & Nakagawa, 2014;Maksymowicz et al., 2018;Yagi et al., 2014). And a number of thrust mechanisms, with significantly rotated planes relative to the plate interface (like the largest foreshock, M6.6; Figure 3), has additionally been interpreted as indication of some degree of preseismic upper plate deformation . Both interplate and upper plate events (foreshocks and aftershocks) form the seismically active wedge volume observed here ( Figure 8a). However, this spread seismicity also reflects very poor event depth resolution, as the depth uncertainties in this area are the highest ( Figure S2).
The incoming Nazca plate features substantial relief, both bending-related horst-and-grabens and spreading fabrics, as well as single prominent seamounts associated with the Iquique Ridge offshore our study area (Figure 1; . Rugged subducting seafloor has been associated with low interplate locking and aseismic creep on megathrusts (Wang & Bilek, 2014), and we further believe that it may have also had a strong influence on some characteristic features observed in the aftershock pattern of the Iquique earthquake. In addition to seamounts visible in the bathymetry, several currently subducting seamounts have been postulated based on the analysis of reflection seismics ( Figure 6; Geersen et al., 2015). Further indications of underthrusting lower plate relief come from large embayments and NW-SE trending antiforms and synforms in the lower continental slope seen in high-resolution bathymetry (Figures 2a and 6; . Embayments are observed in the northern part of the study region, partially at the location of the northernmost proposed seamount (Figure 6). Antiforms and synforms, located southward of 20°S and westward of the shallowest aftershock seismicity, have been interpreted as possible expressions of creep on the plate interface and internal deformation of the upper plate, related to subducting lower plate relief . A second suggested seamount appears to correlate with where we see the major seismicity streaks in our relocated earthquake catalog, whereas the southernmost such seamount locates at the position of the large offshore seismicity cluster C1 in the south, where we did not find repeating earthquakes (Figures 2 and 6). Focal mechanism solutions of the events in the earthquake cluster C1 are significantly more diverse than elsewhere (Figure 3), which indicates complex deformation involving the upper plate, possibly caused by increased stress due to the nearby M7.6 aftershock event activating a network of fractures that may have been generated by the subducting seamount. Correlations regarding subducting seamounts, however, should be considered with care. The positions of seamounts in Geersen et al. (2015) are prescribed by the three analyzed active seismic transects; seamounts may easily also be present elsewhere, where they could not be imaged. Further, the proposed seamounts were suggested based on images of two-way travel time, which can either be explained with plate interface topography or with strong lateral variations of seismic wave speed. As observed in the region updip of the main coseismic rupture, several thrust events, which are oriented~NS and feature dips steeper than the slab surface, likely also indicate the activation of offshore upper plate faults, possibly splay faults forking off from the megathrust (Profile DD′, Figure 2b, may display such a geometry). Even though location uncertainties are largest in this shallowest segment of the megathrust ( Figure S2), this interpretation is supported by the occurrence of two thrust events with rotated planes and steeper dip angles within the suspected splay ( Figure 3). Furthermore, a similar splay was observed in the aftershock sequence of the 2010 Maule earthquake in central Chile based on ocean bottom data (Lieser et al., 2014), that is, not subject to the same location uncertainty, such that aftershock activity on a large splay fault is certainly plausible.
Recently, as has been commonly observed at exhumed fault surfaces (Engelder, 1974;Petit, 1987), a 3D seismic reflection study of the Costa Rica margin has shown for the first time well-defined "corrugations" on a subduction zone megathrust (Edwards et al., 2018). They are expressed as regular shallow plate interface relief with a dominant wavelength, and an orientation deviating slightly from the plate convergence direction. From higher reflection amplitudes along the corrugated portions, Edwards et al. (2018) inferred a higher fluid content/pressure there and suggested that corrugations may act as fluid conduits. These features were seen in a region where seamounts, plateaus, and ridges are being subducted, which has been shown to cause extensive fracturing and faulting of the upper plate (e.g., Sak et al., 2009;Zhu et al., 2009), and also appears to inhibit seismic rupture propagation (DeShon et al., 2003;Wang & Bilek, 2011). Subducting rough topography and resulting fracturing and modification of the plate interface in the direction of convergence have also been linked to the formation of linear seismicity features observed offshore the Aleutian-Alaska subduction zone (Abers et al., 1995). Structural features aligned to the fault slip direction, such as corrugations or blocks of resistant host rock, have also been proposed as possible causes of seismicity streaks in strike-slip faulting environments (Rubin et al., 1999).
We envision that the megathrust offshore Iquique is not only being deformed and fractured by the effect of underthrusting rough relief, but it may also be corrugated as imaged in the Costa Rica subduction zone. Numerous fractures produced by subducting seamounts may allow the escape of fluids from the overpressured oceanic crust to the plate interface, breaking the sealed plate boundary (Audet et al., 2009). Similar to a proposed possible mechanism for generating slip-parallel streaking of tremors in the Cascadia subduction zone (Ghosh et al., 2010), the released fluids at the subduction interface may flow through corrugated conduits along the megathrust. By altering the effective normal stress due to the pressure exerted on the conduit walls, fluids may eventually cause shear failure on the interface and trigger earthquakes along the conduits, forming seismic structures such as the~EW striking streaks in our relocated catalog. Some of these ruptures may give rise to the repeating earthquakes along the streaks that we observe, in case a single fault patch is sequentially activated several times. . Synoptic model summarizing the proposed frictional heterogeneity along the plate interface and the seismotectonic segmentation of the continental wedge in the study area. The shallowest portion of the plate interface is aseismic (velocity-strengthening), underlying a highly fractured and deformed frontal prism (FP) in the outer wedge. Further downdip, a transitional zone is characterized by seismically unstable (velocity-weakening, along seismicity streaks) and conditionally stable (updip of mainshock asperity) regions embedded in aseismically creeping (velocity-strengthening) areas. Here fracture networks caused by subducting seamounts and repeating earthquakes hint at afterslip loading unstable/conditionally stable regions. Above this segment, the outer wedge forms a transitional zone (TZ) approximately below the middle continental slope, where fracturing of the upper plate likely gradually decreases landward. Deeper downdip, the seismogenic (velocity-weakening) segment of the plate interface correlates with the highest coseismic slips, below the inner continental wedge. This seismically unstable region changes in the downdip direction to a frictionally heterogeneous behavior, where we observe an anticorrelation between aftershock seismicity (presumably located in conditionally stable areas) and the two (velocity-strengthening) areas with the highest afterslip. Finally, seismicity clusters are observed in the deepest portion of the plate interface. Features not described in the legend are plotted as in Figures 2 and 6.
In section 4.2, we noted that the streaks are aligned with a preferential eastward rupture direction. Folesky et al. (2018) speculated that directivity may be caused by the bimaterial nature of the subduction process, which predicts downdip rupture direction, due to the higher compliance of the underlying plate material. They estimated that this effect (Weertman, 1980) would be enhanced due to the geometrical limitation of possible rupture directions given the along-dip elongated shape of the observed repeater asperities. We believe that the existence of corrugations is a plausible structural explanation for the observed~EW oriented seismicity streaks and thus for the associated elongated repeater asperities. Besides, the corrugations may behave as a pathway along the megathrust through which favored rupture directions could be reinforced.
The presence of corrugations is suggested by the earthquake mechanisms observed in the region of the major streaks. Although most of them indicate~NS striking thrust faulting occurring at the plate interface, a number of them are small (M ≤ 4) thrust events with dips steeper than the slab surface (Figures 4a, S5a, and 7c), indicating local deviations in the rupture planes that could reflect structural features such as corrugations. Some shallow events located at the edges of both major streaks, with thrust focal mechanisms striking along planes rotated with respect to the plate interface orientation, may additionally hint at upper plate deformation taking place in areas surrounding the observed streaks and presumed megathrust corrugations (Figures 4a, S5a, and 7c). Nevertheless, the scarce number of focal mechanisms in this region is far from being sufficient, and additional observations are required to confirm the existence of such features.

Conclusions
We conducted a detailed seismicity study focused on the aftershock sequence following the 2014 M8.1 Iquique earthquake. Our results allow us to illuminate small-scale heterogeneity on the plate interface and the seismotectonic behavior of the overlying continental wedge. We present our key results in a conceptual model in Figure 9. The continental forearc in the study region can be divided into an outer and an inner wedge. The outer wedge is composed of (1) a highly fractured and deformed frontal prism (FP) overlying the shallowest aseismic segment of the plate interface and (2) a transitional region (TZ) where fracturing of the upper plate gradually decreases landward, located above a transitional frictional segment of the megathrust. The inner wedge overlays the seismogenic segment of the megathrust, where we observe the highest coseismic slip (Figures 8a and 9).
The FP may have acted as a barrier for the updip propagation of the mainshock rupture. Further downdip, the plate interface is characterized by velocity-weakening and conditionally stable regions embedded in a primarily velocity-strengthening area, where aseismic afterslip driven by increased static Coulomb stresses takes place. Afterslip is further favored by subducted seafloor relief, which may impose conditions such as fracturing of plate interface rocks, upper plate deformation, and reduction of plate coupling. The presence of earthquakes with rupture planes significantly different from those expected for megathrust events may indicate upper plate deformation and the activation of splay faults. A prominent cluster of seismicity (C1) may be linked to the fracture network developed by a subducting seamount. Numerous repeating earthquakes updip of the main coseismic rupture and inside seismicity streaks support the occurrence of afterslip, despite it not being detectable geodetically (Figure 9).
We interpret the earthquake streaks and their repeating earthquakes as markers of creep that takes place in the surrounding aseismic regions of the plate interface. The streaks may indicate a change in mechanical properties of the fault and consequently could have played an important role in limiting the updip and along-strike extents of the ruptures of both the Iquique mainshock and its largest aftershock. We speculate that the streaks may have a structural origin, possibly being triggered by shear failure induced by fluid transported along conduits formed due to presumed corrugations present on the plate interface.
In the central portion of the megathrust, aftershock seismicity surrounds the asperities of the mainshock and the M7.6 aftershock, which can mostly be explained by the Coulomb stress redistribution triggered by these events in a region dominated by mixed frictional behavior. Moreover, maximum afterslip is located in two patches flanking and below the main coseismic rupture, which correlate with areas of low aftershock seismicity. These areas are thus interpreted as dominated by a velocity-strengthening frictional regime (Figure 9.).
The southern low seismicity area connects downdip to a depth interval of sparse aftershock seismicity at~40to 50-km depth, which supports previous observations that hint at a possible depth segmentation of the interface seismicity along the Chilean margin. Downdip of this region, the deepest interplate earthquakes reach~55-to 65-km depth in two individual clusters beneath the coastal cordillera, which feature a significant difference in dip from each other. At this depth, the observation of repeating earthquakes and moderate afterslip could imply similar frictional conditions as proposed for the shallow megathrust transitional zone (Figure 9.).