Evaluating earthquake-induced rockfall hazard by investigating past rockfall events : the case of Qiryat-Shemona adjacent to the Dead Sea Transform , northern Israel

We evaluate rockfall hazard for the town of Qiryat-Shemona, northern Israel, situated alongside the Dead Sea 10 Transform, at the foot of the Ramim escarpment. Boulders of 1 m3 to 125 m3 are scattered on the slope above town, while historical aerial photos reveal that before town establishment, numerous boulders had reached the town premises. For the hazard analysis we first mapped the rockfalls, their source and their downslope final stop-sites, and compiled the boulder size distribution. We then simulated the probable future rockfall trajectories using the field observed data to calibrate the simulation software by comparing simulated vs mapped boulders stop-sites along selected slopes while adjusting model input parameters 15 for best fit. The analysis identified areas of high rockfall hazard at the south-western quarters of the town and also indicates that in the studied slopes, falling blocks would stop after several tens of meters where the slope angle is below 10°. OSL age determination of several past rockfall events in the study area suggests that these rockfalls were triggered by large (M >6) historical earthquakes. Nevertheless, not all large historical earthquakes triggered rockfalls. Simulations show that downslope reach of the blocks is not significantly affected by the magnitude of seismic acceleration. Considering the size 20 distribution of the past rockfalls in the study area and the reoccurrence time of large earthquakes in the region, the probability to be affected by a destructive rockfall within a 50 year time-window is of less than 5%. 25 Nat. Hazards Earth Syst. Sci. Discuss., https://doi.org/10.5194/nhess-2018-250 Manuscript under review for journal Nat. Hazards Earth Syst. Sci. Discussion started: 24 September 2018 c © Author(s) 2018. CC BY 4.0 License.


Introduction
Rockfalls are a type of fast mass movement process common in mountainous areas worldwide (Dorren, 2003;Flageollet and Weber, 1996;Mackey and Quigley, 2014;Pellicani et al., 2016;Strunden et al., 2015;Whalley, 1984).In this process, a fragment of rock detached from the bedrock along pre-existing discontinuities (e.g., bedding, fractures) slides, topples or falls along a vertical or nearly vertical cliff.Individual fragments travel downslope by bouncing and flying or by rolling on talus or debris slopes (Crosta and Agliardi, 2004;Cruden and Varnes, 1996;Varnes, 1978;Whalley, 1984, Wei et al., 2014).They travel at speeds of a few to tens of meters per second, and range in volume up to thousands of cubic meters.Due to their high mobility, and despite their sometimes small size, rockfalls are particularly destructive mass movements, and in several areas they represent the primary cause of landslide fatalities (Evans, 1997;Evans and Hungr, 1993;Guzzetti, 2000;Keefer, 2002;Guzzetti et al., 2003;Guzzetti et al., 2005).
The town of Qiryat-Shemona (northern Israel; population 25,000) is located in the northern Hula Valley (Fig. 1), one of a series of an extensional basins developed along the active left lateral fault system of the Dead Sea Transform (DST) (Freund, 1965;Garfunkel, 1981;Quennell, 1958).The town is built at the foot of the fault-controlled Ramim escarpment, which rises 800 meters above the west part of the town.New quarters of the town are being planned and built below the escarpment and up the slopes above the town.These slopes are dotted with cliff-derived boulders with measured volumes of more than 100 m 3 , which have apparently traveled down the slope by rockfall mechanism (Fig. 2).Pre-town aerial photos reveal additional rock blocks with similarly estimated volume-range within the now built town premises.Thus, the field observations and aerial photos interpretation suggest that the western neighborhoods of Qiryat-Shemona, located at the escarpment base, are subjected to rockfall hazard.
The current study aims to evaluate the rockfall hazard for the town of Qiryat-Shemona by (a) studying the past rockfall events, (b) identifying the rockfall triggering conditions and their reoccurrence interval, (c) computer-modeling the most probable down-slope rockfall trajectories which outline the hazard-prone area, and (d) by calculating the expected kinetic energy of the blocks at the town borders.Particular attention is given to the calibration of the computer modeling using mapping of past rockfall events and extracting geometrical and mechanical parameters needed for the simulations from the field observations.The age of the studied past rockfall events and their possible relation to known historical earthquakes are constrained by optically stimulated luminescence (OSL) dating (Wintle, 2008).

Study Area
The rock sequence outcrops in the lower part of the slopes, west of the town (hereafter 'the study area'; Fig. 1b) consists of Lower Cretaceous rocks (Glikson, 1966;Kafri, 1991;): The sandstone of Hatira Formation outcrops at the base of the slope, overlain by limestone and marl of Nabi Said Formation.Further up-slope, about 350 m above the town, outcrops the biomicritic limestone of Ein-El-Assad Formation, creating a 40 m high sub-vertical cliff (Fig. 2).This cliff is the source for rockfalls in the study area (see below).The Ein El Assad Formation is overlaid by a ~700 m Lower to Upper Cretaceous carbonate rocks  (Sneh and Weinberger, 2003a, b).Colluvium and rock-mass movement deposits were mapped on the slopes near Qiryat-Shemona (Shtober-Zisu, 2006, Sneh andWeinberger, 2003a).The slopes are generally covered with up to a few meters of soil.
The studied area is located along a primary fault zone of the DST (Weinberger et al., 2009).Here the western-border fault of the Hula basin branches into several faults towards the north.

Rock block inventory
A catalog of the past-rockfall derived boulders was constructed from two data sources: 76 blocks mapped and measured in the field and 200 blocks mapped using pre-town aerial photos (dating to 1946 and 1951).58 out of the 200 blocks mapped using the aerial photos were identified and measured in the field (green rectangles in Fig. 3).These 58 blocks, which were identified both on the aerial photos and measured in the field, were used to fit a correlation curve between field measured and aerial photo estimated block volumes.The correlation was used for volume estimation of the blocks that were removed from the area during the establishment of the town, but were mapped on the aerial photos predating the establishment (142 blocks out of 200).
In summary, the catalog host a total of 218 boulders, which were mapped and their volumes were estimated.This rock block inventory is the basis for the prediction of probabilities for different block sizes for the calculation of rockfall hazard and its mitigation.

Rockfall simulations
The down-slope trajectory of a rock-block (or the energy dissipated as it travels) is affected by slope geometry and surface material properties and by the rock-block geometry and material properties (Agliardi and Crosta, 2003;Guzzetti et al., 2002;Guzzetti et al., 2004;Guzzetti et al., 2003;Jones et al., 2000;Pfeiffer and Bowen, 1989;Ritchie, 1963).Parameters that quantify these measures are used as input for computer-simulation of rockfall trajectories.Several computer programs have been developed and tested to simulate rockfall trajectories: Guzzetti et al. (2002); Dorren, 2003 and references therein;Giani et al., (2004).The current study uses the 2D Colorado Rockfall Simulation Program, CRSP, v4 (Jones et al., 2000) to analyze two significant aspects of rockfall hazard in the studied area: First, the expected travel distance of rock-blocks along the studied slopes, which signifies the urban area prone to rockfall hazard.Second, the statistical distributions of block travel velocities and kinetic energy, which serves as an input for engineering hazard reduction measures.For the current analysis the model input parameters are the topographic profile of the slope (extracted from 5 m elevation contours GIS database and verified in the field), surface roughness (S), slope rebound and friction characteristics, (R n : normal coefficient of restitution; R t : tangential coefficient of frictional resistance) and block morphology.S was measured in the field according to Jones et al., (2000) and Pfeiffer and Bowen, (1989), where R n and R t were estimated via a calibration process (see below).
The CRSP algorithm simulates rockfall as a series of rock-block bounces, and calculates the changes in the block velocity after each impact with the slope surface, taking into consideration the rock and slope geometric and mechanical properties.Model output is a statistical distribution of velocity, kinetic energy and bounce height along the downslope trajectory, including stopping distances of the blocks (Jones et al., 2000).
For the current study, we simulated the rockfalls characteristics along topographic profiles extending from the Ein-El-Assad Formation, identified as the source for the rockfalls, downslope towards the town.25 topographic profiles covering the study area were extracted for the rockfall hazard analysis (Fig. 4), with high spatial density (30-100 m intervals) where the source for rock-blocks is exposed above the town and lower spatial density (150-500 m intervals) further southwards.A single simulation run (along each profile) modeled 100 rock-blocks, thus allowing the statistical analysis (Jones et al., 2000).CRSP results for each profile were later integrated spatially to compile rockfall hazard maps and other hazard properties as detailed below.

CRSP calibration
The first step in hazard analysis using a computerized model is calibration of the model input parameter.Following Katz et al. (2011), calibration was performed by comparing calculated traveling distance of rock blocks of a given size to field observed ones, while adjusting the assigned model-parameters until best-fit was obtained, i.e., back-analysis.
In the current work, CRSP calibration using back-analysis was performed along four slopes with a high number of field mapped rock-blocks (pink lines in Fig. 3).As an index for calibration quality, we used the difference between the field-observed downslope maximal travel distance along a selected slope and the simulated maximal travel distance along this slope (hereafter ∆MD in meters), for a given block-size bin.We considered the model parameters as calibrated when ∆MD = ±60 m (about 10% of average profile length).80 simulation runs, modeling the largest blocks with diameters of D of 5.8 m and 6.2 m along the four profiles were used for calibration (Determined S value is 0.5 for D=5.8 and 6.2 m).These simulations resulted in the following coefficient value ranges: Rn =0.2-0.25;R t =0.7-0.8, which are in agreement with suggested values for bedrock or firm soil slopes according to Jones et al. (2000).These values were further revised and refined following the initial velocity sensitivity analysis (detailed in the following).
The predicted seismic peak ground acceleration (PGA) for the studied area is 0.26 g (Shapira, 2002).Assuming a PGA (a) of 0.3 g with frequency (f) of 1 Hz (Scholz, 2002), the calculated initial horizontal velocity (V x ) of the rock block is 3 m/s (V x = a/f).Sensitivity analysis was performed for two end members of V x of 0 m/s (simulating a-seismic triggering of a rockfall) and 3 m/s (simulating seismic triggering).Where, R n = 0.12 yielded ∆MD = -90 m and -80 m for V x = 0 m/s and V x = 3 m/s, respectively; and R n = 0.25 yielded ∆MD = +160 m and +150 m for V x = 0 m/s and V x = 3 m/s, respectively.Thus, we infer that initial velocity has no significant effect on travel distance.An exponential regression curve was fitted for ∆MD vs. R n values (at V x = 0 m/s), which yields ∆MD = 0 m at R n = 0.22.Thus, calibration is satisfied for R n = 0.22.CRSP output is less sensitive to changes in the tangential coefficient R t in comparison to changes in the normal coefficient R n .Hence R t value was determined to 0.7 using our calibration.This value is also recommended by Jones et al (2000) for firm soil slopes.
Further simulation runs along the four calibration profiles were performed for all block sizes (D=1.3-6.2 m), using the best-fit coefficient values: S=0.1, 0.3, 0.4 for block diameters D=1.3, 2.7, 4.6 m respectively, and 0.5 for D=5.8 and 6.2 m (measured in the field per block diameter); R n =0.22; R t =0.70.All slope cells were given the same values to maintain model simplicity.
The fit between observation and simulation is plotted in Fig. 5.These results can be divided into two behavior patterns: (a) mid-size and large blocks (D ≥3 m; green, orange and red circles): the observed and simulated results are close to the 1:1 ratio; for large blocks (D ≥4 m), simulated travel distance is a little longer, which yields a more conservative result.(b) small blocks (D< 3 m): some blocks are close to the 1:1 ratio, while others demonstrate significantly longer observed travel distances.This longer observed than simulated travel distance of smaller blocks (D< 3m) may be explained in a few ways: First, smaller blocks may be more subject to creep, being more affected by water runoff and slope material movement due to their lower weight.Therefore they may travel further down after the rockfall event took place.Another possible interpretation is that the construction of town has created a different topographical setting than the slope at the time of rockfall events in the past.To circumvent this discrepancy for the hazard analysis, we use only the larger blocks.
Sensitivity analysis for block shape resulted in an insignificant difference between simulations done with sphere, disc or cylinder rock-block shapes.Accordingly, we used sphere shape rock-blocks in the prediction simulations because they yield maximum volume for a given radius and thus tend toward a worst-case scenario analysis (Giani et al., 2004;Jones et al., 2000).

OSL age determinations
For OSL age determinations of rockfall events, colluvium or soil material from immediately underneath the rock blocks was sampled.This approach constrains the time since last exposure to sunlight before burial under the blocks (following Becker and Davenport, 2003).For sampling we excavated a ditch alongside the rock block to reach the contact with the underlying soil using a backhoe, then manually excavated horizontally under the block and sampled the soil below its center.The Sampling of soil was performed under a cover to prevent sunlight exposure of the soil samples.A complementary sediment sample was taken from each OSL sample location for dose rate measurements.Locations of sampled blocks are marked in Fig. 3. Rockfall OSL age determination was based on the assumption that the sampled blocks did not creep or remove from their initial falling location.Thus, only very large blocks between 8 and 80 m 3 , weighing tens to hundreds of tons, were sampled.OSL equivalent dose De was obtained using the single aliquot regeneration (SAR) dose protocol, using a range of preheats (220 -260°C) and a cutheat 20° below preheat (Murray and Wintle 2006).The gamma and cosmic dose rates were measured in the field using a calibrated gamma scintillator.Alpha and beta dose rates were calculated from the concentration of U, Th, and K in the complementary sediment sample.

Size distribution of rockfalls
Rock-blocks, a result of rockfall events, are commonly observed along the slope west of Qiryat Shemona, at the foot of the Ein-El-Assad Formation.Their volume varies, from the smallest pebbles to boulders tens of cubic meters in volume.In places the blocks form grain-supported piles, revealing impact deformations on their common faces such as chipped corners and imbricated blocks separated along previous fracturing surfaces.To determine rockfall hazard and risk, information on the frequency-volume statistics of individual rockfalls is necessary (Guzzetti et al., 2004;Guzzetti et al., 2003).
We used the field mapped blocks to determine their volume distribution.In total, we consider this field catalog complete for block size >1 m 3 and consists of 76 blocks ranging in volume up to 125 m 3 (mode=56.25m 3 ).Following Malamud et al. (2004), the volume distribution of the mapped blocks can be determined using the probability density, p, of a given block volume Eq. ( 1): where N is the total number of blocks, dN is the number of blocks with volume between V and V+dV, and α is the scaling exponent.Our results show that the volume of the individual rock blocks from the studied area exhibits a distinct negative power law behavior, with a scaling exponent of the right tail of α = -1.17(Fig. 6).This conforms to what was found by others who examined natural rockfalls with observed α ranging: -1.07 --1.4,e.g., Guzzetti et al. (2003) Malamud et al. (2004) Brunetti et al. (2009).The scaling exponent is also similar to the value α = -1.13obtained experimentally by Katz and Aharonov (2006), while Katz et al. (2011) found a larger scaling exponent, α = -1.8.
In accordance, the probability density function (PDF) for rockfall volume (p) may be presented as a power law of the form (Dussauge-Peisser et al., 2002;Dussauge et al., 2003;Guzzetti et al., 2003;Malamud et al., 2004) Eq. ( 2): where V is the given block volume.The power law is -1.17 and R 2 = 0.72 (for the 76 field mapped blocks plotted in Fig. 3).
To simplify the hazard evaluation, the block volumes were binned into size scales of 10, 50, 100, 125 m 3 , with corresponding block diameters of 2.7, 4.6, 5.8 and 6.2 m, respectively (assuming spherical block geometry).Cumulative frequencies were used to derive cumulative probabilities for each block size (Table 1).The probability values per block diameter (Table 1) were fitted a regression curve (R 2 = 0.97), yielding the probability (p D ) for a block of given diameter (D) or smaller following Eq.( 3):

Simulations of block trajectories
For the hazard analysis, we ran computer simulations along 25 profiles including the four profiles used for the calibration (Fig. 3), using the calibrated parameters as the model input.A total of 100 computer runs, each run simulating the fall of 100 individual blocks, were performed (four runs on each of the 25 profiles: using block dimeters of 2.7, 4.6, 5.8, and 6.2 m separately, totaling 10,000 simulated single block trajectories).These results were used to analyze the hazard in the study area.

Stop angle and stop swath
The 'x% stop angle' is defined as the slope angle of the profile cell at which cumulated x% of simulated blocks stop and, in accordance, the 'stop swath' is defined as the distance (m) along the profile cell with the corresponding to stop angle that the simulated blocks covered until stopped.For a 100% stop angle, all traveling blocks will stop within a stop swath distance.The 50% and 100% stop angle and stop swath data were extracted from the CRSP simulation analysis (Fig. 7).100% stop angles for all profiles (red circles in Fig. 7) have a mean of 7.7° with SD = 2.3° (range of 3°-12°); 50% stop angles (blue triangles) have a mean of 10° with SD = 5.3°.All other cell slope angles (gray circles) vary widely between 7°-88°, among them very few are less than 10°.Stop swath distances range between 8 m to 105 m, having a mean of 38 m with SD = 24 m.Only in two profiles (out of 25) did the stop swath distance exceed 65 m.In both these cases, 100% stop angle is steeper than in most other profiles (10°-11°).No significant correlation was found between 100% stop angle and stop swath distance.

Rockfall hazard
Rockfall hazard map for Qiryat-Shemona is presented in Fig. 8.The hazard map was compiled from the simulated maximal travel distance (where 100% of blocks stop) of the largest blocks (D >4.6 m, V >50 m 3 ) with the probability of occurrence, pD = 11% (Eq.3).The calculated block trajectories cross the town border and mark the town premises that are subject to rockfall hazard along 8 out of 25 simulated profiles (#8-#14 and #16, marked by † in Fig. 9).The area subjected to rockfall hazard is about 1.55 km 2 , currently including several houses (according to the last updated google Earth image from Nov 2014).For D=4.6 m, block impact velocity varies between 9.5-13.7 m/s and kinetic energy between 7,400-16,300 kJ (Table 2).CRSP simulated maximal travel distance and CRSP velocity and kinetic energy analysis points at town border impact locations are plotted in Fig. 9.For details of the kinetic analysis see locations marked by profile indexes and refer to Table 2.

OSL age determination of rockfalls
OSL ages were determined for eight rock blocks.The locations of these blocks is marked in red circles in Fig. 3.These ages range from 0.7 to 8.1 ka, with uncertainties of 9-40% (Table 3).

Triggering conditions for rockfalls in the studied area
We interpret the field observed grain supported structure of aggregations of blocks of various sizes, with impact deformations (e.g., chipping) on their common faces as evidence for catastrophic events, involving numerous blocks.Long-term erosion which results in single sporadic block failures would have resulted in matrix-supported blocks and not in the evidence observed here.We conclude that the rockfalls were mainly triggered by discrete catastrophic events such as earthquakes or extreme precipitation events.The question of a triggering mechanism in the case of a catastrophic rockfall event is an important one when attempting to evaluate the temporal aspect of rockfall hazard.The recurrence time of an extreme winter storm or a large earthquake may give some constraints on the expected recurrence time of rain-induced or an earthquake-induced rockfall, respectively.Furthermore, it might suggest a periodical probability for the next rockfall to occur when hazard is calculated.
The correlation of rockfall events to historical extreme rainstorms events is limited due to the lack of long enough historical rainstorm record.Noticeably, even following the extremely rainy winters of 1968/69 and 1991/92, in which annual precipitation in northern Israel was double than the mean (IMS, 2007), no significant rock-mass movements and rockfalls are reported in the study area.Thus, rainstorms are ruled out as a favorable triggering mechanism.A possible correlation between the dated rockfall events and historical earthquakes is analyzed below.

Non-random temporal distribution of rockfalls and correlation to earthquakes
The wide range of OSL ages, between 700 and 8100 years before present (Fig. 10 and Table 3), rules out the possibility of a single rockfall event.Given the rich historical earthquake record in the vicinity of the studied area, the positive correlation between rockfall events and historical earthquakes may shed light on the triggering mechanism of the rockfalls.A similar approach was used by Matmon et al. (2005), Rinat et al. (2014) and Siman-Tov (2009), who dated rockfall events ~ 30 km SW from the studied area, where he found a positive correlation between rockfall events and historical earthquakes, dated 749 AD and 1202 AD.To perform this correlation, we overlaid the 8 OSL ages with a set of 9 large historic earthquakes that occurred in the vicinity of the study area (Table 4, Fig. 10).This set contains earthquakes that: (a) occurred within the time spans of the OSL ages; (b) maximum estimated intensity is at least 'IX' on an EMS macroseismic local intensity scale and/or their estimated moment-magnitude is 6 or larger; (c) the distance between the study area and affected localities reported does not exceed 100 km (following Keefer, 1984).
The 8 OSL ages (Table 3; Fig. 10) span over the past 8000 years with a mean of one age per 1,000 years.The validation of OSL age clustering was obtained performing a binomial distribution test, which gives the discrete probability distribution P(k,p,n) of obtaining exactly k successes out of n trials.The result of each trial is true (success) or false (failure), given the probability for success (p) or failure (1-p) in a single trial.The binomial distribution is therefore given by Eq. ( 4): The selected nine historical earthquakes, each with a ±50 year time window, span over 900 years out of the given 2,800 years period (each event ±50 year time window spans 100 years).Therefore, the probability p for a single random earthquake to occur within this period is p = 900/2800=0.32.The number of trials n is the number of earthquakes n = 9.In five cases (199 BC, 502 AD, 551 AD, 1033 AD, 1202 AD), a success (match between an earthquake and an OSL age) is obtained, therefore the number of successes is k = 5.This fit between OSL ages and earthquakes is detailed in Fig. 11.Accordingly, the binomial distribution is P(k, p, n)= P(5, 0.32, 9) = 0.09, i.e., there is a 9% probability to obtain such a distribution of events in time randomly.Hence, we suggest that the OSL age distribution is significantly clustered around dates of the discussed historical earthquakes, with 91% confidence level, thus suggesting a likelihood of seismic triggering for the rockfalls in the studied area.
Assuming that magnitude six earthquake or larger is needed for rockfall triggering (following Keefer, 1984) and based on the recurrence time of 550 years given for these earthquakes (Hamiel et al., 2009), we predicted a ~550 years recurrence time for rockfalls in the studied area.
Not all historic earthquakes are represented in our OSL data set, such as the 1759 AD (Ambraseys and Barazangi, 1989;Marco et al., 2005) and 1837 AD (Ambraseys, 1997;Nemer and Meghraoui, 2006) earthquakes, which both induced extensive damage to cities not far from the studied area (Katz and Crouvi, 2007).This lack of evidence could be explained by OSL undersampling, or because earthquakes only trigger rockfalls that were on their verge of stability (Siman-Tov et al., 2017).
Based on the above analysis we correlate the past rockfalls to historic earthquakes as follows: a. QS-4 (0.89±0.16 ka) fits the historical earthquakes of 1202 AD or 1033 AD.Since the 1202 AD earthquake accounts for severe damage in other places in northern Israel (Marco et al., 1997;Wechsler et al., 2006), we find it a better candidate for triggering a rockfall event than the 1033 AD event.
b. QS-3 (1.5±0.13 ka) fits the historical earthquakes of 502 AD and 551 AD.The 551 AD earthquake is reported at more localities along the DST than the 502 AD (reported on shoreline localities only).Hence we find it a better rockfall triggering candidate.
e. QS-5 (4.0±0.70 ka) and QS-12 (4.1±09 ka) fit the 2050/2100 BC earthquake (or two separate events) suggested by Migowski et al. (2004).This also fits the findings of Katz et al. (2011) andYagoda-Biran et al. (2010), who found evidence for earthquake and earthquake-induced slope failure east of the Sea of Galilee (ca.50 km south from study area, along the DST) with OSL ages of 5.0±0.3 and 5.2±0.4ka, respectively, suggesting the area had experienced one or more strong earthquakes.f.QS-1 is (7.0±1.1 ka) may correlate with the Mw=7 earthquake derived from slope movement evidence dated by Yagoda-Biran et al. ( 2010) to 6.0±0.4 ka east of the Sea of Galilee.

5.3
Rockfall hazard for the town of Qiryat-Shemona

Area subject to rockfall hazard
The nature of the analyzed past rockfall events in the studied area can be used to constrain the possible characteristics of the expected future rockfall events and promote hazard mitigation.The predicted probabilities P D for specific rock fall with a given block diameter or smaller, derived from the regression curve (Eq.3), are presented in Table 1.P D (2.7), the cumulative probability for a block of D = 2.7 m (V = 10 m 3 ) or smaller is 0.67.Consequently, the probability for traveling blocks of 2.7< D< 6.2 m (10< V< 125 m 3 ), is 1-P D (2.7) = 0.33 or 33%.The occurrence of larger, more destructive blocks amongst these (D = 4.6 -6.2 m or respectively, V = 50 -125 m 3 ) is P D = 11%.Despite their lower probability, these blocks would reach the farthest distances; and hence pose the largest hazard to town.
The area subjected to rockfall hazard for travel distances of the large blocks (D = 4.6-6.2m, V >50-125 m 3 , 11% probability) appears in Fig. 8.About 50,000 m 2 (0.05 km 2 ) of the westernmost, currently inhabited and built, urban area (out of 1.55 km 2 total area mapped as subject to rockfall hazard, which also includes the above town slope outside the urbanized area) is mapped under direct rockfall impact hazard.This hazard mapping may be used to plan mitigation actions and also as a basis for future urban planning.
We note that the main road connection the town southwards, which can serve as evacuation route, is not under high rockfall hazard.We also note that some smaller blocks (D≤ 3 m) were mapped from the historical aerial photos predating town establishment further down slope below the simulated 100% blocks stop (Fig. 3; blue circles in Fig. 6).As suggested above these might be blocks that traveled farther downslope by creeping after the rockfall event.Another possible explanation is that these now inhabited parts of the slope were altered and even levelled by the construction works.Hence the simulated profile extracted from current topography is different from the slope topography on which these smaller blocks traveled before the construction of town.No detailed topography maps pre-dating town establishment is available to verify this.
The stop angle results (Fig. 7) indicate that most blocks (>50%) keep traveling down-slope until the slope angle decreases to 10°-15°.All blocks stop where the slope angle decreases to 5.5°-10.0°.Combined with the stop swath distance results (Sect.

Rockfall hazard probability
We can calculate the probability of a rockfall occurrence P EQ in the next 50 years, assuming magnitude M=6 as the threshold for rockfall: 50/550 = ~0.09or 9%.Considering the probabilities for given block sizes detailed above, the probability for rockfall hazard per specific block size (H R ) may be predicted as Eq. ( 5): where P D is the cumulative probability per block diameter D (Table 1) and P EQ is the rockfall occurrence probability calculated above to be 9%.Accordingly, predicted H R for the next 50 years for D between 2.7 and 6.2 m is H R ~3% and for larger blocks only, D between 4.6 and 6.2 m is H R ~1%.We do not present a time-dependent earthquake recurrence interval calculation because the time passed since the last large earthquake is not well constrained.

Conclusions
In this work, we studied rockfall hazard for the town of Qiryat Shemona to demonstrate computer-simulation based hazard evaluation in cases where the study area is residing below a cliff in an apriori exposure to rockfall hazard, but no documentation of recent or past rockfalls is available.To overcome this lack of observations, we derived the needed geometrical and mechanical parameters for the computer hazard analysis from field study of past events.We located areas subject to current rockfall hazard at the south-western quarters of town.
OSL age determination of several past rockfall events in the study area suggests that these rockfalls were triggered by large (M>6) historical earthquakes and in accordance, the estimated rockfall recurrence interval is hundreds of years.Nevertheless, we conclude that not all historical large earthquakes triggered rockfalls.Additionally, we infer that the downslope travel distance of the blocks is not significantly affected by the magnitude of seismic accelerations.
The field calibrated simulation results indicate that in the studied slope, falling blocks would stop after covering a distance of several tens of meters once the slope angle decreases to around 5°-10°.Considering the size distribution of the past rockfalls in the study area and the reoccurrence time of large earthquakes, the probability to be affected by a destructive rockfall within a 50 years time-window is of less than 5%.

Figure 8
Rockfall hazard map of the study area.The area subject to rockfall hazard is defined from the source escarpment to 100% stop line.
Map compiled from maximal travel distance of 25 rockfall simulation profiles performed using CRSP (green lines in Fig. 4).

Figure 9
Rockfall hazard map for the town of Qiryat Shemona.Yellow line represents the CRSP 100% stop line calculated for large blocks (D is 5.8 m and 6.2 m).Orange-black hazard triangles mark simulated stop line and town border impact location for each profile; profile numbers in yellow refer to simulated profile numbers; indices with sword label ' †' mark locations of rockfall impact at town border; inset: black rectangle outlines enlarged map area.
Figure 10 OSL ages plotted in rank order and suggested rockfall triggering earthquakes.OSL age results for the past 8000 years in black circles with error bars; selected earthquakes suggested as rockfall triggers using data from earthquake catalogs and paleoseismic data in vertical gray lines, chronologically labeled at top axis; see text for details.

4. 2
.1 above), considering the means and SD's, CRSP results indicate that falling blocks would stop after covering a distance of 14-62 m on a slope angle of 5°-10° at the point of stopping.This conclusion may help when considering rockfall mitigation design.Nat.Hazards Earth Syst.Sci.Discuss., https://doi.org/10.5194/nhess-2018-250Manuscript under review for journal Nat.Hazards Earth Syst.Sci. Discussion started: 24 September 2018 c Author(s) 2018.CC BY 4.0 License.

Figure 1 (
Figure 1 (a) Location map of the study area on the backdrop of major faults of the DST; Inset shows the plate tectonic setting of the DST; (b) Orthophoto map of the study area (black rectangle).The rockfall source, Ein-El-Assad Formation (EEA), is marked by a blue line; town border in red dashed line.

Figure 2 A
Figure 2 A block at its stop site and the source for rockfalls: the cliff of Ein El Assad Formation at the back.The large boulder, ~41 m3, is block no 016 (Table3), under which OSL sample QS-11 was excavated.

Figure 4
Figure 4Location of 25 simulation profiles of rockfall trajectories.Faults traces are fromSneh and Weinberger (2003a).The source for the rock-falls (Ein-el-Assad formation) is marked with blue line.

Figure 5
Figure 5 Field observed distance from the source and maximal simulated travel distances of rock blocks.Block diameters are both size-and color-coded.The 1:1 line (x=y) is plotted in gray.

Figure 7
Figure 7Slope gradients at different stop angles.Tangential denotes simulated profile number 1 to 25. Radial denoted the slope angles between 0 and 90 of all cells along the profile.For example: the cells along profile 20 have slope angles that vary between 8°-36°.Red circles are 100% stop angles (slope angle of the profile cell at which cumulated 100% of simulated blocks stop); blue triangles are 50% stop angles; gray circles are all other cells in each profile.

Figure 11 Figures
Figure 11Clustering of OSL ages around historical earthquake dates.Five 'successes' for binomial distribution test defined using a ±50 years time-window (see text for detail); vertical lines are ±50 years time-window around earthquake dates (red solid lines and blue dashed lines are used to distinct between overlapping time windows; black circles are OSL ages with ±50 years error bars (black) and their lab reported error range (gray).

Table 4 . Selected historic earthquakes, candidates as possible rockfall triggers in the studied area*
Keefer, 1984)andidate rockfall triggering earthquakes, which contains historical earthquakes that: (a) occurred within the time spans of the OSL ages; (b) maximum estimated intensity is at least 'IX' on an EMS macroseismic local intensity scale and/or their estimated moment-magnitude from previous paleoseismic studies is 6 or larger; (c) the distance between the study area and affected localities reported does not exceed 100 km (followingKeefer, 1984).Nat.Hazards Earth Syst.Sci.Discuss., https://doi.org/10.5194/nhess-2018-250Manuscript under review for journal Nat.Hazards Earth Syst.Sci. Discussion started: 24 September 2018 c Author(s) 2018.CC BY 4.0 License.
Table 3), under which OSL sample QS-11 was excavated.Block mapping areas: green outline -field mapping area, where 76 blocks >1 m3 (green squares) were mapped; blocks sampled for OSL age determination are marked in red circles; black rectangle outline -area where blocks were mapped from 1946 and 1951 aerial photos (blocks in yellow squares).Size of squares denotes block diameter bins, see legend.Pink lines represent slope profiles on which CRSP calibration simulations were run.