Linear Relationship Between Aftershock Productivity and Seismic Coupling in the Northern Chile Subduction Zone

The aftershock productivity is known to strongly vary for different mainshocks of the same magnitude, which cannot be simply explained by random fluctuations. In addition to variable source mechanisms, different rheological properties might be responsible for the observed variations. Here we show, for the subduction zone of northern Chile, that the aftershock productivity is linearly related to the degree of mechanical coupling along the subduction interface. Using the earthquake catalog of Sippl et al. (2018, https://doi.org/10.1002/2017JB015384), which consists of more than 100,000 events between 2007 and 2014, and three different coupling maps inferred from interseismic geodetic deformation data, we show that the observed aftershock numbers are significantly lower than expected from the Båth's law. Furthermore, the productivity decays systematically with depth in the uppermost 80 km, while the b value increases. We show that this lack of aftershocks and the observed depth dependence can be simply explained by a linear relationship between the productivity and the coupling coefficient, leading to Båth law only in the case of full coupling. Our results indicate that coupling maps might be useful to forecast aftershock productivity and vice versa.


Introduction
An important empirical feature of aftershock sequences is the so-called Båth's (1965) law, which states that the magnitude difference Δm between a mainshock and its largest aftershock does not depend on the mainshock magnitude M and is on average close to Δm = 1.2. Since earthquake magnitudes are known to generally follow the well-known Gutenberg-Richter distribution and are essentially not correlated to the earthquake activity rate, the expected value of the maximum aftershock magnitude is related to the total number of triggered aftershocks. Based on the Gutenberg-Richter distribution, one aftershock of M − Δm or larger is expected if the number of aftershocks above the completeness magnitude M c is (1) However, both Δm and the productivity N a for equal mainshock magnitudes are found to vary more strongly than the random variability expected from a Poisson process (Marsan & Helmstetter, 2017;Page et al., 2016). The reason might be the variability of the mainshock source mechanisms and rheological properties. In particular, it has been shown by Bouchon and Karabulut (2008) that supershear ruptures lead to depleted on-fault aftershocks, which likely affects the total number of triggered aftershocks. Another source effect is proposed by Marsan and Helmstetter (2017), who analyzed the aftershock productivity in southern California and suggested a relation of the observed variability with the stress drop variability of the mainshocks. It can also be expected that further source characteristics such as the slip roughness affect the number of short-term aftershocks (Helmstetter & Shaw, 2006;Marsan, 2006). In addition to source effects, the stress state and rheological properties of the source region may be influential. For example, a dependence on the faulting style has been recognized with highest productivity values found for thrust events (Tahir et al., 2012;Tahir & Grasso, 2015), which might be related to higher background stresses. Finally, increased heat flow is expected to favor viscous instead of brittle deformation (Ben-Zion & Lyakhovsky, 2006), which is in agreement with decreased aftershock productivity found in areas of increased heat flow in southern California (Yang & Ben-Zion, 2009). The latter explanation is also in accordance with the observation that oceanic transform faults trigger relatively weak aftershock sequences (Boettcher & Jordan, 2004;McGuire et al., 2005), indicating that seismic coupling may be an important factor to explain productivity changes. The seismic coupling coefficient C refers to the fraction of the seismic slip related to the total slip on faults (Pacheco et al., 1993). In fault zones without any aseismic processes, all permanent deformation is related to earthquakes, and C is one. In contrast, in purely creeping regions, C is zero. In general, seismic coupling is inferred by the ratio between the recorded average seismic moment rate and the geodetic moment rate (Brune, 1968;McCaffrey, 1997). However, in regions dominated by a single fault such as subduction zones, it can be directly inferred from geodetic data by inverting backslip on the fault interface during interseismic periods assuming that deformation is limited to the interface. The C value is then evaluated by the locking degree of the interface (Savage, 1983). Geodetic observations of interseismic and postseismic deformation indicate heterogeneous frictional properties with seismically slipping (velocity weakening) asperities embedded in a creeping (velocity strengthening) environment (Lay et al., 2012). C can be thus also seen as an approximation of the fraction of the velocity weakening area on the fault. Based on these considerations, the aftershock productivity can be expected to be directly proportional to C, because the density and number of potential aftershock sites increases with C as illustrated in the schematic Figure 1. The same relation can also be expected for 3-D-distributed deformation, where C represents the fraction of the volume with brittle rheology. While this is an intuitive expectation, it has not been proven for observational data so far.
In this study, we analyze a recently published 8-year-long earthquake catalog of the northern Chile subduction zone with event magnitudes down to two, which comprehensively maps upper and interplate deformation and outlines structures within the downgoing Nazca plate in unprecedented detail. For the same region, coupling maps based on geodetic measurements have been recently published by different groups. This offers the opportunity to analyze the aftershock productivity in relation to seismic coupling in detail. Our results presented in the remainder of the paper indicate that both quantities are indeed linearly related.

Data Selection
We use the recently published catalog of Sippl et al. (2018), which is based on continuous recordings acquired in the framework of Integrated Plate Boundary Observatory Chile IPOC (http://www.ipoc-network.org/; GFZ German Research Centre for Geosciences; Institut des Sciences de l'Univers-Centre National de la Recherche CNRS-INSU, 2006) and contains more than 100,000 double-difference relocated earthquakes with absolute location errors of approximately 5 km and relative errors of about 2-3 km within the first 100-km depth. The spatial distribution of events is shown in Figure 2. Here we use the classification of Sippl et al. (2018), differentiating between the three bands of the downgoing slab, namely, interface (P1), upper band (P2), and lower band (P3) events, as well as intermediate depth (ID) and the upper plate activity (UP). The remaining events are in the category NN. Figure 3a shows the frequency-magnitude distribution for the whole data set as well as for each subset separately. The distributions can be well fitted by Gutenberg-Richter distributions considering a detection probability which is given by a cumulative function of the Normal distribution (Ogata & Katsura, 1993, dashed curves). The parameters resulting from maximum likelihood fits are summarized in Table 1. For all subduction zone classes (P1, P2, P3, and ID), the b value is lower than one, with a smallest value of 0.58 for P1. Apart from the NN class, M c is found to be less than 2.7, which is used as conservative threshold for the following analysis.

Aftershock Selection
For the aftershock selection, we use a recently established scheme based on nearest-neighbor distances (Baiesi & Paczuski, 2004Zaliapin et al., 2008;Zaliapin & Ben-Zion, 2013). It is a purely statistical method which does not rely on any particular triggering mechanism, such as static/dynamic coseismic stress changes or afterslip. The method quantifies the correlation between an event i and a preceding event j by its magnitude-weighted space-time distance x, and M being the time, location, and magnitude of the events, respectively. d is the fractal dimension of the hypocenter distribution, which would be 2 for a planar distribution and 3 for a homogeneous three-dimensional distribution. Among all events j preceding i, the identification of the (most likely) trigger of i results from selecting that event with the lowest n ij value. To distinguish between triggered and background activity, a threshold value of n c is set and only events with n ij ≤ n c are considered as plausible mainshock-aftershock pairs. By means of epidemic-type aftershock sequence simulations, the applied detection method has been previously demonstrated to be robust with respect to (1) changes of the involved parameters of the method, (2) catalog incompleteness, and (3) location errors (Zaliapin & Ben-Zion, 2013).  Ogata and Katsura (1993) with parameters provided in Table 1. (b) The temporal decay for aftershocks in the subset P1 (points), P2 (squares), and P3 (crosses) for different mainshock magnitude ranges. Here the squares and crosses are shifted slightly in time to enhance visibility. The points and error bars refer to the mean and the 90% confidence interval related to alternative aftershock selection parameters, while the dashed black line refers to a Omori p value of 1.0. The numbers of mainshocks for the standard selection parameters Here we use hypocentral distances and consider alternative parameter selections and subsets in order to analyze the stability of our results. In particular, we are interested in the events on or near the plate interface which should be most affected by variations of coupling. Thus, we analyze the interface events (P1, b = 0.58) separately, as well as the combination (P12) of subsets P1 and P2 due to their spatial closeness and partial overlap (with b = 0.60) and, for comparison, the whole data set together (all, b = 0.77). Furthermore, we use a fractal dimension of d = 2.3, which is in agreement with the fractal dimensions for these events estimated by means of the correlation integral (Grassberger, 1983; see supporting information Figure S1). However, we test in all cases the alternative values of d = 2.1 and 2.5. The values for the threshold is log(n c ) = −2.3 (using units of kilometers and years) with a tested range between −2.7 and −1.9 for P1 and P12, while log(n c ) is set to −4.4 [−5.0, −3.8] for all events (see Figure S2).
The selection procedure leads to separated sequences of earthquakes connected by n ij ≤ n c . These consist of one earthquake (isolated background event) or more earthquakes (background event with directly or indirectly triggered subsequent events). For all resulting earthquake sequences, the largest event is defined as mainshock, and all following earthquakes are defined as aftershocks. For illustration, the Figure S3 shows Table 1 The Maximum Likelihood Estimates of the b Value and Completeness Parameters Using the Model of Ogata and Katsura (1993)  Note. and are the mean and standard value used in the cumulative Normal function describing the detection probability, while M c is the completeness magnitude related to a completeness level of 99%. The errors refer to the standard deviation derived from bootstrapping.
the temporal and spatial distribution of the P12 mainshocks which are selected by the standard parameters. The intense aftershock sequences of the two largest events (2007 M7.8 Tocopilla and 2014 M8.1 Iquique) cluster in space and time around the rupture zone (Hoffmann et al., 2018;Schurr et al., 2012Schurr et al., , 2014, while the activity related to mainshocks in the range between 4 and 6, which is analyzed in section 3, is almost evenly distributed. Note that many small mainshocks have no recorded aftershocks; nonetheless, their aftershock productivity (N a = 0) is taken into account when the average aftershock numbers and rates are calculated.
The stacked and averaged aftershock rates are shown as function of time and mainshock magnitude in Figure 3b. On average, the aftershock rates are found to follow the well-known Omori-Utsu function, which is a power law with a time-offset parameter c where t is the time relative to the mainshock and the exponent p is usually in the range 0.8-1.2 (Utsu et al., 1995). This is also valid in the case of the selected aftershocks which show a power law decay with p ≈ 1. The c value is mostly related to the short-term incompleteness of empirical catalogs. It has been shown that sophisticated reprocessing of the recorded seismograms leads to the detection of many early aftershocks which are missed by routine detection procedures (Kagan, 2004;Peng et al., 2007) and that the empirically observed delayed onset of the decay is in agreement with the expected incompleteness due to overlapping waveforms (Hainzl, 2016a(Hainzl, , 2016b. This short-term incompleteness is expected to be most prominent for large mainshocks (relative to the cutoff magnitude M c ), while it should be almost negligible for small-to moderate-sized mainshocks. This is also observed in our case, where the onset of the aftershocks for mainshocks with magnitude M ≥ 6 are significantly delayed, particularly in the case of the largest events, the 2007 M7.8 Tocopilla and the 2014 M8.1 Iquique mainshocks. In contrast, a delayed onset of the power law is not obvious for mainshocks with M < 6, indicating that the incompleteness issue is insignificant for mainshocks below M6.
In the following, we do not care about the timing of the aftershocks and only count the total number N a of aftershocks following each mainshock in the first 100 days. The cutoff time is introduced to further minimize the effects of a contamination by background events (see Supporting Information S2). We also test the robustness of our results by a shorter cutoff T max of 10 days. Furthermore, we focus on the small mainshock magnitudes. Besides the completeness issue mentioned above, the reason is twofold: First, smaller mainshocks are more frequent, which leads to more robust results, and second, smaller mainshocks have a smaller rupture dimension and thus provide more localized information.

Coupling Data
Coupling maps are based on records of deformation acquired at spatially distributed locations during the interseismic period. The corresponding uncertainties of the coupling values are rather high. To account for those uncertainties, we use three alternative coupling maps which have been previously published for Northern Chile by different research teams using different inversion schemes. In particular, we use the coupling map of Béjar-Pizarro et al. (2013)  In the following, we restrict our analysis of the coupling data to the range between latitude 19.0 • S and 23.5 • S, which is the range of the recorded slab-related seismic activity in subsets P1 and P2. The three coupling maps are shown in Figure S4 together with the recorded seismicity. The comparison shows that the maps have similar first-order features but significant differences in details.
The depth-averaged coupling values are approximately constant as function of latitude, while the trench-parallel averaged values strongly depend on depth. In particular, the coupling values at the mainshock locations vary strongly for a given latitude but are well constrained for a given depth value (see Figure  S5). Therefore, we also focus on trench-parallel averaged coupling values later on.

Event Density
The time-averaged total slip rate is expected to be constant for the subduction megathrust, because the sum of the seismic and aseismic motion should be equal to the convergence rate. Because the coupling coefficient is the ratio between the seismic and total slip, the average seismic moment release is expected to be proportional to the coupling coefficient. Such a linear relation should also hold for the number of earthquakes, if the frequency-magnitude distribution is constant in space. However, the event density is found to be only partly correlated to the depth dependence of the coupling values. This is shown in Figure 4b, where the number of m ≥ 2.7 events in subset P12 is shown as function of the interface depth. For that purpose, we first projected the earthquake hypocenters to the plate interface. The event distribution is compared to the average coupling in the corresponding depth level, which was first calculated for each coupling map separately. The peak of the event distribution is clearly shifted to deeper depth values compared to that of the coupling data. This discrepancy can be at least partly explained by an increase of the Gutenberg-Richter b value with depth, which is shown in Figure 4a. Here symbols refer to the maximum likelihood estimation of the b value (Aki, 1965) in nonoverlapping depth bins of 5 km. The result can be approximated by the linear dependence b(z) = 0.45 + 0.0042z with depth z in units of kilometers. The systematic increase of b with depth is confirmed by a recently developed Bayesian approach to detect statistical significant change points Thus, the average seismic moment release per event is smaller for larger b values at depth. Consequently, the peak of the seismic moment release is expected to be at smaller depth compared to that of the event number distribution. In principle, this could be tested directly by summing the observed seismic moments in different depth bins. However, the observation time is too short to sample the whole magnitude distribution, in particular the largest magnitudes which dominate the seismic moment release. Thus, the empirical distribution is dominated by random fluctuations. Instead, we calculated the expected mean seismic moment release per event based on the assumption of Gutenberg-Richter distributed magnitudes between M c and a fixed upper limit of M max = 8 (equation 13 of Zakharova et al., 2013). Here we account for the depth-varying b(z) and multiply, in each depth level, the observed number of events with the calculated mean seismic moment per event. The resulting depth dependence of the calculated seismic moment release is shown in Figure 4c. It is found to be close to that of the coupling coefficient as theoretically expected.

Aftershock Productivity
At first, we use the full spatial information of the mainshock location and associate each mainshock with the nearest coupling value of the analyzed coupling map. We use the whole data set but ignore mainshocks with distance larger than 10 km to the closest coupling value to ensure an appropriate association to the local coupling values. Then we sort the sequences according to their coupling value and calculate the mean aftershock number for bins of the coupling value with a bin width of 0.1. The corresponding uncertainties are calculated by repetitions for alternative aftershock selections with varying fractal dimensions d and threshold value n c (see section 2.1). In Figure 5, the final results are shown for mainshocks in the magnitude range between 4 and 5. The results for the different coupling maps show a large variability, which is related to the uncertainties of the coupling values inverted from interseismic deformation. However, the overall trend indicates that mainshocks occurring in areas of larger coupling C produce on average more aftershocks, while those in weakly coupled areas are almost lacking any aftershocks. The trend can be fitted by the linear relation N a = 0.7 C. and T max = 10 and 100 days. Note that the legend in (a) also holds for (b) and that a few of the points are out of the ordinate range for T max = 100 days (one for P1 and 5 for P12; see Figure S6 for the full range of results). The bullet points refer to the average aftershock number calculated in depth bins of 10 km in the case that all data are used. The error bars indicate the mean ± one standard deviation for alternative combinations of d and n c used for aftershock selection in this case and the small green points refer to the corresponding result for early aftershocks (T max = 10 days).
The coupling values are well constrained for given depth (see Figure S5). Thus, we now focus on the analysis of the relation between averaged aftershock productivity and coupling as function of depth, because this is expected to lead to more robust results due to the additional averaging in trench-parallel direction before comparison. For a proper comparison to the coupling depth values, we project each mainshock hypocenter to the plate interface and use the corresponding depth value of the interface as sequence depth. For the analysis of the largest data set (all), we used equidistant depth bins of 10 km. In this case, the uncertainties are again estimated by repetitions for different fractal dimensions d and threshold value n c . However, to ensure sufficient statistics in the case of the subsets P1 and P12, we use variable, nonoverlapping bins to include 40 or 10 mainshocks in the case of mainshocks with magnitudes 4 ≤ M < 5 and 5 ≤ M < 6, respectively. In the latter case, the results are associated with the average depth of the mainshocks in the corresponding window. The results are shown in Figure 6 for alternative aftershock selection parameters. The aftershock productivity shows a clear decaying trend with increasing depth for all data sets. While mainshocks in the range 4 ≤ M < 5 (average magnitude 4.3) have on average approximately 1.3 aftershocks Figure 7. Comparison of the depth dependence of the observed (symbols) and expected (lines) aftershock productivity based on the Båth law as function of the interface depth for mainshocks in the range (a) 4 ≤ M < 5 and (b) 5 ≤ M < 6. In both plots, the thin lines represent the expected value without consideration of the coupling coefficient, while the thick line refers to the theoretical value when the reduced and variable coupling is considered (cyan: T max = 10 days; black: T max = 100 days). Here the gray shaded area indicates the range of results (T max = 100 days) using the minimum and maximum values of the three alternative coupling maps. The symbols are the same as in Figure 6.
at a depth of 20 km, the number systematically decreases to less than 0.2 in a depth of 80 km. Similarly, a decay from an average of six to less than one aftershock is observed for the 5 ≤ M < 6 mainshocks (average magnitude 5.3) in the same depth interval. The observed depth dependence of the aftershock productivity can be compared to that of the coupling values. For that purpose, we calculated the average coupling value at a certain depth for each of the three coupling models independently. In Figure 6, the dashed line and the gray shaded area represents the average, respectively, the range between the minimum and maximum, of those three values with separated scale on the right of the plot. Also, the coupling coefficients show a strong decay from values of about 0.7 at 20-km depth to values around 0.1 at 60-km depth with a shape very similar to the productivity data. The simultaneous decay indicates a simple linear relation between coupling and aftershock productivity.
Note that we also repeated the same analysis for the isolated P3 subset, which represents the lower plane of the so-called double Benioff zone. The results show no clear correlation to the coupling coefficient (see Figure S8). This is reasonable, because this activity is not expected to be correlated (within the current mechanical knowledge) to the coupling between the sliding plates.

Relation to Båth's Law
Based on the empirical Båth's (1965) law, the average number of aftershocks N a is expected to be given by equation (1) as function of the b value and the difference between mainshock and completeness magnitude, M − M c . The Båth's law was previously derived for large mainshocks for which aftershocks were usually selected within longer time intervals; for example, Tahir et al. (2012) considered aftershocks within 1 year finding Δm = 1.2. In our case, we selected the aftershocks for our smaller mainshocks only in limited time intervals of T max = 10 or 100 days. Assuming an Omori law (p = 1), we can expect that only the fraction (1) occur within T max . Here we assume c = 0.001 days, which leads to (10) = 0.72 and (100) = 0.90.
As shown in Figure 4a, the b value is increasing approximately linearly with depth according to b(z) = 0.45 + 0.0042z. Therefore, the expected number of aftershocks should also increase with depth, which is indicated by the thin lines in Figure 7. However, the theoretically expected absolute value, as well as its depth trend, is in strong contrast to the observation. As seen in Figure 7, the expected aftershock activity significantly exceeds the observed one, and the aftershock productivity is expected to increase with depth z but actually decreases. A simple explanation for this discrepancy is that we have so far ignored the seismic coupling coefficient in this consideration. The coupling coefficient C represents the fraction of the total slip which is related to earthquakes. Assuming that the Båth's law holds for fully coupled faults, the aftershock numbers should be only a fraction C of the expected value provided in equation (1). This leads to an expected number of which is shown in Figure 7 as bold lines. The result is found to match the observations very well for both independent data sets of 4 ≤ M < 5 and 5 ≤ M < 6 mainshocks. This means that the observed depth-dependence and absolute values of the average aftershock productivity can be explained based on the empirical b-value trend and the geodetically inverted coupling values alone.

Discussion
A linear correlation between coupling, respectively, locking degree, and aftershock productivity can be explained based on simple considerations. The number of aftershocks is expected to scale with the stored available elastic energy in the area affected by the mainshock, if the rupture properties of aftershocks are constant. The locking degree is itself a measure of the elastic strain loading rate. Thus, areas with higher locking degree will have, on average, a higher level of stored elastic energy and can thus produce more aftershocks when affected by mainshock-induced stress changes. Considering the asperity model, the locking degree is connected to the fractional area of asperities, which represent the fully locked parts embedded in a creeping environment. In this case, elastic strain buildup is mainly limited to asperities, and the number of aftershocks is thus expected to scale with the fractional area of asperities (see Figure 1).
The observed correlation between aftershock productivity and seismic coupling is in agreement with previous findings for the aftershock zones of the Maule 2010 M w 8.8, South Chile, and the Tohoku-oki 2011 M w 9.0, Japan, megathrust earthquakes. Using a modified epidemic-type aftershock sequence model to identify secondary aftershocks (triggered by preceding aftershocks), Zakharova et al. (2017) found, among other things, a spatial correlation between the seismic coupling value and the secondary aftershock productivity.
Although their data analysis differs significantly both in the aftershock association procedure and the data type, the agreement of both studies indicate that the correlation between seismic coupling and aftershock productivity is not restricted to Northern Chile.
Additionally, the comparison of the absolute values of the aftershock productivity with the expectations related to the empirical Båth's law yields interesting insights. Our results show that the average aftershock productivity is in agreement with the aftershock activity associated to the Båth's law with Δm = 1.2 only in the limit of fully coupled fault zones. The average maximum aftershock magnitude is in general a function of the coupling coefficient and is, on average, M − 1.2 + log 10 (C)∕b, according to our results (equation (3)). However, this result is obtained so far only for one specific subduction zone, and it is an open question whether or not the relation also holds for different tectonic regimes.
It is interesting to note that we observe a clear decrease of the aftershock productivity with depth, while the average stress drop Δ increases with depth in subduction zones, strongest in the first 20 km (Lay et al., 2012). Based on simple considerations, Marsan and Helmstetter (2017) related the total number of aftershocks to the stress drop according to N a = Δ 1−d∕3 M d∕3 0 , where M 0 is the seismic moment of the mainshock. For the same moment magnitude and a constant fractal dimension of d ≈ 2.3, this relation would thus predict an increase of the aftershock productivity with depth, which is in contradiction to our observation. A possible explanation is that the fractal dimension is not constant and changes with depth, but our analysis of the correlation integral in different depth intervals does not indicate a significant change of d. Therefore, the stress drop effect on the aftershock productivity, if existing, seems to be completely masked in the first 60-km depth by the effect of the strongly decreasing coupling values. However, the increase of the aftershock productivity at depths around 80 to 120 km might be related to higher stress drops of these intermediate depth events (Derode & Campos, 2019).
Furthermore, our observed linear relationship between the aftershock productivity and the coupling degree might explain that slow slip events are less efficient in aftershock triggering. Pollitz and Johnston (2006) observed that slow slip events trigger significantly fewer aftershocks than earthquakes of the same seismic moment in the San Juan Bautista, California, area. The authors interpreted this observation as evidence that dynamic triggering plays a dominant role in aftershock triggering. However, slow slip events likely occur in fault sections with lower coupling, while earthquakes occur generally in more strongly coupled fault segments. Thus, even if only static stress triggering takes place, slow slip events are expected to trigger fewer aftershocks given the same amount of static stress transfer.
It should be noted that seismicity data have usually a higher resolution than geodetically inverted coupling data. Thus, small mainshocks might sample heterogeneities of the coupling which are not resolved in estimated coupling maps. To obtain comparable aftershock productivity estimations for small mainshocks, the aftershock productivity has to be averaged over many mainshocks, which randomly sample those heterogeneities. Vice versa, for larger mainshocks, the extent of the aftershocks might exceed the resolution length scale of coupling maps, and the inverted coupling values might vary within the aftershock zone. In this case, the total aftershock productivity will be related to the average coupling value in the whole region, where the mainshock induced positive stresses, rather than to the coupling value at the mainshock hypocenter location. Thus, the total aftershock number cannot be explained by a single coupling value for very large earthquakes. However, the local aftershock density is expected to scale linearly with the local coupling values at the aftershock locations in the case of large mainshocks, as indicated by positive correlation coefficients in the case of the aftershocks of the M w 8.8 Maule and M w 9.0 Tohoku-oki mainshocks (Zakharova et al., 2017).
Our results indicate that coupling maps could be potentially useful for aftershock forecasting. However, even in the case that coupling and aftershock productivity are in general linearly related as indicated by our case study for northern Chile, potential forecasts of future aftershock activity would depend on the assumption that coupling is stationary in time, which is questioned by recent findings (Loveless & Meade, 2016;Melnick et al., 2017).

Conclusions
The high-resolution seismicity data set of Northern Chile allows a detailed analysis of the seismicity and its relation to deformation on the fault. Our analysis of the activity in the vicinity of the plate interface shows that coupling, event rates, seismic moment release, b value, and the aftershock productivity are all depth dependent and clearly related to each other. The coupling is found to be linearly related to the total seismic moment release as well as to the aftershock productivity, which can be simply explained by considering coupling as the relative density of asperities (see Figure 1). Where asperities are denser, the coupling is higher, and more aftershocks can occur. Furthermore, we observe that the aftershock productivity is unexpectedly low in deeper parts of the subduction zone, if compared to expected values based on the empirical Båth's law. However, this apparent inconsistency vanishes, and observations are well fitted, if Båth's law is assumed to hold for fully coupled faults, while coupling strongly decays with depth in the subduction zone. All findings for northern Chile indicate a simple linear relation between coupling and aftershock productivity. If generally valid, our results consequently suggest that the average local aftershock productivity might be estimated from geodetically derived coupling coefficient and vice versa. All analyzed data have been previously published; in particular the analyzed earthquake catalog is available as electronic supplement of the manuscript by Sippl et al. (2018). We are grateful to Blandine Gardonio, the anonymous Associate Editor, and the anonymous second reviewer for their helpful comments and suggestions.