Inhomogeneous spatio-temporal epidemic-type aftershock sequence model incorporating seismicity-triggering slow slip events
Isaías Bañales, Tomoaki Nishikawa, Yoshihiro Ito, Vladimir Kostoglodov, Ekaterina Kazachkina, José Santiago
TL;DR
This paper introduces a new model to study how slow fault slips trigger earthquakes, improving understanding of seismic activity patterns.
Contribution
The paper extends the ETAS model by incorporating slow slip events to better capture background seismicity changes.
Findings
Background seismicity increases during slow slip events in Guerrero, Mexico and Boso Peninsula, Japan.
The model uses GNSS data and EM algorithm to infer seismicity changes and earthquake genealogy.
The approach allows for a more comprehensive understanding of the relationship between slow and fast earthquakes.
Abstract
Clarifying the relationship between regular earthquakes and slow fault slip is essential for understanding the mechanisms behind seismic activity. We hypothesize that the background seismic activity is partially triggered by interplate slow-slip events (SSEs). Consequently, we present an extension of the spatio-temporal epidemic-type aftershock sequence (ETAS) model, which incorporates background seismicity as a piecewise constant function over time based on recent advances in the inference of space–time inhomogeneous point processes. In this study, Global Navigation Satellite System (GNSS) data is employed to identify the occurrence periods of SSEs, thereby delineating the intervals during which changes in background seismicity may occur. Due to the technical complexity of performing inference with an inhomogeneous ETAS model, this work employs a maximum likelihood inference method…
Genes, proteins, chemicals, diseases, species, mutations and cell lines named across the full text — each resolved to its canonical identifier and authoritative record.
Click any figure to enlarge with its caption.
Figure 10
Figure 11
Figure 12
Figure 13
Figure 14
Figure 15
Figure 1
Figure 2
Figure 3
Figure 4
Figure 5
Figure 6
Figure 7
Figure 8
Figure 9
Figure 16- —https://doi.org/10.13039/501100001691Japan Society for the Promotion of Science
- —https://doi.org/10.13039/501100009037Science and Technology Research Partnership for Sustainable Development
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Taxonomy
Topicsearthquake and tectonic studies · Geological and Geophysical Studies Worldwide · Seismology and Earthquake Studies
Introduction
At the boundaries between tectonic plates, two types of spontaneous and episodic fault slip phenomena occur: fast (regular) earthquakes and slow earthquakes. These fault slip phenomena are closely related. Among slow earthquakes, those with a relatively large magnitude (approximately Mw 5 or greater) that can be detected geodetically are referred to as slow slip events (SSEs)^1^. It has been observed that SSEs often trigger small to moderate earthquakes in the Sagami Trough subduction zone in Japan^2^. Furthermore, SSEs have preceded and possibly triggered megathrust earthquakes at several subduction plate boundaries^3,4^.
Global Navigation Satellite System (GNSS) networks provide detailed daily information on crustal deformation and allow to detect and describe SSEs^5–8^. SSEs have also been discovered and studied in various regions worldwide, including, Middle America Trench^9^, the Japan Trench^10,11^, the Hikurangi subduction zone in New Zealand^12^ , and Peru^13^. In these regions, SSEs occurring at the plate boundary are thought to significantly impact seismic activity, and therefore quantifying the impact of SSEs is essential for improving the accuracy of earthquake forecasts.
The modeling and forecasting of fast-earthquake activity in a stochastic context has been widely accepted by the seismological community since the presentation of the seminal article by Ogata^14^, in which the epidemic-type aftershock sequence (ETAS) model is defined using the Hawkes process. Moreover, Zhuang et al.^15^ made an important extension to the ETAS model by the inclusion of a spatio-temporal component in the intensity function (i.e., seismicity rate). In addition, Li et al.^16^ have studied in detail the inference of the background intensity function (i.e., the background seismicity rate in seismological terms) in spatio-temporally inhomogeneous point processes, this work was a pillar in the development of the model presented in “Methodology” section.
The detection of aseismic transients and their relationship to seismicity have been extensively studied. There have been several studies of seismicity changes due to aseismic transients considering only the time domain^17–21^. For example, Okutani & Ide^22^ investigated the impacts of SSEs on seismic activity using the temporal ETAS model. They proposed a model called the boxcar model, in which the background seismicity rate increases in a boxcar-like manner during the slow slip period estimated from geodetic observations. Their approach is similar to that one presented by Mattews & Reasenberg^23^, who investigated the quiescence of microearthquakes through a temporally inhomogeneous Poisson process, using a piecewise constant function.
Nishikawa & Nishimura^24^ presented a variant of the ETAS model that explicitly links an increase in background seismicity to detected SSE. Although their research has made significant advances in the modeling and forecasting of seismic activity associated with SSEs, it does not account for the spatio-temporal changes that may occur in background seismic activity during SSE periods.
As described above, extensive research has been conducted on time-domain analysis. In this work, however, the focus is on spatio-temporal seismicity modeling.
Llenos & McGuire^25^ proposed a complex model that combined the ETAS model and the rate- and state-dependent friction seismicity model^26^ to detect seismicity rate changes induced by aseismic transients. They adopted a tricky approach to subtract the coseismically triggered seismicity rate estimated by the conventional spatio-temporal ETAS model from the total seismicity rate and related the residual seismicity rate to the rate- and state-dependent friction seismicity model.
It is important to highlight the contributions of Marsan et al.^27^ and Reverso et al.^28^ to the spatio-temporal modeling of seismicity associated with aseismic transients. Both papers propose a way to model the evolution of the seismicity using a mesh over space and time. They used the conventional spatio-temporal ETAS model throughout the entire period as a null model, and for each earthquake occurring at time \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t_i$$\end{document} in the location \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(x_i,y_i)$$\end{document} , they fit a locally elevated background intensity using earthquake records where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(t_j,x_j,y_j)$$\end{document} satisfy
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} |t_j-t_i|&<\frac{\tau }{2},\\ |x_j-x_i|&<\frac{\mathcal {L}}{2},\\ |y_j-y_i|&<\frac{\mathcal {L}}{2}, \end{aligned}$$\end{document}for all n, where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {L}$$\end{document} are parameters that control the size of the spatio-temporal window. If the locally estimated background intensity significantly differs from that of the null model, as determined by the criterion specified in each study, they replace the background intensity of the null model by the locally estimated value within the vicinity. Additionally, Reverso et al.^29^ have presented a pioneering work relating the ETAS model with SSE, following the ideas presented in^28^.
Possible improvements to the methods of Marsan et al.^27^ and Reverso et al.^28^ include the utilization of geodetic observations. In their methods, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {L}$$\end{document} are subjectively chosen (e.g., \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau$$\end{document} is set to 1 day, 40 days, or 100 days). However, particularly for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau$$\end{document} , the duration of an aseismic transient can sometimes be estimated based on geodetic data, which can then be used as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau$$\end{document} .
Furthermore, they use the small spatio-temporal windows for each earthquake one by one to estimate the local background intensity via maximum likelihood estimation, sequentially comparing it to the null model. However, this is an approach adopted for simplicity. Ideally, the background intensities of multiple windows should be varied simultaneously to estimate the set of background intensities that maximize the likelihood.
In light of the aforementioned studies, this research proposes a new modification of the spatio-temporal ETAS model that incorporates the impact of SSEs on seismic activity (“Model” section). Our model determines the periods in which the background intensity changes based on GNSS observations and simultaneously estimates the spatial distribution of the background seismicity rate for each period. This model is mathematically grounded by Li et al.^16^.
We apply the new model to Mexican earthquakes in the Middle America subduction zone and elucidate the impact of SSEs on the background seismicity within this subduction zone, which is a topic that has not been addressed from the perspective of statistical modeling to date. In addition, the model was applied to the Boso Peninsula, located in the Sagami Trough subduction zone in Japan, a region that has been previously studied and has exhibited substantial changes in the seismic activity during SSE^22,30^. Our new model will be a useful tool in the future for elucidating the characteristics of seismic activity associated with SSEs worldwide.
Methodology
Introduction to spatio-temporal ETAS Model
The spatio-temporal ETAS model is a marked branching point process for earthquake occurrences, and its behavior can be completely defined through its conditional intensity function given by
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned}&\mathbb {P}(\text {an event in }[t,t+dt]\times [x,x+dx]\times [y,y+dy]\nonumber \\&\quad \times [M,M+dM]|\mathcal {H}_t)=\lambda (t,x,y,M|\mathcal {H}_t)dtdxdydM+o(dtdxdydM), \end{aligned}$$\end{document}where M is the magnitude, (x, y) denotes the spatial coordinates, t represents the elapsed time, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {H}_t$$\end{document} denotes the space-time and magnitude occurrence history of the earthquakes up to time t^15^. In particular, based on its assumptions, the spatio-temporal ETAS model has
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \lambda (t,x,y|\mathcal {H}_t)=\mu (x,y)+\sum _{\{k:t_k<t\}} k(M_k)g(t-t_k)f(x-x_k,y-y_k|M_k), \end{aligned}$$\end{document}where the subindex k refers to the values of the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k-$$\end{document} th event considered, with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k=1,2,...,n$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(x_k,y_k)$$\end{document} are its spatial coordinates, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t_k$$\end{document} is the time of occurrence, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$M_k$$\end{document} is its magnitude. Without loss of generality, in this work it is assumed that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t_k$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k=1,2,...,n$$\end{document} are ordered in increasing order. Furthermore,
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} k(M)&=Ae^{\alpha (M-M_0)} \end{aligned}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} g(t)&=(p-1)c^{p-1}(t+c)^{-p}\mathbbm {1}(t>0)\end{aligned}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} f(x,y|M)&=\frac{1}{2\pi \sqrt{d_1d_2}e^{\alpha (M-M_0)}}\exp \Big \{-\frac{1}{2}\frac{1}{e^{\alpha (M-M_0)} }\Big (\frac{x^2}{d_1}+\frac{y^2}{d_2}\Big ) \Big \}, \end{aligned}$$\end{document}where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$M_0$$\end{document} denotes a reference magnitude, typically the minimum observed value, although in some cases a higher cutoff magnitude is adopted to reduce computational load^31^, and
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \mathbbm {1}(t>0)={\left\{ \begin{array}{ll} 1, & \text {if } t>0\\ 0, & \text {otherwise} \end{array}\right. }. \end{aligned}$$\end{document}Different authors^15,32^ have used the Gaussian probability density function in equation (5). In this work, we assume that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d_1=d_2$$\end{document} for the sake of parsimony, since the aftershock classification does not show significant differences compared to the case \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d_1 \ne d_2$$\end{document} for the Mexican data, as can be seen in Figs. 6 and 14a.
Another popular option is to use a Pareto distribution^33–35^, given by
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} f_{\text {Pareto}}(x,y|M)=\frac{q-1}{\pi \sqrt{d_1d_2} e^{\gamma (M-M_0)}}\Big (1+\frac{1}{e^{\gamma (M-M_0)}}\Big (\frac{x^2}{d_1}+\frac{y^2}{d_2}\Big )\Big )^{-q}, \end{aligned}$$\end{document}which can also be seen as a particular case of the bivariate t-distribution. The advantage of using the Pareto distribution is to avoid overestimate the background seismicity function. Nevertheless, for our Mexican data, the heavy tails of the Pareto distrubtion classify as aftershocks earthquakes that are unrealistically far from their respective mainshocks as it can be seen in Fig. 14b in Appendix App. 1.
By defining \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu (x,y)=\nu u(x,y)$$\end{document} and assuming stationarity, Zhuang et al.^15^ propose the estimator of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu$$\end{document} as
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \hat{\mu }(x,y)=\frac{1}{T} \sum _j (1-\rho _j) \frac{1}{2\pi d_j^2}\exp \Big \{-\frac{x^2+y^2}{2 d_j^2} \Big \}. \end{aligned}$$\end{document}where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d_j$$\end{document} is a bandwidth that depends on how many earthquakes are close to the event j, and 1 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-\rho _j$$\end{document} is the probability that the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$j-$$\end{document} event is an immigrant (i.e., background event).
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{\mu }$$\end{document} is fitted in an iterative two-step procedure, in the first iteration the vector \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta =(\nu ,A,\alpha ,c,p,d)$$\end{document} is fitted, in the second step the vector \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta$$\end{document} is taken as known and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{\mu }$$\end{document} is updated until the convergence of the log-likelihood
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \ell (\eta ):=\sum _{k=1}^n \log (\lambda _{\eta }(t_k,x_k,y_k|H_{t_k}))-\int _0^T \iint _S \lambda _\eta (t,x,y|H_t)dxdydt, \end{aligned}$$\end{document}is reached, where the analysis time is [0, T], and S is the analysis region.
It is worth mentioning the novel nonparametric Bayesian approaches to model \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu$$\end{document} . Ross & Kolev^32^ also assume that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu$$\end{document} fulfills \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu (x,y)=\nu u(x,y)$$\end{document} , where u is a probability density function and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nu$$\end{document} is a positive real number. This allows the use of a mixture of Dirichlet processes (MDP)^36^ as the prior for u(x, y). On the other hand, Molkenthin et al.^35^ assume that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu$$\end{document} can be written as
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \mu =\frac{v}{1+e^{-w(x,y)}}, \end{aligned}$$\end{document}where a Gaussian Process (GP) prior is used for w. The main advantage of using the MDP approach is that the function u always integrates 1, since it is a probability density function. In the case of GP, the integral can not be solved analytically; nevertheless, using the GP approach avoids the need for a finite approximation of the infinite mixture required in the MDP case.
It is important to mention that Veen & Schoenberg^37^ discuss in detail the numerical stability problems of maximizing the likelihood of the spatio-temporal ETAS model directly. They proposed using the Expectation-Maximization (EM) algorithm to improve the inference performance, this idea was one of our motivations to develop our model (“Model” section).
Extending the ideas presented by Veen & Schoenberg^37^, Fox et al.^38^ modeled the background intensity rate as a piecewise constant function, which allows a numerically efficient way to realize the approach proposed by Veen & Schoenberg^37^ with spatial inhomogeneities. The modeling of the background seismicity as a piecewise function is also used by different authors^39,40^ to infer the intensity function. Our approach (“Model” section) intends to advance the method proposed by Fox et al.^38^ by considering temporal inhomogeneities that will act as triggering effects due to SSEs.
Another advantage in following the approach presented by Veen & Schoenberg^37^ and Fox et al.^38^ is that the structure of the branching process is inferred in addition to estimating the intensity function. Therefore, we can easily distinguish which earthquakes are background events and which are aftershocks, as shown in the results obtained in “Results” section.
Since the ETAS model was introduced by^14^, the completeness of catalogs remains crucial for the correct estimation of the model. This is the reason that authors such as Seif et al.^33^ have made significant efforts to analyze the consequences of missing recorded aftershocks on the biases of the estimators. In their work, they argue that a possible solution to avoid incompleteness in aftershock sequences after a strong earthquake is to use a large cutoff magnitude. However, this can result in some seismic activity as foreshocks are not being observed and it is important to consider that the inclusion of small earthquake reduce the bias in the ETAS model estimators.
Short-term aftershock incompleteness has been extensively studied in sismology^41,42^.This is a problem of blindness, whereby small earthquakes are undetectable when the signal is saturated by a strong earthquake. Empirical relationships have been propose to model the STAI effect^43^. In the context of the ETAS algorithm, in recent years Hainzl^31^ has proposed the ETASI model, which incorporates the STAI phenomena into the ETAS model by adding a temporary dependence on the number of expected events. The ETASI model has been recently extended by Asayesh et al.^44^, where the spatial kernel of the ETASI model was modified to incorporate information from the stress scalars.
Although STAI is an important phenomenon to take into account, the cutoff magnitudes used in this study are not low enough to appreciate this phenomenon as it can be seen in Fig. 12 in Appendix App. 1. Therefore, we chose to keep b-values and cutoff magnitudes constant throughout the entire period.
We also wish to point out that, in order to model spatial variation in aftershock activity parameters, authors such as Ogata^45^ and Ueda et al.^46^ extended the spatio-temporal ETAS model^15^, allowing parameters such as p and the productivity K to become spatially varying functions. However, the issue of identifiability in their approach required the imposition of a smoothness penalization on the functions p, and K. Nevertheless, extending \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu$$\end{document} , K, and p to functions of (x, y, t) leads to identifiability problems, since changes in seismicity induced by strong earthquakes can occur abruptly in both space and time, rendering smoothness penalization ineffective. For this reason, our work only considers temporal inhomogeneity through \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu$$\end{document} .
Model
Following the idea of Nishikawa and Nishimura^24^, we model the inhomogeneities in the ETAS model through changes in the background seismic activity. The reason for this is that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu$$\end{document} is usually related to tectonic loading and the relative velocity of plate motion^47^. If we regard slow slip as an increase in tectonic loading or interplate slip rate, it is natural to model it through changes in the background seismic activity. However, we cannot deny the impact of SSEs on aftershock activity. Addressing this issue is an important direction for future work.
Since the aim of this work is to model the triggering effect of SSEs, the assumption of stationarity may not be realistic. In addition, Veen & Schoenberg^37^ have discussed that the approach introduced by Zhuang et al.^15^ is numerically expensive and unstable because it is necessary to optimize (7).
To describe the triggering effect due to SSEs, this work proposes to use
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \lambda (t,x,y|\mathcal {H}_t)=\mu (x,y,t)+\sum _{\{k:t_k<t\}} k(M_k)g(t-t_k)f(x-x_k,y-y_k|M_k), \end{aligned}$$\end{document}which is an extension of (2) because it allows a time dependent \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu$$\end{document} function. However, due to the complexity of working with an arbitrary form of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu$$\end{document} , in this study, it is defined as
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \mu (x,y,t)= \sum _{i=0}^m \mu ^i(x,y)\mathbbm {1}(t\in S_i,E_i), \end{aligned}$$\end{document}where m is the total number of SSEs and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$S_i,E_i$$\end{document} are the start and end times of the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$i-$$\end{document} th SSE, with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$i=1,...,m$$\end{document} . In the case of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$i=0$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$S_i=0,E_i=T$$\end{document} , i.e. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu ^0$$\end{document} represents the background seismicity over the entire period, regardless of the presence of an SSE, and the remaining \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu ^i$$\end{document} with i=1,2,..., m represent the increase of the background seismicity with respect to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu ^0$$\end{document} . Thus, the intensity is given by
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \lambda (t,x,y|\mathcal {H}_t)=\sum _{i=0}^m \mu ^i(x,y)\mathbbm {1}(t\in S_i,E_i) +\sum _{\{k:t_k<t\}} k(M_k)g(t-t_k)f(x-x_k,y-y_k|M_k). \end{aligned}$$\end{document}As presented by Veen & Schoenberg^37^ and discussed in detail by others^48,49^, the Hawkes process could be defined as a marked Poisson cluster process where there are two kind of events, immigrants (background events) and offspring (aftershocks), which allows the use of the EM algorithm to maximize the likelihood.
McLachlan & Krishnan^50^ discussed the EM algorithm in detail and showed that it is useful when there is missing information which, if known, makes the maximization of the complete likelihood easier. The EM algorithm is an iterative algorithm that maximizes the likelihood using a two-step procedure, the first step (Expectation) replaces the unknown information by the expected one with the help of an initial value for all the parameters, and the second step (Maximization) uses the expected values obtained in the previous step to optimize the likelihood defined by the augmented information, these two steps are repeated until the convergence of the augmented likelihood.
The development of the model in this section is mathematically based on that proposed by Li et al.^16^, but the definition of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda (t,x,y|\mathcal {H}_t)$$\end{document} by them differs from that presented in (9). In their work they assumed
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \mu (x,y,t)=\alpha u(x,y)v(t), \end{aligned}$$\end{document}where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha$$\end{document} is in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathbb {R}^+$$\end{document} , and u and v are positive functions. Since SSEs are spontaneous large fault slip events with variable slip evolution, the spatial and temporal contributions to background seismicity cannot be assumed to form a product over the entire period. Therefore, we preferred to use the expression in (8).
It is also important to mention that they are not working with a marked point process and that their intensity function does not incorporate information about event magnitudes. However, they use the kernel of a gamma distribution for the decay in the aftershock activity over the time with parameter \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha =1$$\end{document} :
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} g^L(t)=\beta e^{-\beta t} \end{aligned}$$\end{document}where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta$$\end{document} is in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathbb {R}^+$$\end{document} . Their expression differs from (4), proposed by Zhuang et al.^15^, which is based on Omori’s law^51^ and is also adopted in our study.
From the observed data, it is not known whether an earthquake is a background event or an aftershock, and if it is a background event, it is also not known whether it comes from the process with intensity \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu ^0$$\end{document} or from \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu ^i$$\end{document} when the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$i-$$\end{document} th SSE occurs. Therefore, the following random variables are defined:
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \chi ^s_{ii}&={\left\{ \begin{array}{ll} 1, & \text {if earthquake }i\text { is a background event, and it is produced by } \mu ^s\\ 0, & \text {otherwise} \end{array}\right. }\end{aligned}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \chi _{ij}&={\left\{ \begin{array}{ll} 1, & \text {if earthquake }i\text { is an aftershock of }j\\ 0, & \text {otherwise} \end{array}\right. }, \end{aligned}$$\end{document}where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s=0,...,m$$\end{document} .
The random variables defined in (10) and (11) are an extension of those presented in the supplementary material by Fox et al.^38^, where only one \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\chi _{ii}$$\end{document} is introduced for the background intensity over the entire period, while here the triggering effect of multiple SSEs is modeled in a similar way to the model by Li et al.^16^. In the present work, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu ^i(x,y)$$\end{document} is defined as a piecewise constant function for all i:
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \mu ^i(x,y)=\sum _{u=1}^{n_y} \sum _{v=1}^{n_x} \mu _{uv}^i\mathbbm {1}((x,y)\in D_{uv}), \end{aligned}$$\end{document}where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$D_{kl}=((u-1)\Delta x, u \Delta x)\times ((v-1)\Delta y, v \Delta y)$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta x$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta y$$\end{document} are the step size in the partition of the x and y axes, also \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n_x$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n_y$$\end{document} are the number of grids for each axis. The advantages of defining \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu ^i(x,y)$$\end{document} in this way is to facilitate the integral over the space in (17) to recover a closed expression in the maximization step, as it can be seen in (19). The grid applied in this study is shown in Fig. 1.
The total amount of parameters to be estimated by (12) are \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n_y n_x$$\end{document} for each \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu ^i$$\end{document} , with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$i\in 0,1,...m$$\end{document} . Furthermore, because of (3) to (5), we have the five aftershock parameters \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(A,\alpha ,c,p,d)$$\end{document} to be estimated. Additionally, we must estimate the whole branch structure given by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n(n-1)/2$$\end{document} elements in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\chi _{ij}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sum _{i=1}^m n_i$$\end{document} elements in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\chi ^s_{ii}$$\end{document} , where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n_i$$\end{document} is the total number of earthquakes that occurred during \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(S_i,E_i)$$\end{document} .
In this work, the array that contains all the parameters of the previous paragraph is denoted by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta$$\end{document} , i.e.
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \theta =( \{\chi _{ii}\}_{i\in I},\{\chi _{ij}\}_{ij\in I\times I},A,\alpha ,c,p,d ). \end{aligned}$$\end{document}where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$I=\{1,2,...,n\}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\times$$\end{document} denotes the cartesian product. It is important to note that the enormous amount of parameters will produce a not well-posed problem, the optimization of the likelihood has been constrained, following^52^ according to the expression (20). This study assumes that the parameters controlling the aftershock intensity \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(A,\alpha ,c,p,d)$$\end{document} do not change over time.
If the branch structure were known, the log-likelihood of the complete information ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\ell _c$$\end{document} ) would be defined as
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \ell _c(\theta )=\ell ^*_\text {O}(\theta )+\ell ^*_\text {I}(\theta ), \end{aligned}$$\end{document}where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\ell _\text {O}(\theta )$$\end{document} is the log-likelihood due to the offspring, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\ell _\text {I}(\theta )$$\end{document} is the log-likelihood due to immigrants. They are given by
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \ell ^*_\text {I}(\theta )&= \sum _{s=0}^m \Big [ \sum _{i=1}^n \chi _{ii}^s\log (\mu _\theta ^s(x_i,y_i)\mathbbm {1}(t_i\in (S_s,E_s))\nonumber \\&\quad -\int _0^T \iint _S \mu _\theta ^s(x_i,y_i)\mathbbm {1}(t\in (S_s,E_s)) dxdydt \Big ] \end{aligned}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \ell ^*_\text {O}(\theta )&= \sum _{j=1}^n \Big [ \sum _{i>j} \chi _{ij} \log \Big (k_{\theta }(M_j)g_{\theta }(t_i-t_j)f_{\theta }(x_i-x_j,y_i-y_j|M_j)\Big ) \nonumber \\&\quad - \int _{t_j}^T \iint _S k_{\theta }(M_j)g_{\theta }(t-t_j)f_{\theta }(x-x_j,y-y_j|M_j) dxdydt \Big ], \end{aligned}$$\end{document}where the subindex \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta$$\end{document} in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k_{\theta },g_{\theta },f_{\theta }$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu _\theta ^s$$\end{document} means that the functions k, g and f are defined using the parameters \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(A,\alpha ,c,p,d)$$\end{document} given by the array \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta$$\end{document} where the log-likelihood is evaluated.
It is important to note that (10) and (11) are not independent random variables, the sum of all entries in the vector \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(\chi _{ii}^0,...,\chi _{ii}^m,\chi _{i,1},...,\chi _{i,i-1})$$\end{document} is always 1 for all \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$i=1,...,n$$\end{document} , i.e., it is a vector with a multinomial distribution. This means that if we have a set of parameters \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta _r$$\end{document} , the expected value of the entries for all i in 1, 2, ..., n is given by
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \mathbb {E}(\chi _{ii}^s|\theta _r)&=\hat{p}_{ii}^s=\frac{\mu ^s(x_i,y_i)}{\lambda _{\theta _0(t_i,x_i,y_i|H_{t_i}})} \mathbbm {1}(t_i\in S_s,E_s) \quad \text {with } s=0,1,...,m \end{aligned}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \mathbb {E}(\chi _{ij}|\theta _r)&=\hat{p}_{ij}=\frac{k(M_j)g(t_i-t_j)f(x_i-x_j,y_i-y_j|M_j)}{\lambda _{\theta _0(t_i,x_i,y_i|H_{t_i}})} \quad \text {with } j=1,...,i-1. \end{aligned}$$\end{document}The equations (15) and (16) correspond to the expectation step in the EM algorithm. The next step is to replace the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\chi$$\end{document} values in (13) and (14) by their expected values. Using \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta _r$$\end{document} as the initial parameter value,
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \ell _\text {I}(\theta _r)&= \sum _{s=0}^m \Big [ \sum _{i=1}^n \hat{p}_{ii}^s\log (\mu _{\theta _r}^s(x_i,y_i)\mathbbm {1}(t_i\in (S_s,E_s))\nonumber \\&\quad -\int _0^T \iint _S \mu _{\theta _r}^s(x_i,y_i)\mathbbm {1}(t\in (S_s,E_s)) dxdydt \Big ], \end{aligned}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \ell _\text {O}(\theta _r)&= \sum _{j=1}^n \Big [ \sum _{i>j} \hat{p}_{ij} \log \Big (k_{\theta _r}(M_j)g_{\theta _r}(t_i-t_j)f_{\theta _r}(x_i-x_j,y_i-y_j|M_j)\Big ) \nonumber \\&\quad - \int _{t_j}^T \iint _S k_{\theta _r}(M_j)g_{\theta _r}(t-t_j)f_{\theta _r}(x-x_j,y-y_j|M_j) dxdydt \Big ]. \end{aligned}$$\end{document}It is important to note that since the real branch structure (genealogy) of the earthquakes is unobservable, the equations (13) and (14) cannot be evaluated. However, when the real values are replaced by their expected values as in (17) and (18), they can be computed.
Another advantage of the EM approach is that all parameters concerning \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu ^i$$\end{document} are only in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(17)$$\end{document} , while the parameters \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(A,\alpha ,c,p,d)$$\end{document} are only in (18), therefore the optimization can be done separately. To solve (18) it is important to note that
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \int _{t_j}^T \iint _S k_{\theta _r}(M_j)g_{\theta _r}(t-t_j)f_{\theta _r}(x-x_j,y-y_j|M_j) dxdydt \end{aligned}$$\end{document}can be rewritten as
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} k_{\theta _r}(M_j) \int _{t_j}^T g_{\theta _r}(t-t_j)dt \iint _S f_{\theta _r}(x-x_j,y-y_j|M_j) dxdy. \end{aligned}$$\end{document}The utility of taking the expressions in (4) and (5) is that the integrals can be solved easily, since
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \int _{t_j}^T g_{\theta _r}(t-t_j)dt= 1-c^{p-1}(T-t_j+c)^{-p+1} \end{aligned}$$\end{document}and the spatial integral
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \iint _S f_{\theta _r}(x-x_j,y-y_j|M_j) dxdy, \end{aligned}$$\end{document}which corresponds to the double integral of a bivariate Gaussian distribution, which can be solved efficiently using Monte Carlo methods. As an additional note,^53^ suggests the following approximation with different k, g and f functions
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \int _{t_j}^T \iint _S k_{\theta _r}(M_j)g_{\theta _r}(t-t_j)f_{\theta _r}(x-x_j,y-y_j|M_j) dxdydt\approx k_{\theta _r}(M_j), \end{aligned}$$\end{document}this approximation holds if almost all aftershock activity occurs in the observed space and time, which is a very easy assumption to satisfy and it can lead to produce important reductions in computing time. We mention this approach because it has been used by different authors^32,38^, nevertheless, the implementation of algorithm 1 in this work does not use it because, as mentioned in^54^, it can generate inaccuracies in the estimation process.
Regarding (17), it is easy to find that the solution of
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned}&\frac{{\partial }}{{\partial } \mu ^j_{uv}} \sum _{s=0}^m \Big [ \sum _{k=1}^n \hat{p}_{ii}^s\log (\mu _{\theta _r}^s(x_k,y_k)\mathbbm {1}(t_k\in (S_s,E_s))\\&\quad -\int _0^T \iint _S \mu _{\theta _r}^s(x_k,y_k)\mathbbm {1}(t\in (S_s,E_s)) dxdydt \Big ] =0 \end{aligned}$$\end{document}is
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \hat{\mu }^j_{uv}=\frac{\sum _{k=1}^n \hat{p}^j_{ii} \mathbbm {1}((x_k,y_k)\in D_{uv}) \mathbbm {1}(t_k\in (S_j,E_j)) }{(E_j-S_j)\Delta x \Delta y}. \end{aligned}$$\end{document}The main advantage of using the piecewise expression of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu ^i$$\end{document} is to retrieve a closed expression in (19), which enables a fast updating of the parameters. However, the main weakness of this approach is that only the earthquakes inside of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$D_{kl}$$\end{document} update the value of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu ^i_{kl}$$\end{document} , and a smooth estimation of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu ^i$$\end{document} is not achieved.
To address the not well-conditioned problem we restrict the optimization problem, the parameter A must be less than 1, i.e. we are assuming that the smallest earthquake have less than one expected aftershock. Regarding the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha$$\end{document} value we followed the idea of Ross & Kolev^32^ of use the Helmstetter et al.^52^ inequalities.
Assuming the Gutenberg-Richter law^55^, the random variable ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$M-M_0$$\end{document} ) is distributed as an exponential random variable^56^ with rate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$b\log (10)$$\end{document} , where b is the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$b-$$\end{document} value of the Gutenberg-Richter law. Then, we can calculate the average number of offspring created per event, which is defined by
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} r:=\int _{M_0}^\infty Ae^{\alpha (M-M_0) } b\log (10)e^{-b\log (10)(M-m_0)}dM. \end{aligned}$$\end{document}In order to guarantee that r is finite and it is less than one, the following conditions must be fulfilled^52^
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \alpha<b\log (10), \quad \frac{Ab\log (10)}{b\log (10)-\alpha }<1. \end{aligned}$$\end{document}Consequently, in this work we allow only parameters of A and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha$$\end{document} that satisfy the inequalities in (20).
During this study, we assumed a constant b-value over time and across the entire study region. This value was estimated using the maximum likelihood method^57^. Tests were also performed in Mexico using the b-positive estimator^42^ and the novel estimator presented by Lippiello & Petrillo^58^ called b-more-positive; however, the ETAS model estimators did not change significantly. Furthermore, Table 1 in Appendix App. 1 reports estimations assuming b-value lower and higher than those obtained via MLE, b-positive method and b-more-positive and the parameters did not show any important change with respect to those presented in “Mexico” section.
Although the choice of the MLE estimator or the family of b-positive estimators does not generate significant changes in our ETAS model estimators, we note that Nandan et al.^59^ introduced a variant of the ETAS model that allows different b-values to be considered. They speculate that the magnitude of the mainshock may affect the b-value for its triggered earthquakes. Also, Ito & Kaneko^60^, who, using continuum models of fully dynamic earthquake cycles with fault frictional heterogeneities, have observed that the b-value of foreshocks decreases with time prior to the mainshock.
Furthermore, as noted by^61^, the b-value also exhibits a spatiotemporal distribution, which we plan to analyze in detail in future work.
For c and p parameters, we restrict the acceptable values to be between [0, 5] and [1, 2]. We are being very flexible in our boundaries, since Holschneider et al.^62^ have observed that c and p are weakly identifiable.
In the case of d, with units deg \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^2$$\end{document} , we restrict the optimization algorithm to be in (0,1), the idea of take 1 as the upper limit is be noninformative.
The maximization step in the EM algorithm consists in solving
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \theta _{r+1}=\underset{\theta }{\operatorname {argmax}}\quad \ell _\text {I}(\theta _r) +\ell _\text {O}(\theta _r) \end{aligned}$$\end{document}Finally, by replacing r by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${r+1}$$\end{document} the Expectation and Maximization steps must be iterated until \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$|\theta _r-\theta _{r+1}|<\varepsilon$$\end{document} , where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\epsilon >0$$\end{document} is the convergence criterion. In this work, the optimization was done using the implementation of the L-BFGS-B algorithm in python using the implementation in scipy^63^.
As a summary of the EM algorithm that is used in “Results” section to detect the earthquake triggering effect of SSEs, the pseudo code is presented in the algorithm 1.
Algorithm 1Expectation maximization algorithm of inhomogeneous spatio-temporal ETAS.
The initial values of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{p}^s_{ii}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{p}_{ij}$$\end{document} used in Algorithm 1 are defined by the following matrices
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} P^O&=\begin{pmatrix} 0 & 0 & 0& 0& \hdots & 0\\ \frac{1}{2} & 0 & 0 & 0& \hdots & 0\\ \frac{1}{3} & \frac{1}{3} & 0& 0& \hdots & 0\\ \frac{1}{4} & \frac{1}{4} & \frac{1}{4}& 0& \hdots & 0\\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ \frac{1}{n} & \frac{1}{n} & \frac{1}{n} & \frac{1}{n} & \hdots & 0 \end{pmatrix}, \\ \quad P^I&= \begin{pmatrix} \frac{\mathbbm {1}(t_1 \in (S_0,E_0))}{\sum _{i=0}^m\mathbbm {1}(t_1 \in (S_i,E_i))} & \frac{\mathbbm {1}(t_1 \in (S_1,E_1))}{\sum _{i=0}^m\mathbbm {1}(t_1 \in (S_i,E_i))} & \hdots & \frac{\mathbbm {1}(t_1 \in (S_m,E_m))}{\sum _{i=0}^m\mathbbm {1}(t_1 \in (S_i,E_i))}\\ \frac{\mathbbm {1}(t_2 \in (S_0,E_0))}{2\sum _{i=0}^m\mathbbm {1}(t_2 \in (S_i,E_i))} & \frac{\mathbbm {1}(t_2 \in (S_2,E_2))}{2\sum _{i=0}^m\mathbbm {1}(t_1 \in (S_i,E_i))} & \hdots & \frac{\mathbbm {1}(t_2 \in (S_m,E_m))}{2\sum _{i=0}^m\mathbbm {1}(t_2 \in (S_i,E_i))}\\ \vdots & \vdots & \ddots & \vdots \\ \frac{\mathbbm {1}(t_n \in (S_0,E_0))}{n\sum _{i=0}^m\mathbbm {1}(t_n \in (S_i,E_i))} & \frac{\mathbbm {1}(t_n \in (S_2,E_2))}{n\sum _{i=0}^m\mathbbm {1}(t_n \in (S_i,E_i))} & \hdots & \frac{\mathbbm {1}(t_n \in (S_m,E_m))}{n\sum _{i=0}^m\mathbbm {1}(t_n \in (S_i,E_i))}, \end{pmatrix}, \end{aligned}$$\end{document}taking \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{p}^s_{ii}$$\end{document} as the entry (i, s) of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P^I$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{p}_{ij}$$\end{document} as the entry (i, j) of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P^O$$\end{document} . The idea of the above matrices is to try to reflect the unknowns about the branching structure, starting with uniform distributions between being an aftershock or a background event, and uniformity among \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu ^s$$\end{document} that could generate the event if it is a background event.
All the codes used to fit the above model and generate the images in this work are available in https://github.com/isaiasmanuel/ETAS. The code used to perform Algorithm 1 is available in EM2.py. Also, the figures in this study were produced using the code GNSSETAS_Figures.py.
Hypothesis testing
To verify if there is an improvement of our model, i.e., by adding \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu ^i$$\end{document} for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$i = 1, 2,..., m$$\end{document} with respect to the model which only has \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu _0$$\end{document} , a hypothesis test must be performed to analyze the significance of the results. Due to the complexity of the model presented in Algorithm 1, it is not straightforward to derive the theoretical joint distribution of the estimators and determine the rejection regions. For this reason, we propose a Likelihood Ratio Test (LRT) based on parametric bootstrapping^64^.
Since bootstrap is a resampling approach that requires simulating synthetic earthquake catalogs, we will simulate the synthetic catalog in a similar way to Ross & Kolev^32^, based on the cluster process representation of the Hawkes process^65^, the advantage of this approach is avoiding the methodology of thining as in Fox et al.^38^ that requires more computational time. The simulation of synthetic catalogs is as follows
- We simulate our synthetic mainshocks for each square \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$D_{uv}$$\end{document} defined in the equation (12), and for each interval \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$i=1,2,...,m$$\end{document} , from a Poisson distribution with mean \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu ^i_{uv}(E_i-S_i)\Delta x\Delta y$$\end{document} , the occurrence time of this event is simulated from a uniform distribution in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(S_i,E_i)$$\end{document} and the magnitude is simulated from the Gutenberg-Richter law with the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$b-$$\end{document} value estimated via maximum likelihood from the real catalog^57^.
- For each simulated earthquake j that it is direct offspring have not been calculated previously, we simulate the number of direct offspring from a Poisson distribution with mean \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Ke^{\alpha (M_j-M_0)}$$\end{document} , and the offspring locations and times are simulated directly from random variables with densities given by equations (4) using \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t=t_j$$\end{document} and (5) using \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(x_j,y_j)$$\end{document} . If the new earthquakes have times larger than \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$E_0$$\end{document} or if they are outside of our study area, they are discarded.
- Repeat the previous step until no new offspring is generated. The code to simulate synthetic earthquake catalogs is available in Hypothesis_testing.py.
The hypothesis to test in the LRT are
- \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$H_0$$\end{document} : \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta$$\end{document} has \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu ^i$$\end{document} for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$i=1,2,..$$\end{document} equal to 0
- \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$H_1$$\end{document} : \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta$$\end{document} has at least one \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu ^i$$\end{document} for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$i=1,2,..$$\end{document} different to 0 From the real catalog \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\{x_i,y_i,t_i,m_i\}_{i=1}^n$$\end{document} we can obtain the maximum likelihood estimators (MLE) of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta$$\end{document} using the complete model and the model that only has \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu ^0$$\end{document} , which will be denoted by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{\theta }$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tilde{\theta }$$\end{document} respectively, then our test statistic is defined by
where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L_1$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L_0$$\end{document} denote the likelihoods of the complete and reduced model. In order to define the rejection region for our hypothesis testing, we will generate B synthetic catalogs \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\{x_i,y_i,t_i,m_i\}_{i=1}^{n_b}$$\end{document} , with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$b=1,...,B$$\end{document} , using \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tilde{\theta }$$\end{document} from the reduced model. Once the synthetic catalogs are obtained, we estimate the MLE for both models to each synthetic catalog to obtain a pair \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{\theta }_b$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tilde{\theta }_b$$\end{document} , which allows to define
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} T_{LR,b}= \frac{L_{1,b}(\hat{\theta }_b)}{L_{0,b}(\tilde{\theta }_b)}, \end{aligned}$$\end{document}where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L_{0,b},L_{1,b}$$\end{document} denote the Likelihood functions of the reduced and complete model using the synthetic catalogue b. It is important to note that the idea of the bootstrap methodology is to obtain an approximation of the distribution of our estimators through empirical distributions obtained by resampling^64^.
Finally we will reject \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$H_0$$\end{document} in favor of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$H_1$$\end{document} at a confidence level \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma$$\end{document} if the following expression is fulfilled
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \frac{\#\{T_{LR,b}:T_{LR,b}>T_{LR}\}}{B}\le \gamma , \end{aligned}$$\end{document}where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\#$$\end{document} denotes the cardinality of the set.
While the LRT allows us to determine whether there is evidence that the full model is preferable to the reduced model, it is not informative in terms of seeing which \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu ^i_ {uv}$$\end{document} are contributing to that conclusion. Taking this into account and the fact that we already calculate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{\theta }_b$$\end{document} for all synthetic catalogs, we can examine marginally which \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu ^i_{uv}$$\end{document} from \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{\theta }$$\end{document} , denoted by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{\mu }^i_{uv}$$\end{document} , are significantly greater than the corresponding \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu ^i_{uv}$$\end{document} from \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{\theta }_b$$\end{document} , denoted by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu ^{i,b}_{uv}$$\end{document} . i.e., we can define
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} T_{uv}=\frac{\#\{\mu ^{i,b}_{uv}:\mu ^{i,b}_{uv}\ge \hat{\mu }^i_{uv}\}}{B}, \end{aligned}$$\end{document}and conclude that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{\mu }^i_{uv}$$\end{document} is significantly higher than \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu ^i_{uv}$$\end{document} under the reduced model at level \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma$$\end{document} if \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T_{uv}\le \gamma$$\end{document} .
Note that the value \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu _{uv}^i$$\end{document} for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$i=1,2,...$$\end{document} can be 0 if no earthquake was observed in the square \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$D_{uv}$$\end{document} during the interval \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(S_i,E_i)$$\end{document} or if the observed earthquakes were classified as offspring of an earthquake that ocurred prior to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(S_i,E_i)$$\end{document} .
In the “Results” section this hypothesis testing methodology is applied for records in Guerrero, Mexico and the Boso Peninsula, Japan.
Additionally, since we are obtaining synthetic catalogs with a sample size \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n_b$$\end{document} , we select the size of our mesh in equation 12 as the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n_x=n_y$$\end{document} that minimize
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \frac{1}{B} \sum _{b=1}^B (n-n_b)^2. \end{aligned}$$\end{document}In other words, the mesh size is selected through minimize the mean squared error of the number of earthquakes. It is worthy to mention that Asim et al.^66^ have proposed the idea of use a multiresolution mesh to improve the computational time, nevertheless due the size of our region a non multiresolution mesh is still computationally affordable.
Data
Mexico
In Fig. 1, all the epicenters in the study are shown in blue, on a mesh with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n_x=n_y=14$$\end{document} , following the notation in equation (12). The earthquake catalog was obtained from the Advanced National Seismic System (ANSS) ^67^, with dates from 2000/01/01 to 2016/12/31, and magnitudes greater than or equal to 4.3. Only the earthquakes related to plate subduction (i.e., on the landward side of the trench) were considered in this study. The magnitude completeness of the ANSS catalog in the study region was assessed using the Maximum Curvature Method^68^, and the magnitude of completeness (Mc) was found to be M4.0. However, as this method is known to systematically underestimate Mc by a fraction of a unit^69^, M4.3 was used as the threshold \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$M_0$$\end{document} in equations (3) and (5) in the present work.
On September 8, 2017, an earthquake of M8.2 occurred in Mexico and some studies, such as^70^, have suggested that it may have broken the entire subducted Cocos lithosphere and significantly altered the seismicity in the subduction system. For this reason, we limit our study period to the end of 2016.Fig. 1. Seismicity data from Guerrero, Mexico. Earthquakes epicenters^67^ are shown in blue. Green circles are magnitude scale. The red squares define the grid used in Algorithm 1. The pink line denotes the Middle America Trench^71^. The gray lines show the political division of Mexico^72^, the blue, orange, fuchsia contours are the 4 cm slip contours curves digitalized from^73^ of the 2001-2002, 2006 and 2009-2010 SSE respectively, the black curve is the 15 cm slip contour of the 2014 SSE digitalized from^4^. The yellow cross correspond to the location of Cayaco GNSS station. The contour lines in purple are the slab model from^74^. We also present the slip distribution contours at 1-meter intervals from 1 to 6 meters for the earthquakes with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$M>7$$\end{document} , using data from^67^. These earthquakes occurred on 2003/01/22, 2012/03/20, and 2014/04/18, with epicenters \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-104.1040\hbox {E}^{\circ }$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$18.770\hbox {N}^{\circ }$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-98.2310\hbox {E}^{\circ }$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$16.493\hbox {N}^{\circ }$$\end{document} , and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-100.9723\hbox {E}^{\circ }$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$17.397\hbox {N}^{\circ }$$\end{document} and magnitudes 7.6, 7.4 and 7.2, respectively.
In Fig. 2, the North-South component of the GNSS time series of the Cayaco (CAYA) continuous GPS monitoring station is presented. CAYA is located in Guerrero state, Mexico, at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-100.2672\hbox {E}^{\circ },17.0485\hbox {N}^{\circ }$$\end{document} . The duration of SSEs can be defined as the periods of southward displacements at this station^9^, and the SSEs appear to have a periodicity of approximately four years. This GNSS data is available in https://github.com/isaiasmanuel/ETAS as caya_ns.
To estimate the start and end times ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$S_i$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$E_i$$\end{document} ) of the SSEs, b-splines^75^ are fitted to the raw data. The duration of the SSEs is defined by the periods when the spline has a negative slope, as shown in Fig. 2. The data for this work was provided by the Servicio Sismológico Nacional (SSN) of Mexico and by the Department of Seismology, IGEF UNAM.Fig. 2. Global navigation satellite system data at Cayaco station. North-to-south component at Cayaco station (blue). Fitted b-spline polynomial (orange). The periods when the data slope is negative is indicated by vertical red lines, and the purple vertical lines are when an earthquake with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$M>7$$\end{document} ocurred.
The SSE slip contour curves with 4 cm of slip spacing (from 2 cm to 18 cm) of the first 3 SSEs (the 2001-2002, 2006, and 2009-2010 SSEs) have been digitalized from Radiguet et al.^73^. The 15 cm slip contour of the 2014 SSE was digitalized from^4^ and is presented in Figure 1. The expected result of the algorithm 1 is that if an earthquake triggering effect of SSEs exists, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu ^i$$\end{document} with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$i=1,2,3,4$$\end{document} should have an increase near the SSE slip contours.
Japan
To test the flexibility and capability of our model, we additionally explored the Boso Peninsula, located in the Sagami Trough subduction zone in Japan. The seismicity-triggering effect of SSEs and swarm activity produced during SSE periods in the Boso Peninsula have been previously studied by Okutani & Ide^22^ and Fukuda^30^. For this analysis, we used data from the Japanese Meteorological Data^76^ from 2001/01/01 to 2009/01/01 and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$M_0=3$$\end{document} . Accoding to the geodetic analysis by Fukuda^30^, during this period, two SSE ocurred from 2002/10/01 to 2002/10/19 and from 2007/08/12 to 2008/08/25. These periods were used to define \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu ^1$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu ^2$$\end{document} , respectively.
Figure 3 shows the seismicity during our study period, and the mesh used in Algorithm 1, in this example \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n_x=n_y=6$$\end{document} . The daily slip rate distributions of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$4\frac{\text {m}}{\text {year}}$$\end{document} for the 2002 and 2007 SSEs are also presented. The digitized data correspond to the slip rates on 2022/10/08 and 2007/08/16 from^30^, which represent the days with the highest slip rates.Fig. 3 Seismicity data from the Boso Peninsula, Japan. The colors are the same as in Fig. 1, except for the SSEs. The fuchsia and red contours denote the daily slip rate distributions of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$4\frac{\text {m}}{\text {year}}$$\end{document} for the 2002 and 2007 SSEs, respectively, digitized from^30^.
It is important to note that, while the SSEs in Mexico occurred over several months, those in the Boso Peninsula lasted only several weeks. This substantial difference in time scale is a key reason for analyzing the performance of our model in both regions.
Results
Mexico
Using \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu ^1$$\end{document} , we observe a substantial increase in the background seismicity of 0.18 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{\text {Events}}{\text {deg}^2 \text {day}}$$\end{document} , compared to the stationary background rate of 0.01 given by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu ^0$$\end{document} , during the 2001-2002 SSE in the square with a vertex at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-100.11\hbox {E}^{\circ }$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$17.18\hbox {N}^{\circ }$$\end{document} (Fig. 4). It is important to note that this area is located close to the regions with the largest slip gradient in the slip contours of Fig. 4.The significant influence of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu ^1$$\end{document} on the background seismicity is further supported by Fig. 5, which highlights the squares with a significant \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T_{uv}^i$$\end{document} at the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$90\%$$\end{document} confidence level using \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$B=100$$\end{document} . Additionally, the LRT rejects the null hypothesis (p-value = 0.09) in favor of our model.
Since SSE slips release stored elastic energy, the regions with the largest SSE slips have experienced a reduction in cumulative stress. In contrast, stress increases just outside the margins of the SSE patches, where an increase in background seismicity can be expected due to the triggering effects. This is consistent with our estimation in Fig. 4.
Regarding \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu ^2$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu ^3$$\end{document} , we observe values greater than zero next to the epicenter of the 2003/01/22 earthquake and in the margins of the 2006 and 2009-2010 SSEs. Nevertheless, these values are too small and dispersed to conclude that a triggering effect could be visible. Furthermore, as shown in Fig. 5, the observed values in nearly all the squares are not statistically significant.
Radiguet et al.^4^ have discussed that the 2014/04/18 M7.3 Papanoa earthquake was triggered by the 2014 SSE. Consistently, our model shows a significant increase in the background seismicity rate during the SSE period in the square that contains the epicenter, as well as in the square with vertex at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-100.73\hbox {E}^{\circ }$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$16.89\hbox {N}^{\circ }$$\end{document} where, according to^4^, the SSE caused a slip of approximately 0.05 m.
We also found a significant value of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu ^4$$\end{document} in squares with vertices at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-98.88\hbox {E}^{\circ }$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$16.03\hbox {N}^{\circ }$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-98.26\hbox {E}^{\circ }$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$16.89\hbox {N}^{\circ }$$\end{document} which are next to the epicenter of the 2012/03/20 earthquake. Although we have focused this work on discussing changes in seismicity during periods of SSEs, it is important in future work to consider the background seismicity variations around large earthquakes to account for possible tectonic stress rearrangements.
While \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu ^1$$\end{document} reached the highest values among all \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu ^i$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu ^2$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu ^3$$\end{document} exhibited much lower levels of activity. This may be related to differences in the source characteristics of each SSE (i.e., moment magnitude, duration, and slip distribution).
In terms of moment magnitude, the 2001–2002 SSE was the largest (Mw 7.65), while the 2006, 2009–2010, and 2014 SSEs had smaller magnitudes of Mw 7.49, Mw 7.54, and Mw 7.60, respectively, according to Radiguet et al.^4,73^. The duration of the 2001–2002 SSE was 436 days, while the other three SSEs lasted 314, 510, and 378 days, respectively.
In particular, the 2009-2010 SSE exhibited relatively small moment magnitude and longest duration, indicating a lower moment rate. This may explain our observation that this SSE did not produce a detectable triggering effect for earthquakes of M4.3 or larger.Fig. 4 Estimated background seismicity for the Mexican dataset. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu ^i$$\end{document} with i=0,1,2,3,4 are estimated using Algorithm 1. The white contours in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu ^1$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu ^2$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu ^3$$\end{document} represent the slip distributions with 4-cm slip intervals for the 2001-2002, 2006, and 2009-2010 SSEs, respectively, digitized from^73^. In \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu ^4$$\end{document} , the slip contour of the 2014 SSE is shown.Fig. 5. Significant \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T^i$$\end{document} values for the Mexican dataset. The \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T^i$$\end{document} values with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$i=1,2,3,4$$\end{document} that are significantly different from zero are shown.
To show the validity of the estimation presented in this work, the genealogy of earthquakes defined by the estimators obtained using Algorithm 1, is presented in Fig. 6. In this figure, all the earthquakes in our dataset are shown with arrows pointing to their aftershocks, in the figure only arrows where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P^O_{ij}\ge 0.5$$\end{document} are presented.
An anomalous arrow was observed originating from the earthquake with coordinates \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-102.540\hbox {E}^{\circ }$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$18.026\hbox {N}^{\circ }$$\end{document} . This may be due to the fact that the time between both earthquakes was only 6.86 hours, and the fact the second earthquake is in a square where only 3 events were recorded throughout the entire study period. This arrow may be an unreliable result.Fig. 6. Offspring estimation for the Mexican dataset. Map of earthquakes with arrows pointing from events to their offspring according to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P^O$$\end{document} .
To clarify the results in Fig. 6, a zoomed-in view of indices 183 to 192 is presented in the Fig. 7a. This window was selected because it contains the largest earthquake (event 183) in our data set, with a magnitude of 7.6. As it can be seen in the Fig. 7a, earthquake 190, with an epicenter at -104.37E \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^\circ$$\end{document} ,18.54N \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^\circ$$\end{document} , could be an offspring of the earthquakes 183,184, 188 or 189. Therefore, the 4 arrows pointing to it are lighter in color than the arrow between earthquakes 183 and 184, which have an associated probability close to 1.
Figure in 7b allows not only to see which earthquakes are background events and which are offsprings, but also to recover an estimated entire genealogy—an idea that has been previously explored in the works of Rasmussen^48^ and Fox et al.^38^.Fig. 7. Genealogy of events 183-192. In (a), the matrix \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P^O$$\end{document} for indices 183 to 192 is shown. (b) Presents the same events as a genealogy.
The b-value for our dataset is 1.254 and the estimated vector of aftershock parameters \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(A,\alpha ,c,p,d)$$\end{document} is (0.118 events, 1.112, 0.021 days, 1.363, 0.0048 deg \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^2$$\end{document} ). The expected total value of earthquakes using the data in this section and Algorithm 1 is
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \int _{0}^T \iint _S \lambda (t,x,y|\mathcal {H}_t)=774.15, \end{aligned}$$\end{document}using our sample of 794 earthquakes. Figure 8a shows a histogram of the values (i.e., the probability of each earthquake being a background event):
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} (1-\rho _j)=\sum _{s=0}^4 \hat{p}_{ii}^s. \end{aligned}$$\end{document}In this figure, we observe a bimodal density concentration near 0 and 1, which is similar to the one presented by Zhuang et al.^15^. The probability of being an aftershock is considerably lower. This could be due to the use of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$M_0=4.3$$\end{document} , because, as shown in Fig. 8b, for the Boso Peninsula in Japan ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$M_0=3$$\end{document} ), there are more aftershocks in proportion than in the Mexican case.Fig. 8. Background probability histogram. Histograms of probabilities of being a background event for all earthquakes in our dataset are presented. (a) Mexican dataset (b) Japanese dataset.
Japan
According to^30^, the magnitudes of the observed SSEs are 6.67 and 6.65 for the 2002 and 2007 events, respectively. Although they had shorter durations than the SSEs in Guerrero, Mexico, the SSEs in the Boso Peninsula, Japan, have exhibited notable increases of the background seismicity rate, as it can be observed in Fig. 9. In both cases, the square with the highest value has a vertex \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$140.5\hbox {E}^{\circ }$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$35.75\hbox {N}^{\circ }$$\end{document} , where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu ^1$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu ^2$$\end{document} reach values of 2.65 and 3.02, respectively, while \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu ^0$$\end{document} for the same square was 0.03.
For both SSEs, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T^i$$\end{document} for the squares associated with the maximum value of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu ^1$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu ^2$$\end{document} was statistically significant, as shown in Fig. 10, and the the LRT rejects the null hypothesis (p-value = 0.01) in favor of our model. Also, as expected, the increase in the background seismicity rate was observed in the margins of the SSE patches. This is consistent with the Mexican resultsFig. 9Estimated background seismicity for the Japanese dataset. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu ^i$$\end{document} with i=0,1,2 estimated using Algorithm 1 are shown. The white contours represent the daily slip rates presented in Fig. 3 for the 2002 and 2007 SSE.Fig. 10. Significant \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T^i$$\end{document} values for the Japanese dataset. The \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T^i$$\end{document} with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$i=1,2$$\end{document} that are significantly different from zero are presented.
The b-value for this data is 0.709, and the estimated vector of aftershock parameters \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(A,\alpha ,c,p,d)$$\end{document} is (0.365 events, 1.03, 0.002 days, 1.029, 4.101e-05 deg \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^2$$\end{document} ). And we obtain
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \int _{0}^T \iint _S \lambda (t,x,y|\mathcal {H}_t)=554.10, \end{aligned}$$\end{document}for our sample of 575 earthquakes. As in “Mexico” section, Fig. 11 presents the estimated genealogy.Fig. 11. Offspring estimation for the Japanese dataset. Map of earthquakes with arrows pointing from events to their offspring according to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P^0$$\end{document} .
Discussion and conclusions
We have successfully developed an inhomogeneous spatio-temporal ETAS model. Our approach extends the ideas of Fox et al.^38^ by introducing a piecewise constant spatial background seismicity over time, which contrasts with their non-time-dependent approach. Furthermore, this approach can be regarded as the spatio-temporal extension of the model proposed by Okutani & Ide^22^, which considers only the temporal components.
We also extended the work of Li et al.^16^ by employing a marked point process. Furthermore, we did not assume that spatial and time contributions of the background seismicity interact solely through multiplication, and we accounted for differences in the spatial distribution of increased background seismicity for each SSE.
Unlike Marsan et al.^27^ and Reverso et al.^28^, which rely on subjectively chosen local spatio-temporal windows, we used the estimation of the SSE durations from GNSS data and utilized them in our seismicity modeling. Moreover, rather than using an approximate approach that sequentially optimizes changes in the background seismicity rate across numerous small spatio-temporal windows, our method simultaneously performs a maximum likelihood estimation of spatio-temporal variations in the background seismicity rate over the entire analysis region and period, which represents a key technical advance of our model.
Furthermore, for the first time, we quantitatively evaluated the relationship between seismic activity and SSEs in the Mexican subduction zone using a statistical seismicity model that accounts for not only regular earthquakes but also SSEs. Our results from Mexico demonstrate statistically significant increases in the background seismicity accompanying Guerrero SSEs, highlighting the seismological novelty of this study.
Our model also successfully detected statistically significant increases in background seismicity associated with SSEs in the offshore Boso Peninsula region of Japan. Despite the substantial differences in the time scales of SSEs between the Boso Peninsula and Guerrero, Mexico, consistent results were obtained—namely, an increase in background seismicity near the margins of the SSE patches. This demonstrates the adaptability and broad applicability of our model to different datasets.
The assumptions made for our model were flexible enough to be applied in different regions, and we expect to continue improving the model and exploring other plate boundaries where SSEs occur, such as in New Zealand and Peru^12,13^. Applying our model to these regions will enable a deeper understanding of how SSEs trigger fast earthquakes.
In future work, we plan to extend our algorithm to allow dependencies between mainshock and aftershock magnitudes. Currently, we model magnitudes independently. However, it is important to explore possible dependencies between magnitudes, as Nandan et al.^59^ found evidence of such dependencies in California, USA.
We also intend to improve our algorithm by replacing the piecewise estimators of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu ^i$$\end{document} with smooth functions in the space domain. In the context of the ETAS model, interesting ideas have recently been explored for this purpose, such as the use of Dirichlet processes^32^ (DP), or Gaussian processes^35^ (GP). An additional advantage of using DP or GP is that they allow inference to be performed within a Bayesian framework.
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1Ozawa, S. et al. Characteristic silent earthquakes in the eastern part of the Boso Peninsula, Central Japan. Geophys. Res. Lett.30 (2003).
- 2Kostoglodov, V. et al. A large silent earthquake in the Guerrero seismic gap, Mexico. Geophys. Res. Lett.30 (2003).
- 3Schwartz, S. Y. & Rokosky, J. M. Slow slip events and seismic tremor at circum-pacific subduction zones. Rev. Geophys.45 (2007).
- 4Hainzl, S. & Ogata, Y. Detecting fluid signals in seismicity data through statistical earthquake modeling. J. Geophys. Res. Solid Earth 110 (2005).
- 5Llenos, A. L. & Mc Guire, J. J. Detecting aseismic strain transients from seismicity data. J. Geophys. Res. Solid Earth 116 (2011).
- 6Antoniak, C. E. Mixtures of Dirichlet processes with applications to Bayesian nonparametric problems. Ann. Stat. 1152–1174 (1974).
- 7Daley, D. J. & Vere-Jones, D. An Introduction to the Theory of Point Processes: Volume I: Elementary Theory and Methods (Springer, 2003).
- 8Helmstetter, A., Sornette, D. & Grasso, J.-R. Mainshocks are aftershocks of conditional foreshocks: How do foreshock statistical properties emerge from aftershock laws. J. Geophys. Res. Solid Earth 108 (2003).