(this html report is better visualized on an horizontal screen, notepad or monitor)
1 Abstract
Chilean seismic activity is among the strongest ones in the world. As already shown in previous papers, seismic activity can be usefully described by a space-time branching process, like the ETAS (Epidemic Type Aftershock Sequences) model, which is a semiparametric model with a large time scale component for the background seismicity and a small time scale component for the induced seismicity. The large-scale component intensity function is usually estimated by nonparametric techniques, specifically in our paper we used the Forward Likelihood Predictive approach (FLP); the induced seismicity is modelled with a parametric space-time function. In classical ETAS models the expected number of induced events depends only on the magnitude of the main event. From a statistical point of view, forecast of induced seismicity can be performed in the days following a big event; of course the estimation of this component is very important to forecast the evolution, in space and time domain, of a seismic sequence. Together with magnitude, to explain the expected number of induced events we also used other covariates.
According to this formulation, the expected number of events induced by event \(E_i\) is a function of a linear predictor \(\boldsymbol{ \eta}_i=\mathbf{x}_i^{\mathsf{T}} \boldsymbol{ \beta}\), where \(\mathbf{x}_i\) is the vector of covariates observed for the \(i\)-th event (the first covariate is usually the magnitude \(m_i\)), and \(\boldsymbol{ \beta}\) is a vector of parameters to be estimated together with the other parametric and nonparametric components of the ETAS model.
We obtained some interesting result using some covariates related to the depth of events and to some GPS measurement, corresponding to earth movement observed until the time of main events. We find that some of these models can improve the description and the forecasting of the induced seismicity in Chile, after a subdivision of the country in two different spatial regions.
We used open-source software (R package etasFLP, Chiodi and Adelfio (2017)) to perform the semiparametric estimation of the ETAS model with covariates. How to cite: Chiodi M. (2021)
2 Whole Chile 2000-2020: Circles extension proportional to magnitude; red: recent, blu: older
2.1 Distribution of earthquakes by month and magnitude
2.2 Candidate covariates present in the catalog
The catalog contains the ordinary seismic variables and some more information on other variables related to each observed event, as detailed later, in a further section (clickgo).
2.3 Two areas
The choice of the subdivision in these two areas is due to the fact that the Nazca plate is being subducted beneath South America at a rate of 90 mm/yr in the north and approximately 40mm/yr in the center-sur of Chile. Additionally, the considered areas includes the major earthquakes of the last 30 years. The Northern area include the M8.2 Iquique earthquake (Abril, 2, 2014) and the second area include the 2010 Chile earthquake of 8.8 magnitude occurred off the coast of central Chile on Saturday, 27 February 2010 which provoked a rupture of more 1000 km length, and the M8.4 Coquimbo earthquake occurred on September 16, 2015.
2.4 Simple completeness analysis (log transform of back cumulated frequencies)
On the basis fof these information a value of magnitude \(m_0 = 2.7\) has been taken as a magnitude threshold for both areas, and only events from 2006 have been taken, because an empirical completeness analysis of catalogs before 2006 shows that events before 2006 probabily are not identified with the same completeness magnitude.
(Nicolis, Chiodi, and Adelfio 2015, 2017)
2.5 Distribution of depth
We decided to consider only events at a depth less 70 km in both areas. Earthquakes deeper than 70 km have been excluded because they occur within the continental plate in the crust. They are mainly due to the deformations generated by the convergence between the Nazca plate and the continental plate. These deformations caused the rise of the Andes Mountains and they can reach magnitudes greater than 7 (Barrientos et al. (2004))
3 ETAS model for the analysis fo seismicity in Chile
Based also on previous papers, we tentatively fitted two ETAS models to the two areas obtaining some interesting results.
This triggering function is factorized into separate effects of marks, time, and relative location:
\[ {\lambda_{\tilde{\boldsymbol{ \theta}}}}(t,\textbf{s}|\mathcal{H}_t ) \ = \ \mu f(\textbf{s}) \ + \ \sum_{t_j<t}\frac{\kappa_0 \exp{(\eta_j)} }{(t-t_j +c)^p} \left[(\textbf{s}-\textbf{s}_j)^2+d\right]^{-q} \]
where \((t_j , \textbf{s}_j )\) is the time and location of individual occurrence \(j\), \(\eta_j = \boldsymbol{ \beta}^{\mathsf{T}} \mathbf{Z}_j\) is a linear predictor, with \(\mathbf{Z}_j\) the external known covariate vector, including the magnitude (usually coinciding with the first covariate), acting in a multiplicative fashion on the base risk and \(\boldsymbol{ \theta}=(\mu, \kappa_0, c, p, d, q, \boldsymbol{ \beta})^{\mathsf{T}}\), with \(\boldsymbol{ \beta}\) a \(k\)-component vector, to be estimated together with the other 6 components of \(\boldsymbol{ \theta}\) (Adelfio and Chiodi (2020)); \(f(\textbf{s})\) is estimated with a non-parametric technique based on the FLP method (Forward Likelihood Predictive) proposed by Adelfio and Chiodi (2015).
More in details, in the usual ETAS model (Ogata 1988, 1998), \(k=1\), \(\mathbf{Z}_{j1}=m_j-m_0,\) and \(\beta_1=\alpha\). In this model formulation, for an easier correspondence with the ETAS parametrization, in the \(\boldsymbol{ \beta}\) vector an intercept term is not included, because of the presence of the parameter \(\kappa_0\) in the model; in the spatial component \((\textbf{s}-\textbf{s}_j)\), we used simple euclidean distances.
3.1 Possible models and covariates
We tried different models based on covariates to be added to magnitude to explain the induced seismicity. The variables present in the catalog are:
- time: Time of the event, in days with decimal fractional part.
- long: longitude of the event.
- lat: latitude of the event.
- magn1: magnitude of the event.
- z: depth.
- HDist: Haversine distance from the seismic event to the Nazca plate in meters.
- Station: GPS station closest to the seismic event.
- ELon and ELat: longitude and latitude of the GPS station in the day of the seismic event.
- Angle and speed: angle and speed of the GPS station in a 5-day interval until the day of the seismic event.
- HVQuk: Haversine distance from the GPS station to the seismic event in meters.
- Mt1D: displacement of the station in meters, in an interval of one day before the seismic event
- Mt5D: displacement of the station in meters, in an interval of five days before the seismic event
| lat | long | z | magn1 | time | HDist | Estacion | ELon | ELat | Angulo | Velocidad | HVQuk | Mt1D | Mt5D | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | -38.05 | -75.94 | 15.0 | 4.2 | 43952.91 | 89334.6 | UDEC | -72.34 | -37.47 | 0.34 | 0 | 322980.65 | 0 | 0.00 |
| 2 | -38.11 | -75.79 | 10.0 | 3.7 | 43809.51 | 87946.5 | ANG8 | -72.69 | -37.80 | -1.07 | 0 | 274209.58 | 0 | 0.01 |
| 3 | -37.69 | -75.65 | 15.0 | 3.6 | 43778.63 | 86240.0 | ANG8 | -72.69 | -37.80 | -0.03 | 0 | 261134.05 | 0 | 0.02 |
| 16737 | -30.79 | -71.84 | 28.2 | 3.0 | 42596.75 | 71036.2 | EMAT | -71.66 | -31.15 | -0.20 | 0 | 42677.98 | 0 | 0.01 |
| 16738 | -29.96 | -71.84 | 38.8 | 2.7 | 42621.78 | 8928.9 | BTO1 | -71.49 | -30.26 | -1.14 | 0 | 47518.86 | 0 | 0.01 |
| 16739 | -34.71 | -71.84 | 11.6 | 2.5 | 42669.43 | 209647.4 | HLN2 | -71.93 | -35.01 | -0.98 | 0 | 34861.38 | 0 | 0.01 |
| 33476 | -33.70 | -70.00 | 3.7 | 2.9 | 41157.22 | 237126.1 | BN06 | -70.74 | -34.17 | -0.74 | 0 | 86222.26 | 0 | 0.00 |
| 33477 | -33.16 | -70.00 | 2.6 | 2.3 | 39244.62 | 59788.3 | CMBA | -71.00 | -31.19 | -1.16 | 0 | 239300.50 | 0 | 0.01 |
| 33478 | -31.82 | -70.00 | 40.4 | 2.7 | 40013.74 | 281486.3 | PEDR | -70.69 | -30.84 | 1.45 | 0 | 127690.34 | 0 | 0.01 |
3.1.1 Software: the open source package etasFLP
We used the open source R package etasFLP with the new version vith covariates (Chiodi and Adelfio 2021) on CRAN repository (R Core Team 2020). It performs a semiparametric estimate of the components of the model, with many statistical and computational options (Chiodi and Adelfio 2017). An updated version of the software will be published on CRAN in few weeks, but could be obtained asking the first Author.
3.2 Choice of covariates
The strategy followed to choose covariates in both areas is a forward selection: after a fitting with the first covariate magnitude, other covariates are added one by one, and the gain in AIC is computed. The best result, in terms of AIC, was obtained adding the variable depth, \(z\). So the depth \(z\) is introduced in the model, and a third variable, chosen from the candidate covariates, is added one by one, and the gain computed. The best result is obtained with the variable Mt5D, which is the third covariate introduced in the model for aftershocks. As a final step, a fourth covariate,Mt1D, gave some further improvement in AIC value. No other covariates led to appreciable improvements of AIC.
The ranking of variables and the order of inclusion of covariates, related to induced seismicity, was the same in both areas, even with different values of parameters estimates, so that for both areas the vector of covariate was: \(\left\{magnitude, z (the depth), Mt5D, Mt1D\right\}.\)
3.3 Estimate of \(p\).
The only difference between models used in the two areas, concerned the parameter \(p\):
Indeed a first procedure for both areas has been made with a fixed value of \(p\), very close to 1. However the estimation of \(p\), even if useful in term of AIC gain in northern area, led to values very close to 1, even if significantly different from 1 for north area.
The value of \(p\) is significantly different for 1 in the north area, and below a plot of the profile log-likelihood is reported for the northern area. The Akaike Information Criterion (AIC) of the model with p included is sensibly better than the AIC value for the model with p fixed.
Instead, in the south area the value of \(p\) is not significantly different from 1, and the best AIC value is obtained in the model with \(p\) fixed to a value very close to 1.
3.4 Report for Northern area
Table of AIC values for four models fitted to the data of the north catalog, with \(n=8049\) valid observations:
| covariates | |
|---|---|
| magnitude | 168471.2 |
| magnitude+depth | 168287.5 |
| magnitude+depth+Mt5D | 168048.1 |
| magnitude+depth+Mt1D+Mt5D | 168029.8 |
| covariates | |
|---|---|
| magnitude | 168524.3 |
| magnitude+depth | 168375.7 |
| magnitude+depth+Mt5D | 168129.9 |
| magnitude+depth+Mt1D+Mt5D | 168098.6 |
The gain in AIC values using a value of \(p\) estimated is quite evident.
3.5 Report for Southern area
In South area, with \(n=27780\) valid observations, a fixed value of \(p\) was generally preferred on the basis of the AIC.
| covariates | |
|---|---|
| magnitude | 554814.9 |
| magnitude+depth | 554385.3 |
| magnitude+depth+Mt5D | 554319.3 |
| magnitude+depth+Mt1D+Mt5D | 554294.3 |
| covariates | |
|---|---|
| magnitude | 554762.6 |
| magnitude+depth | 554397.4 |
| magnitude+depth+Mt5D | 554303.1 |
| magnitude+depth+Mt1D+Mt5D | 554283.3 |
In south area the situation is generally in favour of a fixed value of \(p\), especially for the model with 3 and 4 covariates.
Again, the big gain in AIC is obtained at the second step, with the insertion of \(z\) (the depth) in the model for induced seismicity. Of course the main contribution is always given by the magnitude.
3.6 Comparison of results in the two Chileean zone.
3.6.1 North area (model with 4 covariates):
Formula for covariates of the triggered components: ~magn1+z+Mt1D+Mt5D-1. \(p\) has been estimated from data together with other parameters.
| Estimates | std.err. | |
|---|---|---|
| \(\mu\) | 0.6432 | 0.0147 |
| \(\kappa_0\) | 0.0066 | 0.0015 |
| \(c\) | 0.0183 | 0.0021 |
| \(p\) | 1.0417 | 0.0083 |
| \(d\) | 13.2706 | 1.1801 |
| \(q\) | 1.6717 | 0.0330 |
| magn1 | 0.7303 | 0.0267 |
| \(z\) | -0.0118 | 0.0017 |
| Mt1D | 1.5592 | 0.2551 |
| Mt5D | 1.3036 | 0.1341 |
3.6.2 South area (model with 4 covariates, p fixed):
Formula for covariates of the triggered components: ~magn1+z+Mt1D+Mt5D-1. \(p\) has been kept fixed since there was no gain in AIC.
| Estimates | std.err. | |
|---|---|---|
| \(\mu\) | 1.2686 | 0.0232 |
| \(\kappa_0\) | 0.0319 | 0.0045 |
| \(c\) | 0.0142 | 0.0006 |
| \(p\) | 1.0001 | 0.0000 |
| \(d\) | 37.2547 | 1.7399 |
| \(q\) | 1.7826 | 0.0199 |
| magn1 | 0.6843 | 0.0145 |
| \(z\) | -0.0104 | 0.0007 |
| Mt1D | -0.1050 | 0.0273 |
| Mt5D | 0.2030 | 0.0129 |
3.7 Plots for North Area (4-covariates model)
Here we report some plot for the ETAS model fitted to the north area with 4 covariates, and a value of \(p\) estimated:
3.7.1 South Area
## Comparison between model residuals: North area.
We compared different models also by mnean of residual analysis (Schoenberg (2003)). briefly here we report the tranformed time residuals: for the right model they should follow the bisector line. In the plot the comparison between tranformed time residuals for two models is reported: a black line for the 4-covariates model and a green line for the 1-covariate model. The red line is the (theoretical) bisector line.
Black line seems to behave better, expecially after a first departure from theoretical model.
3.8 Comparison between model residuals: South area.
The same comparison for the south area gave not similar results, in the sense that residual for the 4-covariates and for the 1-covariate models seem to behave similarly.
References
Adelfio, G., and M Chiodi. 2015. “Alternated Estimation in Semi-Parametric Space-Time Branching-Type Point Processes with Application to Seismic Catalogs.” Stochastic Environmental Research and Risk Assessment 29: 443–50.
Adelfio, Giada, and Marcello Chiodi. 2020. “Including Covariates in a Space-Time Point Process with Application to Seismicity.” Statistical Methods & Applications, July. https://doi.org/10.1007/s10260-020-00543-5.
Barrientos, S., E. Vera, P. Alvarado, and T. Monfret. 2004. “Crustal Seismicity in Central Chile.” Journal of South American Earth Sciences 16 (8): 759–68. https://doi.org/https://doi.org/10.1016/j.jsames.2003.12.001.
Chiodi, M., and G. Adelfio. 2017. “Mixed Non-Parametric and Parametric Estimation Techniques in R Package etasFLP for Earthquakes’ Description.” Journal of Statistical Software 76 (3): 1–28. https://doi.org/10.18637/jss.v076.i03.
Chiodi, Marcello, and Giada Adelfio. 2021. EtasFLP: Mixed Flp and Ml Estimation of Etas Space-Time Point Processes for Earthquake Description.
Chiodi M., Adelfio G., Nicolis O. 2021. “ETAS Space Time Modelling of Chile Induced Seismicity Using Covariates.” In. EGU General Assembly 2021, online. https://doi.org/10.5194/egusphere-egu21-3523.
Nicolis, Orietta, Marcello Chiodi, and Giada Adelfio. 2015. “Windowed Etas Models with Application to the Chilean Seismic Catalogs.” Spatial Statistics 14: 151–65. https://doi.org/https://doi.org/10.1016/j.spasta.2015.05.006.
———. 2017. “Space-Time Forecasting of Seismic Events in Chile.” In Earthquakes, edited by Taher Zouaghi. Rijeka: IntechOpen. https://doi.org/10.5772/66339.
Ogata, Y. 1988. “Statistical Models for Earthquake Occurrences and Residual Analysis for Point Processes.” Journal of the American Statistical Association 83 (401): 9–27.
———. 1998. “Space-Time Point-Process Models for Earthquake Occurrences.” Annals of the Institute of Statistical Mathematics 50 (2): 379–402.
R Core Team. 2020. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. https://www.R-project.org/.
Schoenberg, F. P. 2003. “Multi-Dimensional Residual Analysis of Point Process Models for Earthquake Occurrences.” Journal American Statistical Association 98 (464): 789–95.