Fault System-Based Probabilistic Seismic Hazard Assessment of a Moderate Seismicity Region: The Eastern Betics Shear Zone (SE Spain)

Including faults as seismogenic sources in probabilistic seismic hazard assessments (PSHA) has turned into a common practice as knowledge of active faults is improving. Moreover, the occurrence of earthquakes in multi-fault ruptures has evidenced the need to understand faults as interacting systems rather than independent sources. We present a PSHA for the Southeastern Spain obtained by including the faults of a moderate seismicity region, the Eastern Betics Shear Zone (EBSZ) in SE Spain, as the main seismogenic sources in two separate source models, one considering background seismicity. In contrast with previous studies in Spain, earthquake occurrence of the EBSZ system is modeled considering different hypotheses of multi-fault ruptures at the whole fault system scale and weighted in a logic tree. We compare the hazard levels with those from an area source PSHA and a previous fault-based approach. The results show a clear control of the EBSZ faults in the seismic hazard for all return periods, increasing drastically the hazard levels in the regions close to the fault traces and influencing up to 20 km farther with respect to the area source PSHA. The seismic hazard is dependent on the fault slip rates as peak ground accelerations and territorial extension of the fault influence appear higher around the Alhama de Murcia and Carboneras faults, while lower slip rate faults (Palomares Fault) show minor contribution to the hazard. For the return period of 475 years and near-fault locations, our models are more consistent with the ground motion values reached in the 2011 Mw 5.2 Lorca event than the building code or national seismic hazard map, which suggest that our fault system-based model performs more accurate estimations for this return period. Fault data, mainly slip rates, and its uncertainties have a clear impact on the seismic hazard and, for some faults, the lack of detailed paleoseismic studies can compromise the reliability of the hazard estimations. This, together with epistemic uncertainties concerning the background seismicity, are key discussion points in the present study, having an impact on further research and aiming to serve as a case example for other low-to-moderate seismicity regions worldwide.


INTRODUCTION
The use of geological data in probabilistic seismic hazard assessment (PSHA) has become increasingly important, since seismic catalogs do not usually extend for more than a few centuries. This is a very short time window for characterizing the occurrence of major events in regions of moderate seismic activity, even in areas where the catalogs extend for more than 1,000 years such as SE Spain or Portugal (García-Mayordomo, 2005;Vilanova and Fonseca, 2007, respectively). Traditionally, geological data and fault data have been used to characterize the occurrence of major, maximum or characteristic events (Wesnousky, 1986;Schwartz and Coppersmith, 1984). Maximum events usually have a very low impact in seismic hazard for the return periods considered in conventional building (e.g., 475 years) and so this type of fault-based PSHA has been regarded in studies interested in long return periods (e.g., critical infrastructures). However, there is increasing awareness of the fact that low-to-moderate magnitude events (Mw 4.0-6.0) are also able to produce significant damage, especially at sites located on top or very close to fault sources (e.g., 2011 Mw 5.2 Lorca earthquake in SE Spain; Alguacil et al., 2014). In this context, the identification and characterization of active faults as seismogenic sources able to produce all the range of earthquake magnitudes become critical tasks in the PSHA of moderate active regions.
For this reason, in accordance with the increasing availability and refinement of paleoseismic and geological studies, the incorporation of fault data in PSHA has become an important part of the source modeling aiming to estimate earthquake activity beyond the coverage of seismic catalogs. The use of faults in PSHA is a well-established practice worldwide (Frankel, 1995;WGCEP, 2003;Stirling et al., 2012;Woessner et al., 2015). In most of these cases, faults are modeled as independent sources following a characteristic earthquake model (Youngs and Coppersmith, 1985;Wesnousky, 1986) or a Gutenberg-Richter distribution (Gutenberg and Richter, 1945;Bungum, 2007).
During earthquakes, faults can rupture in very complex configurations involving several structures simultaneously and even with contrasting kinematics or geometries. Earthquakes such as the 2010 Mw 7.2 Cucapah (Fletcher et al., 2014), the 2012 Mw 8.6 Sumatra (Zhang et al., 2012) or the 2016 Mw 7.8 Kaikoura (Little et al., 2018) are good examples. Accordingly, in the recent years seismic hazard modelers and practitioners have developed new approaches for modeling faults as complex and interacting sources in seismic hazard. Including this complexity has demonstrated to be a turning point to model fault processes in PSHA. Approaches like UCERF-3 (Field et al., 2014) or the recent SHERIFS (Chartier et al., 2019) are emerging methodologies whose principle is modeling earthquake occurrence considering scenarios of multi-fault ruptures for source characterization in PSHA. In both cases fault rupturing is treated as an aleatory variable, whose occurrence is assumed to be linked to the randomness of the seismic process.
Most fault-based PSHA approaches are usually carried out in high seismicity regions with advanced paleoseismic fault knowledge to perform such studies (e.g., Stirling et al., 2012;Field et al., 2014 for New Zealand and California, respectively). Conversely, the scarcity of paleoseismic and geological slip rate data on faults has been one of the major challenges faced by PSHA analysts in low-to-moderate seismicity regions, such as Western Europe (e.g., Spain). In Spain, PSHA has been classically modeled using area sources following Cornell (1968) and McGuire (1976) methodology, as it is the case of the latest national seismic hazard map (IGN-UPM, 2013). In that work, fault information has been included as a complementary criterion for defining area sources and, particularly, the fault length, which is used to set the upper bound of the magnitude-frequency distribution in each one. In this map, faults have not been included as seismogenic sources in the calculations due to the incompleteness of the fault data across the country. However, that is not the case of SE Spain, the most seismically active area in the country (García-Mayordomo et al., 2007) and where one of the longest active fault systems is located: the Eastern Betics Shear Zone (EBSZ).
The EBSZ is a ∼400 km long transpressive tectonic corridor of mainly reverse and left-lateral strike-slip kinematics with a characteristic NE-SW sigmoidal trend (Bousquet, 1979;De Larouzière et al., 1988;Silva et al., 1993) ( Figure 1B). It absorbs most part of the NW-SE convergence from the Eurasian and Nubian plates in Iberia (4-6 mm/yr; DeMets et al., 2015). Although it is a moderate seismicity region with rates lower than other European regions (e.g., Italy or Greece), many destructive earthquakes with intensities up to IX-X (I EMS98 ) have affected this region since historical times (e.g., 1518 VIII-IX Vera, 1522 VIII-IX Almería, 1674 VIII Lorca, 1804 VIII-IX Dalías-Berja, 1829 IX-X Torrevieja; Silva et al., 2014) ( Figure 1B). More recently, moderate magnitude earthquakes have taken place (e.g., 1993-1994  Over the last two decades, the EBSZ tectonic context has motivated many geological, geodetic and paleoseismic studies, especially since the 2011 Lorca earthquake (e.g., Moreno, 2011;Ortuño et al., 2012;Echeverria et al., 2013;Insua-Arévalo et al., 2015;Martín-Banda et al., 2015;Ferrater, 2016;Martínez-Díaz et al., 2018;Masana et al., 2018). Accordingly, PSHA studies focused here have successfully incorporated faults as independent sources either considering a characteristic earthquake or Gutenberg-Richter model (García-Mayordomo, 2005;García-Mayordomo et al., 2007;Rivas-Medina et al., 2018). However, no PSHA study has yet modeled fault sources from the perspective of multi-fault ruptures, so the present constitutes the first one in this matter. The incorporation of this variable representing fault behaviors widely observed in nature is seen as a step toward more realistic approximations of fault modeling and seismic hazard in the EBSZ.
The aim of our paper is to explore a fault-based PSHA in a moderate seismicity region as it is the EBSZ in SE Spain, using a state-of-the-art fault source modeling approach that implements fault data as inputs, considers multi-fault rupture possibilities at the whole fault system scale and treats them as aleatory uncertainty. We use the source models computed and published in a previous article (Gómez-Novell et al., 2020). There, an emergent methodology (SHERIFS; Chartier et al., 2019) is applied to perform a source modeling of the EBSZ by computing magnitude-frequency distributions (MFDs) at the scale of the whole EBSZ fault system considering different scenarios of fault ruptures. The adequacy of the models is evaluated and ranked in terms of their agreement with the regional seismic catalog and paleoseismic data. For the present study, we adapt and re-compute these models, weight them in a logic tree and use them as inputs for the hazard calculations. We compute seismic hazard using the OpenQuake engine (Pagani et al., 2014), evaluate the influence of the selected fault source models on the seismic hazard and compare them with previous PSHA approaches in SE Spain. We also place the interest of our study for other low-to-moderate seismicity regions worldwide.

GEOLOGICAL SETTING
The EBSZ is a main active fault system in SE Spain. The faults affect Quaternary sedimentary deposits and generate geomorphological and morphotectonic features (e.g., fluvial channel offsets, fault scarps) along a left-lateral transpressive corridor, evidencing tectonic activity during the recent Quaternary (De Larouzière et al., 1988;Silva et al., 1993). The onset of the neotectonic activity under the nowadays compressional transpressive regime in the Eastern Betics occurred during the Late Miocene (e.g., Rodríguez-Fernández et al., 2012). From SW to NE these faults are ( Figure 2A (Bousquet, 1979;Silva et al., 1993). Characterizing the role and activity degree of these faults during the Quaternary has been the aim of many geomorphological, paleoseismic and geodetic studies in the last two decades, making the EBSZ one of the most studied fault systems in the Iberian Peninsula so far.
The EBSZ is one of the most seismically active regions of the Western Mediterranean and it has hosted many damaging earthquakes since the historical period. Most damaging earthquakes in the area have been attributed to the main faults forming the EBSZ. Important examples are the 1518 Vera earthquake, attributed to PF (e.g., García-Mayordomo, 2005), the 1522 Almería earthquake to CF (e.g., Reicherter and Hübscher, 2007), the 1674 Lorca earthquake to AMF (Martínez-Díaz et al., 2018), the 1829 Torrevieja earthquake to BSF (e.g., Alfaro et al., 2012) or, in recent times, the 2011 Lorca   (IGME, 2015b;García-Mayordomo, 2015), whose parameters are indicated in Table 3. A dashed purple line delineates the mapping area where the seismic hazard is computed. The location of the main faults of the EBSZ from Figure 2A within its buffer area (#55) are indicated.
Frontiers in Earth Science | www.frontiersin.org December 2020 | Volume 8 | Article 579398 4 earthquake to AMF (Martínez-Díaz et al., 2012a) ( Figure 1B). Some of these historical earthquakes are known to have ruptured the surface, as recently unveiled for the 1674 Lorca earthquake (Martínez-Díaz et al., 2018). Importantly, most of the instrumental seismicity in the EBSZ occurs relatively at shallow depths and focal mechanisms are consistent with the strike-slip and reverse character of the main faults (Stich et al., 2007, Stich et al., 2010. Paleoseismic trenching and detailed geomorphological studies have demonstrated the capability of the system to generate large magnitude earthquakes (Mw > 6.0) and allowed to constrain critical seismic parameters (slip rate) in several of the main faults: AMF, CF, CAF or LTF ( Table 1). The slip rate ranges between 1.0 and 1.7 mm/yr for AMF, according to estimates from Ferrater et al. 2017; 1.0 and 1.3 mm/yr for CF (Moreno, 2011); 0.37 ± 0.08 mm/yr for CAF (Martín-Banda et al., 2015) and 0.12-0.17 mm/yr for LTF (Insua-Arévalo et al., 2015). Conversely, PF, BSF and the northern termination of AMF have not been studied using detailed paleoseismic trenching analysis and their slip rate estimations are based mainly on long-term displacements of geological markers (e.g., Herrero-Barbero et al., 2020), which may be affected by larger uncertainties.
Continued tectonic activity in the EBSZ is also inferred from the current geodetic data. The obtained velocities of AMF-PF combined (1.5 ± 0.3; Echeverria et al., 2013) are consistent with the paleoseismic and geomorphologic slip rates of AMF, suggesting a main role of this fault into the deformation. Similarly, geodetic velocities of CF are in agreement with its geological slip rates (1.3 ± 0.2 mm/yr; Echeverria et al., 2015). In the NE end of the EBSZ GPS data indicates a partitioning of the deformation into lateral and reverse components, and compatibility with the kinematics of CAF and BSF, respectively (Borque et al., 2019).

Source Modeling
In this work, we consider two types of seismogenic sources that are combined: faults and area sources ( Figure 2). Fault sources are modeled within the EBSZ and area sources are modeled outside of it. The resulting composite source models are not a classical hybrid approach where faults accommodate earthquakes from a cutoff magnitude and area sources accommodate background seismicity (e.g., García-Mayordomo et al., 2007); instead, in our models both fault and area sources are modeled following a Gutenberg-Richter (GR) MFD for all the magnitude range (Mw ≥ 4.0), as we explain in the following sections.

Fault Sources
Fault sources are modeled exclusively for the EBSZ. The source modeling performed is based on the SHERIFS methodology as published in Gómez-Novell et al., 2020, where the first step of a source modeling for PSHA is performed at the EBSZ: the feasibility of multi-fault rupture scenarios is analyzed and earthquake occurrence on faults is modeled at the whole fault system scale following a GR distribution.
SHERIFS (Chartier et al., 2019) is a method of fault source modeling that uses geological fault data as inputs to calculate fault earthquake rates. Among its advantages, 1) it allows to model earthquake occurrence in fault systems as a whole and considering sets of multi-fault ruptures; 2) is well adapted for regions, such as the EBSZ, where geological data is the prime source of information for characterizing faults. In this sense, SHERIFS requires: i) a 3D geometry of the fault system, ii) slip rate of each fault source, iii) a target shape MFD of the fault system (e.g., GR) which is based on published studies or seismic databases and iv) a scaling relationship to estimate the upper bound of the MFD distribution (Mmax). The slip rate of each fault is treated as a budget and it is consumed following an iterative process to fit the target MFD, in which the different rupture scenarios are computed separately until the budget of the limiting faults is exhausted. The resulting MFDs are representative of the whole fault system and are weighted in a logic tree according to their consistency with data from the regional seismic catalog and paleoseismic studies (see Gómez-Novell et al., 2020 for details on how SHERIFS is applied to compute and weight the fault system MFDs).
Only the faults of the EBSZ have been considered for the SHERIFS modeling, because they are the ones absorbing most of the plate convergence motion in SE Spain (Bousquet, 1979;Silva et al., 1993) and because they show the best geological and paleoseismic data available in the region for seismogenic source modeling. The main fault sources are shown in Figure 2A. Faults outside of the EBSZ or secondary faults have not been considered mainly due to lack of enough quality data (e.g., Crevillente Fault, Torrevieja Fault, San Miguel de Salinas Fault, El Andarax Valley Fault, El Alquián Fault System, among others; Figure 2A).
For the fault source modeling we adopted the fault data and characteristics used in Gómez-Novell et al., 2020, obtained from available and mainly published data. Faults are sectioned ( Figure 2A) according to segmentation in literature based on their geometry, geomorphic expression, seismicity, kinematics and activity evidence. Offshore segmentation of CF is adopted from Moreno (2011) Table 2). Fault traces are simplified with respect to the ones depicted in the Quaternary Active Faults database of Iberia (QAFI v.3) (IGME, 2015a; García-Mayordomo et al., 2017) for the purposes of this study and multi-branch sections are simplified to a single one. This is consistent with studies that suggest that fault traces likely link at depth (e.g., CF: Moreno, 2011;AMF: Martínez-Díaz et al., 2012b). Dip, kinematics and seismogenic depths (Table 1) are extracted from the QAFI database, which compiles data from published studies on these faults. Fault slip rates are directly extracted from published works and, in selected cases, from unpublished data or expert opinion, as specified in   Figure 2A for the location of each fault section mentioned in the Gómez-Novell et al., 2020 for specific details on the criteria followed for slip rate selection, especially for unpublished work). The fault-system source model considered is spatially constrained to the limits of area source #55 of the ZESIS national scale area source model (García-Mayordomo, 2015), available at the ZESIS database (IGME, 2015b), which was delineated as the envelope of the surface projection of the EBSZ fault planes and characterized by the seismic record within (IGN-UPM, 2013). In this present approach, earthquake occurrence in this area (hereinafter called EBSZ buffer area) is completely modeled on faults from the geological fault data shown in Table 1. However, we also consider the alternative possibility that a small portion of the seismic moment release in the EBSZ is due to unknown faults (e.g., blind, unmapped) different from the considered ones, by modeling background seismicity in one of our models. Although the EBSZ buffer area is a narrow area constrained to the surface projection of well-known faults, it cannot be precluded that some activity may occur on unknown active faults not considered in this study (section 3.1.1.2).
Therefore, two alternative fault source models are envisaged within the EBSZ buffer area, one exclusively considering the faults' activity (FM model) and another one considering a fraction of background seismic activity overlapped with the faults (FM_bg model).

Fault Sources: FM Model
For this model, earthquake rates in the EBSZ buffer area are modeled only on the EBSZ fault system, assuming that all seismicity over Mw 4.0 occurring within its limits is related to these faults. This fault-system source model is based on the results of the modeling performed in Gómez-Novell et al., 2020 without any additional modifications. Four scenarios or hypotheses of fault and multi-fault ruptures involving from single fault sections to nearly the whole system (hypothesis 1 to 4; Table 2) are modeled using SHERIFS and the fault data in Table 1 to obtain different GR distributions for the whole EBSZ system. In order to incorporate the epistemic uncertainty in the input data, the GR of each scenario is composed by 20 GR curves resulting from 20 computations exploring random values of input parameters, i.e., fault slip rates, b-values and Mmax of GR distribution. Slip rates are explored within the ranges provided in Table 1 and the b-values within the 0.8-1.2 range (Figure 3), whose mean (1.0) is in agreement with the one of the EBSZ in the ZESIS model (García-Mayordomo, 2015). The Mmax is set following Wells and Coppersmith, 1994 for rupture area and all kinematics. The resulting earthquake occurrences of each scenario (80 GR curves for all four scenarios) are expressed as the seismicity rates of each fault-rupture involved in that scenario (see Gómez-Novell et al., 2020 for details on the criteria applied to define the fault rupture scenarios, the input parameters used for the modeling with SHERIFS and the calculations). Furthermore, in the previous paper, the agreement of the resulting MFDs with the earthquake distribution from the catalog and paleoseismic record is analyzed and a qualitative ponderation of each hypothesis for PSHA is proposed, which we use here to set our logic tree weighting scheme (section 3.3 and Figure 3). All these outputs are used to build the source models in hazard calculations.

Fault Sources Including Background Seismicity: FM_bg Model
We introduced, by means of a background ratio, a modification respect to the previous FM model in order to consider that part of the seismic moment release in the EBSZ could be due to unknown active faults apart from the main ones. The rest of the input parameters are the same as mentioned for the FM model.
The background ratio is a parameter from SHERIFS that allows sharing a fraction (percentage) of the seismicity rate of a given magnitude between the faults and the background or the area directly surrounding the faults. For instance, a background ratio of 70% for Mw 4.0 means that 70% of the total seismicity rate of Mw 4.0 is modeled on faults and 30% as background seismicity. Therefore, the occurrence rates of the selected magnitude range on the faults are lowered with respect to the FM model consistently with the ratio adopted, avoiding double-counting of seismicity and respecting the GR distribution (see Chartier et al., 2019 for details on the 2 | Fault sections involved in the maximum fault and multi-fault ruptures envisaged in each hypothesis of the fault source modeling. See Figure 2A for the location of the mentioned fault sections. The maximum magnitude ranges of each hypothesis of the FM and FM_bg models are indicated. These ranges result from the random sampling of different fault ruptures allowed within a same hypothesis by the SHERIFS methodology.

Hypothesis/fault rupture scenario
Maximum expected multi-fault ruptures Maximum magnitude range (Mw) 1 Only single fault section ruptures are allowed 6.6-7.
December 2020 | Volume 8 | Article 579398 background ratio). In our case, the background area corresponds to the EBSZ buffer area, whose seismic potential is then not only defined by the faults contained but also by the background seismicity as an area source. Chartier et al., 2019 recognize that the user may rely on expert opinion to set this ratio due to its difficulty but propose to use the distance between earthquake epicenters and faults as possible criteria to set it. We cannot follow this approach due to large location errors at the EBSZ: 67% of the Mw ≥ 4.0 earthquakes are previous to year 1962 (IGN-UPM, 2013) and hence show epicentral errors larger than 10 km or even 20 km (García-Mayordomo, 2005). We then set the background ratios following expert opinion but based on criteria from seismicity data of the region as we explain below. These ratios are 70%, 80%, 90%, 95% for Mw 4.0, 4.5, 5.0, and 5.5, respectively, which represent the percentage of the total EBSZ seismic moment that is actually modeled on faults for each 0.5 magnitude bin. The remaining is modeled as background in the area of about 10 km around the faults.
The selected percentages are consistent with the fact that most of the instrumental earthquakes of Mw ≥ 4.0 in the EBSZ have occurred in the modeled faults. This is supported by the location of recent moderate earthquakes (e.g., 2011 Lorca seismic series related to AMF; Martínez-Díaz et al., 2012b) and by the fact that focal mechanisms of seismicity are in agreement with the main faults of the region (Stich et al., 2007, Stich et al., 2010. On the other hand, the GR distributions modeled with SHERIFS for the EBSZ in Gómez-Novell et al., 2020 predict good fits with the regional seismic catalog for most of the rupture scenarios in the Mw < 6.0 magnitude range. This supports that known faults might be responsible of most of the seismicity within the EBSZ. Consequently, we select ratios that do not modify significantly the GR distributions, hence that do not invalidate the proposed logic tree weighting in that study and comparison can be directly done with the FM model. Only events in the range of Mw 4.0-5.5 are affected by the background ratio in our study, because earthquakes this size rarely cause surface rupture (Biasi and Weldon, 2006) and thus, they may occur in unrecognized faults and unadvertised by the paleoseismic studies. Mw > 6.0 earthquakes though, are considered to be produced in larger faults with geomorphic expression so they are modeled exclusively on the known EBSZ faults.

Area Sources
Outside the bounds of the EBSZ buffer area we consider the area sources of ZESIS model (García-Mayordomo, 2015). The area sources considered here are ( Figure 2B): #29, #30, #31, #34, #35, #36, #37, #38, #39, #40, #41, #43, #44, #45 and, for the FM_bg model only, #55. These areas are all located within a radius of 50 km from the external boundary of the mapping area in order to account for the influence of seismogenic sources adjacent to the EBSZ buffer area (#55). The potential influence of farther sources is considered negligible for the purposes of this study. The earthquake occurrence in these areas is modeled according to GR distributions obtained in the frame of the project to update the national seismic hazard map of Spain (IGN-UPM, 2013), and their parameters (a and b-value, Mmin, Mmax; Table 3) are adopted here from the ZESIS database (IGME, 2015b). In the case of area #55, earthquake occurrence is controlled by the seismic moment assigned by the background ratios in each source model and following the target MFD of the faults (section 3.1.1.2). The area sources are added to the fault sources modeled with SHERIFS to build the different source models for the hazard calculations.
The characterization of the area sources is further completed defining a distribution of nodal planes and hypocentral depths in each area (OpenQuake manual v.3.7.1; GEM, 2019). The former controls the strike, dip and rake of earthquake ruptures inside the FIGURE 3 | Logic tree schemes followed in this study for both FM and FM_bg source models. The weights (over 1) of each branch in the logic tree are indicated in brackets. Hyp.: hypothesis; WC94: Wells and Coppersmith, 1994;CB14: Campbell and Bozorgnia, 2014. Frontiers in Earth Science | www.frontiersin.org December 2020 | Volume 8 | Article 579398   Table 3). Eventually, a set of different nodal planes is defined for each area source representing each one the same probability of occurrence. Similarly, two equally probable hypocentral depths are defined for each source area ( Table 3): the upper is set at 5 km and the lower between 8 and 10 km depending on the minimum common seismogenic depth of all faults contained in the area. The maximum seismogenic depth of the faults in the QAFI database defines the seismogenic depth for each area source.

Final Composite Seismogenic Source Models for Hazard Calculations
Two resulting composite source models are used for the calculations: FM + area sources and FM_bg + area sources. The "composite" term is used here to remark that each source model is composed by a fault source model and an area source model ( Figure 3). The former sets the occurrence rates of all the fault ruptures envisaged in each fault rupture hypothesis modeled with SHERIFS within the EBSZ buffer area. The latter sets the occurrence rates of the different area sources considered outside the EBSZ buffer area and the background seismicity in the case of FM_bg model. The final models are referred by the distinctive fault source model name henceforth (i.e., FM and FM_bg). Because each of the four multi-fault rupture scenarios is composed of 20 different GR curves, a total of 80 input source Frontiers in Earth Science | www.frontiersin.org December 2020 | Volume 8 | Article 579398 models for OpenQuake compose both FM and FM_bg source models ( Figure 3). The fault source part (e.g., earthquake rates, Mmax) is different for each 80 models because it depends on a) on the rupture scenario considered and b) the values of the input parameters explored randomly during the SHERIFS calculation. Each FM and FM_bg are computed separately in a different hazard calculation in order to better test and highlight the influence of including a background ratio, and to compare results more directly.

Ground Motion Prediction Equation
We selected the ground motion prediction equation ( We acknowledge that considering a set of GMPEs is a normal practice in the attempt to capture the epistemic uncertainty related to ground motion. However, we consider just one GMPE to facilitate the detection of the impact of fault incorporation into PSHA and to be able to perform direct comparisons with

Logic Trees Explored
Each source model is computed separately. Thus, two logic trees are explored in this study as depicted in Figure 3. For each logic tree, the four different hypotheses of multi-fault ruptures are considered in order to incorporate the epistemic uncertainty of earthquake occurrence on faults. Each rupture hypothesis constitutes a branch of the logic tree, weighted based on the analysis of the agreement with the seismic catalog and paleoseismic data performed in Gómez-Novell et al., 2020.
That study concluded that hypotheses 2 and 3 were the ones with better agreement with the regional seismic catalog and paleoseismic data, followed by hypothesis 1 and lastly, hypothesis. 4. A better agreement implies a more suitable hypothesis for the study region. Accordingly, we assigned a weight of 35% to both hyp. 2 and 3, 20% to hyp. 1 and 10% to hyp. 4. Each hypothesis is composed by 20 GR curves, hence 20 input source models for OpenQuake that include a fault source model and an area sources model. The weight of each of the 20 branches or models that form a hypothesis equals 1/20 of the weight of the whole branch ( Figure 3).

Hazard Calculations
Hazard calculations are performed using the OpenQuake engine (Pagani et al., 2014). We compute a set of mean hazard maps in terms of PGA as function of probabilities of exceedance of 10, 5, 2 and 1% in a time window of 50 years, which correspond to return periods (RPs) of 475-, 975-, 2,475-and 4,975-years, respectively (Figures 4, 5). The region considered is depicted in Figure 2B and follows a 2 × 2 km grid, accounting for a total of 10,556 calculation points for each return period. This region is defined to comprise completely the EBSZ buffer area and include important cities of the region (Murcia, Lorca, Almería, Vera, Torrevieja, Alicante and Cartagena). To analyze fault influence, in the cities adjacent to the EBSZ faults (Murcia, Lorca, Almería, Vera, Torrevieja) we produce mean hazard curves in terms of PGA ( Figure 6) and hazard disaggregation for distance-magnitude (Figure 7 and Supplementary material). Soil conditions are always considered to be rock or hard soil (Vs30 760 m/s). The maximum distance considered for the hazard calculations is 100 km in order to better reproduce the near-fault related hazard of the EBSZ and because farther sources have negligible impact in the study area (see disaggregation results in section 4.1). Source input files of the hazard calculations including source models, GMPE, logic tree and calculation specifications are included in the Supplementary material of the article.

RESULTS
The results of the calculations are presented in mean hazard maps (Figures 4, 5) and mean hazard curves for five cities of the region (Figures 6). Maps were obtained after interpolation using the natural neighbor approach from the OpenQuake outputs. A common colored scale has been set for all maps to facilitate comparison. Mean hazard curves for the cities of Murcia, Lorca, Almería, Vera, Torrevieja are also depicted in terms of PGA on rock ( Figure 6). To test the influence of the different fault sources and the seismic background on the hazard, we computed hazard disaggregation of the FM_bg model for the magnitude-distance relationship. In Figure 7, we depict representative disaggregation results in three cities and for two RPs (475 and 2,475 years). For complete disaggregation, figures of all the cities considered in Figure 6 and all four RPs are available the Supplementary material.
The results show a clear control of the EBSZ fault over the seismic hazard of SE Spain for all RPs, especially in the vicinity of faults (10-20 km). Slip rate values directly affect fault influence into the hazard because highest accelerations are located around faults with higher slip rates of the system, while lower slip rate faults have much less influence. The impact of the background into the models is very limited given the similarity of results from both FM and FM_bg.

Fault Influence
The mean hazard maps from both FM and FM_bg models (Figures 4, 5) produce very similar results in terms of absolute hazard levels and distribution and increments between them are very low (Figure 8). The modeled EBSZ fault system exerts a clear control into the seismic hazard of SE Spain for all RPs, especially In both models, the maximum acceleration values range from 0.5-0.6 g for the 475 years RP to 1.2-1.3 g for the 4,975 years RP (Figures 4, 5) and are located along the faults with higher slip rates. That is the case of the Southwestern termination of AMF, near Lorca (AMF-1 and 2) or CF near Almería, which have the highest slip rates of the system (>1.0 mm/yr; Table 1 and Figure 2A). On the other hand, lower slip rate faults have minor influence in the hazard levels for all RPs, also depending on their slip rates. Mean hazard maps from both FM and FM_bg models (Figures 4, 5, respectively) show that PF is the fault with less footprint on the hazard for all RPs, which is consistent with its slip rate being the lowest of all the system according to the available estimations (<0.1 mm/yr; Table 1).
The hazard disaggregation of the FM_bg model (Figure 7) shows that all the contributions to the seismic hazard in a specific site are from sources located within a 20-30 km radius, suggesting that farther sources do not influence the hazard. Higher slip rate faults contribute the most to the seismic hazard values of the closer zones, while lower slip rate faults have very small contributions, even in the sites in their closest vicinity. In the cities of Lorca and Almería ( Figures 7A-D), the magnitudedistance relationships that influence most the hazard coincide with the location of AMF-1 and 2, and CF-3, respectively, which are the fault sections with higher slip rates in the system ( Table 1). Lorca is located <1 km from AMF in the limit between AMF-1 and 2, which explains the maximum influence from sources located between 1-5 km. These maximum AMF-1 and 2 contributions are in a tight Mw 6.5-6.75 range for both RPs ( Figures 7A, B), which might suggest the controlling earthquakes in this city. Almería is located at ∼15 km from CF-3, consistent with the distance of the source's maximum contribution. The higher distance of the fault from the city causes its contribution to be considerably lower than AMF-1 and hence other non-fault sources influence the hazard as well ( Figures 7C, D). The contributions of CF-3 comprise a wider magnitude range (Mw 5.0-7.0) in the 475 years RP ( Figure 7C), while they concentrate around Mw 6.5-7.0 in the 2,475 years RP ( Figure 7D). Conversely, in zones such as Vera, located next to very low slip rate faults (i.e. PF), the main contributions to the seismic hazard are not from the closer faults, but from other farther faults ( Figures 7E, F). In this city AMF-1 is the main contribution despite this fault being ∼25 km farther, while PF is the least even for long RPs ( Figure 7F).
A hanging wall effect can be observed in non-vertical dipping faults such as AMF, CAF or BSF (Figures 4, 5). Such near-fault effect is more visible at longer RPs (i.e., 2,475 and 4,975 years), where PGAs for two selected sites at the same distance of a fault trace but on opposite sides, are higher on the hanging wall. A clear example is the transect perpendicular to AMF near the city of Lorca ( Figure 9E). This is a significant observation for seismic hazard assessment in near-fault zones as we discuss in section 5.1.1.

Background Influence
The similarities between the mean hazard maps from the FM and FM_bg models suggest that adding a seismic activity background ratio in the FG_bg model produces minor changes to the final hazard results within the EBSZ buffer for any RP (Figure 8). In general, the background contributes to raise the hazard values with respect to the FM model, while the decreases are depreciable. PGA increments, especially for  Figure 9D.
Frontiers in Earth Science | www.frontiersin.org December 2020 | Volume 8 | Article 579398 the shorter RPs, are mild ( Figures 8A-C), but for the 4,975 years RP the hazard values are incremented up to a maximum of ∼0.1 g in some zones (e.g., E of BSF; Figure 8D). For the faults, PGA values remain practically identical for all RPs, which means that earthquake occurrence on most faults is higher than the one generated by the background. However, the background contributes to slightly increment the hazard within the EBSZ buffer area in the zones farther from the traces and influence of the faults (e.g., E of BSF-4) or next to faults with slip rates much lower than the average of the system (∼0.5 mm/yr; based on the mean values in Table 1); mainly PF. In these zones, the poor influence of faults allows the background of the FM_bg model to overcome the occurrence rates (i.e., probabilities of exceedance) of the hazard from the FM model, thus higher values are found in the FM_bg model within a same RP. The maximum hazard increments caused by the background take place E of BSF-4 (Figure 8), where the Bajo Segura faultsystem terminates. This is because fault-related seismic hazard influence is larger and wider in the transversal sections across faults rather than in their longitudinal ends, as it can be seen in CF or the Southwestern termination of AMF (Figures 4, 5). In addition, BSF-4 has a lower slip rate (∼0.2 mm/yr; Table 1) than the average of the faults in the system and the closest area sources in this zone, i.e., #31, #39 and especially #41 ( Figure 2B), have very low activity rates (IGME, 2015b). This particular situation causes low occurrence rates of the hazard from the FM model in this zone and allows the background (FM_bg) to overcome such rates more remarkably than any other zone.
These observations are consistent with the disaggregation results of the FM_bg model as in general the background shows low contributions on the hazard, especially where faults dominate the hazard (Figure 7 and Supplementary material). Only in the case of Vera, the seismic background has a significant contribution due to the lower slip rates, hence activity rates, of the PF ( Figures 7E,F). This causes the background to overcome the contribution of PF in the Mw 4.0-5.5 range and to contribute significantly to the seismic hazard levels of the city, especially in longer RPs ( Figure 7F).

Area Sources Influence
From the mean hazard maps of both FM and FM_bg models it can be observed that the area sources have a much minor footprint in the seismic hazard distribution of the EBSZ than the faults. This is mainly because the activity rates of the area sources modeled are either very low (#31, #36, #39 or #41; IGME, 2015b) or are located too far (>100 km) to have a significant impact in the hazard of the mapped region, such as the North-African source areas (#43, #44 and #45) ( Figure 2B). Nevertheless, area sources such as #37 and #38 have a recognizable influence on the hazard results as their shapes can be recognized in the hazard maps especially for the 2,475 and 4,975 years RPs (Figures 4, 5). These area sources have high activity rates comparable to the EBSZ buffer area (IGME, 2015b) and contain faults known to be capable of producing Mw > 5.0 earthquakes, such as the Crevillente Fault or the Andarax Valley and El Alquián fault systems (#37 and #38, respectively; Figure 2).
The hazard disaggregation shows the contribution of area source #38 on the hazard of Almería for RPs of 2,475-years RP ( Figure 7D) as the city is located inside this area source. Contrarily, neither Lorca or Vera have contribution of area sources into their hazard ( Figures 7A-F, respectively) as they are located ∼8-9 km from the closest area source limit. Although these areas are at relatively close distances, the higher proximity of AMF-1 in Lorca and the low activity rates of the closer area #36 in Vera prevent them to contribute to the hazard in these cities. Figure 6 shows the mean hazard curves obtained for five cities located close to the EBSZ faults (Almería, Lorca, Murcia, Vera and Torrevieja) for both the FM and FM_bg models. It is clear that setting a background ratio for the EBSZ buffer area does not imply dramatic changes in the hazard curves of these cities, which are all very close to active faults. In the case of the Almería, Lorca and Murcia, the respective hazard curves of both models predict practically identical PGA values for any RP. Conversely, in Torrevieja and, especially, in Vera, the FM_bg model predicts higher PGAs than the FM model, although these differences are very small: less than <0.1 g for the ∼10,000 years RP. The reason for these differences between models is the lower slip rates of the BSF and PF, respectively for Torrevieja and Vera, compared to the other faults close to the other studied cities.

Comparison With Previous Probabilistic Seismic Hazard Assessment in SE Spain: Implications of the Modeling
The results presented in this study highlight a clear influence of the faults of the EBSZ in the seismic hazard values of SE Spain, as shown in Figures 4-7. Both of the FM and FM_bg source models applied for the hazard calculations show very similar results, meaning that the application of a background ratio has minor influence on the hazard, especially in those zones closer to faults (Figure 8). Hence, for the comparison analysis that follows we will focus on results of the FM model, though practically all the observations are also applicable to the FM_bg model results. In the same way and due to the low impact of the background on the hazard, the fault influence observations from the disaggregation of the FM_bg model (Figure 7) are valid for the FM model. We compare our results with a PSHA produced only considering the ZESIS area sources model, as well as with a recent fault-based PSHA done in SE Spain (Rivas-Medina et al., 2018).

Comparison With a Probabilistic Seismic Hazard Assessment Using the ZESIS Area Source Model
The PSHA results of our FM model are compared with a PSHA reached here by using exclusively the area sources of the ZESIS model ( Figure 2B) and the same GMPE as our calculations (Campbell and Bozorgnia, 2014 are found between both models within the limits of the EBSZ, where the faults are modeled ( Figures 9A-D); outside these limits, the hazard values do not experience any change since all use the same area sources model. The FM model increases the PGA with respect to the areas approach in the zones strictly around the faults. Hence, faults concentrate the hazard highlighting the controlling effect of the faults on the hazard levels. The curves in Figure 10A depict this effect as Lorca, which is located nearly on top of AMF (Figure 2A), shows the highest increment in the FM model. Next, Murcia and Almería show progressively lower PGA increments with respect to the areas model as the distance to the fault increases. These hazard increments in the near-fault zones are achieved by reducing the hazard in the zones farther from the faults with respect to the areas approach and are more visible in longer RPs (975 years RP or larger; Figures 9B-D), where the faults have more influence.
The influence of faults on the hazard (net PGA and areal extension) respect to the areas approach, strongly depend on the fault slip rates. The cross sections A-A′ and B-B' (Figure 9E), show that higher slip rate faults in the system such as CF or AMF-1, respectively (1-2 mm/yr; Table 1) influence the hazard up to 10-20 km far from the fault trace for all return periods. This causes significant increments in the PGA values within 20-40 km-wide influence zones of ∼70-80% for all RPs in CF-3 and ∼130-150% for all RPs in AMF-1 ( Figure 9E). Such increments are slip rate dependent, thus higher for AMF-1 because it has the highest slip rate. Despite this, the disaggregation results of the FM_bg model show that faults can control the hazard in regions farther than the influence zone of the faults up to 20-30 km from its traces, as it is the case of AMF-1 in Vera ( Figures 7E,F). Farther than this distance limit however, no source has significant influence in any part of the studied area. Fault influence zones are also dependent on the dip direction in non-vertical dipping faults (AMF, BSF), wider in the hanging wall ( Figure 9E), implying a hanging wall effect.
Contrarily, lower slip rate faults have much less impact on the hazard with respect to the areas approach. Although no exhaustive slip rate analysis has been performed, it can be observed that the faults that have slip rates <0.1 mm/yr do not increase the seismic hazard and it is overcome by the areas approach, as it is the clear case of PF (Figures 9, 10A in Vera). Even so, faults with slightly higher slip rates (∼0.1-0.4 mm/yr; AMF-4, LTF or BSF) do influence the hazard, but their impact is more limited than AMF-1, 2 or CF, even for long RPs ( Figure 9D). The cross section in BSF-3 (C-C'; Figure 9E) depicts one of these latter examples. This fault section has a slip rate range 4-14 times smaller than AMF-1 and CF (0.12-0.3 mm/yr; Table 1), which causes its influence zone to be over 2-3 times smaller (<10 km wide). In this case the near- fault PGA increments are negligible for short RPs of 475-and 975-years and very small (<0.1 g) for longer RPs ( Figure 9E). The small influence of the fault causes the FM's PGA values to be overcome by the ones of the areas approach in less than 5 km from the peak of maximum fault influence. This is also depicted by the Torrevieja curves ( Figure 10A) and the maps of Figure 9. The implications of lower slip rate faults on the hazard are discussed in section 5.2. The consideration of near-fault hazard in zones in the vicinity of the EBSZ faults as both FM and FM_bg models do, including effects such as the hanging wall, provides more realistic estimations of the seismic hazard. This is relevant because during earthquakes, zones close to active faults are exposed to higher accelerations than farther zones due to these near-source effects (e.g., 2011 Lorca earthquake; Alguacil et al., 2014). Their location respect to the dip direction of the fault is also relevant for the ground motion, usually higher in the hanging wall (e.g., Yu and Gao, 2001). Conversely, area source approaches do not reproduce these effects, but uniform PGA distributions within an area, as depicted by the nearly identical mean hazard curves of the cities within the EBSZ buffer area ( Figure 10A). A detailed and representative analysis of the relationship between PGA, distance to the fault and influence of fault slip rate is a necessary step to be taken in further hazard studies in the region, as it could be useful to define amplification coefficients in earthquake resistant provisions to account for near-fault effects at sites located next to active faults.

Comparison With Rivas-Medina et al., 2018 Fault-Based Approach
We compare the results from the FM model with the fault-based study of Rivas-Medina, 2014, recently published in Rivas-Medina et al., 2018. In that study a PSHA of Southern Spain is performed considering faults as independent, non-interacting sources, following a hybrid methodology (HM model) that avoids setting an arbitrary cutoff magnitude for modeling seismicity in the background (see details in Rivas-Medina et al., 2018). The PSHA was performed for a larger territory of SE Spain than the present study and for all the active faults of the region as compiled in the old QAFI v.2 database back in 2012. The areal model for this study (GM10) is different from the ZESIS model (aka GM12; IGN-UPM, 2013) as it does not identify the EBSZ buffer area (area #55; Figure 2B). The GMPE used is the same as in our study. We compare the mean hazard curves of the FM model with the HM model, including Alicante ( Figures 10B, 10C). The hazard curves for the cities of Lorca, Vera and Torrevieja were not computed specifically by the authors, hence we inferred visually the PGA values at these cities from pixel screening of the hazard maps for each of the four RPs calculated (475, 975, 2,475 and 4,975 years; Figure 10C).
The main contribution of our work with respect to Rivas-Medina et al., 2018 is the consideration of multi-fault ruptures in the source modeling, which causes differences in the hazard estimations mainly for long RPs. Regarding the EBSZ as a whole in PSHA is a novelty in Spain and might constitute a more realistic approach to model the seismic hazard. The main point is that multi-fault ruptures are regarded as a frequent phenomenon, demonstrated by numerous morphogenetic earthquakes in the last decades (e.g., Choi et al., 2012;Hamling et al., 2017;Villani et al., 2018).
In general, the study by Rivas-Medina et al., 2018 provides hazard levels in the range of the ones obtained by the FM model, however there are noteworthy differences. These are mainly due to 1) the consideration of multi-fault ruptures in the FM model, 2) differences in the fault data used in each source modeling, 3) the background seismicity of the HM model and 4) the number of faults considered. The FM model treats the fault system as a whole by allowing the possibility for faults to rupture simultaneously in multi-fault ruptures including large rupturing scenarios. This consideration allows a larger number of fault rupture possibilities to host a certain magnitude and therefore hazard values are generally higher for this model, as shown by the curves of Lorca, Almería and Alicante (Figures 10B,C). These differences are more notable in RPs over ∼1,000 years because multi-fault ruptures produce higher magnitudes that control the longer RPs. Moreover, the discrepancies in the slip rate values and segmentation between both FM and HM models are also responsible for the hazard differences. For instance, the higher values of the FM in Lorca are because updated slip rate data for AMF-1 and AMF-2 are used (1.65 and 1.0 mm/yr, respectively; Table 1), contrary to the HM model (0.3 and 0.5 mm/yr). The higher values in Almería, on the other hand, are because CF is divided in three segments with a mean slip rate of 1.1 mm/yr (Table 1) in the FM model, while in the HM it is divided in two segments, one of which uses an outdated slip rate (0.04 mm/yr).
In the cities of Murcia, Vera and Torrevieja, the curves of the HM model predict higher values than the FM model ( Figures  10B,C). In Murcia this is only for RPs over ∼500 years, as for lower RPs the tendency is the same as the mentioned for Almería and Lorca. We attribute this to the fact that CAF, which controls the hazard of this city (see disaggregation in the Supplementary material), is considered as a single-segment fault in the HM model (CAF; 0.54 mm/yr), while in the FM as two segments with different and lower mean slip rates (0.37 and 0.53 mm/yr for CAF-1 and 2, respectively; Table 1). In RPs over 500 years this implies that probably a big part of the magnitudes that in the HM model are achieved independently by CAF, in the FM model are achieved by multi-fault ruptures involving CAF-1 and 2. As the slip rates of CAF-1 are lower than CAF-2, the occurrence rates of these multi-fault ruptures are limited by the former (see slip rate treatment in SHERIFS on section 3.1.1) and hence are lower than the HM model in this area. On the other hand, the fact that the fault is divided in two segments increases the earthquake rates of lower magnitudes and consequently the PGA values for low RPs (<500 years) in this city with respect to the HM model ( Figure 10B). As such, linked ruptures between segments, i.e., higher magnitudes, are less frequent in low RPs.
In Vera and especially Torrevieja the higher values of the HM model are visible in all four RPs represented ( Figure 10C). In the former, the hazard is not controlled by the nearest faults due to the low activity rates of PF, so the background seismicity of the HM model overcomes the ones from the FM. Remarkably, performing this comparison with the curves from the FM_bg model would produce practically the same results, as the background seismicity of the HM model is modeled for a larger magnitude range (Mw 4.0-6.9) and implies higher seismicity rates. In the latter, the higher HM values are because in this area, additional fault sources apart from the BSF are considered in the modeling (Torrevieja and San Miguel de Salinas faults; Figure 2A). While the incorporation of more fault sources on the hazard might suggest that the results from the HM model in the Torrevieja zone are more accurate than the FM, it should be noted that including more faults on the hazard does not necessarily imply more accurate results. These additional faults are not included in our models because are secondary faults to the BSF system (Alfaro et al., 2012) and hence are not part of the main structures in the EBSZ. Moreover, neither their relationship with the EBSZ system and their geological parameters have been studied in detail. The inclusion of such faults on the hazard should account for these uncertainties as they could imply misguided hazard estimations.

Implications for Moderate Magnitude Earthquakes
The hazard results from the disaggregation show that moderate magnitudes in the Mw 4.5-6.0 range, although not being the most remarkable, have a non-depreciable contribution to the final hazard in the selected studied sites even in long RPs (e.g., Almería and Vera; Figures 7D,F, respectively). Moreover, according to Gaspar-Escribano et al., 2008, in SE Spain, lowto-moderate   , being the 2011 Mw 5.2 Lorca earthquake the latest one (Martínez-Díaz et al., 2012a).
The near-source effects of the 2011 Mw 5.2 Lorca earthquake produced exceptionally high accelerations around 0.37 g in the city of Lorca (Cabañas et al., 2011) due to: 1) the proximity of the epicenter to the city, 2) the shallow depth of the earthquake and 3) the directivity of the rupture propagation toward Lorca (Alguacil et al., 2014). The hazard map of the current Spanish Building Code (NCSE-02, 2002) for the 500 years RP and the updated national seismic hazard map for the 475 years RP (IGN-UPM, 2013) predict PGA values in Lorca (0.12 and 0.19 g, respectively) that underestimate the PGA reached in the Lorca earthquake. Conversely, the results of the FM model, but also FM_bg, predict values that differ less than 0.1 g from the ones reached in the 2011 Lorca earthquake: ∼0.45 g for the 475 years RP. Rivas-Medina et al., 2018 predict slightly lower values of 0.3 g in Lorca for the same RP, but that are also close to the measured ones.
The fact that fault-based PSHAs in SE Spain predict PGA values for the 475 years RP coherent to the ones of the Lorca earthquake, whose moderate magnitude produced severe damage, may suggest that fault-based approaches perform more realistic hazard predictions than classical area source approaches, particularly for sites located on top of faults or in their vicinity. The PSHA such as the one presented in this study or the one from Rivas-Medina et al., 2018 approximate better the hazard for short RPs due to the effect of faults. However, other site effects such as topography, directivity and changes in soil properties should be incorporated for a more accurate hazard estimation, especially in seismic microzonation.

Implications of the Fault Slip Rates
Slip rate data of faults is a controlling factor in the hazard values obtained in this study as it has a direct impact on the probabilities of exceedance in a given observational period. Higher slip rate faults (e.g., AMF-1 and 2 or CF) imply higher PGA values and influence than lower slip rate faults, some of which have slip rate values so low (<0.1 mm/yr; PF) that are not able to influence significantly the hazard at any RP. This is also observable in the slip rate values used by the source models that compose each branch of the logic tree in relationship with the resulting PGA for each model (Figure 11). Each source model is built by picking a random slip rate value from the fault uncertainty ranges provided as inputs (Table 1). Thus, aside from the fault rupture hypotheses considered, the models that use higher slip rate values within the uncertainty range show an overall tendency to produce higher PGA estimations at closer sites within the same hypothesis (e.g., AMF-1 in Lorca; Figures 11A,B). Nevertheless, the data in the plots from Figure 11 show high dispersion, which suggests that the slip rate uncertainty is not the only parameter affecting seismic hazard variability. In fact, the proximity of the analyzed site to the fault, the rupture hypothesis in question or the RP might control the slope of the mentioned tendency for each hypothesis. Closer sites, i.e., Lorca, result in steeper slopes, similarly to higher RPs ( Figures 11A,B).
PGA variations as a function of slip rate uncertainty are more remarkable in higher slip rate faults ( Figures 11A,B), despite their uncertainties being smaller than for lower slip rate faults ( Figures 11C,D). AMF-1 has an uncertainty that represents only the ∼3% of the mean, while it is ∼100% for PF-1 (Table 1), however the hazard variations linked to the slip rate variations are considerably higher for the former in Lorca than the latter in Vera. This is because lower slip rates have less impact on the return periods studied than higher slip rate faults, which is consistent with the low influence of PF in our hazard calculations ( Figures 7E,F).
The identification of the method and source of information upon which slip rate data is based on for each individual fault ( Table 1) is a crucial step, as it affects the reliability of such data and ultimately the hazard results. At the EBSZ, this is one of the main challenges today as most paleoseismic studies focus only on specific fault sections or single strands of faults with multiple branches, while others remain poorly studied. The faults that concentrate most part of these studies and thus have a better characterization are AMF-1 and 2, CAF and CF. Although most of these focus only on the major fault branches, their slip rate estimations are robust and consistent with overall geodetic data (see section 2 and references of Table 1). Conversely, PF, AMF-3 and 4 or BSF have not been object of detailed paleoseismic research to date and fault data is inferred based on qualitative geomorphological indicators and from geodetic data. These less studied faults at the EBSZ are the ones that have lower slip rate estimations ( Table 1) and contribute poorly to the hazard, which could imply a bias of such hazard towards well-studied faults.
Frontiers in Earth Science | www.frontiersin.org December 2020 | Volume 8 | Article 579398 Considering the clear sensitivity of the seismic hazard to the fault slip rates, the acquisition of new data in the less studied faults is critical to increase the quality and reliability of the current available data compared to other faults in the EBSZ. Updated parameters can have important implications in the hazard results around them, especially PF. On the other hand, in the zones farther from the mentioned faults and controlled by faults with more paleoseismic research, seismic hazard estimations are more reliable and critical seismic parameters are well constrained (see references in Table 1).

Implications of the Ground Motion Prediction Equation Selection
The selection of GMPEs for PSHA is also a difficult task as it has a great influence on the hazard results obtained for a particular study region and sometimes can mask or exaggerate the influence of a seismic source. In this study, we selected a NGA GMPE that is well adapted for fault-based PSHA (Campbell and Bozorgnia, 2014; section 3.2) as it is sensitive to near-fault effects. Therefore, the strong contribution of the faults into our hazard models could be seen as a result of applying such equation. In order to test the sensitivity of the hazard models to the GMPEs we performed the hazard calculation of the FM model using the equation by Akkar et al., 2014 (henceforth Akk14). We used this equation because is not as specifically designed for near-fault effects as the one by Campbell and Bozorgnia, 2014 (henceforth CB14) and the source characteristics of earthquakes are exclusively represented by magnitude and style of faulting. The mean hazard maps of the FM model using Akk14 are shown in Figures 12A,B for the 475and 4,975-years RPs. The differences with the CB14 equation are depicted in the maps of Figures 12C and 12D and the cross sections of Figure 12E for the same RPs. Only these two extreme RPs are shown in order to point out the maximum variations with the CB14 equation.
The sensitivity test performed with the Akk14 equation evidences that the influence of the EBSZ faults on the hazard results performed in this paper is independent of the GMPE used. Therefore, the marked influence of the faults into the PSHA is not a consequence of the GMPE. The resulting mean hazard maps of Akk14 ( Figures 12A,B) are very similar to the ones obtained with the CB14 equation (Figure 4) showing a similar hazard distribution. Faults have a clear influence on the hazard for both RPs explored in the near-fault zones and the highest PGA values are found along the traces of the faults, especially the most active according to their slip rates (e.g., AMF-1 and CF). Hanging wall effects are also recognized in the map using Akk14 for non-vertical dipping faults (AMF, CAF and BSF; Figures 12A  and 12B).
Despite this, due to its less sensitivity to near-fault hazard estimations, Akk14 computes a more diffuse hazard estimation around the faults, with lower PGA values in general and lower FIGURE 11 | (A) Relationship between the slip rate values of AMF-1 used in each 80 source models that compose the logic tree of this study and the respective PGA (g) values obtained in the city of Lorca for the return periods of 475 years and; (B) 2,475 years; (C) Same relationship for PF-1 and the PGA (g) obtained in the city of Vera for the return periods of 475 years and; (D) 2,475 years. Note that the faults considered are the closest to the cities analyzed and so that can highlight a better effect. The source models belonging to each fault rupture hypotheses that compose the logic tree of this study are depicted in different colors.
Frontiers in Earth Science | www.frontiersin.org December 2020 | Volume 8 | Article 579398 PGA increments as a function of distance to the fault trace ( Figure 12E). Therefore, CB14 produces higher PGA values than Akk14 for most of the mapped area and especially around faults, accentuating near-fault effects on the hazard levels. Maximum variations are located close to the traces of faults with higher slip rates, importantly around CF or AMF-1, where these are up to 0.1 and 0.2 g, respectively for the 4,975 years RP and <0.1 g for the 475 years RP (A-A′ and B-B′, respectively; Figure 12E). The mean hazard curves from Lorca, Almería, Murcia and Vera also depict PGA increments for the CB14 equation in RPs over ∼100 years ( Figure 10D). In other lower slip rate faults such as BSF, the use of different GMPEs implies practically negligible PGA variations for all RPs (C-C′ cross section along BSF-3; Figure 12E). In this case, for the 475 years RP, the use of Akk14 provides slightly higher PGA values (∼0.05 g) due to the less sensitivity of the equation to the fault slip rates, up to one order of magnitude lower than AMF-1 or CF (Table 1).
In non-vertical dipping faults, hanging wall effects are more accentuated with the CB14 equation (Figures 12C-E). In the zones further away from the fault trace on the hanging wall, hazard decreases more abruptly with CB14 with respect to Akk14, especially in long RPs. The higher PGA gradient around the fault trace causes a rapid drop on the hazard values with respect to the lower gradient of Akk14 as distance to the fault increases. AMF, CAF and BSF are clear examples of this effect and, in particular, BSF is the one affected by highest reductions on the hanging wall for all RPs, because the surface projection of its fault plane is larger ( Figures 12C,D). The curve of Torrevieja represents this effect related to the BSF ( Figure 10D).

Overview: Limitations and Perspectives
The high hazard values in the near-fault zones, the high contribution of faults on the sites hazard and the agreement of the 475 years RP PGA values with the 2011 Lorca earthquake represent more coherent results with the observed values than classical area source methods. Examples of these classical methods are the current Spanish Building Code (NCSE, 2002) and the updated national seismic hazard map (IGN-UPM, 2013). The inclusion of multi-fault ruptures into fault-based PSHA is a novelty for Spain and yet unusual in other low-to-moderate seismicity regions worldwide. Generally, the lack of fault data prevents such kind of analysis in these regions, but as seen, it can imply very important hazard (PGA) changes respect to classical approaches. Moreover, considering multi-fault ruptures is a step forward towards a more realistic seismic hazard assessment, as the occurrence of multi-fault ruptures in nature is becoming an increasingly identified phenomena (e.g., Langridge et al., 2018;Quigley et al., 2019). We therefore aim that our study could serve as a case example for other studies focused in low-to-moderate seismicity regions, such as the EBSZ.
Nevertheless, it is necessary to assess the uncertainties of our study that affect the results presented and need to be considered when interpreting them. One remarkable topic is the background ratio selection of our study. The ratios modeled are an assumption that is established following expert opinion, based on available geological and seismological data from earthquakes in the region. This means that the ratios are not necessarily defining of the study area, but an approximation. The option selected here considers that a big part of the earthquakes in the Mw 4.0-5.5 range are related to the faults modeled because it is consistent with observations (section 3.1.1.2). However, we acknowledge the existent implicit epistemic uncertainties within this assumption, as exemplified by the main shocks of the 2019 Ridgecrest sequence (Mw 6.4 and 7.1), which developed in non-previously mapped faults (Ross et al., 2019). In this sense, other background ratios could and should be explored in a logic tree in further research to account for the epistemic uncertainty of earthquakes happening outside of the known faults. In the line of Chartier et al., 2019, we also recognize the difficulty of setting this parameter and encourage the discussion among researchers on what parameters could help its determination. This is a key issue in SE Spain considering the aforementioned large uncertainties affecting the seismicity epicenter locations of most Mw ≥ 4.0 earthquakes.
Another limitation of our results is the fact that the present hazard calculations may be underestimating the contribution of other fault sources in SE Spain outside of the EBSZ fault system. There are many Quaternary-active faults within the mapping area selected in our study that have not been considered in the source modeling mainly because of a lack of quality data, contrary to the faults within the EBSZ. Fault systems such as El Alquián and Andarax Valley, Polopos or San Miguel de Salinas (Figure 2A) are some examples of active faults in the area with evidence of related seismicity (QAFI database; IGME, 2015a). These faults could also contribute to the seismic hazard, especially outside the EBSZ buffer area, although none of them have the dimensions (length) of the ones in the EBSZ considered in our study and thus their contribution should not be expected to be that relevant.
Further research is then encouraged to focus on incorporating such uncertainties in the light of a more complete and accurate PSHA analysis. As such, in the zones where there is a lack of paleoseismological studies, obtaining well-constrained and updated fault parameters, especially slip rate values, is critical because they can compromise the reliability or bias of the hazard results. This is not only important for the EBSZ faults but also for the ones outside of it and not included in this study. Updated paleoseismic and geodetic data on faults and the inclusion of more faults as sources would necessarily lead to a reassessment of the source modeling and the hazard in the region. In the meanwhile and given the need of updated and more realistic hazard assessments, it would be useful to consider the near-fault hazard estimations in our models as they could help to set nearfault amplification coefficients to be applied at sites located on top of active faults or in close vicinity. Further research should also focus on exploring the use of other parameters rather than PGA to quantify the seismic hazard such as the spectral acceleration or peak ground velocity and displacement. These parameters are more indicative of the damage generated by an earthquake than PGA as they highlight possible amplification effects in buildings and structures at a given site that ultimately are of interest for seismic engineering purposes. In this sense, and additionally to the use of GMPEs, exploring approaches of spatially correlated ground motion intensity measures (e.g., Park et al., 2007)  interesting for further research in order to better characterize expected ground motion, damage and loss estimates for seismic hazard and risk assessment.

CONCLUSIONS
This paper presents a fault system-based PSHA of SE Spain in which earthquake occurrence in the EBSZ is modeled at the scale of the whole fault system and using geological fault data. We devised two alternative models for characterizing the seismicity in the EBSZ, one assuming that all seismic activity is due to known faults of the EBSZ (FM) and another in which part of the seismicity (Mw 4.0-5.5) also occurs as background, related to unknown active faults (FM_bg). The main conclusion is that, for both models analyzed, faults dominate the seismic hazard in the areas located around faults up to 20-30 km from them, even for short return periods of 475 years. This is remarkable for SE Spain as most of the important cities or towns in this area are commonly found on top or very close to Quaternary active faults (e.g., Lorca, Murcia, Totana, Alhama de Murcia, etc.). Fault effects increase PGA values compared to previous PSHAs based on classical area source methods, but their influence zone is strongly dependent on their slip rates. Around the faults with the higher slip rates of the system, such as the Carboneras Fault or the southwestern end of the Alhama de Murcia Fault (>1 mm/yr), the PGA values are increased up to 150% respect to area source methods for all return periods and the fault influence extends up to 10-20 km from the fault traces. On the other hand, lower slip rate faults such as Bajo Segura Fault (∼0.2 mm/yr) produce minor PGA increments, only significant for long return periods, and their influence zones extend up to 5 km from the trace. Around non-vertical dipping faults (e.g., Alhama de Murcia and Bajo Segura faults), the higher hazard values are found in the hanging wall, implying hanging wall effects. Lastly, the faults with the lowest estimated slip rates of the system such as Palomares Fault (<0.1 mm/yr) have no significant impact in the hazard whatsoever. The very similar results from both models evidence that modeling a small portion of the Mw 4.0-5.5 seismicity rate as background (FM_bg) has a minor impact on the seismic hazard of the EBSZ, because fault activity generally overcomes the rates of the background.
The fault-based approach presented seems to provide more accurate hazard estimations than classical area source approaches in sites located close to the faults, at least for the return period of 475 years typically dominated by moderatemagnitude earthquakes. Our hazard results in Lorca for this return period appear to be more consistent with the instrumental PGA values recorded during the 2011 Mw 5.2 earthquake than the ones predicted in the Spanish Building Code or in the updated national seismic hazard map. Also, the consideration of multi-fault ruptures in the source modeling causes significant differences with previous fault-based PSHA in SE Spain and it is a step forward to account for the complexity observed in nature. This kind of approach is a novelty in Spain and might serve as a case example for other low-to-moderate seismicity regions worldwide.
Fault slip rate and its uncertainties have an important influence on the seismic hazard, especially in higher slip rate faults that have more impact in the studied return periods. However, the lack of detailed paleoseismic studies in some lower slip rate faults of the EBSZ (Palomares Fault, Bajo Segura Fault and northeastern termination of Alhama de Murcia Fault) compromises the reliability of the hazard in their surroundings. We identify thus that further research should focus primarily on the obtention of new data to characterize these faults. In the meanwhile, we suggest the use of a near-fault amplification coefficients to be applied to the basic design acceleration provided in the current seismic code for the buildings near Quaternary active faults in this region.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.