# Earthquake productivity law

P. N. Shebalin, C. Narteau, S. V. Baranov

► To cite this version:

P. N. Shebalin, C. Narteau, S. V. Baranov. Earthquake productivity law. Geophysical Journal International, 2020, 222, pp.1264-1269. 10.1093/gji/ggaa252. insu-03748811

HAL Id: insu-03748811

https://insu.hal.science/insu-03748811

Submitted on 10 Aug 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

EXPRESS LETTER

Earthquake productivity law

P.N. Shebalin<sup>1</sup>, C. Narteau<sup>2</sup> and S.V. Baranov<sup>1,3</sup>

Accepted 2020 May 18. Received 2020 May 17; in original form 2020 March 11

SUMMARY

Mechanisms of stress transfer and probabilistic models have been widely investigated to explain earthquake clustering features. However, these approaches are still far from being able to link individual events and to determine the number of earthquakes caused by a single event. An alternative approach based on proximity functions allows us to generate hierarchical clustering trees and to identify pairs of nearest-neighbours between consecutive levels of hierarchy. Then, the productivity of an earthquake is the number of events of the next level to which it is linked. Using a relative magnitude threshold ΔM to account for scale invariance in the triggering process, we show that the ΔM-productivity attached to each event is a random variable that follows an exponential distribution. The exponential rate of this distribution does not depend on the magnitude of triggering events and systematically decreases with depth. These results could now be used to characterize active fault systems and improve epidemic models of seismicity.

**Key words:** Probabilistic forecasting; Earthquake interaction, forecasting, and prediction; Seismicity and tectonics; Statistical seismology.

# 1 INTRODUCTION

The earthquake productivity is a key parameter in statistical seismology as it gives the number of events within a space–time–magnitude window of observation. In earthquake-size distributions, the productivity determines the level of seismic activity independently of the ratio between the number of large and small events. However, an important feature of seismicity is the occurrence of space–time clusters, which unequivocally demonstrate that earthquakes interact with each other. Hence, by focusing on the way in which a sequence of earthquakes develops over space and time, the productivity can also be regarded as the number of events resulting from the perturbation of the state of stress induced by another earthquake. Such a productivity has first been used to develop appropriate models of aftershock occurrence considering the empirical Omori-Utsu law (Utsu 1969, 1970). The number of aftershocks in a given time interval can then be computed according to the magnitude *m* of the mainshock and the properties of the power-law aftershock decay rate (Utsu *et al.* 1995; Shcherbakov *et al.* 2004; Shcherbakov & Turcotte 2004; Holschneider *et al.* 2012; Davidsen *et al.* 2015).

Since the advent of epidemic models of seismicity (i.e. ETAS models), the productivity has become a major issue because it is the main ingredient that determines the increase in seismicity rate

after each earthquake (Kagan & Knopoff 1981; Ogata 1989; Helmstetter & Sornette 2002). In all these models, the number of events triggered by a magnitude *m* earthquake is considered to vary as a Poisson process of rate $\langle N(m) \rangle \times 3009 \times 232A; = K10^{am}$. The reported α-values vary from 0.5 to 2 (Console *et al.* 2003; Helmstetter 2003; Zhuang *et al.* 2004, 2005; Hainzl & Marsan 2008; Werner & Sornette 2008; Wang *et al.* 2010; Hainzl *et al.* 2013), but they are often close to the observed *b*-value, the slope of the earthquake-size distribution. Nevertheless, these estimates remain uncertain due to the difficulty of isolating the relative contributions of successive events in a sequence.

Despite the diversity of declustering methods implemented in the past, the research dealing with causal links within cascades of triggered seismicity is still in its early stage and there is no definitive classification as yet. A first approach is to separate the branching structure of earthquake sequences from the background rate using an iterative algorithm based on maximum likelihood estimation of the parameters of an epidemic model of seismicity (Zhuang *et al.* 2002). Without an *a priori* model, a second approach is to identify, directly and indirectly, triggered events assuming linear contributions of each earthquake to the overall seismicity rate (Marsan & Lengliné 2008). Another approach is dedicated to the identification of earthquake clusters using proximity functions in time–space–magnitude domains (Baiesi & Paczuski 2004; Zaliapin *et al.* 2008). All these methods confirm the dependency of the productivity on

<footnote n="1">¹Institute of Earthquake Prediction Theory and Mathematical Geophysics, 84/32 Profsouznaya, Moscow 117997, Russia</footnote>

<footnote n="2">²Université de Paris, Institut de physique du globe de Paris, CNRS, F-75005 Paris, France. E-mail: narteau@ipgp.fr</footnote>

<footnote n="3">³Kola Branch, Federal Research Center Geophysical Survey, Russian Academy of Sciences, 14 Fersmana, Apatity 184209, Russia</footnote>
the magnitude *m* of the triggering event. However, less attention has been given to the overall variability in the number *N* of triggering events in the seismic catalogues (Marsan & Helmstetter 2017). Here we examine this specific issue using a proximity function and clusters of nearest-neighbours earthquakes.

## 2 METHODS

For each pair of earthquakes {i, j}, we compute the proximity function (Baiesi & Paczuski 2004),

$$ \eta_{ij} = \begin{cases} t_{ij}(r_{ij})^{dt} 10^{-bm_i} & \text{for } t_{ij} > 0, \\ +\infty & \text{for } t_{ij} \le 0, \end{cases} \quad (1) $$

where $t_{ij} = t_j - t_i$ is the interevent time, $r_{ij}$ the spatial distance between the epicentres, $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. Hierarchical clustering trees of nearest-neighbour events are then constructed using a threshold $\eta_0$ for the proximity function (Zaliapin & Ben-Zion 2013, 2016). Here, this threshold value is computed using a non-parametric method, without fitting a model to the distribution of $\eta$-values. Instead, we generate a random version of the same catalogue by shuffling the earthquake times. Assuming that the distribution of $\eta$-values of this random catalogue provides the non-clustered reference distribution for the real catalogue, we select the $\eta_0$-value for which the two types of errors compensate each other: same probability of having causally related events with $\eta > \eta_0$ and independent events with $\eta < \eta_0$ (Note 1 in the Supporting Information).

Here all earthquakes are considered as triggering events (parents). Using the threshold $\eta_0$-value between two consecutive levels of hierarchy, each earthquake is also considered as a triggered event (child) if it is linked to the higher level of hierarchy. Considering only the first nearest-neighbour in the time-space-magnitude domain (i.e. the smallest $\eta$-value), a triggered event is associated with only one triggering event. According to these rules, isolated earthquakes with no parent and no child can occur. All hierarchical clustering trees are built from a primary triggering event, which has no parent. From any earthquake, there is a single path to a unique primary event. Level 0 of the hierarchical clustering trees is only composed of primary triggering events (no parent). Level 1 is composed of earthquakes triggered by events from level 0 and so on. Interestingly, a triggered event may have a larger magnitude than its triggering events. Thus, we study seismicity patterns beyond the traditional mainshock-aftershock sequences.

For each triggering earthquake, we count the number of triggered events at the lower hierarchical level using a relative magnitude threshold $\Delta M$ to account for scale invariance in the triggering process (i.e. $M_{triggering} - M_{triggered} < \Delta M$). This number of triggered events is defined as the $\Delta M$-productivity (hereafter the productivity). The distribution of the number of triggered events for an earthquake population is defined as the productivity distribution with a mean denoted $\Lambda_{\Delta M}$ and called the clustering factor.

## 3 EARTHQUAKE PRODUCTIVITY IN WORLDWIDE AND REGIONAL CATALOGUES

We select events worldwide using the ANSS Comprehensive earthquake catalogue ComCat and a completeness magnitude $M_c = 4.5$ from 1980 to 2018. We only estimate the productivity for the 1490

<figure><img src="image_3.png" /><figcaption>Figure 1. Earthquake productivity in the worldwide catalogue. Dots show the distribution of the number of triggered events for $M \ge 6.5$ earthquakes using a relative magnitude threshold $\Delta M = 2$. The solid line is the exponential law with parameter $\Lambda_2$, the mean number of triggered events derived from the data. The histogram shows the Poisson distribution with parameter $\Lambda_2$. Inset shows the cumulative productivity distributions for primary and secondary triggering events.</figcaption></figure>

$M_{triggering} \ge 6.5$ earthquakes. They are associated with 6437 triggered events using a relative magnitude threshold $\Delta M = 2$. The choice of these values is not critical as long as $M_{triggering} \ge M_c + \Delta M$ and the number of triggering and triggered events is sufficient for statistical robustness (Note 2 of the Supporting Information).

Fig. 1 shows the distribution of the number of triggered events as well as the exponential and the Poisson distributions with the same rate parameter $\Lambda_2 = 6437/1490$. This comparison indicates that the productivity distribution is consistently described by an exponential law using the mean number of triggered events as a single rate parameter. Earthquakes triggered by primary events can themselves be secondary triggering events and so on until the final branches of the hierarchical clustering trees (see examples in Note 3 of the Supporting Information). The inset of Fig. 1 shows that the productivity distributions are invariant to the level of the triggering event. Thus, there is a consistent cascade of triggering and the exponential function appears to govern the productivity of all $M \ge 6.5$ earthquakes.

The productivity distribution remains exponential when the relative magnitude threshold $\Delta M$ increases from 1 to 2.6 (Fig. 2a). As expected, the mean values $\Lambda_{\Delta M}$ is decreasing according to the b-value of the earthquake-size distribution (Fig. 2b). In addition, keeping a constant relative magnitude threshold $\Delta M$, the distribution of triggered events and its mean value $\Lambda_{\Delta M}$ are almost the same regardless of the magnitude of triggering earthquakes (Figs 2c and d). Hence, both the exponential shape of the productivity distribution and the clustering factor $\Lambda_{\Delta M}$ can be considered as generic properties of earthquakes whatever their size.

The exponential shape of the productivity distribution is observed in the entire seismogenic layer and Fig. 3 shows the dependence of the clustering factor $\Lambda_2$ on depth (examples of productivity distributions using different depth ranges are shown in Note 4 of the Supporting Information). The rapid decay rate of the $\Lambda_2$-value with
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Dependence of earthquake productivity on magnitude ranges in the worldwide catalogue. (a) Distribution of the number of triggered events of $M \ge 6.5$ earthquakes using a relative magnitude threshold $\Delta M \in \{1, 1.2, \dots, 2.6\}$. (b) The mean number of events $\Lambda_{\Delta M}$ with respect to the relative magnitude threshold $\Delta M$. (c) Distribution of the number of triggered events with respect to the magnitude $M_{\text{triggering}}$ of the triggering event using a relative magnitude threshold $\Delta M = 2$. We take $M_{\text{min}} \le M_{\text{triggering}} < M_{\text{max}}$ with $M_{\text{max}} = M_{\text{min}} + 0.2$ and $M_{\text{max}} \in \{6.5, 6.7, \dots, 8.1\}$. (d) The average number $\Lambda_{\Delta M}$ of triggered events with respect to the magnitude of the triggering event for $\Delta M \in \{1, 2, 2.5\}$.</figcaption></figure>

respect to depth from 10 to 100 km indicates that the productivity is impacted by rock properties and environmental parameters governing the rupture process.

In order to validate these observations made worldwide, we analyse the productivity distributions in seven areas using regional seismic catalogues (Fig. 4). In these catalogues, lower completeness magnitudes allow us to investigate the productivity distribution in smaller magnitude ranges. The exponential behavior is ubiquitous
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Earthquake productivity with respect to depth in the worldwide catalogue. Average number of triggered events with respect to depth for $M \ge 6.5$ earthquakes using a relative magnitude threshold $\Delta M = 2$. Triggering events are chosen with respect to depth using an overlapping sliding window of 100 events with a step of 50 events. Horizontal and vertical errorbars show the depth interval and the 90 per cent credibility intervals of the likelihood function, respectively.</figcaption></figure>

in all areas but the mean number $\Lambda_2$ of triggered events varies from one region to another (Fig. 4a). It could be a consequence of the earthquake depth distribution, but also of other seismogenic properties such as creep, stress regime or heat flux (Narteau et al. 2009; Zaliapin & Ben-Zion 2016). In subduction zones along Japan, Kamchatka and New Zealand the $\Lambda_2$-value decreases with the depth of the triggering events from 10 to 100 km (Figs 4a-c). In California, where strike-slip faulting prevails, the productivity also drops with depth and the decrease in $\Lambda_2$-value is systematically observed in the brittle layer from 5 to 13 km, where most of the earthquakes occur (Fig. 4d). Finally, to determine the proportion of foreshocks, we take $\Delta M = 0$ (i.e. $M_{triggered} \ge M_{triggering}$) to compare the average productivity of $M \ge 4.5$ earthquakes in the different areas (Note 5 in the Supporting Information). The largest value is observed in Kamchatka (0.12), the smallest in New Zealand and in the worldwide catalogue (0.06).

## 4 CONCLUDING REMARKS

To explore the underlying dynamics of earthquake clusters in space and time, the basis of our approach is to attach to each event a productivity, that is a number of directly triggered events. In addition to the Omori-type decay in seismicity rate which is generally described by non-stationary Poisson process, the number $N$ of events triggered by a single earthquake follows a Poisson distribution with mean $\lambda$

$$p_p(N = k) = \frac{\lambda^k \exp(-\lambda)}{k!} \quad (2)$$

We propose here that this mean (expectation) number is not deterministic. Instead, it is a random variable that follows an exponential distribution

$$p_e(\lambda) = \frac{1}{\Lambda_{\Delta M}} \exp\left(-\frac{\lambda}{\Lambda_{\Delta M}}\right). \quad (3)$$

Its mean $\Lambda_{\Delta M}$ is the clustering factor. This $\Delta M$-productivity law is independent of the magnitude of triggering events and characterizes the earthquake clustering process across scales according to the relative magnitude threshold. Combining eqs (2) and (3), the unconditional probability density function of $N$ writes

$$\begin{aligned} p(N = k) &= \int_{0}^{\infty} \frac{\lambda^k \exp(-\lambda)}{k!} \times \frac{1}{\Lambda_{\Delta M}} \exp\left(-\frac{\lambda}{\Lambda_{\Delta M}}\right) d\lambda \\ &= \frac{1}{k! \Lambda_{\Delta M}} \times \left(1 + \frac{1}{\Lambda_{\Delta M}}\right)^{-(k+1)} \Gamma(k+1) \quad (4) \\ &= \frac{1}{1 + \Lambda_{\Delta M}} \left(\frac{\Lambda_{\Delta M}}{1 + \Lambda_{\Delta M}}\right)^k. \end{aligned}$$

where $\Gamma(x) = \int_0^\infty t^{x-1} \exp(-t) dt$ is the Gamma function. Then, the number $N$ of triggered events follows a geometric distribution with the mean $\Lambda_{\Delta M}$, which naturally explains the exponential behaviors observed worldwide and in various active tectonic settings (Figs 1, 2 and 4a).

Using hierarchical clustering trees, the persistence of the exponential behaviour of the productivity law across the different levels of hierarchy (Fig. 1) and over a wide range of earthquake size (Fig. 2) demonstrates the relevance and validity of the hypothesis underlying epidemic models of seismicity. Nevertheless, when considering
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Earthquake productivity in regional catalogues. (a) Distribution of the number of triggered events in different seismic regions. (b–d) Average number of triggered events with respect to depth in Kamchatka (b), California (c), New Zealand (d), Italy (e) using a relative magnitude threshold $\Delta M = 2$.</figcaption></figure>

the variability in the number of triggered events, an exponential distribution has never been considered before, either for theoretical predictions or in probabilistic model inversion. Our results suggest that such a law could improve the likelihood of epidemic models and allow a better understanding of the clustered nature of earthquake sequences.

In taking this step, we analyse synthetic sequences produced by these epidemic models of seismicity to check the validity of our statistical procedure. Numerical tests demonstrate that our hierarchical declustering method is able to recover both the predefined Poisson and exponential (geometric) distributions of the productivity (Note 6 in Supporting Information). Hence, we can reject the hypothesis that the exponential behaviour observed in real catalogues is an artefact of our declustering procedure.

Our investigations show similar results regardless of the nearest-neighbours selection method (Note 7 in the Supporting Information). Indeed, the exponential behavior emerges as soon as the threshold distance between statistically related events remains small enough (Note 8 in the Supporting Information). However, large

threshold value are required to distinguish between different regimes and we show here that the best compromise value can be obtained without fitting approximation by comparing the clustering property of real and random catalogues.

Our findings indicate that the perturbation induced by an earthquake produces only a small number of events. They even show that the most likely number of events triggered is zero. This counterintuitive observation must be considered with respect to the dependence of the productivity distribution on $\Delta M$ (eq. 3). As the clustering factor $\Lambda_{\Delta M}$ increases as $10^{b\Delta M}$ (Fig. 2b), the productivity distribution keeps its exponential shape but is rescaled by the same factor. The single exponential regime observed in worldwide and regional catalogues can then be considered as reflecting the overall strength of the rock material, regardless of the magnitude of the perturbation. Beyond the individual contribution of each earthquake to the triggering of future events, the branching process within cascades of seismicity is responsible for the density and spatio-temporal patterns of seismic clusters. In all cases, there is no reason to assume a Poisson distribution to take into account the variability in the number of earthquakes triggered by an event.

We show that the clustering factor $\Lambda_{\Delta M}$ decreases with depth (Figs 3 and 4b–e) and can now be used as a new statistical tool for the characterization of major fault systems, as it is currently done using earthquake-size distributions or aftershock sequences (Narteau *et al.* 2009; Shebalin & Narteau 2017). In addition, the behaviours observed in different regions of the world clearly indicate that the relative number of earthquakes with no triggered events increases as the transition from brittle to ductile behavior is approached (Figs 4b–e). The earthquake productivity and the clustering factor might then be exploited to explore faulting mechanics at depth where continuous and discontinuous deformations coexist.

## ACKNOWLEDGEMENTS

We thank Agnès Helmstetter for its help concerning the spatio-temporal ETAS model. We also thank Jiancang Zhuang and Shyam Nandan for their constructive evaluations that helped to substantially improve the manuscript. Peter Shebalin and Sergey Baranov acknowledge financial support from state contracts within research projects AAAA-A19-119011490127-6 and 075-01304-20. Clément Narteau acknowledges financial support from the UniversityEarthS LabEx program of Sorbonne Paris Cité (grants ANR-10-LABX-0023 and ANR-11-IDEX-0005-02) and the French National Research Agency (grant ANR-17-CE01-0014/SONO).

## REFERENCES

* Baiesi, M. & Paczuski, M., 2004. Scale-free networks of earthquakes and aftershocks, *Phys. Rev. E*, **69**, 066106.

* Console, R., Murru, M. & Lombardi, A.M., 2003. Refining earthquake clustering models, *J. geophys. Res.*, **108**(B10), doi:10.1029/2002JB002130.

* Davidsen, J., Gu, C. & Baiesi, M., 2015. Generalized Omori-Utsu law for aftershock sequences in southern California, *Geophys. J. Int.*, **201**(2), 965-978.

* Hainzl, S. & Marsan, D., 2008. Dependence of the Omori-Utsu law parameters on main shock magnitude: observations and modeling, *J. geophys. Res.*, **113**(B10), doi:10.1029/2007JB005492.

* Hainzl, S., Zakharova, O. & Marsan, D., 2013. Impact of aseismic transients on the estimation of aftershock productivity parameters, *Bull. seism. Soc. Am.*, **103**(3), 1723-1732.

* Helmstetter, A., 2003. Is earthquake triggering driven by small earthquakes?, *Phys. Rev. Lett.*, **91**(5), 058501.
Helmstetter, A. & Sornette, D., 2002. Subcritical and supercritical regimes in epidemic models of earthquake aftershocks, *J. geophys. Res.*, **107**(B10), 2237.

Holschneider, M., Narteau, C., Shebalin, P., Peng, Z. & Schorlemmer, D., 2012. Bayesian analysis of the modified omori law, *J. geophys. Res.*, **117**(B6), B06317.

Kagan, Y.Y. & Knopoff, L., 1981. Stochastic synthesis of earthquake catalogs, *J. geophys. Res.*, **86**(B4), 2853-2862.

Marsan, D. & Helmstetter, A., 2017. Single-link cluster analysis as a method to evaluate spatial and temporal properties of earthquake catalogues, *J. geophys. Res.*, **122**(7), 5544-5560.

Marsan, D. & Lengliné, O., 2008. Extending earthquakes' reach through cascading, *Science*, **319**(5866), 1076-1079.

Narteau, C., Byrdina, S., Shebalin, P. & Schorlemmer, D., 2009. Common dependence on stress for the two fundamental laws of statistical seismology, *Nature*, **462**, 642–645.

Ogata, Y., 1989. Statistical model for standard seismicity and detection of anomalies by residual analysis, *Tectonophysics*, **169**(1-3), 159-174.

Shcherbakov, R. & Turcotte, D.L., 2004. A damage mechanics model for aftershocks, *Pure appl. Geophys.*, **161**, 2379-2391.

Shcherbakov, R., Turcotte, D.L. & Rundle, J.B., 2004. A generalized Omori's law for earthquake aftershock decay, *Geophys. Res. Lett.*, **31**(L11613), doi:10.1029/2004GL019808.

Shebalin, P. & Narteau, C., 2017. Depth dependent stress revealed by aftershocks, *Nat. Commun.*, **8**(1), 1-8.

Utsu, T., 1969. Aftershocks and earthquake statistics (i): some parameters which characterize an aftershock sequence and their interrelations, *J. Faculty Sci., Hokkaido Univ., Ser. VII (Geophys.)*, **2**, 129-195.

Utsu, T., 1970. Aftershocks and earthquake statistics (ii): further investigation of aftershocks and other earthquake sequences based on a new classification of earthquake sequences, *J. Faculty Sci., Hokkaido Univ., Ser. VII (Geophys.)*, **3**, 197-266.

Utsu, T., Ogata, Y. & Matsu'ura, R., 1995. The centenary of the Omori formula for a decay law of aftershocks activity, *J. Phys. Earth*, **43**, 1-33.

Wang, Q., Schoenberg, F.P. & Jackson, D.D., 2010. Standard errors of parameter estimates in the etas model, *Bull. seism. Soc. Am.*, **100**(5A), 1989-2001.

Werner, M.J. & Sornette, D., 2008. Magnitude uncertainties impact seismic rate estimates, forecasts, and predictability experiments, *J. geophys. Res.*, **113**(B8), doi:10.1029/2007JB005427.

Zaliapin, I. & Ben-Zion, Y., 2013. Earthquake clusters in southern California I: identification and stability, *J. geophys. Res.*, **118**(6), 2847-2864.

Zaliapin, I. & Ben-Zion, Y., 2016. A global classification and characterization of earthquake clusters, *Geophys. J. Int.*, **207**(1), 608-634.

Zaliapin, I., Gabrielov, A., Keilis-Borok, V. & Wong, H., 2008. Clustering analysis of seismicity and aftershock identification, *Phys. Rev. Lett.*, **101**(1), 018501.

Zhuang, J., Ogata, Y. & Vere-Jones, D., 2002. Stochastic declustering of space-time earthquake occurrences, *J. Am. Stat. Ass.*, **97**(458), 369-380.

Zhuang, J., Ogata, Y. & Vere-Jones, D., 2004. Analyzing earthquake clustering features by using stochastic reconstruction, *J. geophys. Res.*, **109**(B5), doi:10.1029/2004JB003157.

Zhuang, J., Chang, C.-P., Ogata, Y. & Chen, Y.-I., 2005. A study on the background and clustering seismicity in the Taiwan region by using point process models, *J. geophys. Res.*, **110**(B5), doi:10.1029/2004JB003157.

## SUPPORTING INFORMATION

Supplementary data are available at GJI online.

### Supplementary Information.pdf

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.