Mechanical origin of aftershocks

Aftershocks are the most striking evidence of earthquake interactions and the physical mechanisms at the origin of their occurrence are still intensively debated. Novel insights stem from recent results on the influence of the faulting style on the aftershock organisation in magnitude and time. Our study shows that the size of the aftershock zone depends on the fault geometry. We find that positive correlations among parameters controlling aftershock occurrence in time, energy and space are a stable feature of seismicity independently of magnitude range and geographic areas. We explain the ensemble of experimental findings by means of a description of the Earth Crust as an heterogeneous elastic medium coupled with a Maxwell viscoelastic asthenosphere. Our results show that heterogeneous stress distribution in an elastic layer combined with a coupling to a viscous flow are sufficient ingredients to describe the physics of aftershock triggering.

Aftershocks are the most striking evidence of earthquake interactions and the physical mechanisms at the origin of their occurrence are still intensively debated. Novel insights stem from recent results on the influence of the faulting style on the aftershock organisation in magnitude and time. Our study shows that the size of the aftershock zone depends on the fault geometry. We find that positive correlations among parameters controlling aftershock occurrence in time, energy and space are a stable feature of seismicity independently of magnitude range and geographic areas. We explain the ensemble of experimental findings by means of a description of the Earth Crust as an heterogeneous elastic medium coupled with a Maxwell viscoelastic asthenosphere. Our results show that heterogeneous stress distribution in an elastic layer combined with a coupling to a viscous flow are sufficient ingredients to describe the physics of aftershock triggering.
The first empirical law for aftershock organisation in time dates back to Omori 1,2 and states that the number of aftershocks n(t) decays as a power law with the time t from the mainshock, ( ) ( + ) − n t t c p . Many explanations for the Omori law have been proposed 3 but a complete understanding of its origin is still lacking. The improvement in data acquisition and elaboration has contributed to identify 4 the dependence of the characteristic time c in the Omori law on the rake angle λ . This angle indicates the direction of slip on the fault plane and can be related to the local level of differential stress σ D 5 . In particular under some assumptions, such as that faulting follows Mohr-Coulomb theory, σ D is larger for λ ∈ , ]. Within these hypotheses [4][5][6] , the rake angle can be, therefore, used to infer information on the differential stress acting on seismic faults, a quantity very difficult to measure directly. A similar dependence on λ has been previously observed 6 for the parameter b in the Gutenberg-Richter (GR) law 7 , stating that the number of magnitude m earthquakes, ( ) N m , exponentially decreases with m, ( ) ∝ − N m 10 bm . These results offer new perspectives in earthquake-hazard analysis even if a precise physical interpretation of their origin is still lacking. Moreover, a clear identification of the mechanisms responsible for the b and c dependence on the differential stress might also contribute to a better understanding of the physics behind aftershock triggering, an issue still debated and controversial [8][9][10] . In this letter we show that also the size of the aftershock area depends on λ and we develop a coherent framework able to explain aftershock organisation in time, space and magnitude.

Results
We first consider experimental data from the Southern California region and restrict the study to intermediate mainshock magnitudes ∈ . , . m [2 5 4 5] M . We identify aftershocks as events occurring within 10 min and in a circle of radius 3.3 km centered in the mainshock epicenter (see Methods).
We then define L a as the average main-aftershock epicentral distance normalized by the typical size of the aftershock area 11 ( ) = . × . l m 0 01 10 M m 0 5 M km. This choice, as shown in the Supplementary Information, ensures that the evaluation of L a is not affected by variations of the b value. We evaluate L a for all main-aftershock couples and finally stack sequences according to the rake angle λ of the mainshock in overlapping intervals of amplitude δλ =  5 . Only λ intervals containing at least 5 main-aftershock couples are included in the study. We have verified that results are not significantly affected by the interval value. We also evaluate the c-value for mainshock magnitudes ∈ . , . Results for the average value of c and L a as function of λ are plotted in Fig. 1. We also plot the dependence of b on λ , following the same procedure of ref. 6 without discrimination between aftershocks and mainshocks. Parametric plots of c vs b and c vs L a indicate ( Fig. 2) a proportionality among these quantities. This behaviour is also recovered for larger mainshocks and other geographic areas ( > . , supporting positive correlations among b, c and L a as a stable feature of seismicity. We wish to stress that the observed behaviour is not a spurious effect related to aftershock incompleteness 12,13 (see Supplementary Information).
The dependence of L a on λ leads to a better understanding of previous results on the c and b values 4,6 . Indeed, assuming that the fault area is proportional to the aftershock area, the seismic moment is proportional to σ ∆ L a 3 , where Δ σ is the stress drop due to the mainshock. Therefore, for a given mainshock magnitude m M (or seismic moment) a smaller value of L a implies a larger Δ σ. Conversely, in regions where Δ σ is smaller, the same m M can be only recovered if the stored elastic energy is distributed over a wider area. In the latter case, it is more probable to find several unstable regions scattered in space. This leads to a larger fraction of small aftershocks (a larger b value) and a longer temporal delay for the relaxation of all instabilities (a larger c value). The above description, therefore, predicts positive correlations among b, c and the spatial extent L a of the aftershock area, experimentally observed. Since, under similar fault conditions, it is also reasonable to expect larger Δ σ in regions with larger σ D , the above argument also provides an explanation for the dependence of measured quantities on λ .
In the following we will show that the experimental statistical features of aftershocks are recovered in a model for a single seismic fault. The model can be extended to describe more realistic fault networks including secondary faults and different orientations with respect to the mainshock fault. Off-fault aftershocks can be expected in this case but their number would not be so relevant to affect the observed statistical results. The model implements three main ingredients. Ingredient 1 is the assumption that the fault plane is an elastic medium modelled as blocks interconnected by springs and subject to a stress with constant rate σ  ext caused by the tectonic drive. More precisely we consider a tilted square lattice of spacing a. As soon as the local stress σ ij exceeds the local static friction σ ij th , the block slips and stress is distributed to nearest neighbor blocks with the dynamic friction coefficient μ D set to zero. An earthquake is represented by the ensemble of subsequent slips and its magnitude can be obtained from the size S of the slipping region. Ingredient 2 is the introduction of a heterogeneous local friction [14][15][16] assuming that σ ij th follows a quenched Gaussian distribution with mean σ A and standard deviation drawn, after each slip, from a Gaussian distribution of zero mean and standard deviation δσ p = σ . 0 6 A . The precise value of σ A does not affect our results and the only free parameters are the standard deviations in the stress drops and in the friction levels. The hypothesis of randomness in stress relaxation reflects the existence of asperities in the fault plane leading to irregular local slips with the possibility that more blocks are simultaneously unstable. An earthquake starts at the most unstable site and involves neighboring blocks. Unstable blocks not involved in the event keep their local stress value that will be, eventually, relaxed at subsequent times. Ingredient 3 is the postseismic relaxation caused by the coupling between the elastic lithosphere of thickness H l with a Newton viscous asthenosphere of thickness H a (Fig. 3). We neglect vertical variations of the local strain and carry out a force balance for a given element in the lithosphere as in ref. 17. Approximating the viscous flow in the asthenosphere as a linear Couette flow we obtain the following equation for the local stress evolution 17,18 expressed in terms of the lithospheric Young modulus Y and the asthenospheric viscosity η. Model parameters and further details can be found in the Methods.
In Supplementary Fig. 5 we plot the temporal evolution of a typical synthetic catalog. Aftershock sequences with patterns very similar to experimental data are clearly visible and the GR law, with a realistic value b = 1.1, is recovered under the assumption of local stress conservation. In order to identify the mechanisms for aftershock production we monitor the response of the system to a shear stress per- . Starting from an initially stable configuration and applying this stress perturbation, the excess of stress is relaxed via a mainshock whose magnitude is tuned by Δ 0 and R 0 . In the simplest version of the model (only ingredient 1) all the external stress is relaxed by the mainshock and aftershocks are not produced. The introduction of spatial heterogeneities (ingredient 2) leads to blocks that, not involved in the energy redistribution process during the main event, are still unstable after the mainshock occurrence. These blocks relax their energy at subsequent times and, as a consequence, aftershocks are triggered. Their activity is substantially constant in time and abruptly stops after a time delay c depending on the spatial extent of the perturbed region (R 0 ). Viscoelastic relaxation (ingredient 3) leads to aftershock activity continuing after c. In this case the aftershock number decreases as a power law of time with an exponent close to p = 1.1 (Fig. 4a) in very good agreement with the Omori law of real seismic data. In Fig. 5 we investigate the influence on our results of the two free parameters: the standard deviations of the local friction distribution (δσ th ) and the standard deviation in the value of the local stress after the stress-drop (δσ p ). The variance δσ th can be related to the number of asperities as well as to their size distribution within the fault. In our study we consider mainshocks with magnitude m = 6.3 and plot the number of events with m > 3.5 as a function of time from the mainshock for different choices of δσ p (left panel) and δσ th (right panel). We observe that results are substantially unaffected by δσ p whereas different values of δσ th lead to different results. More precisely, we observe that for δσ σ = .
0 1 th A , aftershocks follow the Omori decay up to a given time when their number abruptly  decreases to zero and the power law regime is no longer observed. On the other hand, for large values of δσ th a constant, roughly stationary, ( background) seismic activity is superimposed to the aftershock decay rate. In this case the power law decay regime reduces to less than one decade before approaching the constant background rate with a quite stable exponent p. We have also explored the influence of different values of μ D > 0 which only affects the level of the background rate, becoming larger for larger μ D . Conversely, the aftershocks decay is not affected by μ D and the Omori parameters p and c are μ D independent.
We wish to stress that many spring-block models, based on ingredients 1 and 2, have been proposed in the literature even implementing more complex, time dependent or state dependent, friction laws [19][20][21][22][23] . Even if these models exhibit non-trivial temporal patterns, they are not able to reproduce aftershock occurrence in agreement with experimental data. This observation, together with our findings (Figs 4  and 5), of numerical aftershocks following the Omori law, confirm previous results concerning the central role of viscous coupling (ingredient 3) for aftershock triggering 19,[24][25][26] . Similar results can be also recovered by the Jagla model [27][28][29] where δσ ij is constant, the friction thresholds σ ij th are randomly updated after each slip and a different equation for stress relaxation is implemented.
To explore the role of the level of differential stress in the aftershock organisation, we analyse different values of ∆ σ , ∈ . , .
[0 38 0 59] A 0 , keeping the mainshock magnitude = . m 6 7 M fixed. As a consequence, the value of R 0 is changed accordingly, with larger Δ 0 corresponding to smaller R 0 . In Fig. 4a we plot the number of events ( − ) n t t 0 with > . m 3 5 as function of time from the mainshock. Each curve is obtained by averaging over 10 different initial random configurations. The Omori law is observed for each value of Δ 0 and the c value (indicated by vertical arrows) decreases for increasing stress levels (larger Δ 0 ). We also find that the b value in the GR law is a decreasing function of Δ 0 (Fig. 4b), leading to positive correlations between the parameters c and b. Fig. 4 also indicates that larger c values correspond to larger values of R 0 , predicting a positive correlation between c and the size of the aftershock area L a . The parametric plots (insets of Fig. 4) reproduce the same linear trends observed in experimental data and therefore, by assuming a given relationship between Δ 0 and λ , the experimental results in Fig. 1 are reproduced by the numerical model. The agreement between experimental and numerical results indicates that the heterogeneous stress distribution in an elastic layer combined with a viscous coupling are necessary and sufficient ingredients to describe aftershock occurrence.

Methods
Mainshocks and aftershocks identification. We apply a space-time window criterion to discriminate between mainshocks and aftershocks 8 : An event is identified as a mainshock if a larger earthquake does not occur in the previous y days and within a distance L. In addition, a larger earthquake must not occur in the selected area in the following y 2 days. We use typical values L = 100 km, y = 3 and y 2 = 0.5. Aftershocks are all events with magnitude larger than m a = 2.4 occurring in the subsequent time interval ∈ , t t t [ ] 1 2 and within a circle of radius R from the mainshock epicenter. The sets of parameters t 1 , t 2 and R are listed in Table 2 of the Supplementary Information. For each sequence the c value is obtained by means of a maximum likelihood maximization routine keeping p = 1.1 fixed for all sequences. We finally average c over all sequences belonging to a given λ interval. In the case of large mainshocks > .  The spring-block model. We represent the fault plane as an elastic medium made by blocks on a tilted square lattice of spacing a interconnected by springs. Blocks are under the action of a uniform tectonic drive σ  ext in the x-direction, and are coupled to a Maxwell viscoelastic layer (Fig. 3). In the hypothesis that the slip is much smaller than a and that the stress redistribution after the slip is instantaneous, the system evolution can be expressed only in terms of σ ij . The simulation proceeds as follows: We randomly assign a quenched threshold σ ij th and initial condition σ ( = ) t 0 ij , at each site. Local stresses are then updated according to Eq.(1) whose discretized form, including the tectonic drive, reads i j i j ij , obeying local stress conservation. If at least one of these blocks is unstable, a further stress relaxation occurs and the process is iterated. The redistribution of stress stops as soon as no further nearest neighbor block is unstable. The whole process is considered instantaneous and afterwards the temporal evolution is iterated according to Eq. . Stress is expressed in units of the average value σ A whose value is irrelevant.