Reconstruction of Surficial Rock Blocks by Means of Rock Structure Modelling of 3D TLS Point Clouds: The 2013 Long-Chang Rockfall

Terrestrial laser scanning (TLS) enables the digital representation of rock outcrops with unprecedented resolution and accuracy. Currently, extracting full-scale block structures from raw TLS data (e.g. for assessing rockfall source areas) requires special techniques. This contribution describes processing methods for reconstructing three-dimensional representations of blocky rock masses from raw point cloud data. The developed workflow involves extracting structural and topographical details, fitting location-dependent discontinuities bounding the rock blocks, and reconstructing the block array embedded in three-dimensional geomorphology. The method is applied to a case study of the 2013 Long-Chang rockfall located in the Guizhou Province of China. Based on these initial results, the workflow will be applied to rockfall source evaluations throughout the karstified Kaili region within Guizhou.


Introduction
Slope engineering in blocky rock masses requires a thorough understanding of the rock mass discontinuity structure, and how the discontinuities intersect with topography to form kinematically removable blocks. The set-based discontinuity orientation, location and 3-dimensional extent of block-forming discontinuities are of particular importance. Traditionally, the necessary data for rock slope engineering are measured using contact-analog methods, where a compass is placed on the rock surface to perform scanline and window mapping. However, hazardous conditions in rugged alpine areas can preclude direct contact measurements. In such cases, remote outcrop characterization using terrestrial laser scanning (TLS) offers significant advantages. Recent reviews (e.g. Sturzenegger and Stead 2009;Abellán et al. 2014) have shown successful applications of TLS for characterizing rock slopes and performing rock slope monitoring.
To detect planar discontinuities, some researchers (e.g. Slob et al. 2005;Kemeny and Donovan 2005;Lato and Vöge 2012) interpolate and reconstruct the TLS point cloud data as a meshed model. Through geometric analysis of this 3-dimensional mesh and stereographic plotting of the facet orientations, the orientation of discontinuities sets can be efficiently delineated. Efforts have also been made to detect non-planar discontinuity surfaces directly from point clouds (e.g. Jaboyedoff et al. 2007;Ferrero et al. 2009;Riquelme et al. 2014). These approaches seek to estimate discontinuity parameters such as the set-based means of orientation, spacing, and persistence. These statistical data then serve as inputs for the generation of discrete fracture networks (e.g. Ferrero et al. 2009;Fekete and Diederichs 2013;Sturzenegger et al. 2014).
There are also challenges in rock engineering where the exact spatial structure of a discontinuity system should be established to evaluate the main risks. One of the examples is rockfall source areas evaluation that requires the accurate characterization of block systems, from which the kinematics and stability conditions for individual blocks can be established. This paper demonstrates a processing procedure for reconstructing location-dependent 3-dimensional representations of blocky rock masses using the TLS datasets of the 2013 Long-Chang rockfall in karstified limestone. Structural and topographical details were extracted directly from the raw TLS points, and in conjunction with Block Theory (Goodman and Shi 1985), removable blocks were identified, analyzed in terms of 3-dimensional stability and finally visualized in the 3-dimensional digital terrain. Based on these initial results, the elaborated processing methods will be applied to rockfall source evaluations throughout the geographic area of Kaili, located in the Guizhou Province of China.

Long-Chang Rockfall
In February 2013, two consecutive large-scale rockfalls occurred from a calcareous cliff above the Yudong River of the local community Long-Chang in the Kaili area of China. The rockfall deposit has a cumulative volume of approximately 300,000 cubic meters and dammed the Yudong River. Five people remain unaccounted for as a result of the disaster, and 79 residents residing in 21 houses were evacuated.
Kaili and its surroundings are located in the bare karst area of southern China (Yuan et al. 1995). The ubiquitous talus deposits blanketing lower portions of the steep limestone hills record a history of rockfalls. Although studies regarding karst hydrogeology, karst geomorphology and abrupt karst collapse in southern China have been made, the influence of karstification in triggering rockfalls and other forms of slope instability has not been evaluated in detail.
Field investigations performed during the summer of 2013 included high-resolution TLS of the 2013 rockfall cliff of Long-Chang, and field inventory of rockfall prone cliffs in the greater Kaili area. This resulted in the identification of eight additional rockfall study sites, which were subsequently surveyed with TLS. During the field investigations, representative rock samples were collected and subjected to laboratory testing to determine the physical and mechanical properties of the Kaili limestone.
Karstified terrains are a distinctive landform of southern China, with the South China Karst mainly covering the provinces of Yunnan, Guizhou and Guangxi (Yuan et al. 1995). The area of continuously distributed carbonate rocks amounts to 500,000 km 2 with a mean thickness of 1500 m (Zhao et al. 2012). Kaili is located 200 km west of Gouyang in the province of Guizhou and has a subtropical humid climate with mean annual rainfall of about 1300 mm. The mean annual temperature ranges between 15 and 21 °C.
Kaili is located along the eastern margin of the Yunnan-Guizhou Plateau. The area was uplifted during the Yanshan Movement of Middle Jurassic time (Li et al. 2014). Fold axes and major faults in the Kaili area trend to the NNE, with faults dips to the WNW and ESE at a mean dip angle of more than 65°. Since the Cenozoic, intermittent episodes of uplift resulted in deep incisions into the karst terrain by surface and subterranean water flow. The resulting landscape exhibits the form of karst hills and depressions, with a mean relief of 300 m (Lu et al. 2013). The Kaili carbonates were deposited between Carboniferous and Permian times. The karstified calcareous strata are Baizuo Formation (C 1b ) of the Carboniferous as well as Qixia (P 2q ) and Maokou (P 2m ) Formations of the Permian (Dong et al. 2015).
The 2013 Long-Chang rockfall (Fig. 1a) involved two consecutive events. The failure on the left side (Area I) occurred first, and 1 day later the failure on the right side (Area II) occurred. The entire failure surface (Areas I and II in Fig. 1b) has a height and width of approximately 190 and 130 m. Bedrocks exposed in the headscarp area are flatly bedded limestones of Qixia and Maokou Formations. The limestones contain several thin-bedded interlayers of calcareous shale (Dong et al. 2015). The uppermost part of the headscarp is overhanging, and this topographic detail appears to be associated with the collapse of the karst cavity shown in Fig. 1c, d. The karst cavity extends from the ground surface to a shale layer at a depth of approximately 30 m and it has an average diameter of about 15 m. According to Sauro (2003) this type of karst cavities represents an inception doline that developed from the concentration of surface runoff and subsurface flow concentrated along the vertical fractures that are bounded by an impermeable shale layer at the bottom of the doline. Field inspections of the headscarp area revealed sixteen additional doline structures.

Methods
Digital reconstruction of the Long-Chang Rockfall involved data collection, followed by the extraction of 3-dimensional rock structure and computation of block molds. On the basis of this reconstruction, combined with the results of uniaxial compressive strength (UCS) and direct shear tests, backcalculation of stability conditions and the failure sequence was performed. The details of digital reconstruction are enumerated below, together with results of the UCS and direct shear testing. Results of the reconstruction and details of the stability analyses are summarized in Sect. 4.

Data Collection
The Long-Chang Rockfall site is well-suited for making remote engineering surveys, as the rugged topography precludes safe personnel access. The study area was surveyed with an Optech ILRIS terrestrial laser scanner. Overlapping scans obtained from all positions were merged using shapefitting algorithms of the software PolyWorks (InnovMetric 2018). After the merging, each laser point was referenced in the scan project Cartesian coordinate system (x, y, z) with X = E (East), Y = N (North) and positive Z = up.

Extraction of 3-Dimensional Rock Structure
Because of limited sensor range, high light absorption and occlusions during scanning, any 3-dimensional rock surface measured with TLS is non-uniformly sampled and may contain undersampled domains Kieffer 2011, 2012). The variable sampling density affects the correctness of meshing (Fabio 2004;Berger et al. 2016) of rock surfaces. As a result, to spatially identify how discontinuities organize to form distinct removable blocks on a natural rock slope, raw point clouds are utilized rather than a derivative meshed surface.
In a point cloud, rock discontinuities and their spatial pattern can be discerned, but the data often includes obstructions such as vegetation. Because automated (algorithmbased) cloud cleaning is still a challenge (Sotoodeh 2006), and to ensure high confidence in the digital mapping results, the points falling on and defining a discontinuity surface are selected manually. The x, y and z coordinates of the selected points are then used to fit 3-dimensional discontinuities (e.g. Liu and Kieffer 2012). Another advantage of manual localization and discrimination of a geological discontinuity is that the interpreter can incorporate engineering geologic perspective and judgements that cannot be replicated by existing algorithms. The main drawbacks of manual procedures relate to efficiency and subjectivity (reproducibility). For this reason, semi-automated algorithms are considered ideal.
The key to modeling the rock structure directly with a point cloud is to estimate the direction of each point in space, which is the normal vector perpendicular to the surface on which the point is located. Consequently, the points on the same discontinuity can be extracted in space due to their identical or similar direction values. To compute normals on a point cloud, it is necessary to consider the relations of a point to its neighbors. Jaboyedoff et al. (2007) describe an approach using the eigenvalue of the covariance matrix that is defined by the nearest neighbors of the queried point. Riquelme et al. (2014) presented a semi-automatic calculation of orientations based on a neighboring point coplanarity test. Our goal is to delineate the block-forming discontinuities and identify how various discontinuities are spatially organized to form blocks. An ideal algorithm can identify planar surface regions and correctly reconstruct edge intersections. Several algorithms have been tested, including the best fitting plane approach (e.g. Hoppe et al. 1992), local 2D triangulating (e.g. Dey and Goswami 2004), quadric surface approximation (e.g. Groshong et al. 1989;Yang and Lee 1999) and the Hough transformation Marlet 2012, 2016). As an illustrative example, consider three discontinuities terminating as shown in Fig. 2a. The discontinuity pattern results in three sharp intersecting edges and a single corner where the three planes intersect. The discontinuity dips/dip directions are 90°/140° (F1), 70°/50° (F2), and 20°/230° (F3). The three square planes have the same side length of 1 m. Figure 2b shows a point cloud model of the three correlated planes, which were discretized using the same regular grid of 0.025 m × 0.025 m in each plane. The normal of each point is then estimated using four different algorithms and the results are depicted in Fig. 3. Each algorithm properly reconstructs the point normals within the plane area, but only the Hough transform implemented in Boulch and Marlet's method can correctly reconstruct the normals on the intersecting edges as shown in Fig. 3d, where the normal of an edge point is consistent with one of the normals of the two intersecting planes. In Boulch and Marlet's method, the starting point for computing is randomly selected, so the normal of a point on a sharp intersecting edge can be any one of the two intersecting planes and the normal assignment to any one of the two intersecting planes is also correct. Other estimation techniques tend to smooth the sharp features ( Fig. 3a-c). As a result, the Boulch's algorithm (Boulch 2018) is considered most suitable for computing discontinuity normals. The next step is to convert the calculated mathematical normals to a perceptible 3-dimensional image of rock structures. Point cloud processing software aids an interpreter in efficient filtering and rendering the 3-dimensional rock structure. The software COLTOP-3D (Jaboyedoff and Couture 2003;Jaboyedoff et al. 2007) is perhaps the earliest rendering technique to analyze rock-slope relief using a digital elevation model (DEM) and 3-dimensional point clouds (Metzger et al. 2009). The approach is to assign a unique color to a certain spatial topographic cell quantified using the normal direction to the cell. In this way, the points on a continuous planar structure will be rendered as the same color. This automated pre-processing step greatly facilitates discrimination of discontinuities and objective interpretations of the 3-dimensional rock structure. In our approach, each point which has Hough's normal orientation is colored using the HSV-rendering technique (Liu and Kaufmann 2015). As shown in Fig. 4a, the HSV-rendering is represented by a cylindrical color space that can be modelled as a cone. The color hue (H) and saturation (S) are defined by the bearing direction and plunge magnitude of the normal, respectively. A pure color is fully saturated (S = 1) when the normal lies at the edge of the hemisphere projection circle (i.e. vertically dipping discontinuities). Conversely, white color (S = 0) represents horizontal discontinuities that have their respective normals plotting at the center of the projection hemisphere. The value (V) of a color describes how dark the color is. In our 3-dimensional structure models, a consistent lightness of V = 0.75 has been adopted. As an example, Fig. 4a depicts six typical pure colors with full saturation (S = 1). In Fig. 4b each normal of 1° resolution, as poles in equal-angle and low-hemisphere projection, is allocated to a unique HSV-color. We use this HSV-color wheel herein for rendering 3-dimensional rock structures. Because a certain orientation is rendered only with a fixed color, the points on the same plane will appear in the same color (Fig. 4c). The points with the same color are then fitted to a discontinuity plane (Fig. 4d). The values of rendering and fitting are given in Table 1. The points on the same plane were rendered with the same color (column 2 of Table 1) and then fitted to the discontinuity plane using least-squares fitting (column 3 of Table 1). Column 4 and 5 are the positions and area size of the fitted planes.
The final step of our 3-dimensional structure extraction is to assign each fitted discontinuity a half-space code (HSC) according to Block Theory (Goodman and Shi 1985). As an example, the last column of Table 1 shows the HSC of each discontinuity in Fig. 4, in which the digit 1 and 0 represent the half-space below and above a plane, respectively. In this example, the regions of intersections of three half-spaces are defined by the joint pyramid (JP) 110, which corresponds to the block formed on the lower half-space of F1, the lower half-space of F2 and the upper half-space of F3 (Fig. 4d). Since F1 has a dip of 90°, its left side is defined here as below space. The JP identification based on HSC is decisive so that the in situ block system can be characterized using the fitted planes.

Computation of 3-Dimensional Block Molds
A block mold represents the negative space that remains after a block detaches from the rock outcrop. Detailed mapping and analysis of block molds provide important information regarding actual block failures that have occurred and conditions that favor the development of unstable keyblocks.
Our computations concerning block molds include the location, shape, volume of the mold, followed by assessments of the corresponding block's kinematic removability, failure mode, and limiting equilibrium conditions.
To find a removable block intersecting a slope surfaces the orientations of the slope faces must be specified to define the excavation pyramid (EP) according to Block Theory. As an example, the synthetic block mold depicted in Fig. 4d represents a block having three free planes E1, E2, and E3 (Table 2). We suggest that the surficial geometry (i.e. EP) of this synthetic rock slope is shaped by discontinuities having the same orientation of planes forming the block mold. The free (or excavation) planes are almost parallel to and 1 m from the respective discontinuities (Table 2 and Fig. 4d). The HSC of this EP has the digits 001. The closed rock  (Hoppe et al. 1992). b 2D triangulating (Dey and Goswami 2004). c Quadric (Yang and Lee 1999). d Hough transform (Boulch and Marlet 2016). To show the normals of the intersection edges clearly, the intersection corner has been rotated to the reader in all of the four figures, i.e. the convex side of the intersection is pointing to reader block referred as the block pyramid (BP), is the intersection (∩) of the JP and EP (i.e. BP = JP∩EP = 110∩001). The complementary space region of EP is referred as space pyramid (SP), i.e. 110. According to Shi's Theorem, removable blocks satisfy the criteria that the JP is nonempty, and the JP is completely contained within the SP (Goodman and Shi 1985). These two requirements can be readily checked via stereographic projection: (1) every JP is nonempty if it appears in the stereographic projection; and (2) the JP is removable if it is located entirely within the SP. Figure 5a shows the projection of the six boundaries of the BP defined in Table 2, using whole sphere upper-focalpoint stereographic projection. The dotted line represents the reference circle, and solid circles represent the three discontinuity planes (F1, F2, and F3). Dashed circles are the free or excavation planes (E1, E2, and E3). The lower and upper half spaces of a plane lie inside (expressed as digit 1) and outside (expressed as digit 0) its corresponding circle, respectively. For instance, JP 110 is the shaded region on the left side of Fig. 5a, which is simultaneously inside the circle of F1, inside the circle of F2 and outside the circle of F3. Because JP 110 is completely within the SP, it is kinematically removable. The digit "3" in parentheses below the JP digits 110 of Fig. 5a indicates the failure mode, in this    Fig. 4d The dips of the three EP planes (E1, E2, and E3) are 1° different from their corresponding discontinuity planes so that the great circles of their stereographic projection can be clearly differentiated from those of the discontinuity planes in Fig. 5a Plane case planar sliding along F3. Other possible modes in Block Theory include "0" for a block lifting or falling without having contact to any bounding plane, a single digit for sliding along the respective joint plane, e.g. 3 for sliding along F3, and double digits for wedge sliding along the intersection of the respective two joint planes, e.g. 13 along the intersection of F1 and F3. The stability of JP 110 depends on the effective friction angle on the actual sliding plane F3. The friction angle required for limiting equilibrium of JP 110 can be readily assessed in 3-dimension via stereographic projection (Goodman and Shi 1985;Goodman 1989). In Fig. 5b, the location of the gravitational force R reveals that the required friction angle for JP 110 is 20°. Additional block computations are performed to determine block volumes and surface areas.

Laboratory Tests and Results
The Kaili limestone bedrock exhibits varying grades of weathering. Selected representative rock samples were subjected to uniaxial compressive testing and direct shear testing. The goal was to assign representative strength parameters according to weathering grade, and scanning electron microscopy was utilized to facilitate correlations. Figure 6a shows boxplots of the unconfined compressive strength (σ c ). The slightly weathered limestone has a mean σ c of 74 Mpa and for moderately weathered limestones the mean σ c is 33 Mpa. Here the class "slightly weathered" means that the limestones are just discolored along micro fractures, some of which may have been opened slightly (Fig. 6b). The slightly weathered samples were taken in the interior of the failed blocks. Moderate weathering, on the other hand, is the collectively weathered or altered part in the vicinity of highly weathered joints surfaces, on which the rocks have brownish-red colors (Fig. 1c, d) and a significant portion of calcites has been removed by dissolution (Fig. 6c). The brownish-red colored surfaces (Fig. 1b, d) evidence the water circulation passing through open joint system. Figure 7 shows the peak shear strength parameters of the intact limestones having different weathering, which help us to appropriately find post-peak residual shear strength on discontinuities.  Table 2 to examine the removability of JP110. Sliding modes are indicated with parentheses. b 3D limit equilibrium analyses for JP 110

The Roughness and Persistence of Discontinuities at the Natural Scale
The area size of the main failure planes of the Long-Chang rockfalls (Fig. 1b) varies between 5000 and 9000 m 2 . The mean dips of the failure planes are all steeper than 80°. The total volume of the first collapse body is approximately 200,000 m 3 and the second about 40,000 m 3 (Dong et al. 2015). The dry unit weight of Kaili limestones changes from 2.50 to 2.73 g/cm 3 . Due to the identifiable solution features (Fig. 1b-d) we take the low limit of 2.5 g/cm 3 to represent the weathered blocks. Thus, the mean normal stress (σ n ) of the failure planes would be less than 0.15 MPa. Barton's ratio (Barton 1973(Barton , 1976) of unconfined compressive strength σ c to the normal stress σ n , i.e. σ c /σ n , would be much higher than 100, which indicates that the shear strength of discontinuities under low normal stress is frictional and consisted of dilation angle (d n ) due to joint roughness and the residual frictional angle (ϕ r ) of weathered limestone. Compared to the peak shear strength of 35.05° (Fig. 7), a residual friction angle of 30° was assumed to represent karstified discontinuities of Kaili limestone. As enumerated in the results section of this paper, block kinematics, failure modes and stability are controlled by very steeply dipping discontinuities and are thus rather insensitive to the selection of residual shear strength.
Using the high-resolution point cloud data, the roughnessrelated dilation angle of a sliding plane is determined. All points of the undulating surface are first selected and the orientations of their Hough's normals are plotted in equal area stereographic projection. The peak dilation angle is estimated as the angle between the cone axis (mean join plane) and the surface formed by the maximum scatter of the surface normals. As an example the grey discontinuity surface of Fig. 8a is discretized by the point cloud, on which the Hough's normal of each point is represented as a short black line. This surface has a mean attitude of 3°/289°, shown as a great circle (global mean plane) in Fig. 8b and the maximum scatter around the mean pole (87°/109°) indicates a peak dilation angle of 14°.  Fig. 8 a The undulation of the grey discontinuity surface is discretized by 289 points whose normals can quantify surface roughness. b The peak dilation angle is the deviation of the local normal orientations from the mean orientation Field observations (e.g. Goodman and Kieffer 2000) of rock slope failures indicate that structurally controlled failures can occur when existing discontinuities have limited persistence, but this requires progressive rupturing of intact rock bridges. Identifying and documenting the location of rock bridges is thus important for understanding the conditions promoting the initiation and progression of failure.

HSV-Colored Cliff for 3-Dimensional Rock Structure
The TLS cloud of the Long Chang cliff includes 17,7333,810 points (Fig. 9a). The portion of the point cloud corresponding to the 2013 rockfall is marked with a white frame in The intact limestone is generally brighter than the scree deposits and vegetation. To elucidate the 3-dimensional rock structure the normal of each point was computed by Hough's method, and the two normal orientations of a point were then used for visualization, i.e. the outerward pointing normal is rendered as white and the inward pointing normal as black. As shown in Fig. 9b, such rendering greatly facilitates the differentiation of discontinuities and their spatial organization. Location-dependent rock structure is quantified by assigning a specific HSV color based on the normal orientation of the respective point. Figure 9c shows the HSV colored rock structure of the cliff. The HSV rendering shows that the surface relief of the cliff is dominated by systematic structural elements, including joint surfaces forming the overhanging (green) and non-overhanging (violet) areas of the main slope face, pervasive bedding planes which are the flat-lying boundaries between green and violet areas, and the two additional sets of persistent joints (blue and yellow) intersecting the cliff surface. In the following, we use the part of 2013 rockfall (Fig. 10a) to show the HSV-color based extraction of in situ rock structure.

Filtering Joint Sets by HSV-Colors
The distribution of HSV colors in Fig. 10a indicates the existence of seven discontinuity sets: the pink and violet areas are clustered into one set; the green is the second set; the blue is the third; yellow the fourth; aqua the fifth; and red the sixth. These six sets are all joints. The seventh set, i.e. the very flat-lying bedding surfaces, is traceable by means of the undulation of the main slope faces (green and violet). Figure 10b-h shows each of the individual discontinuity sets that were extracted from Fig. 10a.

Set orientation and Kinematic Analysis
To reduce sampling bias due to slope orientation, approximately 35 small areas are randomly selected in each of the filtered HSV-colored sets (Fig. 10b-h) and the points of each selected area are fitted to a plane to obtain the average set orientation over the domain sampled. The orientations of Fig. 11 Poles and contours of the discontinuity sets plotted using the code Dips 7.0 (Rocscience 2018a). The set-based attitude is displayed on the right side Table 3 Results of kinematic analyses in code Dips7.0 *Critical percentage is the ratio between the number of poles within the potential failure zone and the number of the total poles involving in the analysis; OT means oblique toppling, DT direct toppling  248 sampled areas are plotted as poles in Fig. 11. The cliff has a mean slope orientation of 84°/106°. Table 3 shows the results of kinematic analysis. With the assumed residual friction angle of 30°, discontinuity sets J 11 , J 21 and J 32 (dipping out of the cliff) tend to fail in-plane and wedge sliding modes, while the sets J 12 , J 22 and J 31 (dipping into the cliff) are prone to toppling failure. To assess how failure mode tendencies change according to the frictional resistance of discontinuities, the effective friction angle was systematically increased from the residual value up to 65°. As shown in Table 3, the failure modes remain constant over this range of friction angle.

Reconstruction of Long-Chang Rock Fall
The failure scenario analysis included: (1) global stability analysis of the entire failure area assuming a single large block was involved; and (2) separate analysis of block molds to reconstruct the first and the secondary collapse events.  is delimited by set J21, and the back boundary (BB, stippled transparently in Fig. 12b), is the mean of sets J11 and J12. The entire front plane (FP) is modelled by simply connecting the front edges of TP, LB, BB, and BT as shown in Fig. 12c. Table 4 shows the input parameters for computing the block. Block Theory analyses show that to be removable, the whole block must have two free boundaries, i.e. bottom (BT) and front boundary (FT). The JP 100 is the common space of below top (TP), above left boundary (LB) and above back boundary (BB). The computed block is shown in Fig. 12d. The failure mode involves wedge sliding along the intersection line of LB and BB. The required friction is 80°, indicating the global block could have been stable only by virtue of rock bridges.

Global Analysis of the Entire Failure Area
Using the technique described in Fig. 8, we concluded that the failure plane LB has a dilation angle (d n ) of 18° and the failure plane BB has d n = 20° (Fig. 13). With the frictional shear strength (d n + ϕ r ) estimated at 48° for LB of set J 21 and 50° for BB of set J 11 , which would give a factor of safety (FOS) of 0.2291 based on the stability analysis of the code SWEDGE (Rocscience 2018b). This indicates that an additional component of shear strength, formed by rock bridges as shown in Fig. 14, had been necessary for stability. Back-analyses in SWEDGE indicate that a cohesion of about 0.12 MPa must be combined with the friction angles mentioned above, to maintain equilibrium. This low cohesion of rock bridges needed for equilibrium is in good agreement with the results relating to the general amount of rock bridge cohesions of other researchers (e.g. Kemeny 2003Kemeny , 2005Eberhardt et al. 2004;Paronuzzi and Serafini 2005;Frayssines and Hantz 2009;Paronuzzi et al. 2016).
Through systematic block theory analyses, a plausible mechanism of block failure involves wedge sliding. However, considering the slender tabular form of the block, it is also imperative to check its susceptibility to toppling. This requires knowledge of block geometry and the friction angle (ϕ) between the block and the surface on which it is resting (Hudson and Harrison 1997). The block geometry is defined by the ratio b/h of base length (b) and height (h) of the investigated block (Fig. 15a). The limiting equilibrium condition of direct toppling is that the location of the line of action of the force due to gravity, which passes through the center of gravity of the block, will coincide with the Fig. 16 a Fitting a waving discontinuity to homogeneous segments. For clarity the HSV-colored point models of J 11 and J 12 have been switched off. b A mold circumvented by discontinuities in detail. c Zooming in on boundaries and form of the mold. The front face is displayed in raster ▸ 1 3 lower apex of the block, i.e. b/h = tan(ψ) where ψ is the dip of the surface on which the block is resting. Toppling will not occur if b/h > tan(ψ), and will occur if b/h < tan(ψ). Due to the irregular shape of the actual block and with respect to the sliding surfaces LB and BB, the actual block has a mean base length of b = 155 m (Fig. 15b) and a mean height of h = 31 m (Fig. 15c). The ψ angle is now the dip of the intersection between LB and BB, i.e. ψ = 80°. Because of ψ (= 80°) > ϕ (=50°) and b/h (= 5) < tan(ψ) (= 7.12), the whole block would slide and topple simultaneously.

Failure Sequence
To reconstruct the first and the second collapse events, wavy discontinuities are divided into several planar segments, e.g. the left blue big joint in Fig. 16a. In Fig. 12 this joint is fitted as a single plane LB; however, detailed segmentation of this structure reveals details of block-forming discontinuities. As an example, Fig. 16b shows a block mold having the HSC 1000 that involved plane sliding mode along the lower blue joint in Fig. 16c. This failure plane is depicted in Fig. 17a as plane 4, and it is a member of J 21 . The required friction angle for this block is 65°. Since the actual friction of set J 21 is about 48°, the JP 100 would have also remained stable only by virtue of rock bridges. Figure 17b shows the determination of the block b and h in relation to the sliding plane to analyze the kinematic of direct toppling. The simultaneous fulfilment of criteria ψ (= 65°) > ϕ (=48°) and b/h (= 1.312) < tan(ψ) (=2.143) confirms that this block would also rotate about the lowest contact edge of its sliding plane during the failure.
The extent of the two blue and almost vertical joints in the middle of Fig. 16a, here referred as middle boundary (MB), represents the border of the first and the second collapse events. This boundary is also observable both in the point cloud models of Fig. 12a and in right photograph of Fig. 14. Using the fitting technique described above, the rest of the side discontinuities delimiting the block molds are precisely tracked and digitally mapped in 3-dimension (Fig. 18a). Totally 38 mold planes have been identified. The failed blocks are then reconstructed separately. A total of 15 blocks, in the first collapse event 7 and in the second 8, have been computed and visualized in the HSV colored cloud (Fig. 18b).
The results of block computing, kinematic and stability analyses are summarized in Table 5. In the first collapse, except the topmost block B_16 whose failure mode is free fall, other six blocks have the kinematics of simultaneous sliding and toppling. Table 5 shows that the slope face (SF) was the only free space that block B_9 (Fig. 18b) required to move into, i.e. the B_9 failed first. Details on the computation of B_9 are given in Fig. 17. The initialization of the whole first collapse requires that the block B_9 failed at first, which would on the one hand make the space into which B_8 can move and on the other hand switch on the toppling of the overhead blocks. The movement of B_8 again would function as the initiator of movements of its overhead blocks, etc. Here, we do not examine the exact locationdependent relationship between sliding and toppling using 3-dimensional stress-strain analysis, but the modelling of the in situ rock structure based on remote sensing data has been justified to be very practical for the assessment of rockfall sources in the high and steep limestone slopes of Kaili area. Another special feature of 2013 Long-Chang failures are the effect of karstification on the mechanical behavior of the middle boundary (MB). It can be seen from Table 5 Fig. 17 a Block computing and visualizing, background planes in white. b Estimating b and h to examine conditions for sliding and/or toppling that if B_6 and B_4 are removable, the blue MB in Fig. 18b must be free and open. The locations of B_6 and B_4 are exactly where the doline is. Indeed, the set J 21 is the most significant structure that affects the failure pattern of the whole Long-Chang calcareous cliff, as shown in the HSVstructural model of Fig. 9c. The aerial photograph after the failure also shows the spatial extension of this steep and highly persistent joint set J 21 (Fig. 1a), along which inception dolines have developed. This confirms the general link between slope instability and the karstic processes observed by other researchers (e.g. Santo et al. 2007;Gutierrez et al. 2014). The second collapse also starts at the lowermost block B_17 and progressively expands upwards. The blocks of the middle area (B_14, B_13, and B_12) still undergo sliding and toppling simultaneously, but the blocks at the top just fall free down. In particular, the block B_10, which has a volume of 9924 m 3 , falls down due to lack of support from the left side and from the base. The only resistance during a failure is the cohesion of the small and weathered rock bridges between B_10 and the surface of the joint set J 11 and J 12 . A video recording the right collapse confirms our reconstruction of the failure sequence.

Conclusions
Remotely measured point clouds, such as obtained with high-resolution TLS, are of great value where direct contact measurement is precluded by terrain or safety constraints. The data sets are well suited for making quantitative measurements of rock structures and for performing efficient rock slope stability analysis. The HSV-color representation of points permits a rapid 3-dimensional structural analysis, which is not achievable in most Geographical Information Systems due to the true 3-dimensional nature of point clouds. The results presented herein indicate that the methodology developed for structure extracting and block computing provides an auditable and efficient means for assessing rock fall hazards. The workflow involves: (1) computing Hough's normal for sets of laser points; (2) HSV-coloring of the normal orientations; (3) extracting set-based discontinuity points in space based on HSV filtering; (4) fitting 3-dimensional discontinuity planes; and (5) identification, computing and visualization of rock blocks. The practicality of this developed modelling approach is verified by reconstructing the 2013 Long-Chang rockfall event that developed in karstified limestone. Based on these initial promising results, the methodology will also be applied to rockfall source evaluations throughout the geographic area of Kaili, located in the Guizhou Province of China.