Approach for combining fault and area sources in seismic hazard assessment: application in south-eastern Spain

This paper presents a methodological approach to seismic hazard assessment based on a hybrid source model composed of faults as independent entities and zones containing residual seismicity. The seismic potential of both types of sources is derived from different data: for the zones, the recurrence model is estimated from the seismic catalogue. For fault sources, it is inferred from slip rates derived from palaeoseismicity and GNSS (Global Navigation Satellite System) measurements.

Distributing the seismic potential associated with each source is a key question when considering hybrid zone and fault models, and this is normally resolved using one of two possible alternatives: (1) considering a characteristic earthquake model for the fault and assigning the remaining magnitudes to the zone, or (2) establishing a cut-off magnitude, M_c, above which the seisms are assigned to the fault and below which they are considered to have occurred in the zone. This paper presents an approach to distributing seismic potential between zones and faults without restricting the magnitudes for each type of source, precluding the need to establish cut-off M_c values beforehand. This is the essential difference between our approach and other approaches that have been applied previously.

The proposed approach is applied in southern Spain, a region of low-to-moderate seismicity where faults move slowly. The results obtained are contrasted with the results of a seismic hazard method based exclusively on the zone model. Using the hybrid approach, acceleration values show a concentration of expected accelerations around fault traces, which is not appreciated in the classic approach using only zones.

1 Introduction

Active faults are the main earthquake sources in the crust. However, their incorporation in seismic hazard assessment is not straightforward since there are not enough data available to adequately model them. This leads to a limited use of faults as independent sources in seismic hazard analyses and to an extended use of seismic zones that cover a significant portion of the crust, assuming uniform seismic characteristics within each source.

This situation has begun to change in recent years, as more studies on active tectonics, palaeoseismicity and fault deformation rates derived from GNSS and other measurements become available. These recently available studies constrain fault parameters such as rupture plane geometry, predominant sense of slip, slip rates, etc. (e.g. Dixon et al., 2003; Langbein and Bock, 2004; Papanikolaou et al., 2005; Walpersdorf et al., 2014; Metzger et al., 2011).

Taking fault type rather than zones into consideration in seismic hazard studies requires addressing two factors: the 3-D geometry of the source and the data required to characterize its seismic potential. In most practical cases, the seismic potential of faults is characterized based on the slip rate using characteristic earthquake models proposed by Wesnousky (1986) (for instance, Field et al., 2014; Akinci and Pace, 2017) instead of Gutenberg–Richter recurrence models (Parsons and Geist, 2009). Other approaches such as extracting the seismic parameters of every single fault from the earthquake catalogue are not always viable, especially in areas with slow-moving faults. Additionally, the period considered in the catalogue may be too short compared with the recurrence time of the fault to provide an unbiased estimation of fault seismic parameters.

In principle, modelling all existing active faults as independent entities could be conceived as the most accurate source model for seismic hazard assessment. However, this vision is still rather idealistic. A more realistic view would include only a limited number of active faults (those with the highest seismic activity) as independent sources. Accordingly, small faults that generate low-magnitude events or slow faults that produce rare events cannot be properly characterized. To prevent a possible deficit in the seismic source model for a given region, the use of faults as seismic sources may be completed with zones that account for the seismic potential associated with these small or slow faults or simply with unknown faults that cannot be characterized independently. Hence, we propose considering a hybrid source model composed of faults and zones: the first modelled as independent sources and the second including residual seismicity.

Adequately establishing the distribution of seismic potential using a model that combines zones and faults poses a challenge, since these are derived from different data sources. For zones, the recurrence model is calculated based on the seismic catalogue, whereas for faults, the recurrence model is derived from fault geometries and slip rate estimates based on GNSS-measured deformation rates. The problem is that some of the events contained in the catalogue may be associated with the faults and may have already been included when calculating the seismic potential of the faults based on the slip rate estimates. If all events are assigned to the zone, the events associated with the faults would be counted twice, leading to an overestimation of the total seismic potential (for both faults and zones).

Some authors assign initial β values to seismic sources (e.g. Bungum, 2007) or propose a simple way of distributing the seismic potential based on a uniform magnitude value, M_c, assigning events with a magnitude lower than M_c to the zone and events with magnitude higher than M_c to the faults (Frankel, 1995; Woessner et al., 2015). The question is, how is this M_c value determined? Why can the fault not generate events with a magnitude below the M_c value? In a study region such as southern Spain, with slow faults and maximum magnitudes around 6.5–7.0, it is difficult to choose an M_c in a non-arbitrary manner.

The approach presented in this paper addresses the challenging question of how to estimate the anticipated ground motion exceedance rate using a short period of earthquake observations and limited geological data (with significant uncertainties). This challenge is common to all probabilistic seismic hazard models (Kijko et al., 2016). The purpose of this study is to approach this challenge proposing a model that contains different types of seismic sources (faults and zones) and adequately distributes the seismic potential, preventing double counting and taking completeness periods into account.

Figure 1 Diagram showing the distribution of the seismic potential of a region, expressed as the sum of the seismic potential of the faults and the seismic potential of the zone.

An application of the approach presented is carried out in SE Spain, the area with the highest seismic hazard in Spain. Most of the previous work that partly or wholly addresses this area includes zones only (García-Mayordomo et al., 2007; Benito et al., 2010; Mezcua et al., 2011; IGN-UPM working group, 2013; Salgado-Gálvez et al., 2015) or is based on zoneless methods (Peláez and López-Casado, 2002; Crespo et al., 2014). A first attempt to combine faults and zones was carried out by García-Mayordomo (2005), who developed a zone model for the area taking into account the use of the characteristic earthquake model for faults.

2 Source (zones and faults) hybrid approach to hazard estimation

The hybrid model proposed is composed of fault-type sources and zone-type sources. In addition, the term “region” is defined as the geometric container for both source types. Thus, the region presents the same geometry as the zone and its seismic potential (seismicity rate and seismic-moment rate) is the sum of the potentials of the two types of sources (faults and zone). The zone is used to represent the seismic potential of events that cannot be associated with specific faults. Although there is a geometrical equivalence between region and zone, their seismic potential is very different, as the seismic potential of the region equals the seismic potential of the zone plus the seismic potential of the faults contained within the region (Fig. 1).

Figure 2 Completeness analyses of the seismic catalogue. Panel (a) shows the (normalized) cumulative number of earthquakes per year for different magnitude intervals. Solid circles indicate the inflection point that marks the lower limit of the completeness period for the respective magnitude interval CP(m). Panel (b) shows the CP(m) corresponding to each magnitude interval. Note that the CP(m) is not well constrained for magnitudes above the M_MaxC value (note dashed and dotted curves in both figures).

The problem is then how to distribute the seismic potential of the region between the zone and the faults without counting some faults twice. The following considerations were taken into account:

• The seismicity rate of the region is derived from the seismic catalogue after excluding the events that lie outside their respective completeness periods, CP(m). This period is defined for a given magnitude M as the period during which a catalogue of events of magnitudes M and higher is complete. This is comparable to assuming that all events of a given magnitude, M, that have actually occurred are effectively contained in the catalogue within the period CP(m=M) (and not outside of this period).

• The completeness periods, CP(m), for different magnitudes up to a maximum-completeness magnitude value, M_MaxC, are lower than the observation period (OP) of the catalogue.

• Magnitude values above M_MaxC present recurrence periods higher than the catalogue OP. These values usually constitute a sample that does not include a high enough number of records to clearly establish the recurrence period, as this makes it increasingly difficult to constrain rates for rarer events.

By representing the number of recorded events for different magnitude intervals as a function of time it is possible to identify the reference years, RY(m), for different magnitude intervals using the slope method (Hakimhashemi and Grünthal, 2012), also known as the temporal course of earthquake frequency (TCEF) (Nasir et al., 2013). This method consists of plotting the cumulative number of earthquakes of a given magnitude range over time and estimating the year, presenting a significant gradual change in slope (Fig. 2). Consequently, CP(m) and M_MaxC values can be calculated for $m<=M_{\mathrm{MaxC}}$. It is then possible to estimate seismicity rates (Eq. 1) and the seismic moment in each magnitude interval (Eq. 2) as follows:

$$\begin{array}{}\text{(1)}& {\displaystyle }{\displaystyle }\dot{n}\left(m\right)={\displaystyle \frac{n\left(m\right)}{\mathrm{CP}\left(m\right)}},\end{array}$$

where $\dot{n}\left(m\right)$ is the annual rate of events with magnitude (m) and n(m) is the number of recorded events with magnitude (m) in the catalogue in the completeness period, CP(m).

$$\begin{array}{}\text{(2)}& {\displaystyle }{\displaystyle }{\dot{M}}_{\mathrm{o}}\left(m\right)=\left(m\right)\cdot M_{\mathrm{o}}\left(m\right),\end{array}$$

where M_o(m) is the seismic moment released by events of magnitude m, obtained using the equation proposed by Hanks and Kanamori (1979) ($\log \left(M_{\mathrm{o}}\right)=\mathrm{1.5}\cdot M_{\mathrm{w}}+\mathrm{16.1}$).

Finally, the cumulative rates in the interval [M_Min, M_MaxC] can be estimated, where M_Min is the minimum-magnitude value used to compute seismic hazard, as shown in Sect. 2.1. This is illustrated with an example in Fig. 2, with [M_Min, M_MaxC]=[4.0, 5.9].

Although faults are capable of generating earthquakes with magnitude m>M_MaxC, the distribution of seismic potential is carried out in the completeness period [M_Min, M_MaxC]. In this way, we avoid using magnitudes with long recurrence periods that have not been recorded in the catalogue within the completeness periods. The computation of the seismic potential of the fault in the interval [M_MaxC, M_Max(fault)], where M_Max(fault) is the maximum-magnitude value of events that can be generated in a fault, is constrained with other geological criteria (see below).

The seismic potential is represented by the total rate of earthquakes (${\dot{N}}_{\min }$) and the cumulative rate of seismic moment (${\dot{M}}_{\mathrm{o}}$), for the magnitude range [M_Min, M_MaxC], in the completeness period CP(m). Details on how to determine ${\dot{N}}_{\min }$ and ${\dot{M}}_{\mathrm{o}}$ for the entire region, the corresponding zone and faults are explained in the following section.

2.1 Seismic potential of the region

The ${\dot{N}}_{\min }$ and ${\dot{M}}_{\mathrm{o}}$ values representing the seismic potential of the region are derived from the seismic catalogue of the magnitude interval [M_Min, M_MaxC] for the completeness periods, CP(m).

$$\begin{array}{}\text{(3)}& {\displaystyle }& {\displaystyle }{\left.{\dot{N}}_{{\min }_{\mathrm{region}}}\right|}_{M_{\mathrm{Min}}}^{M_{\mathrm{MaxC}}}=\sum _{M_{\mathrm{Min}}}^{M_{\mathrm{MaxC}}}\dot{n}\left(m\right)\text{(4)}& {\displaystyle }& {\displaystyle }{\left.{\dot{M}}_{{\mathrm{o}}_{\mathrm{region}}}\right|}_{M_{\mathrm{Min}}}^{M_{\mathrm{MaxC}}}=\sum _{M_{\mathrm{Min}}}^{M_{\mathrm{MaxC}}}\dot{n}\left(m\right)\cdot M_{\mathrm{o}}\left(m\right),\end{array}$$

with $\dot{n}\left(m\right)$ the annual rate of events with magnitude (m) recorded in CP(m) and M_o(m) the seismic moment for magnitude m. The notation $X{|}_{M_i}^{M_j}$ represents the magnitude interval (M_i, M_j) in which variable X is computed.

2.2 Seismic potential of faults

The cumulative-moment rate of the faults is estimated assuming that the fault planes are accumulating energy evenly and using the equation proposed by Brune (1968):

$$\begin{array}{}\text{(5)}& {\displaystyle }{\displaystyle }{\dot{M}}_{{\mathrm{o}}_{\mathrm{fault}}}=\mathit{\upsilon }\cdot \mathit{\mu }\cdot A,\end{array}$$

where υ is the slip rate, μ is the shear modulus and A is the area of the fault plane.

The slip rate υ and the area of each fault plane can be derived from specific studies based on paleoseismic analyses and GNSS measurements. There are also some databases available to search for these data, including the EDSF for Europe (Basili et al., 2013), the DISS for Italy (DISS Working Group, 2018) and the QAFI for Spain and Portugal (Garcia-Mayordomo et al., 2012). The shear modulus may be estimated from values close to $\mathit{\mu }=\mathrm{3.2}\times {\mathrm{10}}^{\mathrm{10}}$ Pa (Walters et al., 2009; Martínez-Díaz et al., 2012).

This moment rate represents the average annual seismic moment accumulated in each fault that will be released by earthquakes of different magnitudes m=0 up to the maximum magnitude of the fault, M_Max(fault). The value M_Max(fault) can be evaluated from a geometrical aspect of the fault planes using empirical relationships proposed in the literature, such as Wells and Coppersmith (1994), Stirling et al. (2002) and Leonard  (2010). Thus, ${\dot{M}}_{{\mathrm{o}}_{\mathrm{fault}}}$ can be expressed as follows (Eq. 6):

$$\begin{array}{}\text{(6)}& {\displaystyle }{\dot{M}}_{{\mathrm{o}}_{\mathrm{fault}}}=\underset{Mm=\mathrm{0}}{\overset{M_{\mathrm{Max}\left(\mathrm{fault}\right)}}{\int }}\dot{n}\left(m\right)\cdot M_{\mathrm{o}}\left(m\right)dm,\end{array}$$

where $\dot{n}\left(m\right)$ can be estimated applying a recurrence model, for instance, the modified GR model shown in Eq. (7).

$$\begin{array}{}\text{(7)}& {\displaystyle }{\displaystyle }\dot{n}\left(m\right)={\dot{N}}_{{\min }_{\mathrm{fault}}}\cdot \mathit{\beta }\left({\displaystyle \frac{e^{-\mathit{\beta }m}}{e^{-\mathit{\beta }}\left(M_{m=\mathrm{0}}\right)-e^{-\mathit{\beta }\left(M_{{\mathrm{Max}}_{\left(\mathrm{fault}\right)}}\right)}}}\right),\end{array}$$

and the seismic moment M_o(m), can be estimated from the Hanks and Kanamori (1979) relation expressed in exponential terms, $M_{\mathrm{o}}\left(m\right)=e^{dm+c}$ with $c=\mathrm{16.1}\cdot \ln \left(\mathrm{10}\right)$ and $d=\mathrm{1.5}\cdot \ln \left(\mathrm{10}\right)$, where m is the moment magnitude M_w (Anderson and Luco, 1983).

Substituting the previous relations in Eq. (6), solving the integral and reordering the equation for ${\dot{N}}_{{\min }_{\mathrm{fault}}}$, we get Eq. (9).

$$\begin{array}{}\text{(8)}& {\displaystyle }{\displaystyle }{\dot{N}}_{{\min }_{\mathrm{fault}}}={\displaystyle \frac{{\dot{M}}_{{\mathrm{o}}_{\mathrm{fault}}}\cdot \left(d\mathit{\beta }\right)\cdot \left(e^{-\mathit{\beta }\left(M_{m=\mathrm{0}}\right)}-{\mathrm{e}}^{-\mathit{\beta }\left(M_{{\mathrm{Max}}_{\left(\mathrm{fault}\right)}}\right)}\right)}{\mathit{\beta }\cdot \left[e^{-\mathit{\beta }{M}_{\mathrm{Max}}{M}_{\mathrm{o}}\left(M_{\mathrm{Max}\left(\mathrm{fault}\right)}\right)}-e^{-\mathit{\beta }{M}_{m=\mathrm{0}}{M}_{\mathrm{o}}\left(M_{m=\mathrm{0}}\right)}\right]}},\end{array}$$

where M_m=0 is the minimum magnitude that may be generated at a fault rupture (here taken as m=0), M_Max(fault) the maximum magnitude, and ${\dot{M}}_{{\mathrm{o}}_{\mathrm{fault}}}$ the seismic-moment rate accumulated in the fault (Eq. 5).

The total seismic-moment rate for each fault (${\dot{M}}_{{\mathrm{o}}_{\mathrm{fault}}})$ and seismicity rate (${\dot{N}}_{{\min }_{\mathrm{fault}}}$) can be formulated as the sum of the seismic-moment rate released at different magnitude intervals; thus, it follows

$$\begin{array}{}\text{(9)}& {\displaystyle }{\displaystyle }{\dot{N}}_{{\min }_{\mathrm{fault}}}={\left.{\dot{N}}_{\min }\right|}_{M_{M=\mathrm{0}}}^{M_{\mathrm{Min}}}+{\left.{\dot{N}}_{\min }\right|}_{M_{\mathrm{Min}}}^{M_{\mathrm{MaxC}}}+{\left.{\dot{N}}_{\min }\right|}_{M_{\mathrm{MaxC}}}^{M_{\mathrm{Max}\left(\mathrm{fault}\right)}}\end{array}$$

$$\begin{array}{ll}{\displaystyle }{\dot{M}}_{{\mathrm{o}}_{\mathrm{fault}}}& {\displaystyle }=\sum _{M_{M=\mathrm{0}}}^{M_{\mathrm{Min}}}\dot{n}\left(m\right)\cdot M_{\mathrm{o}}\left(m\right)+\sum _{M_{\mathrm{Min}}}^{M_{\mathrm{MaxC}}}\dot{n}\left(m\right)\cdot M_{\mathrm{o}}\left(m\right)\\ \text{(10)}& {\displaystyle }& {\displaystyle }+\sum _{M_{\mathrm{MaxC}}}^{M_{\mathrm{Max}\left(\mathrm{fault}\right)}}\dot{n}\left(m\right)\cdot M_{\mathrm{o}}\left(m\right).\end{array}$$

By implementing a recurrence model, it is possible to derive the seismicity rate and the moment rate in the interval [M_Min, M_MaxC] (see example in Fig. 3 with [M_Min, M_MaxC]=[4.0, 6.9]).

In this approach it is considered that all faults included in the same region will present the same β value and different seismicity rates (${\dot{N}}_{{\min }_{\mathrm{fault}}}$), as this parameter depends on the seismic-moment rate of each fault (${\dot{M}}_{{\mathrm{o}}_{\mathrm{fault}}}$).

Figure 3 (a) Seismicity rate (cumulative number of events per year vs. magnitude) and (b) moment rate (cumulative seismic moment per year vs. magnitude) plots. The different magnitude intervals mentioned in the text are marked.

2.3 Seismic potential of the zone

The parameters representing the zone are initially unknown. They can be calculated for the interval [M_Min, M_MaxC] given that

$$\begin{array}{ll}{\displaystyle }\mathrm{Seismic}{\mathrm{Potential}}_{\mathrm{zone}}& {\displaystyle }=\mathrm{Seismic}{\mathrm{Potential}}_{\mathrm{region}}\\ \text{(11)}& {\displaystyle }& {\displaystyle }-\mathrm{Seismic}{\mathrm{Potential}}_{\mathrm{faults}},\end{array}$$

or specifically,

$$\begin{array}{}\text{(12)}& {\displaystyle }{\displaystyle }{\left.{\dot{N}}_{{\min }_{\mathrm{zone}}}\right|}_{M_{\mathrm{Min}}}^{M_{\mathrm{MaxC}}}={\left.{\dot{N}}_{{\min }_{\mathrm{region}}}\right|}_{M_{\mathrm{Min}}}^{M_{\mathrm{MaxC}}}-\sum {\left.{\dot{N}}_{{\min }_{\mathrm{fault}}}\right|}_{M_{\mathrm{Min}}}^{M_{\mathrm{MaxC}}}\end{array}$$

$$\begin{array}{}\text{(13)}& {\displaystyle }{\displaystyle }{\left.{\dot{M}}_{{\mathrm{o}}_{\mathrm{zone}}}\right|}_{M_{\mathrm{Min}}}^{M_{\mathrm{MaxC}}}={\left.{\dot{M}}_{{\mathrm{o}}_{\mathrm{region}}}\right|}_{M_{\mathrm{Min}}}^{M_{\mathrm{MaxC}}}-\sum {\left.{\dot{M}}_{{\mathrm{o}}_{\mathrm{fault}}}\right|}_{M_{\mathrm{Min}}}^{M_{\mathrm{MaxC}}}.\end{array}$$

In principle, there are two equations with two unknowns related to the zone:

$${\displaystyle }{\displaystyle }{\left.{\dot{N}}_{{\min }_{\mathrm{zone}}}\right|}_{M_{\mathrm{Min}}}^{M_{\mathrm{MaxC}}},$$

and

$${\displaystyle }{\displaystyle }{\left.{\dot{M}}_{{\mathrm{o}}_{\mathrm{zone}}}\right|}_{M_{\mathrm{Min}}}^{M_{\mathrm{MaxC}}}.$$

Regarding the faults, ${\dot{N}}_{{\min }_{\mathrm{fault}}}{|}_{M_{\mathrm{Min}}}^{M_{\mathrm{MaxC}}}$ and ${\dot{M}}_{{\mathrm{o}}_{\mathrm{fault}}}{|}_{M_{\mathrm{Min}}}^{M_{\mathrm{MaxC}}}$ are derived using an initial (not definitive) β value. Regarding the region, ${\dot{N}}_{{\min }_{\mathrm{fregion}}}{|}_{M_{\mathrm{Min}}}^{M_{\mathrm{MaxC}}}$ and ${\dot{M}}_{{\mathrm{o}}_{\mathrm{region}}}{|}_{M_{\mathrm{Min}}}^{M_{\mathrm{MaxC}}}$ are known, as they were extracted from the catalogue (Eqs. 1 and 2). A new additional equation is obtained relating ${\dot{N}}_{{\min }_{\mathrm{zone}}}{|}_{M_{\mathrm{Min}}}^{M_{\mathrm{MaxC}}}$ and ${\dot{M}}_{{\mathrm{o}}_{\mathrm{zone}}}{|}_{M_{\mathrm{Min}}}^{M_{\mathrm{MaxC}}}$ using Eq. (8) for the interval [M_Min, M_MaxC] in the zone, resulting in the following:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\left.{\dot{N}}_{{\min }_{\mathrm{zone}}}\right|}_{M_{\mathrm{Min}}}^{M_{\mathrm{MaxC}}}=\\ \text{(14)}& {\displaystyle }& {\displaystyle }{\displaystyle \frac{{\left.{\dot{M}}_{{\mathrm{o}}_{\mathrm{zone}}}\right|}_{M_{\mathrm{Min}}}^{M_{\mathrm{MaxC}}}\cdot \left(d-{\mathit{\beta }}_{\mathrm{zone}}\right)\cdot \left(e^{-{\mathit{\beta }}_{\mathrm{zone}}\left(M_{\mathrm{Min}}\right)}-e^{-{\mathit{\beta }}_{\mathrm{zone}}\left(M_{\mathrm{MaxC}}\right)}\right)}{{\mathit{\beta }}_{\mathrm{zone}}\cdot \left[e^{-{\mathit{\beta }}_{\mathrm{zone}}{M}_{\mathrm{MaxC}}}{M}_{\mathrm{o}}\left(M_{\mathrm{MaxC}}\right)-e^{-{\mathit{\beta }}_{\mathrm{zone}}{M}_{\mathrm{Min}}}{M}_{\mathrm{o}}\left(M_{\mathrm{Min}}\right)\right]}}..\end{array}$$

Notice that Eqs. (8) and (14) are similar: they differ in the type of source and computation interval. Equation (8) is for faults and it is computed in the magnitude interval [M_m=o, M_Max]. Equation (14) is for zones, and the magnitude interval is restricted to [M_Min, M_MaxC]. Also note that the β value of the zone in this equation can be equal to the β value of the region, as both sources present similar seismic natures.

Table 1 COV coefficient associated with seismic-moment rate obtained using synthetic catalogues.

Figure 4 Graph extrapolating the recurrence model of the fault up to the maximum expected magnitude value, as deduced from geological criteria.

With this third equation, it is possible to solve the system and obtain a new β value for the faults (second iteration) that balances the three equations. The result is the distribution of seismic potential between the zone and the faults in the interval [M_Min, M_MaxC].

Considering that the faults may generate events with magnitudes larger than M_MaxC, the corresponding distribution of seismic potential in the interval (M_MaxC , M_Max(fault)] is calculated by extrapolating the recurrence model with the last adjusted β value (Fig. 4).

Regarding the M_Max value expected for the zone (M_Max(zone)), this can be considered equal to M_MaxC or extended to a higher-magnitude value if it is assumed that bigger events can occur in other unidentified sources (such as blind faults).

2.4 Analysis of uncertainty

The proposed approach strongly relies on computing seismicity, earthquake rates and moment rates, within the magnitude interval [M_Min, M_MaxC] of the seismic catalogue that contains the complete record of events that have occurred in the entire region.

In order to capture the variability of seismic-moment rates calculated from the earthquake catalogue, a sensitivity analysis of three key factors is conducted. These factors are (1) the number of records used to compute moment rates, (2) the magnitude range covered by the complete catalogue and (3) the proportion of earthquakes of different magnitude (b value).

Figure 5 Three-dimensional view of the seismic sources considered for hazard calculation, including faults (red) and zones (brown).

Synthetic catalogues derived from GR-modified recurrence models are generated for this purpose. Earthquake rates are computed using different numbers of events, magnitude intervals and b values that could be representative of areas of low-to-moderate seismic activity.

The procedure consists of five steps:

• generating 2000 synthetic catalogues for different combinations of earthquake rates, magnitude intervals and b values;

• calculating earthquake rates for different magnitude values for each synthetic catalogue (Eq. 6);

• calculating moment rates for different magnitude values for each synthetic catalogue;

• calculating the sum of moment rates for different magnitude values in order to obtain the cumulative-moment rate for each synthetic catalogue (Eq. 7);

• computing the mean and the standard deviation of the distribution of calculated seismic-moment rates.

Table 1 shows the coefficient of variation (COV; the standard deviation/mean) associated with each combination: number of events, magnitude interval and b value. As can be seen, the greater the number of records in the sample and the lower the magnitude range, the lower the uncertainty associated with the rate of seismic moment calculated. The b value presents a different trend, recording the greatest variability for b values between 1.0 and 1.5. This table is useful for estimating the uncertainty of the seismic-moment rate calculated from the seismic catalogue as a function of the number of earthquakes, magnitude interval and b value.

It is also important to consider the uncertainty associated with the slip rate and the area of the fault, as these are propagated into the distribution of seismic-moment rates of the fault in proportion to the deviation of the area or slip rate value. The uncertainty of the slip rate value is more relevant for low slip rate values than for large slip rate values (a similar trend can be deduced for low and high area values). For instance, a deviation of ±1 mm yr⁻¹ in a slip rate of ±2 mm yr⁻¹ represents an uncertainty of 50 %, leading to a COV value of 0.5 at the moment rate of the fault. However, the same deviation (±1 mm yr⁻¹) for a fault with a slip rate of ±10 mm yr⁻¹ represents an uncertainty of 10 %, leading to a COV coefficient moment rate of only 0.1 for the fault.

3 Application of the approach in south-eastern Spain

The approach described above is applied in south-eastern Spain, the most seismically active area in the country. The tectonic deformation and seismicity is related to the north-western boundary between the Eurasian and African plates (e.g. Kiratzi and Papazachos, 1995), with an approximate shortening rate of about 4 mm yr⁻¹ (Argus et al., 1989) in a roughly NNW–SSE direction. Crustal deformation is accumulated over a broad area in which seismicity is diffuse (Benito and Gaspar-Escribano, 2007).

Assigning earthquakes to specific faults is not an easy task, partly due to errors in earthquake location and to the existence of blind, unknown faults: whereas earthquakes can be clearly associated with a rupture, such as the 2011 M 5.2 Lorca event generated in the Alhama de Murcia fault system (Cabañas et al., 2011), other events have occurred in areas where there are no mapped active faults, for instance the 2007 M_w=4.7 Pedro Muñoz and 2015 M_w=4.7 Ossa de Montiel earthquakes, both located in central Spain (QAFI database, García-Mayordomo et al., 2012).

3.1 Source input data

The seismogenic source model considered for SE Spain is composed of 12 regions that contain a total of 95 faults (Supplement) Active fault data are taken from the QAFI database (v2.0) (García-Mayordomo et al., 2012), which includes information about fault segmentation, geometry and slip rate (see Fig. 5). The maximum expected magnitude in each fault is derived from the rupture length using Stirling et al. (2002) equations derived from the instrumental data set. These equations are chosen because they are also the ones used in the QAFI database and ensure consistency with said database. Moment rates accumulated in the faults are estimated using the fault plane area and the slip rate value according to the formula proposed by Brune (1968). A value of $\mathit{\mu }=\mathrm{3.2}\times {\mathrm{10}}^{\mathrm{10}}$ Pa is assumed for the shear modulus (Walters et al., 2009; Martínez-Díaz et al., 2012).

Table 2 Seismic rate and seismic-moment rate recorded in the different regions for two magnitude intervals (M_Min–M_Max) and (M_Min–M_MaxC) obtained from the seismic catalogue. The table includes the ratio of seismic-moment rate of the two intervals, indicating what percentage of the seismic movement rate liberated from M_min to M_max is contemplated in the magnitude intervals over which hazard is distributed (M_Min–M_MaxC). Note that no faults have been catalogued within regions 28, 29, 33 and 40, which is why no values have been assigned (–).

The zone model proposed by García-Mayordomo et al. (2010) is used to obtain the geometries of the 12 regions (and thus of the zones) that account for the seismicity that cannot be ascribed to faults (see Fig. 5). All the regions considered in this model contain fault sources, with the exception of regions 28, 29, 33 and 40. In these cases, the seismic potential of the corresponding region is assigned to the zones.

The seismic moment released in the region is estimated from the seismic catalogue of Spain homogenized to M_w (IGN-UPM Working Group, 2013; Cabañas et al., 2015). This catalogue contains 1,496 earthquakes, with magnitudes ranging from 4.0 to 6.6. The uncertainty assessment of the catalogue used in this study is explained in Gaspar-Escribano et al. (2015). According to the completeness analysis, a M_MaxC of 5.9 is estimated for SE Spain (although not every region reaches this maximum-magnitude value). The recurrence periods for magnitudes higher than 6 are too long to allow us to establish completeness periods for these magnitude ranges (see Fig. 2).

Table 2 shows the seismic potential for each region, calculated in the magnitude intervals [M_Min, M_MaxC] and [M_Min, M_Max(region)]. It is observed that the seismic potential in the first interval up to M_MaxC, constitutes at least a 60 % of the seismic potential in the second interval, up to M_Max(region).

Subsequently, a recurrence model (GR-mod) is assigned to all regions, obtaining the corresponding b values and COV coefficients (see Table 3). Note that zone 30 lacks a COV estimate because the sample of records (only 7) is very limited, and the [M_Min, M_MaxC] interval is very narrow, resulting in increased uncertainty in the hazard estimates for this region. A GR-mod recurrence model is also assigned to the faults. Finally, the distribution of seismic moments among all seismic sources is carried out (Table 4). As may be observed, the seismic-moment rate associated with the zone has a strong influence on the estimated seismic hazard of the region. This is due to the limited number of known active faults that can be modelled as independent sources, a common situation in areas with low and moderate seismic activity. However, it is worth noting that the seismic potential of regions 35, 36 and 38 is dominated by the seismic potential of faults.

Table 3 Parameters extracted from the seismic catalogue for each region used to estimate the COV coefficient for Table 1, regions 28, 29, 33 and 40 have been excluded because they contain no faults.

The seismic hazard calculation is carried out using the software CRISIS2012 (Ordaz et al., 2013), considering the strong motion equation of Campbell and Bozorgnia (2014), which makes it possible to include the fault geometry and the faulting style. The ground motion parameters predicted include peak ground acceleration (PGA) and 15 spectral accelerations within the period range (0.05–10 s), all obtained in hard soil ($V_{{\mathrm{s}}_{\mathrm{30}}}=\mathrm{760}$ m s⁻¹) conditions.

Table 4 Seismic potential distribution of faults and zones. The last column includes the percentage of regional seismic potential assigned to each source within the region.

3.2 Results

Seismic hazard results obtained with the proposed hybrid model (HM) and with the classical method based in zone (CM) are shown in Fig. 6a and b. Only the geometry of the zone model differs in the two analyses: the ground motion prediction equation (GMPE) and the other calculation parameters are the same in both approaches. The definition of seismic zones applied in the classic method is explained with detail in IGN-UPM Working Group (2013).

PGA estimates for the return period of 475 years using the zone approach (CM) reach maximum values in Granada, Almeria and the Murcia region, around 0.20 g. Minimum PGA values are obtained in Jaén, with values as low as 0.06 g.

Figure 6a shows the seismic hazard map resulting from applying our approach (HM). It can be seen that the largest accelerations are estimated around the Carboneras Fault and the fault set of Granada, (0.38 g), followed by the Alhama de Murcia and La Viña faults systems (0.30 g) and, to a lesser extent, by the Venta de Zafarraya, Carrascoy, Bajo Segura, Baza, Mijas and Cartama fault systems.

The seismic hazard map obtained using the HM displays more spatial variability than the one obtained with the CM, showing maximum values along fault sources that decrease sharply away from the faults. This trend reflects a source proximity effect, implying higher acceleration values for the surface projection of the fault rupture plane that rapidly decrease away from the fault (by one half at a distance of about 15 km).

The differences between the expected maximum acceleration obtained with the two methods, CM and HM, for return periods of 475 and 4975 years appear in Fig. 7a and b. The trend presented in both maps is very similar for the two return periods. A different case is found in region 30 (Case Lietor Fault), a very complex region with scarce seismic activity and large faults with low slip rates (see Supplement). Here, the HM gives higher seismic hazard than the CM only for long return periods. For this region, the magnitude range [M_Min, M_MaxC] is very small and it is necessary to extrapolate the model to a larger scale, given the high uncertainty shown in Table 3. However, the results reflect that, for longer periods, these slow faults play a relevant role in the seismic hazard of the region (see Fig. 8), where the HM hazard curve reflects a substantial increase in hazard for long return periods.

To clarify how faults are conditioning the final seismic hazard in our model, the seismic hazard curves showing a partial contribution of different sources in Murcia, Almeria and Granada are shown in Fig. 9 for PGA and SA (1.0 s). For each city, black lines show the total seismic hazard curve and coloured lines show the seismic hazard curve associated with different sources (zone and faults) for each city.

In Murcia, seismic hazard for short return periods is associated with multiple sources (zone and faults), but for return periods exceeding 475 years (an exceedance probability of 0.1 or lower in 50 years) the seismic hazard is dominated by the Carrascoy Fault. This effect is very similar for PGA and SA (1.0 s).

In Almeria, only two sources, zone 38 and the Carboneras Fault, contribute significantly to seismic hazard. In PGA both sources combine equally to give the seismic hazard for the city, but, for SA (0.1 s), the Carboneras Fault predominates, especially for return periods of more than 475 years.

In Granada, there are many sources contributing to seismic hazard for the city. This is because there are many known faults in its vicinity. Seismic hazard is controlled by zone 35 for PGA and SA (1.0 s) and shorter return periods. This trend changes for return periods greater than 975 years.

Figure 6 PGA for return period of 475 years derived from (a) the proposed hybrid approach and (b) classic zone methodology. Note the fault proximity effects in (a) for these faults: (1) Carboneras Fault, (2) Alhama de Mucia Fault, (3) Carrascoy Fault, (4) Bajo Segura Fault, (5) Mijas and Cartama faults, (6) Zafarraya Fault and (7) Baza Fault.

Figure 10 shows the uniform hazard spectra obtained for four cities in the study area. These graphs can be used to compare the maximum accelerations predicted with the CM and HM in different spectral ordinates, evidencing that the trend observed in PGA persists throughout the entire spectrum.

4 Discussion

We present a hybrid method (HM) for determining a seismic source model that combines zones and faults as independent sources. The HM is based on the distribution of seismic potential among different sources and does not impose any restriction with respect to the type of recurrence model assigned to seismic sources. Moreover, the HM does not require defining a fixed cut-off magnitude M_c that separates the magnitude range in which faults and zones produce earthquakes, as in the works of Frankel et al. (1996).

Figure 7 Comparison of seismic hazard results from the two models (HM and CM) for return periods of (a) 475 and (b) 4975 years.

Figure 8 Seismic hazard curve (with HM and CM) of a site close to the Lietor Fault.

This marks a difference with other approaches that model the seismic potential of faults using two alternatives: (1) single-magnitude rupture models such as Wesnousky's (1986) characteristic earthquake model, as seen in Field et al. (2014) and Akinci and Pace (2017), and (2) models that set a fixed cut-off magnitude and assume that the biggest magnitudes take place in the faults and the smaller ones occur in the zone. In contrast, the formulation of the HM considered above uses a GR-type recurrence model for faults, in line with the proposals made by some other authors (Woessner et al., 2015).

Figure 9 Seismic hazard curve (with HM) of Murcia, Almería and Granada considering all the seismic sources involved. The black lines show the total seismic hazard curve and the coloured lines show the seismic hazard curve associated with different sources (zone and faults).

One strong point of the HM is that it ensures a distribution of seismic potential between faults and zones that prevents double counting of seismicity. This is achieved by computing the seismic potential of faults and zone using the events contained in the completeness period of the catalogue for different magnitude ranges. Identifying the magnitude interval used to distribute seismic potential between zone and faults is of fundamental importance. Specifically, determining M_maxC is crucial for adequately limiting this distribution: a low M_maxC value leads to a notable extrapolation of the recurrence model for faults with large rupture planes; a high M_maxC value may not ensure the complete record of all events of that magnitude. In applying the method to SE Spain the M_maxC value identified is 5.9, which may seem too low. However, this value is consistent with the low–moderate level of seismicity in the study area. In fact, the IGN seismic catalogue (http://www.ign.es/web/ign/portal, last access: 20 October 2018) does not contain any shallow event with magnitude equal to or higher than 6.0 in the instrumental period.

In addition, for the purposes of hazard calculations, the distribution of seismic potential for the entire magnitude interval (between the M_min and the M_max values expected for each source) requires and is dependent on the selected recurrence models to represent fault and zone activities. In this regard, different types of recurrence models may be used: modified Gutenberg–Richter, truncated Gutenberg–Richter (Gutenberg and Richter, 1944), or the models proposed by Main and Burton (1981), Chinnery and North (1975) and others. The use of one model versus another depends on each application and on the available data.

The results of applying the HM are compared with the results of the CM in terms of expected accelerations. A single GMPE is used for both calculations. We have not used any other GMPE (or combination of GMPEs through a logic tree) to simplify the calculations and allow a direct comparison of hazard results.

Figure 10 Uniform hazard spectra obtained in four cities with CM and HM for three return periods.

The results obtained with the HM show an increment of expected accelerations near fault traces (in a factor of 2) in relation to the results of the CM approach. This is consistent with observations of very high ground motions in the epicentral areas of recent earthquakes, such as the 2009 L'Aquila and 2011 Lorca events (Akinci et al., 2010; Cabañas et al., 2014). This increment is achieved at the expense of decreasing expected accelerations in areas located farther away from faults. This is a consequence of the redistribution of seismic potential in the region, which is not increased, but redistributed in several sources (zone and faults).

5 Conclusions

An approach for combining zones and faults in a seismic source model is formulated in this paper.

It is based on the distribution of seismic potential among different sources under certain conditions to prevent counting seismicity twice. Two points of the methodology are critical and must be carefully assessed: the analysis of completeness and the choice of recurrence model used to represent the seismic activity of either source. They are determined by the data available (composition of the seismic catalogue, fault slip rates and geometries, etc.) in the study region, and hence not easily automatized and extendible to other areas. Thus, the approach followed for the application in SE Spain should be reevaluated when it is applied to a different area. For instance, it is to be expected that implementing this approach in a region with rapidly moving faults would produce significantly different results, requiring further adjustments. The higher fault slip rates would imply that the faults consume a larger proportion of the seismic potential available, compromising the convergence of the iterative method to obtain the zone β value.

An initial assumption of the approach is that the seismic-moment potential accumulated in active faults is released only seismically. This condition can be easily modified in the formulation presented above. Additional data informing about other ways of releasing seismic energy, such as slow slip events or aseismic transients, would help to constrain this point.

The seismic hazard map obtained with the HM presents a more heterogeneous aspect compared to the CM seismic hazard map, which assigns a uniform seismic potential to each region. In the HM hazard map, the accelerations expected along fault traces increase and decrease farther away from fault traces, thus keeping the seismic potential budget of the region in balance. This effect can be useful for applications in which the effects of being near a fault must be emphasized, such as urban seismic risk studies for cities located atop active fault planes.

As a final conclusion, we identify some points that require further development and are the focus of an interesting line of research. Specifically, these include (1) determining catalogue completeness for different time-magnitude intervals in the study area; (2) selecting the recurrence model assigned to fault sources according to the data available, and (3) determining the proportion of seismic potential accumulated in faults that is released through earthquakes.