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

The aftershock productivity law is an exponential function of the form K∝exp(αM), with K being the number of aftershocks triggered by a given mainshock of magnitude M and α≈ln (10) being the productivity parameter. This law remains empirical in nature although it has also been retrieved in static stress simulations. Here, we parameterize this law using the solid seismicity postulate (SSP), the basis of a geometrical theory of seismicity where seismicity patterns are described by mathematical expressions obtained from geometric operations on a permanent static stress field. We first test the SSP that relates seismicity density to a static stress step function. We show that it yields a power exponent q = 1.96 ± 0.01 for the power-law spatial linear density distribution of aftershocks, once uniform noise is added to the static stress field, in agreement with observations. We then recover the exponential function of the productivity law with a break in scaling obtained between small and large M, with α=1.5ln (10) and ln (10), respectively, in agreement with results from previous static stress simulations. Possible biases of aftershock selection, proven to exist in epidemic-type aftershock sequence (ETAS) simulations, may explain the lack of break in scaling observed in seismicity catalogues. The existence of the theoretical kink, however, remains to be proven. Finally, we describe how to estimate the solid seismicity parameters (activation density δ₊, aftershock solid envelope r_∗ and background stress amplitude range Δo_∗) for large M values.

1 Introduction

Aftershocks, one of the most studied patterns observed in seismicity, are characterized by three empirical laws, which are functions of time, such as the modified Omori law (e.g., Utsu et al., 1995), space (e.g., Richards-Dinger et al., 2010; Moradpour et al., 2014) 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

$$\begin{array}{}\text{(1)}& K\left(M\right)=K_{\mathrm{0}}\exp \left[\mathit{\alpha }\left(M-m_{\mathrm{0}}\right)\right],\end{array}$$

with m₀ 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, 1965), respectively:

$$\begin{array}{}\text{(2)}& \left\{\begin{array}{l}N\left(\ge m\right)=A\exp \left[-\mathit{\beta }\left(m-m_{\mathrm{0}}\right)\right]\\ N\left(\ge M-\mathrm{\Delta }{m}_{\mathrm{B}}\right)=\mathrm{1}\end{array}\right.,\end{array}$$

with N the average number of events above magnitude m, A a seismic activity constant, β the magnitude size ratio (or $b=\mathit{\beta }/\ln \left(\mathrm{10}\right)$ in base-10 logarithmic scale) and Δm_B the magnitude difference between the mainshock and its largest aftershock, such that

$$\begin{array}{ll}{\displaystyle }K\left(M\right)& {\displaystyle }=N\left(\ge m_{\mathrm{0}}\left|M\right.\right)\\ \text{(3)}& {\displaystyle }& {\displaystyle }=\exp \left(-\mathit{\beta }\mathrm{\Delta }{m}_{\mathrm{B}}\right)\exp \left[\mathit{\beta }\left(M-m_{\mathrm{0}}\right)\right],\end{array}$$

with $K_{\mathrm{0}}=\exp \left(-\mathit{\beta }\mathrm{\Delta }{m}_{\mathrm{B}}\right)$ and α≡β. Equation (3) was only implicit in Utsu (1970a) and not exploited in Utsu (1970b), where K₀ 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

$$\begin{array}{}\text{(4)}& S\left(M\right)={\mathrm{10}}^{M-\mathrm{4}}=\exp \left[\ln \left(\mathrm{10}\right)(M-\mathrm{4})\right]\end{array}$$

(Kanamori and Anderson, 1975; Yamanaka and Shimazaki, 1990; Helmstetter, 2003), 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₀ and α relate to the underlying physical parameters.

The aim of the present article is to describe the Utsu aftershock productivity equation (Eq. 1) in terms of a geometrical theory of seismicity coined “solid seismicity”, where the Eq. (4) scaling is parameterized using the solid seismicity postulate (SSP). The SSP 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, 2016a). The theory is applied here for the first time to describe aftershocks.

2 Physical expression of the aftershock productivity law

2.1 Demonstration of the productivity law by geometric operations

“Solid seismicity”, a geometrical theory of seismicity, is based on the following postulate (Mignan et al., 2007; Mignan, 2008, 2012, 2016a).

Solid seismicity postulate: Seismicity can be strictly categorized into three regimes of constant spatiotemporal densities δ – background δ₀, quiescence δ₋ and activation δ₊ (with ${\mathit{\delta }}_{-}\ll {\mathit{\delta }}_{\mathrm{0}}\ll {\mathit{\delta }}_{+}$) – occurring respective to the static stress step function:

$$\begin{array}{}\text{(5)}& \mathit{\delta }\left(\mathit{\sigma }\right)=\left\{\begin{array}{ll}{\mathit{\delta }}_{-},& \mathit{\sigma }<-\mathrm{\Delta }{o}_{\ast }\\ {\mathit{\delta }}_{\mathrm{0}},& \mathit{\sigma }\le \left|\pm \mathrm{\Delta }{o}_{\ast }\right|\\ {\mathit{\delta }}_{+},& \mathit{\sigma }>\mathrm{\Delta }{o}_{\ast }\end{array}\right.,\end{array}$$

with σ the static stress (stress unit), Δo_∗ the background stress amplitude range (stress unit), a stress threshold value separating two seismicity regimes, and δ the spatial density of events (number of events per unit of volume) per seismicity regime.

We mean by “strictly categorized” that any seismicity population is either part of the background, quiescence or activation regime (or class), with no other regime or class possible (i.e., a sort of hard labeling). 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 (2016a) 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 clock-change model (Hainzl et al., 2010; Fig. 1a). For the testing of the SSP on the observed spatial distribution of aftershocks, see Sect. 2.2. The power of Eq. (5) is that it allows seismicity patterns to be defined in terms of “solids” described by the spatial envelope $r_{\ast }=r\left(\mathit{\sigma }=\pm \mathrm{\Delta }{o}_{\ast }\right)$, where r is the distance from the static stress source (e.g., mainshock rupture) and r_∗ is the distance r at which there is a change of regime (quiescence–background at $\mathit{\sigma }=-\mathrm{\Delta }{o}_{\ast }$ or background–activation at $\mathit{\sigma }=\mathrm{\Delta }{o}_{\ast })$. The spatiotemporal rate of seismicity is then a mathematical expression defined by the density of events δ times the volume characterized by r_∗ (see previous demonstrations in Mignan et al., 2007 and Mignan, 2011, 2012, 2016a where simple algebraic expressions were obtained).

In the case of aftershocks, we define the static stress field of the mainshock by

$$\begin{array}{}\text{(6)}& \mathit{\sigma }\left(r\right)=-\mathrm{\Delta }{\mathit{\sigma }}_{\mathrm{0}}\left[{\left(\mathrm{1}-{\displaystyle \frac{c^{\mathrm{3}}}{(r+c{)}^{\mathrm{3}}}}\right)}^{\frac{-\mathrm{1}}{\mathrm{2}}}-\mathrm{1}\right],\end{array}$$

with Δσ₀ < 0 the mainshock stress drop, c the crack radius and r the distance from the crack. Equation (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³ 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 Sect. 3). 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).

Figure 1 Definition of the aftershock solid envelope in a permanent static stress field: (a) event density stress step-function δ(σ) (Eq. 5) of the solid seismicity postulate (SSP) in comparison to the linear clock-change model; (b) static stress σ versus distance r for different effective crack radii c and rupture stress drops Δσ₀ (Eq. 6); (c) linear relationship between effective crack radius c and aftershock solid envelope radius r_∗ for different $\mathrm{\Delta }{\mathit{\sigma }}_{\ast }/\mathrm{\Delta }{\mathit{\sigma }}_{\mathrm{0}}$ ratios (Eq. 7); (d) relationship between mainshock magnitude M and effective crack radius c for different seismogenic widths w₀ (Eq. 8).

Figure 2 Possible static stress fields and inferred aftershock spatial distribution: (a) right-lateral Coulomb stress field for optimally oriented faults, where the mainshock relieves all of the regional stresses σ_r=10 bar, with $\mathrm{\Delta }{\mathit{\sigma }}_{\mathrm{0}}\approx -Gs/L\approx -\mathrm{10}$ bar ($G=\mathrm{3.3}\times {\mathrm{10}}^{\mathrm{5}}$ bar the shear modulus, s=0.6 m the slip, L=20 km the fault length and w=10 km the fault width); (b) radial static stress field computed from Eq. (6) with $\mathrm{\Delta }{\mathit{\sigma }}_{\mathrm{0}}=-\mathrm{10}$ bar and $c=\sqrt{\left(Lw\right)/\mathit{\pi }}$ for consistency with panel (a); (c) aftershock distribution of the largest strike-slip events in the southern California relocated catalog, identified here as all events occurring within 1 day of the mainshock (see data Sect. 3.1); (d) right-lateral Coulomb stress field for optimally oriented faults, where the mainshock relieves only a fraction of the regional stresses σ_r=100 bar with $\mathrm{\Delta }{\mathit{\sigma }}_{\mathrm{0}}=-\mathrm{10}$ bar (same rupture as in panel a) – the black contour represents 1 bar in panels (a), (b) and (d) and a 10 km distance from rupture in panel (c). Coulomb stress fields of panels (a) and (d) were computed using the Coulomb 3 software (Lin and Stein, 2004; Toda et al., 2005).

For $r_{\ast }=r\left(\mathit{\sigma }=\mathrm{\Delta }{o}_{\ast }\right)$, Eq. (6) yields the aftershock solid envelope of the following form:

$$\begin{array}{}\text{(7)}& r_{\ast }\left(c\right)=\left\{{\displaystyle \frac{\mathrm{1}}{{\left[\mathrm{1}-{\left(\mathrm{1}-\frac{\mathrm{\Delta }{\mathit{\sigma }}_{\ast }}{\mathrm{\Delta }{\mathit{\sigma }}_{\mathrm{0}}}\right)}^{-\mathrm{2}}\right]}^{\frac{\mathrm{1}}{\mathrm{3}}}}}-\mathrm{1}\right\}c=Fc\end{array}$$

function of the crack radius c and of the ratio between background stress amplitude range Δo_∗ and stress drop Δσ₀ (Fig. 1c). With Δσ₀ independent of earthquake size (Kanamori and Anderson, 1975; Abercrombie and Leary, 1993) and Δo_∗ assumed constant, r_∗ is directly proportional to c with proportionality constant, or stress factor, F (Eq. 7). Geometrical constraints due to the seismogenic layer width w₀ then yield

$$\begin{array}{}\text{(8)}& c\left(M\right)=\left\{\begin{array}{ll}{\left({\displaystyle \frac{S\left(M\right)}{\mathit{\pi }}}\right)}^{\frac{\mathrm{1}}{\mathrm{2}}},& S\left(M\right)\le \mathit{\pi }{w}_{\mathrm{0}}^{\mathrm{2}}\\ w_{\mathrm{0}},& S\left(M\right)>\mathit{\pi }{w}_{\mathrm{0}}^{\mathrm{2}}\end{array}\right.,\end{array}$$

with S the rupture surface area 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₀ instead of w₀∕2) comes from the free surface effect (e.g., Kanamori and Anderson, 1975; Shaw and Scholz, 2001).

Figure 3 Geometric origin of the aftershock productivity law: (a) sketch of the aftershock solid for a small mainshock rupture represented by a disk; (b) sketch of the aftershock solid for a large mainshock rupture represented by a rectangle; (c) relative role of the two terms of Eq. (9), here with w₀=10 km and $\frac{\mathrm{\Delta }{\mathit{\sigma }}_{\ast }}{\mathrm{\Delta }{\mathit{\sigma }}_{\mathrm{0}}}=-\mathrm{0.1}$ (to first estimate c and r_∗ from Eqs. 8 and 7, respectively); (d) aftershock productivity law (normalized by δ₊) predicted by solid 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 $c\left(M\right)+r_{\ast }\left(M\right)=w_{\mathrm{0}}/\mathrm{2}$ and $S\left(M\right)=\mathit{\pi }{w}_{\mathrm{0}}^{\mathrm{2}}$, respectively.

The aftershock productivity K(M) is then the activation density δ₊ times the volume V_∗(M) 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 r_∗(M) from the rupture plane geometric primitive, i.e., a disk or a rectangle for small and large mainshocks, respectively. This is illustrated in Fig. 3a–b and can be generalized by

$$\begin{array}{}\text{(9)}& V_{\ast }\left(M\right)=\mathrm{2}{r}_{\ast }\left(M\right)S\left(M\right)+{\displaystyle \frac{\mathit{\pi }}{\mathrm{2}}}{r}_{\ast }^{\mathrm{2}}\left(M\right)d,\end{array}$$

where d is the distance traveled around the geometric primitive by the geometric centroid of the semicircle of radius r_∗(M) (i.e., Pappus's Centroid Theorem), or

$$\begin{array}{}\text{(10)}& d=\left\{\begin{array}{ll}\mathrm{2}\mathit{\pi }\left(c\left(M\right)+{\displaystyle \frac{\mathrm{4}}{\mathrm{3}\mathit{\pi }}}{r}_{\ast }\left(M\right)\right),& c\left(M\right)+r_{\ast }\left(M\right)\le {\displaystyle \frac{w_{\mathrm{0}}}{\mathrm{2}}}\\ \mathrm{2}{w}_{\mathrm{0}},& c\left(M\right)+r_{\ast }\left(M\right)>{\displaystyle \frac{w_{\mathrm{0}}}{\mathrm{2}}}\end{array}\right..\end{array}$$

For the disk, the volume (Eq. 9) corresponds to the sum of a cylinder of radius c(M) and height 2r_∗(M) (first term) and of half a torus of major radius c(M) and minus radius r_∗(M) (second term). For the rectangle, the volume is the sum of a cuboid of length l(M) (i.e., rupture length), width w₀ and height 2r_∗(M) (first term) and of a cylinder of radius r_∗(M) and height w₀ (second term; see red and orange volumes, respectively, in Fig. 3a–c). Finally inserting Eqs. (7), (8) and (10) into Eq. (9), we obtain

$$\begin{array}{}\text{(11)}& K\left(M\right)={\mathit{\delta }}_{+}\left\{\begin{array}{l}\left[{\displaystyle \frac{\mathrm{2}F}{\sqrt{\mathit{\pi }}}}+F^{\mathrm{2}}\sqrt{\mathit{\pi }}\left(\mathrm{1}+{\displaystyle \frac{\mathrm{4}}{\mathrm{3}\mathit{\pi }}}F\right)\right]S^{\frac{\mathrm{3}}{\mathrm{2}}}\left(M\right),\\ S\left(M\right)\le {\left({\displaystyle \frac{w_{\mathrm{0}}\sqrt{\mathit{\pi }}}{\mathrm{2}(\mathrm{1}+F)}}\right)}^{\mathrm{2}}\\ {\displaystyle \frac{\mathrm{2}F}{\sqrt{\mathit{\pi }}}}{S}^{\frac{\mathrm{3}}{\mathrm{2}}}\left(M\right)+F^{\mathrm{2}}{w}_{\mathrm{0}}S\left(M\right),\\ {\left({\displaystyle \frac{w_{\mathrm{0}}\sqrt{\mathit{\pi }}}{\mathrm{2}(\mathrm{1}+F)}}\right)}^{\mathrm{2}}<S\left(M\right)\le \mathit{\pi }{w}_{\mathrm{0}}^{\mathrm{2}}\\ \mathrm{2}F{w}_{\mathrm{0}}S\left(M\right)+\mathit{\pi }{F}^{\mathrm{2}}{w}_{\mathrm{0}}^{\mathrm{3}},\\ S\left(M\right)>\mathit{\pi }{w}_{\mathrm{0}}^{\mathrm{2}}\end{array}\right.\end{array}$$

which is represented in Fig. 3d. Considering the two main regimes only (small versus large mainshocks) and inserting Eq. (4) into (11), we get

$$\begin{array}{ll}\text{(12)}& {\displaystyle }& {\displaystyle }K\left(M\right)={\displaystyle }& {\displaystyle }{\mathit{\delta }}_{+}\left\{\begin{array}{l}\left[{\displaystyle \frac{\mathrm{2}F}{\sqrt{\mathit{\pi }}}}+F^{\mathrm{2}}\sqrt{\mathit{\pi }}\left(\mathrm{1}+{\displaystyle \frac{\mathrm{4}}{\mathrm{3}\mathit{\pi }}}F\right)\right]\exp \left[{\displaystyle \frac{\mathrm{3}\ln \left(\mathrm{10}\right)}{\mathrm{2}}}\left(M-\mathrm{4}\right)\right],\\ \text{small }M\\ \mathrm{2}F{w}_{\mathrm{0}}\exp \left[\ln \left(\mathrm{10}\right)\left(M-\mathrm{4}\right)\right]+\mathit{\pi }{F}^{\mathrm{2}}{w}_{\mathrm{0}}^{\mathrm{3}},\\ \text{large }M\end{array}\right.\end{array}$$

which is a closed-form expression of the same form as the original Utsu productivity law (Eq. 1). Note that K and δ₊ are both, implicitly, functions of the selected minimum aftershock magnitude threshold m₀.

Here, we predict that the α value decreases from $\mathrm{3}\ln \left(\mathrm{10}\right)/\mathrm{2}\approx \mathrm{3.45}$ to ln (10)≈2.30 when switching regime from small to large mainshocks (or from 1.5 to 1 in a 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. Hainzl et al. (2010) simulated aftershocks using the clock-change model where events were advanced in time by the static stress change produced by a mainshock in a 3-D medium. They explained the scaling break observed in simulation as a transition from 3-D to 2-D scaling regime when the mainshock rupture dimension approached w₀, which is compatible with the present demonstration. 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).

2.2 Testing of the SSP on the aftershock spatial distribution

The SSP predicts a step-like behavior of the aftershock spatial density for an idealized smooth static stress field (Fig. 4a–b), which is in disagreement with real aftershock observations. A number of studies have shown that the spatial linear density distribution of aftershocks ρ is well represented by a power law, expressed as

$$\begin{array}{}\text{(13)}& \mathit{\rho }\left(r\right)\propto r^{-q},\end{array}$$

with r the distance from the mainshock and q the power-law exponent. This parameter ranges over $\mathrm{1.3}\le q\le \mathrm{2.5}$ (Felzer and Brodsky, 2006; Lipiello et al., 2009; Marsan and Lengliné, 2010; Richards-Dinger et al., 2010; Shearer, 2012; Gu et al., 2013; Moradpour et al., 2014; van der Elst and Shaw, 2015). Although Felzer and Brodsky (2006) suggested a dynamic stress origin for aftershocks, their results were later questioned by Richards-Dinger et al. (2010). Most of the studies cited above suggest that the q value is explained from a static stress process. As for the examples of aftershocks shown to be dynamically triggered (e.g., Fan and Shearer, 2016), they are too few to alter the aftershock productivity law and too remote to be consistently defined as aftershocks in cluster methods.

Figure 4 Spatial distribution of aftershocks following the SSP: (a) Smooth static stress field as a function of distance r from the mainshock, with $\mathrm{\Delta }{\mathit{\sigma }}_{\mathrm{0}}=-\mathrm{10}$ bar and c=10 km (Eq. 6); (b) step-like aftershock spatial linear density ρ(r) with ${\mathit{\delta }}_{+}=\mathrm{1000}$ event km⁻¹, δ₀=1 event km⁻¹ and $\mathrm{\Delta }{\mathit{\sigma }}_{\ast }=-\mathrm{0.3}\mathrm{\Delta }{\mathit{\sigma }}_{\mathrm{0}}$ (ad hoc ratio yielding $r_{\ast }=\mathrm{3.5}$ km; Eq. (7) – event distances sampled from the δ(r) distribution, repeated 100 times). Such distribution is not observed in nature; (c) same as panel (a) but with random uniform noise representative of spatial heterogeneities added to the regional stress field; (d) power-law-like aftershock spatial linear density ρ(r) with power exponent MLE estimate q=1.96, representative of real aftershock observations (see Fig. 5a), due to the addition of uniform noise to the static stress field.

In a more realistic setting, the static stress field must be heterogeneous (due to the occurrence of previous events and other potential stress perturbations). We therefore simulate the static stress field by adding a uniform random component bounded over $\pm \mathrm{\Delta }{o}_{\ast }$ following Mignan (2011) (see also King and Bowman, 2003). Note that any deviation above Δo_∗ would be flattened to Δo_∗ over time by temporal diffusion (the so-called “historical ghost static stress field” in Mignan, 2016a). Figure 4c shows the resulting stress field and Fig. 4d the predicted aftershock spatial density. Adding uniform noise blurs the contour of the aftershock solid, switching the aftershock spatial density from a step function (Fig. 4b) to a power law (Fig. 4d). We fit Eq. (13) to the simulated data using the maximum likelihood estimation (MLE) method with $r_{min}=r_{\ast }$ (Clauset et al., 2009) and find $q=\mathrm{1.96}\pm \mathrm{0.01}$, in agreement with the aftershock literature. This result alone is, however, insufficient to prove the validity of the SSP.

3 Observations and model fitting

3.1 Data

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, considering only the first aftershock generation). Only events with magnitudes greater than m₀=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).

3.2 Aftershock spatial density distribution

Figure 5a represents the spatial linear density distribution of aftershocks ρ(r) for the four largest strike-slip mainshocks in southern California: 1987 M=6.6 Superstition Hills, 1992 M=7.3 Landers, 1999 M=7.1 Hector Mine and 2010 M=7.2 El Mayor. The distance between mainshock and aftershocks is calculated as $r=\sqrt{(x-x_{\mathrm{0}}{)}^{\mathrm{2}}+(y-y_{\mathrm{0}}{)}^{\mathrm{2}}}$, with (x, y) the aftershock coordinates and (x₀, y₀) the coordinates of the nearest point to the mainshock fault rupture (as depicted in Fig. 2c). The dashed black lines shown in Fig. 5a are visual guides to q=1.96, showing that the SSP is compatible with real aftershock observations.

Figure 5 Estimating the solid seismicity parameters from the spatial distribution of aftershocks: (a) spatial linear density distribution ρ(r) of aftershocks for the four largest strike-slip mainshocks in southern California (with first-generation aftershocks only; the density distribution comprising all aftershocks generated by the Landers mainshock is represented by the dotted curve to illustrate the type of spatial heterogeneity, such as the Big Bear cluster, not considered in the present study – see also Fig. 2c). The solid seismicity parameters $r_{\ast }=\mathrm{1}$ km and ${\mathit{\delta }}_{+}(m_{\mathrm{0}}=\mathrm{2})=\mathrm{1.23}$ event km⁻³ can be retrieved from the observed plateau ρ (r < r_∗), in agreement with the SSP (see Fig. 4d). Note that the spatial power-law decay at high r is similar to the one expected by the SSP in the case of a static stress field with additive uniform noise (expected q=1.96 represented by the dashed black lines); (b) aftershock productivity K for M > 6. The curves represent the productivity law as defined by solid seismicity (Eq. 17) for different w₀ values (first term only corresponds to w₀=0; Eq. 18).

Comparing Fig. 5a to Fig. 4d suggests that r_∗ can be roughly estimated from the spatial linear density plot, being the maximum distance r at which the plateau ends, here leading to $r_{\ast }\approx$ 1 km. This parameter is constant for different large M values since both w₀ and Δσ₀ are constant while Δσ_∗ is also a priori a constant. We can then estimate the ratio $\mathrm{\Delta }{\mathit{\sigma }}_{\ast }/\mathrm{\Delta }{\mathit{\sigma }}_{\mathrm{0}}$ from Eq. (7). However, the result is ambiguous due to uncertainties in the width w₀. For $w_{\mathrm{0}}=\left\{\mathrm{5},\mathrm{10},\mathrm{15}\right\}$ km, we get $\mathrm{\Delta }{\mathit{\sigma }}_{\ast }/\mathrm{\Delta }{\mathit{\sigma }}_{\mathrm{0}}=\left\{-\mathrm{0.54},-\mathrm{1.01},-\mathrm{1.38}\right\}$.

As for the plateau value ρ (r < r_∗), it provides an estimate of the aftershock activation density δ₊, with

$$\begin{array}{}\text{(14)}& {\mathit{\delta }}_{+}={\displaystyle \frac{\mathit{\rho }\left(M,r<r_{\ast }\right)}{\exp \left[\ln \left(\mathrm{10}\right)(M-\mathrm{4})\right]}}\end{array}$$

a volumetric density, i.e., the linear density ρ normalized by the mainshock rupture area (Eq. 4). Due to the fluctuations in ρ (r < r_∗), δ₊ will be estimated from the productivity law instead (see Sect. 3.3), and ρ (r < r_∗) will then be estimated from Eq. (14) (horizontal dashed colored lines), as detailed below.

It should be noted that we consider only the first-generation aftershocks to avoid ρ heterogeneities from secondary aftershock clusters occurring off-fault. An example of such heterogeneity and anisotropy is illustrated by the Landers–Big Bear case (Fig. 2c; dotted colored curve in Fig. 5a). Those cases are not systematic and therefore not considered in the aftershock productivity law. However, they are also due to static stress changes (e.g., King et al., 1994) with the anisotropic effects explainable by solid seismicity through the concept of historical ghost static stress field (Mignan, 2016a).

3.3 Aftershock productivity law

The observed number n of aftershocks of magnitude m≥m₀ produced by a mainshock of magnitude M (for a total of N mainshocks) in southern California is shown in Figs. 5b (for large M≥6) and 6a (for the full range M≥m₀). We fit Eq. (1) to the data using the MLE method with the log-likelihood function

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\mathrm{LL}\left(\mathit{\theta };X=\left\{n_i;i=\mathrm{1},\mathrm{\dots },N\right\}\right)=\\ \text{(15)}& {\displaystyle }& {\displaystyle }{\sum }_{i=\mathrm{1}}^N\left[n_i\ln \left[K_i\left(\mathit{\theta }\right)\right]-K_i\left(\mathit{\theta }\right)-\ln \left(n_i\mathrm{!}\right)\right]\end{array}$$

for a Poisson process, representing the stochasticity of the count K of aftershocks produced by a mainshock at any given time. Inserting Eq. (1) in Eq. (15) yields

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\mathrm{LL}\left(\mathit{\theta }=\left\{K_{\mathrm{0}},\mathit{\alpha }\right\};X\right)=\ln \left(K_{\mathrm{0}}\right){\sum }_{i=\mathrm{1}}^N{n}_i\\ {\displaystyle }& {\displaystyle }+\mathit{\alpha }{\sum }_{i=\mathrm{1}}^N\left[n_i\left(M_i-m_{\mathrm{0}}\right)\right]-K_{\mathrm{0}}{\sum }_{i=\mathrm{1}}^N\exp \left[\mathit{\alpha }\left(M_i-m_{\mathrm{0}}\right)\right]\\ \text{(16)}& {\displaystyle }& {\displaystyle }-{\sum }_{i=\mathrm{1}}^N\ln \left(n_i\mathrm{!}\right)\end{array}$$

(note that the last term can be set to 0 during LL maximization). For southern California, we obtain α_MLE=2.32 (1.01 in log₁₀ scale) and K₀=0.025 when considering large (M≥6) mainshocks only to avoid the issues of scaling break and data dispersion at lower magnitudes. This result, represented by the black solid line in Fig. 5b, is in agreement with previous studies in the same region (e.g., Helmstetter, 2003; Helmstetter et al., 2005; Zaliapin and Ben-Zion, 2013; Seif et al., 2017) and with $\mathit{\alpha }=\ln \left(\mathrm{10}\right)\approx \mathrm{2.30}$ predicted for large mainshocks in solid seismicity (Eq. 12). Moreover we find a bulk β_MLE=2.34 (1.02 in log₁₀ scale) (Aki, 1965), in agreement with α≡β.

Figure 6 Aftershock productivity defined as the number of aftershocks K(m₀=2) per mainshock of magnitude M: (a) observed aftershock productivity in southern California with aftershocks selected using the nearest-neighbor method; (b) seismicity time series with distinction made between background events and aftershocks, observed (“obs”, in black) and ETAS-simulated (“sim”, colored); (c) true simulated aftershock productivity with kink, defined from Eq. (20); (d) retrieved simulated aftershock productivity with aftershocks selected using the nearest-neighbor method – data points in panels (a), (c) and (d) are represented by grey dots; the model MLE fits are represented by the dashed and dotted black lines for M ≥ 6 and M≥m₀, respectively; dashed and dotted grey lines are visual guides to $\mathit{\alpha }=\mathrm{3}/\mathrm{2}\ln \left(\mathrm{10}\right)$ and ln(10), respectively.

Let us now rewrite the solid seismicity aftershock productivity law (Eq. 12) by only considering the large M case and injecting $r_{\ast }=F{w}_{\mathrm{0}}$ (by combining Eqs. 7–8). We get

$$\begin{array}{}\text{(17)}& K\left(M>M_{\mathrm{break}}\right)={\mathit{\delta }}_{+}\left\{\mathrm{2}{r}_{\ast }\exp \left[\ln \left(\mathrm{10}\right)(M-\mathrm{4})\right]+\mathit{\pi }{r}_{\ast }^{\mathrm{2}}{w}_{\mathrm{0}}\right\}.\end{array}$$

The role of w₀ is illustrated in Fig. 5b for different values (dashed and dotted curves) and shown to be insignificant for large M values. Therefore Eq. (17) can be approximated to

$$\begin{array}{}\text{(18)}& K\left(M>M_{\mathrm{break}}\right)\approx \mathrm{2}{\mathit{\delta }}_{+}{r}_{\ast }\exp \left[ln\left(\mathrm{10}\right)(M-\mathrm{4})\right].\end{array}$$

By analogy with Eq. (1), we get

$$\begin{array}{}\text{(19)}& {\mathit{\delta }}_{+}={\displaystyle \frac{K_{\mathrm{0}}\exp \left[\ln \left(\mathrm{10}\right)\left(\mathrm{4}-m_{\mathrm{0}}\right)\right]}{\mathrm{2}{r}_{\ast }}}.\end{array}$$

With $r_{\ast }\approx \mathrm{1}$ km estimated from ρ(r) (Sect. 3.2) and K₀=0.025, we obtain ${\mathit{\delta }}_{+}=\mathrm{1.23}$ events km⁻³ for m₀=2. We then get back the plateau ρ (r < r_∗) for different M values from Eq. (14), as shown in Fig. 5a (horizontal dashed colored lines). Although based on limited data, this result suggests that the activation parameter δ₊ is constant (at least for large M) in southern California. Note that if ρ (r < r_∗) was well constrained, it could have been estimated jointly with r_∗ from Fig. 5a to predict the aftershock productivity law of Fig. 5b without further fitting required (hence removing K₀ from the equation, K₀ having no physical meaning in solid seismicity).

4 Role of aftershock selection on productivity scaling break

We tested the following piecewise model to identify any break in scaling at smaller M, as predicted by Eq. (12):

$$\begin{array}{ll}\text{(20)}& {\displaystyle }& {\displaystyle }K\left(M\right)={\displaystyle }& {\displaystyle }\left\{\begin{array}{l}{K}_{\mathrm{0}}{\displaystyle \frac{\exp \left[\ln \left(\mathrm{10}\right)(M_{\mathrm{break}}-m_{\mathrm{0}})\right]}{\exp \left[\frac{\mathrm{3}}{\mathrm{2}}\ln \left(\mathrm{10}\right)(M_{\mathrm{break}}-m_{\mathrm{0}})\right]}}\exp \left[{\displaystyle \frac{\mathrm{3}}{\mathrm{2}}}\ln \left(\mathrm{10}\right)(M-m_{\mathrm{0}})\right],\\ M\le M_{\mathrm{break}}\\ K_{\mathrm{0}}\exp \left[\ln \left(\mathrm{10}\right)(M-m_{\mathrm{0}})\right],\\ M>M_{\mathrm{break}}\end{array}\right.\end{array}$$

but with the best MLE result obtained for M_break=m₀, suggesting no break in scaling in the aftershock productivity data, as observed in Fig. 6a. Final parameter estimates are α_MLE=1.95 (0.85 in log₁₀ scale) and K₀=0.141 for the full mainshock magnitude range M≥m₀ (dotted line), subject to high scattering at low M values.

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

$$\begin{array}{}\text{(21)}& \left\{\begin{array}{l}\mathit{\lambda }\left(t,x,y\right)=\mathit{\mu }\left(t,x,y\right)+{\sum }_{i:t_j<t}K\left(M_i\right)\\ f(t-t_i)g\left(x-x_i,y-y_i\left|M_i\right.\right)\\ f\left(t\right)=c^{p-\mathrm{1}}(p-\mathrm{1})(t+c{)}^{-p}\\ g\left(x,y\left|M\right.\right)={\displaystyle \frac{\mathrm{1}}{\mathit{\pi }}}{\left(d{e}^{\mathit{\gamma }\left(M-m_{\mathrm{0}}\right)}\right)}^{q-\mathrm{1}}\\ {\left(x^{\mathrm{2}}+y^{\mathrm{2}}+d{e}^{\mathit{\gamma }\left(M-m_{\mathrm{0}}\right)}\right)}^{-q}(q-\mathrm{1}).\end{array}\right.\end{array}$$

Aftershock sequences are defined by power laws, both in time and space (for an alternative temporal function, see Mignan, 2015, 2016b; the spatial power-law distribution is in agreement with solid seismicity in the case of a heterogeneous static stress field – see Sect. 2.2). The southern California background seismicity, as defined by the nearest-neighbor method (with same t, x, y and m), is denoted by μ. We fix the ETAS parameters to $\mathit{\theta }=\left\{c=\mathrm{0.011}\right.$ day, p=1.08, d=0.0019 km², q=1.47, γ=2.01, β=2.29, K₀=0.08}, following the fitting results of Seif et al. (2017) for the southern California relocated catalog and m₀=2 (see their Table 1). However, we define the productivity function K(M) from Eq. (20) with M_break=5. Examples of ETAS simulations are shown in Fig. 6b for comparison with the observed southern California time series. Figure 6c allows us to verify that the simulated aftershock productivity is kinked at M_break, as defined by Eq. (20).

We then select aftershocks from the ETAS simulations with the nearest-neighbor method. Figure 4d represents the estimated aftershock productivity, which has lost the break in scaling originally implemented in the simulations (with an underestimated α_MLE=2.07 as observed in the real case for M≥m₀). Note that a similar result is obtained when using a windowing method (Gardner and Knopoff, 1974). 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 a bias is possible, it 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.

One other possible explanation for the lack of scaling break is that our demonstration assumes moment magnitudes while the southern California catalogue is in local magnitudes. Deichmann (2017) demonstrated that while M_L∝M_w at large M, M_L∝1.5M_w at smaller M values. This could in theory cancel the kink in real data. However, the scaling break predicted by Deichmann (2017) occurs at several magnitude units below the geometric scaling break expected by solid seismicity, invalidating this second option for mid-range magnitudes M.

5 Conclusions

In the present study, a closed-form expression defined from geometric and static stress parameters was proposed (Eq. 12) to describe the empirical Utsu aftershock productivity law (Eq. 1). This demonstration is similar to the previous ones made by the author to explain precursory accelerating seismicity and induced seismicity (Mignan, 2012, 2016a). In all these demonstrations, the main physical parameters remain the same, i.e., the activation density δ₊ (also δ₋ and δ₀), the background stress amplitude range Δo_∗ and the solid envelope r_∗ which describes the geometry of the “seismicity solid” (Fig. 3a–b). Further studies will be needed to evaluate whether the δ₊ and Δo_∗ parameters are universal or region-specific and if the same values apply to different types of seismicity at a same location.

Although the solid seismicity postulate (Eq. 5) remains to be proven, it is so far a rather convenient and pragmatic assumption to make to determine the physical parameters that play a first-order role in the behavior of seismicity. The similarity of the SSP-simulated and observed values of the power-law exponent q of the aftershock spatial density distribution shows that the SSP is consistent with large aftershock observations once uniform noise is added to the stress field (Figs. 4d–5a). The impact of other types of noise on q has yet to be investigated. The SSP 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 to better understand the behavior of seismicity. Its current limitation in the case of aftershock analysis consists of assuming that the static stress field is radial and described by Eq. (6) (e.g., 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. 6a), such as overpressure due to trapped high-pressure gas for example (Miller et al., 2004 – see also Mignan, 2016a, for an overpressure field due to fluid injection). 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).

Finally, the disappearance of the predicted scaling break in the aftershock productivity law once declustering is applied (Fig. 6) indicates that more work is required in that domain. Only a declustering technique that does not dictate a constant scaling at all M will be able to identify whether a scaling break really exists or not.