Articles | Volume 18, issue 11
Research article
30 Oct 2018
Research article |  | 30 Oct 2018

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

Alicia Rivas-Medina, Belen Benito, and Jorge Miguel Gaspar-Escribano

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, Mc, 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 Mc 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, Mc, assigning events with a magnitude lower than Mc to the zone and events with magnitude higher than Mc to the faults (Frankel, 1995; Woessner et al., 2015). The question is, how is this Mc value determined? Why can the fault not generate events with a magnitude below the Mc 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 Mc 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 1Diagram 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 2Completeness 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 MMaxC 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, MMaxC, are lower than the observation period (OP) of the catalogue.

  • Magnitude values above MMaxC 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 MMaxC values can be calculated for m<=MMaxC. It is then possible to estimate seismicity rates (Eq. 1) and the seismic moment in each magnitude interval (Eq. 2) as follows:

(1) n ˙ ( m ) = n ( m ) CP ( m ) ,

where n˙(m) 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).

(2) M ˙ o ( m ) = ( m ) M o ( m ) ,

where Mo(m) is the seismic moment released by events of magnitude m, obtained using the equation proposed by Hanks and Kanamori (1979) (log(Mo)=1.5Mw+16.1).

Finally, the cumulative rates in the interval [MMin, MMaxC] can be estimated, where MMin 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 [MMin, MMaxC]=[4.0, 5.9].

Although faults are capable of generating earthquakes with magnitude m>MMaxC, the distribution of seismic potential is carried out in the completeness period [MMin, MMaxC]. 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 [MMaxC, MMax(fault)], where MMax(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 (N˙min) and the cumulative rate of seismic moment (M˙o), for the magnitude range [MMin, MMaxC], in the completeness period CP(m). Details on how to determine N˙min and M˙o for the entire region, the corresponding zone and faults are explained in the following section.

2.1 Seismic potential of the region

The N˙min and M˙o values representing the seismic potential of the region are derived from the seismic catalogue of the magnitude interval [MMin, MMaxC] for the completeness periods, CP(m).


with n˙(m) the annual rate of events with magnitude (m) recorded in CP(m) and Mo(m) the seismic moment for magnitude m. The notation X|MiMj represents the magnitude interval (Mi, Mj) 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):

(5) M ˙ o fault = υ μ A ,

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 μ=3.2×1010 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, MMax(fault). The value MMax(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, M˙ofault can be expressed as follows (Eq. 6):

(6) M ˙ o fault = M m = 0 M Max ( fault ) n ˙ ( m ) M o ( m ) d m ,

where n˙(m) can be estimated applying a recurrence model, for instance, the modified GR model shown in Eq. (7).

(7) n ˙ ( m ) = N ˙ min fault β e - β m e - β M m = 0 - e - β M Max ( fault ) ,

and the seismic moment Mo(m), can be estimated from the Hanks and Kanamori (1979) relation expressed in exponential terms, Mo(m)=edm+c with c=16.1ln(10) and d=1.5ln(10), where m is the moment magnitude Mw (Anderson and Luco, 1983).

Substituting the previous relations in Eq. (6), solving the integral and reordering the equation for N˙minfault, we get Eq. (9).

(8) N ˙ min fault = M ˙ o fault ( d β ) e - β M m = 0 - e - β M Max ( fault ) β e - β M Max M o M Max ( fault ) - e - β M m = 0 M o M m = 0 ,

where Mm=0 is the minimum magnitude that may be generated at a fault rupture (here taken as m=0), MMax(fault) the maximum magnitude, and M˙ofault the seismic-moment rate accumulated in the fault (Eq. 5).

The total seismic-moment rate for each fault (M˙ofault) and seismicity rate (N˙minfault) can be formulated as the sum of the seismic-moment rate released at different magnitude intervals; thus, it follows

(9) N ˙ min fault = N ˙ min M M = 0 M Min + N ˙ min M Min M MaxC + N ˙ min M MaxC M Max ( fault )


By implementing a recurrence model, it is possible to derive the seismicity rate and the moment rate in the interval [MMin, MMaxC] (see example in Fig. 3 with [MMin, MMaxC]=[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 (N˙minfault), as this parameter depends on the seismic-moment rate of each fault (M˙ofault).

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 [MMin, MMaxC] given that


or specifically,

(12) N ˙ min zone M Min M MaxC = N ˙ min region M Min M MaxC - N ˙ min fault M Min M MaxC

(13) M ˙ o zone M Min M MaxC = M ˙ o region M Min M MaxC - M ˙ o fault M Min M MaxC .

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




Regarding the faults, N˙minfault|MMinMMaxC and M˙ofault|MMinMMaxC are derived using an initial (not definitive) β value. Regarding the region, N˙minfregion|MMinMMaxC and M˙oregion|MMinMMaxC are known, as they were extracted from the catalogue (Eqs. 1 and 2). A new additional equation is obtained relating N˙minzone|MMinMMaxC and M˙ozone|MMinMMaxC using Eq. (8) for the interval [MMin, MMaxC] in the zone, resulting in the following:


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 [Mm=o, MMax]. Equation (14) is for zones, and the magnitude interval is restricted to [MMin, MMaxC]. 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 1COV coefficient associated with seismic-moment rate obtained using synthetic catalogues.

Download Print Version | Download XLSX

Figure 4Graph 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 [MMin, MMaxC].

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

Regarding the MMax value expected for the zone (MMax(zone)), this can be considered equal to MMaxC 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 [MMin, MMaxC] 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 5Three-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−1 in a slip rate of ±2 mm yr−1 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−1) for a fault with a slip rate of ±10 mm yr−1 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−1 (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 Mw=4.7 Pedro Muñoz and 2015 Mw=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 μ=3.2×1010 Pa is assumed for the shear modulus (Walters et al., 2009; Martínez-Díaz et al., 2012).

Table 2Seismic rate and seismic-moment rate recorded in the different regions for two magnitude intervals (MMinMMax) and (MMinMMaxC) 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 Mmin to Mmax is contemplated in the magnitude intervals over which hazard is distributed (MMinMMaxC). Note that no faults have been catalogued within regions 28, 29, 33 and 40, which is why no values have been assigned (–).

Download Print Version | Download XLSX

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 Mw (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 MMaxC 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 [MMin, MMaxC] and [MMin, MMax(region)]. It is observed that the seismic potential in the first interval up to MMaxC, constitutes at least a 60 % of the seismic potential in the second interval, up to MMax(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 [MMin, MMaxC] 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 3Parameters 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.

Download Print Version | Download XLSX

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 (Vs30=760 m s−1) conditions.

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

Download Print Version | Download XLSX

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 [MMin, MMaxC] 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 6PGA 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 Mc that separates the magnitude range in which faults and zones produce earthquakes, as in the works of Frankel et al. (1996).

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


Figure 8Seismic 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 9Seismic 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 MmaxC is crucial for adequately limiting this distribution: a low MmaxC value leads to a notable extrapolation of the recurrence model for faults with large rupture planes; a high MmaxC value may not ensure the complete record of all events of that magnitude. In applying the method to SE Spain the MmaxC 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 (, 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 Mmin and the Mmax 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 10Uniform 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.

Data availability

The available data of active faults are in the Supplement.


The supplement related to this article is available online at:

Author contributions

All authors contributed to the preparation of this paper.

Competing interests

The authors declare that they have no conflict of interest.


We would like to thank Mario Ordaz Schroeder for his time and support during a research stay carried out by ARM at the Instituto de Ingeniería, UNAM, and we acknowledge financial support from Vicerrectoría de Investigación y Desarrollo (VRID), UdeC, “216.419.003-1.0IN”.

Edited by: Maria Ana Baptista
Reviewed by: five anonymous referees


Akinci, A. and Pace, B.: Effect of Fault Parameter Uncertainties on PSHA explored by Monte Carlo Simulations: A case study for southern Apennines, Italy, in: AGU Fall Meeting Abstracts, 11–15 December 2017, New Orleans, 2017. 

Akinci, A., Malagnini, L., and Sabetta, F.: Characteristics of the strong ground motions from the 6 April 2009 L'Aquila earthquake, Italy, Soil Dynam. Earthq. Eng., 30, 320–335, 2010. 

Anderson, J. G. and Luco, J. E.: Consequences of slip rate constraints on earthquake occurrence relations, B. Seismol. Soc. Am., 73, 471–496, 1983. 

Argus, D. F., Gordon, R. G., DeMets, C., and Stein, S.: Closure of the Africa–Eurasia–North America plate motion circuit and tectonics of the Gloria Fault, J. Geophys. Res., 94, 5585–5602, 1989. 

Basili, R., Tiberti, M. M., Kastelic, V., Romano, F., Piatanesi, A., Selva, J., and Lorito, S.: Integrating geologic fault data into tsunami hazard studies, Nat. Hazards Earth Syst. Sci., 13, 1025–1050,, 2013. 

Benito, M. B. and Gaspar-Escribano, J. M.: Ground Motion Characterization in Spain: Context, Problems and Recent Developments in Seismic Hazard Assessment, J. Seismol., 11, 433–452, 2007. 

Benito, M. B., Navarro, M., Vidal, F., Gaspar-Escribano, J., García-Rodríguez, M. J., and Martínez-Solares, J. M.: A new seismic hazard assessment in the region of Andalusia (Southern Spain), B. Earthq. Eng., 8, 739–766, 2010. 

Brune, J. N.: Seismic moment, seismicity, and rate of slip along major fault zones, J. Geophys. Res., 73, 777–784,, 1968. 

Bungum, H.: Numerical modelling of fault activities, Comput. Geosci., 33, 808–820, 2007. 

Cabañas Rodríguez, L., Carreño Herrero, E., Izquierdo Álvarez, A., Martínez Solares, J. M., Capote del Villar, R., Martínez Díaz, J., Benito Oterino, B., Gaspar Escribano, J., Rivas Medina, A., García Mayordomo, J., Pérez López, R., Rodríguez Pascua, M. A., and Murphy Corella, P.: Informe del sismo de Lorca del 11 de mayo de 2011, Instituto Geográfico Nacional, Madrid, p. 138, 2011. 

Cabañas, L., Alcalde, J. M., Carreño, E., and Bravo, J. B.: Characteristics of observed strong motion accelerograms from the 2011 Lorca (Spain) Earthquake, B. Earthq. Eng., 12, 1909–1932, 2014. 

Cabañas, L., Rivas-Medina, A., Martínez-Solares, J. M., Gaspar-Escribano, J. M., Benito, B., Antón, R., and Ruiz-Barajas, S.: Relationships between Mw and Other Earthquake Size Parameters in the Spanish IGN Seismic Catalog, Pure Appl. Geophys., 172, 2397–2410,, 2015. 

Campbell, K. W. and Bozorgnia, Y.: NGA-West2 ground motion model for the average horizontal components of PGA, PGV, and 5 % damped linear acceleration response spectra, Earthq. Spectra., 30, 1087–1115,, 2014. 

Chinnery, M. A. and North, R. G.: The frequency of very large earthquakes, Science, 190, 1197–1198, 1975. 

Crespo, M. J., Martínez, F., and Martí, J.: Seismic hazard of the Iberian Peninsula: evaluation with kernel functions, Nat. Hazards Earth Syst. Sci., 14, 1309–1323,, 2014. 

DISS Working Group: Database of Individual Seismogenic Sources (DISS), Version 3.2.1, A compilation of potential sources for earthquakes larger than M 5.5 in Italy and surrounding areas, Istituto Nazionale di Geofisica e Vulcanologia,, 2018. 

Dixon, T. H., Norabuena, E., and Hotaling, L.: Paleoseismology and Global Positioning System: Earthquake-cycle effects and geodetic versus geologic fault slip rates in the Eastern California shear zone, Geology, 31, 55–58, 2003. 

Field, E. H., Arrowsmith, R. J., Biasi, G. P., Bird, P., Dawson, T. E., Felzer, K. R., Jackson, D. D., Johnson, K. M., Jordan, T. H., Madden, C., Michael, A. J., Milner, K. R., Page, M. T., Parsons, T., Powers, P. M., Shaw, B. E., Thatcher, W. R., Weldon, R. J., and Zeng, Y.: Uniform California Earthquake Rupture Forecast, Version 3 (UCERF3) – The Time-Independent Model, B. Seismol. Soc. Am., 104, 1122–1180,, 2014. 

Frankel, A. D., Mueller, C., Barnhard, T., Perkins, D., Leyendecker, E., Dickman, N., Hanson, S., and Hopper, M.: National seismic-hazard maps: documentation June 1996, US Geological Survey, Denver,, 1996. 

García-Mayordomo, J.: Caracterización y Análisis de la Peligrosidad Sísmica en el Sureste de España, PhD Thesis, University Complutense of Madrid, Madrid, Spain, 2005. 

García-Mayordomo, J., Gaspar-Escribano, J. M., and Benito, B.: Seismic hazard assessment of the Province of Murcia (SE Spain): analysis of source contribution to hazard, J. Seismol., 11, 453–471, 2007. 

García-Mayordomo, J., Insua-Arévalo, J. M., Martínez-Díaz, J. J., Perea, H., Álvarez-Gómez, J. A., Martín-González, F., González, A., Lafuente, P., Pérez-López, R., Rodríguez-Pascua, M. A., Ginez-Robles, J., and Azañón, J. M.: Modelo integral de zonas sismogénicas de España, in: Contribución de la Geología al Análisis de la Peligrosidad Sísmica, edited by: Insua-Arévalo, J. M. y Martín-González, J. J., Sigüenza, Guadalajara, España, 193–196, 2010. 

García-Mayordomo, J., Insua-Arévalo, J. M., Martínez-Díaz, J. J., Jiménez-Díaz, A., Martín-Banda, R., Martín-Alfageme, S., Álvarez-Gómez, J. A., Rodríguez-Peces, M., Pérez-López, R., Rodríguez-Pascua, M. A., Masana, E., Perea, H., Martín-González, F., Giner-Robles, J., Nemser, E. S., and Cabral, J.: The quaternary active faults database of iberia (QAFI v.2.0)/La base de datos de fallas activas en el cuaternario de iberia (QAFI v.2.0), J. Iber. Geol., 38, 285–302, 2012. 

Gaspar-Escribano, J. M., Rivas-Medina, A., Parra, H., Cabañas, L., Benito, B., Barajas, S. R., and Solares, J. M.: Uncertainty assessment for the seismic hazard map of Spain, Eng. Geol., 199, 62–73,, 2015. 

Gutenberg, B. and Richter, C. F.: Frequency of earthquakes in California, B. Seismol. Soc. Am., 34, 185–188, 1944. 

Hakimhashemi, A. H. and Grünthal, G.: A statistical method for estimating catalog completeness applicable to long-term nonstationary seismicity data, B. Seismol. Soc. Am., 102, 2530–2546, 2012. 

Hanks, T. C. and Kanamori, H.: A moment magnitude scale, J. Geophys. Res., 84, 23480–23500, 1979. 

IGN-UPM working group: Actualización de mapas de peligrosidad sísmica de España 2012, Centro Nacional de Información Geográfica, Madrid, 267 pp., 2013. 

Kijko, A., Smit, A., and Sellevoll, M. A.: Estimation of earthquake hazard parameters from incomplete data files. Part III. Incorporation of uncertainty of earthquake-occurrence model, . Seismol. Soc. Am., 106, 1210–1222, 2016. 

Kiratzi, A. A. and Papazachos, C. B.: Active crustal deformation from the Azores triple junction to the Middle East, Tectonophysics, 243, 1–24, 1995. 

Langbein, J. and Bock, Y.: High-rate real-time GPS network at Parkfield: Utility for detecting fault slip and seismic displacements, Geophys. Res. Lett., 31, L15S20,, 2004. 

Leonard, M.: Earthquake Fault Scaling: Self-Consistent Relating of Rupture Length, Width, Average Displacement, and Moment Release, B. Seismol. Soc. Am., 100, 1971–1988,, 2010. 

Main, I. G. and Burton, P. W.: Rates of crustal deformation inferred from seismic moment and Gumbel's third distribution of extreme magnitude values, in: Earthquakes and Earthquake Engineering: The Eastern United States, edited by: Beavers, J. E., Ann Arbor Science, Michigan, 937–951, 1981. 

Martínez-Díaz, J. J., Bejar-Pizarro, M., Álvarez-Gómez, J. A., de Lis Mancilla, F., Stich, D., Herrera, G., and Morales, J.: Tectonic and seismic implications of an intersegment rupture: The damaging May 11th 2011 Mw 5.2 Lorca, Spain, earthquake, Tectonophysics, 546, 28–37, 2012. 

Metzger, S., Jónsson, S., and Geirsson, H.: Locking depth and slip-rate of the Húsavík Flatey fault, North Iceland, derived from continuous GPS data 2006–2010, Geophys. J. Int., 187, 564–576, 2011.  

Mezcua, J., Rueda, J., and Blanco, R. M. G.: A new probabilistic seismic hazard study of Spain, Nat. Hazards, 59, 1087–1108, 2011. 

Nasir, A., Lenhardt, W., Hintersberger, E., and Decker, K.: Assessing the completeness of historical and instrumental earthquake data in Austria and the surrounding areas, Aust. J. Earth. Sci., 106, 90–102, 2013. 

Ordaz, M., Martinelli, F., D'Amico, V., and Meletti, C.: CRISIS2008: A Flexible Tool to Perform Probabilistic Seismic Hazard Assessment, Seismol. Res. Lett., 84, 495–504,, 2013. 

Papanikolaou, I. D., Roberts, G. P., and Michetti, A. M.: Fault scarps and deformation rates in Lazio–Abruzzo, Central Italy: Comparison between geological fault slip-rate and GPS data, Tectonophysics, 408, 147–176, 2005. 

Parsons, T. and Geist, E. L.: Is there a basis for preferring characteristic earthquakes over a Gutenberg–Richter distribution in probabilistic earthquake forecasting?, B. Seismol. Soc. Am., 99, 2012–2019, 2009. 

Peláez, J. A. and Lopez-Casado, C.: Seismic Hazard estimate at the Iberian Peninsula, Pure Appl. Geophys., 159, 2699–2713, 2002. 

Salgado-Gálvez, M. A., Cardona, O. D., Carreño, M. L., and Barbat, A. H.: Probabilistic seismic hazard and risk assessment in Spain, Centro Internacional de Métodos Numéricos en Ingeniería, Barcelona, 2015. 

Stirling, M., Rhoades, D., and Berryman, K.: Comparison of earthquake scaling relations derived from data of the instrumental and preinstrumental era, B. Seismol. Soc. Am., 92, 812–830,, 2002. 

Walpersdorf, A., Manighetti, I., Mousavi, Z., Tavakoli, F., Vergnolle, M., Jadidi, A., Hatzfeld, D., Aghamohammadi, A., Bigot, A., Djamour, Y., Nankali, H., and Sedighi, M.: Present-day kinematics and fault slip rates in eastern Iran, derived from 11 years of GPS data, J. Geophys. Res.-Solid, 119, 1359–1383, 2014. 

Walters, R. J., Elliott, J. R., D'agostino, N., England, P. C., Hunstad, I., Jackson, J. A., Parsons, B., Phillips, R. J., and Roberts, G.: The 2009 L'Aquila earthquake (central Italy): A source mechanism and implications for seismic hazard, Geophys. Res. Lett., 36, L17312,, 2009. 

Wells, D. L. and Coppersmith, K. J.: New empirical relationships among magnitude, rupture length, rupture width, rupture area, and surface displacement, B. Seismol. Soc. Am., 84, 974–1002, 1994. 

Wesnousky, S. G.: Earthquakes, Quatemary faults, and seismic hazard in California, J. Geophys. Res., 91, 12587–12631, 1986. 

Woessner, J., Laurentiu, D., Giardini, D., Crowley, H., Cotton, F., Grünthal, G., Valensise, G., Arvidsson, R., Basili, R., Demircioglu, M. B., Hiemer, S., Meletti, C., Musson, R. W., Rovida, A. N., Sesetyan, K., Stucchi, M., and The SHARE Consortium:: The 2013 European seismic hazard model: key components and results, B. Earthq. Eng., 13, 3553–3596, 2015. 

Short summary
We present an approach for seismic hazard assessment that considers a hybrid source model composed of faults and zones containing the remaining seismicity. The seismic-moment rate is used to distribute seismic potential, avoiding double counting. The approach is applied in SE Spain, a region of low-to-moderate seismicity. Results show a concentration of expected accelerations around fault traces using the hybrid approach, which is not appreciated in the classic approach using zones exclusively.
Final-revised paper