Simulation of fragmental rockfalls detected using terrestrial laser scans from rock slopes in south-central British Columbia, Canada

Rockfall presents an ongoing challenge to the safe operation of transportation infrastructure, creating hazardous conditions which can result in damage to roads and railways, as well as loss of life. Rockfall risk assessment frameworks often involve the determination of rockfall runout in an attempt to understand the likelihood that rockfall debris will reach an element at risk. Rockfall modelling programs which simulate the trajectory of rockfall material are one method commonly used to assess potential runout. This study aims to demonstrate the effectiveness of a rockfall simulation prototype which uses the Unity 3D game engine. The technique is capable of simulating rockfall events comprised of many mobile fragments, a limitation of many industry standard rockfall modelling programs. Five fragmental rockfalls were simulated using the technique, with slope and rockfall geometries constructed from high-resolution terrestrial laser scans. Simulated change detection was produced for each of the events and compared to the actual change detection results for each rockfall as a basis for testing model performance. In each case the simulated change detection results aligned well with the actual observed change in terms of location and magnitude. An example of how the technique could be used to support the design of rockfall catchment ditches is shown. Suggestions are made for future development of the simulation technique with a focus on better informing simulated rockfall fragment size and the timing of fragmentation.

1 Introduction

Rockfall is a mass-movement hazard often found in mountainous environments, posing risk to human lives and the safe operation of infrastructure. In Canada, rockfall hazard is particularly problematic for transportation infrastructure, where traffic corridors have been constructed in steep river valleys, adjacent to natural and cut slopes. In many of these valleys, right-of-way space is limited, forcing operators to construct highways and railways with minimal clearance to potentially unstable rock masses. The proximity and extensive presence of fractured and weathered rock slopes, combined with the seemingly unpredictable nature of rockfall events, makes it difficult for operators to manage this hazard.

Like any geohazard, we are interested in knowing how frequently we can expect events to occur, where they are most likely to come from, and whether or not they will reach and cause harm to our elements at risk. For rockfall this could involve identifying lithologies which are more susceptible to rockfall (e.g. Rosser et al., 2005; van Veen, 2016) or correlating rockfall frequency to triggering events like severe rain storms (e.g. D'Amato et al., 2016; Pratt et al., 2018). In order to build these relationships, we first need knowledge of previous rockfall in the area. This is one of the main challenges in the study of rockfall. Due to challenging terrain in mountainous areas, and safety concerns where rockfall activity is high, site access to the slopes of interest is often limited, prohibiting direct measurement.

In order to circumvent these obstacles, modern remote sensing techniques, such as lidar, have seen widespread use in the study of landslide hazards like rockfall (Jaboyedoff et al., 2012). Terrestrial laser scanning (TLS) in particular has been increasingly applied in recent years to the study of rock slope instabilities (Abellán et al., 2014; e.g. Lato et al., 2009, 2015; Kromer et al., 2017). TLS provides a detailed 3-D model of the slope geometry, supporting a range of geotechnical and geomechanical analyses. From a single scan we can extract slope angle measurements, map rock mass structure and look for evidence of past rockfall (Telling et al., 2017). Using multiple successive scans of the same slope, progressive changes to the slope surface can be measured. Multi-temporal scanning enables researchers to detect rockfall events over time and build detailed magnitude–frequency relationships for slopes, supporting the study of processes such as coastal cliff erosion (e.g. Rosser et al., 2007; Williams et al., 2018), as well as hazard and risk assessments for linear infrastructure such as railways (e.g. van Veen et al., 2017). The identification of rockfall events using sequential scans also permits the back analysis of rockfall-triggering factors and failure mechanisms. A study by Kromer et al. (2015) used TLS change detection to investigate the occurrence of a 2600 m³ failure above a section of the Canadian National (CN) Railway in western Canada, including an analysis of structural constraints, pre-failure deformation, and precursor rockfall leading up to failure.

While the application of TLS to rock slope monitoring is often focused on determining how likely future rockfall is or where the rockfall may come from, it is also important to determine how likely it is that the rockfall material will reach an element at risk should a fall occur. In order to answer that question, rockfall modelling programs can be used to simulate the runout of falling rock material using numerical models (Turner and Duffy, 2012). A rockfall simulation requires a representation of the slope surface and rockfall mass being modelled in 2-D, 2.5-D, or 3-D. TLS can support rockfall simulation by providing a high-resolution 3-D model of the slope surface and, in cases where a rockfall event has been identified using change detection between multiple scans, the volume and location of the rockfall mass. In order to make use of the quality and quantity of 3-D point cloud data being collected as part of a rock slope monitoring program in western Canada, a 3-D rockfall simulation technique using game-engine technology has been developed (Ondercin, 2016; Sala, 2018). The technique uses the video game development platform Unity 3D (Unity Technologies, 2018) and is capable of simulating rockfall runout using fully 3-D meshes built from TLS point cloud data.

One of the strengths of this simulation technique is its ability to simulate rockfall runout using numerous interacting bodies. This capability allows us to produce simulations of fragmental rockfall events, which are defined by the presence of multiple mobile fragments of rock during runout (Hungr and Evans, 1988). Conventional rockfall modelling programs (e.g. RocFall, Rocscience Technologies, 2016; Rockyfor3D, Dorren, 2015; RAMMS:ROCKFALL, Bartelt et al., 2016) simulate single boulder trajectories at a time and therefore are unable to model fragmental rockfall runout. Efforts to model fragmental rockfall processes, including the disaggregation of falling rock masses along pre-existing discontinuities, or the breakage of intact blocks during impact, have been demonstrated by previous authors. Cuervo et al. (2015) utilized a discrete element technique to model the disaggregation of a 1000 m³ rockfall event in southern France comprised of numerous mobile fragments. Wang and Tonon (2011) developed a discrete element code which models the fragmentation of a falling rock at impact, taking into consideration the effect of impact velocity, ground condition, energy loss, and fracture properties such as persistence. The GIS-based tool RockGIS, presented in Matas et al. (2017), incorporates both breakage along pre-existing discontinuities and the generation of new fragments during impact into a 3-D runout modelling program and has been utilized to model a 10 000 m³ fragmental rockfall in the eastern Pyrenees mountains. The goal of the work presented in this paper is to demonstrate the capability of our novel game-engine-hosted simulation technique to model a series of rockfall events detected using TLS change detection at two rock slopes in south-central British Columbia. Emphasis is placed on the ability of the technique to be used for fragmental rockfall runout simulation, including a discussion of how these types of simulation could support mitigation design.

2 Study sites

The rock slopes of focus for this paper are part of the Ashcroft subdivision of the CN Railway. The subdivision follows sections of the Thompson and Fraser River valleys, located between the towns of Ashcroft and Lytton, British Columbia. This region serves as an important transportation corridor for Canada, with the presence of the CN Rail mainline, as well as sections of the Canadian Pacific Railway, and the Trans-Canada Highway. A previous study by Piteau (1977) identified that rock cuts in this region are subject to slope instability issues due to a combination of lateral erosion from river activity, over-steepening from blasting during the construction of the railways, and a lack of adequate rockfall catchment areas.

Two sites in the Ashcroft subdivision will form the basis for this research, White Canyon and Goldpan. The location of these sites along the Thompson River and proximity to the town of Lytton, BC, can be seen in Fig. 1. The monitoring of these slopes is part of the Railway Ground Hazard Research Program, a collaborative research initiative focused on the characterization and assessment of mass movement hazards affecting Canadian railways.

Data collection in the Ashcroft subdivision first took place in 2012, and regular monitoring at both sites has been ongoing since 2014, with data collection taking place on a seasonal basis, approximately every 3 months. Point cloud data collection consists of TLS scans acquired using an Optech ILRIS 3D-ER (2012–2017) or Riegl VZ-400i (2018) lidar system. Additionally, aerial laser scan (ALS) data coverage for both sites, with an average point spacing of 0.3 m, was acquired in 2014 and 2015. Site photos are collected using a Nikon D700 or similar camera with a 135 mm prime lens. A series of photos of the slope are taken using a GigaPan EPIC Pro robotic camera mount and stitched together into a single high-resolution site panorama using GigaPan Stitch (GigaPan Systems, 2013).

Figure 1 Aerial imagery showing the location of the White Canyon and Goldpan field sites situated along the Thompson River, near the town of Lytton in south-central British Columbia. This region is approximately 160 km northeast of the city of Vancouver. (Image source: © Google Earth; Digital Globe 2018.)

Goldpan is located on the north side of the Thompson River, approximately 26 km upstream from the town of Lytton. Present at the base of the slope is a section of the CN Rail main line. The rock slope spans 800 m, rising up to 65 m vertically above track level. The average slope angle at the site is 55–60^∘. TLS data are collected from scanning vantage points across the river from the slope, accessed via Goldpan Provincial Park. Scan distances range between 170 and 230 m, producing an average point spacing of approximately 6 cm.

White Canyon is located on the north side of the Thompson River, approximately 5 km upstream from the town of Lytton. Present at the base of the slope is a section of the CN Rail main line. The rock slope spans 2.4 km, rising up to 375 m vertically above track level. The average slope angle at the site is 40–45^∘. TLS data are collected from scanning vantage points across the river from the slope. Scan distances range between 400 and 600 m, producing an average point spacing of approximately 10 cm.

Goldpan and White Canyon are located in the Intermontane belt of the Canadian Cordillera. Both sites belong to the Quesnellia volcanic arc terrane, which is composed of Carboniferous to mid-Jurassic volcanic, sedimentary, and plutonic rocks (Struik, 1987; Monger and Nokelberg, 1996). The rock mass at Goldpan is largely massive and is predominantly composed of volcanics of the mid-Cretaceous Kingsvale–Spences Bridge Group. The main rock types of this group are basaltic to andesitic flows interspersed with volcaniclastic sandstones, shales, and conglomerates (Brown, 1981). The comparatively large rock slope of White Canyon belongs to the Mt Lytton plutonic complex (Greig, 1989). Locally, the dominant rock type is an amphibolite grade quartzofeldspathic gneiss. Amphibolite banding is also present and gneissic layering is crosscut throughout the canyon by intrusive phases of gabbro, tonalite, and granodiorite (Brown, 1981). Preferential weathering around these intrusions often results in the formation of vertical rock spires which act as source zones for rockfall.

Mitigative measures at both sites have been installed in response to frequent rockfall activity impacting railway operations. At Goldpan this includes rockfall wire mesh draped over parts of the slope, extensive sections of shotcrete, and four concrete rock sheds. In the eastern half of White Canyon there are two timber rock sheds and one concrete shed, as well as wire mesh rockfall nets and concrete lock blocks adjacent to the track. The western half of White Canyon also has wire mesh rockfall nets and lock block retaining walls, as well as an additional three concrete rock sheds and one timber rock shed. Slide detector fences are present at both sites, comprised of horizontal wires strung between upright telephone poles, and provide warning by switching the track signal to stop when broken.

3 Rockfall events

Five rockfall events were selected from a database of rockfalls which have been identified using change detection at White Canyon West, White Canyon East, and Goldpan since 2012 (Kromer et al., 2015). These events have masses ranging from 2 to 170 m³, exhibiting a combination of structurally controlled failure modes including wedge sliding, planar sliding, toppling, and overhanging blocks. Images of the five slope sections where the rockfall events of interest took place can be seen in Fig. 2.

Figure 2 Photos of five rock slope sections, prior to failure, from the Goldpan (a) and White Canyon (b–e) sites adjacent to the CN Railway. The red regions indicate the source zone for the rockfall events discussed in this paper.

Each event showed material accumulation below the source zone in the change detection results. This is essential information for comparison with our simulation results. Each event is also close to track level, with the highest fall occurring 46 m vertically above the track. These events were selected because the shorter distance from source to a notable accumulation of material presents a simpler trajectory for back analysis and reduces potential confusion in interpreting whether the material gain is due to the selected rockfall or another mass-movement nearby. In each case, change was detected using the Multiscale Model to Model Cloud Comparison (M3C2) point–point distance calculation (Lague et al., 2013) in CloudCompare (2018), between the pre- and post-fall TLS scans. A summary of each of the rockfall events including change detection results, and site photos before and after the event, is provided in Figs. 3–7. Discontinuity and slope angle measurements used for the stereonet failure mode analysis shown in each figure were completed using the Compass tool in CloudCompare.

Figure 3 Site A. Visual overview of the 170 m³ wedge sliding rockfall event detected at the Goldpan site between July and October 2016. The event involved a large triangular slab of weathered and jointed rock mass which failed approximately 35 m above a rock shed which is protecting the track at Goldpan. Site photos (a, b) of the source zone pre- and post-failure are shown. Change detection results of the event are shown (c) with cool colours indicating material loss and warm colours indicating material accumulation. Material from the rockfall accumulated on the bench of the mid-slope gully, as well as on top of the rock shed, with the majority of the material running out over the shed and leaving the slope. The rockfall hull of the event extracted from the pre- and post-failure meshes is shown (d), as well as a stereonet representation of the wedge sliding failure mode.

Figure 4 Site B. Visual overview of the 24 m³ wedge sliding rockfall event detected at the White Canyon East site between May and July 2016. The event involved a large pseudo-cubic block of weathered and jointed rock mass which failed approximately 46 m above track level. Site photos (a, b) of the source zone pre- and post-failure are shown. Change detection results of the event are shown (c) with cool colours indicating material loss and warm colours indicating material accumulation. A small portion of the rockfall volume was retained in the gully leading down to track level, with the majority of the accumulation taking place in the track-side ditch. The rockfall hull of the event extracted from the pre- and post-failure meshes is shown (d), as well as a stereonet representation of the wedge sliding failure mode. It should be noted that parts of the source rock mass also exhibited overhanging sections.

3.1 Event fragmentation

While the selected rockfall events were not observed directly, each rockfall is believed to have been comprised of multiple mobile fragments. This interpretation is based on observations of the fractured state of the source zones pre- and post-fall, as well as the size and distribution of visible rock fragments in the post-fall areas of accumulation, below the source zone and adjacent to the track. An example of this can be seen in Fig. 8 for White Canyon overhanging wedge event. In these photos it is clear that the source rock mass is heavily jointed and that the accumulation of material generated from the event is not a few large blocks but rather a deposit of coarse granular material.

Figure 5 Site C. Visual overview of the 32 m³ flexural toppling event detected at the White Canyon East site between June and August 2016. The event involved a large slab of weathered and heavily jointed rock mass which failed approximately 14 m above track level. Site photos (a, b) of the source zone pre- and post-failure are shown. Change detection results of the event are shown (c) with cool colours indicating material loss and warm colours indicating material accumulation. Change detection results indicate the bulk of material from the rockfall event was retained in the track-side ditch. The rockfall hull of the event extracted from the pre- and post-failure meshes is shown (d), as well as a stereonet representation of the flexural toppling failure mode. It should be noted that parts of the source rock mass also exhibited overhanging sections.

Figure 6 Site D. Schematic overview of the 15 m³ planar sliding rockfall event detected at the White Canyon West site between October 2015 and February 2016. The event involved an irregularly shaped slab of weathered and heavily jointed rock mass which failed approximately 9.5 m above track level. Site photos (a, b) of the source zone pre- and post-failure are shown. Change detection results of the event are shown (c) with cool colours indicating material loss and warm colours indicating material accumulation. Change detection results indicate that the bulk of the material from the event was retained by the track-side ditch. The rockfall hull of the event extracted from the pre- and post-failure meshes is shown (d), as well as a stereonet representation of the planar sliding failure mode.

The presence of multiple mobile fragments means that these events may be classified as fragmental rockfalls (Hungr and Evans, 1988; Hungr et al., 2014). Hungr and Evans (1988) originally proposed that for the case of fragmental falls, the movements of the most mobile fragments in the fall are independent of each other. This suggests that the modelling of these events as a volume of fractured material is unnecessary, and instead single design blocks with a specified average or maximum fragment volume could be used. This is in contrast to the idea that larger volume rock slope failures such as rock avalanches should be modelled as granular flows rather than independent ballistic trajectories (Bourrier et al., 2013). The distinction between these two types of motion is often discussed in relation to the volume of material mobilized as part of the event, with larger volumes (> 10³–10⁴ m³) suggested to show stronger interaction between individual fragments. A discussion of the various volumes and classifications of rock slope failure relevant to the transition between these two styles of motion is presented by Corominas et al. (2017) and suggests that volume thresholds and terminology for these types of events is not yet consistent in the literature.

Figure 7 Site E. Schematic overview of the 8 m³ planar sliding event detected at the White Canyon West site between July and October 2016. The event involved the failure of two distinct sections of the jointed source rock mass, with the upper (6 m³) and lower (2 m³) failures releasing from 10 and 6 m above track level respectively. Site photos (a, b) of the source zone pre- and post-failure are shown. Change detection results of the event are shown (c) with cool colours indicating material loss and warm colours indicating material accumulation. Interpretation of the post-fall imagery indicates that the lower failure was inside the slide path of the upper rockfall event. It is our interpretation that these events likely occurred simultaneously, with the lower failure occurring as a result of impact from the upper failure. Change detection results indicate that the bulk of the material from the two events was retained by the track-side ditch. The rockfall hull of the events, extracted from the pre- and post-failure meshes is shown (d), as well as a stereonet representation of the planar sliding failure mode for the upper, larger event.

Figure 8 Before and after photos of the 24 m³ wedge sliding event from White Canyon East, Site B. The top row of photos shows the pre-fall debris present in the track-side ditch (a), as well as the pre-fall source zone (b). The bottom row shows the post-fall accumulation in the track-side ditch (c), as well as the post-fall source zone (d). From these images we can see the heavily jointed state of the rockfall mass pre-fall and the rockfall back scarp post-fall. A notable increase in the quantity and size of debris fragments in the track-side ditch is visible post-fall.

While rockfall events < 1000 m³, such as those considered in this study, may be described as having limited interaction between mobile fragments, the use of a single design block is not effective in cases where sufficiently large fragmental falls might overwhelm ditches at the base of a slope. A schematic overview of this process can be seen in Fig. 9. Material which builds up in the ditch or other retaining structures may allow trailing rockfall debris to roll out over the newly formed surface. Additionally, a trailing fragment of rock may impact the accumulated pile of debris with enough force to push some fragments out onto the track. Simulation of the entire fragmental rockfall volume at once, as a moving mass of many fragments, allows for important slope-stopping features such as benches, gullies, and ditches to be filled, impacting the runout of trailing rockfall material. Snapshots from a video of a recorded fragmental rockfall, which filled a ditch and impacted a section of railway in western Canada, can be seen in Fig. 10. In this case, the leading portion of the rockfall event was pushed forward by trailing material before it fully came to rest in the ditch. Additionally, subsequent individual rock fragments were able to run out into the track region as a result of the catchment ditch being full.

Figure 9 Schematic of potential fragmental rockfall runout behaviour. Here leading material from the event has filled up a catchment ditch. Trailing rockfall material impacting the back of the deposit has the potential to push the leading material forward and out of the ditch (black arrows) reaching the track area. Additionally, the initial material has created a surface over which trailing rock fragments can move (grey trajectory), reducing the effectiveness of the ditch.

Figure 10 Snapshots from a video of a fragmental rockfall event occurring above a section of railway in western Canada. Elapsed time between frame (1) and frame (6) is approximately 2 s.

The simulation of these events as single block falls would produce unrealistic runout due to interaction with retaining structures like ditches. In addition, the mobility of larger falls is much different than smaller fragments, with significantly higher potential energies at release, as well as larger moments of inertia affecting rotation during runout. At track level the impact energy of a single coherent block would be larger, resulting in an overestimate of the potential for damage to the track. Additionally, a single block simulation may underestimate the spatial extent of track interaction, as individual fragments from a multi-block fall are free to disperse from each other laterally, impacting multiple points on the track at once. From a design perspective, this would influence the type of mitigation selected for the site. The isolated high-energy impacts of a large single block may require a fixed, high-strength system like a mechanically stabilized earth wall or rock shed. In the case of the numerous, small impacts from a fragmented source volume, wire mesh draping to divert the fragments into the ditch may suffice.

4 Runout simulations

For each of the selected rockfall events from our field sites, we can set up a simulation in the Unity game engine. The main input parameters for simulation setup are

• coefficient of friction,

• coefficient of restitution,

• coefficient of viscoplastic ground drag,

• slope geometry,

• rockfall geometry, and

• rockfall release position and orientation.

Collision detection, contact resolution, and the application of material-specific impact coefficients are handled by the built-in physics system in the Unity game engine. The resolution of linear and angular velocity per contact along a 3-D fragment surface is dependent on the simple friction and restitution equations shown below. The friction coefficient (μ) is used during collision to determine the proportion of the normal force (N) between the simulated rockfall and slope surfaces which is to be used for frictional resistance to movement. The friction force (F_f) acts in the opposite direction of the downslope component of the incoming velocity vector. The force calculated per contact is applied over the simulation time step (as an impulse) and used in the solution of the new translational and rotational velocities which occur as a result of collision. In the Unity environment, surface collisions are resolved into two contact anchor points, resulting in the frictional force applied being twice as strong. Therefore, the effective friction angle (ϕ) is equal to the arctangent of double the Unity friction coefficient. The equation for the frictional force applied per contact anchor point, as well as the conversion between the coefficient of friction and friction angle, is shown below. Unlike the coefficient of dynamic friction, which is applied to the tangential component of velocity, the coefficient of restitution operates only in the normal direction. The coefficient is equal to the ratio of the per-contact relative velocity between the two colliding objects before and after the collision. The third parameter, viscoplastic ground (vpg) drag, is assigned to regions of the slope which are expected to have deformable surface materials. The equation used here to model the force applied due to vpg drag is based on the ground drag equation from the RAMMS::ROCKFALL model and is given below (Bartelt et al. 2016). The drag force (F_v) is proportional to the mass of the rock block and the square of the block velocity at impact ($V_{\mathrm{s}}^{\mathrm{2}}$). Essentially, impacts with higher kinetic energy (heavier, faster moving blocks) are expected to result in a larger amount of dampening as they cause more surface deformation of the slope material, resulting in an increased dragging effect on the rockfall block.

$$\begin{array}{}\text{(1)}& {\displaystyle }{F}_{\mathrm{f}}=\mathit{\mu }\ast N\text{(2)}& {\displaystyle }\mathit{\varphi }=\arctan \left(\mathrm{2}\mathit{\mu }\right)\text{(3)}& {\displaystyle }{R}_{\mathrm{v}}={\displaystyle \frac{v_{\mathrm{incident}}}{v_{\mathrm{rebound}}}}\text{(4)}& {\displaystyle }{F}_{\mathrm{v}}=-{\displaystyle \frac{m}{\mathrm{2}}}{C}_{\mathrm{v}}{V}_{\mathrm{s}}^{\mathrm{2}}\end{array}$$

A more detailed description of the material coefficients and their influence on the rotational and translational velocity of simulated rock fragments, along with a parametric investigation into the effects of simulated block shape and release orientation, can be found in Sala (2018).

The geometry used to build simulations is comprised of slope and rockfall elements. The slope geometry is a 3-D mesh produced from the pre-fall TLS scan of the slope, generated using Poisson reconstruction in CloudCompare. The rockfall geometry is a 3-D mesh produced from the pre- and post-fall TLS scans of the source region, generated using Poisson reconstruction and connected into a single hull in Blender (Blender Foundation, 2017). By extracting the rockfall source volume directly from the same TLS data used to generate the slope geometry, the release position and orientation of the source volume relative to the slope surface is well defined. This 3-D hull is then segmented into a collection of convex hull fragments, using Voronoi fracturing in Blender. This method subdivides the source rockfall volume into a specified number of fragments selected by the user, with control over the size and shape of fragments produced. An example of various Voronoi fracturing results, including recursive fracturing, and fracturing with preferential shape constraints can be seen in Fig. 11. The end result is a simulated rockfall source volume which reflects the shape and fractured nature of the rockfall event detected from the field.

Figure 11 Showing the results of different Voronoi fracturing methods. Recursive fracturing in (b) results in the additional fracturing of the initial convex polyhedra in (a). Shape constraints in (c) result in the elongation of the convex polyhedra, creating a more sheet-like fracture network.

A detailed explanation of the steps necessary to move from field remote sensing data to a functioning simulation in Unity, including the workflow from CloudCompare to Blender to Unity, can be found in Sala (2018). One of the key strengths of this modelling technique is its ability to simulate the collision between multiple moving rigid bodies at once. Support for multi-body collisions allows us to simulate the disaggregation of heavily jointed source rockfall volumes. This enables us to study the types of rockfall runout expected to occur at our field sites, taking into consideration the implications of fragmental runout discussed above. Snapshots of an example fragmental rockfall simulation for the 170 m³ Goldpan event (Site A) can be seen in Fig. 12.

When developing or attempting to calibrate any new simulation technique, it is necessary to compare simulation outputs to real observations of the phenomena you are modelling in order to test how well the simulation performs. Where direct observation of rockfall events is possible, such as in rockfall drop tests, (e.g. Ushiro et al., 2006; Vick et al., 2015; Volkwein et al., 2018), this information could include translational and rotational velocity, mapped runout trajectories, or pass height above the slope. For the five rockfall cases selected from our field sites, this information is not available as the events were not observed directly. Instead the location and magnitude of change present in the change detection results for each event serve as the basis for simulation comparison.

Simulations of each of the five selected rockfall events were run and change maps of the simulation results were produced. In order to generate change detection maps for the simulations, the positions of the post-simulation fragments were merged with the slope geometry to create a post-simulation 3-D mesh. Similarly, the fragmented rockfall geometry, pre-simulation, was merged with the slope geometry to create a pre-simulation 3-D mesh. Each mesh was then converted into a point cloud in Blender. The conversion of the mesh data to a point cloud allows us to use the M3C2 point–point distance calculation for both the actual and simulated rockfall events. The transition from simulation mesh to point cloud is done by selecting only visible vertices in the mesh, using a similar vantage point to the direction of the original TLS scan. The selection of only visible vertices in the mesh is necessary as the 3-D rockfall fragments, pre- and post-simulation, have vertices across their entire surfaces, both front and back facing. This results in multiple layers of points in the pre-fall source zone and post-fall accumulation zone, adversely affecting point cloud normal calculations and producing errors in the distance measurements between what should be distinct pre- and post-fall surfaces. The extraction of visible vertices from the mesh geometry, using a specified vantage point, ensures that only a single layer of points is present in the simulated point cloud.

Figure 12 Screenshots of a fragmental rockfall simulation (1000 fragments) for the 170 m³ rockfall event at Goldpan.

Figure 13 A comparison of the actual and simulated change detection results from the 170 m³ wedge sliding event at Goldpan, Site A. The simulation produced a good comparison with the actual change results, yielding a similar spatial distribution, and magnitude of change. In both cases accumulation takes place on the mid-slope bench and on top of the rock shed, with the majority of the material running out over the shed and off the slope. The percentage of initial source volume retained at various locations on the actual and simulated slope is displayed and shows a strong agreement. Areas of additional loss in the actual change detection, on top of the rock shed and at the mid-slope bench, are not captured in the simulation results. The simulation technique uses a fixed, rigid slope surface and is not able to capture smaller slope failures which may have occurred as a result of talus impacts during the initial 170 m³ event.

Figure 14 A comparison of the actual and simulated change detection results from the 24 m³ wedge sliding event from White Canyon East, Site B. The simulation produced a good comparison with the actual change results, yielding a similar spatial distribution, and magnitude of change. The percentage of initial source volume retained at various locations on the actual and simulated slope is shown. In this case the majority of the material in both the actual and simulated analysis ends up in the track-side ditch, with a smaller component of the volume stopped inside the gully. In the simulated change a small number of fragments do reach track level, representing approximately 2 % of the simulated volume.

The results of the change detection analysis for each simulated rockfall event are shown in Figs. 13–17. In each case the game-engine-based simulation prototype was able to produce rockfall runout which compared well in terms of the location and magnitude of the measured actual change. The percentages of source rockfall volume retained at different locations on the slope is shown. Volume measurements using the actual event point cloud data were completed using the 2.5-D Volume tool in CloudCompare. Simulated volume percentages were determined using the 3-D mesh data for each of the convex hull fragments simulated in Unity 3D. For each event, the simulated volume percentages for each of the slope locations was within 5 % of the actual percentage, with the majority of the rockfall material ending up in the track-side ditches for every rockfall event except the Goldpan wedge slide. In the Goldpan case, Site A, we see two distinct areas of accumulation, on the mid-slope bench and on top of the rock shed. In this case the simulation results produced similar accumulations in both locations, with the majority of the material leaving the slope over the edge of the rock shed.

Figure 15 A comparison of the actual and simulated change detection results from the 32 m³ flexural toppling event from White Canyon East, Site C. The simulation produced a good comparison with the actual change detection results, yielding a similar spatial distribution, and magnitude of change. The percentage of initial source volume retained at various locations on the actual and simulated slope is shown. The majority of the material in both cases is retained in the track-side ditch. In the simulated case a small amount of the simulated volume comes to rest on the slope (0.5 %) and runs out into the track area (1.5 %).

Figure 16 A comparison of the actual and simulated change detection results from the 15 m³ planar sliding event from White Canyon West, Site D. The simulation produced a good comparison with the actual change detection results, yielding a similar spatial distribution, and magnitude of change. The percentage of initial source volume retained at various locations on the actual and simulated slope is shown. The majority of the material in both cases is retained in the track-side ditch. In the simulated case, 1 % of the volume reaches the track area.

Figure 17 A comparison of the actual and simulated change detection results from the 8 m³ planar sliding event from White Canyon West, Site E. The simulation produced a good comparison with the actual change detection results, yielding a similar spatial distribution and magnitude of change. The percentage of initial source volume retained at various locations on the actual and simulated slope is shown. The majority of the material in both cases is retained in the track-side ditch. In the simulated case a small amount of the simulated volume comes to rest on the slope (1.5 %) and runs out into the track area (1.5 %).

4.1 Simulation with multiple parameterizations

In each of the above change-based simulation comparisons, a single model iteration which produced a positive comparison to the observed change was shown. The input parameters used to produce each comparison are not identical, with variation between the five sites. For a single simulation, a parameter set consisting of coefficients of dynamic friction, restitution, and viscoplastic ground drag that best simulated the actual results was chosen. The parameter sets used for each of the five site comparisons are shown in Table 1.

Table 1 The material coefficients used for each of the five rockfall event simulations. The fragmentation levels selected for the 15 simulation iterations run for the White Canyon East wedge sliding event are also included. Fric indicates Friction and Rest indicates Restitution.

While a single parameter set can be selected, producing results which align with the material accumulations seen in the change detection, we cannot be certain that this parameter set will best reflect conditions at any given site in the area. The friction, restitution, and damping coefficients used are a simplification of a suite of interconnected and complex processes which govern the loss of energy during collision between the rockfall and slope. For example, variations in moisture content for a section of soil affect slope deformation and therefore the transfer of energy which takes place during collision (e.g. Vick, 2015). From a forward-modelling perspective, while advances in rock slope monitoring have shown that we are now able to detect pre-failure deformation for unstable sections of slopes (e.g. Royan et al., 2013; Kromer et al., 2015, 2017), in some cases providing time-to-failure estimates, we are not able to predict the exact time at which a failure will occur. Therefore, in terms of material parameters, we cannot know the exact state of the slope when runout will take place. For this reason, it makes sense that a range of potential material coefficients be used to capture the variability in these parameters.

Figure 18 Results of the 15 simulations run for the 24 m³ wedge sliding event at White Canyon East, Site B. Each point represents the end point of a simulated fragment in one of the 15 simulations. The envelope of red points indicates the runout of 95 % of the simulated volume, with the purple and blue portions representing envelopes containing the additional 4 % and 1 % of the material, respectively.

Similarly, for forward modelling of events based on pre-failure deformation or precursor rockfall, the exact volume of a failure would not be known prior to the event. While volume estimates based on structure in the source rock mass are possible (e.g. Lambert et al., 2012; e.g. Salvini et al., 2013), they require assumptions about joint persistence and shape. As a result, it would be advisable to estimate a range of potential volumes using the least and most conservative persistence estimates. Additionally, while we may be able to identify potential failure events, and even estimate block volume from visible structure, we won't know the degree of fracturing expected as part of the disaggregation of the source volume. With this in mind, it is important that we also use a range of fragmentation values as well, in this case termed coarse, medium, and fine. The use of a suite of material coefficients, volumes, and levels of fragmentation allows us to produce an envelope of potential rockfall runout rather than a single deterministic result.

An example suite of simulations has been produced for the 24 m³ overhanging wedge failure at White Canyon East, Site B. Varying source volumes were not used in this case as the rockfall has already occurred, and the source volume is therefore known. The material coefficients used were based on the five parameter sets from Table 1, each of which produced a positive comparison at one of the five sites. These five material parameterizations were each run for three levels of fragmentation (250 fragments, 500 fragments, 1000 fragments) produced using Voronoi fracturing in Blender. The distribution of potential rockfall runout for the event was produced for each of these 15 simulations. Each point shown on the slope surface in Fig. 18 is the end position of a fragment from one of the 15 simulations. The red points represent the runout of 95 % of the simulated volume. This boundary illustrates that the majority of the material from the event, based on 15 different simulation parameterizations, is retained by the track-side ditch. This result compares well with the actual change detection from the event, in which all of the rockfall material appears to have come to rest in the gully above track or in the track-side ditch. The 99 % (purple) and 100 % (blue) runout boundaries are also shown, illustrating that a portion of the additional 5 % of simulated material does come to rest in the track area.

4.2 Volume comparisons and ditch design

From an engineering design perspective, one application of this simulation technique is the assessment of mitigation performance. The simulation of a fragmental rockfall event as a single large block overestimates potential impact energies and underestimates the spatial distribution of impacts possible when rockfall volumes run out as hundreds or thousands of individual, interacting fragments. The simulation of rockfall volumes as a collection of free-moving fragments allows us to produce more realistic material accumulations in retaining ditches adjacent to the track.

Using the 15 simulations run for the 24 m³ White Canyon East wedge sliding event, Site B, we can analyze the distribution of accumulated material along the slope surface. Figure 19 shows the same end point locations as in Fig. 18 but in section view. Also included is a plot displaying the volume of material accumulated along the slope surface. From this distribution we observe a notable peak in volume occurring at the location of the track-side ditch. This style of volume accumulation curve allows us to visualize the effectiveness of countermeasures at retaining material. In preparation for an expected rockfall event, or during the construction of new ditches along a section of railway or highway, this analysis could be used to evaluate the effectiveness of different ditch shapes or the use of simple retaining structures like lock blocks. Additionally, the retention capacity and the effectiveness of countermeasures, as material progressively accumulates, could be tested.

Figure 19 A section view of the 15 simulation iterations run for the 24 m³ wedge sliding event at White Canyon East, Site B. (a) shows the accumulation of simulated fragments across the section. The notable peak visible at 270 m is in alignment with the location of the track-side catchment ditch, as shown in (b).

In the initial 15 simulations run for the White Canyon East wedge sliding event, the pre-failure slope surface was used in order to examine the match between the simulated and actual conditions on the slope during runout. We can also re-run these 15 simulations using the post-failure slope surface. In this case, the ditch area is full of debris and large rockfall fragments from the 24 m³ event, representing a ditch near full retention capacity. The simulation results using the post-failure surface conditions can be seen in Fig. 20, compared to the initial pre-failure analysis. From this we can see that the 95 % volume-runout boundary has moved closer to the track for the full ditch scenario, with some end points now touching the track boundary. A comparison of the two volume accumulation profiles for these different ditch scenarios can be seen in Fig. 21. Two vertical lines indicate the location of the 95 % volume-runout boundary for the two curves, with the red line indicating the start of the track area. The use of the post-failure, full ditch slope model resulted in the 95 % boundary moving forward 1.2 m, illustrating a decrease in the effectiveness of the ditch at stopping material when full. From a back-analysis perspective this effect on runout illustrates the importance of using the true pre-failure slope conditions for our simulation comparisons. From a forward analysis perspective, it may not be possible to know the state of the ditch during the runout event, and therefore different ditch conditions should be simulated in order to assess the range of possible outcomes when evaluating countermeasure options.

Figure 20 Comparison of the 15 simulation iterations of the 24 m³ wedge sliding event at Site B using the pre-fall (empty) and post-fall (full) ditch geometries. In the full ditch scenario the 95 % and 99 % runout boundaries have shifted further forward into the track area, indicating a decreased effectiveness in the ditch at stopping material from reaching the track.

Figure 21 The distribution of simulated volume for the pre-fall and post-fall ditch geometries is shown. The vertical grey and black lines indicate the 95 % volume-runout boundaries for the two cases. Simulation using the post-fall (full) ditch geometry results in the 95 % volume-runout boundary shifting 1.2 m closer to the inside rail.

A comparison of the percentage of rockfall volume retained in different areas of the model for the empty and full ditch scenarios is shown in Fig. 22. Simulation using the post-fall ditch geometry results in almost twice as much material reaching track level. Additionally, we can look at the size distribution of the fragments which end up in the track area, as shown in Fig. 22 for the full ditch scenario. The use of volume-based runout envelopes takes into consideration both the size and number of fragments which reach the track. Different-sized fragments have different implications from a hazard management perspective. For example, in the CN Rockfall Hazard Risk Assessment framework (Abbot et al., 1998a, b), rockfall fragments with maximum dimensions of 0.3–1 m are noted as the sizes most likely to cause derailments due to their potential to get wedged underneath a train car. Of the 11 % of simulated fragments which reach the track in the 15 simulations of the full ditch scenario, approximately 57 % of them fall within the 0.3–1 m size range.

Figure 22 A schematic overview of the volume percentages stopping at different locations for the 15 simulation iterations using the empty and full ditch geometries. Simulation using the full ditch geometry resulted in nearly twice as much material reaching track level. The distribution of the grain sizes which reached track level in the full ditch scenario is also shown.

5 Limitations

The objective of this work was to demonstrate the ability of our rockfall simulation workflow to model complex rockfall source geometries comprised of many moving fragments. While this process is a step forward in the simulation of fragmental rockfall events, a discussion of current limitations is important. The two main limitations of this work are fragmentation timing and the distribution and size of rockfall fragments used.

5.1 Fragmentation timing

The events modelled in this paper were not observed directly, and therefore it is not possible, using the available TLS data and site photographs, to understand at what point fragmentation took place during the event. While the post-fall accumulations visible in the photos and change detection indicate that the final state of the source volumes were deposits of coarse rock fragments, the transition from coherent rock mass to thousands of fragments could have occurred immediately at the time of detachment, or as a result of impact during the first or subsequent collisions between the source volume and slope surfaces.

The current methodology simulates the release of rock from the source zone as a gravity-induced, rapid unravelling of highly jointed material (Fig. 23a). It is possible that this unravelling of jointed rock mass could have instead occurred slowly over time as small isolated events (Fig. 23b). Alternatively, the separation of the source mass into individual blocks could have occurred as a result of loading during impact (Fig. 23c). At the same time, the decision to model the events as a rapid disaggregation of jointed rock was not arbitrary. This style of fragmentation is akin to the disaggregation-without-breakage rockfall mechanism described by Ruiz-Carulla et al. (2016a), in which the rockfall body separates along discontinuities present in the source rock mass, with no breakage of individual, discontinuity-bounded blocks. This appears to be the dominant behaviour in the railway rock slope video shown in Fig. 10. This type of rock mass behaviour is also observed in low-stress underground excavations, in which structurally controlled failures of highly jointed or crushed rock will unravel from excavation ceilings under the influence of gravity (Palmström, 1995).

Figure 23 Comparison of potential timing of rock mass fragmentation, which could have taken place for the events discussed. Red blocks indicate material which breaks away from the initial source volume at a given time.

The choice to model the events as a rapid unravelling of the entire source volume (Fig. 22a) was also inherently influenced by the current capabilities of the modelling technique. In Sala (2018), the possibility to connect simulated fragments in the source mass using strength-based object-to-object connections was briefly discussed. Each simulated connection between objects can be assigned a strength value, with the connection broken when this threshold is exceeded. The use of these connections could facilitate future simulations in which the source mass does not disaggregate immediately but instead as a result of force from impact. While this is one potential method for simulating break-up of the source volume during runout, field data to support the calibration of these connections are not currently available.

5.2 Fragment size distribution

One of the goals of modelling rockfall source volumes as a collection of rock fragments was to better capture the size of mobile fragments involved in the rockfall events rather than using a single large block. While this effect was achieved using our modelling workflow, it should be noted that an evaluation of actual versus simulated grain size distribution has not been completed.

At this time, Voronoi fracturing is used for source volume fragmentation. In 3-D this method partitions the source mesh into disjointed convex polyhedra, based on seed points distributed in the mesh volume (Ledoux, 2007). Several parameters can be set in the fracturing algorithm including the recursive partitioning of larger fragments and control over fragment shape. While the use of this method is capable of rapidly creating fragmented 3-D volumes, ready for use inside the Unity game engine, the fracture network produced is not based on actual discontinuities in the rock mass. For rock masses which are less jointed, discontinuity planes can be used in Blender to separate the source volume into blocks based on discrete joints mapped from available imagery or 3-D models of the slope. In the case of the White Canyon and Goldpan events, where hundreds to thousands of fragments were visible in the post fall accumulations, the mapping of discrete joint planes was not practical. For these cases the use of Voronoi fracture networks, instead of discrete joint planes, permitted us to model the sheer number of fragments present in the rockfall deposits.

The size of the fragments deposited during the five rockfall events is clearly not unimodal, with a variety of grain sizes present. Grain size distributions for these events were not created due to a lack of adequate data. Point spacing in the TLS data is too sparse to be able to isolate individual grains. Additionally, vantage points from our scan sites often result in notable occlusions at the track level, resulting in missing data in the post-fall debris piles. Photos of the events do capture the post-fall debris, but due to the oblique angle of the images, and homogenous colour of the material, accurate image-based grain segmentation is challenging.

A select few fragments measured for the White Canyon East 24 m³ case, Site B, show approximate sizes ranging from 1.6 to < 0.05 m. Particles smaller than 5 cm are clearly present but cannot be measured as the images are too noisy at this scale due to resolution limitations. Size ranges for the coarsest models (250 fragments) and finest models (1000 fragments) of the 24 m³ event were 0.09–0.8 m and 0.02–0.6 m respectively. While we do not have a full grain size distribution for the actual event, we can see from these measurements that the grain sizes produced using the current Voronoi fracture parameters underestimate the largest block sizes and result in an overall smaller range of values. In terms of the smaller grains, the 2 cm minimum size of the 1000 fragment model is sufficient to capture the lower size limit of the 24 m³ event. Ultimately more investigation needs to be done into fine-tuning the parameters of the Voronoi fracturing process in order to produce grain sizes with a broader range of values and that more accurately reflect the sizes produced by events in our study area. In order to do this, complete grain size distributions for events in our rockfall database, and from other events in the literature, will need to be generated.

6 Conclusions

Five rockfall simulations were run based on TLS change detection inputs from events identified at our field sites in south-central British Columbia. The simulations utilized high-resolution 3-D meshes generated from 6–10 cm TLS point clouds, and complex rockfall volumes, extracted from pre- and post-event scans. The fracturing of source rockfall volumes enabled the modelling of fragmental rockfall runout, taking into consideration the interaction of fragments with each other as well as slope-stopping features such as gullies and ditches. The results of this work demonstrate the ability of our rockfall modelling prototype to produce runout material accumulations which agree well with observed change from the actual events. With no direct observations of rockfall runout, or on-site measurement, the conversion of post-simulation mesh data in order to produce simulated change detection presents a novel way to perform rockfall runout comparisons in the absence of other data types.

The application of the technique to mitigation design was also discussed, demonstrating the potential of the model to be used for investigating the effectiveness of rockfall catchment ditches. This type of analysis, visualizing the accumulation of rockfall runout volume downslope, could also be applied to the assessment of other countermeasure options such as lock block retaining walls or rock sheds. The ability to perform these types of analyses is limited in industry standard rockfall modelling programs, as they are often able to simulate only a single moving block at a time. While these models are very useful for single block falls, and regional runout assessments, they are unable to capture the build-up of rockfall material on the slope, which can affect the path of subsequent rockfall fragments and the effectiveness of retaining structures. The simulation of these falls as multiple moving bodies in our technique is able to capture this interaction between rockfall material and the finite capacity of the structures used for rockfall protection.