It has recently been found that the number of aftershocks of large earthquakes in the world and the number of direct aftershocks of earthquakes in different regions of the world, considered in a fixed magnitude range relative to the main shock, obeys an exponential distribution (Shebalin et al., 2020a; Shebalin et al., 2020b). This law, called earthquake productivity law, was established for different magnitudes of the main shocks, different ways of identifying direct aftershocks, in a wide range of parameters of the algorithm for identifying aftershocks (Shebalin et al., 2020b). The data used made it possible to consider ranges of aftershock magnitudes with a lower threshold within the range of up to 2.5 units of magnitude less than the magnitude of the main shock. A further increase in the range in the reviewed catalogs was impossible, since the entire range of considered magnitudes should be above the general completeness threshold. Therefore, an increase in the range is possible only by increasing the magnitude threshold for the main shocks, which leads to their number being too small.

The exponential shape of the distribution means that the most probable number of direct aftershocks (the mode of the distribution) is 0. This property of earthquake productivity seems so counterintuitive that testing whether this shape of the distribution persists as the range of aftershock magnitudes increases has an important independent meaning. The importance of this verification is also reinforced by the fact that the exponential form of productivity contradicts one of the key elements of the ETAS stochastic model (Ogata, 1998) widely used for modeling seismicity. It is assumed in ETAS model that the productivity of earthquakes is a function of their magnitude. Under this assumption, productivity should have the form of a Poisson distribution with a non-zero mode. The same assumption is used in stochastic methods for declustering earthquake catalogs, in particular, the method based on the ETAS model (Zhuang et al., 2002) and in the model-independent MISD method (Marsan and Lengline, 2008).

The aim of this work is to investigate whether the exponential form of the distribution of earthquake productivity will be preserved with a significant expansion of the range of aftershock magnitudes. This can only be done in a region with a dense network of seismic stations and, at the same time, a high frequency of large earthquakes. We chose the land part of Japan, where the representative magnitude for crustal earthquakes since 2000 is about 1.0, and earthquakes of magnitude 6 and above occur several times a year.

We follow the definition of earthquake productivity adopted by Shebalin et al. (2020b). In the earthquake flow, each event is considered as a potential “parent” of subsequent earthquakes, and vice versa, each event can have a “parent”, in which case it is considered an “offspring”. The offsprings can be interpreted as immediate aftershocks. We use the nearest neighbor scheme by Zaliapin and Ben-Zion (2013). In this scheme, each parent (triggering event) can have multiple offsprings (triggered events), but each offspring can only have one parent. The parent of a given event is determined by the minimum of the proximity function (“nearest neighbor”). The proximity function determines the extent to which “parent” and “offspring” can be considered related. In this paper, following Zaliapin and Ben-Zion (2013), we use the proximity function proposed by Baiesi and Paczuski (2004): η i j = t i j r i j d f 1 0 − b m i , t i j > 0 , + ∞ , t i j ≤ 0 , (1) where t _ij = t _j − t _i is the interevent time, r _ij the spatial distance between the epicenters, m _i the magnitude of event i, d _f the fractal dimension of the epicenter distribution and b the slope of the earthquake-size distribution.

Zaliapin and Ben-Zion (2013) have shown that the distribution of the nearest neighbor proximity function usually has a bimodal distribution, in which small values correspond to related events (aftershocks, foreshocks, swarms), and large values correspond to independent events. To separate related and independent events, a threshold value for the proximity function is introduced. Zaliapin and Ben-Zion (2013) proposed to approximate this bimodal distribution by the sum of two log-normal distributions. However, in reality, the form of distribution can be more complex. Therefore, Shebalin et al. (2020b) proposed to approximate the right side of the distribution using a randomized catalog of earthquakes. We believe that the distribution of the nearest neighbor proximity function in the randomized catalog is the same as for independent events. The distribution of related events is then defined as the difference between the overall distribution and the distribution for the randomized catalog. The threshold is defined in such a way as to equalize the error rates (related events with a proximity function above the threshold and independent events with a proximity function below the threshold).

After introducing the threshold, it turns out that some of the events do not have a “parent” because the proximity to the nearest neighbor exceeds the threshold. We call such events “background”. If the proximity is below the threshold the event has a “parent”. We call such events offsprings. The productivity of an earthquake then is the number of its offsprings. Here, as in (Shebalin et al., 2020b), we will use ΔM-productivity. We count the number of offsprings of magnitude M _a ≥ M _m − ΔM, where M _m is the magnitude of the parent and M _a is the magnitude of the offspring. Note that the magnitude of the offspring may be greater than the magnitude of the parent.

The scheme used here differs from the traditional mainshock-aftershock outline, which usually assumes that aftershocks are weaker than the mainshock. Another difference is that an aftershock sequence consists of a hierarchical tree of parent-offspring sequences with several levels of hierarchy (offspring of offspring, etc.). In the traditional scheme (Zaliapin and Ben-Zion, 2013), the total number of aftershocks is calculated. We count the number of offsprings at one level of the hierarchy, which can be interpreted as direct aftershocks. In this scheme, each earthquake is characterized by its ΔM-productivity value. It makes sense to analyze the values of ΔM-productivity for earthquakes with magnitude M ≥ M _c + ΔM. In this case, the offsprings are counted above the magnitude of complete registration M _c .

For the territory under consideration, estimates were made of the values of the proximity function parameters (1): b = 0.86, d _f = 1.68 and log₁₀ η ₀ = −1.46. b-value is determined by Bender (1983) method in magnitude interval m ≥ 3; d _f is determined by Grassberger (1983) method also for events of magnitude m ≥ 3. The η ₀ threshold was determined using the method from (Shebalin et al., 2020b). To avoid the possible influence of the Tohoku earthquake on 11 March 2011 with M = 9.1, all three parameters were estimated using data for 2000–2010.

For each earthquake with M _m ≥ 6.0, we calculated the productivity: the number of offsprings with magnitude M _a ≥ M _m − ΔM, ΔM = 5. It turned out that 39,464 triggered events were associated with 56 parent events. The frequency-productivity graph is shown in Figure 2. This plot is similar to the commonly used frequency-magnitude cumulative graph, in which frequencies are summed starting at higher values. The rectilinear form of the graph, as in the case of the Gutenberg-Richter law, indicates the exponential form of the distribution. The difference is that the productivity of each earthquake is an integer, so the resulting distribution is more correctly interpreted as a geometric distribution, for which the cumulative frequency plot, starting from large values of the argument, also has a linear form on a logarithmic scale. The slope of the graph is uniquely related to a single distribution parameter.

Thus, we may conclude that the exponential form of the distribution of ΔM-earthquake productivity established by Shebalin et al. (2020b) for ΔM ≤ 2.5 is also confirmed for a very large range of magnitudes for ΔM = 5.

In the presented analysis, the earthquake productivity was not separated according to different levels of the hierarchy due to small number of M _m ≥ 6.0 events. Using worlwide stastistics of productivity, it was shown by Shebalin et al. (2020b) that the exponential form of the distribution and the value of its parameter change little for different levels of the hierarchy. It is this observation that gave grounds to simultaneously analyze the productivity of all earthquakes, regardless of whether they are main shocks or aftershocks in the traditional terminology. We repeated this check for earthquakes of the land part of Japan with M _m ≥ 5.0 and ΔM = 4. This property is maintained: for at least eight levels of the hierarchy, starting with background events (hierarchy level 0), the frequency-productivity plots remain linear on a logarithmic scale and have approximately the same slope (Figure 3).

We verified whether the observed exponential distribution of productivity is a property of the data, or, alternatively, the result of the choice of the proximity function. We selected only background earthquakes (events of the hierarchy level 0) from the catalog and, assuming η ₀ = ∞, repeated the procedure for finding for each event its parent (the closest neighbor) and calculated the productivity. The resulting distribution of the productivity of background earthquakes has a pronounced non-zero maximum (Figure 4) and, thus, is not geometric. Comparison of the distributions for background and clustered seismicity demonstrates completely different spatiotemporal structure. Thus, it is confirmed that the exponential form of the productivity distribution is a property of clustered seismicity.

The main result obtained here—the exponential form of the distribution of the ΔM-productivity of earthquakes [the productivity law (Shebalin et al., 2020b)] is preserved even at a very large value of ΔM. The density of the geometric distribution—an integer version of the exponential one—has a maximum value at 0. This makes this form of distribution counterintuitive. It is difficult to imagine that the most probable number of direct aftershocks is 0. It is even more difficult to believe this when direct aftershocks are considered, the magnitude of which is 5 units less than the magnitude of the earthquakes that caused them. Our analysis, however, shows that this is possible, and it becomes clear why the productivity law can remain valid even for large ΔM. The point is that for large magnitude ranges under consideration, the total number of offsprings is orders of magnitude greater than the number of parent earthquakes. In the considered example (Figure 2), with M _m ≥ 6 and Δ _M = 5, the productivity reaches 3,000, but the total number of parent events is only 56. Under these conditions, the probability of realizing the productivity value exactly 0 is small: p = 0.0014, and expected number of such earthquakes is 0.08, e.g., zero, because number of earthquakes is integer. Nevertheless, the smallest productivity values still predominate. If the number of analyzed events is much greater than the average productivity, then the number of events without direct aftershocks (productivity 0) is indeed large compared to the number of events of any other productivity (see Supplementary Material).

It is usually assumed that the productivity of earthquakes depends mainly on their magnitude. This, in particular, is used in the ETAS (Zhuang et al., 2002) and MISD (Marsan and Lengline, 2008) stochastic declustering methods. The expected ΔM-productivity distribution for such models is a Poisson distribution having a pronounced non-zero mode, this is actually predefined in these models. However, Shebalin et al. (2020b) showed, and here it is confirmed for large values of ΔM, that in reality this distribution is rather a monotone geometric one with a mode at 0. This gives additional advantages to the Zaliapin and Ben-Zion (2013) declustering method, which does not use any specified assumption. The spatiotemporal version of the ETAS model (Ogata, 1998; Zhuang et al., 2002) is also often used for seismicity modeling. The Poisson productivity distribution embedded in the model contradicts the observations. The earthquake productivity law confirmed here makes it possible to correct the ETAS model.

In a nearest neighbor schema, each parent can have multiple offsprings, but each offspring can only have one parent. This makes the productivity averaging procedure meaningful, since in such a procedure each offspring is taken into account only once. We have also shown that productivity depends little on the level of the hierarchy. This means that there is no need to distinguish between main shocks and aftershocks during averaging. ΔM-productivity averaged over some space-time region, we call after Shebalin et al. (2020b) the clustering factor. Note that in ETAS model the clustering factor is called a branching ratio. It should be less than 1 to avoid diverging aftershocks occurrence rate (Helmstetter and Sornette, 2002). For large ΔM we obtain clustering factor much larger than 1. But there is no apparent contradiction here, because the ETAS model assumes a strong relationship between the number of direct aftershocks and the main shock magnitude (Helmstetter, 2003), which was refuted by Shebalin et al. (2020b) and here again. In ΔM-analysis, the number of offsprings turns out to be weakly dependent on the magnitude of the parent, and its statistical distribution has exponential form with clear maximum at 0. This explains why large values of the clustering factor do not cause a diverging aftershock sequence. Of course, if we consider all offsprings with magnitudes above a certain threshold, their number increases with the magnitude of the parent. But this is controlled by Gutenberg-Richter law.

ΔM-productivity, similar manner to magnitude, can be seen as a property inherent in every earthquake. In this case, this value does not have to be an integer. Let’s denote this value λ. The observed value in this case is the concrete realization of the “potential” productivity of λ. It can be assumed that a particular sample is a Poisson random variable with rate λ. One more assumption can be made: the parameter λ of earthquakes, similarly to the magnitude, has an exponential distribution of the form: p e λ = 1 Λ Δ M exp − λ Λ Δ M , (2) where Λ_ΔM is a parameter.

It is easy to show (Shebalin et al., 2020b) that the distribution of a sample of the productivity of an arbitrary earthquake in this case has the form of a geometric distribution. The estimate of the parameter Λ_ΔM is the average value of λ over the ensemble. Thus, the clustering factor has the meaning of the Λ_ΔM parameter and its estimate. The clustering factor is a convenient and a simple parameter to characterize the productivity in a set of earthquakes, for example, earthquakes in a certain space-time volume.