Evaluation of Maximum Rockfall Dimensions Based on Probabilistic Assessment of the Penetration of the Sliding Planes into the Slope

There is intrinsic difficulty in the investigation of the largest volume of rockfalls that is expected in an area, which lies in the small number of large events, in registrable times. The maximum credible rockfall size has been associated with the properties of the rock mass discontinuities, as they delimit detachable rock blocks, and in particular with the penetration of those discontinuities that comprise rockfall sliding planes. In highly fractured rock masses, the evaluation of the penetration remains an issue. A probabilistic methodology is proposed, to measure the penetration of potential sliding planes into the interior of a rocky slope. The main hypothesis of the method is that the sliding plane persistence is interrupted along its two directions, at the intersection with two lateral discontinuity sets, as the latter displaces the former. Due to the displacement, the sliding planes are formed by quasi-planes that contain a maximum number of spacings of the intersecting joints, hence their size is restricted. The methodology requires as an input the spacing of the intersecting joint sets. Its application to a granodiorite slope confirms that for the study site, there is a maximum volume of rockfalls, excluding the possibility of large stepped failures and rupture of rock bridges. The maximum calculated persistence for the two existing sliding planes in the study site is, respectively, 28.0 m and 48.5 m. The maximum calculated sliding plane surfaces are, accordingly, 282.5 m2 and 289.3 m2. These results are compared against the observed scar dimensions at the study site, which have been retrieved alternatively, by processing a LiDAR point cloud. The results from the two alternative sources are similar, indicating that the methodology can be efficiently used to assess the sliding plane persistence and the expected maximum rockfall magnitude at the study site.


Introduction
Rockfall risk assessment and mitigation require data for the magnitude of previously or potentially mobilized rock masses from a slope. The rockfall magnitude, expressing the total volume of the released mass or the area of the source, is necessary to estimate the destructive energy of an event. For fragmental rockfalls, the number and size of blocks depend on the initial rockfall volume with a subsequent effect on their intensity and run out . The design of protection measures is usually implemented for the most likely expected scenarios and it entails a residual risk associated with big but rare rockfall events, for which full mitigation is not feasible (Corominas et al. 2005;Nicot et al. 2001). The rockfall size has additionally been used as a discriminating value to mark the transition from solid particle type downslope movement, for fragmental rockfalls, to flow type propagation for rockfall avalanches, with a subsequent effect on the motion properties and propagation distance (Evans and Hungr 1993;Hungr 2001;Corominas et al. 2018). The evaluation of the expected rockfall magnitude with a focus on the biggest expected sizes is of major importance in the rockfall analysis.
De facto methods for assessing quantitatively the rockfall hazard using frequency-magnitude relation are based on techniques that find application from statistical physics and complexity theory to natural hazards, with the prevalence of the Gutenberg-Richter power law. In the latter, the probability p(x) of an event of magnitude x or greater occurring is given by the equation p(x) ~ x -b , where b is a constant (Malamud and Turcotte 2006;Malamud et al. 2006). The use of the Gutenberg-Richter power law for the expression of the rockfall magnitude-frequency relation has been applied to a variety of rocky slopes and geological settings, in numerous and well-documented studies of natural slopes and rock cuts (Dussauge et al. 2003;Guzzetti et al. 2003). Power law relations have been indicated to fit well medium and large size landslides (Hungr et al. 1999).
For the occurrence of earthquakes, it is well-known that for large seismic moments the Gutenberg-Richter power law distribution has to be modified due to energy conservation and geometrical reasons (Sornette and Sornette 1999). This might take place in terms of a second power law for which a larger b value is applied beyond a cross-over magnitude, or based on either a "hard" or a "soft" magnitude cutoff using an exponential taper. Similarly, for rockfalls, the extrapolation of the same power law as for smaller rockfall sizes has been questioned (Corominas and Moya 2008). For large magnitudes, deviations from the power law line have been registered, demonstrating a heavy tailed, exponential behaviour, calling into question the scale invariant character of the rockfall magnitude-frequency relation (Guthrie and Evans 2004;Pelletier 1997). The characteristics of this heavy tailed behaviour have received minor attention for landslides and rockfalls compared to earthquakes (Tanyaş et al. 2019).
There is intrinsic difficulty in the investigation of the largest volume of rockfalls that is expected in an area, which lies in the small number of large events, in registrable times. Besides the large statistical uncertainties from which large events suffer, the problem of defining the tail of the rockfall frequency-magnitude distribution and of assessing potential outlier events remains an issue that few researchers have tackled so far (Brideau et al. 2009;Corominas et al. 2018).
The rockfall size has been associated with the properties of the rock mass discontinuities that bound blocks of intact rock mass, and define detachable rock masses, and in particular with their spacing (Palmstrom 2005;Lambert et al. 2012;Ferrero et al. 2011). Elmouttie and Poropat (2012) used discontinuity fracture networks to study how variations of discontinuity spacing and orientation affect the block size. Kim et al. (2007) introduced the effect of disrupted discontinuities. Rosser et al. (2013), Royán et al. (2015), and Stock et al. (2012) argued that big rockfall events take place through progressive failure of the intact rock bridges in the rock mass. Therefore, instabilities are controlled by the length of alternative paths for joint propagation, through the intact rock, and the strength of the rock material (Einstein et al. 1983). Intensively fractured rock masses, with small length rock bridges, are prone to produce numerous small size failures, leaving scars that might eventually lead to large continuous surfaces. Continuous surfaces are formed either as a result of isolated events, by coalescence, or in a single event (Mavrouli et al. 2015). On the opposite, when extended rock bridges are present, they increase the overall rock mass resistance and rockfall events are less frequent. However if the rock mass resistance is exceeded, large failures can be expected. Rock mass classifications like the Geological Strength Index (GSI) are based on the rock mass fracture degree (Palmström 2009). These studies provide a valuable input for the analysis of the rock block size, when the penetration of discontinuities in the rock mass is known. However, although the penetrability of discontinuities into the slope has been identified in several works as a key factor for the rockfall size , there is a certain difficulty for its direct assessment. In selected cases, ground-penetrating radar or drilling methods have been used for this purpose (Deparis et al. 2011). Still, these methods are not applicable for extensive and highly fractured rock masses.
Slope topography as intersected by the main discontinuity sets, with given properties (dip, dip direction and spacing) is the primary geometrical constraint that marks an upper limit for the credible rockfall volume. Brideau et al. (2009) after field work and numerical analysis of four sites suggested that the geometry of the rock slope failures is strongly influenced by the presence and location of tectonic structures and they discussed the role of tectonic damage and brittle rock fracture in the development of large slope failures. Corominas et al. (2018) argued that there exist geological constraints on the maximum size of failure, which is controlled by the fracture pattern, for a slope in Andorra. In that study, mutually interrupting highly persistent joint sets prevent the formation of continuous sliding surfaces and large slope failures. As therein mentioned, the volume restriction can be overcome, to some extent, either by coalescence of basal planes or through step-path failures involving the breakage of rock bridges, which in any case involve smaller volumes than in the case of fully persistent basal joints. While there is no evidence of large broken rock bridges in that study area, small stepped path failures and breakage of short rock bridges of up to 20 cm have been observed on images obtained by helicopter flights and UAV-based close-range digital photogrammetry. Corominas et al. (2018) argued that the observed step corresponds to the systematic displacement of a discontinuity set by an intersecting posterior joint set, which moved across the former. Such displacements represent common tectonic effects leading to the genesis and evolution of current fracture patterns.
The present work focuses on the same study site and builds upon the work of Corominas et al. (2018) and Mavrouli and Corominas (2017). It uses as a starting point the afore-mentioned systematic tectonic displacement of the discontinuity sets to investigate the frequency and extent of their coalescence. The discontinuity sets correspond to the predominant rock sliding planes during the rockfall release. The aim of this study is to assess probabilistically the formation of continuous and large sliding surfaces, and consequently the potential for large volume failures, or the existing constraints. Given that the systematic assessment of the discontinuity penetration into the rocky slope is practically unfeasible, we alternatively propose a statistical method for evaluating the length and width of (un)interrupted sequences of discontinuity planes. The method uses data for the discontinuity properties that have been previously acquired by measurements on a LIDAR point cloud of the area (Santana et al. 2012).

The Solà d´Andorra Study Site: Background and Data Collection
The study site is the slope above the urban area of Andorra la Vella and Santa Coloma (Fig. 1). It is a slope very active in rockfalls with an average frequency of events once every 2 years (Corominas and Moya 2010). The predominant material is highly fractured granite gneiss. Demographic pressures and low availability of space for urban expansion have led to the extension of built areas to the lowest part of the slope, with a substantial risk for properties and people. High dissipative steel fences, with a capacity exceeding 5000 kJ, have been installed. After the installation of the protection barriers, Corominas et al. (2005) calculated the residual risk for blocks with a size larger than 10 m 3 . In 2008, a block of 30 m 3 overtopped the barrier and damaged a workshop. In 2013, further blocks reached the nearby buildings . Known rockfall volumes in the area reach up to 300 m 3 (Pica talus slope, 2003). Santana et al. (2012) estimated the statistical frequency of historical rockfall volumes, assuming that the actual relief Fig. 1 (Left) the couloir of the Santa Coloma on the Solà d´Andorra; (right) four principal discontinuity sets that form the detachable blocks from the slope of the slope face has been shaped as a consequence of the release of rock masses. The main hypothesis of their work has been that the rockfall volume distribution corresponds to the missing mass from the existing scars. Through a stochastic procedure, maximum rock blocks missing from the scars were indicated to be of a few thousands of cubic metres, with the maximum surface measured for a basal plane of rupture equal to 213 m 2 . In that work, the possibility of large step-path failures was not taken into consideration. No large rupture niches have been identified on the slope by the analysis of high-resolution 3D topographic profiles (Corominas et al. 2018). Mavrouli et al. (2015) assessed the volume of large potential instabilities through an approximate GIS-based methodology, considering the interaction of the topographical surface with the discontinuity patterns serving as potential basal planes. They assumed cubic or parallelepiped rockfall shapes. Maximum volumes were indicated to be 50,000 m 3 and 25,000 m 3 , respectively. Corominas et al. (2018) used a high-resolution 3D model of the slope, acquired by UAVbased images and close-range digital photogrammetry, and they delimited on the slope face big rock spurs up to 20,200 m 3 , with a sliding surface area of 3268 m 2 . They argued that, for that slope, there exist geological constraints to the maximum credible rockfall volume. They suggested that a random distribution of large rockslides and rock avalanches is not to be expected, as the disruption of the basal planes restricts the formation of very large detachable rock masses. Given the ample background of studies in the area, and the extensive field data collection along decades, this slope was selected to explore the effect of joint disruption on the formation of large rockfalls, in quantitative terms.
Eight discontinuity sets have been identified via field data collection and scanlines, as well as by the processing of a LiDAR obtained point cloud (Santana et al. 2012). Amongst them, the joint sets F3 (dip dir/dip: 157°/56°) and F5 (182°/47°) are basal sliding planes, which intersect with the two lateral cracks F1 (54°/59°) and F7 (141°/89°), as seen in Fig. 1. The lateral cracks are almost perpendicular to the former. Along their height, they include various spacings of the basal planes. The predominant type of failure is plane failure along the basal planes F3 and F5. Lateral scar surfaces are mostly plane.
The present relief does not provide evidence of large steppath failures, as remaining extended continuous surfaces cannot be observed. Only short steps corresponding to the breakage of intact rock bridges, smaller than 20 cm have been observed (Fig. 2), where the largest rockfall scars are. The formation of the V-shaped channel of Santa Coloma has been interpreted as the consequence of gradual and independent rockfalls instead of a large event. This is supported by the fact that no evidence of historical large failure has been recorded or found by deposits in the narrow valley at the toe of the slope or in the Valira River, which crosses the valley (Corominas et al. 2018). A closer look of the discontinuities and movement indices indicates that the joint sets F3 and F5 have been displaced by the joint sets F1 and F7. An example of joint displacement is shown in Fig. 3.

Probabilistic Assessment of Penetration of Basal Planes into the Slope
The methodology presented in this section aims at investigating the persistence and penetration of basal sliding surfaces into the slope and consequently, the rockfall size constraints, in the study site. To evaluate how the displacement of the basal surfaces by the lateral joint sets affects their size in the two directions of length and width, the probability of the basal quasi-planes preserving or losing their continuity at each intersection with the lateral joints was calculated. From now on, with the term basal quasi-planes or quasiplane surfaces we refer to the coalesced surfaces, which can be, in principle, separated by small steps up to 20 cm as a result of displacement. Figure 4 shows a general representation of a simplified rock mass fracture pattern with three discontinuity sets, the basal joint set A and the lateral joint sets B and C. The threedimensional rock mass is composed by columns i = {1…n}, in two directions, and bounded by the two vertical joint sets B and C. The columns contain k = {1…m} series of planes  Rock mass profile with three discontinuity sets for the calculation of the probability mass function (PMF) of generating a continuous plane A containing i = 1…n spacings, in the direction of length or width. At their intersection with either B or C faults, the planes of set A are displaced at a distance d ik . If d ik < 20 cm. The successive sections of plane A form a single continuous basal quasi-plane of the joint set A. The continuity of the sections of plane A is controlled at each intersection with the planes B or C. The spatial extension of the basal quasi-planes depends on the amount, as well as the length or width of the individual sections, between consecutive spacings of the joint sets B and C. By definition, the length is assumed to be along the sliding direction and the width along the transversal. Through this process, the size of the surfaces of the set A can grow by merging i = {1…n} sections, if two criteria are fulfilled: (a) the surfaces are continuous and (b) the surfaces' relative height permits sliding. Two consecutive sections are considered to be continuous in length or width, if the perpendicular distance between the two is lower than a given threshold. The establishment of this threshold can be made upon field observations as more extensively discussed in Sect. 4. To guarantee the role of the surfaces of set A as possible sliding surfaces, the second condition is translated into this: starting from the exterior and moving towards the interior of the rock mass, the surfaces of set A grow only if the successive sections are progressively on a higher level than the former (Fig. 5). This constraint is applied only along the sliding direction, to eliminate quasi-planes with rock recesses that obstruct sliding.
The proposed procedure requires as an input a range of discrete spacing values for the joint sets A, B and C. It consists of the three steps described in Sects. 3.1, 3.2 and 3.3.

Generation of the Rock Mass Sample
For a given slope, the generation of the rock mass according to Fig. 4, considering the intersection of the three discontinuity sets A, B, and C takes place as following. Sample columns (sections corresponding to spacings of the intersecting joints) numbered i = {1…n}, of k = {1…m} joint series within each column are generated. Between the planes k = m-1 and k = m, there is random spacing of distance s, that takes values from the discrete spacing set S (s ϵ S) of A. The length of each section is random and takes values from the discrete spacing set C. Accordingly, the width takes values from the discrete spacing set B. The length or width of each column is independent from its neighbouring ones. The number of joint columns and series in the generated sample representing the rock mass should be sufficiently large to achieve convergence of the sample mean to the true values of the probabilities calculated according to Sect. 3.2, and thus stability of the results. Frey 2010 estimated the required sample size of Monte Carlo simulations using confidence intervals. Hahn 1972 andOberle 2015 analysed the effect of Monte Carlo iterations and the accuracy or error in the estimation of the mean of the probability distribution for the quantity of interest, using the central limit theorem and the Wald confidence interval. They proposed the calculation of the sample size in function of the maximum percentage error of the mean, when the standard deviation of the probability distribution for the quantity of interest is known. As in the literature, as far as the authors know, there are no clear recommendations for the a priori calculation of the sample size when the probability distribution is not known; an empirical approach is followed for the selection of the sample size, with the objective to capture maximum joint set persistence, under the computing capacity limitations of Excel.

Calculation of the Probability Mass Function (PMF) of a Continuous Quasi-Plane Surface Being Composed by i = {1…n} Sections Along a Direction
First, the vertical distances y ik of all the generated planes with i = {1…n} and k = {1…m} from a common reference plane are calculated. The reference plane intersects the joint set A at the point O (0, 0). This permits the identification of the minimum distance d ik between each joint of the column i and the joints of the column i + 1. The afore-mentioned criteria, (a) guarantying continuity and (b) permitting sliding, are then applied to assess the amount of spacings in the sampled rock mass, which form basal quasi-plane surfaces, along each direction. The criterion (b) is applied only for the sliding direction. The probability mass function (PMF) of a plane being constituted by i = {1…n} consecutive spacings along the sliding direction is calculated using Eq. (1): (1) P s (i) = N As (i)∕m, for i = {1, 2 … n} where N As (i) is the number of series (spacings) of the rock mass fulfilling the criteria (a) and (b) for i successive spacings, m is the total number of series k in the first column, along length.
Along the transversal direction, the respective probability is calculated using Eq. (2): where N At (i) is the number of series (spacings) of the rock mass fulfilling criterion (a) for i successive spacings, m is the total number of series k in the first column, along width.
The PMF is highly dependent on the spacing of the joint set A. When this is smaller (denser), then the probability of having more extensive A surfaces increases.

Calculation of the Total Length, Width, and Area of the Basal Quasi-Plane Surfaces
The length and width of each section are determined by the spacing of the intersecting and displacing joints B and C. The total length L A and width W A of a continuous quasiplane surface are given by Eqs. (3) and (4): where L A is the total length of quasi-plane, along the sliding direction, W A is the total width of quasi-plane, along the transversal direction, s B and s C are the spacing s ϵ S of the intersecting set B or C, respectively. Random samples of L A and W A are generated through a Monte Carlo simulation, applying Eqs. (3) and (4), for i being a variable following the PMF calculated in Sect. 3.2. The s B and s C take random values from the spacing sets B and C.
Last, random samples of the area of the quasi-plane surfaces A A are generated through a Monte Carlo simulation, using Eq. (5). For the sake of simplicity, rectangular basal areas are assumed.
where A A is the total area of basal quasi-plane for the set A.

Dimensions of Quasi-Plane Basal Surfaces in the Study Site
The procedure of Sect. 3 was applied to the study site, to assess the continuity of the basal planes that belong to the joint sets F3 and F5, and are intersected and displaced by the joint sets F1 and F7. For the planes F3, the sliding direction and length are parallel to the direction of F1. Its width is parallel to F7. Vice versa, for the plane F5, the sliding direction and length are parallel to the F7 and its width is parallel to the F1 (Fig. 4). In the study area, the intersecting joints are almost perpendicular to the basal planes, forming parallelepipeds.
The sliding surfaces present certain undulation and roughness (Santana et al. 2012). As a result, neighbouring sections may lay at a different height either because the discontinuity surface is undulated or because they belong to different parallel planes of the same discontinuity set. The undulation of the discontinuity surfaces at the study site has usually amplitudes less than 10 cm and exceptionally up to few decimetres. On the other hand, the minimum observed spacing in the field and with TLS images is 10 cm. For values between 10 and 25 cm, there exists some overlapping between undulation and spacing, thus we selected a threshold of 20 cm, to distinguish between spacing and undulation. This means that two contiguous parallel surfaces of the same discontinuity set, which are separated less than 20 cm will be considered to be in the same quasi-plane. Therefore, the value of 20 cm was selected as the threshold for the application of the criterion (a), on the continuity between successive sections. The criterion (b) was checked just along the sliding direction for F3 and F5. Spacing data were extracted from the LiDAR point cloud. For this, a representative sample of planes was identified for each joint set and the perpendicular distance between adjacent planes was measured using the software Rhinoceros ® . The spacing distributions for each joint set are shown in Fig. 6.
Following the steps of Sects. 3.1, 3.2, and 3.3, ten sample columns (sections) with 50,000 series each were generated, in an excel environment. The high number of 50,000 simulations was chosen after trials, to guarantee the stability of the calculations.
Median and standard deviation values of the randomly sampled spacings were compared with measured spacings to test their representativeness. The median value error was 0 and the maximum standard deviation error was 4%.
The vertical distance (height) of each respective joint from the reference point O (0,0) was then extracted. The height variation within a joint was not considered, hence equal height for both joint ends was assumed. Eventually, the minimum vertical distance between the closest joint ends of successive columns was calculated.
For both joint sets F3 and F5, following the steps of Sect. 3.2, the probabilities P s (i) and P t (i) of a given number of spacings i along the sliding and the transversal direction were obtained from Eqs. (1) and (2). The calculated probabilities for the 50,000 random samples were afterwards proportionally downscaled to the number of the observed F3 and F5 surfaces that were identified on the LiDAR point cloud. These are 4760 for F3 and 3920 for F5. The results are shown in Table 1. Table 1 also summarizes the proportional to the observed scars number of planes, N, that are expected to contain i spacings along the sliding and transversal direction. N 3s (i) and N 5s (i) are, accordingly, the number of planes along the sliding and N 3t (i) and N 5t (i) along the transversal direction, for the two joint sets.
The calculated probabilities, downscaled to the observed scar number, indicated, for the F3, a maximum number of spacings i equal to 5 along the length and 6 along the width. For the F5, the numbers are 6 and 9.
To proceed with the steps of Sect. 3.3, for the F3 basal planes, it was considered that the length and width of each quasi-plane surface depend on the number and size of spacings of F7 and F1 and vice versa for the basal planes F5. To calculate statistically the total length and width, L A and W A , of the continuous quasi-plane surfaces, random samples of lengths and widths were created. Equations (6)-(9) were adapted from Eqs. (3) and (4), for the two basal planes. Equations (10) and (11) for the area calculation were also adapted from Eq. (5).  Table 1 Probabilities P s (i) and P t (i) of an F3 or an F5 quasiplane surface containing exactly i spacings, and N number of quasi-planes with i spacings, for the 4760 and 3920 observed planes of F3 and F5 s, sliding direction, t, transversal direction - where L F3 is the total length of quasi-plane, along the sliding direction, for the set F3. L F5 is the total length of quasi-plane, along the sliding direction, for the set F5. W F3 is the total width of quasi-plane, along the transversal direction, for the set F3. W F5 is the total width of quasi-plane, along the transversal direction, for the set F5. s F1 and s F3 are the spacing of the intersecting set F1 or F3, respectively. A F3 and A F5 are the total area of the basal quasi-plane for the set F3 or F5, respectively.
It is assumed that each quasi-plane contains a number of fixed spacings s F1 and s F7 , although in reality each quasiplane might contain different spacings. As over the slope there are homogeneous areas, with denser or sparser spacing, the variability of the spacings within each quasi-plane can be considered low and a fixed value can be assumed. For the generation of random samples for the quasi-plane length and width, the number of spacings i follow the PMF of Table 1.
The results of the median and maximum dimensions, and standard deviations, are shown in Table 2. As a validation, the results from the application of this procedure were compared with the dimensions of the planes F3 and F5, on the actual relief, representing historical rockfall scars. Their real dimensions were observed and measured on the LiDAR point cloud. They are summarized in Table 2.
The maximum calculated length, width and area from the probabilistic simulations are for the F3: 28.0 m, 48.5 m, and 282.5 m 2 , which closely approximate the respective values of 27.1 m, 32 m and 236 m 2 , measured with LiDAR by Santana et al. (2012). For the F5, the differences are slightly higher, with the calculated dimensions larger than the observed. The former was found to be 48.5 m, 50.4 m, and 289.3 m 2 in comparison with the latter which was 14.7 m, 19.5 m and 144 m 2 . The calculated and observed median and standard deviation values for the length and width are similar. The exception is the median and standard deviation of the areas, which have a difference of one order of magnitude between the calculated and observed surfaces.
The theoretical maximum values were also calculated, applying Eqs. (6)-(8) for the maximum identified number of spacings i and their size (Table 3). Given the sufficiently large number of samples in this step (100,000) the calculated maximum from the simulations is equal to the calculated theoretical maximum, for the lengths and widths. However, the maximum areas from the simulations are one order of magnitude less than the theoretical. This is due to the random combinations of widths and lengths, which despite the large number of generated samples do not reproduce the theoretical maximum areas. Considering the sample size, the statistical probability of the latter is less than 10 -5 .  Table 2 shows that for the set F3, the calculated dimensions of the generated basal surfaces areas are very similar to those observed with the LiDAR. This confirms the hypothesis argued by Corominas et al. (2018), that the presence of the F1 and F7 restricts the extent of the F3 sliding surfaces and bounds the size of the slope failures. Additionally, as the calculated maximum rock block dimensions fit the size of the maximum observed scars, it is implied that the principal hypothesis of the calculations performed here is certain and that the displacement of the basal planes F3 is the principal cause of rockfall size restriction in the study site. Nevertheless, for the set F5, the calculated length, width as well as areas present larger differences compared to the observed. The size of the calculated areas, despite being of the same order of magnitude, is almost double of the observed. This suggests that slightly larger failures than the historical ones are possible, still in the order of few hundreds of square metres.

Discussion
The differences of order of magnitude in the median and standard deviation of the areas between the calculated and observed surfaces imply that their size distribution is not the same. This is attributed to the polygon shape of the real surfaces, instead of the rectangular form assumed in this process.
The effect of the criterion (b) (requirement for sliding of progressively higher spacings as moving from the exterior towards the interior of the slope) was evaluated, comparing the number of spacings and the maximum expected length of the F3 and F5 quasi-plane surfaces, first applying the criterion and then without applying it. For the set F3, the application of this criterion leads to a difference of just one spacing, while for the F5 the effect is stronger and involves three spacings. In general, as expected for the calculated surfaces, the widths, where the restriction of the criterion (b) is not applied, are longer than the lengths, where the criterion (b) is applied. All the same, this is valid for the observed basal planes, which implies that there is indeed a negative effect of the downward steps for the extent of the basal planes, in particular for the F5 set.
Considering the afore-mentioned restrictions in the extent of the rockfall basal planes, and assuming that each quasi-plane surface corresponds to at least one event, the maximum rockfall volume is restricted. The exceedance of this maximum could occur in the event of substantial breakage of intact rock bridges or in the event of the creation of an extensive stepped path failure surface. In both cases, various quasi-plane surfaces would merge into a larger basal plane. However, as afore-mentioned there is no field evidence supporting the occurrence of large slope failures (larger than 100,000 m 3 ) in the Solà d' Andorra, at least during the last 10.000 years (Corominas et al. 2018). Similarly, the breakage of substantial rock bridges cannot be evidenced by field observations, suggesting that for the given geological structure, the intact rock resistance is sufficiently high to prevent extended bridge failures. Besides, the intensely fractured, nevertheless prominently stepped, rock pattern favours the separate detachment of rock masses from the individual quasi-plane surfaces, rather than their merging into a massive simultaneous rock mobilization. Nonetheless, further studies, beyond the work presented here, are required to assess the credibility of a large rockfall scenario due to progressive bridge failure, incorporating the mechanical properties of the fractured rock mass into the analysis of slope stability.

Conclusions
The proposed probabilistic procedure for the calculation of the size of the basal quasi-planes, constituting sliding surfaces for the rockfall detachment, provides coherent results with the size of the observed scars on the slope in terms of order of magnitude. The advantage of the developed approach is that, using as a starting point the displacement of the basal planes by the lateral intersecting joints, it effectively indicated the extent of penetration of the basal planes into the rock mass, overcoming the difficulties for its assessment.
In the study site, it has been argued that the presence of the F1 and F7 joint sets restricts the extent of the F3 and F5 sliding surfaces, bounding the size of the potential slope failures (Corominas et al. 2018). The results of this work bring additional evidence for this restriction, with quantitative information for the expected penetration and size of the basal planes. They indicate that the maximum joint penetration of the set F3 is 28.0 m and 48.5 m, accordingly in the sliding and transversal direction. For the set F5, the maximum joint penetration is 48.5 m and 50.4 m. The probabilistic simulation indicated basal areas of 282.5 m 2 , for the F3, and 289.3 m 2 , for the F5. This is one order of magnitude lower than the maximum theoretical areas (1539.1 m 2 and 2446.4 m 2 ), which were assessed considering the maximum number and size of spacings. However, considering the sample size, the statistical probability of the latter is lower than 1 × 10 -5 .
The application of the proposed procedure for the calculation of the penetration of the basal planes into the slope, was based on the random generation of spacings of the basal planes with respect to a reference point O (0,0), provided that no detailed data for the displacement of the basal planes by the intersecting joints existed. In that respect, the real displacement has not been integrated into the model. This is a limitation to be improved if possible, after systematic collection of data, through field observations or using remote sensing methods (e.g. UAV-based high resolution image processing). Given that the methodology is sensitive to the selection of the threshold to be applied for the continuity criterion (a), and as a future step, close field observations of the rock surface steps and detailed stability analysis, to provide an insight into the breakage of bridges, are suggested for a more systematic evaluation of it.
As a further development, the proposed methodology could set the basis for assessing maximum rockfall volume, using as a starting point the area of the sliding planes. In that case, the calculation of the probability distribution of the heights, which is related to the location of the respective quasi-place surface along y ik (Fig. 4), would be additionally required.
Acknowledgements This work has been performed within the framework of the research project RockModels, funded by the Spanish Ministry of Economy and Competitiveness, and co-funded by the Agencia Estatal de Investigación (AEI) and The European Regional Development Fund (ERDF or FEDER in Spanish) on the framework of the State Plan of Scientifico-Technical Research and Innovation with reference code BIA2016-75668-P (AEI/FEDER.UE).
Funding This study was funded by the Spanish Ministry of Economy and Competitiveness, and co-funded by the Agencia Estatal de Investigación (AEI) and The European Regional Development Fund (ERDF or FEDER in Spanish) on the framework of the State Plan of Scientifico-Technical Research and Innovation with reference code BIA2016-75668-P (AEI/FEDER.UE).

Compliance with ethical standards
Conflict of interest The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.