Dependence of seismic hazard assessment on the observation time interval: insights from a synthetic earthquake catalogue in southeastern Spain

Traditionally Probabilistic Seismic Hazard Assessment (PSHA) is mainly based on historical and instrumental catalogues. The seismic catalogue in SE Spain covers a period of almost 1 kyr. Although this catalogue can be considered long, its time span is not enough to record the complete seismic cycle of the slow moving regional faults, for example, the ones which conform the Eastern Betics Shear Zone (EBSZ). To assess whether the seismic hazard at the EBSZ depends on the time interval in which the earthquake catalogue has been recorded, this study performs a PSHA using synthetic seismicity. The synthetic seismicity in the EBSZ consists of a 1 Myr catalogue generated in previous studies using the RSQSim earthquake simulator, applied to the fault system that forms the EBSZ. This catalogue has been divided into ten thousand sub-catalogues of the same duration as the historical and instrumental one and randomly distributed over time. The magnitude-frequency distributions of the synthetic sub-catalogues show significant variability in the maximum reached magnitude, in the slope of the Gutenberg-Richter relationship and in the annual rate of earthquakes. A hundred sub-catalogues have been selected to perform individual PSHA and have been compared with the results derived from historical and instrumental seismicity. Using R-CRISIS software, this study obtains seismic hazard curves for the main cities in the region, showing the estimated return period for different values of Peak Ground Acceleration (PGA). The hazard curves reveal that each sub-catalogue leads to different return period values. The obtained variability ranges from 11 % to 21 % for a PGA=0.04 g, and from 25 % to 58 % for a PGA=1 g. Our results show that there is a dependence between seismic hazard and the observation time interval in which an earthquake catalogue is recorded.

1 Introduction

Consequences of large earthquakes affect all aspects of human life (Daniell et al., 2012; Marano et al., 2010), both by direct and indirect effects (Gokhale et al., 2004). Characterization of seismic hazard allows us to prevent and mitigate the devastating effects of earthquakes. The development of a Probabilistic Seismic Hazard Assessment (PSHA) through a statistical analysis of the occurrence of earthquakes, leads to the determination of exceeding probabilities for different ground motion levels (Cornell, 1968; McGuire, 2008). The results derived from a PSHA are used in the definition of seismic actions in building codes (e.g. Bisch et al., 2012; Barbat et al., 2005), which makes them of crucial importance in the protection of human lives and the limitation of damages.

Recorded seismicity is the fundamental input data of PSHA, in which the so-called Gutenberg-Richter relationship (Gutenberg and Richter, 1944) and the maximum magnitude of the earthquake catalogue describe the seismic activity of a region. Gutenberg-Richter relationship is calculated fitting to a log-linear equation the cumulative magnitude distribution of an earthquake catalogue (Eq. 1).

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

where N(M) represents the number of earthquakes with magnitude equal to, or higher than M, a is the annual activity rate of earthquakes and b is slope of the function, which depends on the relationship between the number of large and small events.

The probability density function associated with the Gutenberg-Richter relationship (f_M(M), Eq. 2) describes a dependence with the minimum magnitude. This magnitude can be considered as the magnitude of completeness, M_c, from which all seismicity is assumed to have been recorded.

$$\begin{array}{}\text{(2)}& f_M\left(m\right)=\mathrm{1}-e^{-\mathit{\beta }(M-M_{\mathrm{c}})}{M}_{\mathrm{c}}\le M\le \mathrm{\infty }\end{array}$$

In Eq. (2), β represents another way of expressing the slope of the Gutenberg-Richter relationship (Eq. 3).

$$\begin{array}{}\text{(3)}& \mathit{\beta }=b\cdot \ln \left(\mathrm{10}\right)\end{array}$$

Most of the implementations of the Gutenberg-Richter equation on PSHA are by means of the truncated, or modified, relationship (Cornell and Vanmarcke, 1969; Cosentino et al., 1977), in which probability density function ($f_M^t\left(M\right)$, Eq. 4) depends also on M_max. In this case, M_max is considered to be the maximum possible value of the earthquake magnitude. This approximation avoids the overestimation of the mean return period of the largest earthquakes of the general Gutenberg-Richter exponential relationship.

$$\begin{array}{}\text{(4)}& f_M^t\left(m\right)={\displaystyle \frac{\mathrm{1}-e^{-\mathit{\beta }(M-M_{\mathrm{c}})}}{\mathrm{1}-e^{-\mathit{\beta }(M_{\max }-M_{\mathrm{c}})}}}\end{array}$$

At slow-moving fault systems, with tectonic deformation rates of less than 2 mm yr⁻¹ and seismic cycles of thousands of years, recorded seismicity is usually moderate (Echeverria et al., 2013; Shi et al., 2014; Anderson and Luco, 1983; Perea et al., 2006; Alfaro et al., 2012; Gómez-Novell, 2021). An example is the case of SE Spain. In this region, seismicity is low to moderate, even though the EBSZ fault system has the potential to generate large earthquakes (Perea et al., 2012; Gómez-Novell et al., 2020; Herrero-Barbero et al., 2021).

Paleoseismological studies have shown that the time interval between major earthquakes at the EBSZ single faults is of the order of several thousands of years (Ortuño et al., 2012; Martínez-Díaz et al., 2012b, 2019; Masana et al., 2010, 2018; García-Mayordomo, 2005; Canora-Catalán et al., 2016; Ferrater et al., 2015; Moreno, 2011; Roquero et al., 2019; Insua-Arévalo et al., 2015; Martín-Banda et al., 2016). Meanwhile, the historical and instrumental record used in hazard studies at SE Spain, barely covers a thousand years: from 1048 to the present (IGN and UPM, 2013; González, 2017). This period is shorter than the seismic cycle of the faults. For example, in the case of the Carrascoy fault, the interval between large earthquakes is estimated to be of three thousand years (Martín-Banda et al., 2016). Consequently, historical and instrumental record covers just one-third of the seismic cycle at this fault. Meaning a narrow part of the seismic cycle, earthquake catalogues neither may record the entire seismogenic potential of faults nor changes in their behaviour. An example of the latter is the study carried out by Martín-Banda et al. (2021) in which fault slip rate variability has been proposed during the last 210 kyr in the case of the Carrascoy fault at the EBSZ. This variability can be related to accelerations of the system and interactions between adjacent faults, which can lead to the coexistence of short seismic cycles with large ones (Sieh et al., 2008; Benedetti et al., 2013).

Characterization of slow-moving fault systems uses geological data to complement the scarce information about seismic activity. A better understanding of fault systems has turned into the inclusion of them as independent seismogenic sources in hazard studies around the world (Field et al., 2015; Stirling et al., 2012; Valentini et al., 2017; Peruzza et al., 2011; García-Mayordomo et al., 2012a). According to García-Mayordomo (2007), at cities located near active faults with recurrence intervals of 2 kyr, the inclusion of faults as seismogenic sources may represent a contribution of 10 %–30 % of hazard for return periods of 500 yr. Another approach has modeled the EBSZ as seismogenic source zone. To include the contribution of its main faults to a classical PSHA, a hybrid approach was used by defining an area source formed by the fault system (IGN and UPM, 2013).

Synthetic seismicity appears as another tool to address the limitations of seismic record (Robinson et al., 2011; Dieterich, 1994; Ward, 2000; Rundle et al., 2006). Computer models based on earthquake physics allow the creation of synthetic catalogues with the same characteristics as the instrumental ones but extended in time. Time expansion makes it possible to reproduce several seismic cycles (Robinson et al., 2011; Console et al., 2018) and to evaluate different models of seismogenic processes (Wilson et al., 2017). By using rate- and state- dependent friction numerical approach (Dieterich, 1979), it is also possible to simulate the physics of rupture processes of earthquakes (Richards-Dinger and Dieterich, 2012). All of these features make synthetic seismic models useful for gaining a better understanding of earthquake occurrence, as illustrated in this study.

In this study, a classical PSHA methodology (Esteva, 1967; Cornell, 1968) has been applied to 100 different 1 kyr long seismic catalogues extracted randomly from a 1 Myr synthetic catalogue of the EBSZ (Herrero-Barbero et al., 2021; Álvarez-Gómez et al., 2023). Depending on the time interval of analysis, each catalogue is supposed to cover a different part of the seismic cycle of faults. To compare the results obtained, a PSHA is also carried out with data from the historical and instrumental catalogue. The results derived from these analyses will show how the calculation of seismic hazard depends on the time interval in which the seismic catalogue is recorded.

2 Seismotectonic context

Southeastern Spain is characterized by weak to moderate seismicity (Stich et al., 2003; Rodríguez-Escudero et al., 2014). Nevertheless, it is considered the most seismically active region of the Iberian Peninsula (García-Mayordomo et al., 2007), due to the convergence between the African and Eurasian plates (DeMets et al., 1990; Argus et al., 1999). This study focuses on the Eastern Betic Shear Zone (EBSZ), located at the Betic Cordillera, specifically, at its Internal Zones. Betic Cordillera conforms the northern termination of the Rif-Betic Alpine orogenic belt (De-Larouziere et al., 1988). The EBSZ is an active deformation zone that conforms mainly left-lateral strike-slip faults oriented N-S to ENE-WSW (Ortuño et al., 2012). The main faults that constitute the system are (Fig. 1, Table 1), from south to north: the Carboneras Fault (CF), the Palomares Fault (PF), the Alhama de Murcia Fault (AMF), the Los Tollos Fault (LTF), the Carrascoy Fault (CAF) and the Bajo Segura Fault (BSF). The mean slip rate of the system is estimated to be 1.5 ± 0.3 mm yr⁻¹ (Echeverria et al., 2013) unevenly distributed across different faults.

Martínez-Díaz (1998)(Ferrater et al., 2016; Ortuño et al., 2012)(Martínez-Díaz et al., 2012a)Alonso-Henar et al. (2020)(Ferrater et al., 2016)(Herrero-Barbero et al., 2020)Silva et al. (2003)(Moreno, 2011)(Moreno et al., 2015)Echeverria et al. (2015)Booth-Rea et al. (2004)(IGN, 2019)Insua-Arévalo et al. (2015)Martín-Banda et al. (2016)(Martín-Banda, 2020; Sanz de Galdeano et al., 1998; Silva, 1994)Alfaro et al. (2012)Borque et al. (2019)Alfaro et al. (2012)Perea et al. (2012)

Table 1 Model input parameters for fault sections in the EBSZ obtained from Herrero-Barbero et al. (2021).

The EBSZ is characterized by a moderate seismic activity, although it has produced many damaging earthquakes (M_w≥6.0) in historical times: 1048 Orihuela (Silva et al., 2015), 1482 OHerrerorihuela (García-Mayordomo, 2005), 1487 Almería (IGN and UPM, 2013), 1518 Vera (García-Mayordomo et al., 2012a), 1522 Alhama de Almería (Reicherter and Hübscher, 2007), 1673 Orihuela (García-Mayordomo, 2005), 1674 Lorca (Martínez-Díaz et al., 2019), 1804 Dalías (Murphy, 2019), 1829 Torrevieja (García-Mayordomo and Martínez-Díaz, 2006) and 1910 Adra (Rodríguez-Escudero et al., 2014). The most recent damaging earthquake was the Lorca event of 2011, with M_w=5.1, which claimed nine fatalities and huge economic losses (Martínez-Díaz et al., 2012b; López-Comino et al., 2012). All these major earthquakes have been associated with the faults of the area (Yazdi et al., 2023; Yazdi and García-Mayordomo, 2024). For producing a PSHA of the region, the EBSZ has been modeled previously as a seismogenic zone (García-Mayordomo et al., 2012a; IGME, 2015) which gathers the seismic features of the fault system.

Figure 1 Main faults which constitute the EBSZ (García-Mayordomo et al., 2012b, 2017): the Carboneras Fault (CF), the Palomares Fault (PF), the Alhama de Murcia Fault (AMF), the Los Tollos Fault (LTF), the Carrascoy Fault (CAF) and the Bajo Segura Fault (BSF). The representation of GNSS vectors (Echeverria et al., 2013; Borque et al., 2019) is made with fixed Eurasia. The earthquakes are sourced from the catalogue used in this study (IGN, 2022).

3 Materials and methods

In order to gain insight into the influence of the observation time interval on seismic hazard results at the EBSZ, we have resorted to the use of synthetic seismic catalogues. From a 1 Myr-long computed catalogue (Herrero-Barbero et al., 2021; Álvarez-Gómez et al., 2023), we randomly selected ten thousand 1 kyr sub-catalogues to characterise the seismic hazard within a classical PSHA framework. These 1 kyr sub-catalogues were extracted to match the duration of the historical catalogue (spanning from 1048 AD and therefore approximately 1 kyr long). Each sub-catalogue represents a possible realisation of the long-term seismic process, allowing us to compare the hazard estimated from multiple stochastic 1 kyr realisations with that derived from the observed historical and instrumental seismicity, which can likewise be interpreted as a single realisation.

Synthetic seismicity consists of a million-year catalogue generated using RSQSim earthquake simulator. RSQSim simulations require a fault system model that includes a representation of the geometry, kinematics and frictional parameters of the fault segments. The implementation of the EBSZ in RSQSim (Herrero-Barbero et al., 2021) is based on a fully-interacting 3D model built externally to the simulator. This model is composed of triangular cells of 1 km² resolution, allowing the simulation of events with magnitudes greater than M_w=3.0 and approximating some geometric complexities. Rate- and state-dependent friction parameters (Dieterich, 1979) control the sliding resistance and are selected in this model following realistic relations between rupture area and earthquake magnitude. Seismic cycles are divided into several states (healing, nucleation and seismic rupture) that are passed by fault elements throughout the simulation. The evolution of stress, slip variation and state variable are controlled by analytical expressions that use event-driven steps. Mean slip rates of faults are used by the simulator to apply tectonic loading to the faults by the backslip method (Savage, 1983). The result of the simulations applied to this model is an earthquake catalogue characterized by a Magnitude Frequency Distribution (MFD) that tends to fit well with the historical and instrumental features. For more details on the physics-based model and the validation procedure of the best-fit catalogue, see Herrero-Barbero et al. (2021, 2023), Álvarez-Gómez et al. (2023), and Gómez-Novell et al. (2024).

The historical and instrumental earthquake catalogue in SE Spain is compiled by the Instituto Geográfico Nacional (IGN, National Geographic Institute). It can be downloaded freely online (IGN, 2022). The earthquake catalogue provides information on the date, the location, and the magnitude of earthquakes. To conduct this study, it has been necessary to process the catalogue collected until 2022, since the last hazard studies (IGN and UPM, 2013) were carried out in 2012. Originally, earthquakes were reported on heterogeneous magnitude scales. In order to obtain the size of earthquakes on the same magnitude scale, a homogenization process is needed. Using the equation proposed by Cabañas et al. (2015), from a Reduced Major Axis Regression (RMA), the magnitudes m_bLg and m_b are converted to moment magnitude, M_w. Then, earthquakes with magnitudes smaller than M_w=3.0 are removed from the catalogue, since the equations of Cabañas et al. (2015) do not work in this range of values and the completeness of the catalogue cannot be guaranteed. In addition, earthquakes occurring at depths greater than 35 km are removed as they do not impact the seismic hazard (IGN and UPM, 2013).

According to the Poisson hypothesis (Cornell, 1968; Cornell and Winterstein, 1988), earthquakes must be considered as independent phenomena that occur randomly in time. For this purpose, a declustering process (Gardner and Knopoff, 1974) has been applied to the earthquake catalogue using a set of equations developed ad-hoc for the region (IGN and UPM, 2013) (Eqs. 5, 6, and 7).

$$\begin{array}{}\text{(5)}& \log \left(T\right)=\mathrm{0.6}\cdot M_{\mathrm{w}}-\mathrm{0.7}\end{array}$$

for M_w<5.9

$$\begin{array}{}\text{(6)}& \log \left(T\right)=\mathrm{0.04}\cdot M_{\mathrm{w}}+\mathrm{2.6}\end{array}$$

for M_w≥5.9

$$\begin{array}{}\text{(7)}& \log \left(L\right)=\mathrm{0.2}\cdot M_{\mathrm{w}}+\mathrm{0.4}\end{array}$$

where T represents the time interval between earthquakes; M_w, the maximum magnitude of each cluster and L the linear distance that separates two earthquakes. M_w is determined by means of an iterative process (Wiemer, 2001) in which the largest earthquake is selected from a cluster of related events. The declustering process allows differentiation between mainshocks and fore- after-shocks. Mainshocks are selected as independent events with information that is representative of the cluster. This procedure has been applied to both historical-instrumental and synthetic catalogues as is the common procedure on classical PSHA analyses.

The characterisation of the EBSZ as a seismogenic source is represented by the Gutenberg-Richter relationship and the maximum magnitude of each catalogue. Different scenarios for seismic hazard assessment have been generated using the Mersenne-twister random number generator (Matsumoto and Nishimura, 1998) and selecting intervals of 1 kyr from the general synthetic catalogue of 1 Myr. Each sub-catalogue is subjected to a PSHA, which results in an individual hazard curve for the selected sites.

The input data for the hazard calculation consist of the annual earthquake rate, the slope of the Gutenberg-Richter relationship and the maximum magnitude of each catalogue at the EBSZ. The seismic parameters of the seismogenic zones adjacent to the EBSZ zone are also required for the model. These data have been adopted from the ZESIS database (Table 2). These parameters were obtained within the framework of the National Seismic Hazard Map project, following a process of declustering and magnitude homogenization of the earthquake catalogue (IGN and UPM, 2013). In that work, different zoning models (García-Mayordomo et al., 2012b) and procedures for calculating the Gutenberg-Richter parameters were tested and eventually incorporated into the hazard calculations through a logic-tree approach. The maximum magnitude parameter was modelled using a normal probability distribution. The mean maximum magnitude and its standard deviation were derived for each zone based on criteria that balanced the available information on the presence and size of active faults and the magnitude of the largest historical events. Active fault parameters were obtained from the QAFI database (IGME, 2022) which at that time was version 2. Several procedures were evaluated to convert fault-size parameters into moment magnitude using empirical relationships from the literature (IGN and UPM, 2013). Table 2 shows the values referenced for the ZESIS model (IGME, 2015), which we adopt here.

For each scenario, a MFD has been created and adjusted to a Gutenberg-Richter power relationship using the maximum likelihood fitting technique (Eq. 8, Utsu, 1964; Aki, 1965; Sandri and Marzocchi, 2007).

$$\begin{array}{}\text{(8)}& b={\displaystyle \frac{{\log }_{\mathrm{10}}e}{(\stackrel{\mathrm{‾}}{M}-M_{\mathrm{c}})}}\end{array}$$

where $\stackrel{\mathrm{‾}}{M}$ describes the mean magnitude of earthquakes with a magnitude at least M_c, representing the magnitude of completeness of the catalogue.

The magnitude of completeness has been calculated using the Goodness-of-fit test (GMF, Wiemer and Wyss, 2000). Maximum magnitude has been calculated, following the procedure of IGN and UPM (2013), as an average between the geological maximum and the maximum magnitude of each catalogue.

Table 2 Seismic parameters of the seismogenic zones adjacent to the EBSZ zone: β represents the slope of the Gutenberg-Richter relationship; τ(4.0) is the annual rate of earthquakes with magnitude M_w≥4.0; and “MM” is the mean maximum magnitude with its standard deviation (σMM) (after IGN and UPM, 2013; IGME, 2015).

In the development of seismic hazard studies, a set of Ground Motion Prediction Equations (GMPE) is often used in combination with a logic tree. In this case, in which we are interested in the variation of results linked to the observation time window, we use a single GMPE, to avoid including more parameters that might condition the variability of the results. GMPE has been selected as the one derived by Campbell and Bozorgnia (2014) because it includes parameters that describe the behaviour of a fault system, such as the style of faulting and the depth of nucleation of earthquakes. Other studies of seismic hazard on the Iberian Peninsula had also chosen Campbell and Bozorgnia (2014) as GMPE (Rivas-Medina et al., 2018; Gómez-Novell et al., 2020).

Hazard curves reveal the estimated return period for different ground motion levels. Using three types of seismic parameters for the EBSZ: the ones that came from (1) synthetic seismicity, (2) historical and instrumental seismicity updated to the present and (3) ZESIS database, this study conducts a PSHA for six of the most populous cities in SE Spain (Fig. 2, Table 3). Hazard curves have been computed using R-CRISIS software (Ordaz and Salgado-Gálvez, 2020), which calculates seismic hazard using a fully probabilistic approach. Using seismic and geometric features of seismic sources, R-CRISIS performs a spatial integration by subdividing the original source. This procedure is based on the assumption that seismicity is evenly distributed over the area source (Ordaz et al., 2021). Considering the selected GMPE, R-CRISIS evaluates the exceedance probability at different values of ground motion at the selected sites. The probabilities of all the magnitudes and seismic sources are summed to obtain the total exceedance values. Exceedance values have been transformed into return periods. The results obtained consist of 100 curves from 100 synthetic catalogues, one curve for the historical and instrumental catalogue updated to the present, and one curve for ZESIS data. A hundred synthetic scenarios have been considered enough to show whether the variability within the time interval of analysis exists and to quantify it.

Table 3 Populous cities at SE Spain that have been selected to perform a PSHA based on synthetic seismicity.

Figure 2 Seismogenic sources adjacent to EBSZ (IGME, 2015): Central Guadalquivir (CG), Cazorla-Segura & Albacete pre-Betics (CS-ApB), Valencian Platform & Alicante pre-Betics (VP-ApB), Western Internal Betics (WIB), Granada basin (Gb), Sierra Nevada-Filábride & Guadix-Baza basin (SNF-GBb), Murcian pre-Betics (MpB), Central Internal Betics (CIB), Eastern Internal Betics (EIB), Western Alboran sea (WAs), Argelian-Balear basin (ABb), Alborán-Rif Central crest (ARCc), Eastern Rif (ER) and Orán (O) and cities where seismic hazard analyses have been performed: Alicante (ALI), Torrevieja (TOR), Murcia (MUR), Lorca (LOR), Vera (VER) and Almería (ALM).

4 Results

4.1 Magnitude Frequency Distribution

MFD displays the rate of earthquakes over a range of magnitudes. The curve of accumulated events is described by the Gutenberg-Richter relationship (Fig. 3, in blue). In this study, the regression resulting from a maximum likelihood curve fitting technique (Eq. 8) for the historical and instrumental catalogue, has coefficient values of a=3.8 and b=1.0. On PSHA, the annual rate of earthquakes, τ, higher or equal to a given magnitude is a needed parameter. Following IGN and UPM (2013), the selected magnitude for this study is M_w≥4.0. In the case of the EBSZ, we obtained a value of $\mathit{\tau }(M_{\mathrm{w}}\ge \mathrm{4.0})=\mathrm{0.63}$, which means that, on average, less than one earthquake with a magnitude equal to or higher than four is recorded in this region per year.

Figure 3 MFD of the EBSZ obtained from the historical and instrumental record since 1048 to 2022. Discrete rate of annual earthquakes, in pink; cumulative annual rate, in black; Gutenberg-Richter relationship (Eq. 1) using maximum likelihood fitting technique (Eq. 8, Utsu, 1964; Aki, 1965), in blue; and the magnitude of completeness (M_c), in red.

In the case of synthetic sub-catalogues, a variety of frequency–magnitude distributions (Fig. 4) with their corresponding Gutenberg-Richter parameters (Fig. 5) can be observed. The distribution of a (Fig. 5a) shows two maxima between the values a=3.0–3.5 and a=4.5–5.5. In the case of b (Fig. 5b), a similar pattern can be found. The slope of the Gutenberg-Richter relationship presents both maxima at b=0.8–0.9 and b=1.2–1.3. If we evaluate the frequency distribution of a and b using ten thousand sub-catalogues as a Gaussian distribution, we obtained that it is centred at a=4.5 and b=1.1. These mean values are significantly close to the values which are obtained using the historical and instrumental seismic record, especially in the case of b.

Figure 4 Frequency magnitude distribution for the 100 synthetic catalogues (grey lines) and the historical-instrumental catalogue (red line). The dashed line represents the magnitude of completeness for the full length of synthetic catalogues according to Herrero-Barbero et al. (2021). The magnitude of completeness for the historical-instrumental catalogue is M_w 3.5 (IGN and UPM, 2013).

Figure 5 Statistical variability of (a) a – parameter (Gutenberg-Richter relationship), (b) b – parameter (Gutenberg-Richter relationship) and (c) M_max for 10 000 synthetic sub-catalogues and (d), (e) and (f) for 100 synthetic sub-catalogues.

Maximum magnitude values (Fig. 5c) shows a bimodal pattern with its main maximum at M_max=7.3–7.4 and a secondary maximum over M_max=6.3–6.4. Given that we are considering the fault system as a whole, this bimodal pattern could be related to differences in the behaviour of faults. A segmentation of the synthetic sub-catalogues by faults has been carried out towards the evaluation of earthquake nucleation at each fault (Fig. 6). Obtained results reveal that, on the one hand, the Alhama de Murcia fault (AMF, Fig. 6) and the Carboneras fault (CF, Fig. 6) keep the same bimodal distribution as the whole system. On the other hand, the Carrascoy fault (CAF, Fig. 6) and the Bajo Segura fault (BSF, Fig. 6) only present a normal distribution centred in M_max=5.6 for both faults. The nucleation of large earthquakes in the model, with M_max>7.1 in the CAF and the BSF has not been notified in these sub-catalogues. According to previous research in the region (Herrero-Barbero et al., 2021), nucleation of large earthquakes is possible on the BSF, but less likely than in other faults of the system, such as the AMF and the CF.

Figure 6 Distribution of maximum magnitudes for 10 000 sub-catalogues divided by faults: the Alhama de Murcia Fault (AMF), the Carboneras Fault (CF), the Carrascoy Fault (CAF) and the Bajo Segura Fault (BSF). The “count” represents the number of times the 10 000 sub-catalogues show the corresponding magnitude as the maximum magnitude on the indicated fault.

4.2 Seismic Hazard Assessment

Hazard curves have been calculated for Murcia, Lorca, Alicante, Vera, Torrevieja and Almería (Fig. 7) using 100 synthetic subcatalogues. We consider a representative sample of the full set of 10 000 catalogues generated. The parameter distributions obtained from the 100 sub-catalogues reproduce the central tendency and dispersion of the reference set of 10 000 sub-catalogues (Fig. 5d, e, f). No systematic bias is observed in the distributions of the a, b, and maximum magnitude values; therefore, the sample of 100 sub-catalogues is statistically representative of the variability of the dataset. The differences are thus attributable to expected fluctuations in finite samples. This number of catalogues offers a balance between statistical robustness and computational feasibility, ensuring that the variability of seismic behaviour is sufficiently captured while maintaining compatibility with the operational constraints of R-CRISIS. These hazard curves show the variation of return period for the same Peak Ground Acceleration (PGA). The hazard curves obtained using historical-instrumental and ZESIS data show less hazardous scenarios for almost all PGA values. Hazard curves are plotted on a logarithmic scale, which hinders the correct perception of variability. Evaluating specific return periods with seismic design interest (i.e. 100, 500, 1000, 2000, 5000 and 10 000 yr) for the cities of Murcia, Lorca, Alicante, Vera, Torrevieja and Almería (Fig. 8) it can be shown that variability in PGA is not uniform. From one side, for small return periods, the range of values is tiny and variability seems to be modest. On the other side, for return periods longer than 5000 yr, the range of PGA value increases, which translates into higher variability. For large return periods, there is an asymmetry in the dispersion of data with a predominance of higher PGA.

Figure 7 Seismic hazard curves obtained from 100 sub-catalogues of a thousand years build upon a million-year synthetic catalogue (Herrero-Barbero et al., 2021), in grey; upon updated to 2022 catalogue, in red; upon ZESIS database (IGME, 2015), in blue at cities: (a) Murcia, (b) Lorca, (c) Alicante, (d) Vera, (e) Torrevieja y (f) Almería.

Figure 8 PGA variability for different return periods at the cities of (a) Murcia, (b) Lorca, (c) Alicante, (d) Vera, (e) Torrevieja y (f) Almería. The median is plotted in red; percentiles 25 % and 75 %, in blue; minimum and maximum values, in black and anomalous values as a “+” in red.

A comparison of the mean hazard curves between cities (Fig. 9a) reveals that shorter return periods are found for the city of Murcia for all values of PGA. In contrast, Alicante is the location where return periods are longer for all PGA with a significant difference. Meanwhile, the cities of Lorca, Vera, Torrevieja and Almería have quite similar values. Figure 9b represents the return period variability coefficient for each PGA at studied sites. What stands out in this figure is that variability presents a rapid decrease for small values of PGA, till a tipping point located at PGA=0.05 g for Murcia, PGA=0.05 g for Lorca, PGA=0.03 g for Alicante, PGA=0.05 g for Vera, PGA=0.04 g for Torrevieja and PGA=0.02 g for Almería. Henceforth, the variability coefficient presents a continual growth that reaches the following values for PGA=1 g: 58 %, for Torrevieja; 57 %, for Murcia; 54 % for Vera; 53 %, for Lorca, 40 % for Alicante and 25 %, for Almería. It is notorious that, even though Murcia presents the larger values of PGA, the variability is slightly higher in the case of Torrevieja.

Figure 9 Comparison among the cities studied in this work: (a) mean hazard curves and (b) variability coefficient.

5 Discussion

The main objective of this work is to examine the variability between subcatalogues with analogous characteristics derived from a single synthetic catalogue. This approach ensures that the epistemic uncertainties associated with their generation are common to all subcatalogues, making the differences observed in their comparison independent of those uncertainties. Since our study focuses specifically on the variability introduced by the observation time window, we have deliberately decided not to address other sources of epistemic uncertainty. These uncertainties have been previously explored in Herrero-Barbero et al. (2021), from which the synthetic seismic catalogue used here originates.

5.1 Implications for fault seismogenic behavior

A matter of particular attention in the discourse around physical earthquake simulators, such as RSQSim (Richards-Dinger and Dieterich, 2012), is the MFD obtained with them. These simulations depend on the frictional properties of fault systems, which are difficult to estimate. They control the generation of large earthquakes and, with that, the configuration of the MFD.

The model shows a decrease in the seismic rate between M_w=5.0–6.0, suggesting a lack of moderate earthquakes compared to large ones. This scarcity may be related to the lack of secondary faults, such as the Torremendo (Alfaro et al., 2012) or Hinojar (Baena-Escudero, 1993) faults, in the model. It should be noted that this model only represents the main seismogenic sources in the region, which are responsible for the most dangerous earthquakes.

The description of the earthquake catalogue using a Gutenberg-Richter relationship does not contemplate the characteristic earthquake behaviour in the fault system (Youngs and Coppersmith, 1985; Wesnousky, 1994). The choice between Gutenberg-Richter relationship or the characteristic earthquake model depends on the knowledge of faults and data available (Valentini et al., 2017). The characteristic earthquake model suggests that earthquakes at large magnitudes occur with a higher frequency than expected if a potential relationship is considered. According to Ben-Zion et al. (2003), the cyclic repetition of large earthquakes can be described as the dynamic regime of intermittent criticality caused by the spontaneous evolution of stress heterogeneities. The statistical distribution of earthquake magnitudes obtained by Ben-Zion et al. (2003) on a mature fault with positive dynamic weakening is consistent with the characteristic earthquake model. The repeated presence of local maxima at high magnitudes in the synthetic MFDs compatibles with the dimensions of the coherent brittle zones of faults (Ben-Zion and Rice, 1993) may suggest that EBSZ can also be described as an intermittent criticality dynamic system. Although these large events do not necessarily dominate hazard at engineering return periods in classical PSHA frameworks, their intermittent occurrence strongly influences the variability of hazard estimates when different observation time windows are considered. Large earthquakes in this region could be caused by multi-fault ruptures involving several fault segments which are kinematically compatible (Álvarez-Gómez et al., 2018; Herrero-Barbero et al., 2021; Herrero-Barbero, 2022). The results suggest that these multi-fault ruptures are more common than previously thought. The appearance of a bimodal pattern in the maximum magnitude frequency evaluation with a main maximum (Fig. 5c) at large magnitudes appears to support this assumption. The maximum magnitude count appears at M_max=7.2. According to the simulation, these magnitudes would be associated with large earthquakes nucleated mainly on the Alhama de Murcia Fault (AMF, Fig. 6) and on the Carboneras Fault (CF, Fig. 6).

5.2 Seismic Hazard Assessment

The hazard curves calculated in this study (Fig. 7) show different return periods for each catalogue. This variation results from the different values of b, M_max and τ(4.0) calculated for each catalogue, reflecting the variability of seismicity over finite time windows. Each 1 kyr catalogue represents a possible realisation of the long-term seismic process, and therefore produces a distinct hazard curve. Consequently, the observed differences highlight that seismic hazard calculation depends on the time interval in which an earthquake catalogue is recorded. The use of 100 catalogues for the calculation of the hazard curves is justified by the similarity in the mean parameters of the Gutenberg-Richter relation in both datasets (Table 4).

Table 4 Statistical parameters for earthquake catalogues using 100 and 10 000 sub-catalogues.

The notified variability of the hazard curves (Fig. 8) increases for long return periods and PGA values, both associated with large magnitude earthquakes. For a 500-year return period (Table 5), which is often used in seismic building codes for structures of normal importance, the values obtained for PGA show a range of scenarios that implies significant uncertainties that could affect the definition of seismic actions.

Table 5 PGA associated to RP=500 yrs for the cities considered in the PSHA calculations.

Figure 10 Comparison of the return periods obtained in the hazard curves for the synthetic scenarios and the historical-instrumental, HI, one (grey squares) in the cities of: (a) Murcia, (b) Lorca, (c) Alicante, (d) Vera, (e) Torrevieja and (f) Almería. The red line indicates the hypothetical scenario in which the results were identical for both datasets.

Figure 11 Frequency distribution of cumulative seismic moment for 100 sub-catalogues including the historical and instrumental result (in lilac).

The variability of hazard among cities (Fig. 9a) shows the different influence that the EBSZ faults have on the calculation of seismic hazard in each location. The variability observed for the city of Torrevieja is higher, suggesting that there may be an important dependence between the hazard at Torrevieja and the seismic parameters of the Bajo Segura Fault at the NE extreme of the EBSZ. The return period variability coefficient (Fig. 9b) ranges from 11 % to 21 % for a PGA=0.04 g, and from 25 % to 58 % for a PGA=1 g. This range of variability means that the standard deviation is high in relation to the mean hazard curve. This fact illustrates the impact that the use of different seismic catalogues can have on the quantification of seismic hazard at specific locations.

5.2.1 Synthetic vs. historical-instrumental hazard curves

Hazard curves (Fig. 7) obtained using historical and instrumental records, as well as ZESIS database, represent one of the less conservative results for almost all the values of PGA. If return periods are compared (Fig. 10), it can be noticed that the values obtained through synthetic seismicity are, in general, lower than those derived from historical and instrumental data. The lower hazard obtained from historical and instrumental catalogues likely results from a limited sampling of large-magnitude earthquakes occurred within the observation window, making them representative of a less conservative seismic scenario. This would explain the longer return periods compared to those derived from synthetic catalogues. Still, the possibility that the synthetic model overestimates seismic activity rates cannot be entirely ruled out.

The seismic moment released in the different sub-catalogues is calculated (Fig. 11) to assess the similarity between the synthetic scenarios and the real one. The total cumulative seismic moment for the historical and instrumental catalogue is $M_{\mathrm{0}}=\mathrm{1.779}\times {\mathrm{10}}^{\mathrm{19}}$ N m. Synthetic values range between $M_{{\mathrm{0}}_{\mathrm{m}}}=\mathrm{9.8}\times {\mathrm{10}}^{\mathrm{18}}$ N m to $M_{{\mathrm{0}}_M}=\mathrm{4.1}\times {\mathrm{10}}^{\mathrm{20}}$ N m in the case of 100 catalogues. The mean cumulative seismic moment for 100 synthetic catalogues is $\stackrel{\mathrm{‾}}{M_{\mathrm{0}}}=\mathrm{8.09}\times {\mathrm{10}}^{\mathrm{19}}$ N m. The median $M_{{\mathrm{0}}_{\mathrm{med}}}=\mathrm{3.2}\times {\mathrm{10}}^{\mathrm{19}}$ N m is slightly lower than the mean and it is closer to the historical and instrumental value. The histogram of values of the cumulative seismic moment for the synthetic catalogues shows a maximum within the range that contains the cumulative moment of the historical and instrumental catalogue. However, this histogram spreads to values higher than the historical and instrumental one. This finding may suggest that the cumulative seismic moment released as earthquake occurrence could be overestimated. The assumption of total coupling likely inflates the cumulative seismic moment and, consequently, the estimated hazard. Introducing partial coupling or aseismic slip in future models could reduce this bias.

Our findings are consistent with the conclusions of Beauval et al. (2008), highlighting that short observation windows can significantly underestimate seismic hazard. By providing a broader sampling of earthquake variability, our synthetic catalogues underscore the importance of incorporating long-term fault behavior into hazard assessments, especially where the historical record is limited or incomplete.

5.3 Applicability of the model

The methodology applied in this study uses synthetic seismicity as data source. The comprehensive approach of this work, which requires evaluating a large amount of data spread out over time, hinders to replicate the experience with real data. Consequently, the model proposed is subordinated to the constraints associated with synthetic seismicity. The main limitations are related to the absence of secondary faults, the lack of dynamic interactions and the uncertainties of the data input in the model. For a detailed revision see Herrero-Barbero et al. (2021).

The observed variability is likely not an artifact of the model but rather a reflection of the intrinsic complexity of the fault system's seismic cycle. The EBSZ has been proposed as the selected region for this study because of the wide knowledge that we have of this region. Importantly, the methodology is transferable and could be reproduced on other fault systems. It will be an exercise of great interest to develop the same kind of study in a region with faster slip rates and compare the differences with a region of slower ones as proposed in this study.

6 Conclusions

Seismic hazard estimates depend strongly on the characterization and parametrization of the seismogenic sources model used. When the source model characterization relies entirely on the historical earthquake catalogue, the incompleteness of the data becomes a major limitation, as fault seismic cycles are usually longer. On slow-moving fault systems, such as the EBSZ, the historical and instrumental record only represents a small part of the seismic cycle and the occurrence of large earthquakes could have not been recorded yet. In this study, we have shown how the variability of the time interval of analysis leads to different results in seismic hazard assessment. For the selected calculation sites, the return period for different PGA values showed a strong dependency on the time window in which the catalogue has been recorded. The variation coefficient ranges from 11 % to 21 % for a PGA=0.04 g, and from 25 % to 58 % for a PGA=1 g. Using synthetic seismicity, this study evidences that there is a dependence in the quantification of seismic hazard with the time interval in which the earthquake catalogue spans.

The evaluation of a statistically representative set of 10 000 sub-catalogues reveals the recurrent appearance of large-magnitude events (M_w>7.2) associated with multi-fault ruptures. Although infrequent, these events reflect the long-term seismogenic behaviour emerging from the synthetic fault system, typically nucleating on the AMF and CF, which exhibit the highest slip rates within the EBSZ. Rather than implying a systematic dominance of large earthquakes in seismic hazard at engineering return periods, their presence highlights how rare events may significantly influence hazard variability when different observation time windows sample the seismic cycle unevenly.

At places where the geometry and behaviour of faults are well-known, synthetic models of seismicity can be built with confidence. This approach is especially useful in cases where systems have slow rates of deformation and seismicity is just low or moderate. In these cases, synthetic catalogues can reveal the inner variability of such systems and clarify the seismogenic potential of faults. Nevertheless, they are constrained to the quality of the geological knowledge available for the active faults present in the region. To better understand the models and their impact on seismic hazard, further investigations should focus on issues such as fault segmentation, structural development, frictional properties, and degree of coupling (aseismic deformation and secondary faulting).