Comparison of postseismic afterslip models with aftershock seismicity for three subduction-zone earthquakes : Nias 2005 , Maule 2010 and Tohoku 2011

D. Lange,1 J. R. Bedford,2 M. Moreno,2 F. Tilmann,2 J. C. Baez,3 M. Bevis4 and F. Krüger5 1GEOMAR Helmholtzzentrum für Ozeanforschung, Kiel, Germany. E-mail: dlange@geomar.de 2Helmholtz-Zentrum Potsdam, Deutsches GeoForschungsZentrum GFZ, Potsdam, Germany 3Centro Sismológico Nacional, Universidad de Chile, Santiago de Chile, Chile 4School of Earth Sciences, Ohio State University, OH 43210, USA 5Universität Potsdam, Institut für Erdund Umweltwissenschaften, Potsdam, Germany

difficulty, only a few studies successfully combined afterslip and viscoelastic flow to explain time series of postseismic displacements with multiple mechanisms (Freed et al. 2006;Bruhat et al. 2011;Rousset et al. 2012).
The temporal decay of aftershock sequences follows the modified Omori's law (Utsu 1961): where n is the number of aftershocks in a time period t, t is the time after the main shock and K S , C S and p are constants. K S is a measure of the aftershock productivity of the sequence and the exponent p is usually between 0.9 and 1.5 with a median of about 1.1 (Reasenberg & Jones 1989;Utsu et al. 1995). The variable C S is often a fraction of a day (e.g. Utsu et al. 1995;Enescu et al. 2009) but it is unclear whether this has a physical reason or simply reflects decreased catalogue completeness early in the aftershock sequence (Utsu 2002). Integrating eq. (1) from 0 to t using p = 1 and C s > 0 gives the cumulative number of events since the time after the mainshock (Jonsdottir et al. 2006): (2) The postseismic displacement can be modelled using a rateand state-dependent friction law (e.g. Dieterich 1979Dieterich , 1994 and the displacement (D) versus time in the velocity-strengthening frictional regime can be approximated with (Marone et al. 1991;Wennerberg & Sharp 1997): where K G and C G are constants. Savage (2007) shows that the conventional, 1-D, spring-block model with a rheology compatible with ordinary transient creep leads to the same temporal decay of slip rate as described by eq. (3). This logarithmic law is thus commonly used to fit afterslip displacements (e.g. Hsu et al. 2006;Kreemer et al. 2006;Freed 2007;Savage & Svarc 2009). In case long postseismic or preseismic data are available eqs (2) and (3) can be expanded with an additional linear term (at) representing the interseismic deformation. However, in Chile interseismic GPS rates of coastal stations are of the order of 3-4 cm yr −1 , nearly an order of magnitude less than the displacement in 1 yr of postseismic deformation (Fig. S1), such that we opt to not include this term. Comparing eq. (2) with (3) shows that the cumulative number of events calculated from the modified Omori law results in time dependent functions of the same form but with possibly different parameters.
If surface displacements and aftershocks show the same decay with time, that is C S ≈ C G , the triggering processes should be mechanically linked. An agreement in the decay function would result in a linear correlation between the cumulative number of aftershocks and the cumulative displacement and such a correlation was indeed observed for continental strike-slip (e.g. Savage & Yu 2007;Murray-Moraleda & Simpson 2009;Wang et al. 2009;Yang & Ben-Zion 2009) and other subduction zone events (Perfettini & Avouac 2004a;Hsu et al. 2006;Ozawa et al. 2012). Wennerberg & Sharp (1997) show that the cumulative moment from seismicity is proportional to the cumulative number of aftershocks multiplied by the integral of the frequency-magnitude distribution assuming that the seismicity rate is high enough that the Gutenberg-Richter distribution is 'filled in' for any time interval. Hence, for a Gutenberg-Richter distribution which is constant in time the cumulative moment is proportional to the cumulative number of aftershocks which in turn is proportional to the accumulated strain release. An underlying assumption in this chain of arguments is that the nucleation time of the aftershocks governed by the time evolution of the stress is negligible (e.g. Hsu et al. 2006;Helmstetter & Shaw 2009).
Previous studies have shown that the total moment released by aftershocks is smaller than the total equivalent moment released by afterslip as constrained by geodetic observations. The equivalent total moment release from afterslip is often found to be around 20-30 per cent of the co-seismic moment release (Melbourne et al. 2002;Hsu et al. 2006;Chlieh et al. 2007;Perfettini et al. 2010;Ozawa et al. 2012), while the cumulative moment from aftershocks is generally much less (often around 1-10 per cent of the postseismic deformation) (Hsu et al. 2006; Barbot et al. 2009;Perfettini et al. 2010;Ozawa et al. 2012). Therefore, aftershocks were suggested as being triggered by the stress perturbations induced by the afterslip (Perfettini & Avouac 2004a;Savage 2010). An alternative is that seismicity groups in areas of stress concentration marking the transition between largely homogeneous seismogenic and aseismic regions . The distribution of aftershock activity also shows complex spatio-temporal behaviour that may correlate with zones of high postseismic strain (Das & Henry 2003). Often, the aftershocks cluster around the coseismic rupture area due to elevated differential stresses (Mendoza & Hartzell 1988). However, a detailed understanding of the spatio-temporal relation between aseismic slip and distribution of aftershocks remains elusive.
Here, we combine geodetic and seismological data of the M w 8.8 Maule 2010 earthquake and other recent great subduction earthquakes (2005Nias, 2011 to explore the temporal and spatial relationship between seismic and aseismic slip. First, we calculate maps of cumulative seismic slip (CSS) estimated from earthquake catalogues using scaling relations (e.g. Wells & Coppersmith 1994;Blaser et al. 2010). These relations allow us to estimate the rupture size on the basis of earthquake moment magnitudes (M w ) for individual events, which in turn allows us to calculate the average slip in the rupture area (Kanamori 1977). The resulting maps of CSS are then compared to afterslip models. To gain a more global perspective, we also apply the method to the 2011 Tohoku and 2005 Nias earthquakes. Secondly, for earthquakes where a spatial-temporal afterslip model is available we suggest a method to map the relation between cumulative number of aftershocks and total afterslip (TA), which can be regarded as a proxy for the seismic generation efficiency (e.g. a measure of the fraction of seismicity relative to TA). This in turn will allow us to infer frictional properties of the plate interface and better study the transition from seismogenic to aseismic behaviour.

T H E 0 1 0 F E B RUA RY 7 M w 8 . 8 M AU L E E A RT H Q UA K E A N D I T S P O S T S E I S M I C S E Q U E N C E
The Chilean margin shows right-lateral oblique subduction of the Nazca Plate beneath South America. The present-day subduction rate in the Maule region of Chile is about 62.5-66.0 mm yr −1 with a bearing of N78 • -80 • E (Angermann et al. 1999;Kendrick et al. 2003), which implies an obliqueness of ∼20 • . The M w 8.8 Maule earthquake on 2010 February 27 ruptured the megathrust segment north of the great 1960 Chile earthquake (Plafker & Savage 1970) and filled a well known seismic gap (Kelleher 1972) which last broke in 1835 (Caldcleugh 1836). Before the 2010 Maule earthquake, high GPS velocities indicated that the plate interface was mostly not slipping (highly locked) (Ruegg et al. 2002). Therefore, Blue triangles indicate seismic station from the International Maule Aftershock Deployment (IMAD) and red triangles indicate GPS stations which were operating during the mainshock. Co-seismic slip distribution of the 2010 and 1960 Chile earthquake (5 m slip contours, Moreno et al. 2009Moreno et al. , 2012 are shown with blue and grey contours, respectively. Rupture zones of other previous earthquakes are shown with black lines: 1835 event from Barrientos (1988); 1985 earthquake from Beck et al. (1998). Yellow star indicates the hypocentral location of the Maule 2010 earthquake (Servicio Sismológico Nacional de Chile). Double-couple focal mechanisms from the gCMT catalogue (lower hemisphere projection, M ≥ 6.9) indicate the 2010 mainshock and large aftershocks. this area had accumulated a high slip deficit, which apparently was totally released by the Maule earthquake (Melnick et al. 2012).
The aftershock sequence was dominated by shallow thrust events on or near the plate interface with a band of seismicity downdip of the coseismic peak slip (Lange et al. 2012;Rietbrock et al. 2012). A concentration of thrust-faulting aftershocks was also observed between the main patches of co-seismic slip (Agurto et al. 2012). In the northern part of the rupture zone, Cumulative number of events for the whole series (red) and for the events in the northern (blue) and southern (green) rupture area. Continuous lines show fits for the modified Omori law for the first 300 d (with parameters). Strong earthquakes (M ≥ 6.5) are shown with black squares. Lower panel shows the residuals of the modified Omori law fit to the first 300 d of the aftershock sequence.
intense seismicity also occurred in the overriding plate in the Pichilemu sequence, which started on 2010 March 11 with a M w = 6.9 and M w = 7.0 normal-faulting aftershock doublet (Farías et al. 2011;Ryder et al. 2012;Lieser et al. 2014). Notably, the largest aftershock (M w = 7.4) occurred in the outer rise adjacent to the southernmost limit of the rupture approximately 90 min after the mainshock. The largest aftershock on the megathrust (M w = 7.1) for the main time period under consideration (2010 March 2 to 2011 April 15) occurred on 2011 January 2nd in the southern part of the rupture area ( Fig. 1).

Earthquake data
For seismicity we use the NEIC catalogue between 2010 February 27 and 2011 June 16 and magnitudes larger than 4.5. This value is above the magnitude of completeness (M c ) of 4.3 to 4.4 reported for the NEIC catalogue (Sipkin et al. 2000;Legrand et al. 2012) . The largest aftershock (M w = 7.1) for the time period under consideration occurred on 2011 January 2nd in the southern part of the rupture area ( Fig. 1). Additionally, we use the local seismicity catalogue of Lange et al. (2012), which is derived from a dense local network of seismic landstations. This catalogue contains 20 205 local aftershocks between 2010 March 15 and September 30 of which 2177 events are above the magnitude of completeness (M c = 3.5). 151 focal mechanisms were taken from the gCMT catalogue and are used to select thrust events which are related to the plate interface. Out of those we selected 111 thrust type aftershocks (2010 March 2 until 2011 May 12) with differences between the inclination of the SLAB 1.0 model (Hayes & Wald 2009) and dip less than 20 • and differences between the slip vector to the trench normal direction of less than 30 • .
The temporal behaviour of the aftershock seismicity closely follows the modified Omori Law (eq. 1) for the 300 d after the Maule 2010 mainshock (Fig. 2). The data were fitted to the modified Omori law (eq. 1) by minimizing the negative log-likelihood function using the procedure from Ogata (1983). The modified Omori law describes the first 300 d of the aftershock series as can be seen by the small residuals during this initial period (Fig. 2, lower panel). After postseismic day 300 the modified Omori law does not give a at UB Kiel on November 20, 2014 http://gji.oxfordjournals.org/ Downloaded from good fit because of the occurrence of major earthquakes. The exponent p of the modified Omori Law is 1.02 for the seismicity south of 36 • S and 1.192 in the northern part of the rupture, and C S is less than 1 d.

GPS observations and afterslip model
A dense network of continuous GPS (cGPS) stations was installed in the Maule region after the main shock (Bevis et al. 2010;Vigny et al. 2011) densifying the existing cGPS network (Bevis et al. 2001;Ruegg et al. 2009;Vigny et al. 2009). The large number of GPS stations allows an unprecedented spatio-temporal resolution of the rapidly time-varying postseismic response in near, medium and far-fields following a large subduction earthquake. In this study, we use GPS data from 16 stations that were operating before and after the Maule earthquake, and from 42 stations that were installed days to weeks after the event (Figs 1 and S2). Additionally, we included all regional IGS stations south of the equator for the processing. Then, we estimated daily station coordinates and their uncertainties using the Bernese GPS Software (Dach et al. 2007) between 2010 February 27 until 2011 December 31. The postseismic phase was isolated by starting the day after mainshock on 2010 February 28 (e.g. the postseismic displacement is zero at the day after the mainshock). More information about the estimation of the daily geodetic displacements can be found in Bedford et al. (2013) and in the supplementary material.
The main feature of the cGPS time-series is the postseismic signal after the Maule earthquake, evident by an ongoing trenchward (westward) movement (Fig. 1) and a temporal decay. We note that the displacements are not proportional to ln (t) (Fig. S1) and that high values of C G of the order of several days (Fig. S3) are needed to fit the displacements with eq. (3). Alternatively, we tried to fit a power law function (d = At B ) to the displacements but eq. (3) (with large values of C G ) results in an overall better fit because of smaller residuals. The cumulative postseismic displacements of the Maule earthquake indicate significant lateral variations on the magnitude of the deformation (Fig. S1). Maximum westward displacement is found at MOCH, the station on Mocha Island, which moved ∼73 cm in the year after the earthquake (Fig. S4).
We use the spatio-temporal afterslip model from Bedford et al. (2013) for the time between 2010 March 2 and 2011 April 15 based on the Principal Components Analysis Inversion Method (PCAIM, Kositsky & Avouac 2010). The model uses a geophysically constrained plate interface geometry (Tassara & Echaurren 2012) where the interface is made of triangular patches with an average patch area of 190 km 2 and extending from the trench to a maximum depth of 100 km between 35 • and 39 • S. One major assumption of the model is that all slip is on the plate interface. The model assumes homogeneous elastic parameters and a Poisson's ratio of 0.25 and a shear modulus of 35 GPa. The model resolution is best around the coast and the Coastal Cordillera and decreases towards the trench. In order to minimize the influence of the major crustal normal-faulting event in the region of Pichilemu, the displacement occurred between March 10 and 12 was subtracted from the cumulative afterslip for all later times.

M A P P I N G A F T E R S L I P
First, we want to establish two ways to compare the seismic and aseismic slip following a larger earthquake. In the interest of brevity and clarity we define the following terms: (i) seismic afterslip (SA) represents the cumulative slip caused by all aftershocks and can be estimated from seismological observations. We use the term afterslip in SA since it is seismic slip that occurs during the postseismic period. (ii) pure aseismic afterslip (AA) is that part of the postseismic slip on the plate interface not caused by aftershocks. It cannot be inferred directly but can be deduced from TA and SA (see below). (iii) Total afterslip (TA) is the sum of SA and AA. TA can be estimated by inverting GPS displacement observations for slip on a predefined fault plane and therefore depends on assumptions such as the plate geometry and shear modulus.
The main relation between these quantities is that TA implicitly contains the displacements from the aftershocks (SA) and the AA, that is TA=SA+AA and therefore SA ≤ TA. Furthermore, when considering the total equivalent moment (e.g. averaging over the whole fault surface) most of the studies observe SA AA (e.g. Hsu et al. 2006;Ozawa et al. 2012). In order to relate SA with an afterslip model (TA) we first need to add the contributions of all aftershocks.

Representation of seismicity using the cumulative seismic slip (CSS)
SA for individual events can be calculated using the scaling relations which relate the seismic moment with rupture lengths and areas (Fig. S5). We calculate the rupture areas (A) of aftershocks using the scaling relations of Blaser et al. (2010) for subduction zone environment and thrust events based on an event range from M w = 6.1 to 9.5. These scaling relations are similar to the scaling relations of Wells & Coppersmith (1994). In order to estimate the slip for the individual aftershocks from their moment magnitude we use the equation from Kanamori (1977) to relate moment magnitude (M w ) to seismic moment (M 0 ). The slip s of the fault can then be calculated with We assume a constant shear modulus μ = ρ * v 2 s = 35 GPa calculated from v s = 3385 m s −1 and ρ = 3050 km m −3 derived from tomography and gravity studies (Tassara et al. 2006;Haberland et al. 2009;Hicks et al. 2012). Together with the scaling relations this results in a unique mapping between magnitude, slip and area ( Fig. S5). For simplicity we assume circular rupture areas with uniform slip. While the abrupt transition at the edge of the uniform slip region is not physically realistic, any tapered rupture model or more complex rupture geometries would introduce additional unconstrained parameters into the calculation without substantially effecting results. The largest event for which slip is mapped this way is the M w 7.1 thrust event in the southern part of the rupture area ( Fig. 1). At the scales examined in this study, the assumption of a roughly circular region should be satisfied for an event of that magnitude. We calculate CSS by adding the slip of all aftershocks. The slip of overlapping aftershocks will add up in case of overlapping source regions (Fig. S6). CSS has some similarity with the graphical representation of slip by Nadeau & McEvilly (1999) where slip of repeating event clusters is added within circles of constant diameter. Here, we additionally consider the size of the events from scaling relations. Fig. 3, panel A shows the CSS together with the TA from Bedford et al. (2013) for the Maule 2010 event in the same colour scale. For the Maule 2010 event this representation shows clearly that the aftershocks only cover a small portion of the fault zone, implying that large parts of the subduction interface show no seismic energy release during the aftershock sequence. This finding supports the hypothesis of aftershocks being driven by afterslip and occurring on seismogenic patches within a largely aseismic background (Perfettini & Avouac 2004b). However, the scaling relations are based on an empirical global compilation of earthquakes and the underlying assumption is that the scaling relations are valid for smaller magnitudes as well. Nevertheless, in order to fill the whole plate interface the scaling relations would need to differ significantly which we think is unlikely given the self-similar behaviour of seismicity down to very small magnitudes (e.g. Kwiatek et al. 2010).

Method 1: relating CSS to TA and the different spatial resolution of both observations
In the next step we calculate the percentage of slip from the CSS in relation to the TA model from Bedford et al. (2013). The result is shown in Fig. 3, S5) imply that single aftershocks with magnitudes larger than M w ≥ 5 exceed ∼20 cm slip which is in most places larger than at UB Kiel on November 20, 2014 http://gji.oxfordjournals.org/ Downloaded from the TA estimated by inverting GPS data. In the northern part of the 2010 Maule rupture a significant number of aftershocks occurred in the region of the southwest dipping focal planes of the Mw 6.9 normal faulting event on 2010 March 11. These events occurred in the crust of the overriding plate, so the SA model in Fig. 3 in this region will not represent faithfully the SA on the plate interface. However, south of 35.5 • S the slip of still a significant number of events exceeds the afterslip from the GPS in a region where largely aftershocks related to the plate interface were located (Lange et al. 2012). Agurto et al. (2012) compared the seismic slip of the first 12 d with the TA model from Vigny et al. (2011) and similarly found regions where the slip from seismicity exceeds the TA model. As the TA on the plate interface must be inferred from observations of surface displacements the TA model is severely limited in spatial resolution (Hetland & Simons 2010), particularly offshore, and effectively averages small slip patches like those caused by moderate and even strong events over larger regions. Therefore, CSS is much more heterogeneously distributed compared to TA models which is reflected by the spatially heterogeneous CSS map compared to the TA model.
In order to better quantify the smearing effect that spreads out SA in the TA inversion we forward calculated displacement data from slip derived from scaling relations (for the Maule 2010 case). Then the displacement data were inverted back to an afterslip model in order to directly compare between slip predicted by the scaling relation and by the geodetic inversion. The regularization, station distribution and model constraints for this inversion are the same as used for the inversion of TA (more details are found in Bedford et al. 2013). The advantage is that the SA is subject to the same model resolution and smearing as TA. Therefore, we believe the comparison of TA and SA for Maule 2010 case is as fair as possible. We made a number of tests and found that the seismic moment is recovered for different smoothing values. We found no significant lateral offset during these recovery tests. The input slip model and the inverted model show the tendency of the inversion to distribute the inverted slip into a larger region while reducing the amplitude of the signal in order to produce the same moment (Fig. 4). The pattern of the inverted slip model consists effectively only of two larger anomalies and small slipping areas are not reproduced, particularly close to the trench. Because of this smearing out of strongly peaked seismic slip, the influence on the cumulative TA model is so small in most areas, that the TA model can be viewed as a (somewhat degraded) proxy for the unobservable AA model. By comparing the TA model with the CSS model we can thus directly make inferences about the relative local importance of SA and AA. The implicit assumption is that the true AA is much smoother than the slip caused by earthquakes. This assumption is physically reasonable and would be expected from rate-state behaviour in the velocity-strengthening regime but cannot easily be tested directly.

Application to 2011 Tohoku and 2005 Nias earthquakes
We carry out an analogue analysis using published afterslip models for the M w = 9.0 2011 Tohoku (Ozawa et al. 2011) and the M w = 8.7 2005 Nias (Hsu et al. 2006) earthquake. During the 14 d of afterslip of the 2011 Tohoku earthquake only five aftershocks exceed 2.5 times the TA model (Fig. 5) which can also be explained with the different spatial resolutions of CSS and TA. A much larger difference between the afterslip estimated from seismicity and the TA model (covering 11 months of afterslip) is found for the Nias earthquake (Fig. 6). For the Nias earthquake, Tilmann et al. (2010) attempted a first estimate of the spatial variability of the relationship between TA and SA and found that the ratio between cumulative seismic and equivalent afterslip total moment differs by a factor of ∼20 between four segments of the Nias rupture. Here, we are examining this relationship at higher spatial resolution and find that the largest aftershocks of the Nias 2005 event occur in regions of very low afterslip (<0.15 cm) so that the cumulative slip from a number of stronger aftershocks exceeds the afterslip model (Hsu et al. 2006) more than 12-fold (events indicated in black Fig. 6,  panel B).
The finding that the total moment release of the aftershocks locally exceeds by far the afterslip (as in the Nias 2005 case) is difficult to explain just by the spatial smearing of small local afterslip patches into larger ones of the afterslip model. As by definition, the true TA must be larger than or equal to the true SA, we need to discuss the extent to which the calculation of SA can be erroneous.
The scatter of the data set on which the scaling relations are based indicates that rupture areas for earthquakes of the same magnitude can vary by up to a factor of 5 in the moment magnitude range between 6 and 7 (Blaser et al. 2010) and the estimated slip values will vary accordingly. Since the scaling relations are the result of a regression analysis, which results in a curve roughly in the middle of this range, this would only account for an exceedance of ∼2.2 times the afterslip if we assume that the aftershocks of the Nias 2005 events are exclusively on the upper bound of the magnitude-area distribution. One possibility to explain the discrepancy might be related to strong local heterogeneity due to the presence of material with strongly increased shear modulus, which would result in a smaller slip estimate for the individual aftershocks. Such small regions of elevated shear modulus might be a possible explanation for the 2010 Maule and 2011 Tohoku earthquake.
For the Nias 2005 earthquake, however, this explanation would result in unrealistically high values of the shear modulus. Even assuming rocks with extreme high shear modulus such as Peridotite or Dunite (shear modulus of 64 GPa) would influence the slip only by a factor of two because, for a given scalar seismic moment, the shear modulus is related linearly to slip. Another potential explanation might be due to a systematic shift of the hypocentres from NEIC for the Nias region. For the gCMT catalogue Tilmann et al.  Fig. 6) the moment release in the postseismic phase is mostly characterized by aftershocks and only to a small degree, if at all, by AA, in contrast to the large scale pattern where AA dominates. To further test this assertion we calculated averaged moment release from CSS and the 11-month afterslip model (Hsu et al. 2006) for circles of variable radius and centred on the 2005 Nias mainshock epicentre; aftershocks up to M w = 6.9 occurred in this region, which is characterized by reduced slip in the afterslip model. Fig. 7 suggests that for circles smaller than 50 km radius the CSS exceeds the afterslip model. CSS needs to be smaller than TA by definition, and we interpret this discrepancy as originating from the limited resolution of the TA model as discussed above. In spite of uncertainties, we conclude that in the region of the Nias 2005 hypocentre a number of stronger aftershocks locally exceed  the AA and that these aftershocks therefore are not solely driven by afterslip. It is highly likely that postseismic deformation is locally dominated by SA even though on the scale of the whole rupture AA is much more significant. This suggests that the ratio of SA and AA varies significantly when looking at spatial scales of 50-100 km.

Method 2: calculation of afterslip to aftershocks ratio (ASAR)
Based on an observed linear relation between cumulative number of aftershocks and displacements and on rate-state concepts, Perfettini & Avouac (2004a)  triggered by the stresses imposed by afterslip with the rate of triggering being proportional to the rate of stressing. In order to test whether this linear relation holds also true locally for the Maule postseismic sequence we plot the cumulative number of aftershocks with a magnitude larger than a threshold M c within a radius of 50 km around each GPS station against postseismic displacement (Fig. 8). The search radius of 50 km was chosen by trial and error and represents a reasonable compromise between available GPS data and the need to ensure a sufficient number of events in the neighbourhood of most stations for the NEIC and the local catalogue (Fig. S7). A smaller radius results in a smaller number of events within the search radius. Before the trial and error procedure we have excluded GPS sites outside the Maule 2010 rupture zone where only sparse seismicity occurred. This was done by only considering GPS sites where seismicity occurred on more than 20 d within a 50 km search radius in the period covered by each catalogue.
As discussed in Section 1, the integrated version of the modified Omori law for p = 1 (eq. 2) and the formula for GPS displacements for a steady-state sliding system (eq. 3) have similar time-dependent functions. If the time constants governing the decay are similar, then it was concluded that the temporal evolution of aftershocks is governed by the time evolution of the stress induced by afterslip (Perfettini & Avouac 2004a,b), and that the delay introduced by the nucleation time of the aftershocks has a negligible effect. In this case C G and C S should have similar values as observed for the Nias 2005 earthquake (supplementary material of Hsu et al. 2006). However, for the Maule 2010 case C S is on the order of half a day and C G of several days indicating that aftershock and afterslip are related in a more complicated way. The first-order relation between both observations is still linear (Fig. 8). We note that station MOCH located on Isla Mocha is characterized by a small number of aftershocks in relation to its displacement. This station is located closer to the trench compared to other stations ( Fig. 1  and S8). The station was not operational from 2010 July 16 until 2011 March 3. During this time the station was displaced 36.8 cm to the east, presumably due to the strong aftershock (M w = 7.1) on 2011 January 2 (Fig. 1). As a result the linear relation between cumulative number of events and displacements and equally the temporal behaviour of the displacements of station MOCH cannot be approximated by simple time-dependent functions until 2011 March 8 (136 d after the 2010 mainshock, Fig. S1).
While for the calculation of seismicity only limited assumptions are made (e.g. velocity model), the method assumes implicitly that the major part of afterslip and the seismicity occurs on or close to  (Hsu et al. 2006) and aftershocks shown for the same 11 months period after the mainshock. Grey lines indicate co-seismic slip distribution (2.5-m slip contours, Konca et al. 2007).

Figure 7.
Cumulative moment release (shown as equivalent moment magnitude, M w ) determined by spatial averaging of CSS and afterslip model (Hsu et al. 2006) for circles of different radii around the Nias 2005 mainshock hypocentre (yellow star, Fig. 6B). the subduction interface. In fact, this is true only for the southern part of the Maule 2010 rupture area where most of the aftershocks are related to the plate interface, whereas the northern part of the Maule rupture area also hosts the aftershock seismicity of the normal faulting sequence of Pichilemu in the overriding plate. Although we still observe a linear relation between aftershocks and displacements in this area (north of 35.5 • S), a significant portion of the observed displacement is presumably influenced by deformation in the overriding plate.
In the next step, we calculated the cumulative number of events in relation to the Bedford et al.'s (2013) temporal afterslip model. This is possible because the temporal decays of the individual nodes of the inverted afterslip model follow similar time-dependent decays (compare Fig. S1 and S9). Instead of measured surface displacements from GPS we use the average displacement in the afterslip model within 50 km search radius around each gridpoint on the megathrust. Because the megathrust is not horizontal these cells will appear to be ellipses in map view. For the seismicity we consider only events from the same ellipse. We consider the dip of the slab by calculating ellipses using the dip of the megathrust from the SLAB1.0 model beneath the cells of the afterslip model. The cartoon in Fig. 9 shows how the values of afterslip to aftershock ratio (ASAR) were calculated and the result for the Maule 2010 earthquake is shown in Fig. 10, panel A. Due to the seismicity in the overriding plate mentioned above we repeated the analysis but excluded the events from the local seismicity catalogue with vertical distances of more than 10 km to the SLAB1.0 model (Fig. 10,  panel B).
Following the argument in the introduction ASAR can be seen as a proxy of the variability of the ratio of equivalent moment release of the total displacement estimated from GPS to the moment released by seismicity. We prefer the cumulative number of events rather than the seismic moment as the moment sums are more strongly influenced by the largest earthquake (Frohlich 2007 (Lange et al. 2012) occurring within 50-km radius around the GPS stations and for events with magnitudes larger than 3.5. Note that the local seismicity catalogue begins 21 d after the mainshock and ends 2010 September 30: the cumulative displacement was set to zero on the first day that an event occurred within the search radius of the station in order to make it easy to compare the slopes of the different lines.
deeper part of the afterslip model slip along the whole rupture due to the small number of events at more than 50 km depth. For the shallow part (offshore) the results shows lateral changes alternating between moderate and small values (compared to the down-dip end of the megathrust) which is suggested to reflect spatial changes of the properties of the subduction interface. The linear regression fits show larger errors of the estimated slope parameter for regions near the outer limits of the modelled interface (Fig. 10, panel C). We explain these with the small afterslip near the boundaries of the afterslip model together with a smaller seismicity rate which results in poorly constrained slope parameters of the regression. The larger deviation from linearity outside the main rupture (Fig. 10, panel C) is due to the smaller number of aftershocks outside the main rupture and reduced resolution of the temporal afterslip model rather than due to a breakdown of the linear relation between aftershocks and the cumulative number of events.
The absolute values of the linear regressions are dependent on the smallest magnitudes incorporated in the calculation. We repeated the same calculation by only considering events with magnitudes larger than 4 (Fig. S10) which results in a smaller number of aftershocks taken into account. As a result the absolute values increase by an average of 31 per cent with a standard deviation of 4 per cent. Clearly, the values of ASAR are dependent on M c of the seismicity catalogue, so only relative differences of ASAR should be interpreted. Additionally, strong spatial changes of the b-value of the Gutenberg-Richter distribution will distort the result and have to be considered. For the Maule earthquake this is in particular the case in the region of the Pichilemu aftershock cluster which has a large number of small events (Lange et al. 2012). The b-values of aftershocks of the 2010 Maule rupture vary moderately between 0.9 and 1.2 using the NEIC catalogue and M c = 4.5 . We tested the robustness of the results for different M c and find that the crustal seismicity of the Pichilemu sequence (large number of small magnitude events) influences the outcome of the method for the northern part of the rupture, while the relative values of ASAR south of 35.5 • S still indicate the same anomaly pattern when considering the relative changes of neighbouring grid nodes. Hence, these variations of the b-value will not significantly change the results south of 35.5 • S.
The values of ASAR reflect the ratio of aftershocks to geodetic afterslip and help in mapping the average frictional properties of the megathrust in the postseismic period. We expect that the values of ASAR are sensitive to changes of material (e.g. Kopp 2013), temperature (Oleskevich et al. 1999) or overpressured fluids (Husen & Kissling 2001). We observe a heterogeneous distribution of ASAR both in down-dip and lateral directions. Down-dip of the co-seismic rupture of the Maule 2010 event elevated values of ASAR are found because of the reduced seismicity at larger depths (>50 km) with still observable afterslip (Fig. 10). Because of the crustal seismicity in the region north of 35.5 • and low resolution of focal depths in the offshore regions the values of ASAR differ in these regions for a modified data set where events with depths more than 10 km from the plate interface have been removed (compare Fig. 10 A and Fig. 10 B). We restrict the interpretation of the values of ASAR to the regions of robust and reliable ASAR defined by the diagonal of the resolution matrix from Bedford's et al. afterslip model (shown with a green line in Fig. 10) and for regions south of 36 • S in order to minimize the influence of crustal seismicity. Along the coastline at 37 • S we find reduced values of ASAR due to reduced afterslip in the Bedford et al.'s (2013) model. Additionally, the region is characterized by stronger aftershocks (Fig. 3) and by its location between the two co-seismic slip maxima of the Maule 2010 event (e.g. Vigny et al. 2011).
We compared this minimum as the most robust and best resolved feature in the ASAR distribution with the distribution of co-seismic slip, Bouguer anomaly  and interseismic coupling . We cannot find a simple relation between these measures (Fig. S12) which suggests that the distribution of ASAR for the Maule 2010 earthquake is not strongly controlled by these parameters. There might be a weak correlation between coseismic slip and ASAR in that the minimum in coseismic slip is associated with a local maxima of ASAR, both in the region around 36 • S, and at the edges of the rupture. At larger depths of around 50 km we find elevated values of ASAR reflecting increasingly aseismic behaviour of the plate interface.

C O N C L U S I O N S
We focus on the relation of afterslip models derived from GPS displacements (implicitly including aseismic and SA) and aftershock seismicity for the Nias 2005, Maule 2010 and Tohoku 2011 subduction-zone earthquakes: (i) We calculate the cumulative contribution of the aftershock seismicity (considering slip and areas of events), estimated from scaling relations by summing up the contributions of the individual events from a seismicity catalogue. The resulting spatial grid represents the cumulative seismic slip (CSS) of the aftershock sequence. By comparing this cumulative SA with afterslip inverted from GPS data we find that aftershocks only cover a small part of the rupture area as expected from the much smaller moment release of the at UB Kiel on November 20, 2014 http://gji.oxfordjournals.org/ Downloaded from Figure 9. Cartoon showing the counting of aftershocks within a radius of the megathrust assigned to the centre of the cell. The afterslip value is calculated by averaging the afterslip over the same circle. Since the megathrust surface is not horizontal, these cells will appear to be ellipses in map view.
aftershock seismicity compared to the total moment released by the afterslip seen from GPS. However, the slip estimated from scaling relations from aftershock seismicity exceeds locally the afterslip model estimated from GPS. This must be a consequence of both the smearing inherent in geodetic inversions and the uncertainty in converting moment to slip from both the variability of scaling relations and the value of the shear modulus. For the Tohoku 2011 and the Maule 2010 earthquake the overshoot is moderate (up to factor 2) but still implies that locally the displacement due to the largest aftershocks 'runs ahead' of the AA. This effect is much more extreme for some strong aftershocks of the Nias 2005 earthquake where we observe regionally (around 100 km diameter) a very high fraction of seismic moment release equivalent to a moment magnitude of M w = 7.1. Taking into account the uncertainties of the cumulative slip from seismicity and the afterslip model this moment release is comparable with the moment release from the afterslip model indicating a very small amount of aseismic slip in this area. Although such an exceedance is only very localized such a behaviour is incompatible with aftershocks driven purely by afterslip and suggests that the ratio of SA and AA varies significantly when looking at the spatial scale of aftershocks.
(ii) We investigate the first order linear relation of cumulative number of events for smaller regions around individual geodetic stations for the Maule M w 8.8 earthquake in Central Chile on 2010 February Maule 27. The first order temporal evolution of aftershocks is governed by the time evolution of the stress induced by afterslip but we find different values for the parameters involved in the modified Omori law and temporal behaviour of the postseismic displacements (assuming a similar functional form for postseismic displacements C S and C G in eqs 2 and 3, respectively) suggesting that for the Maule 2010 earthquake the postseismic displacements and the cumulative number of aftershocks are not completely linearly related, in particular when considering several hundred days of afterslip. Based on the first-order linear relation between the afterslip model and the cumulative number of aftershocks for regionalized subsets of events we suggest a method to map the ratio between afterslip and aftershocks (ASAR) which we use to map the relation between both processes along the Maule 2010 rupture. The values of ASAR show significant lateral changes and reduced values at the coast at 37 • S. For the downdip end of the seismogenic zone (∼50 km depth) elevated values of ASAR suggest an increase of aseismic faulting with increasing depths.

A C K N O W L E D G E M E N T S
We thank the German Research Foundation (DFG) for the funding of the MARISCOS bundle project, grant LA 2970/1-1 and MO2310/ 1-1. We thank all persons and organizations who deployed, serviced and provided seismic stations for the seismic network along the Maule 2010 rupture zone International Maule Aftershock Deployment (IMAD). We acknowledge the U.S. National Science Foundation for funding an emergency deployment of ∼30 new cGPS at UB Kiel on November 20, 2014 http://gji.oxfordjournals.org/ Downloaded from stations in the area of the M w 8.8 Maule event via a RAPID award, and for a second grant from its Continental Dynamics program which has considerably extended the lifetime of this geodetic network. We thank UNAVCO, Inc. for its help during the construction of this cGPS network in 2010, especially with regards to data telemetry. The processing of the GPS data was supported by FONDECYT (grant 1101034). We thank the Servicio Sismológico Nacional de Chile (SSN) and the University of Santiago de Chile for logistical support during the deployment. We thank three anonymous reviewers, Sylvain Barbot and the editor Duncan Agnew for their thoughtful reviews of this manuscript. Figures were generated using GMT (Wessel & Smith 1998).

S U P P O RT I N G I N F O R M AT I O N
Additional Supporting Information may be found in the online version of this article: Figure S1. Total postseismic GPS displacements plotted versus time after the mainshock in linear-linear scale (panel A) and semi-log scale (panel B). Displacements are plotted such that the displacement is zero at the day after the mainshock. Only stations where data is available for the first day after the mainshock are shown. Displacement errors are shown as vertical lines (partly obscured by black points). Station MOCH located on Isla Mocha was not operational from 16 July 2010 until 3 March 2011. During this time the station was displaced 36.8 cm to the east, presumably due to the strong aftershock (M w = 7.1) on 2 January 2011. More detailed timeseries of individual GPS displacements are shown in the supplementary material of Bedford et al. (2013). Figure S2. Map showing the distribution of seismic stations from the IMAD deployment (inverted green triangles) and the continuous GPS stations with circles (labeled with station name). Coseismic slip distribution of the 2010 Maule earthquake shown with blue lines (5 m slip contours, Moreno et al. 2012). Plate convergence after Angermann et al. (1999). Figure S3. Residuals of total displacements after fitting equation 3 for stations which produced a reliable fit. 9 stations clearly outside the co-seismic rupture (NIEB, VALN, LNDS, VAPL, PORT, ROBL, SANT VNEV, LLFN) with small displacements are not shown. Additionally, stations with many gaps or influenced by strong aftershocks (MOCH) are not included. The resulting variables K G and C G of the individual fits are given in the list. Stations BAVE, DGF1 and CONZ (shown in blue) result in a poor fit and we fixed C G to 12.5 d. Figure S4. Aftershocks (NEIC catalogue) and displacements for the station MOCH which cannot be approximated with simple time dependent functions. Displacements (left) and cumulative number of events (centre) plotted versus day after main event. On the right the total displacement versus the number of aftershocks is shown. Figure S5. Scaling relations from (thick lines, Blaser et al. 2010, optimized for subduction zones) and (thin lines Wells and Coppersmith 1994, thrust events). The scaling relation is shown for the calculated slip (using a shear modulus of μ = 35GPa) and the total rupture area. Figure S6. Cartoon illustrating the calculation of CSS by spatially summing up individual seismic afterslip contribution (estimated from scaling relations) for each aftershock. Figure S7. GPS displacements after the Maule 2010 earthquake plotted versus the number of local events occurring within 50 km radius around the GPS stations. A: Displacement versus the local aftershock catalogue (Lange et al. 2012) for events with (M ≥ 3.5). Same as Fig. 8 but only using events with depth differences of less than 10 km to the SLAB1.0 reference model (Hayes & Wald 2009). Because the seismicity catalogue begins 21 days after the mainshock, we set for the first pair of earthquake and displacements the values to zero, such that all lines start at the origin of the plot making a direct comparison of the slope of lines more feasible. Therefore, only the slope of the lines are meaningful. B: Like Fig. 8 but the y-axis is reduced with a reduction value of 1 mm/event (y = y−(x * 1 mm/event)). Additionally this panel does not show stations north of 35.5 • S in order to exclude an influence of events from the crustal normal faulting sequence of Pichilemu. C: Like a but using events from the NEIC catalogue (M ≥ 4.5, 28 February 2010 until 31 December 2011).We set for the first pair of earthquake and available displacement data the values to zero. Figure S8. Slope of lines shown in Fig. 8A using the NEIC catalogue (panel A, 2 March 2010 until 14 June 2011) and the local aftershock catalogue (panel B, Lange et al. 2012) in mapview. Linear regression for different time windows and catalogues leads to similar relative values. The main difference between both maps is larger number of data points for the deep part of the rupture zone using the local catalogue due to resolution of the downdip events more inland which is not resolved by the NEIC catalogue. In the region of the crustal aftershock cluster of Pichilemu (stations PTPC, LEMU and PMQE) the local catalogue resolves a significant large number of events which are not resolved by NEIC resulting in smaller values. The yellow star indicates the main shock epicentre from the SSN catalogue. Co-seismic slip distribution is shown with blue lines . Volcanoes are plotted as upright grey triangles (Siebert & Simkin 2002). Figure S9. GPS displacements from afterslip model (Bedford et al. 2013), inverted onto the plate interface, versus time after the mainshock in log scale. Each line represents the temporal behaviour of an individual node of the afterslip model. The time interval of the Pichilemu event which was not included in the afterslip model is shown with the blue lines. Figure S10. Same as Fig. 10 but only considering events with magnitudes larger than M ≥ 4 to show the moderate influence of b-value changes in the central part of the model. A: Result using the local seismicity catalogue with M ≥ 4. Blue lines show the cumulative afterslip model (0.5 m contour lines) for the same time period. The green line encloses the region with good resolution of the afterslip model from Bedford et al. (2013). Only nodes with more than 20 pairs of aftershocks and total slip inverted on the plate interface are shown. Black line indicates the location of the profile shown in Fig. S12. B: Plot showing the difference (in %) using M c ≤ 3.5 (shown in Fig. 10A) and M c = 4 (panel A). We observe a sensitivity of the method which is related to strong spatial changes in the b-value of the Gutenberg-Richter distribution. In particular the region close to the city of Pichilemu (∼34 • S) is characterized by an pronounced aftershock series of swarm-like behaviour which can be seen by differences of the linear regressions using different values of M c . The inlay shows the histogram (mean: 39%, standard deviation: 9%) of the difference of both calculations. Figure S11. Forward calculated displacements of normal faulting events estimated by the formula from Okada (1992). A: Seismicity from the gCMT catalogue (28 February 2010 until 15 April 2011) colored according to the type of thrusting (red: normal faulting events, green: thrust type events, blue: strike-slip events). Magnitudes for the normal faulting events used for the forward calculation of slip in B and C are labeled with magnitudes. B: Calculated East-West displacements for the normal faulting events shown in A. C: Vertical displacements.  Moreno et al. (2010). Location of the profile is shown in Fig. 10A. Blue lines show ASAR values using the whole seismicity catalogue (Fig. 10, panel A), red lines show ASAR values of events which distances of less than 10 km to the SLAB1.0 model (Fig. 10, panel B). Region with crustal seismicity is indicated with black arrows. Bouguer anomaly is plotted for regions onshore and the free-air anomaly offshore . (http://gji.oxfordjournals.org/lookup/suppl/doi:10.1093/gji/ ggu292/-/DC1).
Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.