under a Creative Commons License. Nonlinear Processes in Geophysics Loading rates in California inferred from aftershocks

We estimate the loading rate in southern California and the change in stress induced by a transient slip event across the San Andreas fault (SAF) system in central California, using a model of static fatigue. We analyze temporal properties of aftershocks in order to determine the time delay before the onset of the power law aftershock decay rate. In creep-slip and stick-slip zones, we show that the rate of change of this delay is related to seismic and aseismic deformation across the SAF system. Furthermore, we show that this rate of change is proportional to the deficit of slip rate along the SAF. This new relationship between geodetic and seismological data is in good agreement with predictions from a Limited Power Law model in which the evolution of the duration of a linear aftershock decay rate over short time results from variations in the load of the brittle upper crust.


Introduction
In the last decades, geodetic measurements have considerably improved the description of spatio-temporal properties of strain accumulation and release along faults (Savage and Burford, 1973;Sauber et al., 1986;Langbein et al., 1990;Bennett et al., 1996;Peltzer et al., 2001;Fialko, 2006). Overall, these new sets of data may now provide a range of information over time scales that approach the characteristic times of loading and discharge along faults. An important result is that, complementary to the deformation accommodated by earthquakes, aseismic deformation such as post-seismic slip (Langbein, 1990), slow earthquakes (Dragert et al., 2001) and creep (Simpson et al., 2001), may accommodate an important part of the deformation. A major challenge for seismic hazard assessment remains in coupling these different modes of deformation with different patterns of seismicity Correspondence to: C. Narteau (narteau@ipgp.jussieu.fr) McEvilly, 1999, 2004;Schmidt et al., 2005;.
Aftershocks are earthquakes of smaller magnitude that occur after an event in the region of seismogenic rupture. The aftershock rate decays with time according to a power law (Omori, 1894), and, this rate of decay following the mainshock is essential in determining physical mechanisms ruling the transition from the dynamic rupture to the relaxation phase (Kagan and Houston, 2005). Unfortunately, it is extremely difficult to evaluate the exact aftershock frequency over short times due to problems in: the identification of aftershocks in coda waves of the mainshock, overlapping aftershock records, catalog compiler overload, absence or malfunction of seismic stations close to the source zone. In this case, solutions consist of scrutinizing the high-frequency signal (Vidale et al., 2004), or analyzing aftershocks of larger magnitude which are more likely to be observed (Utsu et al., 1995). In both cases, non-power law behavior may remain at the beginning of the aftershock sequence (Narteau et al., 2002;Peng et al., 2006Peng et al., , 2007Enescu et al., 2007). A usual measure of the time delay before the onset of the power law aftershock decay rate is the parameter c of the modified Omori law (MOL), where is the aftershock rate, K is a constant, t is the elapsed time from the mainshock and p is the slope of the power law aftershock decay rate (Utsu et al., 1995). In a vast majority of cases, the c value is determined empirically and its magnitude is only discussed with respect to the artifacts cited above. Inspired by a model of static fatigue suggested by Scholz (1968), Narteau et al. (2002) proposed a physical interpretation to the parameter c of the MOL by relating its magnitude to an upper limit of the overload within the aftershock zone. Narteau et al. (2005) have tested such a prediction in southern California and have shown that, at a regional Published by Copernicus Publications on behalf of the European Geosciences Union and the American Geophysical Union.  and an exponential transition rate λ(σ 0 ) produce a limited power-law (LPL) with an exponent q=1, and two characteristic aftershock rates λ b ∼1/t b and λ a ∼1/t a over short and long times, respectively. We determine N b from ∞ 0 N (σ 0 )dσ 0 . We arbitrarily choose the following numerical values: σ a =10 bars, σ b =125 bars, 1/λ a =1 yr. Note that λ b =λ(σ b ). (b) ∂ ln ( (t)) /∂ ln (t), the local slope of the LPL on logarithmic scales for σ b values ranging from 0 to 140 bars, all the other parameters being kept constant. The time delay c∼t b before the transition from a linear regime to a power law regime decreases exponentially with respect to σ b (see λ(σ 0 ) in (a)). length scale, the evolution of the c value exhibits an asymmetry which may be related to the classical seismic cycle picture (i.e. slow loading and rapid discharge).
In the present paper, we study the central segment of the San Andreas fault (SAF) where, instead of large M>7 earthquakes, the seismicity is characterized by small earthquakes along the creeping segment, and M≈6 earthquakes near Parkfield along the transition zone between the creeping and locked segments of the fault. Most of these M≈6 earthquakes occurs on the SAF, as it was the case in 1857, 1881, 1901, 1922, 1934, 1966and on 28 September 2004(Bakun et al., 2005. Nevertheless, the 1983 M 6.5 Coalinga earthquake in the west, and the 2003 M 6.5 San Simeon earthquake in the east indicate that, across the boundary between the American plate and the Pacific plate, the deformation is distributed on a population of faults and on off-fault structures (Titus et al., 2005). We consider the spatial distribution of seismic events and the systematic occurrence of aftershocks to provide, through a model of static fatigue, an estimation of loading and unloading rates across the SAF system.

A band-limited power law model of aftershock decay rate
In order to model the aftershock decay rate, let us describe briefly the relaxation mechanism which is fully developed in Narteau et al. (2002) and generalized to any stress history in Narteau (2007).
In the aftershock zone, we consider a finite population of domains which are able to fail under continued static overload. Each failure results of aging accelerated by stress perturbations induced by the mainshock and changes in the strength of the rock. Assuming that each domain can produce a single aftershock, the characteristic time of this event is given by a relationship between the overload σ 0 and failure rate λ. Classical expressions of λ(σ 0 ) are based on subcritical crack growth experiments (Atkinson and Meredith, 1987) and empirical relationships between the stress intensity factor and the crack velocity. In all cases, the failure rate is an increasing function of the stress perturbation with either an exponential (Charles and Hillig, 1962;Wiederhorn and Bolz, 1970) or a power-law (Charles, 1958;Atkinson, 1984) form. In addition to this rate, N(σ 0 , t), the overload distribution over the population of domains is enough to evaluate the aftershock rate (t) in a probabilistic manner: For any single overload value, the aftershock rate is Substituting the solution of Eq. (3) into Eq.
(2) yields where N (σ o )=N(σ o , 0) is the overload distribution just after the mainshock. The aftershock decay rate is therefore a sum of exponential decay rates, and Narteau et al. (2002) have shown that various failure rates λ(σ 0 ) and simplified overload distributions N (σ 0 ) result in the same formula that has been called the Limited Power Law (LPL) In this formula, t is the elapsed time since the mainshock, A is a constant, is the incomplete Gamma function, and λ b and λ a are two characteristic aftershock rates (Fig. 1a): -λ a =λ(0) corresponds to a threshold of crack growth at low stress levels for which strengthening processes may dominate, and prevent the propagation of the rupture (Cook, 1986).
In this model, a large range of power law decay rates result from different shapes of N (σ 0 ) and λ(σ 0 ) (see example with q=1 in Fig. 1, and Narteau et al. (2002) for examples with q<1 and q>1 associated with exponential and power law relationship for λ(σ 0 ), respectively). However, the main characteristic of the LPL is to limit the power law aftershock decay rate by an exponential regime over long times t>t a (t a ∼1/λ a ) and a linear decay rate over short times t<t b (t b ∼1/λ b ). Transition from one regime to another can be related to physical properties of the brittle layer where the aftershock sequence occurs. First, in Narteau et al. (2002) and Narteau et al. (2003), we have shown that the exponential cutoff in the power-law scaling over long time may be related to structural properties of the fractured medium. Second, over short times, the upper limit of the overload distribution can be directly estimated from the time delay before the onset of the power law regime as illustrated in Fig. 1b. It follows that assuming q=p=1 and λ a → 0 in Eqs. (1) and (5) (this is the case in Fig. 1a), we have www.nonlin-processes-geophys.net/15/245/2008/ Nonlin. Processes Geophys., 15, 245-263, 2008 Table 1. Window algorithm for aftershocks. Note that the spatial window is larger than the original window as suggested by Gardner and Knopoff (1974). We verify that the L value and T value do not significantly influence our results. at t=0. In this study, we will work under such conditions to link the LPL to the MOL and to discuss at any stage either the λ b value or the c value in relation with the σ b value (Fig. 1b).
What is the physical meaning of the upper limit of the overload distribution in the brittle crust? Locally, where the rupture stops, compositional, structural and thermal heterogeneities probably play an important role in maintaining the stress level below the yield strength of the rock (Wesnousky, 2006). In addition, aftershocks are expected to occur preferentially on pre-existing fractures distributed over a wide area. Therefore, the aftershock generation process may not be dictated solely by the complexity of earthquake slip and microscopic details of the stress redistribution in the vicinity of the rupture tips. It may also depend on prestress patterns and time-dependent behavior of a discrete population of fractures further away from the rupture. Taking a step function with a maximum threshold ( Fig. 1), we model such a population of seismic sources that have not broken during the mainshock event by following the degradation of their strength over time due to stress. Hence, the rupture is delayed according to a minimum time delay which represents the remaining strength. Such a delay is likely to be highly dependent on the magnitude thresholds for mainshocks and aftershocks. For the specific implementation developed below, this is why we focus only on large aftershocks of intermediate size mainshocks (Utsu et al., 1995;Shcherbakov et al., 2004). The main idea behind our analyzes is that, at a regional length scale, when averaged on a representative sample of aftershock sequences triggered by mainshocks in the same magnitude range, the variation of the upper limit of the overload distribution should reflect variations in the load (i.e. in the mean value of the overload distribution itself). This load being essentially affected by tectonic motions and major seismic or aseismic events, the evolution of a spatially averaged σ b value (expressed by either λ b or c) could allow for a better understanding of the accommodation of deformation along the plate boundary. Again, for large magnitude earthquakes producing higher stress perturbations, temporal variations of the σ b value due to tectonic loading are more likely to be isolated if we study only small magnitude events. Let us now present how we estimate such a macroscopic parameter from seismicity catalogs.

The average aftershock decay rate within the first day
We extract mainshocks from the U.S. Advanced National Seismic System (ANSS) composite catalog according to the algorithm of Gardner and Knopoff (1974) (Table 1). We select M Min M <M<M Max M mainshocks, and events occurring over 10 d prior to an earthquake of greater or equal magnitude at a distance shorter than 50 km are excluded from the analysis as they may represent potential foreshocks. For all the remaining mainshocks, we record the corresponding aftershock sequence within 1 d and a 40 km diameter circle.
In order to avoid artifacts arising from overlapping records, we do not consider large earthquakes and their aftershock sequences. Thus, we analyze events outside what is usually called the aftershock zone of largest events (Wiemer and Katsumata, 1999 is small and aftershocks under consideration are always relatively large in comparison to their mainshocks (Utsu et al., 1995); second, they are not in the zone of highest seismicity when they occur. Most importantly, using intermediate magnitude mainshocks, aftershocks are distributed in the entire seismic zone and the resulting catalog of aftershocks is the best available sampling of the seismicity of an entire area for a given period of time.

Statistical properties of selected aftershocks
In this section, we study the statistical properties of aftershocks that have been selected by our procedure starting with M Min A =1.0. Overall, we try to verify that the frequency-size distributions of the selected aftershocks according to the Gutenberg and Richter (1944) relationship where M and N are the magnitude and the number of earthquakes respectively. The b value is estimated by a maximum likelihood method in a magnitude range which depends on M c , the magnitude of completeness of the catalog. This magnitude of completeness is evaluated according to the procedure suggested by Wiemer and Wyss (2000) using a threshold of 10% for the residual fit between observed and predicted cumulative number of events. This method consists in Figure 3(a) shows for each of these catalogs the cumulative and non-cumulative frequency-size distributions as well as the maximum likelihood fit of the cumulated frequencysize distribution. Figure 3(b) shows the magnitude of aftershocks and the magnitude of completeness for logarithmic time periods with respect to the elapsed time from the mainshock (i.e. time periods with the same width in logarithmic scale). In all cases, M c ≤1.8 despite some fluctuations due to the small number of events over short times. In addition, frequency-size distributions exhibit a power-law behavior which is almost the same in all catalogs. Figure 4 shows the cumulative and non-cumulative frequency-size distributions as well as the maximum likelihood fit for logarithmic time periods with respect to the elapsed time from the mainshock. The power-law regime is persistent over short times, and the b value is in the vast majority of cases stable over the different time periods despite the difference in the number of aftershocks under consideration.
From the comparison between Figs. 3 and 4, we conclude that the catalogs of aftershocks obtained by our selection procedure does not exhibit systematic bias due to the incompleteness of catalogs for M≥1.8. Furthermore, these figures show that the slope of the frequency-magnitude relationship is not only stable in different subregions over various time periods but also stable over logarithmic time periods after the mainshock. Consequently, as in Narteau et al. (2005), M Min A =1.8 and M =2.7 are used as default values in our procedure of aftershock selection. 3.2 Estimation of the time delay before the onset of the power law decay rate Using a M Min A value determined above, our procedure of selection is repeated every two months for past events occurring over a period of two years (i.e. sliding windows of T w =2 yr with time shift of 2 months). Then all aftershock sequences are stacked by sorting each event according to the time interval from its respective mainshock. Artifacts resulting from catalog compiler overload should be significantly attenuated when averaged over such long times. Similarly, administrative and technical artifacts can be considered as noises of different natures which should be reduced by stacking of sequences over long time periods and over large areas. Nevertheless, as a variable at time t is estimated from events occurring in [t− T w ; t], the response to a perturbation may be shifted forward in time, particularly if there is not a significant increase in seismicity associated with such a perturbation (i.e. few new selected aftershocks). Before estimating the parameters of the MOL and the LPL, we emphasize a fundamental feature of the model: the spatial distribution of the selected mainshocks is broad and covers the entire fault system under consideration (Fig. 5); furthermore, before and after a large earthquake, the distribution of selected mainshocks exhibits strong similarities despite a zone of dense seismicity around the rupture (Fig. 5a). Hence, we analyze macroscopic properties on a disperse population of seismic events.
As shown in Sect. 2 for an individual sequence, the power law regime of the LPL results from a large number of exponential decay rates which overlap with one another according to the overload distribution just after the mainshock (see Eq. (4)). On the basis of this summation, combining different aftershock sequences in a stack is a natural and relevant extension of this model for describing of the state of stress and strength of an entire seismic zone.
By selecting aftershocks through our time window approach, we end up with a bimonthly average aftershock decay rate over one day and we investigate the onset of the power law regime. From the stacked catalogs with more than 40 events, a best-fitting procedure using the method of maximum likelihood is devoted to estimating the parameters {K, c, p} of the MOL (Eq. (1)) and the parameters {A, q, λ a , λ b } of the LPL (Eq. (5)). For a sequence with N aftershocks occurring at time t j , j ∈ [1, ..., N ] within a [t 1 , t 2 ] time interval, the maximum likelihood function for Eqs. (1) and (5) is The parameters are estimated via a method of continuous minimization by simulated annealing (Press et al., 1992), which is more likely to converge to the true global maximum than classical gradient methods. As mentioned above, we consider p=q=1 and λ a →0 (Fig. 1) in order to facilitate the evaluation of c values and λ b values (see Eq. (7)) as well as the comparison between each of these parameters over different time periods.

The onset of the power law regime across the San Andreas fault system
As an example of the data we are dealing with, Fig. 6 shows the average aftershock decay rates over six different periods of time and the best fit provided by Eqs. (1) and (5). Before interpretation of these results, the quality of fit of the MOL and the LPL should be assessed. In fact, despite the relationship between c and λ b when t→0, the LPL and the MOL do not behave identically. We compare these models by calculating AIC values, the difference between their Akaike Information criterion where n p is the number of parameters for a given model (3 for the MOL and 4 for the LPL). For all time periods since 1985 in southern California, Fig. 7 shows mostly smaller AIC values for the LPL than for the MOL; the additional parameter in the LPL is already taken into account (Eq. (11)). Nevertheless, the differences in AIC remain small. We take this as an indication that the performance of the two models is approximately the same. Because the AIC is known to favour complex models over simpler ones and applicable for nested models only, a rigorous statistical discrimination is not possible here. However, it is not the aim of the paper to pursue this.
More than the power law regime and the linear regime, the main advantage of the LPL is to provide a better estimation of the aftershock rate during the transition from one regime to another. Such a transition is the cornerstone of the present paper, and, in the following, we will focus on the evolution of λ b and its correlation with patterns of seismicity (Fig. 8) Fig. 6. The average aftershock decay rate within the first day for six periods of time before and after transient slip events along the SAF and the last M6 Parkfield earthquake. Lines represent the best fits provided by the LPL (solid) and the MOL (dashed). Note that the increase in the time delay before the onset of the power decay rate corresponds also to a decrease of the aftershock rate for t→0.

Southern California
Let us first recall some results obtained for earthquakes located in southern California between 31 • and 35 • N, and 240 • and 246 • E over the last 20 yr (Fig. 5). Narteau et al. (2005) observed that λ b decreases suddenly after large earthquakes, and slowly increases at a constant rate during periods of low seismicity (Fig. 8a). This observation was verified at smaller length scales of few hundreds of km in the neighborhood of the Landers, Northridge, and Hector Mine mainshocks. This asymmetric behavior has been related to the seismic cycle picture and, using the correlation with the Benioff strain accumulated over the moving time window (see Fig. 4a in Narteau et al., 2005), characteristic patterns of loading and discharge have emerged: 1. During interseismic periods, tectonic forcing increases the load within the system. As a result, σ b , the upper limit of the perturbations induced by 2.5<M<4.5 earthquakes increases at the same rate in logarithmic scale (log 10 (λ b )∼t/t 0 with t 0 ≈2.78 yr), and the delay before the onset of the power law decay rate decreases (i.e. λ b value increases).
2. After a large earthquake, dissipation, relaxation and stress transfer processes are associated with the rupture and its aftermath (Nur and Booker, 1972;Deng et al., 1999). Almost immediately, σ b , the upper limit of the perturbations induced by 2.5<M<4.5 earthquakes collapses to a smaller value, and the delay before the onset of the power law decay rate increases (i.e. λ b value decreases).
Are these observations similar and these inferences valid to the North where the SAF system and the seismicity exhibit different types of structural and spatiotemporal patterns?

Central California
In central California, in order to avoid complexities arising from slip partitioning between the SAF and the Calaveras fault, we analyze only earthquakes located south of the branching point (the central California zone in Fig. 9a). In this region, Fig. 6 shows that the time delay before the onset of the power law decay rate is not constant, but increases, at least since 1993. For the last 20 yr, the evolution of the λ b value has some characteristic features which can be compared to the Benioff strain in Fig. 8b. In addition, Fig. 9d shows the evolution of λ b value on a smaller region centered on the SAF (the SAF zone in Fig. 9a). The comparison between Fig. 9c and d shows that λ b values behave stably and that this behavior originates from the SAF. Unfortunately, more detailed analysis are impossible elsewhere in central California because of the small number of seismic events outside of the SAF zone. In the appendix, we explore the parameter space of the aftershock selection procedure to evaluate the statistical significance of the behaviors observed in Fig 9. In particular, we show that all the results discussed below are valid for more constraining declustering methods and magnitude thresholds for mainshocks and aftershocks.
Starting from a relative small value, λ b is continuously increasing between 1987 and 1990 (log 10 (λ b )∼t/t 0 with t 0 ≈6.67 yr in Fig. 9c and t 0 ≈11.1 yr in Fig. 9d), before it reaches a plateau. This relatively high value is conserved for 4 yr at a regional length scale. In the SAF zone, the λ b value continues to increase at the same rate. At the beginning of 1995, the λ b value starts to decrease at a constant rate until the end of 1998 (log 10 (λ b )∼−t/t 0 with t 0 ≈4.76 yr in Fig. 9c and t 0 ≈3.23 yr in Fig. 9d) Some of this behavior reinforces the hypothesis that the evolution of the λ b value is strongly dependent on seismic discharge and, in some way, correlated to the unloading rate across the SAF system: 1. the constant rate increase between 1988 and 1990 corresponds to a period of time with very low seismicity (M<4).
2. the plateau between 1991 and 1995 is associated with a series of 3<M<5 earthquakes, especially in the vicinity of Middle Mountain.
3. the collapse of the λ b value is associated with the San Simeon and the Parkfield earthquakes.
However, the constant rate decrease between 1995 and 1998 cannot be associated with any seismic pattern. Is this feature related to another type of discharge mechanism?
In fact, across the SAF system in central California, the moment rate calculated from the slip rate distribution cannot be neglected when compared to the moment rate calculated from seismic events. Furthermore, during the 1990's, if the seismicity is relatively low in central California, the seismic behavior along the SAF is evolving and, following Langbein et al. (1999), transient slip events have been quantitatively identified and measured by different groups.
From the analysis of more than a decade of high quality data, particularly those from the two-color electronic distance meter in the Parkfield area, Gao et al. (2000) suggest that the SAF underwent two transient phases of slip summarized in Table 2: a transient decrease in slip rate of about 1.5 mm/yr between 1991 and 1993; a transient increase in slip rate of about 3.3 mm/yr between 1993 and 1998. The standard deviation associated with the slow down of slip rates makes the observation much more controversial than the subsequent increase. Langbein et al. (1999) and Gao et al. (2000) note that this increase of slip rate occurs just after a sequence of 3<M<5 earthquakes in the vicinity of Middle Mountain, but last for few years afterward.
Recently, such triggering has also been identified by a time-dependent inversion done by Murray and Segall (2005). More importantly, this work has provided new information about spatiotemporal properties of aseismic deformation on the SAF at Parkfield. Spatially, from a fault plane model 40 km long with a seismogenic depth of 14 km, they have shown that the slip rate may have reached 49 mm/yr northwest of Carr Hill. This slip rate is not only higher than the slip rate of 15 mm/yr predicted between 1986 and 1990 but also higher than the long-term geological rate of 39 mm/yr between the American and Pacific plates. Temporally, Murray and Segall (2005) Gao et al. (2000) from the examination of the Parkfield deformation dataset, particularly those from the two color electronic distance meter. V 0 is an apparent velocity estimate and V 1,2 =V 0 −V 1,2 are transient velocities.
October 1992 to July 1996 (Fig. 9b). The moment release during this time period is equivalent to a M w 5.6 earthquake, a major event for this segment of the SAF. Then, they conclude that the transient aseismic event relaxed more slip than required to dissipate the stress perturbation induced by the triggering seismic events. Thus, the transient slip corresponds to an effective release of the strain stored along the SAF. These results inferred from geodetic observations can be correlated in a straightforward manner with the evolution of the λ b value: 1. The acceleration of the slip rate observed by Gao et al. (2000) overlaps almost completely with the constant decay rate of the λ b value ( Fig. 9b and Table 2). Indeed, the λ b values estimated from time windows encompassing only earthquakes that occurred during the period of higher slip rate (i.e. 1993-1998) extend from 1995 to 1998 (i.e. the time period during which the decrease is the most significant). Fig. 9c and d, we can observe than the short period of high slip rate determined by Murray and Segall (2005) triggers a decrease of the λ b value. As above, taking into account the duration of the time window, the transient slip rate is synchronized with the decrease of the λ b value. Similarly, during a short period of low seismicity (i.e. M<4), a constant rate increase of the λ b value between 1987 and 1990 corresponds to a constant slip rate of approximately 15 mm/yr.

As shown in
More generally, from the comparison between Fig. 8a and b we can observe that, on average, the λ b value is higher in central California than in southern California. In addition, the logarithmic slope of the variation of the λ b value is always lower in central California than in southern California. Since the characteristics of fatigue failures may be related to a given rheology, the behavior exposed above in term of load can be translated in terms of stress and strain.

Stress changes and loading rates inferred from aftershocks
In fracture mechanics, classical expressions of near field stress redistribution are proportional to the remote applied stress (Atkinson and Meredith, 1987). Considering a set of pre-existing domains that did not ruptured during the earthquake (i.e. the stress change is not infinity), this positive correlation allows to relate the evolution of a spatially averaged λ b value to variation in the load. Fig. 10 illustrates schematically how a mean value of the overload distribution, i.e.
can be related to the upper bound σ b of the overload distribution at the length scales of the 2.5<M<4.5 mainshocks that we consider. In addition, it shows how the µ σ 0 value evolves according to different modes of deformation in the upper crust of the Earth. Based on Fig. 10, let us discuss the loading rates in the framework of the LPL model (Fig. 1a). If λ(σ 0 ) is an exponential transition rate, and λ b =λ(σ b ), the gradual increase of the λ b value in southern California during interseismic periods can be related to a constant increase of the σ b value Assuming uniform strength over the entire population of domains, this upper limit of the overload distribution is proportional to an absolute level of differential shear stress (i.e. µ σ 0 in Fig. 10). For a linear elastic rheology, it follows that the rate α in Eq. (13) can be related to a strain accumulation rate.
In other words, we conclude that, in stick-slip zone, the time derivative of the logarithm of the λ b value is proportional to the interseismic strain accumulation rate (Fig. 10). Applied to central California, where aseismic slip accommodates a more significant part of the deformation, we find that the maximum strain accumulation rate between 1987 and 1990 is approximately three times lower than in southern California (see the constant slopes of λ b (t) in Fig. 9c  and d). From 1991 to 1993, a plateau of λ b (t) tends to indicate a balance between strain release and strain accumulation rates. Later, strain release dominates: through aseismic deformation between 1995 and 1998, and then through two major seismic events for this region (San Simeon and Parkfield earthquakes). We conclude that, in the creep-slip zone, the time derivative of the logarithm of the λ b value is not only proportional to the strain accumulation rate but also inversely proportional to the rate of strain release (Fig. 10).
Coupled with a model of interseismic deformation, the behavior of the λ b value may provide more quantitative assessments of the stress changes in central California. From the exponential transition rate presented in Fig. 1a, we have On the other hand, for a linear elastic rheology, we can write where dσ/dt,ǫ and G are the variation rate of the differential shear stress, the shear strain accumulation rate, and the shear modulus, respectively. By considering a simple screw dislocation model with a vertical strike-slip fault locked at depth D, Savage and Burford (1973) suggested that, on the free surface in x (the distance perpendicular to the fault) the shear strain accumulation rate can be related to a deficit of slip rate v s bẏ It follows that After the substitution of the time derivative of the differential shear stress (Eq. (17)) into Eq. (14), we obtain This equation allows us to relate the results exposed in Figs. 8 and 9 to different measures of slip rates across the SAF system. Summarized in Table 3, we discuss the 6 points derived from the evolution of the λ b value (see the t 0 value in Sect. 4 and in Fig. 9) and from the slip rate inferred from geodetic data. At a regional length scale, W ≈200 km, three significant trends in the variation of the logarithm of the λ b value can be estimated: 1. During interseismic periods in southern California, d log 10 (λ b )/dt=0.36 yr −1 , while the slip rate is assumed to be zero, the deformation being essentially accommodated by seismic events. As a consequence, the deficit of slip rate is given by the long-term geological slip rate. 2. From 1987 to 1990 in central California, d log 10 (λ b )/dt=0.15 yr −1 , while the maximum slip rate along the SAF has been estimated at 15 mm/yr by Murray and Segall (2005).
3. From 1995 to 1998 in central California, d log 10 (λ b )/dt=−0.21 yr −1 , while the maximum slip rate along the SAF has been estimated at 49 mm/yr by Murray and Segall (2005).
For all the above cases, the long-term geological slip rate is taken as 39 mm/yr, the angular velocity that best describes motion of the Sierra Nevada Great Valley block relative to the Pacific plate (Argus and Gordon, 2001). At a smaller length scale of W ≈1 km in the vicinity of the SAF near Parkfield (see the SAF zone in Fig. 9), two significant trends in the variation of the logarithm of the λ b value can be estimated: 4. From 1987 to 1990, d log 10 (λ b )/dt=0.09 yr −1 . During this time interval, the maximum slip rate along the SAF has been estimated at 15 mm/yr by Murray and Segall (2005).
5. From 1995 to 1998 in central California, d log 10 (λ b )/dt=−0.31 yr −1 . During this time interval, the maximum slip rate along the SAF has been estimated at 49 mm/yr by Murray and Segall (2005).
For all these cases, the long-term geological slip rate is taken equal to 28 mm/yr, as suggested by Titus et al. (2005) from continuous GPS measurements between pairs of sites that flank the creeping segment at intersite distances of 1 km.
Furthermore, we consider the null hypothesis: 6. For slip rates equal to the long-term geological slip rate, the deficit of slip rate is null, and we consider that d log 10 (λ b )/dt is equal to zero.
Gathering all these points together, Fig. 11 shows the relationship between d log 10 (λ b )/dt and the deficit of slip rate v s . There is a clear linear relationship, and the best fit straight line of slope s=0.01 mm −1 , obtained by regression, is also shown on the figure. Such a positive slope confirms that the evolution of the λ b value is likely to result from variations in the load of the brittle upper crust.   Murray and Segall (2005) and (c) the evolution of the λ b value in central California and (d) along the SAF near Parkfield (zones are shown in (a)). In (c) and (d), solid lines limit λ b values calculated from aftershocks that occurred only between October 1992 and July 1998. Dashed lines limit λ b values calculated over a time period that incorporate the maximum of slip rate. Green lines are regression lines that follow log 10 (λ b )∼t/t 0 with t 0 =6.67 yr from 1987 to 1990 and t 0 =−4.76 yr from 1995 to 1998 for (c) and with t 0 =11.1 yr from 1987 to 1990 and t 0 =−3.23 yr from 1995 to 1998 for (d). Table 3. Relationship between the evolution rate of the λ b value in logarithmic scale and the deficit of slip rate v s =v g −v. W is the characteristic length scale of the zone under consideration perpendicular to the SAF plane, v is the maximum slip rate estimated from geodetic measurements (Murray and Segall, 2005), and v g the long term geological slip rate (Titus et al., 2005).

Region
W ( Fig. 10. Schematic representation of the relationship between different deformation modes and λ b (or c) a statistical parameter extracted from catalogs of seismicity. Note thatǫ andǫ G are the strain accumulation rate and the long-term geological strain rate respectively. µ σ 0 is the mean value of the overload at a given time. The σ b value is a function of the λ b value (see Fig. 1). Hence, assuming that σ b is positively correlated with µ σ 0 , the evolution of the λ b value may be connected to loading/unloading rates in the brittle upper crust, especially where aseismic deformation transients occur.
It follows that, if no deformation is accommodated on faults or on off-fault structures, the normalizing stress constant σ a can vary with respect to the distance from the SAF. Then, From this equation and the exponential expression of λ(σ 0 ), we can simply quantify stress variations within the upper crust, σ b , in response to change in the λ b value: For example, for the gradual decrease of the λ b value along the SAF near Parkfield from 1995 to 1998 (Fig. 9d) with G=3 10 5 bars, D=14 km, s=0.01 mm −1 (see Fig. 11).
Far from the fault, our estimation is on the same order of magnitude as Coulomb stress changes induced by M≈6.5 earthquakes like Coalinga or San Simeon earthquakes over the same distances. In the vicinity of the SAF, the change in stress due to the transient slip event is on the order of magnitude of the lower range limit of earthquake stress drop (Hanks, 1977). This stress variation is also consistent but smaller with observation of low stress drop events at a border between locked and creeping fault patches (Sammis and Rice, 2001 11. Relationship between dlog 10 (λ b )/dt and the deficit of slip rate between the long-term geological rate and the slip rate across the SAF. Circle and square symbols corresponds to zones with widths of 200 km and 1 km centered on the SAF (see Figs. 8 and 9). Taking into account the increase of right lateral deformation rate with distance from the fault, we consider geological slip rates of 39 mm/yr for circles and 28 mm/yr for squares (see text). The red line is the best fit straight line obtained by regression (r 2 =0.974). The green lines represent the 95% prediction interval of the regression line.
in stress average over an entire volume within which the aftershocks take place, it is likely that local variations of stress can be much larger especially along pre-existing fractures and discontinuities. However such a value of 0.3 bars is significant: a discharge rate of 0.1 bars/yr (0.3 bars during 3 yr) is comparable with classical estimate of the loading rate along faults. In other words, transient slip events can significantly delay the recurrence time interval for earthquakes. From Eq. (20) and the exponential expression of λ(σ 0 ), we can also quantify the loading rate along the SAF (i.e. x=0): For example, in southern California, where the SAF is not creeping and where d log 10 (λ b )/dt=0.36 yr −1 , we obtain a loading rate equal to 0.117 bars/yr. For this region, paleoseismic and geological observations have shown that the mean recurrence time of major events along the SAF is about 250 yr (Sieh, 1984). Over such an interseimic period, the loading rate estimated above yields to change in stress of about 30 bars. This value is in good agreement with stress drop associated with large interplate earthquakes (Scholz, 1990).
These numerical investigations demonstrate that our method, which is only based on the examination of the aftershock decay rate and on geodetic measurements, can lead to realistic estimates of stress variations and loading rates within the brittle upper crust. Hence, the evolution of the λ b value may provide useful information for recognizing characteristic patterns of strain accumulation and release across seismic and aseismic fault systems.

Discussion and conclusion
It is likely that earthquakes are strongly under-reported during early parts of aftershock sequences, and the c value may be significantly influenced by non-physical effects. Unfortunately, it is impossible to assess the (in)completeness of the catalogs of seismicity todays, despite continuously improving techniques in seismic recording and new types of analyzes of initial phases of seismograms (Peng et al., 2007;Enescu et al., 2007). For example, here, the λ b value is saturated at high frequency over short time (e.g. 7 min and 2 min in southern and central California, respectively). Nevertheless, we consider that the origin and the variation of the c value have to be examined quantitatively in relation with independent observations (e.g. the slip rate along faults). As we are aware of possible artifacts, we try to minimize them by investigating small magnitude events over a large area and stacking different aftershock sequences occurring over a long period of time. This averaging technique is key in providing information at length scales of the entire fault system.
The suggested relationship between the time derivative of the λ b value and the deficit of slip rate may be tested in various types of tectonic settings worldwide. In such analysis, as in this paper, where absolute values being too dependent on magnitude thresholds, only relative variations have to be investigated. We have described a pattern of seismicity that allows the identification of transient slip events. In the future, this could be an opportunity to estimate acceleration or deceleration of slip across remote fault zone where geodetic measurements remains impossible.
Recently, using the same catalogs of seismicity that we use to determine the λ b value,  and  have suggested that the power law exponent of magnitude-frequency distributions (the b value) is directly related to the differential stresses in the Earth's crust. We analyze the aftershock decay rate to infer the same type of relationship between the delay before the onset of the power law regime and a measure of loading and unloading rates across a fault system. In addition, by comparing geodetic measurements to our seismological data, the LPL and a simple screw dislocation model allow for the quantification of stress variations over long times along the SAF.
We have estimated for the non-creeping section of the fault that the loading rate is on the order of 0.117 bar/yr. For the creeping section of the fault, in central California, we have shown that changes in stress for transient slip events are on the same order of magnitude as static stress changes associated with earthquakes. In the vicinity of the SAF the changes in stress tend to the lower limit of earthquake stress drop. Such a discharge is not instantaneous, but it occurs over a period of 3 yr along a 60 km segment of the SAF north-west of Parkfield. An estimated unloading rate of 0.1 bars/yr can compensated for an equivalent time interval of strain accumulation and then potentially delay by more than 6 yr the occurrence of the next M 6 earthquake in Parkfield.
As a conclusion, we suggest that, the time delay before the onset of the power-law aftershock decay rate inferred from catalogs of seismicity could provide an independent constraint on loading/unloading rates across active fault systems. In particular, if the long term geological slip rate is known, one may quantify in real time the strain accumulation and release rates as well as the ratio between aseismic and seismic deformation. Furthermore, stress changes in the upper crust can be evaluated according to a set of observations which do not rely on specific geometrical constraints of the fault population. Carlo method is used to evaluate the uncertainty of the maximum likelihood estimate of λ b (black crosses). The error bar curve corresponds to 16% and 84% quantiles of the maximum likelihood estimates of λ b for 500 synthetic aftershock sequences at each point.

Appendix A
Independent fits of c and λ b Figure 6a shows the evolution of c and λ b obtained by independent fits of the MOL and the LPL respectively in the central California zone from 1984 to 2005. As predicted by Eq. (7), the comparison between the evolution of c and 1/λ b indicates that the MOL and the LPL give similar estimates of the time delay before the onset of the power law decay rate. It is important to note that the aftershock decay rate for both laws can differ significantly from one model to the other, especially during the transition period toward the power law regime (see Eqs. (1) and (5), and Fig. 6).
Uncertainty of maximum likelihood estimates Figure 6b shows the uncertainty of maximum likelihood estimates (MLE) of λ b using a Monte Carlo approach. Practically, at each time step, 500 independent aftershock sequences are generated using a non-stationary Poissonian process with a frequency determined by the LPL and the MLE of K and λ b (remember that q=1 and λ a →0). An individual sequence lasts from t=10 s to t=1 d. Then, for each    Table 1. The black dotted lines indicate the trends estimated in Fig. 9 for periods of 1987-1990. of them, we obtain new MLE of λ b and evaluate the 16% and 84% quantiles of their distribution. Error bars between these quantiles remain small and the envelope curves reflect the same behaviors than the initial MLE of λ b (Fig. 6b). In addition, the uncertainty decreases with the number of observations, and it is much lower in southern California where more aftershocks are stacked together at each time step.

Appendix B
Our analysis is based on a selection procedure of aftershocks that requires few input parameters particularly for the spacetime windows and the magnitude thresholds. Figure 13 shows the effect of these parameters on the number of selected aftershocks and on the variation of λ b in the central California and the SAF zones from 1984 to 2005. For comparison with the results presented in Sect. 4.2, dotted lines indicate the trends calculated in Fig. 9 for the periods of 1987-1990 and 1995-1998. In all figures, the quality of the fit and the stability of the results are highly dependent on the number of staked aftershocks. Tests have shown that reliable estimates are obtained with stacks consisting of a minimum of 40 aftershocks. Note on Fig. 13 that below such a threshold, the signal investigated is still present, but shows a significant increase in noise levels (the black curves).

Magnitude thresholds
In order to analyze the properties of aftershock sequences over short times, the classical procedure is to eliminate events of smaller magnitude, larger events being identified more easily in seismograms. Utsu et al. (1995) suggest that the time delay before the onset of the power law aftershock decay is not an artifact if it converges to a constant value for an increasing M Min A value. Here, given our selection procedure with magnitude thresholds for mainshocks and aftershocks, such a test can be done by decreasing the M value (Eq. (8) Fig. 9c and d where M Min A =1.8, M Max M =4.5) despite strong fluctuations when the number of aftershocks in the stack is too low. These results indicate that there is no significant bias associated with the magnitude thresholds in our aftershock catalogs. We emphasize that it is because we are only using the largest events of intermediate magnitude mainshocks (minimization of M ).
Time window duration Figure 13b shows the evolution of λ b for different time window durations, with time steps of two months. The main characteristic of a time window is to reduce the level of noise by averaging a number of consecutive measurements over time: the shorter the duration of the time window, the higher the level of noise. This is the case here since the number of events in the stack is correlated to the duration of the time window. Practically, we choose T w =2 yr because this is the shortest time window which always gives a number of stacked aftershocks larger than 40. Nevertheless, in all other cases, the increase and decrease rates of λ b remain very similar for the periods of 1987-1990 and 1995-1998, respectively.

Space-time windows for aftershock selection
The effect of the algorithm of Gardner and Knopoff (1974) ( Table 1) is tested by multiplying and dividing the space-time windows by 2 and 3, respectively (Fig. 13c). Increasing windows by a factor 2 diminishes the number of aftershocks in the stacks and consequently increases the variability of the signal. Reducing windows by a factor 3 has no impact on the result. On the other hand, a larger number of events in stacks offers the possibility of increasing the M Min A value (i.e. of decreasing the M value). In this case, the evolution of λ b is almost the same as in Fig. 9c and d, and, interestingly, is more stable than the curve with the same magnitude thresholds in Fig. 13a (i.e. the black curve with M Min A =2.2, M Max M =4.5).

Declustering method
Finally, we modify our technique of selection itself. The algorithm of Gardner and Knopoff have only been used to select mainshocks, but, for the selection of aftershocks, the space windows scale with the magnitude M of the mainshock. We consider a circular area with a radius R=10 β(M−M Min M ) r, where β is a constant, and r is an arbitrary distance. Figure 13d shows the evolution of λ b for a fixed β value and different r values, and for a fixed r value and different β values. It shows also the results from aftershock catalogs obtained by taking the spatial parameters of the Reasenberg (1985) declustering method (β=0.41, r=1.12 km). All curves behave similarly with the exception of the largest r value in the SAF zone, where aftershocks are mixed with significant uncorrelated seismicity along the fault, which affects the temporal decay of the aftershock decay rate over time (green curve in Fig. 13d). The Reasenberg (1985) parameters give a smaller number of events than our default procedure but the shape of the evolution of λ b remains unchanged.
From Fig. 13, we can conclude that, considering large aftershocks of intermediate size mainshocks (small M value), it is possible to capture time variations of the λ b parameter. However, the number of aftershocks in the stacks provides a strong constraint. This number must be larger than 40 to ensure the quality of fit of the MOL and the LPL and to reduce statistical fluctuations.