Computing the time-dependent activity rate using non-declustered and declustered catalogues – a first step towards time-dependent seismic hazard calculations for operational earthquake forecasting

Probabilistic seismic hazard analysis (PSHA) typically requires tectonic b values and seismic activity rates using declustered catalogues to compute the annual probability of exceedance of a given ground motion (for example, the peak ground acceleration or PGA). In this work, we propose a methodology that includes the spatially gridded time-dependent b value and activity rate computation using seismic clusters in PSHA calculations. To account for the spatial variability and the relationship of the earthquakes with the seismic sources, we incorporate the distance from the centre of the grid cell to the closest fault and the epicentre's uncertainty into the smoothing kernel as the average distance and the variance, respectively. To illustrate this methodology, we selected two scenarios as representatives of the high-seismicity region and low-to-moderate-seismicity region. The first one is located in Central Italy, where the L'Aquila earthquake happened, while the other is in south-eastern Spain, where several earthquakes with a moment magnitude (M_w) greater than 4.0 have taken place over the last 30 years, including two earthquakes with magnitude greater than or equal to 5.0. We compared three different seismic activity models based on the parameters considered in the calculations (distance from spatial cells to faults and epicentral distance uncertainty), and we defined and calculated the changes in the annual probability of exceedance for a given background PGA value. The results reveal noticeable changes in the annual probability of exceedance in the proximity of the occurrence of significant events. In the case of Italy, the annual probability of exceedance increases significantly, but in the case of Spain not all the earthquakes have an associated increase in the exceedance probability. However, we have observed how, for moderate- to low-seismicity regions, the use of a non-declustered catalogue can be appropriate when computing time-dependent PSHA, as in the case of Spain.

1 Introduction

Probabilistic seismic hazard analysis (PSHA) has been the basis for seismic engineering design since Cornell (1968) proposed it in order to account for all the possible earthquake scenarios and ground motion levels that can occur in the different seismic sources affecting the site of interest. One of the key points of PSHA is how the uncertainties are incorporated into the ground motion computation, so the results are much more appropriate for use in engineering decision-making for risk reduction. However, the procedure increases in complexity (Budnitz et al., 1997).

PSHA results depend on combining the pertinent input models (those which, according to the scientific and engineering communities, represent the relevant phenomena in an appropriate way). Therefore, the choice of these models will evolve as our knowledge of the seismic activity and occurrence increases.

PSHA determines the probability of exceeding the ground motion level over a specified time period based on the occurrence rate of earthquakes and ground motion prediction equations (GMPEs). The occurrence rate of earthquakes is generally described by the truncated exponential model (Cosentino et al., 1977) and the characteristic earthquake model (Schwartz and Coppersmith, 1984). Additionally, this earthquake occurrence rate or activity rate is assumed constant during the computation process. Therefore, it provides results which can be used for the seismic design. Once the knowledge of the seismic activity and occurrence improves due to the recording of new rare events or new tectonic information and models, the PSHA can be calculated again, and the seismic building codes will be updated if needed.

On the other hand, many authors have begun to focus the PSHA computations from a temporal or real-time perspective, so the term “time-dependent probabilistic seismic hazard” (TDPSHA) is now widely used. These models are based on how the probabilities of large events increase as stress builds up on a fault plane until it reaches the breaking point of the rock (Kanamori and Brodsky, 2004) and also how the probabilities of large aftershocks are a decreasing function after the main large event (Ogata, 1988; Reasenberg and Jones, 1989). However, measuring changes in the stress caused by the main shock is possible only indirectly and with somewhat low precision.

In general, small earthquakes are more frequent than large earthquakes. This is quantitatively stated in the Gutenberg–Richter law (Gutenberg and Ritcher, 1956) (G–R from now on) that can be seen in Eq. (1):

$$\begin{array}{}\text{(1)}& {\log }_{\mathrm{10}}N(M\ge m)=a-b\cdot m,\end{array}$$

where N is the cumulative number of earthquakes with magnitude M above m; the b value is the average size distribution of earthquakes (which expresses the ratio between high-magnitude and low-magnitude earthquakes); and a is the productivity, or more precisely, 10^a is known as the seismic activity rate. As the PSHA results are given as an annual probability of exceedance for a given intensity of the ground motion, the most common way to work with the activity rate is using the annual activity rate, which is obtained by dividing 10^a by the duration in years of the seismic catalogue. So, if we can identify seismic sources a priori, then the seismic data inside each seismic source are used to compute a source-specific (a and b) magnitude frequency distribution. However, it is often challenging to identify the corresponding boundaries and to have enough data allowing a significant statistical fitting.

Therefore, Frankel (1995) instead of specifying spatial borders for each seismic source adopted a boundary-less source model when computing the PSHA for the central and eastern United States. Under this approach, the historical seismicity is spatially smoothed, and activity rates are computed at a grid of locations through the analysis domain. First, they divided the region into a grid, and then they counted for each cell of the spatial grid the number of earthquakes greater than a reference magnitude (M_ref) depending on the occurrence year of the event (1700 for magnitudes greater than 5.0 M_w and 1924 for magnitudes greater than 3.0 M_w). Next, the author obtained a maximum likelihood estimate for 10^a (Weichert, 1980) that they would then smooth using a Gaussian kernel with a correlation distance, c, of 50 km. This normalised smoothed value, $\stackrel{\mathrm{̃}}{n_i}$, was calculated as follows (Eq. 2):

$$\begin{array}{}\text{(2)}& \stackrel{\mathrm{̃}}{n_i}={\displaystyle \frac{\sum _j{n}_j\cdot \exp \left(\frac{-{\mathrm{\Delta }}_{ij}^{\mathrm{2}}}{c^{\mathrm{2}}}\right)}{\sum _j\exp \left(\frac{-{\mathrm{\Delta }}_{ij}^{\mathrm{2}}}{c^{\mathrm{2}}}\right)}},\end{array}$$

where Δ_ij is the distance between the ith and the jth spatial grid's cell centre and then the summation of the counts, n_j, over j is done considering cells within distance equal to 3 times c (with c being the aforementioned correlation distance) from the ith cell's centre.

Later, Woo (1996) proposed an alternative finite-range form for the kernel, based on the fractal dimension of epicentres and shown in Eq. (3):

$$\begin{array}{}\text{(3)}& \left\{\begin{array}{ll}K(M,x)=\frac{D}{\mathrm{2}\cdot \mathit{\pi }\cdot h\left(M\right)}{\left(\frac{h\left(M\right)}{r}\right)}^{\mathrm{2}-D}& \mathrm{if}r\le h\left(M\right)\\ K=\mathrm{0}& \mathrm{if}r>h\left(M\right)\end{array}\right.,\end{array}$$

where M is the magnitude for a location x; r is the radial separation distance; and h(M) is a magnitude-dependent bandwidth parameter which can be parameterised as $h\left(M\right)=H\cdot \exp (k\cdot M)$, where H and k are regionally estimated constants using seismological and geological considerations and D is the fractal dimension of the epicentres.

Subsequently, Helmstetter et al. (2006) proposed a model for the seismicity density calculation by means of an isotropic adaptive kernel (Izenman, 1991) that smoothed the seismicity depending on the number of events (in order to increase or decrease the detail in the seismicity calculations). Here, the parameter used for the smoothing kernel depends on the average distance between all the events around an earthquake but also on the accuracy of the epicentre location in the first instrumental era of the earthquake catalogue.

Hiemer et al. (2014) created a model based on the seismicity and the fault moment release in order to consider the active mechanisms that generate seismicity in a more direct manner in order to smooth the seismicity (i.e. the locations of the earthquakes). They use the kernel defined by Helmstetter et al. (2006), and, similarly, the fault moment rate was smoothed with an isotropic kernel (Eq. 4):

$$\begin{array}{}\text{(4)}& k_i\left(r\right)={\displaystyle \frac{C\left(d\right)}{{\left(r^{\mathrm{2}}+d^{\mathrm{2}}\right)}^{\mathrm{1.5}}}}\cdot M_i,\end{array}$$

where k is the value of the smoothing kernel for a fault point i, M_i is the fault moment at that point, d is the constant smoothing distance, C(d) is the normalisation constant and r is the epicentral distance.

More recently, the 2020 European Seismic Hazard model (ESHM20) has been released (Danciu et al., 2021). The authors combine the smoothing seismicity algorithms with active fault models. In this case, they point out the challenge of avoiding double-counting events around faults when they consider the background seismicity and the one linked to the fault's activity. Another example of this approach is shown in the work by Pandolfi et al. (2023), where the authors combine 3D information of the seismic sources with the data in the seismic catalogue to calculate the seismic rate.

The works cited in the previous paragraph showcase the importance of considering the active seismogenic sources when computing the activity rate. A common assumption within PSHA is that seismicity can be well described by a Poisson process (Cornell and Winterstein, 1987). A fundamental property of Poisson processes is that the instantaneous rate of events is constant and does not depend upon the occurrence of other events located close in either space or time. However, earthquake sequences feature a significant number of aftershocks, and these events are dependent upon the main shock. The purpose of declustering seismicity data is to remove these dependent events so that the underlying long-term average rate of occurrence can be estimated.

Taroni and Akinci (2021) proposed the use of aftershocks and foreshocks in the seismic activity calculation since removing such events from seismic catalogues may lead to underestimating seismicity rates and, consequently, the final seismic hazard in terms of ground shaking. To do this, they used as kernel a simple weight function of the form (Eq. 5)

$$\begin{array}{}\text{(5)}& k_i\left(N\right)={\displaystyle \frac{\mathrm{1}}{N}},\end{array}$$

where N is the number of events in the seismic series in which the event i belongs.

This weight function ensures that the contribution of each event is the same for the activity rate computation, regardless of its association with a seismic series.

With all the exposed factors, we investigate the sensitivity of the activity rate computation model to both the proximity of the spatial cells to the seismic sources and the epicentral uncertainty related to each event of the catalogue and its influence on a time-dependent seismic hazard. Therefore, we evaluate if the obtained values may be used as a decision factor on operational earthquake forecasting (OEF). Additionally, we consider the foreshocks and aftershocks in order to calculate this activity rate by means of a previous clustering process so each main shock and the corresponding foreshocks and aftershocks are grouped in a given cluster, but we also compare the results with the ones obtained by using a declustered catalogue. To do this, we consider two case studies: Central Italy, a high-seismicity area which helps calibrate the models proposed, and south-eastern Spain, a moderate-seismicity area in which different treatments in the catalogue are tested (declustering and using tectonic b value vs. time-dependent b value).

2 Methodology

In this section, the procedure followed to obtain the parameters that are used inside the smoothing kernel is described. The purpose of this kernel is to smooth the gridded seismic activity (for which a spatial grid is previously defined), helping to improve the description of the seismicity in the area. These parameters define the different models to be tested in the different areas.

2.1 Smoothing kernel parameters

For this work, a modification of the kernel proposed by Frankel (1995) has been used to smooth the gridded seismicity (Eq. 6):

$$\begin{array}{}\text{(6)}& f\left(r\right)=A\cdot \exp \left(-{\displaystyle \frac{(r-\mathit{\mu }{)}^{\mathrm{2}}}{\mathrm{2}\cdot {\mathit{\sigma }}^{\mathrm{2}}}}\right),\end{array}$$

where r is the distance between the centre of the spatial grid cells and the centre of the cell in which the seismic activity is being computed, A is the normalisation constant, μ is the parameter that controls the r value at which the maximum of the function is reached, and σ constraints dispersion of the function around the maximum value.

A more detailed review of the smoothing function, including some examples, can be seen in Sect. 2.1.3.

As expressed in the last paragraph of the introduction, we avoid any arbitrary choice in the definition of these parameters (μ and σ) by assuming they have a geophysical meaning.

2.1.1 Geophysical meaning of the parameter μ

The meaning of this parameter, within the context of the seismic activity smoothing, is the distance from a given cell centre to the point(s) at which the probability of having an earthquake is higher.

It is common to find that the value of this parameter is set to zero (Frankel, 1995; Helmstetter et al., 2006; Hiemer et al., 2014), as the maximum probability of having an earthquake is where it has already occurred before. So the smoothing function has its maximum value in the cell in which the seismic activity rate is being smoothed. This constitutes the first option regarding this parameter: μ=0.

An alternative model is proposed, where the maximum probability is set at the location of the nearest seismic sources. For this to be implemented, the minimum distance between the point in which the seismic activity rate is being computed and the location of the nearest seismic source is calculated and referred to in this work from now on as $d_{{\mathrm{f}}_i}$. So, the second option for the parameter value is $\mathit{\mu }=d_{{\mathrm{f}}_i}$.

For areas in which the tectonic structures are only present in part of the region, a hybrid approach may be applied by using cut-off distance. This cut-off distance may be calculated as follows:

$$d_{\mathrm{c}}=\stackrel{\mathrm{‾}}{d}+\mathrm{2}\cdot {\mathit{\sigma }}_{\mathrm{d}},$$

where d_c is the cut-off distance, $\stackrel{\mathrm{‾}}{d}$ stands for the mean value of the distance between all the structures and σ_d is the standard deviation for all these distances.

If the distance from the centre of the spatial grid cell to the nearest fault is higher than the cut-off distance then μ=0. Otherwise, it is set to $d_{{\mathrm{f}}_i}$.

2.1.2 Geophysical meaning of the parameter σ

This parameter accounts for the dispersion of the values of the distribution around the mean value – that is to say, how far one might expect to find earthquakes around the most probable value (of distance). Therefore, we have considered that this second parameter is related to the accuracy of earthquake's epicentre measurement. This means that it would depend on the methodologies and instrumentation used for the calculation of the epicentre and, thus, on both the year and the location of the catalogue.

It should be noted that σ may depend on other geophysical parameters such as the characteristics of the ground, the style of faulting and/or the tectonic stress regime, to cite a few. Nevertheless, in this work only the influence of the uncertainty in the epicentre's location is considered in the smoothing process.

As in the previous section, two different options regarding the epicentre uncertainty, ε, have been considered: either it depends on the year of occurrence (ε₁) or it is constant and computed as the mean value of the epicentral uncertainty for all the events (ε₂).

Three different models have been proposed to account for the variations in these parameters (Table 1), where ε₁ refers to the different epsilon values depending on the period of the catalogue and ε₂ refers to the fixed value for the whole catalogue.

Table 1 Models for the seismic activity smoothing.

2.1.3 Examples of the smoothing kernel implementation

In this section, some examples of how the smoothing kernel works are shown. Here, the seismic sources are faults in a shallow seismicity context, so the trace of such faults has been considered the location with the maximum probability of having an earthquake. This approach has also been considered for the two case studies in this work. Three main scenarios have been considered to showcase the smoothing kernel.

Usual 1D Gaussian filter, μ=0. This is the case when using models 2 and 3 also when the distance from the centre of the spatial grid cell in which the seismic activity rate is computed to the nearest fault is greater than d_c as defined in the Sect. 2.1.1. An example can be seen in Fig. 1a.

Single fault, μ≠0. When the nearest fault is closer than d_c from the centre of the spatial grid cell, then the resulting function provides a ring-shaped smoothed activity, the width of which depends on σ. Only the section of this ring in which the fault is located is used in the smoothing. This can be achieved by considering the n closest points to the spatial grid cell centre and then computing the angles to define the ring arc (Fig. 1b).

Several faults, μ≠0. This case is a generalisation of the former, with the exception that when the spatial grid cell's centre is in between faults and at similar distances, then the full ring is used as smoothing function (Fig. 1c). On the other hand, if the distance to both faults is similar but the spatial grid cell's centre is not in between the faults, then the resulting smoothing is a ring arc (Fig. 1d).

Figure 1 (a) Smoothing function for μ=0. This can happen when either the distance is greater than d_c or when the spatial cell is over the seismic source. (b) Smoothing function for μ≠0 and a single fault, i.e. when only one seismic source is present and within a distance less than d_c. (c) Smoothing function for μ≠0 when several seismic sources are surrounding the spatial grid cell at similar distances (with angular amplitude greater than or equal to 180°). (d) Smoothing function for μ≠0 in the event that several seismic sources are surrounding the spatial grid cell at similar distances (with angular amplitude less than 180°). The blue lines show the fault tracers from the QAFI database (García-Mayordomo et al., 2012; IGME, 2022). In this example d_c equals 48 km. The stars mark the spatial grid cell considered in each case.

2.2 Cluster identification and seismicity smoothing

First, the spatial grid is defined by creating a rectangle spanning the maximum and minimum longitudes and latitudes of the catalogue with the desired resolution. The choice of the resolution can be motivated by similar studies in comparable tectonic settings or the order of the epicentral uncertainty of the earthquakes in the catalogue. For this work, in the case of Spain, the same resolution as in a previous work in the same area by Montiel-López et al. (2023) has been used. In the case of Italy, although Murru et al. (2016) use a 0.025°×0.025° grid and Gulia and Wiemer (2019) use a 2 km spaced grid, we decided to use the same resolution for both case studies – a 0.015°×0.015° grid. All the data used in this work can be found in the “Data availability” section.

Then, all the events of the catalogue must be assigned to each cell. This is done by calculating the minimum distance of each event to all the spatial grid cell's centres.

One of the most important steps regarding the activity rate calculation in this work is the identification of the seismic clusters present in the area for the selected period of time. As indicated in the introduction, we do not pretend to remove the foreshocks and aftershocks but to identify the main event and all related events in the corresponding cluster.

To do so, even though the epidemic-type aftershock sequence (ETAS) model allows us to assign to each event the probability of being an aftershock (Zhuang et al., 2002; Marsan and Lengliné, 2008; Console et al., 2010), in this work we have decided to select a non-stochastic method based on the performance classifying events of a relevant seismic series.

There are several options for this task: (a) using the Reasenberg and Jones (1989) (RJ) algorithm by applying the ZMAP software (Wiemer, 2001) or (b) using the AFTERAN (A) algorithm (Musson, 1999) or the Gardner and Knopoff (1974) (GK74) declustering algorithms by applying the Python libraries included in OpenQuake's Python scripts (Pagani et al., 2014).

Figure 2 Example of earthquake clusters identified in south-eastern Spain. All the cross markers with the same colour indicate events that belong to the same cluster, while the grey circles are events not belonging to any cluster. The fault traces from the QAFI database (García-Mayordomo et al., 2012; IGME, 2022) are represented by red-coloured lines, and the spatial grid cell's limits are represented by green lines. A zoom in on the Lorca cluster is shown in the bottom-right corner, where the events are plotted using purple circles, the size of which depends on the magnitude of the event.

These algorithms flag each event from each cluster with an identifier which is added as a column to the catalogue, and the events that do not belong to any series have a value of zero for this field. In order to decide which algorithm performs best on the data, a comparison between the RJ, A, and GK74 declustering algorithms has been made using default parameters (see Sect. 3.2.1). An example can be seen in Fig. 2. The smoothing procedure explained by Taroni and Akinci (2021) has been adapted to work with as many clusters as can be found in each spatial grid cell.

To do this, all the events belonging to each cluster are counted and define the cluster weight, c_j. If an event does not belong to any cluster (i.e. the cluster label for that event is set to zero), then c_j equals 1. The weighted counts for each spatial grid cell are calculated as the summation of all the events over the different clusters (Eq. 7):

$$\begin{array}{}\text{(7)}& k_i=\sum _{\mathrm{1}}^j\left({\displaystyle \frac{\mathrm{1}}{c_j}}\sum _{\mathrm{1}}^m\mathrm{1}\right),\end{array}$$

where k_i is the weighted count of events inside the cell i; the first sum goes over the number of clusters, j, inside the ith cell; and the second summation goes over each event, m, of the cluster j that is inside the spatial grid cell i.

For instance, if inside a cell there are 20 events that belong to a cluster composed of 100 events in total and 13 events that do not belong to any cluster, the weighted number of events for that cell is $k_i=\mathrm{13}+\mathrm{20}/\mathrm{100}=\mathrm{13.2}$.

2.3 Seismic activity rate computation

After the smoothing process explained before, we consider both the nature and source of the earthquakes and the uncertainty related to the earthquake location.

Thus, the seismic activity rate is calculated as the product of the weighted counts and the smoothing kernel values (Eq. 8):

$$\begin{array}{}\text{(8)}& {\mathit{\lambda }}_i={\mathbf{w}}_i\cdot \mathit{k},\end{array}$$

where w is an n×n matrix, with n being the number of cells of the spatial grid, that for each cell, i, contains the n values of the smoothing kernel associated with the cell. On the other hand, k is a vector containing the weighted count for each cell i as defined in Eq. (7). So, for each cell the vector product between the smoothing kernel (w_i can be seen as a vector) and the weighted count is calculated. This means all cell counts are added in each cell activity rate computation, and the smoothing function works similarly to the correlation distance proposed by Frankel (1995).

2.4 Exceedance probability computation

The annual exceedance probability of a given PGA has been obtained by developing a Python script based on OpenQuake (Pagani et al., 2014). Two models have been tested for the computation of the needed b value: a fixed (time-independent) b value assigned from the tectonic zones of each country and a gridded (time-dependent) b value calculated using the methodology proposed by Montiel-López et al. (2023). In both cases, another Python script was developed to obtain a gridded time-dependent seismic activity rate for moment magnitudes greater than or equal to 4.0 M_w for each cell.

The Akkar and Bommer (2010) empirical equation has been used as a ground motion prediction model since it is appropriate for Mediterranean regions such as Spain, whereas the Akkar et al. (2014) ground motion empirical equation has been used for Central Italy since it was used in the EHSM20 (Danciu et al., 2021). An example of the workflow for this work is represented in Fig. 3.

Figure 3 Workflow diagram for the exceedance probability computation. Examples of the main inputs are given as the three spatial mappings of the (a) b value, (b) seismic activity rate and (c) tectonic zones' b value.

Our goal is to investigate if the changes in the seismic activity and b-value time series can be observed as a trend in the PSHA results. In the event that such trends are observed, this methodology could be used for OEF. For this reason, the temporal evolution of the annual probability of exceedance (PoE) of a background PGA corresponding to a 475-year return period (i.e. 0.002 PoE) has been computed as a time-dependent value. The results have been expressed as a relative change (RC in percentage change, Eq. 9) between background annual exceedance probability (long-term value) and the time-dependent annual exceedance probability. Depending on the country, the background (long-term) PGA value may have been updated in the corresponding seismic hazard studies. This could be due to the occurrence of new damaging earthquakes or improvements in the seismic knowledge of a region, amongst other reasons. In case such changes have been made, a new background PGA value has been computed using the data up until the year the seismic hazard information was updated.

$$\begin{array}{}\text{(9)}& \mathrm{RC}=\mathrm{100}\cdot \left({\displaystyle \frac{\mathrm{PoE}}{\mathrm{0.002}}}-\mathrm{1.0}\right)\end{array}$$

In order to save computation time, the annual exceedance probability is only calculated for a selection of cities located inside the spatial grid.

Additionally, we have computed the annual variation of the RC in the exceedance probability (${\mathrm{RC}}_i-{\mathrm{RC}}_{i-\mathrm{12}}$, with i the computed month) and the monthly variation (${\mathrm{RC}}_i-{\mathrm{RC}}_{i-\mathrm{1}}$) to investigate if any of these metrics is effective as an indicator for OEF.

3 Case studies

As explained before, the goal of our smoothing methodology is to test the viability of producing time-dependent seismic hazard results which may be used for making decisions before the main earthquake. Therefore, now we present and discuss the results obtained for two different regions with different seismic behaviour – Central Italy (high seismicity) and south-eastern Spain (low seismicity). We check if there are significant changes in the defined metrics before the occurrence of important earthquakes carrying out a retrospective validation of how useful the results are.

3.1 Central Italy

3.1.1 Catalogue preparation and parameters for computation

As mentioned before, since Central Italy is a very active region, this case study will help us to decide which of the models (Table 1) performs better. Central Italy (Abruzzo, Campania, Lazio, Marche, Molise, Toscana and Umbria) is a macro-region where several high-magnitude earthquakes and significant seismic series have occurred in the past. The main focus is on L'Aquila, where a 6.29 M_w earthquake (Table 2) struck the area in 2009 and caused 309 deaths and 1500 injuries. Therefore, the city of L'Aquila has been selected as the site for the hazard computation.

Table 2 The L'Aquila earthquake data and distance to hazard computation site.

Figure 4 shows the location of the area of study (a rectangle with longitudes from 11.394 to 15.391° E and latitudes from 40.359 to 44.353° N) and the tectonic zones and main faults as defined by Danciu et al. (2021) to compute the European Seismic Hazard Map. For this area, the Italian HORUS (Lolli et al., 2020) catalogue has been used as it has been homogenised and comprises events from 1960 to 2023. This catalogue is a homogenised instrumental catalogue based on the hypocentral location of earthquakes compiled from the Italian Seismological Instrumental and Parametric Database (ISIDe) (ISIDe Working Group, 2007). A spatial filtering process has been applied to the catalogue to extract the events within the above-mentioned area, so 49 112 events with maximum depth of 30 km and maximum magnitude of 6.81 M_w remain.

Figure 4 The map shows the tectonic zones (green lines) in Central Italy with their acronyms and tectonic b values as can be found in Danciu et al. (2021). The star marks the epicentre of the L'Aquila earthquake (Table 2), and the red lines represent the fault traces (Basili et al., 2022).

The events that are not earthquakes (such as quarry blasts, eruptions and explosions) have been filtered out from 2012 on (as the catalogue has such information). In order to consider the influence of such events prior to 2012, the area that has been selected for this study does not show important changes in the b value according to the results of Taroni et al. (2022). In this case the catalogue has not been declustered, but the clusters have been identified by using the GK74 algorithm, and this information has been used to weight down the influence of the non-independent events towards the seismic parameters' computation.

A spatial cell grid of 0.015°×0.015° spanning the above longitude and latitude ranges has been created (using 70 756 points). The completeness magnitude (Table 3) has been retrieved from Taroni et al. (2021).

Table 3 Completeness magnitude values proposed by Taroni et al. (2021).

Since the σ parameter used in the smoothing kernel computations and based on the location uncertainty aims to account for the physical variability in the location of the earthquakes, three models with different uncertainty values have been tested to showcase the variability in the results as a consequence of increasing or decreasing the uncertainty of the epicentre location. In order to obtain such values, the work from Scudero et al. (2021) gives insight into the variation of the horizontal error (ERH) in Italy as well as a range of mean values for different revision processes of the data (2.2, 3.3 and 13.1 km). Given that the HORUS (Lolli et al., 2020) catalogue, used in this work, has no information on the ERH but the locations of the events are obtained through the ISIDe database, their spatial uncertainty can be deduced from the CPTI15 catalogue (Rovida et al., 2020, 2022). The aforementioned range of mean values for the ERH is coherent with the mean spatial uncertainty obtained from the CPTI15 catalogue.

Therefore, a minimum value of 6 km, in agreement with the previous explanation, and a maximum value of 30 km, following the work of Taroni et al. (2021), have been chosen to characterise the spatial uncertainty. The three models proposed for the seismic activity smoothing are presented in Table 4.

Table 4 Models for the activity smoothing in Central Italy.

^* From Taroni et al. (2021).

3.1.2 Results

Figure 5 (Model 1t) presents a moderate increase in the annual exceedance probability (25 %) 1 month before the L'Aquila earthquake occurred, and not only the annual but also the monthly variations of relative change reach values higher than 35 %. This sudden change is due to the foreshock activity that preceded the main shock, as a 4.1 M_L earthquake occurred on 30 March 2009. Figure 6 (Model 2t) shows a trend in all the metrics that is similar to the previous model with a slightly lower value for the exceedance probability change before the earthquake (22 %) and the annual and monthly variations (32 %). The Model 3t (Fig. 7) provides the lowest values for the metrics (−3 %, 4 % and 3 %, respectively). After this increase the RC slowly decreases over time for all three models.

Figure 5 Model 1t. Relative change (RC) of the annual exceedance probability (a) and its annual and monthly variation (b). The vertical dashed black line marks the occurrence of the L'Aquila earthquake.

Figure 6 Model 2t. Relative change (RC) of the annual exceedance probability (a) and its annual and monthly variation (b). The vertical dashed black line marks the occurrence of the L'Aquila earthquake.

Figure 7 Model 3t. Relative change (RC) of the annual exceedance probability (a) and its annual and monthly variation (b). The vertical dashed black line marks the occurrence of the L'Aquila earthquake.

Given that the main objective is to be able to perform OEF in the area of study, Model 1t (Fig. 5) has been selected since it performs better in terms of exceedance probability change. Moreover, its annual and monthly variations are the highest 1 month before the main shock (when compared with the models 2t and 3t).

3.2 South-eastern Spain

3.2.1 Catalogue preparation and parameters for computation

The south and south-east of Spain are regions with a higher seismic hazard in Spain (IGN-UPM Working Group, 2013; Kharazian et al., 2021), with values reaching 0.23 g for a 10 % probability of exceedance in 50 years (i.e. return period of 475 years). Although Spain is a moderate- to low-seismicity region compared to other European countries such as Italy or Greece, it has been exposed to several damaging earthquakes in the past, with the 1829 Torrevieja earthquake and the 1884 Arenas del Rey earthquake being the most representative, both with a maximum intensity of IX-X. Additionally, in the last 25 years, south-eastern Spain has suffered seven earthquakes with M_w greater than or equal to 4.5 (Table 5 and Fig. 8 present only those classified as main shocks and located in the area of study), with the 2011 Lorca earthquake being the most relevant since it was the most recent earthquake causing damage to buildings and injuries to the population. The seismicity is usually very shallow (mainly lower than 10 km). Three main cities (Murcia, Lorca and Vera from north to south) have been chosen as representative of the region in terms of decreasing seismic hazard values for a 475-year return period.

Figure 8 Map showing the tectonic zones in the south-eastern Spain region with their corresponding b values (black numbers) as computed by García-Mayordomo (2015) and IGME (2015). The red triangles mark the earthquakes and the red numbers the order in which they appear in Table 4. The orange circles show the location of the main sites for the seismic hazard analysis.

Table 5 Damaging earthquakes in the last 25 years and epicentral distance to some chosen cities in the area of study.

In order to compute the seismic activity rate to be used in a PSHA, first we need to compile a homogeneous and complete seismic catalogue in the influence area, which is needed for the chosen locations. This catalogue, obtained from the Instituto Nacional de Geografía web page (IGN, 2022), comprises all the events from 1396 to August 2023 in south-eastern Spain inside the tectonic zones of Eastern Betics Shear Zone (ZCBOR), Eastern Inner Betics (BIOR), Valencian Plateau and Alicante's Prebetic (PVPA), Murcian Prebetic (PM), Sierra Nevada–Filábrides and Guadix–Baza (SNFCGB), Central Inner Betics (BIC), Southern Plateau (MS), Cazorla–Segura and Albacete's Prebetic (CSPA), and Central Guadalquivir and Algerian-Balearic Basin (CAB), as defined by García-Mayordomo (2015) to create the Spanish seismic hazard map (IGN-UPM Working Group, 2013).

The catalogue contains a total of 20 279 events that span from 1396 to August 2023. Their moment magnitudes range from 0.1 to 6.8 after being homogenised using the magnitude correlation equations for this region (IGN-UPM Working Group, 2013). Their depth goes up to 90 km, although in the calculations only the earthquakes shallower than 30 km are considered (which amount to a total of 20 168 events). A spatial grid of 0.015°×0.015° covers the area of study (using 40 401 points).

Table 6 presents the number of clusters and the events in clusters for the whole seismic catalogue. The RJ algorithm classifies a total of 652 clusters in the catalogue, while GK74 detects 1012 clusters. The A algorithm identifies 1245 clusters. As can be seen, despite the three methods relying on windows for their calculations, there are significant differences in the results, not only in the number of clusters but also in the number of events inside each cluster.

Table 6 Comparison of cluster identification and total events inside clusters among three declustering algorithms: an analysis using Lorca's seismic series.

Zaliapin and Ben-Zion (2020) pointed out these problems with the identification of aftershocks and main shocks and proposed an algorithm to discriminate between background and clustered events by randomly thinning a complete catalogue by removing nearest-neighbour earthquakes. Moreover, Anderson and Zaliapin (2023) examine the effect on the hazard estimation when using different declustering thresholds. They conclude that hazard estimates are most sensitive to the catalogue thinning near the aftershock zone and less sensitive elsewhere.

In spite of the difficulties in defining the clusters, Caban̈as et al. (2011) carried out a detailed study on the 2011 Lorca earthquake seismic series. This study is the best definition at the moment, so we use it to validate the best algorithm. Caban̈as et al. (2011) identified 146 events (including the foreshock, the main shock and the aftershocks) that belong to Lorca's series, from 11 May until 19 July 2011. With this information, in order to test the performance of the declustering methods, the confusion matrices for each one have been computed. In the area of study, a total of 249 events have been recorded, which means a total of 103 background events should be identified. For this analysis, all the events classified in a cluster different from the one of Lorca's series have been considered background for simplicity. Figure 9 shows that the GK74 method is the most adequate (with a 94.43 % mean for the metrics compared with the 92.88 % for RJ and a 78.79 % for A) and also the one that is able to identify more events belonging to Lorca's series.

Figure 9 Confusion matrices for the tested declustering methods. Inside each square, the number of events (bold) and some metrics computed using the data are presented (NPV stands for negative predictive value).

The catalogue starts in the historical period (when there was a lack of instrumentation and procedures to accurately locate the epicentres and evaluate the magnitude of the earthquakes) and ends in the present day. This implies that not all the magnitude values are complete in the catalogue (low magnitudes are missing in the historical period) and the location uncertainty also differs depending on the year of detection. First, we characterise the completeness magnitude – the minimum magnitude from which the catalogue is not missing any record – and periods for the Spanish seismic catalogue.

Gaspar-Escribano et al. (2015) defined different threshold magnitudes for different regions around Spain. The class marks of these magnitude intervals for the zone of interest (south-eastern Spain) have been selected as the completeness magnitudes up until 1962. From 1962 on, the completeness magnitude has been computed by spatially averaging the gridded completeness magnitude results available from González (2017) over the area of study. The values used in this work are presented in Table 7.

Table 7 Completeness magnitude for each period according to Gaspar-Escribano et al. (2015) on top and the spatially averaged completeness magnitude for each period using the results from González (2017) on the bottom.

The uncertainty of the epicentral location (ε) varies with time, showing a decreasing behaviour since the techniques and instrumentation have continuously been improved. The appropriate estimation of this uncertainty is very important in order to correctly assign the location of each earthquake to a given seismic source.

Following the research of Peláez and López Casado (2002), the ε values for each period are presented (Table 8). The value corresponding to the period 1990–2023 has been obtained as the average epicentral uncertainty using the data provided by the national seismic network. A second fixed ε value of 7.5 km has been computed as the mean value for all the uncertainties in the catalogue.

Table 8 ε values proposed by Peláez and López Casado (2002).

^* Calculated as the average epicentral uncertainty for the 1990–2023 period events in our catalogue.

The three models to be evaluated are presented in Table 9. The fixed model implies a fixed b value (the b value from the tectonic zone), while the time-dependent model indicates a time-dependent b value.

Table 9 Models for the exceedance probability calculation in south-eastern Spain.

Figure 10 Relative change (RC) of the annual exceedance probability for models 1f, 2f, and 3f. The vertical dashed grey lines mark the earthquakes considered in Table 5 which are closer than 75 km to each one of the sites.

3.2.2 Results

After computing the time-dependent PSHA for the different models shown in Table 9, we have observed that although all the graphs show similar behaviour for the RC, Model 1t provides greater annual and monthly variations of the RC for some of the earthquakes compared to the rest, similarly to Italy's case study. Therefore, with the exception of the fixed models (models 1f, 2f and 3f) where we present a general comparison between them, we present the results of Model 1t in this section along with figures comparing all three models. The stand-alone figures for the models 2t and 3t can be found in Appendix A.

Time-dependent PSHA using tectonic zones' b values (fixed model). The time-dependent PSHA (PGA for a return period of 475 years) has been computed using the proposed methodology, for the compiled non-declustered catalogue, in 1-month increments starting from 1990. The b value is constant and given by the zonation proposed by García-Mayordomo (2015) (models 1f, 2f and 3f). This background PGA value is a long-term PSHA which varies each time that the seismic normative changes in Spain. Our first background PGA corresponds to the PSHA computed using a catalogue with the same length as the one used for the NCSE-94 (España. Ministerio de Obras Públicas, Transportes y Medio Ambiente and Colegio Oficial de Arquitectos de Andalucía Oriental, 1994). This background PSHA value is used from 1990 to December 1998 (since the next code updated the seismic hazard map using a seismic catalogue up to 1999). The second background value corresponds to the PSHA computed with a catalogue of the same length as the one used for the NCSE-02 (España. Dirección General del Instituto Geográfico Nacional, 2002) and is used from 1999 to May 2011 (since that is the year when the seismic hazard map was updated again). Finally, the last background value corresponds to PSHA computed using a catalogue of the same length as the one used for the current seismic hazard map for Spain (IGN-UPM Working Group, 2013). This last value is used from June 2011 to August 2023.

Figure 10 represents the temporal evolution of the results for each tested model. The vertical lines correspond to the main earthquakes from Table 5 and with an epicentral distance less than 75 km from the chosen city, since more distant events would not require forecasting as they are not expected to cause damage. As can be seen, the behaviour is similar for the three models. The exceedance probability decreases continuously since 1990 except for Lorca in 1997 and 2006, Murcia in 1996, and Vera in 1994 and 1999. These variations are due to seismic activity changes, but they do not appear to be related to the occurrence of any of the main earthquakes from Table 5. On the other hand, all the models provide similar changes in the exceedance probability, although in the city of Lorca, Model 3f is the one with the lowest percentage change and Model 2f is the one with the highest percentage change.

Figure 11 Relative change (RC) of the annual exceedance probability and corresponding uncertainty for Model 1t in Lorca, Murcia and Vera (from top to bottom). The vertical dashed black lines mark the earthquakes considered in Table 5 which are closer than 75 km to each one of the sites. A zoom in on the mentioned increase in the RC during 1998 at the Lorca site appears in the upper-left side of the graph.

Figure 12 Mean value of the relative change (RC) of the annual exceedance probability for models 1t, 2t and 3t in Lorca, Murcia and Vera (from top to bottom). The vertical dashed red lines mark the earthquakes considered in Table 5 which are closer than 75 km to each one of the sites.

It seems that the use of a constant b value coupled with a time-dependent seismic activity rate leads to a RC that decreases over time with a constant rate and sudden increases mainly due to changes in the background PGA values. This behaviour is due to the update of this parameter depending on the period of the catalogue, which is increased by a month at a time (as it was the selected minimum time step for this area). We find that this approach is not appropriate towards earthquake forecasting for areas with low to moderate seismicity. This uniform behaviour potentially rules out the possibility of finding any metrics for OEF.

Time-dependent PSHA using the time-dependent b value (time-dependent model). In this section, models 1t, 2t and 3t (Table 9) are tested using the three PGA background values explained previously (and computed in January 1990, December 1998 and May 2011). As can be seen in Fig. 11, the annual probabilities decrease before the Mula earthquake for the Lorca site. However, close to the occurrence of the earthquake, the RC shows a slight increase even at the Vera site, although it is 101.4 km away from the earthquake's epicentre. At the Murcia site, the RC continuously decreases until 5 months before the earthquake, when it shows a sharp increase from −75% to −60 % in the change in exceedance probability. This change is also seen in the annual and monthly variations of the RC (Fig. 11 zoom in). After the Mula earthquake the change in probability exceedance remains higher than 20 % (even increasing until 50 % in the case of Lorca and 100 % in the case of Murcia) for both the Lorca site and the Murcia site until the Lorca earthquake happens. At the Vera site, this parameter oscillates about the baseline. After 2011, the RC steadily increases in Vera, whereas in Lorca and Murcia it stays constant after 2019. On the other hand, Model 2t and Model 3t (Fig. 12) show a similar behaviour to Model 1t. It should be noted that the sudden changes in the RC in January 1999 and 1 month after the Lorca earthquake, i.e. the −60 % to 0 % increase in January 1999 (for both Lorca and Murcia) and the 100 % to 0 % decrease in June 2011 (for both Lorca and Murcia), are artefacts due to the change in the background PGA and cannot be considered in the analysis. However, the annual and monthly variations of the RC shown in Figs. 13 and 14 allow us to see the changes related to the aforementioned earthquake occurrences.

Figure 13 Annual and monthly variations of the relative change in the annual probability of exceedance for Model 1t in Lorca, Murcia and Vera (from top to bottom). The vertical dashed black lines mark the earthquakes considered in Table 5 which are closer than 75 km to each one of the sites.

Figure 14 Monthly variations of the relative change in the annual probability of exceedance for models 1t, 2t and 3t in Lorca, Murcia and Vera (from top to bottom). The vertical dashed red lines mark the earthquakes considered in Table 5 which are closer than 75 km to each one of the sites. The earthquake that could cause the peak in 1994 at the Vera site has also been indicated.

Figure 13 shows a 20 % mean decrease in the annual variation of the RC from January 1991 to March 1993 that could be explained by the RC uncertainty (as can be seen in Fig. 11 for both the Lorca site and the Murcia site, with higher uncertainty for that period). Then, a 15 % increase in the annual variation can be seen in the RC from October 1998 until August 1999 for the Lorca site (3 months before the Mula earthquake and then 6 months after it). This increase can also be seen in the declustered catalogue scenarios (Figs. 16, A5 and A6). At the Vera site, the increased RC variation during 1999 could be due to changes in the b value from April to December (six earthquakes with magnitude from 3.5 to 3.8 M_w occurred). In the case of the Gergal and Bullas earthquakes in 2002, an increase in the variation of the RC cannot be observed. However, it can be seen for the Lorca site and the Murcia site that from July 1999 to May 2001 the annual variations of the RC reach values higher than 15 % with respect the baseline and with a mean value of 10 % over this period. These increased values cannot be related to any close seismic activity greater than or equal to 4.0 M_w. It can also be seen that at the Lorca site the annual variation stays higher than 20 % for 1 year after the Lorca earthquake. Lastly, the peak in the annual and monthly variation in Vera in 2022 appears due to the seismic activity in Turre (a town 14 km south from Vera) where a 4.0 M_w earthquake struck on 31 December 2022.

Similar results are obtained for all three locations for the considered models regarding the monthly variations (Fig. 14). Overall, the monthly variations do not show changes preceding relevant earthquakes for this case study. One of the possible explanations is the lack of foreshocks in most of the main shocks. In the Lorca earthquake, even though there was a 4.5 M_w earthquake almost 2 h before the main shock, the 1-month increments in the computation process are not able to show any change in RC.

The annual variations (Fig. 15), on the other hand, show periods of increased RC before some of the selected earthquakes. An example is seen at the Lorca site, where a 15 % increase is seen before the Mula earthquake from June 1998 (the earthquake occurred in February 1999). Another example can be seen at both the Murcia site and the Lorca site, where a 10 % increase can be seen before Aledo earthquake from May 2004 until the earthquake occurrence in January 2005. For this metric, differences between the three models can be seen. For instance, Model 1t and Model 2t show greater changes after the Lorca earthquake at the Lorca site.

Figure 15 Annual variations of the relative change in the annual probability of exceedance for models 1t, 2t and 3t in Lorca, Murcia and Vera (from top to bottom). The vertical dashed red lines mark the earthquakes considered in Table 5 which are closer than 75 km to each one of the sites.

The most prominent increase on the annual variation occurs after the Lorca earthquake (32.8 %, 36.2 % and 21 % for models 1t, 2t and 3t).

Effect of the declustering on the results. In order to compare the effect of the catalogue declustering on the results, Model 1t has been plotted using both the declustered catalogue (with a total of 13 841 events) and the full catalogue (with the clusters identified and weighted down accordingly).

Figure 16 represents the changes in the annual exceedance probability for Model 1t. As can be seen, the results using a non-declustered catalogue provide lower changes in the exceedance probability for M_w 4.0 at the Lorca site from 1996 until 2011. Then, from 2011 until 2023, the non-declustered catalogue provided a higher RC. At both Murcia and Vera, the results are similar for the non-declustered and declustered catalogues. It should also be noted that the mean uncertainty of the RC is slightly higher for the declustered catalogue (11.21 % for Lorca, −0.41 % for Murcia and 5.58 % for Vera). Since the results are compatible, keeping the foreshocks and aftershocks, i.e. using the non-declustered catalogue, seems to be a better choice if the aim is to perform OEF. Some of the advantages would be a lower uncertainty in the RC and the possibility of using more detailed timescales in case foreshocks are present.

Figure 16 Model 1t. Comparison of the relative change (RC) of the annual exceedance probability and its uncertainty for a non-declustered and a declustered catalogue in Lorca, Murcia and Vera (from top to bottom). The vertical dashed black lines mark the earthquakes considered in Table 5 which are closer than 75 km to each one of the sites.

4 Conclusions

This methodology considers the influence of all the events in the seismic clusters and also the location of the seismic sources (corresponding active faults) for seismic activity rate smoothing and b-value computation, showing that when computing a time-dependent PSHA the use of a non-declustered catalogue provides similar results to using a declustered catalogue with the added benefit of keeping the foreshock activity. Therefore, if we compute the changes in the annual probability of exceedance for a given PGA value (fixed as a background value which may change according to the updates in the seismic normative), we are able to show how this probability is changing with time.

The changes in the annual probability of exceedance (increases and decreases) can be more accurately described using a spatially gridded time-dependent b value instead of a fixed one for each tectonic zone. This can be seen when comparing Fig. 10 with Fig. 12. Therefore, we suggest using spatially gridded b values for the corresponding period (time-dependent) when computing the background PGA value and the corresponding changes in the annual probability of exceedance in the time-dependent PSHA.

Regarding which of the proposed models can be more effectively used to describe these changes, we have to consider several factors. One could be how close the computed PGA values are to the national seismic hazard maps for each country. In the case of Central Italy, models 1t, 2t and 3t provide the following background PGA values: 340.61, 359.72 and 334.28 cm s⁻². The ESHM20 model (Danciu et al., 2021) computes 334.38 cm s⁻². The closest match would be Model 3t followed by Model 1t. However, by looking at Figs. 5 and 7, it can be seen that Model 3t seems to be less affected by changes in the seismic activity than Model 1t, as the monthly and annual RC variations suggest (Model 1t monthly and annual variations are 4.5 times higher than Model 3t variations). With this information Model 1t seems appropriate for the purpose of this work.

In general, this methodology benefits from complete catalogues in zones with increased seismicity – assuring less uncertainty in the b-value computation – and well-defined seismicity sources, where the seismicity smoothing is accurate. Figure 16 shows this result, as the non-declustered catalogue (with weighted-down cluster events) has less RC uncertainty and enables the use of the foreshocks in daily to weekly timescales.

Although our results are not significant enough to relate these changes to the occurrence of a main earthquake for low- to moderate-seismicity areas, the methodology can be useful for other countries with a higher seismicity or in the future if new significant earthquakes occur in the studied region of Spain. As we saw, for Central Italy both the annual and monthly changes in the exceedance probability show important variations related to the foreshock activity preceding the L'Aquila earthquake. This could be useful for OEF.

Finally, in the case of south-eastern Spain, the relative change in the annual probability of exceedance remained high in the region after the Mula earthquake and did not decrease until the occurrence of the Lorca earthquake. However, the continuous increase in this parameter in Vera after the Lorca earthquake cannot be directly related to a potential upcoming earthquake similar to the one from Lorca. Therefore, more time and data are needed to confirm this.

Appendix A: Other results for south-eastern Spain

A1 Time-dependent PSHA using the time-dependent b value (time-dependent model)

Figure A1 Relative change (RC) of the annual exceedance probability and corresponding uncertainty for Model 2t in Lorca, Murcia and Vera (from top to bottom). The vertical dashed black lines mark the earthquakes considered in Table 5 which are closer than 75 km to each one of the sites.

Figure A2 Relative change (RC) of the annual exceedance probability and corresponding uncertainty for Model 3t in Lorca, Murcia and Vera (from top to bottom). The vertical dashed black lines mark the earthquakes considered in Table 5 which are closer than 75 km to each one of the sites.

Figure A3 Annual and monthly variations of the relative change in the annual probability of exceedance for Model 2t in Lorca, Murcia and Vera (from top to bottom). The vertical dashed black lines mark the earthquakes considered in Table 5 which are closer than 75 km to each one of the sites.

Figure A4 Annual and monthly variations of the relative change in the annual probability of exceedance for Model 3t in Lorca, Murcia and Vera (from top to bottom). The vertical dashed black lines mark the earthquakes considered in Table 5 which are closer than 75 km to each one of the sites.

A2 Effect of the declustering on the results

Figure A5 Model 2t. Relative change (RC) of the annual exceedance probability for a non-declustered and a declustered catalogue in Lorca, Murcia and Vera (from top to bottom). The vertical dashed black lines mark the earthquakes considered in Table 5 which are closer than 75 km to each one of the sites.

Figure A6 Model 3t. Relative change (RC) of the annual exceedance probability for a non-declustered and a declustered catalogue in Lorca, Murcia and Vera (from top to bottom). The vertical dashed black lines mark the earthquakes considered in Table 5 which are closer than 75 km to each one of the sites.