ETAS Space–Time Modeling of Chile Triggered Seismicity Using Covariates: Some Preliminary Results

: Chilean seismic activity is one of the strongest in the world. As already shown in previous papers, seismic activity can be usefully described by a space–time branching process, such as 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 triggered seismicity. The use of covariates can improve the description of triggered seismicity in the ETAS model, so in this paper, we study the Chilean seismicity separately for the North and South area, using some GPS-related data observed together with ordinary catalog data. Our results show evidence that the use of some covariates can improve the ﬁtting of the ETAS model.


Introduction
Chile is considered one of the most active seismic countries in the world due to its particular location. In particular, Chile is located on three great tectonic plates: Nazca, South American, and Antarctic. The seismic events in Chile result from the convergence of the oceanic lithosphere of the Nazca plate beneath the South American continental plate at a rate of about 6.5 cm/yr [1]. The strong coupling between the two plates produces the most intense seismicity and the largest earthquakes in the country [2][3][4]. Seismic activity induces various types of ground movements, as seismic waves (P and S-waves), which can cause significant damages to structures, and consequently losses of lives [5,6]. The velocity of the waves can be different according to the depth [7,8].
Seismic activity can be described by a spatio-temporal branching process, such as the ETAS (Epidemic Type Aftershock Sequences) model [9], which is a semiparametric model with a large time-scale component for the background seismicity and a small time-scale component for the triggered 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) approach [10][11][12]; the triggered seismicity is modeled with a parametric space-time function. In classical ETAS models, the expected number of triggered events depends only on the magnitude of the main event. In previous papers [13,14] we already used the classical ETAS model to describe Chilean seismicity, also using different models for different Chilean areas.
From a statistical point of view, forecast of triggered 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 of a seismic sequence in space and time domain.
In the classical ETAS model [9], the average number of future triggered events just depends on the magnitude of the previous event. Although that represents a very comfortable and useful model, in many situations, it may be too simple to describe the aftershocks dynamics, as often highlighted by diagnostic analysis. Furthermore, the high seismicity of Chile, together with the recent increasing GPS data availability, allowed us to consider a much more complex model for describing and analyzing the Chilean triggered activity.
In this paper, we provide a description of the triggered seismicity in Chile after a subdivision of the country into two spatial areas. In particular, based on the methodology proposed by [15], we suggest the use of a specific branching-type model for earthquake description (the ETAS model) in a regression-oriented version, accounting also for external covariates, related to the depth of events and to some GPS measurements. For further examples on applications incorporating the external information in self-exciting models, see [16][17][18][19].
Interesting and encouraging results, though not completely satisfactory, are provided, accounting for covariates related to the depth of events and to some GPS measurements, corresponding to the Earth movement observed until the time of main events. In addition, some of the proposed models can improve the description and the forecasting of the triggered seismicity in Chile, if compared to the classical ETAS model.
All the reported results are obtained using the open-source software developed in R [20,21]. This paper is organized as follows. In Section 2 a description of data is provided jointly to some preliminary analysis. Section 3 reports some details of the theoretical background, recalling some basic theory of space-time point processes and the ETAS model, with an extension of the FLP approach with the inclusion of covariates. The choice of covariates is described in Section 4. A seismic application to the considered data is discussed in Section 5. Section 6 is devoted to discussion of results. Conclusions are reported in Section 7.

Description of the Datasets
In this paper, two datasets are analyzed: the catalog of the Chilean seismic events and the GPS measurements in the area between latitudes [−37.99, −31.00] and longitudes [−74.99, −71.00].

The Seismic Catalog and Chilean Seismicity
The Chilean catalog has been provided by the Chilean National Seismological Center of the University of Chile (http://sismologia.cl (accessed on 10 February 2021)). It contains the ordinary seismic variables, such as latitude, longitude, magnitude (on the Richter scale), depth, and time of occurrence of each seismic event (year, month, day, hour, minute, and seconds).
From Figure 1, it is observed that the largest number of seismic events occurred in the years 2010, 2014, and 2015, which coincide with the occurrence of the highest magnitude earthquakes, of magnitudes 8.8, 8.2 and 8.4 on the Richter scale, respectively. The relationship between month and magnitude is explicitly represented in Figure 2. From this figure it is also observed that the original seismic data are normally collected with a threshold in the minimum magnitudes, since very low magnitudes are difficult to be detected.
Some summary statistics, related to the historical seismic activities in Chile in the last 21 years, are reported in Table 1. This table shows that 2010 has the maximum number of events (second column) and the maximum magnitude (third column) as well. In the years 2000 and 2011, two seismic events, with the maximum magnitude in two different spatio-temporal locations, occurred.
Since the seismicity has a different behavior along the coast of Chile, due to the different velocity that the Nazca Plate ( Figure 3, left panel) is being subducted beneath South America (see, [22]), we considered two distinct areas (see Figure 3, right panel). We consider only events with depths less than 70 km in both areas (Figures 4 and 5). Earthquakes deeper than 70 km, which occur within the subducted oceanic lithosphere [8], are excluded in our analysis for both areas. The reason is that the largest earthquakes occur at shallow depths down to 50 km and near the limit of the intermediate-depth seismic segment at 60-70 km [1], beneath the forearc section of the convergent margin. Additionally, most seismic events in Chile occur at depths between 30 and 40 km (Figures 4 and 5). These 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 [22].   The threshold magnitude m 0 for the ETAS model, i.e., the lower bound for which earthquakes with higher values of magnitude are surely recorded in the catalog, is obtained by analyzing the log-cumulate distribution of the magnitudes recorded in the two areas (see, Figure 6).
Based on this information, a value of magnitude m 0 = 2.7 is taken as the threshold magnitude for both areas. Starting from an empirical completeness analysis of the catalogs before 2006, events before that date are not identified with the same completeness magnitude. Therefore, only events from 2006 are considered in the proposed study [13,14].
The estimated b-values are 0.586 for the North zone and 0.542 for the Southern one, with standard errors 0.0104 and 0.0050, respectively.

The GPS Measurements
For calculating the daily accumulated displacements associated with seismic events, we use the GPS measurements (longitude, latitude and altitude above sea level, a.s.l.), downloaded by the website http://geodesy.unr.edu/NGLStationPages/GlobalStationList (accessed on 22 February 2021). Figure 7 represents the daily coordinates of a GPS station (which name is EMAT) located in the Coquimbo region in the period 2005-2020. The calculation of the accumulated daily horizontal displacements is performed using the GPS measurements of both one day and five days prior to the seismic event, collected by the GPS station closest to the same event. For this calculation, the Euclidean distance, given by (x 1 − x 2 ) 2 + (y 1 − y 2 ) 2 , is used. The GPS altitude component, normally used for evaluating the vertical displacement, has not been considered in this application. As documented in various studies, this variable does not have sufficient accuracy, and it is not largely used to describe seismic events [23,24]. Therefore, in the proposed analysis we have not included this variable as potential covariate. To acquire the accumulated displacements (called Mt1D and Mt5D, for one and five days, respectively) the data were pre-processed by interpolating the missing values using the R package imputTS, [25], and denoising them through the R package wmtsa [26].

Branching Point Processes and ETAS Model
Branching processes are used to model reproduction phenomena. These models have been recently considered for the description of different application fields: biology [27], demography [28], epidemiology [29,30].
Any analytic spatio-temporal point process is uniquely characterized by its associated conditional intensity function (CIF) [31] that represents the instantaneous rate or hazard for events at time t and location s given all the observations up to time t, conditioning on the random past history of the process H t .
The conditional intensity function of the branching model is defined as the sum of a term describing the large time-scale variation (spontaneous activity or background) and one relative to the small time-scale variation due to the interaction with the events in the past (triggered activity or offspring): with θ = (φ, µ) , the vector of parameters of the triggered intensity (φ) together with the parameter of the background general intensity (µ), f (s) the space density, and τ φ (t, s) the triggered intensity, given by: where ν φ (t − t j , s − s j ) is the space-time intensity at (t, s) triggered by a previous event at time t j . The self-exciting component of the model basically provides a description of the intensity at a space-time location (t, s) caused by each previous event.
As introduced above, a branching process for earthquake description, widely used in seismological context is the Epidemic Type Aftershocks-Sequences (ETAS) model [9].
The ETAS conditional intensity function can be written, starting from model (1), as follows: where f (s) is the space density of the background/long term component, stationary in time. The aftershock (triggered) component is the product of the density of aftershocks in time, i.e., the Omori law representing the occurrence rate of aftershocks at time t, following the earthquake of time t j and magnitude m j , and the density of aftershocks in space. In particular, m j is the magnitude of the j-th event and m 0 the threshold magnitude, k is a normalizing constant, c and p characteristic parameters of the seismic activity of the given region; p is useful for characterizing the pattern of seismicity, indicating the decay rate of aftershocks in time; d and q are two parameters related to the spatial influence of the mainshock. For the spatial triggered distribution, conditioned to magnitude of the generating event, the distribution is: relating the occurrence rate of aftershocks to the mainshock magnitude m j , through the parameters α, γ that measure the influence on the relative weight of each sequence. The FLP approach, which is a nonparametric estimation procedure based on the subsequent increments of log-likelihood obtained by adding an observation one at a time, to account for the information of the observations until t k on the next one, has been developed for estimating the background component [11,32].
The package etasFLP [11,33,34] provides tools to implement this mixed approach for a wide class of ETAS models for the description of seismic events, developed in an R environment [20].

ETAS Model with Covariates
The ETAS model is widely used in the framework of statistical modeling of earthquake occurrence; however, in some applications, it may not adequately describe the triggered seismicity. In the classical ETAS model formulation, as in Equation (2), the average number of events triggered by an event at time t j depends only on its magnitude m j . Ref. [15] proposed an ETAS model with covariates in the framework of the survival analysis. This model relates the average number of events triggered by an event at time t j with the value of some covariates, including magnitude, corresponding to the j−th event. From a statistical point of view, the inclusion of new explanatory variables can significantly improve the fitting of the model to real observed data.
These variables may also vary continuously in space, and their effect is incorporated in the triggering part of the model. The background component of the model is estimated using the Forward Likelihood for Prediction (FLP) approach [10][11][12].
As proposed by [35] in a context of infection occurrences, in [15] we incorporate the space-time phenomenological laws of the triggering part of the ETAS model with the effects of covariates. In particular, we model the covariates of the ETAS model as in a GLM framework, such that η j is a classical linear predictor given by η j = β Z j , where Z j is the vector of covariates observed for the j-th event and β is a vector of unknown parameters. Therefore, the triggering function is factorized into separate effects of marks, time, and relative location: where (t j , s j ) is the time and location of individual occurrence j, η j = β Z j is a linear predictor, with 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 θ = (µ, κ 0 , c, p, d, q, β) , with β a k component vector, to be estimated.
More in details, in the usual ETAS model [9,36], k = 1, Z j1 = m j − m 0 , and β 1 = α − γ. In this model formulation, for an easier correspondence with the ETAS parametrization, in the β vector a constant term is not included, since the presence of the parameter κ 0 in the model.
In the seismic context, that would provide a more general formalism for the earthquake occurrence in space and time. Indeed, the main idea is that the effect on the future activity depends not only on the closeness of the previous events but also on other characteristics of the main event, such as magnitude, as usual, and the distance from the faults, or other geological sources.
An extended version of the package etasFLP for generalized offspring component is available at the CRAN [21].

Choosing the Covariates in the ETAS Model Applied to Chile Seismicity
First analysis conducted on Chile catalogs, showed that probably the use of covariates could have improved the performance of the model. Possible choices have also been encouraged by the availability of GPS measurements, which were supposed to be useful in the explanation of triggered seismicity. To explain the intensity of the triggered seismicity, we estimated several models with different sets of covariates, in addition to the magnitude.
In particular, the suitable covariates in both the areas are chosen based on a forward selection strategy: starting from the model with the first covariate magnitude, further covariates are added one by one, and the gain in terms of the Akaike Information Criterion (AIC) is computed. Since the best gain in the AIC values is obtained by adding the variable depth z, it is the second variable introduced in the model. Then, a third variable, chosen from the remaining candidate covariates, is added, 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, gives some further improvement in the AIC value. No further covariates lead to appreciable improvements in terms of AIC.
The ranking and the order of inclusion of the covariates, in the triggered seismicity component, is the same in both the areas, even with different estimated parameters. Then, the selected covariates are: magnitude, z, Mt5D, Mt1D.

Results
The two models used in the analyzed areas differ just for the estimation of the parameter p. Indeed, a first procedure for both the areas is carried out with a fixed value of p, very close to 1. However, for the Northern area, the estimation of p leads to a value very close to (but significantly different from) 1.
In Figure 8 the profile log-likelihood for the parameter p is reported for the Northern area. The AIC of the model estimating also p is sensibly lower than the AIC value of the model with p taken as fixed. Conversely, in the Southern 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.

Comparison of the AIC Values for the Two Areas
The AIC values for the four models fitted to the data of the Northern catalog, with n = 8049 valid observations, are reported in Table 2. The first column reports the AIC values obtained when p is estimated together with other parameters, while the AIC values of the second columns are obtained from models in which p is fixed to 1.0001. The gain in the AIC values is evident in models with p estimated from the data.
Similar AIC values for Southern area (n = 27780) are reported in Table 3. Here, evidence is generally in favor of a fixed value of p, especially for the model with three and four covariates, as already seen from the profile likelihood plot (Figure 8).
The model with four covariates seems to perform better than the simpler models both for the North and South areas.
No further covariates are introduced because no improvement in terms of AIC is observed. Again, the big gain in terms of AIC is obtained at the second step, with the introduction of z (the depth) in the model for triggered seismicity. However, the main contribution is always given by the magnitude.

Comparison of Results in the Two Chilean Zone
For the triggered components for the Northern area, according to the AIC values reported in the previous tables, we use as covariates, besides the magnitude, the depth, Mt1D and Mt5D (see estimation results in Table 4). The parameter p is estimated from data. For the triggered components for the Southern area, according to the AIC values reported in the previous tables, we consider, besides the magnitude, the same set of covariates used for the Northern area, i.e., the depth, Mt1D and Mt5D (see estimation results in Table 5). p is kept fixed, since no gain in terms of AIC is observed.

Discussion of Results
The obtained results suggest that from a simple AIC point of view, models with four covariates must be preferred for both areas, as already underlined in Section 5.1. Moreover, the estimated parameters differently characterize the two areas, as shown in Tables 4 and 5 (note for instance the different sign of the variable Mt1D), reflecting the different behavior of the seismicity along the Chilean coast.
However, we report a diagnostic residual plot to see how the four-covariates models perform well. In general, in diagnostic analysis, if the estimated model is a good approximation of the real model that describes data, then the transformed observations (moving each point on the integrated intensity function) should behave like a homogeneous Poisson process [37]. Deviations from it are caused by some characteristics of data not taken into account by the fitted model. Therefore, in Figure 9 we report the transformed time residual plots for the Northern area, on the left, and for the Southern area, on the right. In both plots, we draw in red the theoretical line (if the model is estimated correctly), in black the line for the transformed residuals obtained fitting an ETAS model with four covariates, and finally in green the line of the transformed residuals fitting the classical ETAS model. From the left panel of Figure 9 it is clear that the four-covariates model can describe the triggering behaviors in the Northern area. Nevertheless, the black and green lines on the right panel of Figure 9 related to the South area are very close to each other, showing that probably the fitted four-covariates ETAS model does not fully describe the complex triggering mechanism of the Southern area.

Conclusions
In this paper, we present an application of the ETAS model with covariates to describe the Chilean seismicity, also using GPS information of the observed data. In particular, the triggered seismicity in Chile is analyzed distinguishing between the North and the South of the country, since the different behavior of the seismicity along the coast of Chile. According to the analyzed areas, different results are obtained, accounting for the main relevant covariates related to the depth of events and some GPS measurements, corresponding to the Earth movement observed until the time of main events. It represents an innovative method, improving the assessment of seismic events in space and time, considering a contagion model within a regression-like framework accounting for external covariates.
Moreover, the introduction of covariates to the classical ETAS model seems to improve the description of the Chilean seismicity, explaining better the overall variability of the studied phenomenon and leading to a decrease in the unpredictable variability. Indeed, in the seismic context, the proposed approach can provide more general formalism for the earthquake occurrence in space and time, such that the effect on the future activity does not depend only on the closeness of the previous events but also on specific characteristics of the main event, such as magnitude, as usual, and further information, such as geological features or some GPS-related data.
The reported results confirm the need for a more flexible model for describing the complexity of the spatio-temporal aftershock activity. Moreover, introducing external information is an innovative and promising perspective of study that could be also relevant for studying different phenomena with epidemic features.
The proposed approach could still be improved and would require more investigation, both in terms of inference and diagnostic results (e.g., accounting for directional effect for a more realistic analysis). Additionally, the proposed model could be used for predicting future events in space and time such as in [14,38].