Topological properties of epidemic aftershock processes

TL;DR

This paper explores the topological properties of epidemic aftershock processes modeled by the ETAS framework, linking branching process theory with seismic event clustering to better understand earthquake triggering mechanisms.

Contribution

It extends the ETAS model from classical branching processes, analyzing how its topological properties depend on key parameters and providing a null-model for empirical seismic clustering analysis.

Findings

Topological properties depend on average branching ratio and exponent ratio.

ETAS model bridges Poisson and scale-free tree behaviors.

Provides a null-model for empirical aftershock sequence analysis.

Abstract

Earthquakes in seismological catalogs and acoustic emission events in lab experiments can be statistically described as a linear Hawkes point process, where the spatio-temporal rate of events is a linear superposition of background intensity and the aftershock clusters triggered by preceding activity. Traditionally, statistical seismology has interpreted this model as the outcome of an epidemic branching process, where one-to-one causal links can be established between mainshocks and aftershocks. Declustering techniques have been used to infer the underlying triggering trees and relate their topological properties with epidemic branching models. Here, we review how the standard Epidemic Type Aftershock Sequence (ETAS) model extends from the Galton-Watson (GW) branching processes and bridges two extreme cases: Poisson sampling and scale-free power-law trees. We report the most essential topological properties expected in GW epidemic trees: the branching probability, the distribution of tree size, the expected family size, and the relation between average leaf-depth and tree size. We find that such topological properties depend exclusively on two sampling parameters of the standard ETAS model: the average branching ratio $N_b$ and the exponent ratio $\alpha/b$ determining the branching probability distribution. From these results, one can use the memory-less GW as a null-model for empirical triggering processes and assess the validity of the ETAS model to reproduce the statistics of natural and artificial catalogs.

Equations

$$N_{AS}(\tau |m_{\mathrm{MS}})\sim \frac{k10^{m_{\mathrm{MS}}}}{(C+\tau {)}^p},$$

$$\mu (t,\mathbf{r},m)={\mu }_0(t,\mathbf{r},m)+\sum _{i|t_i<t}{\Psi }_i(m,t,\mathbf{r}|m_i,t_i,{\mathbf{r}}_{\mathbf{i}}).$$

$${\Psi }_i(m,t,\mathbf{r}|m_i,t_i,{\mathbf{r}}_{\mathbf{i}})={\rho }_m(m)n_b(m_i){\Psi }_t(t-t_i){\Psi }_r(||\mathbf{r}-{\mathbf{r}}_{\mathbf{i}}||).$$

$${\rho }_m(m):=10^{-b(m-m_c)}/(b\log (10)),$$

$$\mathbb{P}(N_i=k|n_b)=\frac{n_b^k{e}^{-n_b}}{k!}.$$

$$n_b(m_i):=k_{c}10^{\alpha (m_i-m_c)},$$

$${\rho }_{n_b}(n_b)=\frac{b}{\alpha }\frac{1}{k_c}{\left(\frac{n_b}{k_c}\right)}^{-\left(\frac{b}{\alpha }+1\right)}$$

$$\mathbb{P}(N_i=k)=\int {\rho }_{n_b}(n)\mathbb{P}(N_i=k|n)dn=\frac{b}{\alpha }\frac{{\left(N_b\left(1-\frac{\alpha }{b}\right)\right)}^{\frac{b}{\alpha }}}{k!}\Gamma \left(k-\frac{b}{\alpha },N_b\left(1-\frac{\alpha }{b}\right)\right).$$

$$p_T(K|n_b)=\frac{(n_{b}K{)}^{K-1}{e}^{-n_{b}k}}{K!}.$$

$$⟨⟨d_l⟩|N_T⟩\propto N_T^{{\gamma }_d}.$$

Full text

\journalname

JGR: Solid Earth

Centre for Mathematical Research, Campus de Bellaterra, Edifici C 08193 Bellaterra, Barcelona, Spain

\correspondingauthor

Jordi Barójbaro@crm.cat

{keypoints}

Epidemic aftershock models are interpreted as branching processes with causal-links.

All topological properties of a simple epidemic model are determined by only two parameters.

If valid, model predictions can be used to classify natural and artificial catalogs.

1 Introduction

The concept of aftershocks is traditionally associated with seismology [Utsu \BOthers. (\APACyear1995)], but similar phenomena have been observed in other natural systems and are common in many mechanical processes in rocks, composites and porous materials [Benioff (\APACyear1951), Hirata (\APACyear1987), Baró \BOthers. (\APACyear2013), Ribeiro \BOthers. (\APACyear2015), Davidsen \BOthers. (\APACyear2017)]. Aftershocks are identified in sequences of point events as a sudden increase of the activity causally linked to a previous event —usually stronger— called a mainshock. The empirical Omori-Utsu law [Utsu \BOthers. (\APACyear1995)] describes the common temporal evolution of the number of aftershocks after time $\tau$ since a mainshock of magnitude $m_{\mathrm{MS}}$ as:

[TABLE]

where $p$ is a power-law exponent usually close to 1 and $k$ a productivity factor. Additionally, aftershocks are spatially clustered, usually according to power-law decay with distance to the mainshock $r:=||\mathbf{r}-\mathbf{r_{\mathrm{MS}}}||$, being $\mathbf{r}$ and $\mathbf{r_{\mathrm{MS}}}$ the locations of the aftershock and the mainshock respectively [Guo \BBA Ogata (\APACyear1995)]. Different sophistications have been proposed accounting for more precise observations such as anisotropic spatial drift of the aftershock production and non-factorizable magnitude dependencies [Ogata \BBA Zhuang (\APACyear2006)], generalized scaling forms [Vere-Jones (\APACyear2005), Saichev \BOthers. (\APACyear2005), Davidsen \BBA Baiesi (\APACyear2016)] or more complex temporal decay forms [Davidsen \BOthers. (\APACyear2017), Baró \BBA Davidsen (\APACyear2017)]. Such details are excluded from the following mathematical and numerical developments but will be recovered in the discussion of the results.

Due to the complexity of the seismogenic mechanisms, statistical seismology considers all earthquakes —mainshocks and aftershocks— as non-isolated stochastic events in space and time, rather than the outcome of deterministic mechanical processes. Mainshocks and aftershock sequences from different mainshocks coexist in the same regions and temporal windows. As consequence, all events are customarily interpreted as point events in a stochastic point process, determined by an intensity $\mu(t,\mathbf{r},m)$ accounting for the instantaneous probability of finding an event defined by a mark —in this case, the earthquake magnitude $m$— at a spatiotemporal ($t,\mathbf{r}$) location. Simple proportional hazard models consider events to be independent, with a space-time dependent $\mu$ [Varotsos \BOthers. (\APACyear1996)]. More sophisticated hazard models take into account the correlations between events caused by aftershock production [Vere-Jones \BBA Ozaki (\APACyear1982), Ogata (\APACyear1988), Michael (\APACyear1997), Zhuang \BOthers. (\APACyear2002), Turcotte \BOthers. (\APACyear2007)]. The most simple aftershock model is to consider a linear Hawkes self-exciting process [Hawkes \BBA Oakes (\APACyear1974)] incorporating the observations of statistical seismology. In a linear Hawkes process, the intensity can be expressed as the linear superposition of a background rate $\mu_{0}$ and the individual contribution of all previous events $\{i\}$ with a triggering term $\Psi_{i}$:

[TABLE]

Notice that the linear Hawkes model assumes an additive contribution for each event in the intensity, meaning that all events can generate aftershocks. This effective stochastic process is often interpreted as the result of an epidemic or branching process [Vere-Jones \BBA Davies (\APACyear1966), Ogata (\APACyear1988), Saichev \BBA Sornette (\APACyear2004), Turcotte \BOthers. (\APACyear2007)]. Fig. 1 shows a schematic representation of a sequence of point events interpreted as the outcome of a history-dependent intensity, or Hawkes process, and its representation as a branching process. The branching process is strictly constituted by two distinct categories of events: background independent events generated entirely by the background rate ($\mu_{0}$ in Eq. (2)), and triggered events, caused by a unique preceding parent event $i$, through its individual and independent contribution to the intensity ($\Psi_{i}$ in Eq. (2)), represented as arrows in Fig. 1. Under the branching process assumption, earthquake catalogs are built as branching forests, linear superpositions of independent topological objects that here we call triggering trees. Each triggering tree is a structure of causally connected events initiated by a background event, the root of the tree (at Depth = 0 in Fig. 1), that can trigger a number $N_{1}$ of events in a first generation of aftershocks (Depth = 1 in Fig. 1). In its turn, each of the $N_{d}$ events in the $d$-th generation of aftershocks can trigger events in the $(d+1)$-th generation of aftershocks, and so on. We call leaves those events with no offspring, extinguishing a branch (dark-green circles in Fig. 1). Singlets are background events which are also leaves, i.e. have no offspring. The triggering tree is extinguished when all events in a maximum generational depth $D$ are leaves.

Assuming the validity of the branching process approach, the identification of triggering trees would provide valuable topological information of the branching and enable the direct measurement of the triggering kernel $\Psi_{i}$ [Zhuang \BOthers. (\APACyear2004), Davidsen \BOthers. (\APACyear2017)]. Hence, this representation of aftershock sequences is a useful approach, whether we can argue if strictly valid. Whilst the stochastic point process resulting from the branching process can be represented as a linear Hawkes process, the linear Hawkes processes do not require the existence of explicit one-to-one causal links between events. In theory, univocal causal links in mechanical interactive systems can be defined from an energy stability point of view. This can be shown in micromechanical models such as the viscoelastic democratic fiber bundle model [Baró \BBA Davidsen (\APACyear2018)]. In field studies, however, we have a limited capability to identify such a deterministic process. The bare statistical analysis of earthquake catalogs do not determine the explicit causal links but provide, instead, an assessment of its point process representation (2) where all terms contribute to $\mu$ with their specific weight. Advanced declustering techniques, either stochastic [Zhuang \BOthers. (\APACyear2004)] or based on the nearest neighbor distance [Baiesi \BBA Paczuski (\APACyear2004), Zaliapin \BBA Ben-Zion (\APACyear2013\APACexlab\BCnt2), Zaliapin \BBA Ben-Zion (\APACyear2013\APACexlab\BCnt1)] can be used to infer the most plausible one-to-one causal structure. This classification is never entirely free of uncertainty [Zhuang \BOthers. (\APACyear2002)] but appears to be reliable when tested against synthetic catalogs [Zaliapin \BBA Ben-Zion (\APACyear2013\APACexlab\BCnt2)]. Notice that the topological concepts discussed in the following sections are only valid under the branching process assumption and have no correspondence to the more general point process description.

Recently, the analysis of seismological catalogs as branching processes revealed significant deviations between the reconstructed clusters or triggering trees, and the branching model expectations in southern California [Zaliapin \BBA Ben-Zion (\APACyear2013\APACexlab\BCnt1)]. The same authors suggested a regional classification based upon such inconsistencies in the topological properties. Thereupon, aftershock sequences or clusters were classified into two distinct categories: bursts and swarms. Burst-like clusters were defined as clusters of events with shallow generational depth ($d$), where most activity is accumulated in the first generation of aftershocks, mostly a consequence of significantly strong mainshocks. The burst-like activity was found to be compatible with the numerical results of the Epidemic Aftershock Sequence (ETAS) model, defined in the next section, by imposing the parameters fitted from field catalogs. Burst-like sequences were linked to tectonic settings with low heat flow. On the contrary, swarm-like clusters designate aftershock sequences with deeper generational depths, usually growing with the size of the swarm. The swarm-like activity was not predicted by the ETAS model with the fitted parameters and was found to match those regions with high heat flow. The results by \citeAZaliapin2013b paved the road for a new analytical methodology based upon the topological statistics of triggering trees. This methodology will potentially lead to a better understanding of the seismogenic mechanisms behind triggering processes and improve the accuracy of stochastic point-process models with potential applicability to hazard assessment [Field \BOthers. (\APACyear2014)]. Overall, the study of topological properties of triggering trees will foreseeably be gaining more popularity in the following years, thanks to the improved refinement of event detection techniques [Shelly \BOthers. (\APACyear2016), Ross \BOthers. (\APACyear2017)]. However, few works [Saichev \BBA Sornette (\APACyear2004), Saichev \BOthers. (\APACyear2005)] have addressed the topological properties expected from these hypothetical branching models, and a global picture of the model predictions remains incomplete.

Here, we will introduce some of the most common topological features used to characterize natural triggering for the standard ETAS model, defined as a particular Galton-Watson branching model. We will revisit the results by \citeASaichev2005 and introduce general predictions regarding the concepts of family size and average leaf-depths. These results can be useful as a benchmark to validate ETAS as a null-hypothesis to natural and synthetic catalogs. As a main result, we will prove that the relation between average leaf-depths and cluster sizes depends exclusively on the probabilistic nature of the individual branching ratios, determined in the ETAS model by the ratio of parameters $b/\alpha$ and the average branching rate. This has a direct implication on the separation between swarm-like and burst-like clusters.

2 Epidemic aftershocks as Galton-Watson branching models

Epidemic Type Aftershock Sequence (ETAS) models are based on field and lab observations such as the empirical Omori-Utsu law (1) and the spatial distribution of aftershocks. Here we discuss the standard ETAS model  [Ogata (\APACyear1988)], simulating a Hawkes process where the triggering kernel is factorized in its dependencies as:

[TABLE]

Both the temporal and spatial kernel are normalized to $\int_{0}^{\infty}\Psi_{t}(t-t_{i})dt=1$ and $\int_{\mathbb{R}}\Psi_{r}(||\mathbf{r}-\mathbf{r_{i}}||)d\mathbf{r}=1$ and, therefore, have no effect on the topology of the triggering trees. The magnitudes ($m$) of the events are independent and identically distributed (i.i.d.) following the Gutenberg-Richter law [Gutenberg \BBA Richter (\APACyear1944)]:

[TABLE]

being the magnitude of completeness $m_{c}$ an effective lower-bound to the distribution. For the sake of simplicity, we consider $m_{c}$ to be also the minimum magnitude able to generate aftershocks. We also consider that the number $N_{i}$ of aftershocks generated by event $i$ is a Poisson number with a characteristic branching ratio $n_{b}:=\mathbb{E}(N_{i})$:

[TABLE]

The branching ratio is given by the aftershock production in Eq. (3):

[TABLE]

reproducing the mainshock-magnitude ($m_{i}$) dependence in the Omori-Utsu law (1) with a productivity exponent $\alpha$ usually found between 0.5 and 1 in field [Utsu \BOthers. (\APACyear1995)] and experiments [Baró \BOthers. (\APACyear2013), Davidsen \BOthers. (\APACyear2017)]. The term $k_{c}$ normalizes the aftershock production for $m_{i}=m_{c}$. Since magnitudes are independent, the number of aftershocks is equivalent to a random sampling of i.i.d. $n_{b}$ for all $m_{i}$ values. Given Eqs. (4,6):

[TABLE]

where the ratio $\alpha/b$ is bound inside the range $(0,1]$. From now on we change the notation from $k_{c}$ to the more convenient average branching ratio $N_{b}:=\langle n_{b}\rangle=k_{c}\frac{b}{b-\alpha}$. The number of first generation aftershocks for all events and all $n_{b}$ is i.i.d. as:

[TABLE]

Therefore, the standard ETAS model is a particular case of a Galton-Watson (GW) process [Pitman (\APACyear2006)], where all individual events $i$, background and triggered, have the same probability to trigger a number $N_{i}$ of events defined by an offspring distribution $\mathbb{P}(N_{i}=k):=p_{1}(k)$. In this case, for large $k$, this distribution can be approximated to a power-law $p_{1}(k)\sim k^{-\gamma_{1}}$ with the exponent value $\gamma_{1}={b/\alpha+1}$ inherited from the $n_{b}$ distribution (7). Considering $0<\alpha\leq b$, this exponent values are constraint to $2\leq\gamma_{1}$. This result agrees with the reconstructed trees in seismology [Baiesi \BBA Paczuski (\APACyear2004), Zaliapin \BBA Ben-Zion (\APACyear2013\APACexlab\BCnt2)], where $b\approx\alpha$ and the distribution of the number of first generation aftershocks —called degree distribution by \citeABaiesi2004— is a power-law with an exponent $\gamma_{1}\approx 2$.

A singular case of the ETAS model is found for $\alpha=0$. In that case, the distribution (8) becomes a Dirac delta around $k_{c}$: $\lim_{r\to 0}\rho(n_{b})=\delta(n_{b}-k_{c})$, i.e. $n_{b}$($=k_{c}=N_{b}$) is unique and all events have the same probability of generating aftershocks given by Eq. (5). The ETAS model for $\alpha=0$ is, thus, equivalent to a Poisson Galton-Watson (P-GW) branching process [Pitman (\APACyear2006)].

The numerical results in the following sections are obtained through Monte-Carlo generation of a number of aftershock sequences up to $N=10^{7}$. We use an arbitrary threshold $m_{c}=1$. The background rate and the spatio-temporal kernel parameters are ignored since they play no role in this study.

3 Tree and family sizes

A fundamental concept for the topological characterization of the branching process is the tree-size ($N_{T}$) defined here as the total number of members in an extinguished tree for $t\to\infty$: $N_{T}:=\sum_{d=0}^{D}N_{d}$. Notice that, for $N_{b}>1$, there is a non-zero probability of finding eternal trees ($D\to\infty$) with a non-defined $N_{T}$. We only consider values $N_{b}<1$ through all the discussion, imposing that all trees are extinguished at a finite time. Again, since all events $i$ are i.i.d., $N_{T}$ is i.i.d. as well, with a given $\mathbb{P}(N_{T}=K):=p_{T}(K)$. In the case of the P-GW, obtained by imposing $\alpha=0$, this distribution is known to be a Borel distribution [Pitman (\APACyear2006)]:

[TABLE]

As a matter of fact, notice that, for large $n_{b}$ and $K$, this distribution can be approximated by the exponentially tapered power-law: $p_{\>T,n_{b}}(K)\sim K^{-3/2}e^{-K(1-n_{b})^{2}/2}$. This exponent $\gamma_{N}=3/2$ is common in the avalanche size distribution of mean-field models with avalanche dynamics [Sethna \BOthers. (\APACyear1993), Zapperi \BOthers. (\APACyear1995), Fisher \BOthers. (\APACyear1997), Vespignani \BBA Zapperi (\APACyear1998), Baró \BBA Davidsen (\APACyear2018)], sometimes regarded as loop-less triggering processes, i.e. P-GW processes.

A decade ago, \citeASaichev2005 found the mathematical expression for the distribution of tree sizes ($N_{T}$). Fig. 2 shows a numerical verification of the results. The insets (Figs. 2.b,c) show in log-log scale the distributions of $N_{T}$ given $N_{b}=0.995$ and $N_{b}=0.30$ respectively for a selected collection of $\alpha/b$ values, compared with the theoretical result for the P-GW. Fig. 2.a shows the power-law exponent ($\hat{\gamma}_{N}$) estimated by maximum likelihood [Baró \BBA Vives (\APACyear2012)] for $N_{b}=0.99$ and $N_{b}=0.30$ as function of $\alpha/b$. As expected, the Borel distribution (9) is recovered for $\alpha=0$. The distribution for $0<\alpha\leq b$ is different and can be summarized in the following observations: (i) For all values $0\leq\alpha/b<1$ the original power-law behavior in (9) is preserved below a chacateristic $N_{T}$ value $N_{c}\approx(1-N_{b})^{1/(1-\alpha/b)}$ [Saichev \BOthers. (\APACyear2005)] but with a general exponent value $\gamma_{N}^{\mathrm{low}}$. (ii) This effective power-law exponent $\gamma_{N}^{\mathrm{low}}$ observed below $N_{c}$ increases from the expected $1.5$ for $\alpha/b<0.5$ towards higher values $\gamma_{N}^{\mathrm{low}}\approx 2.0$ for $\alpha/b\approx 1$ (see Fig. 2.a). The thick gray line in Fig. 2.b represents the Borel distribution for $N_{b}=0.99$, indistinguishable to a power-law with $\gamma_{N}^{\mathrm{low}}$ within the range of observation. Since the characteristic scale is high ($N_{c}>10^{4}$ for all $\alpha/b$), the exponents are used as a proxy for $\gamma_{N}^{\mathrm{low}}$. The present results for $\gamma_{N}^{\mathrm{low}}$ agree with \citeASaichev2005, which predicted a transition from $\gamma_{N}^{\mathrm{low}}=1.5$ for $b/\alpha<0.5$ to $\gamma_{N}^{\mathrm{low}}=1+b/\alpha$ for $0.5\leq\alpha/b\leq 1$ (dashed green line in Fig. 2.a). (iii) The exponential regime for large trees found in the Borel distribution becomes fat-tailed with a power-law exponent $\gamma_{N}^{\mathrm{high}}=1+b/\alpha$ for all the regime $0<\alpha<b$, reminiscent of the large scale regime in Eqs. (7,8). This power-law regime can be observed in the distributions for $N_{b}=0.30$ represented in Fig. 2.c. The exponents for $N_{b}=0.30$ are used as a proxy for $\gamma_{N}^{\mathrm{high}}$ since the selected estimation interval is considerably above $N_{c}$ except for $\alpha=b$. (iv) For $\alpha=b$, both exponents coincide and $N_{c}$ diverges. The power-law behavior has a single exponent ($\gamma_{N}\approx 2.0$) with an infinite domain. This singular scale-free sampling is usually observed in field catalogs and incited the development of the more restrictive Branching Aftershock Sequence (BASS) model introduced by \citeATurcotte2007. In this singular solution, the value of $N_{b}$ only affects the proportion of singlets $\mathbb{P}(N_{T}=0)$ (8), offsetting the distribution for all $N_{T}>0$. This can be shown by comparing the faintest lines in Fig. 2.b and Fig. 2.c.

Other common topological observables can be inferred from the sampling distribution. Eq. (8) establishes a relation between the ETAS parameters and the number of singlets, by imposing $k=0$. Given that the branching process is GW, the same value provides the average fraction of leaves in trees and a good approximation to the average family size $\langle B\rangle$. The family size ($B$) is defined by \citeAZaliapin2013b as ‘the average number of offspring over all earthquakes in the family that have at least one offspring’. This definition is mathematically identical to $B:=\frac{N_{T}-1}{N_{T}-n_{l}}$, being $N_{T}$ the cluster size and $n_{l}$ the number of leaves. Notice that, for large clusters, $B\approx\left({1-\frac{n_{l}}{N_{T}}}\right)^{-1}$. Thus, the expected value of $B$ is directly related to the probability that an event is a leaf ($n_{l}/N_{T}$), which is equivalent to $P(N_{i}=0)$ or $g(0)$ in \citeAZaliapin2013b. Notice that $P(N_{i}=0)$ and $B$ depend on $N_{b}$ and the ratio $\alpha/b$ but are essentially independent of the number of events ($N_{T}$) in the standard ETAS model when $N_{T}$ is large enough.

4 Generational depth and sizes

Some of the most commonly used topological properties of a GW process related to the generational depth of extinguished trees can be derived from the dualities between branching processes and random walks [Bennies \BBA Kersting (\APACyear2000), Pitman (\APACyear2006)]. The concept of Harris path [Harris (\APACyear1951)] is particularly useful for the topological analysis of aftershock sequences since it establishes a link between topological concepts such as size and depth with properties of a one-dimensional stochastic process. The Harris path can be interpreted as a sorted exploration of the branching tree in the following way: Consider an event $i$ in the tree (parameterized as $t-1$ in the Harris path) of depth $d(t-1):=d_{i}$ that can either be a root ($d=0$) or a triggered event from a parent with index $p$. We ask if this event $i$ has a new child $c$. If it does, we repeat the process with that first child event, now with $d(t):=d_{c}=d(t-1)+1$; if it does not, we decrease the generation by one ($d(t):=d_{p}=d(t-1)-1$), move back to the parent and ask if $p$ has another child unexplored by the Harris path. Starting from the mainshock or root of the tree, $d(t=1)=0$, all events $\{i\}$ are explored $\{N_{i}+1\}$ times before reaching $d(T)=-1$ at exactly $t=T\equiv 2N_{T}$. Leaves are identified as local maximums in the profile. We find a leaf $l$ at $t$ if $d({t})=d({t-1})+1$ and $d(t+1)=d({t})-1$, and $l$ is explored exactly once by the Harris path.

Let’s consider, for a moment, that the tree is generated through a GW process with a geometric offspring distribution. This is $p_{1,N_{b}}(k)={N_{b}}^{k}(N_{b}+1)^{-k-1}$ in terms of $N_{b}$. The probability of finding a new child, adding a step $d\to d+1$ in the Harris path, is independent of the number of previous children of the same parent and, therefore, the Harris path is equivalent to a random walk [Harris (\APACyear1951)]. Hence, the depth profile is a diffusion process that scales with time as $t^{1/2}$. All magnitudes associated with charactersitic depths such as the maximum depth $D$ and the average leaf-depth ($\langle d_{l}\rangle$) are expected to scale with the size of the tree as $N_{T}^{1/2}$. Whilst the random walk analogy is only strictly valid in the geometric case, the Harris path of a GW with another offspring distributions is asymptotically equivalent to a random walk for large tree sizes as long as $N_{b}\lesssim 1$ and the variance is well defined ($0<\sigma^{2}<\infty$) [Aldous (\APACyear1991), Pitman (\APACyear2006)]. This approximation is valid for the P-GW process, but not for the ETAS model when $\alpha/b>0.5$, rendering exponents values lower than 3 in Eq. 8 and, hence, infinite variance.

Fig. 3.a shows the numerical validation of the diffusive assumption in the case of the P-GW. On average, the dependence between $d_{l}$ and $N_{T}$ follows the relation $\langle\langle d_{l}\rangle|N_{T}\rangle\sim N_{T}^{0.5}$, which extends to significantly low values of $N_{b}$. Although not shown here, a proportional relationship is also found in other characteristic depths such as the average depth of all events $\langle\langle d\rangle|N_{T}\rangle\sim N_{T}^{0.5}$ and average maximum depth $\langle\langle D\rangle|N_{T}\rangle\sim N_{T}^{0.5}$. Notice that, similarly as how a stochastic process cannot diffuse faster than a ballistic trajectory, the maximum depth cannot be larger than the tree size, forcing the limit $\langle d_{l}\rangle<N_{T}$ which biases the average values for small trees (usually $N_{T}<20$ as seen in Fig. 3).

The topological relations between depths and sizes are more complex in the ETAS model due to the power-law sampling (8). Here we only introduce the numerical results and leave the mathematical derivation, if possible, as an open question for future works. Fig. 4.a shows the conditional average $\langle\langle d_{l}\rangle|N_{T}\rangle$ for different parameters $\alpha/b$. An scatter plot of the data pairs ($\langle d_{l}\rangle$, $N_{T}$) for $\alpha/b=0.6$ and $N_{b}=0.995$ is also shown as example. The distribution of small trees exhibits the aforementioned limit $\langle d_{l}\rangle<N_{T}$, biasing the bivariate distribution for $N_{T}<20$. For larger trees, an average power-law relation is observed beyond statistical fluctuations:

[TABLE]

This power-law relation gets distorted for low values of $N_{b}$, although one is still able to estimate a reliable power-law relation for a shorter range. The inset in Fig. 4.b shows the exponent values $\gamma_{d}$ fitted for trees with $30<N_{T}<1000$. The exponent appears to be consistent with the typical GW process ($\gamma_{d}=0.5$) for the relatively broad range $0<\alpha/b<0.3$ but drops towards lower exponent values for larger ratios. This first regime for low $\alpha/b$ values is consistent with the asymptotic limit since the effective power-law exponent $\gamma_{1}$ in Eq. 8 is higher than 3. For $0.5<\alpha/b<1$, the drift from the diffusive relationship is apparent for all values of $N_{b}$. Close to $\alpha/b=1$, the depth of triggering trees appear to be independent of $N_{T}$. The branching ratio $N_{b}$ changes the range of the distribution in $N_{T}$ as well as the dependence on the bivariate distribution. For low $N_{b}$ values, large events tend to have shallower trees than predicted by the diffusive relationship. As a consequence, the effective exponent $\gamma_{d}$ depends on $N_{b}$ for intermediate values of $\alpha/b$ but almost coincide for the extreme cases $\alpha<0.3b$ and $\alpha\approx 1$ . The exponent transition is steeper for low $N_{b}$ values.

5 Discussion

The analytical and numerical results obtained here from the branching process models highlight three direct measurements that stand out to characterize the topology of epidemic aftershock processes: the direct triggering distribution $p_{1}(k)$, the relation $\langle d_{l}\rangle$, $N_{T}$ and the size distribution $p_{N}(K)$. Their joint analysis can verify whether the relation between $\langle d_{l}\rangle$ and $N_{T}$ corresponds to the expected behavior of the specific branching model hypothesis and an ETAS model with independently fitted parameters $k_{0}$, $b$ and $\alpha$. Scalar measurements such as the fraction of singlets or leaves, the family size $B$, and the bulk average leaf-depth $\langle\langle d_{l}\rangle\rangle$, which entangle both dependencies in $N_{b}$ and $\alpha/b$, can mislead the validation or fitting of a branching model. However, the expected relationships between such scalar measurements, such as the fraction of leaves and singlets, can serve for preliminary rejection tests.

The ETAS model can potentially explain the topological variability of triggering trees reconstructed from seismological catalogs [Zaliapin \BBA Ben-Zion (\APACyear2013\APACexlab\BCnt1)] and acoustic emission experiments [Davidsen \BOthers. (\APACyear2017)]. Provided the validity of the branching model approach, unequivocal relations exist between the topological properties and the parameters of aftershock production, determined by $N_{b}$ and $\alpha/b$ in the ETAS model (8). Assuming independence between the values of $n_{b}$ —guaranteed in the present study because of the independence of $m$— all topological information is contained in the distribution of $n_{b}$ (as Eq. 7 in the ETAS). In general, the same relations can be extrapolated to other GW models with well defined $p_{1}(k)$. Beyond strict GW processes, the same results would be locally valid in models implementing spatio-temporal variations of the ETAS parameters such as the spatially-variating ETAS model (SVETAS) introduced by \citeANadan2017. In that case, the topological properties would depend on the local distribution of $N_{b}$ and $\alpha/b$. A global evaluation would report a spread in the ($\langle d_{l}\rangle$,$N_{T}$) space with a blurred power-law relationship, as reported by \citeAZaliapin2013b in southern California.

The ETAS model is a peculiar case of GW processes where the aftershocks sampling is power-law [Saichev \BOthers. (\APACyear2005)], leading to a natural cluster classification in the terms presented by \citeAZaliapin2013b. The results shown in Fig. 4 validate the hypothesis that swarms, like bursts, can appear as a consequence of event-event triggering processes, i.e. aftershocks, represented as one-to-one causal links in branching processes [Zaliapin \BBA Ben-Zion (\APACyear2013\APACexlab\BCnt1)]. The topological properties of the trees used for the classification of swarms and bursts —and, in particular, the exponent $\gamma_{d}$— differ depending on the parameter ratio $\alpha/b$ and the branching ratio $N_{b}$. Such classification is noticeably sharp in the parameter space for low $N_{b}$ ( Fig. 4.b) but smooth for high branching ratios ($N_{b}\lesssim 1$). The two classes are found in the extreme cases $\alpha=b$ and $\alpha=0$. When $\alpha<0.5b$, the P-GW limit is recovered and trees grow as swarms, forming relatively slender tree structures, with $\gamma_{d}\approx 0.5$. In the opposite case scenario, when $\alpha\approx b$, only strong events are likely to generate aftershocks. Because the branching ratio is fixed, most of the triggered activity for $\alpha\approx b$ is due to the few stronger events, which, as consequence, are more likely to be background events. The first generation offspring of this strong event is unlikely to generate aftershocks of their own, rendering spray-like short tree sequences and star-shaped spatial structures characteristic of burst-like activity [Zaliapin \BBA Ben-Zion (\APACyear2013\APACexlab\BCnt1)]. Although the transition is not sharp in the parameter space, the empirical $\alpha/b$ values are typically close to one in seismological catalogs. Hence the separation between $\alpha/b\approx 1$ and $\alpha/b<1$ is a natural choice to define the classification between burst-like and swarm-like clusters. The GW model does not expect values $\gamma_{d}>0.5$. Significantly higher exponent values might indicate memory in the branching process, which cannot be modeled as a GW.

Overall, the ETAS model establishes a clear relationship between the topological properties and the ratio between the Gutenberg-Richter exponent $b$ and the productivity exponent $\alpha$. Swarms are resulting particularly problematic when validating the ETAS model through this relationship. On the one hand, \citeAZaliapin2013b had shown how the ETAS model parameterized with field estimations fails to predict the observed topology of swarms in southern California, where $\alpha\approx b$ and $\gamma_{T}\sim 1$. The distribution of swarm-like activity in hot areas [Zaliapin \BBA Ben-Zion (\APACyear2016)] is consistent with the power-law behavior in Eq.(8), rendering the expected value for $b/\alpha=1$. Hence, the swarm-like shape of the trees is inconsistent with the ETAS model given the fitted parameters (see Fig. 4). On the other hand, cluster reconstruction techniques are based on the homogeneous Poisson null-hypothesis [Zaliapin \BBA Ben-Zion (\APACyear2013\APACexlab\BCnt2)] and do not take into account spatio-temporal variations due to exogenous geological or anthropogenic processes, i.e. when $\mu_{0}(t,\mathbf{r})$ has a dependence in both space and time. This is the case, for example, of episodic volcanic [Roberts \BOthers. (\APACyear2016)], natural geothermal [Gaeta \BOthers. (\APACyear1998)], or human-induced [Ellsworth (\APACyear2013)] seismicity, or even tectonic seismicity in the presence of seasonal variations [Ueda \BBA Kato (\APACyear2019)]. Spatio-temporal correlations in such settings are not necessarily a consequence of a history-dependence and might disrupt the performance of cluster detection techniques, which will overestimate triggering relationships. Precisely, the random linking of uncorrelated events would lead to the generation of a P-GW process, compatible with the tree shapes reported by \citeAZaliapin2013b. In particular, geothermal systems reproducing P-GW processes might not represent actual triggering, but rather exogenous variations in the background rate. Studies focused on well confined episodes of seismic activity are advised to validate the results.

Discrepancies between model and data might be corrected by modifications of the ETAS models implementing more sophisticated field observations. Magnitude-magnitude correlations [Lippiello \BOthers. (\APACyear2008), Davidsen \BBA Green (\APACyear2011)] or depth-dependent $m$-distributions would, in general, fall outside the GW category. However, this is not a rule of thumb. In particular, recent observations suggest a simple distinction in the Gutenberg-Richter exponent for aftershocks ($b_{AS}$) and mainshocks ($b$) [Gu \BOthers. (\APACyear2013), Davidsen \BBA Baiesi (\APACyear2016), Davidsen \BOthers. (\APACyear2017)]. This specific modification falls still within the GW category. The relation between topological properties would, in that case, depend on $k_{c}$, $\alpha$, $b_{AS}$ and $b$.

Finally, any empirical study on the topological properties should account for the fact that natural catalogs have a limited spatial and temporal range. Even if the branching approach is valid, the concept of leaf and root are ill-defined in spatially confined or finite time series [Zhuang \BOthers. (\APACyear2002)]. Events misclassified as mainshocks might be actually triggered by remote or ancestral events [van der Elst (\APACyear2017)] and apparent leaves might generate aftershocks outside the observational range. The results can be especially biased for the analysis of short catalogs, considering the power-law kernels in $\Psi_{t}$ and $\Psi_{r}$ of the ETAS model.

6 Conclusions

Here, we have considered the branching model as a valid representation of aftershock sequences and other triggering processes. The reconstruction of triggering trees accounting for all event-event correlations opens new perspectives to learn about the seismogenic mechanisms behind aftershocks and improve our current forecasting techniques. We have revisited, and added to, the expected topological properties for the standard Epidemic-Type Aftershock Sequence model (ETAS), which is interpreted as a fat-tailed Galton-Watson process, extended from the Poisson Galton-Watson process (P-GW) which is recovered as a particular case. This list of properties had been proposed in previous works to be helpful to characterize and classify aftershock sequences. Specifically, this analysis serves to distinguish between swarms and bursts [Zaliapin \BBA Ben-Zion (\APACyear2013\APACexlab\BCnt1), Zaliapin \BBA Ben-Zion (\APACyear2016)].

All the singular properties of the ETAS model within the GW category derive from the power-law distribution of branching ratios $n_{b}$ leading to a similar power-law in $p_{1}(k)$. In particular, all topological properties depend only on two parameters: the average branching ratio $N_{b}$ and the ratio between exponents $\alpha/b$. Since the P-GW is recovered as a limiting case of the ETAS model we observe a transition in the distribution of tree sizes ($N_{T}$). The $N_{b}$ controls the characteristic $N_{T}$ marking the transition between two distinct power-law exponents ($\gamma_{N}^{\mathrm{low}}$ and $\gamma_{N}^{\mathrm{high}}$) , which coincide to $\gamma_{N}^{\mathrm{low}}=\gamma_{N}^{\mathrm{high}}=2$ for $\alpha=b$, and recovers the Borel distribution for $\alpha=0$. Characteristic generational depths are strongly dependent on $\alpha/b$ once the power-law tail has a significant statistical weight ($\gamma_{1}<3$). The average leaf-depth of a tree ($\langle d_{l}\rangle$) has a power-law dependence with $N_{T}$ for high branching ratios or low exponent ratios $\alpha/b$. The exponent of this power-law relation coincides with $\gamma_{d}=0.5$, typical in well defined GW processes, for $\alpha=0$ and decreases for values $b\lesssim 2\alpha$, vanishing to zero at the limit $\alpha=b$. For low branching ratios this transition in the exponent gets sharper and occurs at lower values, for example $b/\alpha\sim 0.3$ for $N_{b}=0.5$. This numerical result interprets the separation between bursts and swarms as a phenomenological observation based upon the common $\alpha\approx b$ found in nature. A regional analysis of the ratio is required to validate the exponent dependence on the $\alpha/b$ ratio. In any case, the results in the ETAS model prove that the topological structure of swarms, as bursts, can be explained as an event-event triggering processes, i.e. aftershocks, represented by one-to-one causal links.

In general, this study can be used to validate the ETAS model for the description of the triggering processes associated with tectonic and induced seismicity as well as acoustic emission experiments [Benioff (\APACyear1951), Hirata (\APACyear1987), Baró \BOthers. (\APACyear2013), Ribeiro \BOthers. (\APACyear2015), Costa \BOthers. (\APACyear2016), Davidsen \BOthers. (\APACyear2017)] and micromechanical models [Yamashita \BBA Knopoff (\APACyear1987), Dieterich (\APACyear1994), Hainzl \BOthers. (\APACyear1999), Lyakhovsky \BOthers. (\APACyear2005), Jagla \BBA Kolton (\APACyear2010), Zhang \BBA Shcherbakov (\APACyear2016), Baró \BBA Davidsen (\APACyear2018)]. A rejection of the ETAS model from the topological properties of the triggering trees might indicate more sophisticated epidemic processes, involving magnitude-magnitude correlations or depth dependencies. On the contrary, the validation of the ETAS model would set a step forward in the testing and development of new micromechanical models implementing seismogenic mechanisms of aftershocks and other event-event triggering mechanisms.

Acknowledgements.

Thanks to A. Corral, K. Wiese and I. Zaliapin for fruitful discussions. Thanks to AXA Research Fund for financial support through the project RheMechFail.