Utsu aftershock productivity law explained from geometric operations on the permanent static stress field of mainshocks

The aftershock productivity law, first described by Utsu in 1970, is an exponential function of the form K=K0.exp({\alpha}M) where K is the number of aftershocks, M the mainshock magnitude, and {\alpha} the productivity parameter. The Utsu law remains empirical in nature although it has also been retrieved in static stress simulations. Here, we explain this law based on Solid Seismicity, a geometrical theory of seismicity where seismicity patterns are described by mathematical expressions obtained from geometric operations on a permanent static stress field. We recover the exponential form but with a break in scaling predicted between small and large magnitudes M, with {\alpha}=1.5ln(10) and ln(10), respectively, in agreement with results from previous static stress simulations. We suggest that the lack of break in scaling observed in seismicity catalogues (with {\alpha}=ln(10)) could be an artefact from existing aftershock selection methods, which assume a continuous behavior over the full magnitude range. While the possibility for such an artefact is verified in simulations, the existence of the theoretical kink remains to be proven.


Introduction
Aftershocks, the most robust patterns observed in seismicity, are characterized by three empirical laws, which are functions of time (e.g., Utsu et al., 1995;Mignan, 2015), space (e.g., Richards-Dinger et al., 2010) and mainshock magnitude (Utsu, 1970a;b;Ogata, 1988).The present study focuses on the latter relationship, i.e., the Utsu aftershock productivity law, which describes the total number of aftershocks K produced by a mainshock of magnitude M as with m 0 the minimum magnitude cutoff (Utsu, 1970b;Ogata, 1988).This relationship was originally proposed by Utsu (1970a;b) by combining two other empirical laws, the Gutenberg-Richter relationship (Gutenberg and Richter, 1944) and Båth's law (Båth, 1964), respectively: with β the magnitude size ratio (or b = β/ln(10) in base-10 logarithmic scale) and Δm B the magnitude difference between the mainshock and its largest aftershock, such that with != exp −Δm ! and = .Eq. ( 3) was only implicit in Utsu (1970a) and not exploited in Utsu (1970b) where K 0 was fitted independently of the value taken by Båth's parameter Δm B .The α-value was in turn decoupled from the β-value in later studies (e.g., Seif et al. (2017) and references therein).
Although it seems obvious that Eq. ( 1) can be explained geometrically if the volume of the aftershock zone is correlated to the mainshock surface area S with = 10 !!! = exp ln(10)( − 4) (see physical explanation of Eq. ( 4) in Kanamori and Anderson, 1975), there is so far no analytical, physical expression of Eq. (1) available.Although Hainzl et al. (2010) retrieved the exponential behavior in numerical simulations where aftershocks were produced by the permanent static stress field of mainshocks of different magnitudes, it remains unclear how K 0 and α relate to the underlying physical parameters.
The aim of the present article is to explain the Utsu aftershock productivity equation (Eq. 1) by applying a geometrical theory of seismicity (or "Solid Seismicity"), which has already been shown to effectively explain other empirical laws of both natural and induced seismicity from simple geometric operations on a permanent static stress field (Mignan, 2012;2016).The theory is applied here for the first time to the case of aftershocks.

Solid Seismicity Postulate (SSP):
Seismicity can be strictly categorized into three regimes of constant spatiotemporal densities -background !, quiescence ! and activation !(with !≪ !≪ ! ) -occurring respective to the static stress step function: with Δ * the background stress amplitude range.
Based on this Postulate, Mignan (2012) demonstrated the power-law behavior of precursory seismicity in agreement with the observed time-to-failure equation (Varnes, 1989), while Mignan (2016) demonstrated both the observed parabolic spatiotemporal front and the linear relationship with injection-flow-rate of induced seismicity (Shapiro and Dinske, 2009).It remains unclear whether the SSP has a physical origin or not.If not, it would still represent a reasonable approximation of the linear relationship between event production and static stress field in a simple clockchange model (Hainzl et al., 2010) (Fig. 1a).The power of Eq. ( 5) is that it allows defining seismicity patterns in terms of "solids" described by the spatial envelope * = = ±Δ * .The spatiotemporal rate of seismicity is then a mathematical expression defined by the density of events δ times the volume characterized by * (see previous demonstrations in Mignan et al., 2007;Mignan, 2011;2012;2016 where simple algebraic expressions were obtained).
In the case of aftershocks, we define the static stress field of the mainshock by with Δσ 0 < 0 the mainshock stress drop, c the crack radius and r the distance from the crack.Eq (6) is a simplified representation of stress change from slip on a planar surface in a homogeneous elastic medium.It takes into account both the square root singularity at crack tip and the 1/r 3 falloff at higher distances (Dieterich, 1994) (Fig. 1b).It should be noted that this radial static stress field does not represent the geometric complexity of Coulomb stress fields (Fig. 2a).However we are here only interested in the general behavior of aftershocks with Eq. ( 6) retaining the first-order characteristics of this field (i.e., on-fault seismicity; Fig. 2b), which corresponds to the case where the mainshock relieves most of the regional stresses and aftershocks occur on optimally oriented faults.It is also in agreement with observations, most aftershocks being located on and around the mainshock fault traces in Southern California (Fig. 2c; see section "Observations & Model Fitting").The occasional cases where aftershocks occur off-fault (e.g., Ross et al., 2017) can be explained by the mainshock not relieving all of the regional stress (King et al., 1994) (Fig. 2d).(Lin and Stein, 2004;Toda et al., 2005).
For * = = Δ * , Eq. ( 6) yields the aftershock solid envelope of the form: , function of the crack radius c and of the ratio between background stress amplitude range Δ * and stress drop Δσ 0 (Fig. 1c).With Δσ 0 independent of earthquake size (Kanamori and Anderson, 1975;Abercrombie and Leary, 1993) and Δ * assumed constant, * is directly proportional to c with proportionality constant, or stress factor, F (Eq. 7).Geometrical constraints due to the seismogenic layer width w 0 then yield with S the rupture surface defined by Eq. ( 4) and c becoming an effective crack radius (Kanamori and Anderson, 1975;Fig. 1d).Note that the factor of 2 (i.e., using w 0 instead of w 0 /2) comes from the free surface effect (e.g., Kanamori and Anderson, 1975;Shaw and Scholz, 2001).
The aftershock productivity K(M) is then the activation density !times the volume * () of the aftershock solid.For the case in which the mainshock relieves most of the regional stress, stresses are increased all around the rupture (King et al., 1994), which is topologically identical to stresses increasing radially from the rupture plane (Fig. 2a-b).It follows that the aftershock solid can be represented by a volume of contour * from the rupture plane geometric primitive, i.e., a disk or a rectangle, for small and large mainshocks respectively.This is illustrated in Figure 3ab and can be generalized by where d is the distance travelled around the geometric primitive by the geometric centroid of the semi-circle of radius * (i.e., Pappus's Centroid Theorem), or For the disk, the volume (Eq.9) corresponds to the sum of a cylinder of radius c(M) and height 2 * (first term) and of half a torus of major radius c(M) and minus radius * (second term).For the rectangle, the volume is the sum of a cuboid of length l(M) (i.e., rupture length), width w 0 and height 2 * (first term) and of a cylinder of radius * and height w 0 (second term; see red and orange volumes, respectively, in Figure 3a-c).Finally inserting Eqs. ( 7), ( 8) and ( 10) into (9), we obtain which is represented in Figure 3d.Considering the two main regimes only (small versus large mainshocks) and inserting Eq. ( 4) into (11), we get which is a closed-form expression of the same form as the original Utsu productivity law (Eq.1).Seismicity (Eq.11).This relationship is of the same form as the Utsu productivity law (Eq.1) for large M (see text for an explanation of the lack of break in scaling in Eq. 1 for small M).

Dotted vertical lines represent M for
Here, we predict that the α-value decreases from 3ln(10)/2 ≈ 3.45 to ln(10) ≈ 2.30 when switching regime from small to large mainshocks (or from 1.5 to 1 in base-10 logarithmic scale).It should be noted that Hainzl et al. (2010) observed the same break in scaling in static stress transfer simulations, which corroborates our analytical findings.For large M, the scaling is fundamentally the same as in Eq. ( 4).Since that relation also explains the slope of the Gutenberg-Richter law (see physical explanation given by Kanamori and Anderson, 1975), it follows that ≡ , which is also in agreement with the original formulation of Utsu (1970a;b) (Eq. 3).

Observations & Model Fitting
We consider the case of Southern California and extract aftershock sequences from the relocated earthquake catalog of Hauksson et al. (2012) defined over the period 1981-2011, using the nearest-neighbor method (Zaliapin et al., 2008) (used with its standard parameters originally calibrated for Southern California).Only events with magnitudes greater than m 0 = 2.0 are considered (a conservative estimate following results of Tormann et al. (2014); saturation effects immediately after the mainshock are negligible when considering entire aftershock sequences; Helmstetter et al., 2005).The observed number of aftershocks n produced by a mainshock of magnitude M (for a total of N mainshocks) is shown in Figure 4a.
We first fit Eq. ( 1) to the data using the Maximum Likelihood Estimation (MLE) method with the log-likelihood function for a Poisson process, or, with Eq. ( 1), (note that the last term can be set to 0 during LL maximization).For Southern California, we obtain α MLE = 2.04 (0.89 in log 10 scale) and K 0 = 0.23.It should be noted that this approach does not include the case of mainshocks that produce zero aftershock.Therefore we also compute the MLE for the Zero-Inflated Poisson (ZIP) distribution: where w is a weighting constant.It finally follows that α MLE(ZIP) = 2.13 (0.93 in log 10 scale, with K 0 = 0.15), corrected for zero-values.This result is in agreement with previous studies in the same region (e.g., Helmstetter et al., 2005;Zaliapin and Ben-Zion, 2013;Seif et al., with α = ln(10) ≈ 2.30 predicted for large mainshocks in Solid Seismicity.Moreover we find a bulk β MLE = 2.34 (1.02 in log 10 scale) (Aki, 1965), in agreement with α = β.It should be noted that no significant difference is obtained when computing β MLE for background events or aftershocks alone, with β MLE = 2.29 and 2.35, respectively (0.99 and 1.02 in log 10 scale).
We also tested the following piecewise model to identify any break in scaling, as predicted by Eq. ( 12): but with the best MLE result obtained for M break = m 0 , suggesting no break in scaling in the aftershock productivity data.

Role of aftershock selection on productivity scaling-break
We now identify whether the lack of break in scaling in aftershock productivity observed in earthquake catalogues could be an artefact related to the aftershock selection method.We run Epidemic-Type Aftershock Sequence (ETAS) simulations (Ogata, 1988;Ogata and Zhuang, 2006), with the seismicity rate Aftershock sequences are defined by power laws, both in time and space (for an alternative temporal function, see Mignan, 2015;2016b).µ is the Southern California background seismicity, as defined by the nearest-neighbor method (with same t, x, y and m).We fix the ETAS parameters to θ = {c = 0.011 day, p = 1.08, d = 0.0019 km 2 , q = 1.47, γ = 2.01}, following the fitting results of Seif et al. (2017) for the Southern California relocated catalog and m 0 = 2 (see their Table 1).However, we define the productivity K(M) from Eq. ( 16) with M break = 5, K 0 = 0.23, α = 2.04 and β = 2.3.
Examples of ETAS simulations are shown in Figure 4b for comparison with the observed Southern California time series.Figure 4c allows us to verify that the simulated aftershock productivity is kinked at M break , as defined by Eq. ( 16).
We then select aftershocks from the ETAS simulations with the nearestneighbor method.Figure 4d represents the estimated aftershock productivity, which has lost the break in scaling originally implemented in the simulations.This demonstrates that the theoretical break in scaling predicted in the aftershock productivity law can be lost in observations due to an aftershock selection bias, all declustering techniques assuming continuity over the entire magnitude range.While such an artefact is possible, it yet does not prove that the break in scaling exists.The fact that a similar break in scaling was obtained in independent Coulomb stress simulations (Hainzl et al., 2010) however provides high confidence in our results.

Conclusions
In the present study, a physical closed-form expression defined from geometric and static stress parameters was proposed (Eq.12) to explain the empirical Utsu aftershock productivity law (Eq.1).This demonstration, combined to the previous ones made by the author to explain precursory accelerating seismicity and induced seismicity (Mignan, 2012;2016), suggests that most empirical laws observed in seismicity populations can be explained by simple geometric operations on a permanent static stress field.Although the Solid Seismicity Postulate (SSP) (Eq.5) remains to be proven, it is so far a rather convenient and pragmatic assumption to determine the physical parameters that play a first-order role in the behavior of seismicity.It is also complementary to the more common simulations of static stress loading (King and Bowman, 2003) and static stress triggering (Hainzl et al., 2010).
Analytic geometry, providing both a visual representation and an analytical expression of the problem at hand (Fig. 3), represents a new approach to try better understanding the behavior of seismicity.Its current limitation in the case of aftershock analysis consists in assuming that the static stress field is radial and described by Eq. ( 6) (Dieterich, 1994), which is likely only valid for mainshocks relieving most of the regional stresses and with aftershocks occurring on optimally oriented faults (King et al., 1994).More complex, second-order, stress behaviors might explain part of the scattering observed around Eq. (1) (Fig. 4a).Other σ(r) formulations could be tested in the future, the only constraint on generating so-called seismicity solids being the use of the postulated static stress step function of Eq. (5) (i.e., the Solid Seismicity Postulate, SSP).

Figure 1 .
Figure 1.Definition of the aftershock solid envelope in a permanent static stress field: (a)

Figure 2 .
Figure 2. Possible static stress fields and inferred aftershock spatial distribution: (a) Right-

Figure 3 .
Figure 3. Geometric origin of the aftershock productivity law: (a) Sketch of the aftershock

Figure 4 .
Figure 4. Aftershock productivity, observed and simulated, defined as the number of