Approach for combining faults and area sources in seismic hazard assessment : Application in southeastern Spain

This paper presents a methodological approach for seismic hazard assessment that considers a hybrid source model composed by faults as independent entities and zones (containing the 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 10 catalog. For fault sources, it is inferred from kinematic parameters derived from paleoseismicty and GNSS measurements. Distributing the seismic potential associated to each source is a key question when considering hybrid models of zone and faults, which some authors solve by assigning to the fault only the earthquakes that exceed a fixed magnitude value Mc. In the present approach, instead of restricting the magnitudes of each type of source, the distribution of seismic potential is carried out only for magnitudes below the maximum magnitude value completely recorded in the catalog (Mmaxc). This is 15 derived from a completeness analysis and can be lower than the Mmax generated by the faults, taking into account that their the recurrence period can be higher than the observation period of the catalog. The proposed approach is applied in southern Spain, a region of low-to-moderate seismic where faults move slowly. The results obtained are contrasted with the results of a seismic hazard model using the traditional zone model exclusively. Results show a concentration of expected accelerations around fault traces using the hybrid approach, which is not 20 appreciated in the classic approach using zones exclusively.


Introduction
Active faults are the main earthquake sources in the crust.However, their incorporation in seismic hazard assessment is not straightforward due to the lack of data required to model them appropriately.This leads to a limited use of faults as independent sources in seismic hazard analyses and to an extended use of seismic zones, which cover a significant portion of the crust, assuming uniform seismic characteristics within each source.This situation is changing in the last years, as more studies on active tectonics, paleoseismicity and fault deformation rates derived from GNSS measurements (among others) constraining fault parameters (such as rupture plane geometry, predominant sense of slip, slip rates, recurrence periods, etc.) are available.The consideration of sources as "fault type" instead of zones in seismic hazard studies requires addressing two factors: the 3D geometry of the source and the data needed to characterize its seismic potential.In most practical cases, the characterization of the seismic potential of faults is based on their kinematics.Other approaches, such as extracting the seismic parameters of each single fault from the earthquake catalog is not always viable, especially in areas with slowly moving faults.Additionally, the period considered in the catalogue may be too short in comparison with the recurrence time of the fault to provide unbiased estimated of fault seismic parameters.
In principle, the modeling of 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 quite 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 generating low magnitude events or slow faults producing rare events cannot be properly characterized.To prevent a possible deficit in the seismic source model of a region, the use of faults as seismic sources may be completed with zones that account for the seismic potential linked to these small/slow faults or just unknown faults that cannot be characterized independently.Hence, we propose to consider hybrid source models, composed by faults and zones: the first ones modelled as independent sources and the second ones encompassing the residual seismicity.
The problem of using a model combining zones and faults is establishing the distribution of seismic potential among them, taking into account that they are derived from different data sources.For zones, the recurrence model is calculated from the seismic catalog, and for faults, it is derived from fault geometries and slip rate estimates (based on GNSS-measured deformation rates).The problem is that a part of the events contained in the catalog must be linked to faults and they are already included in the seismic potential of faults derived from slip rate estimates.If all events are assigned to the zone, the events linked to faults would be double-counted, leading to an overestimate of the total seismic potential (related to both faults and zones).The approach presented in this paper, as all probabilistic seismic hazard models, face the challenging question of estimating the expected ground motions with the basis of a short period of observations of earthquake occurrences and limited geological data (with significant uncertainty) to construct recurrence models.The purpose of this work is not to solve this challenge, but rather, to propose a model that contains different types of seismic sources (faults and zones) and distributes the seismic potential appropriately, avoiding double-counting and considering periods of completeness.

Source hybrid approach (zones & faults) for hazard estimation
The hybrid model proposed is composed by fault-type sources and zone-type sources.At the same time, the term region is defined as the geometric envelope of both source types.Thus, the region presents the same geometry as the zone and the seismic potential (seismicity rate and seismic moment rate) of both sources (faults and zone).The zone is defined as the source which seismic potential is residual, excluding the seismic potential of faults.Although there is a geometrical equivalence between region and zone, their seismic potential are very different, as the seismic potential of the region equals the seismic potential of the zone plus the seismic potential of the faults enclosed in the region (Fig. 1).
The problem is then distributing the seismic potential of the region between the zone and the faults trying to avoid doublecounting.The following considerations are made: -The seismicity rate of the region is derived from the seismic catalog using the records which origin time and magnitude are contained in the period for which the catalog can be considered complete, CP(m) (complete period for magnitude m and higher).
-The complete periods PC(m) (for different magnitude up to a maximum magnitude of completeness value, M MaxC. ) are included in the observation period (OP) of the catalog.
-Magnitudes values above M MaxC present recurrence periods higher than the catalog OP or constitute a sample an insufficient number of records to apply statistics.
The representation of the number of recorded events for different magnitude intervals as a function of time allows the identification of the reference years RY(m) for different magnitude intervals using the Stepp (1972) approach.Consequently, the PC(m) and M MaxC values can be calculated for m<= M MaxC .Then, it is possible to estimate seismicity rates and seismic moment in each magnitude interval and cumulative rates in the interval [M Min , M MaxC ], where M Min is the minimum magnitude value used to compute seismic hazard.This is illustrated with an example in Fig. 2, with [M Min , M MaxC ] = [4.0,6.9].
Although faults are capable of generating earthquakes with magnitude m > M MaxC , the distribution of seismic potential is carried out in the complete period [M Min , M MaxC ].This way we avoid miscalculating problems related to recurrence periods larger than OP and the eventual occurrence of events with magnitudes higher than the maximum expected magnitudes in the fault (M Max ) that are not contained in the seismic catalog.The computation of the seismic potential of the fault in the interval [M MaxC , M Max ] (where M Max 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 cumulative rate of earthquakes (Ṅmin) and the cumulative rate of seismic moment (Ṁ), for the magnitude range [M Min , M MaxC ], in the complete period PC(m).Details on how to determine Ṅmin and Ṁ for the entire region, the corresponding zone and faults are explained in the following.

Seismic potential of the region
The Ṅmin and Ṁ values representing the seismic potential of the region are derived from the seismic catalog in the complete periods CP(m), for the magnitude interval [M Min , M MaxC ]. (1) with ṅ(m) the annual rate of events with magnitude (m) recorded in PC(m) and Mo(m) the seismic moment for magnitude m.

Seismic potential of faults
The cumulative moment rate of the fault is estimated from the slip rate (), the rigidity modulus (µ) and the area of the fault plane (A), using the expression of Brune (1968), assuming that the fault plane is accumulating energy evenly: This moment rate represents the average seismic moment accumulated in the fault per year that will be released by earthquakes of different magnitudes (from very low, close to m=0 up to the maximum magnitude of the fault, M Max ).The value of M Max can be evaluated from a geometrical element of the fault plane using empirical relationships proposed in the literature, as Wells and Coppersmith (1994), Stirling et al. (2002) or Leonard (2010) (among others).
Assigning a recurrence model to the fault and assuming an initial b-value (for example, b=1) that will be adjusted iteratively, it is possible to stablish a relation between the seismicity rare and the moment rate by solving the integral proposed by Anderson (1979).Using a modified Gutenberg -Richter recurrence model, N ̇min is calculated as: The total seismic moment rate of each fault (Ṁ  ) and seismicity rate (N ̇min) can be formulated as the sum of the seismic moment rate released in different magnitude intervals, it follows: (5)

Seismic potential of the zone
The parameters representing the zone are initially unknown.The can be calculated for the interval [M Min , M MaxC ] considering that: Or specifically: In principle, there are 2 equations with two unknowns related to the zone:  Considering that the fault may generate events with magnitudes larger than M MaxC , the corresponding distribution of seismic potential in the interval (M MaxC , M Max ] is calculated by extrapolation of the recurrence model with the last b-value adjusted (Fig. 4).
Regarding the M Max value expected for the zone, this can be equalled to M MaxC or extended to a higher magnitude value if it is accepted that bigger events can occur in other unidentified sources (such as blind faults).
unidentified sources (such as blind faults).

Analysis of uncertainty
The proposed approach strongly relies on the computation of the seismicity (earthquake rates and moment rate) within the magnitude interval [M Min , M MaxC ] of the seismic catalog that contains the complete record of events occurred in the entire region.In order to capture the variability of seismic moment rates calculated from the earthquake catalog 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 catalog and 3) the proportion of earthquakes of different magnitude (b-value).
Synthetic catalogs derived from GR-modified recurrence models are generated for this purpose.Earthquake rates are computed using different amounts of events, magnitude intervals and b-values that could be representative from areas of low-to-moderate seismic activity.
The procedure comprises five steps: − Generation of 2000 synthetic catalogs for different combinations of earthquake rates, magnitude intervals and bvalues.
− Calculation of earthquake rates for different magnitude values, for each synthetic catalog (Eq.( 4)).
− Calculation of moment rate different magnitude values, for each synthetic catalog.
− Sum of moment rate different magnitude values, to obtain the cumulative moment rate for each synthetic catalog (Eq.( 5)).
− Computation of the mean and the standard deviation of the distribution of calculated seismic moment rates.
Table 1 shows the coefficient of variation (COV = sigma/mean) associated to each combination (number of events,

Application of the approach in southeast Spain
The approach described above is applied in southeastern Spain, the most seismically active area of the country.The tectonic deformation and seismicity is related with the convergence between the African and Eurasian plates, with approximate shortening rate of 5 mm/year in NNW-SSE direction.Crustal deformation is accumulated over a broad area where the seismicity is diffuse.The assignation of earthquakes to specific faults is not an easy task: whereas earthquakes can be clearly associated to a rupture (such as the 2011 M 5.2 Lorca event generated in the Alhama de Murcia fault system), other events occurred in areas with unknown active faults (as the 2007, Mw 5.1 Pedro Muñoz and 2015, Mw 4.7 Ossa de Montiel earthquakes, both located in Central Spain).
These seismotectonic characteristics makes southeastern Spain a suitable scenario for applying the hybrid approach combining faults and zones proposed in this work.The recently available data on fault geometries and activity parameters in Spain and the seismic catalog that contain records dating back to 1048 are used to constrain the source model parameters.

Sources input data
The seismogenic source model considered for SE Spain is composed by 12 regions that contain a total of 95 faults.Active fault data are taken from the QAFI database (v2.0), which includes information about faults segmentation, geometry and slip rate (see Fig. 5).The maximum expected earthquake magnitude generated in each fault is derived from the fault geometry, applying the empirical relationship of Stirling (2002).Moment rates accumulated in the faults are estimated using the fault plane area and the slip rate value by the formula of Brune (1968).The zones model of 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).Not all the regions considered contain faults sources, i. e., 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 catalog recently homogenized (IGN-UPM Working Group, 2013; Cabañas et al., 2015).It contains 1496 earthquakes, with magnitudes ranging from 4.0 to 6.6.The uncertainty assessment of the catalog is explained in Gaspar-Escribano et al. (2015).A M MaxC is estimated in 5.9 for SE Spain (although not all the regions reach this maximum magnitude value).
Table 2 shows the seismic potential for each region, calculated in the magnitude intervals [M Min , M MaxC ] and [M Min , M Max ].It is appreciated 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.
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 is not significant and the [M Min , M MaxC ] interval is very narrow, resulting in an increased uncertainty in the hazard estimates in this region.A GR-mod recurrence model is also assigned to the faults.Finally, the distribution of seismic moment among all seismic sources is carried out (Table 4).As it can be observed, the seismic moment rate associated to the zone has a strong influence on the estimated seismic hazard of the region.This is due to the reduced amount of known active faults than can be modelled as independent sources, a common situation in areas with low and moderate seismic activity.However, it is noteworthy that the seismic potential of regions 35, 36 and 38 is dominated by the seismic potential of faults.
The seismic hazard calculation is done with the software CRISIS2012 (Ordaz et al., 2013), considering the strong motion equation of Campbell (2013), which allows including the fault geometry and the style of faulting.The ground motion parameters predicted include peak ground acceleration (PGA) and 15 spectral accelerations in the period range (0.05 -10 seconds), all obtained in hard soil (Vs30=760 m/s) conditions.

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 6b, respectively.Only the geometry of the zone model differs from both analyses, remaining the ground motion prediction equations (GMPEs) and the other calculation parameters in all the estimates.
The seismic hazard map resulting from the CM approach is rather homogeneous: there are no big differences in the spatial distribution of expected acceleration values.This result is characteristic of this method because the seismic sources are very big.PGA estimates for return period of 475 years using the zone approach (CM) reach maximum values in Granada, Almeria and the Murcia region, around 0.20 g.The minimum PGA values are obtained in Jaen, with values as low as 0.06 g.
Fig. 6a shows the seismic hazard map resulting from the application 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.30g) and, to a minor extent, by the Venta de Zafarraya, Carrascoy, Bajo Segura, Baza, Mijas and Cartama fault systems.
The seismic hazard map obtained through the HM displays more spatial variability, showing maximum values along fault sources that decrease strongly away from faults.This trend reflects a "source proximity effect", implying higher acceleration values in 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).These results agree with the real observations in recent earthquakes, including Lorca 2011.
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 7b, respectively.The tendency 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 small seismic activity and large faults with low slip rates.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 very much, giving 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.The hazard curves associated with the two models (HM and CM) in this region are shown in In Murcia, seismic hazard for short return periods corresponds to several sources involved (zone and faults), but for return periods exceeding 475 years (probability of exceedance of 0.1 and 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 contribute to the seismic hazard significantly, zone 38 and the Carboneras fault.In PGA both sources combine equally to give the seismic hazard of 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 the seismic hazard of 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 tendency changes return periods greater than 975 years.

Discussion and conclusions
The estimation of expected ground motions for a future time span is the purpose of a seismic hazard assessment.Such assessment is specially challenging in areas with low seismic activity and faults with long recurrence periods.In these areas, the seismic record contains many information gaps, including the estimation of earthquake rates, seismic recurrences and maximum magnitude values.Additionally, the identification and characterization of seismogenic, active faults are not evident.
This paper proposes an approach to prepare a hybrid source model for use in seismic hazard assessment, especially suitable for low seismicity areas.For a given region, the model includes the faults that can be adequately characterized for hazard applications as independent seismic sources.The remaining seismicity, related to other unidentified or unconstrained faults, is taken into account trough seismic zones where the seismicity is uniformly distributed.Both types of seismic sources are modelled with their corresponding recurrence models and maximum expected magnitudes, which are constrained with different kind of data.Fault geometries and slip rates are used to assess the seismic potential of faults.Earthquake rates derived from the seismic catalogue and magnitudes complemented with geological data are considered to constraint the seismic potential of the zone.To prevent double counting, a strong point of this approach is that it assures a distribution of seismic potential between faults and zones within the complete period of the catalog for different magnitude ranges.This distribution of seismic potential depends on the selection of recurrence models to represent fault and zone activities.This is one of the controlling points that have to be carefully addressed considering the specific information available for the study area.The application of the hybrid model in SE Spain takes advantage of the available (still incomplete) fault data and the relatively long catalog (although lacking high magnitude events).The GR-mod model is used for faults and zones and provides reasonable results in the SE Spain.This approach is particularly suitable for areas of moderate seismicity such as SE Spain, where we focused the application.The results show an increment of expected accelerations near fault traces (in a factor of 2).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.This increment is achieved at the expense of decreasing expected accelerations in areas located farther away from faults.Therefore, this approach gives a more important role to active faults than the classic seismic hazard assessment based on the zone model exclusively.Thus, it can be useful for applications in which near fault effects need to be highlighted, such as urban seismic risk studies for cities located atop active fault planes.
The performance of the HM needs to be proven in areas with higher seismic activity, where the budget of observed seismic potential available for distribution between faults and zones can be consumed by the very active faults.In this case it is necessary to analyze different recurrence models, to properly assess catalog completeness and to establish the MmaxC value.
This important task will be the focus of future studies.
Some authors propose a simple way of distributing the seismic potential based on a uniform magnitude value Mc, and assigning the events with magnitude lower than Mc to the zone and the events with magnitude higher than Mc to the faults.The question is: how is this Mc value determined?Why the fault cannot generate events with magnitude below the Mc value in possible blind faults.In a study region such as southern Spain, with slow faults and maximum magnitudes around 6.5-7.0,fixing a Mc value results certainly complicated.Hence, a new approach to distribute the seismic potential is proposed in this paper.
10) Note that the b-value of the zone appears in this equation can be equaled to the b-value of the region as both sources present similar seismic nature.With this third equation, it is possible to solve the system and obtaining a new b value for the fault (second iteration) that balances the three equations.This gives the distribution of seismic potential between the zone and the faults in the interval [M Min , M MaxC ].
magnitude interval, b-value).As can be seen, the greater the number of records in the sample and lower magnitude range, the lower the uncertainty associated with the rate of seismic moment calculated.The b-value presents a different trend, recording the biggest variability for b-values between 1.0 and 1.5.This table is useful to estimate the uncertainty on the seismic moment rate calculated from the seismic catalogue, as a function of the number of earthquakes, magnitude interval and bvalue.Is important also to consider the uncertainty associated with the slip rate and the area of the fault, as they propagate into the distribution of seismic moment rate of the fault proportionally 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/year in a slip rate of ±2 mm/year represents an uncertainty of 50%, leading to a COV value of 0.5 in the moment rate of the fault.However, the same deviation (± 1 mm/year) for a fault with slip rate of ±10 mm/year, represents an uncertainty of 10%, leading to a COV coefficient moment rate of the fault of only 0.1.Nat.Hazards Earth Syst.Sci.Discuss., https://doi.org/10.5194/nhess-2018-28Manuscript under review for journal Nat.Hazards Earth Syst.Sci. Discussion started: 13 March 2018 c Author(s) 2018.CC BY 4.0 License.

Fig. 8 .
Fig.8.The HM hazard curve reflects a substantial increase in hazard for long return periods.

Figure 1 .
Figure 1.Diagram with the distribution of the seismic potential of a region different sources.

Figure 2 .
Figure 2. Graphs with the completeness analyses of the seismic catalog to identify the completeness period for magnitude PC(m) an MMaxC value.

Figure 3
Figure 3 Graphs with the identification of seismicity rate and seismic moment rate of a fault, using a modified Gutenberg-Richter recurrence model.

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

Figure 5 .
Figure 5. 3D view of the seismic sources considered for hazard calculation, including faults and zones.

Figure 7 .
Figure 7. Seismic Hazard Results Comparison of the two models (HM / CM) for return periods of (a) 475 and (b) 4975 years.

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

Figure 9 .
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 colour lines the seismic hazard curve associated with different sources (zone and faults).
Ṅmin zone | M  M  and Ṁo zone | M  M  .Regarding the fault, Ṅmin fault | M  M  and Ṁo fault | M  M  are derived using an initial (not definitive) b-value.An new additional equation is obtained relating Ṅmin zone | M  M  and Ṁo zone | M  M using Eq.(4) for the interval [M Min , M MaxC ] in the zone, giving: