%\documentclass[twoside]{article}
\documentclass[a4paper, 12pt, oneside]{article}
%\usepackage[silent,pdf,par,rus]{elxpaper}
\usepackage{color}
\usepackage[pdfpagemode=UseNone]{hyperref}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage[pdftex]{graphicx}    

%\journalid{BS}{4}{1}{16}
%\citnumber{}{1}
%\paperid{}
%\papertype{RAR}
%\cpright{}{2019}
%\lefthead{}
%\righthead{}
%\received{20}{November}{2019}
%\accepted{25}{November}{2019}
%\published{11}{December}{2019}
%\authaddr{Institute of earthquake prediction theory and mathematical geophysics of Russian academy of sciences (shebalin@gmail.com)}

\begin{document}
	\title{Epidemic type stochastic model of seismicity with exponential distribution of the earthquake productivity}
	
	\author{Baranov S.V.$^3$, Gvishiani A.D.$^{1,2}$, Narteau C.$^4$, Shebalin P.N.$^5$}
	\date{%
		$^1$Geophysical Center of Russian Academy of Sciences, Moscow, Russia.\\%
		$^2$Institute of the Physics of the Earth, Russian Academy of Sciences, Moscow, Russia.\\%
		$^3$Kola Branch, Federal Research Center Geophysical Survey, Russian Academy of Sciences, Apatity, Russia.\\%
		$^4$Institut de Physique du Globe de Paris, Universit\'e de Paris, France\\%
		$^5$Institute of Earthquake Prediction Theory and Mathematical Geophysics, Moscow, Russia.\\[2ex]%
		\today
	}
	%\author{1}{Baranov}{S.V.}{}{1}{}{}
	%\author{2}{Gvishiani A.D.}{S.V.}{}{2}{3}{}
	%\author{3}{Narteau}{C.}{}{4}{}{}
	%\author{4}{Shebalin}{P.N.}{}{5}{}{}
	\maketitle
	%\affil{1}{������������� ����� ���}{������}{}{������}
	
	%\keywd{����������� �����}
	%\keywd{����������� ����������}
	%\keywd{�������������� ����������}
	%\keywd{���������������� �����}
	%\keywd{������������� ���������}
	%\keywd{������������}
	%\keywd{CrossRef}
	%\keywd{eLibrary}
	
	\abstract{
		Studying hierarchical structure of aftershock sequences of the three largest earthquake of the last decade, we show that the earthquake productivity, the number of ''offspring'' events per ''parent'' event, counted in a fixed magnitude band relative to the parent event, unexpectedly follows exponential distribution with a mode at zero.This finding extends previous results about exponential distribution of the number of aftershocks of the largest worldwide earthquakes. The most popular stochastic model of seismicity ETAS, Epidemic Type Aftershock Sequences, predicts Poisson distribution of the earthquake productivity with a pronounced non-zero mode. We construct an alternative model, incorporating the found property of the earthquake productivity, and compare two models for the three aftershock sequences. We estimate parameters of the models using data of the first 48 hours of the sequences. Next, for each of the models, we simulate a set of synthetic earthquake catalogues and compare the average cumulative distributions of the times of events for the two models.The ETAS model in the three cases largely overestimates the number of events in the interval 48 hours to 365 days, while our alternative model gives a satisfactory cumulative time distribution for this interval. We conclude that exponential distribution of the earthquake productivity seems to be an important property of the seismic process, although we are not able yet to suggest a plausible mechanism.}
	
	\section{Introduction}
	An important feature of seismicity is the occurrence of space-time clusters
	demonstrating that earthquakes interact with each other. Focusing on the way in which a sequence of earthquakes develops over space and time, we may consider every event as the trigger for the perturbation of the state of stress in some area around the event. The zones of such a perturbation intersect in time and space, and each new event may be considered as an ''offspring'' of all preceding earthquakes. For a Poisson process such epidemic behavior may be described  as a superposition of the processes, initiated by every event. Since the advent of epidemic models of seismicity, the productivity has become
	a major issue because it determines the increase in 
	seismicity rate after each earthquake\cite{Kaga81,Ogat89,Helm02}. In 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=K10^{\alpha m}$. The value $\alpha$ was estimated in a range from 0.5 to 2 \cite{Cons03,Zhua04,Zhua05,Hain08,Wern08,Wan10,Hain13}. This value is usually close to the observed $b$-value, the slope of the earthquake-size distribution\cite{Helm03}. 
	However, these estimates remain uncertain due to the difficulty of isolating 
	the relative contributions of successive events in a sequence. 
	This task remains in early stage despite the diversity of declustering methods implemented in the past. 
	Two stochastic approaches have been suggested to find causal links within cascades 
	of triggered seismicity. 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 ETAS\cite{Zhua02}. A second approach is model independent. 
	A linear contributions of each earthquake to the overall
	seismicity rate is assumed\cite{Mars08}. Those two approaches suppose that every new earthquakes is an ''aftershock'' of all preceding events, and the goal is to estimate the impact of each preceding event on each subsequent event in terms of probability. An alternative approach \cite{zali08, zali13} goes directly to consider a tree of events in which each event may be a ''parent'' of several later events, but it can be an ''offspring'' of only one earlier event. Technically the ''parent'' is found as a ''nearest neighbor'' using proximity functions in
	time-space-magnitude domains\cite{Baie04,zali08}. Here we shall call ''parent'' events triggering events, and their ''offsprings'' triggered events.
	
	All these methods confirm the dependency of the productivity on the magnitude of the triggering event. However, it was found that within this dependency there is 
	a huge variability in the number of triggering events in the seismic catalogs\cite{Mars17}. Recently it was found that unexpectedly this variability may 
	be described by exponential distribution \cite{shb_Dokl2018} with maximum at 0. This result was obtained for the global statistics of the number of aftershocks from earthquakes of magnitude 6.5 and above. Aftershocks of magnitude above a threshold, defined as the magnitude of the main shock minus 2, were counted. 
	
	Here we study a distribution of 
	the productivity in a tree of aftershocks for three largest earthquakes of the last decade: 11 March 2011, Mw=9.1, Tohoku earthquake; 27 February 2010, Mw=8.8, Chile earthquake; 11 April 2012, Mw=8.6, Sumatra earthquake. We find that for each 
	sequence the distribution of the productivity also tends to an exponential form. Thus, exponential distribution of earthquake productivity seems to be a general property of seismicity. Exponential form of the distribution of the real number of triggered events per triggering event contradicts with an expected form of Poisson distribution. To better understand this issue we compare two models for the three aftershock sequences.
	
	\section{Earthquake productivity}
	
	We define earthquake productivity using ''delta-analysis'' \cite{zali08, zali13}, in a magnitude band of a given width $\Delta M$ relative to the magnitude of each triggering event. Thus, the productivity is a property of each earthquake. 
	
	\subsection {Tree of clustered earthquakes}
	To count the productivity values we decompose the earthquake catalogue into a hierarchical tree of pairwise links triggering-triggered events. For each pair of earthquakes $\{i,\,j\}$, we compute the proximity function\cite{Baie04},
	%
	\begin{equation}
	\eta_{ij}=
	\begin{cases}
	t_{ij}(r_{ij})^{d_{\rm f}}10^{-bm_i} & \text{{\small for} $t_{ij}>0$},\\
	+\infty & \text{{\small for} $t_{ij}\le 0$}.
	\end{cases}
	\label{eq: prox}
	\end{equation}
	%
	where $t_{ij}=t_j-t_i$ is the inter-event time,
	$r_{ij}$ the spatial distance between the epicenters,
	$m_i$ the magnitude of event $i$,
	$d_{\rm f}$ the fractal dimension of the epicenter distribution and
	$b$ the slope of the earthquake-size distribution.
	
	Using the proximity function for each event we find the preceding nearest-neighbor. In case the $\eta$ value exceeds a threshold  $\eta_0$, the event is considered as a background event because it has no triggering event. For each sequence we computed $\eta_0$ using the original technique of \cite{zali08,zali13}, approximating the distribution of the nearest-neighbor values $log(\eta)$ by a mixture model of two Gaussian distributions, one modeling independent events, the other causally-related events.  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$.
	
	Accordingly, all clusters are built from a primary triggering 
	event. There is a single path from any earthquake in a cluster to the corresponding primary event. Primary events by definition are ''background events''. Background events with the largest magnitude in the cluster are mainshocks. But a triggered event may have a larger magnitude than its triggering events. In this case the main shock is not a primary event and nor background event. Those definitions describe an aftershock sequence as a hierarchical tree of events. This is slightly different from a standard definition of the foreshocks -- main shock -- aftershock series in which the majority of the aftershocks (and also foreshocks) are linked directly to the main shock. In order to avoid the excessive influence of the definition of the proximity function on our analysis, for each of the three series we consider not a single hierarchical tree of events, including the mentioned earthquakes, but all earthquakes in a simple spatial region with a form of stadium \cite{bshFZ22017} in a time interval of 365 days after the earthquake.
	
	For each triggering event, we count the number $N$ of triggered events at the lower
	hierarchical level using a relative magnitude threshold $\Delta M$ (i.e.
	$M_{\rm triggering}-M_{\rm triggered}<\Delta M$). This number $N$ is
	defined as the productivity. 
	
	\subsection {Distribution of 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_0(\Delta M)$, we call the clustering factor. The productivity may vary from place to place, and also from sequence to sequence.  Recently it was found \cite{shb_Dokl2018} that in a global scale the productivity has a distribution of the exponential shape with maximum at 0. It was also shown that this shape does not depend on the magnitude of the triggering earthquakes nor on the width of the considered magnitude band. For the exponential distribution the clustering factor is an important parameter, as it is a single parameter of this distribution.
	
	Here we concentrate on the productivity distribution within aftershock sequences of large earthquakes. It is important to know whether the exponential form is also characteristic in much more homogeneous conditions in comparison to the global variability of the aftershock sequences considered earlier \cite{shb_Dokl2018}.
	
	The estimated completeness magnitude for the Sumatra and Chile sequences is 4.5, for Tohoku 5.0. For the analysis we use the value $\Delta M=1.5$. This choice allowed to choose the minimum magnitude of triggering earthquakes $M_tr$ 6.0 for the Sumatra and Chile sequences and 6.5 for the Tohoku sequence. Figure 1 shows the histograms of the productivity in comparison to exponential and Poisson distribution. We note that exponential distribution is the distribution of a real value, and the actual productivity is an integer. We may interpret this challenge by supposing that the productivity is an internal property of each earthquake, similar to its magnitude. Specific realizations of the productivity is an integer. It is natural to suppose that the specific values have Poisson distribution with a rate equal to the ''internal'' productivity. The existing epidemic models of seismicity imply that the internal productivity is constant. The difference between the ''internal'' and the real productivity may explain some distortion of the exponential distribution. However it is clear from the figure that for all the three sequences exponential distribution is preferable in comparison to Poisson distribution.
	
	\begin{figure}
		\noindent\includegraphics[width=\linewidth]{fig1}
		\caption{
			{\bf Earthquake productivity for three aftershock sequences}.  
			a) 27 February 2010, Mw=8.8, Chile earthquake, b) 11 March 2011, Mw=9.1, Tohoku earthquake, c) 11 April 2012, Mw=8.6, Sumatra earthquake. 
			Dots show the number of triggered events for
			$M\geq6.5$ (Tohoku) and  $M\geq6.0$ (Chile and Sumatra) triggering events using a relative magnitude threshold
			$\Delta M=1.5$. The solid line is the exponential law with parameter
			$\Lambda_0$, the clustering factor. The histogram shows the Poisson distribution with parameter
			$\Lambda_0$. 
		}
		\label{fig1}
	\end{figure}
	
	\section {Two alternative epidemic models of seimicity}
	We test here whether the found property is important for modeling the seismicity. We compare two epidemic models, the classic ETAS model and a similar model in which the constant internal productivity is replaced by a random internal productivity with exponential distribution.
	Using an interval of 2 days right after the large earthquakes 
	%%and also data prior to them, we estimate parameters of the two models
	. 
	Using those estimates we construct two versions of a synthetic catalog. Repeating simulations many times, we calculate an average cumulative number of events as a function of time in the interval 0 to 365 days.
	Finally we compare the results for the two models with the real data. 
	
	\subsection{ETAS model}
	The classic ETAS model \cite{Kaga81,Ogat89} considers seismicity as a non-stationary Poisson process with a rate described by the equation 2.
	
	\begin{equation}
	\lambda_{ETAS}(t) = \mu + K_0\sum_{t_i<t} \frac{10^{\alpha(M_i - M_0)}}{(t-t_i+c)^p},
	\label{eq: etas}
	\end{equation}
	
	\noindent where $t_i$ and $M_i$ are the time and the magnitude of the $i$-th earthquake, $\mu$, $K_0$, $\alpha$, $с$ and $p$ are parameters, $M_0$ magnitude threshold for counting events. Parameters $с$ and $p$ describe temporal decay of the triggered events according to the empirical Omori-Utsu law\cite{Utsu65,Utsu70}. Parameter $\alpha$, as discussed above, is usually close to the $b$-value of the Gutenberg-Richter relation. Parameter $\mu$ is the background seismicity rate, which is assumed constant.
	
	To estimate parameters of the model we use a standard maximum likelihood procedure \cite{Ogat89}. Earthquake catalogs are not complete right after large earthquakes. To minimize the impact of this effect we omit in the analysis the interval upto 0.05 days after the earthquake. 
	
	For the synthetic catalog simulations there are two equivalent ways \cite{Zhua04}.  The first is to  treat the ETAS model as a single point process, with the probability of an event at each
	point in time reflected in the conditional intensity which contains a component
	of background and a component of triggered seismicity contributed by all past
	events in the history of the process. The second way is to simulate the background events as a Poisson process, and
	then recursively simulate the aftershocks resulting from each of these background
	events in turn. Finally, all events are combined and put into the correct temporal
	order of occurrence.   The first way is not appropriate for our modification of the ETAS model, because the internal productivity is random. For this reason we apply here the second way. In both methods the magnitude of events is simulated independently as a random number above $M_0$ with distribution defined by Gutenberg-Richter relation.
	
	\subsection{ETAS model with random productivity with exponential distribution}
	We suggest here a simple modification of the ETAS model we shall call EP model. This model modifies the recursive definition of the ETAS model described above. The model represents seismicity as a sum of Poisson processes with the decay according to the Omori law, initiated by background events. Each new event  also initiates a similar decaying process. For the synthetic catalogue simulations for each new event we randomly generate its internal productivity $\Lambda_i$ according to the exponential distribution 
	$$p(\Lambda_i)=\frac{1}{\Lambda_0}e^{\lambda_i/\Lambda_0}.$$ 
	
	Then we simulate each branch of the non-stationary process with the rate defined by the equation \cite{Hols12}:
	
	\begin{equation}
	\lambda_{EPi}(t) = \frac{\Lambda_i}{(1+(t-t_i)/c)^p}.
	\label{eq: ep}
	\end{equation}
	
	
	In simulations we use Bayesian estimates of the parameters $c$ and $p$ with Gaussian priors \cite{bshFZ42019} in the interval (0.05,2) days after the large earthquakes. The clustering factor $\Lambda_0$ is also estimated in this interval simply as the average productivity. 
	
	This number is corrected to the interval of 365 days 
	using a multiplier 
	$$U=\frac{U(0.05,365)}{U(0.05,2)},$$
	where $U(a,b)=\int_a^b{(1+x/c)^{-p}dx}$. Like in ETAS model, magnitudes of events are independently assigned according to the Gutenberg-Richter relation. We use a Bayesian estimate of the  $b$-value with Gaussian priors \cite{bshFZ42019}.
	
	\section{Data}
	
	We used data from ANSS ComCat earthquake catalog provided by USGS. Aftershocks of M8.8 Chile earthquake of 2010 were taken for a year after the mainshock from the circle of radius 900 km surrounding its epicenter. Aftershocks of M9.1 Tohoku earthquake of 2011 were taken for a year after the mainshock from the circle of radius 1000 km surrounding its epicenter. Aftershocks of M8.6 Sumatra were taken for a year after the mainshock from the circle of radius 700 km surrounding its epicenter. We did not apply any other special aftershock selection procedure. 
	
	\section{Results and discussion}
	
	For each of the three sequences we have performed 2500 simulations of the synthetic catalogs using the two models. We estimated  ETAS and EP models parameters using data for 2 days after the mainshocks (Table \ref{tab1}).  The average cumulative number of event has been calculated as a function of time. Results are shown in Figure \ref{fig2}. 
	
	\begin{table}[htp]
		\caption{ETAS and EP models parameters estimated using data for 2 days after mainshock.}
		\label{tab1}
		\begin{tabular}{lll}
			\textbf{Mainshock}  & \textbf{$b$-value} & \textbf{Model parameters} \\
			\hline
			Chili, 2010, M8.8   &$b=1.36$  & ETAS: $\mu=0.65$, $c=0.08$, $p=1.26$, $K_0=0.028$, $\alpha=1.82$    \\
			&          & EP: $\Delta M=1.5$, $\mu_{EP}=0.1$, $\Lambda=0.2$, $c=0.034$, $p=1.01$ \\
			Tohoku, 2011, M9.1  & $b=1.06$ & ETAS:  $\mu=0.81$, $c=0.17$, $p=1.56$, $K_0=0.04$, $\alpha=1.74$  \\                                       &          & EP:  $\Delta M=1.5$, $\mu_{EP}=0.2$, $\Lambda=1$, $c=0.04$, $p=1.02$  \\
			Sumatra, 2012, M8.6 &$b=1.26$  & ETAS:  $\mu=0.43$, $c=0.1$, $p=1.52$, $K_0=0.01$, $\alpha=1.846$  \\
			&          & EP:  $\Delta M=1.5$, $\mu_{EP}=0.1$, $\Lambda=0.7$, $c=0.05$, $p=1.36$ \\
			\hline
		\end{tabular}
	\end{table}
	
	\begin{figure}
		\noindent\includegraphics[width=0.8\linewidth]{fig2}
		\caption{
			{\bf Cumulative number of seismic events  for three aftershock sequence, real and forecasted by models ETAS and EP}.  
			a) 27 February 2010, Mw=8.8, Chile earthquake, b) 11 March 2011, Mw=9.1, Tohoku earthquake, c) 11 April 2012, Mw=8.6, Sumatra earthquake. 
			Circles show the real cumulative number of event of magnitude 
			$M\geq6.5$ (Tohoku) and  $M\geq6.0$ (Chile and Sumatra), dashed line the average for 2500 catalogue simulation of the ETAS model, solid line 
			average for 2500 catalogue simulation of the EP model. Parameters for both models have been estimated in the interval (0.05,2) days afters the corresponding large eartquakes, their value are given in Table.
		}
		\label{fig2}
	\end{figure}
	We see that the ETAS model gives a drastic overestimation of the earthquake rates at later times, while the EP model is quite acceptable in the whole forecasting period from 2 to 365 days. Two major reasons explain this as we can suppose. First, the mean value of the exponential distribution is smaller than the mean of the Poisson distribution with an equal parameter. Thus, the EP model should predict smaller rates in comparison to the ETAS model. Second, in all three cases the background rate of the ETAS model estimated at the beginning of the sequence is much higher than estimated in the interval of 365 days. In the EP model the calculations of the background rate give similar values for the beginning and for the whole sequence. 
	
	One of interesting results of this paper is exponential form of the productivity 
	distribution within single aftershock sequences from large earthquakes. Earlier this property was established in global and regional scales \cite{shb_Dokl2018,shbFZ52019}. 
	
	The obtained results demonstrate advantages of the suggested modification of the ETAS model. It is too early to come to final conclusions about those two model, however we obviously may conclude that exponential shape of the productivity distribution is an important property of the seismic process, although we are not able yet to suggest a plausible mechanism.
	
	
	\section*{Acknowledgments}
AG acknowledges financial support from ... SB developed software for the data analysis with support of Russian Foundation of Basic Researches, project 19-05-00812.

	
%	\bibliographystyle{naturemag}
%	\bibliography{bib_prod.bib}
	\input{shebalin_RJES2019.bbl}
\end{document}







