A fault-based probabilistic seismic hazard model for Lebanon, controlling parameters and hazard levels

The present work develops a comprehensive probabilistic seismic hazard study for Lebanon, a country prone to a high seismic hazard since it is located along the Levant fault system. The historical seismicity has documented devastating earthquakes which have struck this area. Contrarily, the instrumental period is typical of a low-to-moderate seismicity region. The source model built is made of a smoothed seismicity earthquake forecast based on the Lebanese instrumental catalog, combined with a fault model including major and best-characterized faults in the area. Earthquake frequencies on faults are inferred from geological as well as geodetic slip rates. Uncertainties at every step are tracked and a sensitivity study is led to identify which parameters and decisions most influence hazard estimates. The results demonstrate that the choice of the recurrence model, exponential or characteristic, impacts the most the hazard, followed by the uncertainty on the slip rate, on the maximum magnitude that may break faults, and on the minimum magnitude applied to faults. At return periods larger than or equal to 475 years, the hazard in Lebanon is fully controlled by the sources on faults, and the off-fault model has a negligible contribution. We establish a source model logic tree populated with the key parameters, and combine this logic tree with three ground-motion models (GMMs) potentially adapted to the Levant region. A specific study is led in Beirut, located on the hanging-wall of the Mount Lebanon fault to understand where the contributions come from in terms of magnitudes, distances and sources. Running hazard calculations based on the logic tree, distributions of hazard estimates are obtained for selected sites, as well as seismic hazard maps at the scale of the country. Considering the PGA at 475 years of return period, mean hazard values found are larger than 0.3 g for sites within a distance of 20–30 km from the main strand of the Levant Fault, as well as in the coastal region in-between Saida and Tripoli (≥ 0.4 g considering the 84th percentile). The study provides detailed information on the hazard levels to expect in Lebanon, with the associated uncertainties, constituting a solid basis that may help taking decisions in the perspective of future updates of the Lebanese building code.


Introduction
"Earthquakes don't kill people, buildings do", this often repeated saying highlights the importance of building safely in regions exposed to earthquakes. It is not possible to predict the next destructive earthquake, however countries may adopt building codes in order to design structures capable of withstanding given levels of shaking. The probabilistic seismic hazard assessment (PSHA) method was introduced in the late 1960s (Cornell 1968;Esteva 1968) and became the methodology enforced in most regulations worldwide to determine the reference ground-motions (e.g. Petersen et al. 2014;Grünthal et al. 2018;Beauval and Bard 2021). PSHA consists in combining probabilistic earthquake forecasts with ground-motion prediction models, to determine at the sites of interest the probability that a given ground-motion level is exceeded in the future.
Lebanon is a country with high seismic potential since it is located along the Levant Fault System (LFS), a plate boundary that accommodates the northward motion of the Arabian plate relative to the Sinai subplate (Fig. 1). This ~ 1200 km long left-lateral strike-slip fault connects the Red Sea spreading center to the East Anatolian fault system (Garfunkel 2014;Marco and Klinger 2014). Current estimates about the rate at which the fault accommodates horizontal displacement vary between 4 and 6 mm/year, based on a number of geomorphologic, paleoseismic, and geodetic studies (e.g. Gomez et al. 2007a, b;Le Béon et al. 2012;Wechsler et al. 2018). Major destructive earthquakes occurred along this fault system, as demonstrated in numerous historical, archeoseismic  Table 1  and paleoseismic studies (e.g Guidoboni et al. 1994;Daëron et al. 2007;Elias et al. 2007;Ferry et al. 2007;Ambraseys 2009;Wechsler et al. 2014).
In a PSHA study, the source model must include all earthquakes that we believe may occur in the future, with associated probabilities of occurrence. To build this source model, different observations available over different time and space windows may be combined: past earthquakes, e.g. precise instrumental events as well as more uncertain historical events; mapping of active faults with information on the occurrence of historical and paleoearthquakes and on the slip rates, both from geologic and geodetic studies. Considering only the instrumental period, Lebanon would be classified as a country with low-to-moderate seismic potential. In terms of seismic hazard assessment, these low seismicity levels are not representative of the large destructive events that occurred in the past. As a consequence, we believe that in this region, any source model should include faults and that area models relying essentially on earthquake catalogs should not be applied.
The most recent seismic hazard estimates for this area were delivered by a regionalscale project, the Earthquake Model of the Middle East project (EMME), that extends from Turkey to Pakistan . The EMME seismic hazard calculations relied on a source model made of two main branches, an area source model where frequencymagnitude distributions are established from the earthquake catalog, and a fault model combined with smoothed seismicity. At the scale of Lebanon, a previous hazard study was performed by Huijer et al. (2016), including area sources and faults. Earthquake recurrence models have been established from past seismicity, however no details are provided on how earthquakes were associated to faults and on the number of events available for each fault. Besides, a number of other probabilistic seismic hazard studies exist, performed within private consultancies (e.g. Arango and Lubkowski 2012; AMEC 2012; Erdik et al. 2014). None of these past studies explore in depth the uncertainties associated to the different steps in the process. The present work aims at establishing a fault-based source model, based on faults' characterization in the Lebanese region, as well as up-to-date information on past seismicity (newly derived earthquake catalog Brax et al. 2019;GRAL local catalog) and faults characterization in the Lebanese region, applying a strategy for tracking uncertainties, understanding their impact on hazard estimates, and delivering hazard results representative of these uncertainties.
Our source model for Lebanon includes faults and off-fault seismicity. We combine this source model with a set of ground-motion models to evaluate seismic hazard. The first section of the paper describes the past seismicity, the set of active faults in the area, the information on the deformation accumulating along these faults and the potential historical events they generated. The second section details the building of the source model. To forecast off-fault seismicity two alternative smoothed seismicity models are built from the Lebanese instrumental catalog, taking into account uncertainties on the declustering step and on the magnitude range considered to model the recurrence. Then we describe how we build moment-balanced recurrence models for faults based on the available slip rates and assumptions on the geometry of faults. The third section describes the ground-motion models selected. The fourth section aims at identifying which parameters and decisions bearing uncertainties have a strong influence on hazard estimates, at four cities in Lebanon.

3
A logic tree is established for the source model, and the impact of each parameter or decision on hazard is quantified by sampling different branches of the logic tree. To better understand parameters' impacts obtained in Beirut, the capital city, a disaggregation study is conducted. At last, the final section of the paper aims at deriving seismic hazard maps at the scale of the country, exploring the logic tree populated with the controlling parameters.
2 Seismic potential of the area

Past seismicity
Although large historical earthquakes have struck the Lebanese country centuries ago, the seismicity observed since the beginning of the XXth century is moderate. In the instrumental earthquake catalog built from global catalogs by Brax et al. (2019), four earthquakes with a magnitude larger than 5 produced light damages to buildings in some localities, in 1956 (M W 5.3 and 5.5, Fig. 1a), 1997 (M W 5.1) and 2008 (M W 5.1). Much larger magnitude earthquakes occurred in previous centuries. Brax et al. (2019) critically analyzed the published studies available for pre-instrumental events and delivered a list of events with estimated macroseismic epicenters within Lebanon, including the 551 earthquake (M S ~ 7.3), 1170 (M S ~ 7.3), the 1202 earthquake (M S ~ 7.6), the two events in 1759 (M S ~ 6.6 in October 30 and M S ~ 7.4 in November 25) and the 1837 earthquake (M S ~ 7.0) (Fig. 1a).

Active faults
The LFS can be divided into two North-South trending sections connected by the Lebanese restraining bend, the Jordan Valley fault to the south and the Ghab fault to the north. The restraining bend comprises two mountain ranges, Mount Lebanon and Anti-Lebanon, separated by the Beqaa Valley. Within the bend, the fault splays into three major strands: the Roum fault on the west, the Yammouneh fault in the center, and the Rachaya and Serghaya faults on the east (Fig. 2). The majority of the strain accumulation is accommodated by the Yammouneh fault (Daëron et al. 2004;Gomez et al. 2007a, b). Another fault has been identified, the Mount Lebanon thrust (MLT), an off-shore east-dipping thrust Carton et al. 2009).
In the present work, we will consider only major faults affecting the Lebanese region, secondary fault segments that are not well-characterized are not included in the model (Fig. 2). A fault can be included in the hazard analysis if the fault geometry can be defined and if its seismogenic potential can be quantified.

Yammouneh fault
The Yammouneh fault is the throughgoing structure that connects the northern and southern sections of the LFS and forms the eastern boundary of Mount Lebanon. The fault has a left-lateral strike-slip mechanism and extends over ~ 200 km. Daeron et al. (2004) estimated a 25-ka mean slip rate between 3.8 and 6.4 mm/yr, based on cosmogenic dating of offset alluvial fans. Analyzing paleoseismological trenches in lacustrine deposits, Daëron et al. (2005) provide evidence of coseismic slip on this fault corresponding to the A.D. 1202 earthquake (M S 7-7.8, Ambraseys and Jackson 1998) and show that this segment has remained locked since then. Moreover, Daëron et al. (2007) expanded the trench work across the Yammouneh paleolake and identified ten to thirteen paleo-events (M ~ 7.5), extending back more than ~ 12 kyr (the latest being the 1202 event). The paleoseismic study of Gomez et al. (2007b) provides additional cosmogenic age constraints and confirms a slip rate between 3.9 and 6.1 mm/yr. Besides, based on the geodetic velocity field and a 1D elastic dislocation model, Gomez et al. (2007a) quantitatively assess an average slip rate of 4.4 mm/year along the fault.

Missyaf fault
To the North of the restraining bend, the Missyaf segment is ~ 90 km long, limited by the Yammouneh Fault to the south and the Ghab basin to the north . Based on a combined analysis of archaeological excavations and paleoseismic trenches, Sbeinati et al. (2010) establish the occurrence of four faulting events in the last ~ 3500 yr (magnitude 7.3-7.5), the last one being the 1170 destructive earthquake, as well as a range of 4.9-6.3 mm/yr for the slip rate along this segment. The timing of the most ancient event in the sequence is poorly constrained. Assuming that this ancient event corresponds to an historical earthquake in 1365 BC, they propose a range of 4.9-5.5 mm/year for the slip rate (Sbeinati et al. 2010). On this part of the Levant Fault, a significant discrepancy is observed between slip rates based on geological observations and rates estimated from geodetic measurements or modeling, which are much lower, Alchalbi et al. (2010) finds 2.5 ± 0.7 mm/year; Gomez et al. (2020) 2.2 ± 0.5 mm/year.

Roum fault
The Roum fault is a left-lateral strike-slip fault, that bounds the Mount Lebanon range to the south (Fig. 2). Nemer and Meghraoui (2006) provide a detailed mapping of the fault trace from the western edge of the Hula pull-apart up to the Jarmaq basin (latitude 33.4°). To the north of this basin, they acknowledge that the trace is less clear. We consider the trace of Nemer and Meghraoui (2006) up to latitude 33.5°, and extend the fault towards the north west along a constant strike, based on morphological analyses (Elias 2006;. At a paleoseismic site on the southern part of the fault (Fig. 2), Nemer and Meghraoui (2006) establish the occurrence of at least 4-5 surface rupture events with magnitude > 6.5 during the last ~ 10ky. For the most recent event, they provide a lower bound for the date (post A.D. 84-239), and suggest that it might be the 1837 historical earthquake (M S ~ 7.0, Ambraseys 2006). Assuming that the stream offsets measured close to the trench correspond to the accumulated coseismic displacements in the last ~ 10ky, they estimate a range of 0.86-1.05 mm/year for the slip rate.

Mount Lebanon fault
A geophysical survey has shown evidence for submarine seismogenic thrusts off shore central Lebanon that would be responsible for the A.D. 551 earthquake ). This east-dipping thrust system is called Mount Lebanon Thrust (MLT). Its geometry is complex and it is limited to the south and north by oblique lateral ramps. Elias et al. (2007) evaluates a rupture length of at least 100 km, possibly 150 km, corresponding to an earthquake with magnitude 7.4-7.6 according to the scaling relationship of Wells and Coppersmith (1994). From shoreline uplifts, they evaluate an average recurrence time of at least 1500 years for an A.D. 551-type event since the early Holocene, and an average slip rate of 1-2 mm/yr for the fault. Here we use the fault trace used in Danciu et al. (2018) prolongated to the south (inclusion of the Saida lateral ramp) and to the north (inclusion of the Tripoli and Aakkar thrusts on land, Elias et al. 2007).

Serghaya fault
The Serghaya fault is located approximately along the Syrian-Lebanese border. This strikeslip fault follows the Anti Lebanon Range east of the Beqaa plain. Meghraoui et al. (2003) and Nemer and Meghraoui (2008) perform a detailed mapping of the fault. The fault trace taken into account here extends from latitude ~ 33.4° to latitude ~ 34.2° (Fig. 2). We do not include the segment bordering Mount Hermon to the south, as the fault trace there is less clear. Based on trench excavations (Fig. 2), Gomez et al. (2003) identify five surface rupturing events within the past ~ 6500 years, with displacement of ~ 2 m and estimated magnitudes larger than or equal to 7.0. The last event postdates A.D. 1650 and may correspond to the 1759 earthquake (M S ~ 7.4, the November shock, that destroyed many villages in the Beqaa and ruined Baalbeck; Ambraseys and Barazangi 1989;Daëron et al. 2005). From the buried displaced channels in the trenches, Gomez et al. (2003) estimate a Holocene slip rate of 1.4 ± 0.2 mm/yr.

Rachaya fault
The Rachaya fault extends along the western part of the Anti-Lebanon range for approximately 67 km. Nemer and Meghraoui (2008) provided a detailed mapping for this fault segment as well as geomorphic evidence for recent faulting. Here we consider a trace that splays from the eastern bounding fault of the Hula Basin up to latitude ~ 33.5°. Nemer and Megrahoui (2008) excavated a trench that revealed the occurrence of at least two seismic events. The last event belonging to the period 1686-1924 A.D. could thus be correlated to the first 1759 shock (M S ~ 6.6, Ambraseys and Barazangi 1989), consistent with the conclusions of Daeron et al. (2005). This fault constitutes the southern continuation of the Serghaya fault. Both faults connect through the compressional jog that forms the Mount Hermon. Hence, we make the hypothesis that they have the same slip rate.

Jordan valley fault
The fault is limited to the north by the Hula basin and to the south by the Dead Sea pullapart (e.g. Ferry et al. 2011). Many studies suggest that the Jordan Valley fault slipped during the earthquakes of A.D 363 (M S ~ 7.4), 746 (M S ~ 7), 1033 (M S ~ 6) and 1546 (M S ~ 6) (see Fig. 3 and Table S1 in Brax et al. 2019). The Jordan Valley fault has been more studied than the northern part of the LFS and is therefore better characterized. Different segmentations have been proposed (e.g. Ferry et al. 2011;Lefevre et al. 2018). A wealth of paleoseismic data that provided slip rates along the fault (e.g. 4.1 +0.4 / −1 mm/yr north of the Sea of Galilee, over the last 4000 years, Wechsler et al. 2018;or 4.9 ± 0.3 mm/yr between the Sea of Galilee and the Dead Sea, for the last 25 ka, Ferry et al. 2011). Slip rates estimated with geodetic models (e.g. 4.7 ± 0.5 mm/yr in Al Tarazi et al. 2011;4.1 ± 0.8 mm/ yr in Hamiel et al. 2016) are overall consistent with the geological slip rates. We consider a unique long segment of ~ 180 km associated with an average slip rate along the fault that ranges from 4 to 5 mm/yr.

3 3 Earthquake recurrence model
The source model must describe the locations of future earthquakes, together with their magnitudes and probabilities of occurrence over future time windows. In the Lebanese region, the seismicity in the instrumental period can be qualified as typical of a low-tomoderate seismicity area. In terms of seismic hazard assessment, the seismic rates inferred from the instrumental catalog are not representative of the true seismic potential in the area. We believe that it is not meaningful to derive area source models that heavily rely on earthquake catalogs. On the other hand, faults are rather well characterized and their activity has been demonstrated. We propose here to develop a source model made of faults combined with off-fault seismicity. Occurrences of earthquakes on fault are constrained by the seismic moment rates inferred from the available slip rate estimates. Occurrences of offfault earthquakes are based on a smoothed seismicity model that relies on the instrumental earthquakes recorded by the local network.

Modeling off-fault seismicity
To forecast earthquakes in the background, a smoothed seismicity model is developed.
Smoothed seismicity models rely on the assumption that future seismicity occurs in areas close to where past seismicity occurred (e.g. Frankel 1995;Stock and Smith 2002). Within Lebanon, the global catalog homogenized in moment magnitude (Brax et al. 2019) contains 23 events with magnitudes M W ≥ 4.1 (spatial window: longitude 35-36.8 latitude 33-34.7), most of them belong to three earthquake sequences that occurred in 1956, 1997. If this catalog is used to establish the smoothed seismicity model, the model will forecast highest probabilities of occurrence at the locations of these sequences. We decided to use instead the Geophysical Research Arrays of Lebanon (GRAL 2019) local catalog that covers the time window 2006-2019 and includes duration magnitudes down to M D 2.4 (3677 events within the spatial window of Fig. 1b). These low-magnitude events highlight areas of recent seismic activity. Although the time window is short, we assume that earthquake densities in this low-magnitude range may be a good guide for future off-fault seismicity. Equivalence between duration magnitude and moment magnitude needs to be assumed as there are too few events with both magnitude estimates to establish a correlation (considering M W from Global Centroid Moment tensor (CMT), Brax et al. 2019).
Historical events are not included in the earthquake catalog used for smoothing for two reasons: the epicentral location of a large magnitude historical event has a limited meaning (in comparison with e.g. the extension of an intensity isoseismal), and most of these large historical events are believed to have been generated on faults (earthquakes in 551, 746, 1157, 1170, 1202, 1759 and 1837, see Fig. 1 and Sect. 2.2).
Smoothed seismicity rates are estimated from the earthquake catalog applying the algorithm developed by Hiemer et al. (2013Hiemer et al. ( , 2014SEIFA). Seismicity rates account for completeness in time. The density of seismicity in each spatial cell is obtained by smoothing the location of each earthquake with an isotropic variable power-law kernel. We use an adaptive smoothing distance (distance to the second closest neighbor). The b-value is evaluated from the Gutenberg-Richter frequency-magnitude distribution (Gutenberg and Richter 1944) obtained at the scale of Lebanon from the GRAL catalog.
To identify clustered events in the catalog, two well-known and widely applied declustering algorithms are considered. The Gardner and Knopoff (1974, GK74) algorithm relies on simple magnitude-dependent windows in time and space. The Reasenberg (1985, R85) algorithm is a linking algorithm where clusters are allowed to grow in time and space. In time, it is based on the Omori's law, which describes the temporal decrease of the number of aftershocks after the occurrence of the mainshock. In space, the algorithm assumes a circular rupture radius for each event and searches in this area for clustered events.
We test different parameterizations for the application of these algorithms. For Reasenberg, we consider the original parameters used in California (Reasenberg 1985), as well as three alternative values producing a stronger declustering (see Table 2, range of parameters from Christophersen et al. 2011). Four frequency-magnitude distributions are established based on the declustered catalogs obtained (Fig. 3a). All distributions display two different slopes, below and above M D 3.3 (respectively ~ 2.8 and ~ 1.2). In many places, the analysis of the number of events versus hour of the day shows an abnormal distribution, with many more events during working hours, suggesting the presence of non-tectonic events (examples in Figs. 3c and 3d). These swarms are only partially removed in the declustering (see declustering D4 in Fig. 3d, 46% of earthquakes with M D ≥ 3.3 are identified as clustered events). A specific analysis would be required to identify recordings that have a high probability of being anthropogenic events (e.g. quarries and mines); however this is not possible within the timeframe of this project. Besides, we also estimate seismic rates considering only night events (see e.g. Gulia and Gasperini 2021), and observe that the rates above M D 3.3 remain quite stable (Fig. 3b).
The Gardner and Knopoff algorithm is applied considering two alternative sets of windows in time and space, the original windows from GK74 and the more conservative windows proposed by Burkhard and Grünthal (2009, BG09) for Europe (Fig. 4). Gardner and Knopoff algorithm removes many more events than the Reasenberg algorithm. Considering M D ≥ 2.8, the GK algorithm removes 62% earthquakes identified as clustered events with the GK74 windows, and 84% with the BG09 windows. Frequency-magnitude distributions are estimated considering all events and night event only. Except for the declustering with BG09 on night events only, these distributions still display two slopes (Fig. 4a). Considering only night events and the BG09 windows, a unique b-value of 1.3 is found over the magnitude range 2.8-4.4; the rates for earthquakes M D ≥ 3.3 remain rather stable.
Based on these different tests, we decided to consider two declustered catalogs that lead to two earthquake recurrence models and two alternative smoothed seismicity models to be included in the source model logic tree: • The catalog declustered with R85 algorithm, set of parameters D4 ( Table 2). The frequency-magnitude distribution displays a slope equal to 1.2 from magnitudes 3.3 on. The smoothed seismicity model 1 is established from events with M D larger or equal to 3.3 (Fig. 4b).
• The catalog declustered with GK74 algorithm applying the windows of BG09, considering only night events. The frequency-magnitude distribution displays a slope equal to 1.3 from magnitudes 2.8 on. The smoothed seismicity model 2 is established from events with M D larger or equal to 2.8 (Fig. 4b). These two models rely on rather close annual rates in the magnitude range 3.2-4.2, but provide different spatial densities (Fig. 5). The models are built with a grid of 0.05° × 0.05° in space. We allow magnitudes up to 6.5 in the background, which means we assume that such magnitudes can occur anywhere on faults presently unknown.

Earthquake forecast on faults
From the knowledge available on faults, earthquake recurrence models can be derived. We use the average slip rates estimated from geologic studies to determine the seismic moment rates that can be released in earthquakes. Slip rates are determined either from paleoseismic trenches when they document enough earthquakes to allow averaging the slip variability for successive earthquakes (Wechsler et al. 2018), or from offset geomorphological features that provide a slip rate averaged over several thousands of years (Daeron et al. 2004). Paleoseismic studies mostly provide evidence for the large magnitude earthquakes with surface ruptures. Moderate events that rarely produce surface ruptures are often not detected. However, moderate events are more frequent than large events, and they also contribute to seismic hazard assessment (see e.g. Beauval et al. 2006).
To establish a moment-balanced recurrence model, we need to assume how earthquake sizes distribute. For the faults in Lebanon, there are not enough observations available to discriminate between a Gutenberg-Richter exponential decrease of earthquakes with magnitudes and a more characteristic behavior (see e.g. Ishibe and Shimazaki 2012;Stirling and Gerstenberger 2018). The GRAL instrumental catalog (14 years) contains many  Reasenberg (1985) algorithm, 75 events with M D ≥ 3.3 (squares), (b) GRAL catalog declustered with Gardner&Knopoff (1974) algorithms, window parameters from Burkhard and Grünthal (2009), considering only periods at night, 148 events with M D ≥ 2.8 low-magnitude events, but their spatial distribution does not seem aligned with the major fault segments (Fig. 1b). There are not enough historical events associated to the faults (see Fig. 1a) to estimate meaningful mean recurrence times.
We reviewed a number of published seismic hazard models. The decision on which distribution to use for faults is most often not detailed, some studies use the exponential distribution (e.g. Rong et al. 2020;Danciu et al. 2018;Wiemer et al. 2016;Poggi et al. 2020;Meletti et al. 2021); whereas others use a characteristic model (e.g. Stirling et al. 2012;Wang et al. 2016;Demircioğlu et al. 2018). In its simplest form, the characteristic model assumes that each fault source produces a single earthquake magnitude that ruptures the full fault source length (e.g. Stirling and Gerstenberger 2018). Here we consider both recurrence models, namely the exponential form 2 of Anderson and Luco (1983) and the characteristic model of Youngs and Coppersmith (1985).
The seismic moment rate is estimated from the slip rate assumed constant over the fault area: with Ṡ the slip rate per year, μ the rigidity modulus (taken as 3.6 × 10 10 N/m 2 ), and A the surface of the rupture area.
Using the exponential form 2 of Anderson and Luco (1983), the annual rate N of events with magnitude larger or equal to m can be written as: with b the exponential coefficient and M max the maximum magnitude corresponding to the largest earthquake that may break the fault. We use the universal b-value of 1 (Kagan 1999) for two reasons: the b-values inferred from the instrumental GRAL catalog in duration magnitude are abnormally high; using the Brax et al. (2019) catalog at a regional scale we found a b-value close to 1.
For the model to be moment-balanced, the a-value must be estimated as follows (see e.g. Mariniere et al. 2021): with parameters c and d from the moment-magnitude definition ( log 10 Ṁ 0 = cM W + d ; with c = 1.5 and d = 9.1, for moment in units of N m, Hanks and Kanamori 1979).
The Youngs and Coppersmith (1985) characteristic model combines an exponential model for moderate magnitudes and a characteristic model for large magnitudes, with rates of characteristic earthquakes in the magnitude range M max -0.5 to M max much higher than predicted by the moderate magnitude range. Its probability density function can be written as follows (see e.g. Convertito et al. 2006): Then, for the recurrence model to be moment-balanced, the cumulative rate of characteristic magnitude can be calculated according to the following equation (see e.g. Convertito et al. 2006): Fault planes need to be defined. A dip of 45° is assumed for Mount Lebanon Thrust ); whereas all strike-slip faults are assumed vertical with a dip of 90°. As the extension at depth of these strike-slip crustal faults is not well constrained, we consider three alternative maximum depths to account for this uncertainty (10,14 and 18 km), according to the range found for locking depths in elastic block modeling studies (e.g. Gomez et al. 2007a). For Mount Lebanon Thrust, the depth is fixed to 18 km but the uncertainty on the dip is considered (40, 45 and 50°).
Maximum magnitudes considered for each fault are reported in Table 1. The uncertainty on this parameter is significant. We identify the largest historical earthquake that has been associated to each fault, and its magnitude range inferred from macroseismic intensity data (see Brax et al. 2019). The maximum magnitude M max on the fault must be larger or equal to the magnitude of the historical event, accounting for uncertainties. For strike-slip faults, the Leonard (2014) scaling relationships are applied to infer maximum magnitudes from the total length of the fault (Table 1). Then, we consider the mean value predicted by the scaling relationship, as well as the mean plus one standard deviation. For thrust faults, Leonard (2014) scaling relationships imply that the maximum rupture that the fault can host is limited by its width. Therefore, the length of the maximum rupture is first obtained from the width (applying L = ( W 12 ) 1.5 , considering e.g. a width 25 km, dip of 45° and depth 18 km). Then M max is inferred from the length considering the mean as well as mean plus one standard deviation (see Table 1).
To account for the uncertainty on the slip rates, an interval is considered (Table 1). This interval corresponds to the range provided by the paleoseismic studies available for Missyaf (4.9-5.5 mm/yr, Sbeinati et al. 2010), Yammouneh (3.9-6.1 mm/yr, Gomez et al. 2007b), Roum (0.86-1.05, Nemer and Meghraoui 2006), and Serghaya (1.2-1.6, Gomez et al. 2003). As there is no direct information for Rachaya, the slip rate is assumed to be equal to the slip rate evaluated for Serghaya. The Jordan valley fault has been studied extensively and different slip rate estimates are available along the fault. We infer a mean range from these different estimates that we apply to the whole fault (4 to 5 mm/yr, Table 1). Figure 6 displays the recurrence models obtained for the 7 faults, considering the exponential Gutenberg-Richter model (first row) or the characteristic model (second row). In this example, the total seismic moment rate is inferred from the lower bound of the slip rate estimates applied to the fault area. As expected, the characteristic model leads to larger rates in the upper magnitude range, and lower rates in the moderate magnitude range, with respect to the exponential model. The impact on the rates of the M max (5)   Table 2 Four sets of parameters required by the Reasenberg (1985) algorithm (  choice is also shown. For a moment-balanced recurrence model with fixed moment rate, a decrease of M max leads to an increase of the rates in the moderate magnitude range. To build the source model, the fault model is combined with the smoothed seismicity model. To avoid double counting earthquakes, a buffer zone of 5 km around the projection of the rupture plane to the surface is considered. Within this buffer zone the frequencymagnitude distributions of the gridded seismicity are truncated at the minimum magnitude taken into account on the fault (following e.g. Danciu et al. 2018;Demircioğlu et al. 2018).
Here we consider two alternative minimum magnitudes, 5.5 and 6.0. This parameter, which is chosen rather arbitrarily, has an impact on hazard results. This minimum magnitude considered for the faults varies from one study to the other: M W 5.5 in Danciu et al. (2018) in the Middle East EMME project, M W 6.0 in Demiciorglu et al. (2017) in Turkey, M W 6.4 in  Youngs and Coppersmith (1985) model, minimum bound for M max ; d) characteristic Youngs and Coppersmith (1985) model, maximum bound for for M max . For a given fault, the four models displayed correspond to the same seismic moment rate, estimated considering the minimum bound of the slip rate interval and the fault area (maximum depth 14 km for strike-slip faults, dip 45° for MLT). Seismic rates are normalized per km 2 . See Table 1 for faults acronyms Woessner et al. (2015) in Europe, or M W 6.5 in Rong et al. (2020) in China. Earthquakes on faults contribute to the hazard from this minimum magnitude up to the maximum magnitude attributed to the fault.

Ground motion prediction
There are no ground-motion records available in Lebanon to guide the selection of the ground-motion prediction model. Any site within Lebanon is located within 20 km of a fault, either strike slip for faults inland or thrust offshore (Mount Lebanon). Disaggregation studies show that for the return periods larger than or equal to 475 yrs, most contributions come from source-site distances lower than 50 km (see Sect. 5). We select three ground-motion models established from shallow crustal earthquakes, relying on three different datasets, as an attempt to capture the epistemic uncertainty on the prediction of ground motions.
The model of Chiou and Youngs (2014) has been established from the NGA-West2 ground-motion database (Ancheta et al. 2014) complemented by numerical simulations. This database contains mainly California earthquakes, supplemented with recordings from magnitudes M W 6 to 7.9 events that occurred in other active regions (Iran, Italy, Japan, New Zealand, Taiwan, Turkey, China). These recordings have been used to constrain the magnitude scaling and the aleatory variability for large magnitudes. The functional form accounts for styles of faulting effects, and for the position of the top of rupture. The model has additional terms to account for hanging-wall effects, as well as effects of rupture directivity at large magnitudes. It accounts for magnitude-dependent effects of extended ruptures on distance scaling in the rupture distance range of 0 km to approximately 100 km. Equations are provided at spectral periods from 0.01 to 10 s.
The model of Akkar et al. (2014) relies on a pan-European database including strongmotion recordings obtained in the Mediterranean region and the Middle East. It provides equations for shallow (< 30 km depth) crustal earthquakes, at spectral periods from 0.01 to 4 s. Its distance metric is the Joyner and Boore distance (shortest distance to the surface projection of the rupture). In the database, magnitudes up to roughly M W 7 are well represented, particularly for normal and strike-slip faulting. For larger magnitudes there are almost no records from reverse-faulting events, the available data are mainly from three large strike-slip earthquakes (largest magnitude M W 7.6). Reverse-faulting earthquakes are quite poorly represented (in contrast with the NGA model for which reverse earthquakes contribute a large proportion of the database).
The third model considered is the Kotha et al. (2020) ground-motion model developed for shallow crustal earthquakes in Europe, with magnitudes M W up to 7.4. It is based on the new European Strong-Motion datataset (Lanzano et al. 2019). The model includes a region-dependent coefficient of attenuation, and a region-dependent source-scaling factor. However, too few data were available both in the attenuation regionalization and event localization polygons including Lebanon, therefore we apply the ergodic version of the model. Distance measure is Joyner and Boore distance, owing to significant uncertainty in the finite rupture dimensions and hypocentral depths. Events have been categorized into three depth bins (< 10 km, 10 to 20, and > 20 km). Equations are provided at spectral periods from 0.01 to 8 s.

Identification of key parameters/key decisions
We start with building a logic tree to explore uncertainties on the source model and understand which parameters, models or decisions, impact hazard estimates. This logic tree is made of the two alternative smoothed seismicity models, lower and upper bounds for the slip rate accumulating on the fault, two alternative earthquake recurrence models for the faults, two alternative b-values (the universal value 1 and an arbitrary value fixed to 0.85, simply to test the sensitivity of the results with respect to this parameter), two alternative minimum magnitude on faults, lower and upper bounds on the maximum magnitude that can break the fault, and three geometries of fault planes (Fig. 7). For strike-slip faults, three alternative maximum depths of fault planes are considered; whereas for Mount Lebanon Thrust three alternative values for the dip are considered. We combine this source model logic tree with the three ground-motion models equally weighted.
Hazard calculations are led considering a minimum magnitude of M W 4.5, a maximum distance of 150 km, the scaling relationship Leonard (2014) for generating extended rupture planes, and a truncation of the gaussian predicted by the ground-motion models at 3 standard deviations above the mean. The OpenQuake engine, developed by the Global Earthquake Model Foundation, is used (Pagani et al. 2014). Figure 8 displays the results for four cities in Lebanon (Fig. 2), for the PGA at 475 years return period, considering generic rock site (V S30 = 760 m/s). The black distribution is obtained by sampling the full logic tree (2 6 × 3 = 192 source models, combined with 3 Fig. 7 Source model logic tree. In Sect. 5, the logic tree is explored considering all branches as equally probable, delivering 2 × 2 × 2 × 2 × 2 × 2 × 3 = 192 different branches' combinations, combined with 3 ground-motion models, resulting in 576 hazard results. In Sect. 6, the branches in grey are discarded, and 2 × 2 × 2 × 3 = 48 branches' combinations are kept to derive seismic hazard maps (144 hazard values) GMMs, resulting in 576 hazard estimates). In the capital city Beirut, located on the coast on the hanging-wall of Mount Lebanon Thrust, the mean PGA obtained is 0.3 g. In Saida and Tripoli, coastal cities close to the MLT fault but not above it, the mean PGA obtained is around 0.25 g. Higher values are obtained in Zahle in the Beqaa Valley (0.42 g mean value), located along the Yammouneh strike-slip fault. Besides, the smallest variability on hazard levels is found for Saida (0.17-0.3 g, percentiles 16th-84th), whereas the largest variability is obtained in Zahle (0.28-0.57 g) and Beirut (0.17-0.45 g).
Following Beauval et al. (2020), different distributions for the acceleration at 475 years (PGA) are displayed, corresponding to different branches of the logic tree sampled (Fig. 8). We choose this representation rather than an impact analysis (as in Beauval and Scotti 2004) or a Tornado diagram (e.g., Anderson 2018), to determine an absolute impact which is not dependent on a specific set of reference parameters. To quantify the impact of the uncertainty on the slip rate, the 576 values are split into two separate distributions, one that relies on the slip rate lower bound (light blue, 288 values), and the other that relies on the slip rate upper bound (dark blue, 288 values). The difference in the mean values of the two distributions corresponds to the impact on hazard of slip rate uncertainty. Similarly, the impact on hazard of the following uncertainties or parameter choices can be evaluated: uncertainty on the geometry of faults, choice of the recurrence model, uncertainty on the maximum magnitude, the decision on the minimum magnitude on faults, decision on the b-value, alternative smoothed seismicity models. The larger the distance between the two mean values of distributions (or the three in case of geometry), the larger the impact the parameter has on hazard. This analysis shows that the choice of the recurrence model for faults is the decision that impacts the most the hazard values, for all sites. In Beirut, the exponential model leads to hazard values (0.4 g) twice as high as the characteristic model (0.2 g). The uncertainty on the slip rate and on the M max also impacts significantly hazard estimates (~ 50% increase in Beirut). Using a minimum magnitude of 5.5 rather than 6 produces an increase of hazard estimates (~ 20% increase in Beirut). The geometry has a stronger impact in Zahle than in the coastal cities. The uncertainty on the b-value or on the smoothed seismicity model has a negligible impact on hazard. This hierarchy of parameters' impacts is the same for the four cities, it is also stable considering other spectral periods (0.2 s and 1 s) and larger return periods (2475 and 4975 yrs) (see Online Resource, Figs. S1 and S2). of the recurrence model on earthquake rates, exponential (grey) or characteristic (black), for a fixed moment rate, and (d) impact on the UHS. All recurrence models correspond to the same seismic moment rate, estimated from the slip rate upper bound (6.1 mm/yr) applied to the Yammouneh fault that extends down to 10 km depth It is expected that the hazard levels increase when the slip rate increases, or when the fault area is enlarged, as it leads to larger seismic moment rates to be released in earthquakes. Some impacts are less straightforward. At 475 years, hazard values increase when M max values on faults decrease (at the PGA, increase of 0.36 g to 0.5 g in Zahle, selecting the M max lower bound rather than the upper bound, Fig. 8). For a fixed moment rate, a decrease of the maximum magnitude leads to an increase of seismicity rates in the moderate magnitude range (Fig. 9a, Yammouneh fault), and an increase in the amplitudes of the hazard spectrum at 475 yrs at a nearby site (Fig. 9b, Zahle city). Similarly, hazard values decrease when the characteristic model is chosen over the exponential model (Fig. 8, PGA at 475 years, decrease of 0.52 g to 0.31 g in Zahle). For a fixed moment rate, the characteristic model leads to lower seismic rates in the moderate magnitude range with respect to the exponential model (Fig. 9c, Yammouneh fault), and a decrease in the amplitudes of the hazard spectrum at 475 yrs at a nearby site (Fig. 9d, Zahle). However, for the return period 2475 years, at spectral periods larger or equal to 1.0, opposite trends are observed, acceleration levels increase. In this case most contributions come from the upper magnitude range.

Beirut, hazard controlled by Mount Lebanon fault
Beirut is the capital city of Lebanon. It is located on the hanging-wall of Mount Lebanon fault, at ~ 30 km from the Yammouneh fault. At 475 years return period, the exploration of the logic tree leads to a mean value of 0.3 and 0.72 g, respectively at the PGA and spectral period 0.2 s (Fig. 8 and Fig. S3). Considering the percentiles 16 and 84 th , PGA values vary between 0.17 and 0.45 g. The uncertainty obtained exploring the logic tree results in a large variability of hazard estimates. It is important to understand where the contributions to hazard come from. Figure 10 displays the hazard curve obtained in Beirut considering the full model (in red), as well as the hazard curves considering only one source (Mount Lebanon fault, Yammouneh fault, smoothed seismicity model). Considering the exponential Gutenberg-Richter model, for return periods larger than or equal to 475 years (annual rate of exceedance lower or equal to 0.0021), all the contributions come from the Mount Lebanon fault. Considering the characteristic Youngs and Coppersmith model, lower hazard levels are obtained and at 475 years return period, the Yammouneh fault also contributes to hazard in Beirut. The smoothed seismicity has a negligible contribution to hazard. A disaggregation study in magnitude and distance is also useful to understand which earthquakes at which distances contribute to the hazard (Fig. 11). At 475 years, for the PGA, assuming that earthquake occurrence on faults follows an exponential model, most contributions come from sources at less than 20 km, with magnitude earthquakes between 5.5 and M max , on Mount Lebanon rupture plane. Assuming instead that earthquake occurrence on faults follow a characteristic model, there are also contributions from the more distant Yammouneh fault, from magnitude earthquakes both in the moderate and upper magnitude range.
Mount Lebanon is an off-shore fault, much less well characterized than the inland faults ). The slip rate, the fault trace and fault geometry, as well as extension at depth, are less constrained than for the other faults included in the model. The regional Fig. 11 Beirut site, disaggregation in magnitude (left column) and distance (right column), for the PGA at 475 years return period. First row: exponential model (PGA 0.39 g); second row: characteristic model (PGA 0.19 g). The source models and hazard calculations used are identical as in Fig. 10. The Chiou and Youngs (2014) ground-motion model is used 1 3 EMME project Danciu et al. 2018) obtained a higher mean hazard value in Beirut, close to the percentile 84th we obtain (0.41 g, Fig. 12). This difference is due to a combination of decisions. For establishing earthquake recurrences on faults, the exponential Gutenberg-Richter model alone was considered, only the central segment was considered (latitude 33.76° to 34.34°) thus reducing the total length (with a maximum magnitude range 7.3-7.6), and a unique minimum magnitude was considered for faults (M W 5.5). Moreover, a different set of ground-motions has been considered: Chiou & Youngs (2008), Zhao et al. (2006), Akkar et al. (2014) and Akkar and Çaǧnan (2010), respectively with the weights 0.35, 0.1, 0.35 and 0.2.
If we consider only the exponential model in our logic tree, we get the distribution displayed in Fig. 12, the mean PGA at 475 years increases up to 0.4 g. Then, if the segment is reduced to its central part with M max 7.3-7.6, the mean PGA increases up to 0.44 g. Next, considering only a minimum magnitude 5.5 for earthquakes on faults, the mean PGA reaches a mean value around 0.5 g. At last, applying the EMME ground-motion logic tree, the mean PGA reduces to 0.43 g, close to the mean value obtained in EMME. Nonetheless, the EMME source model includes two branches, area source and fault model . Using the xml OpenQuake input files, available on the efehr.org website, we have calculated the uniform hazard spectrum in Beirut, considering the mean as well as independently the EMME area model, and the EMME smoothed seismicity and fault model (Fig. S7). Results show that a deeper analysis would be required to understand the differences between the present study and the EMME study.

Fig. 12
Hazard estimates in Beirut, distribution of accelerations at the PGA, return period 475 years, mean and percentiles 16th and 84th. Calc1: full logic tree (Fig. 7), hazard results as in Fig. 8. Calc2: only the exponential model is considered. Calc3: one modification with respect to Calc2, Mount Lebanon fault is reduced to the central segment with M max 7.3-7.6 as in EMME. Calc4: one further modification with respect to Calc3, only the minimum magnitude 5.5 is considered on faults. Calc5: same as Calc4 with the GMM logic tree used in the EMME project. The dashed line corresponds to the mean PGA value obtained by EMME for Beirut (0.41 g, Sesetyan et al. 2018)

Seismic hazard maps
We aim at providing hazard values for the whole country, so that this work can be used in the perspective of future updates of the Lebanese building code. We exclude from the logic tree the branches that have a negligible impact on hazard. A unique b-value of 1.0 is used for the recurrence model on faults; only the smoothed seismicity model 1 is considered. Our aim is to display the range of hazard levels that are obtained, and to show how some decisions on the source model logic tree may impact seismic hazard levels at the scale of the country.

Hazard maps at the PGA and 475 years, based on our logic tree
A distribution of seismic hazard maps can be obtained from the exploration of the source model logic tree combined with the ground-motion logic tree (48 × 3 = 144 different branches' combination, equally weighted). From this distribution, a mean hazard map as well as hazard maps at given percentiles can be displayed, for any return period or spectral period. Figure 13 displays the mean hazard map (Fig. 13a) and the 84th percentile hazard map (Fig. 13b) obtained at the PGA for the return period 475 years. Mean hazard values are higher than 0.2 g throughout most of the country (Fig. 13a). Highest levels are found at sites along the main strand of the Levant fault. South of Lebanon, the 0.3 g contour is at a distance of ~ 10 km from the Jordan Valley fault trace. North of 33° latitude, the area with hazard values larger than 0.3 g widens, as the contributions of the Rachaya fault and then the Roum fault add onto the contribution of the Yammouneh fault. Hazard values larger than 0.3 g are also found along the coast, from north of Saida to south of Tripoli, due to The slip rate range associated to this fault (4.9-5.5 mm/yr), combined with an M max of 7.1-7.5 inferred from the length of the segment (75 km), lead to rates in the moderate magnitude range that are higher than those estimated for the Yammouneh and Jordan Valley fault ( Fig. 5; higher mean slip rate, shorter segment, smaller M max , with respect to Yammouneh fault). As for hazard values corresponding to the 84th percentile, they are larger than 0.3 g within most of the country (Fig. 13b), similarly with highest levels along the Yammouneh (> 0.6 g) and Missyaf (> 0.7 g) faults. Mean seismic hazard maps have also been produced for the spectral periods 0.2 s and 1 s, for 475 years and 2475 yrs return periods (see Online Resource,Figs. S5 and S6). To appraise the hazard levels at different spectral periods, the uniform hazard spectra for four cities in Lebanon are displayed, for the return periods 475 and 2475 years (Fig. 14). Also, to appraise the shape of the distribution obtained for a given spectral period, the distribution of hazard levels in Beirut for the PGA and spectral period 0.2 s are displayed in the Online Resource (Fig. S8).
The analysis of input parameters' impacts on hazard has shown that the choice of the recurrence model for faults impacts the most the hazard levels (Sect. 5.1). In many countries, in order to deliver national seismic hazard maps, an expert committee is set up and an elicitation process is led to make the final decisions on the exact input parameters to be  Table S1) used in hazard calculations, as well as branches' weights (e.g. Gerstenberger et al. 2020). Such an expert committee might be established in the future in Lebanon. Decisions on branches weights have obvious consequences on hazard. If the Youngs and Coppersmith characteristic model is favored over the Gutenberg-Richter model, and hazard calculations are led applying e.g. a weight of 80% on the characteristic branch with respect to 20% to the exponential Gutenberg-Richter, final hazard values decrease (Fig. 15a). Mean PGA values at 475 yrs decrease all over the country and the 0.3 g contour is now restricted to a strip of 20-40 km along the main strand of the Levant fault. Hazard values are still larger than 0.2 g in most parts of the country. If instead, a larger weight is applied to the exponential model (80%), with respect to the characteristic model (20%), mean hazard values increase, the whole country is within the contour 0.3 g, except for the regions around Tyr and north east of Baalbeck, and hazard values exceed 0.5 g along the main strand of the Levant fault (Fig. 15b).

Geodetic versus geologic slip rates on the Missyaf segment
Whereas slip rates inferred from tectonic and geodetic studies are generally consistent for the Jordan Valley fault (e.g. Lefevre et al. 2018), they differ significantly on the northern part of the Levant fault, along the Missyaf segment. Meghraoui et al. (2003) Fig. 15 Mean seismic hazard maps, PGA at 475 years return period, impact of different decisions on the weights associated to the earthquake recurrence model applied to faults. (a) 80% weight on the characteristic Youngs and Coppersmith (1985) and 20% weight on the exponential Gutenberg-Richter model branch.
(b) 20% weight on the characteristic model branch and 80% weight on the exponential model branch. B Beirut; Ba Baalbeck, Ty Tyr, T Tripoli, Z Zahle, S Saida estimated a slip rate of 6.8-7.0 mm/year along this segment, based on the occurrence of three large earthquakes in the past 2000 years identified with archeoseismicity and trenching. The complementary studies by Sbeinati et al. (2010) extended the time window to the last 3500 years, including a fourth event with a very uncertain age (1500 to 900 BC). Assuming that this ancient event corresponds to an historical earthquake in 1365 BC, they propose a range of 4.9-5.5 mm/year for the slip rate (Table 1). However, geodetically derived slip rates obtained for the northern part of the Levant fault are much lower. Alchalbi et al. (2010) derive an average value of 2.5 mm/year (1.8 to 3.3 mm/year), associated with a locking depth of 10 km, from GPS sites measured in 2000 and 2008. More recently, Gomez et al. (2020) estimated 2.2 ± 0.5 mm/year, associated to a locking depth of 10 km, using a block model with GPS velocities that rely on observations from ~ 2001 to ~ 2012. Such a discrepancy between the geologic slip rate (representing an average over a long period of time) and the geodetic slip rate (representing the current state of the fault) would need to be better understood. Both geodetic and geologic slip rates should eventually be included in the source model logic tree, as they represent the epistemic uncertainty on the amount of the deformation that accumulates along the fault and that can be released in future earthquakes. Figure 16a displays the mean hazard map at 475 years obtained applying a slip rate range of 1.7-2.7 mm/ year on the Missyaf segment. As expected, lower hazard levels are found with respect to the hazard calculation that relies on paleoseismic slip rates, with PGA values lower than 0.4 g for sites along the fault trace.  Fig. 14a: (a) the slip rate range applied to Missyaf segment is now inferred from geodetic studies (1.7-2.7 mm/year); (b) the main strand of the Levant fault is considered as one unique long segment with slip rate range 4 to 5 mm/year

Considering only one long fault for the primary strand of the Levant fault
In our source model, the primary strand of the Levant fault is constituted of three faults, from south to north, the Jordan Valley fault, the Yammouneh fault and the Missyaf fault. The Yammouneh fault stops south of the Hula Basin (Fig. 2, Table 1). This segmentation makes the strong assumption that the Hula Basin and the El Boqueaa pull-apart Basin to the north are major obstacles to a rupture propagation. However, several archeoseismic and paleoseismic studies on the Jordan Gorge section, the northern segment of the Jordan Valley fault, have demonstrated that this segment could rupture during earthquakes originating from the north. One example is the A.D. 1202 event, centered on the Yammouneh fault, that ruptured to the south through this segment (Marco et al. 2005;Wechsler et al. 2014). To understand how this segmentation impacts hazard results, we perform a test where the segmentation is relaxed and a unique long segment is considered that encompasses all segments of the primary strand of the Levant fault. In the hazard calculations, earthquake ruptures are now distributed uniformly along this ~ 450 kmlong segment and we consider a unique slip rate interval of 4 to 5 mm/year. We consider a range of 7.5-7.9 for M max (same range as Yammouneh and JVF faults, Table 1). The resulting mean seismic hazard map is displayed in Fig. 16b (PGA, 475 years return period). With respect to Fig. 13a, hazard levels slightly decrease along the Yammouneh fault, and slightly increase south of latitude ~ 33.3˚. The largest impact is obtained for sites along the Missyaf segment with much lower values considering one unique long segment for the Levant fault, as expected (~ 0.4 g rather than ~ 0.6 g).

Conclusions
We have built a source model for Lebanon and bordering regions, made of a smoothedseismicity model (off-fault forecast), combined with a fault model. The fault model relies on the findings of geologic, paleoseismic and geodetic studies in terms of fault traces, fault geometries and the quantification of the deformation that accumulates along faults. Maximum magnitudes on faults are inferred from the magnitude range of historical events associated to the faults, as well as from the application of scaling relationships to lengths and widths of faults. Two sets of alternative moment-balanced earthquake recurrence models are derived for faults, assuming either an exponential Gutenberg-Richter model or a characteristic Youngs and Coppersmith (1985) model. The off-fault seismicity forecast relies on the Lebanese instrumental catalog. To evaluate seismic hazard levels, the source model logic tree is combined with three ground-motion models equally weighted (Akkar et al. 2014;Youngs 2014 andKotha et al. 2020), as an attempt to capture the epistemic uncertainty on the prediction of ground-motion in this area.
We consider four cities, Saida, Beirut (the capital city), Tripoli on the coast and Zahle in the Beqaa plain, where we perform an extensive sensitivity analysis to understand which source model parameters' uncertainties most impact hazard estimates. This analysis demonstrates that: • The choice of the recurrence model on faults, either exponential or characteristic, impacts the most the hazard levels (up to 100% increase); • The uncertainties on the slip rate and on M max also impact strongly the hazard (up to 50% increase); • The minimum magnitude attributed to faults, a parameter that is poorly constrained, has a lower influence on hazard, but still significant (up to 30% increase of the hazard values, reducing M min from 6.0 to 5.5) • The impact of fault geometry uncertainties depends on the site, it is more significant for sites along the strike-slip faults than for sites close to Mount Lebanon thrust; • The uncertainty on the b-value for faults has a negligible impact; • The uncertainty on the smoothed seismicity model (derived only from the short instrumental GRAL catalog, 2006GRAL catalog, -2019 has no impact on hazard at these four cities; the disaggregation studies show that for most sites in Lebanon, the off-fault model has a negligible contribution to hazard. We compute mean seismic hazard maps based on the source model logic tree obtained (branches equally weighted), including only the uncertainties that control hazard variability. Our calculations lead to mean PGA values at 475 years larger than 0.3 g over twothird of the country, including the coastal region between Saida and Tripoli. Largest hazard values (> 0.5 g) are obtained along the main strand of the Levant fault. The variability obtained for hazard in the capital city Beirut is quite large, the PGA at 475 years return period varies between 0.16 and 0.45 g (16th and 84th percentiles) with a mean around 0.3 g. Disaggregation studies show that the hazard estimation in Beirut is mostly controlled by the Mount Lebanon thrust fault located beneath the city. Future studies aimed at characterizing this off-shore fault will be crucial for the hazard assessment of the city as well as within the coastal region.
Our source model relies on annual earthquake frequencies inferred from annual slip rates estimated from geologic, paleoseismologic and geodetic studies. Meanwhile, we ignore part of the information delivered by these studies. Along the Levant fault system, geologic and historical studies have demonstrated that large earthquakes occur in clusters (e.g. the sequence that lasted between the eleventh and thirteen century on the main strand, including the 1033, 1170 and 1202 events). It is not clear yet how this clustering of large events could be accounted for in seismic hazard assessment, but time-dependent models that account for fault interactions are required.
This work delivers a series of alternative source models for Lebanon, with the understanding of which parameters' uncertainties or decisions control hazard values, as well as which sources contribute for given sites, depending on the spectral period and return period. A range for seismic hazard values is obtained, for selected Lebanese cities and at the country scale. The hazard levels obtained in Beirut are compared to the output of the regional EMME project. These results may constitute a basis for discussion and for taking decisions in the perspective of future updates of the Lebanese building code.