Emergent characteristics of rockfall inventories captured at a regional scale

High‐resolution rockfall inventories captured at a regional scale are scarce. This is partly owing to difficulties in measuring the range of possible rockfall volumes with sufficient accuracy and completeness, and at a scale exceeding the influence of localized controls. This paucity of data restricts our ability to abstract patterns of erosion, identify long‐term changes in behaviour and assess how rockfalls respond to changes in rock mass structural and environmental conditions. We have addressed this by developing a workflow that is tailored to monitoring rockfalls and the resulting cliff retreat continuously (in space), in three‐dimensional (3D) and over large spatial scales (>104 m). We tested our approach by analysing rockfall activity along 20.5 km of coastal cliffs in North Yorkshire (UK), in what we understand to be the first multi‐temporal detection of rockfalls at a regional scale. We show that rockfall magnitude–frequency relationships, which often underpin predictive models of erosion, are highly sensitive to the spatial extent of monitoring. Variations in rockfall shape with volume also imply a systemic shift in the underlying mechanisms of detachment with scale, leading us to question the validity of applying a single probabilistic model to the full range of rockfalls observed here. Finally, our data emphasize the importance of cliff retreat as an episodic process. Going forwards, there will a pressing need to understand and model the erosional response of such coastlines to rising global sea levels as well as projected changes to winds, tides, wave climates, precipitation and storm events. The methodologies and data presented here are fundamental to achieving this, marking a step‐change in our ability to understand the competing effects of different processes in determining the magnitude and frequency of rockfall activity and ultimately meaning that we are better placed to investigate relationships between process and form/erosion at critical, regional scales. © 2020 The Authors. Earth Surface Processes and Landforms published by John Wiley & Sons Ltd


Introduction
Rockfalls are a frequent process initiated when rock blocks become detached from a rock mass under the influence of gravity (Selby, 1982). They exert a first-order control on the rate of rock wall retreat on mountain slopes and on rock cliffs (Moore et al., 2009). Their volumes typically range from~10 -2 to 10 2 m 3 , but in some cases they have been known to reach 10 5 m 3 (for example, Wieczorek et al., 1998;Stock et al., 2012). Rockfall activity is also a chronic hazard (Evans and Hungr, 1993;Guzzetti et al., 2003;Wieczorek et al., 2008), often posing significant risks to transportation corridors Katz et al., 2011;Blais-Stevens et al., 2012;Michoud et al., 2012;Ansari et al., 2014), pipelines Couture et al., 2010), and to areas beneath (sea) cliffs Marques et al., 2013). Rockfall activity has been monitored extensively in these settings, and in some cases this monitoring has been used to provide hazard and risk forecasting Stock et al., 2012;Royán et al., 2013). Much of the research into rockfall activity has drawn upon datasets covering relatively small length scales (~10 1 -10 2 m), often defined by instrument capabilities rather than by any scientific rationale, in order to infer wider behaviour. Two important consequences arise: first, that it is often difficult to capture a statistically significant number of the largest rockfalls while also capturing the smallest events in a single inventory; and second, that data from a single site are likely to reflect site-specific conditions. For example, in a 10-month dataset captured across an 8500m 2 rock face, only 10 of~1.8 × 10 5 rockfall events were >1.0m 3 and none exceeded 10.0 m 3 (Williams et al., 2018), whereas rockfall with volumes up to 2.5 × 10 3 m 3 are known to have occurred at this site in previous years (Rosser et al., 2005a).
Rock slope evolution is not uniquely governed by large, infrequent events; it instead reflects a continuum of change where failures can also be small in magnitude and variable in frequency over large areas (Lim et al., 2010). While the smallest events have been observed to occur at high frequencies, often resulting in near-continuous mass wasting and therefore representing a chronic hazard in some areas, the scars and debris of catastrophic events tend to remain evident in the landscape for longer, controlling long-term rates of landform and landscape evolution (Hovius and Stark, 2006). Monitoring rockfall activity across the full range of potential scales poses a number of challenges, as data collected at scales sufficient to capture the largest events are often not sufficiently high in resolution to capture the smallest events. Upscaling detailed monitoring of rock slopes is difficult, both in terms of capturing and processing data, as the topographic complexity of the area monitored inevitably increases with scale. This can mean moving from a single, near-planar rock face to a more complex series of hillslopes with variable lithology, geometry, and structure. Upscaling is, therefore, very rarely a case of applying local approaches more extensively. Similarly, increasing the likelihood of capturing a large event by prolonging the period of monitoring can be prohibitively costly, and, where monitoring intervals do increase, the data captured are inevitably subject to rockfall coalescence and superimposition, which decreases the likelihood of detecting small events (Williams et al., 2018). Whether or not ergodic reasoning (space-for-time substitution) can be applied also remains to be tested. This has implications for our understanding of rock slope failure. A fundamental uncertainty, for example, is whether monitoring 1.0km 2 of cliff face over oneyear would generate a rockfall inventory that is statistically comparable to that captured from 0.1km 2 over 10years, from a set of apparently uniform cliffs. This is unlikely where the timescales of path-dependent behaviour in rockfall evolution, via brittle fracture growth (Kemeny, 2003) and progressive failure , and over longer timescales via changes to slope-profile form and/or post-glacial or post-incision relaxation (Cordes et al., 2013;Messenzehl and Dikau, 2017), are commensurate with or exceed those of most monitoring campaigns.
One empirical/statistical approach to compensate for the difficulty in capturing regional-scale observations is to use the power law behaviour in rockfall magnitude and frequency to upscale, in both time and space, and model future rockfalls and hence cliff erosion, assuming that what is monitored at a small scale is more widely representative (Lim et al., 2010). These approaches have inherent assumptions and limitations that restrict their application, including (1) that they rely on extrapolating a non-biased, complete portion of an inventory to predict both larger and smaller volume frequencies, (2) the need to apply power laws within limits, in order to avoid generating biased scaling coefficients (Barlow et al., 2012), (3) the implicit assumption that a single underlying mechanism, and hence a single form of power law behaviour, transcends all scales of events under investigation (Brunetti et al., 2009), and that extrinsic controls are essentially constant, and (4) that all rockfalls in an inventory are statistically independent of one another, although it is known that rockfalls exhibit some degree of spatial and/or temporal path-dependency (Rohmer and Dewez, 2015). Fundamentally, this approach loses any site specificity, generating only broad rockfall magnitude probabilities rather than an estimation of what could or will happen at an individual location.
An increasingly viable alternative, enabled by rapid advances in (mobile) three-dimensional (3D) data capture on near-vertical surfaces, is to monitor rockfalls over a larger area while retaining a high spatial resolution (Lato et al., 2009). To achieve this, some airborne LiDAR (light detection and ranging) systems are capable of collecting from oblique as well as vertical view angles, permitting the capture of point cloud data both on near-vertical surfaces and over much larger extents. However, the volume and quality of data that can be collected using airborne LiDAR present their own unique challenges. Such data require methods that are able to retain the 3D character of the data while also being able to measure rockfall volumes that can span over six orders of magnitude, and over spatial extents that can exceed~10 6 m 2 . These settings could include, but are not limited to, a length of coastline (for example, Teixeira, 2006;Rosser et al., 2007;Marques, 2008;Lim et al., 2010;Young et al., 2011;Barlow et al., 2012;Rohmer and Dewez, 2013;Kuhn and Prüfer, 2014;Williams et al., 2018Williams et al., , 2019, cut slopes along a transport corridor (for example, Bunce et al., 1997;Hungr et al., 1999;van Veen et al., 2017), or on montane, alpine or arctic rock walls (for example, Dussauge-Peisser et al., 2002;Malamud et al., 2004;Santana et al., 2012;Messenzehl and Dikau, 2017). Without data obtained at these scales, it remains difficult to assess whether rockfalls are truly scale invariant across all the possible volumes of a given distribution, to put limits on modelled power laws of rockfall magnitude and frequency, and therefore to test whether rockfalls can be considered stochastic phenomena. It therefore remains difficult to define the extent of rock face required to capture sufficiently 'representative' rockfall behaviour a priori.
To address these challenges, we present an approach to monitor both small and large rockfall over large areas of near-vertical rock face. We monitored changes along a 20.5 km stretch of coastal rock cliffs (~40m height) using a unique, helicopter-based LiDAR system that captures data on surfaces both beneath and oblique to the flight path, retaining high 3D positional accuracy and geometrical detail on surfaces of all angles. The resulting data were used in the development and application of a fully 3D workflow for detecting and characterizing rockfall activity, with the aim of deriving a detailed rockfall inventory and assessment of coastal cliff erosion. The resulting inventory comprised >58000 rockfalls recorded over four repeat, approximately annual, surveys. We use this dataset to explore the characteristics of rockfall scaling behaviours that emerge at a regional scale, and the implications of our findings for monitoring and modelling rockfall activity.

Study Site
Our study area extends along the North Yorkshire coast (UK) for 24km between Skinningrove and Sandsend, of which 20.5km is near-vertical rock face ( Figure 1). The cliffs are cut into complex, near-horizontally interbedded Lower Jurassic mudstones, shales, siltstones, limestones and sandstones, much of which is capped by silty glacial tills (Powell, 2010). In places, headlands rise up to form cliffs reaching 150m, while embayments are often characterized by sandy beaches backed by low (<30m) cliffs. Much of the coastline is fronted by a gently sloping (<2°) foreshore platform that extends more than 300m seaward in places and is fully exposed when high atmospheric pressure systems coincide with the lowest astronomical tides. The tidal range along the coast is largely macro-tidal, experiencing two daily tides that cycle between spring and neap highs over a range of~6m. When high spring tides coincide with storm events and high swells, the vertical reach of the tide up the cliffs can exceed 4.2m (Rosser et al., 2013).
The North Yorkshire coast is shaped by marine, subaerial, and anthropogenic (mining exploitation and management intervention) processes, leading to complex patterns of morphological change as the coastline adjusts and retreats (Lim, 2014). The rate of cliff erosion at Staithes (Figure 1b) has been estimated at approximately 0.05myr -1 over the last century, based on the analysis of cliff top position from historic maps and photographs (Agar, 1960). This has been derived from a calculated retreat rate of approximately 0.04myr -1 for headlands and 0.07myr -1 for embayments. However, when considered at this scale and monitoring interval, this rate falls well beneath the minimum achievable mapping precision at any given point. Measurements of coastal erosion using historic 2774 J. BENJAMIN ET AL.
map data neglect processes of undercutting and small scale, iterative failures of localized sections of the cliff face, instead focussing on the overall recession of the cliff top or toe (see Lim et al., 2005;Rosser et al., 2005a). The practical implications of this are that erosion rates determined using these approaches are associated with very high levels of uncertainty.

Methods
We detail a robust workflow for detecting and characterizing regional-scale changes in cliff morphology due to rockfall in full 3D. We applied this workflow to four high-resolution point clouds derived from approximately annual airborne LiDAR surveys of coastline in North Yorkshire, UK.

Data acquisition and processing
Four airborne LiDAR surveys (Table 1) were captured at approximately equal intervals (~320days) between August 2014 and March 2017 using a helicopter equipped with a mobile mapping system. The system (ROBIN + WINGS) comprised a RIEGL VQ-450 or RIEGL VUX-1 laser scanner coupled with an IGI AeroControl III navigation system, which combined a global positioning system (GPS) receiver with an Inertial Measurement Unit (IMU-IIe) in order to measure the position and attitude (pitch, roll and yaw) of the helicopter. Uniquely, the system used a high pulse repetition rate (up to 550kHz) near-infrared laser and a rotating mirror to return a 180°swath of data normal to the flight track (RIEGL, 2014). This enabled both the cliff top, foreshore and near-vertical cliff faces to be scanned simultaneously with the same ground sample distance. A downward-looking 36.3 MP Nikon D-800 camera with a 24mm lens was also used to capture optical imagery. Ground GPS data were recorded at one sample per second by 12 Leica 1200 GPS receivers (Table S1, Supporting Information). The GPS antennas were mounted on tripods placed over targets in 13 locations ( Figure S1, Supporting Information).
Overlapping flight lines were flown in each survey epoch to increase point density, where intentional variations in the attitude of the aircraft gave multiple LiDAR incidence angles onto FIGURE 1. Map of the North Yorkshire coast, showing the cliff areas monitored in this research. The total length of cliff face monitored is approximately 20.5km. A foreshore platform, which in places extends more than 300m seaward, fringes much of the coastline and is fully exposed at the lowest astronomical tides. Light grey dashed lines indicate the 38 blocks that the data were divided into for processing. the same area of the cliff face, minimizing any occlusions due to surface roughness. The system was deployed at a mean flying height of~100m above ground, providing a spatial resolution of~0.01m for the optical imagery ( Figure S1) and a mean point density from~30 to 80 points m -2 , irrespective of surface angle. Each survey was registered to a global coordinate system using the ground control points collected. Point cloud data for each of the four surveys were clipped in plan-view to retain the extent of the vertical cliff faces. For each survey, the data were then divided into 38 blocks, with each block 500m in width at the centreline (Figures 1, S1). To improve the quality of the alignment between surveys, the data were further aligned on a block-by-block basis using an iterative closest point alignment, with a maximum permissible 3D registration error of 0.10m.

Change detection
Changes were measured between four monitoring epochs; 2014-2015, 2015-2016, 2016-2017 and 2014-2017. For each epoch, the point clouds n 1 and n 2 were processed block-by-block using the M3C2 algorithm (Lague et al., 2013). M3C2 was implemented in three stages: (1) the estimation of 3D surface normals (i.e. the direction in which the surface points) within a user-defined normal diameter, D, (2) calculation of the average position of each point cloud within a user-defined neighbourhood diameter or projection scale, d and (3) measurement of the distance between the average position of each point cloud along the normal vector (Lague et al., 2013). A trial-and-error approach was first used to estimate the normal diameter, D. Surface normals must be estimated at a scale that is small enough to capture medium-large scale changes in surface orientation, such as small indentations or coves, or changes in cliff aspect, but large enough to avoid fluctuation of the resulting normals due to small scale changes in roughness, σ, such as overhangs on the rockface. On this basis, a normal diameter, D, of 2m was selected as a compromise. Projection scales where d <1m sample too few points from each point cloud and therefore exaggerate change, while larger projection scales (d >3m) cause data smearing and averages out many of the largest changes recorded (Williams et al., 2018). This effect is most pronounced on blocky and vegetated surfaces. The diameter of the projection cylinder, d, was therefore specified as 1m to ensure that the number of points sampled in each cloud was >30, while minimizing the effect of averaging.
Cloud-to-cloud distances and the associated statistics were projected onto both the pre-and post-event point clouds (Figure 2a), and both point clouds were then filtered to remove areas of deposition and insignificant change, as defined by a spatially variable confidence interval ( Figure 2b). This is also known as the Level of Detection (LoD) threshold for a 95% confidence interval: where reg is the user-defined registration error, which is substituted for the root mean square alignment error for point clouds n 1 and n 2 (Lague et al., 2013). For consistency within and between monitoring intervals, the largest possible registration error between two point clouds was used (~0.10m). It is assumed that this error is isotropic and spatially uniform across the dataset, given the manner in which the data was collected and the generally consistent point distribution and spacing.
Having isolated the areas of erosion, the pre-and post-event rockfall points were then merged into a single cloud for each block, for each monitoring period ( Figure 2c).

Rockfall identification and characterization
Following pointwise change detection, the data belonging to each individual rockfall event were grouped using the clustering algorithm DBSCAN (Density-Based Spatial Clustering of Applications with Noise), developed by Ester et al. (1996). Contiguous points of change were assumed to belong to a single event, although the influence of rockfall scar coalescence has been recognized at intervals beneath the annual timescales used here (Barlow et al., 2012;Williams et al., 2019). DBSCAN is the most commonly used single-scan clustering technique and defines clusters based on the local density of points. The algorithm requires the user to define the minimum number of points (MinPts) within a maximum distance (ε) from each randomly chosen point (p) in the dataset, but does not rely on a user-defined number of clusters. DBSCAN defines a neighbourhood of points, N ε , which falls within the circle of radius ε around a point, p. MinPts is defined as the minimum number of neighbours of point p to consider p as a core object. If N ε contains more than MinPts, the algorithm creates a new cluster with p as the core point, and iteratively collects directly density-reachable points from p. The process terminates when no new points can be added to any cluster. If each point is mapped to the distance from its kth nearest neighbour, the threshold point p in the k-distance graph can be used to define ε = kth-dist(p), where MinPts = k. All points with a higher k-dist value are considered to be noise (Ester et al., 1996). For databases where each point only occurs once, Sander et al. (1998) propose that MinPts is equal to 2 × the number of dimensions.
After merging the filtered pre-and post-event rockfall inventories, each block of points was run sequentially using a parallel DBSCAN algorithm (PDSDBSCAN-S) developed by Patwary et al. (2012) for shared-memory computing. MinPts was set to 6 (2 × the number of dimensions) and ε was determined by plotting the sixth-distance graph for the rockfalls in each block and averaging the distance at the point of inflection, p, for all blocks (ε = 0.45m). The results for each block were manually verified, and any noise objects were filtered out of the dataset. For each monitoring period, this gave each rockfall point a unique class identifier ( Figure 2d).
The use of triangulated surfaces for surface reconstruction and the volumetric characterization of objects is well-established, and some attempts to quantify rockfall volume and shape have been made using the alpha shapes method at a site-specific scale (Guerin et al., 2014;Carrea et al., 2015;van Veen et al., 2017). However, these are often subject to considerable post-processing in order to ensure a watertight mesh and fail to provide error estimates for the calculated volumes. Here, the Power Crust algorithm was used on an event-by-event basis to construct a triangular surface mesh for each rockfall (Amenta et al., 2001; see Figure S2, Supporting Information for a detailed description). The power crust is a watertight boundary of a 3D solid described by the approximate medial axis transform (or 'power shape'), and has the advantage that it eliminates the need for the polygonization, hole-filling and/or manifold extraction post-processing steps required in most surface reconstruction algorithms, including the alpha shapes method (Berger et al., 2014;Lim and Haron, 2014). The algorithm uses each set of rockfall points and a tolerance, defined between 0 and 1, as inputs. The tolerance is used to define the boundary of the power crust; in the majority of cases, a higher tolerance yields a more robust fit, although this varies from mesh to mesh. To ensure the most robust fit in every case, each rockfall was meshed at nine different tolerances and the mesh closest to the average volume of the nine resulting meshes was chosen (Figure 2e). Volume error estimation has not traditionally been accounted for in rockfall magnitude-frequency distributions; uniquely, here we have determined the lower and upper bounds of the calculated volumes for each rockfall by using the smallest and largest possible meshes, respectively. A detailed description of this process is given in Figure S3 (Supporting Information). It should be noted that small rockfalls will, on average, have a lower error in the volume estimate, as there are limited meshing configurations for small numbers of points. If there is only one meshing configuration for a set of points, the error is reported as ± 0.00m 3 .
The volume of each rockfall mesh was calculated using the Divergence Theorem. All rigid bodies, and therefore their parameters, can be expressed in terms of 3D moments (Semechko, 2014). Closed-form expressions for the 3D moments of objects represented by triangular surface meshes can be derived and used to calculate volume, V, which is equal to the zeroth moment: and the position of the centre of mass, which is equal to the value of the first moments of x, y and z, divided by volume: The surface area, width, depth and height of each rockfall mesh was also calculated. This gave four rockfall inventories (2014-2015, 2015-2016, 2016-2017 and 2014-2017) captured across the 24km of the North Yorkshire coastline between Skinningrove and Sandsend. The properties recorded in the rockfall inventories are listed in Table S2 (Supporting Information).

Negative power law estimation
Rockfall magnitude-frequency distributions exhibit a negative power law scaling and can be modelled as: where f (V R ) is the frequency density of a rockfall of magnitude V R , and s and β are empirical constants (for example, Dussauge-Peisser et al., 2002;Guzzetti et al., 2004;Brunetti et al., 2009;Barlow et al., 2012;Bennett et al., 2012; see Table S3, Supporting Information for relationships derived from previous terrestrial monitoring of rockfalls). These provide an overall indicator of the level of activity and relative size distribution in an inventory, respectively. For each inventory, rockfall magnitude-frequency was plotted on logarithmic axes using logarithmically spaced bins. Frequency densities were calculated for rockfalls of differing magnitudes using the formula provided by Malamud et al. (2004): where δN R is the number of rockfalls with volumes that fall within the range of δV R , and δV R is the associated bin width. 2777 CHARACTERISTICS OF REGIONAL SCALE ROCKFALL INVENTORIES Parameter estimation was undertaken using least squares regression (LSR) on the logarithmically transformed data, as LSR is considered a robust estimator provided that the input data is complete through the mid-range (Goldstein et al., 2004). This is consistent with previous research undertaken along the coastline and ensures full comparability with other findings. Due to the high rate of rockfall along the cliffs monitored here (Lim et al., 2005(Lim et al., , 2010Rosser et al., 2005aRosser et al., , 2005bRosser et al., , 2007Rosser et al., , 2013; Barlow et al., 2012;de Vilder et al., 2017), we also assessed the effect of superimposition and coalescence of rockfall scars on the form of the magnitude-frequency distribution by comparing the scaling coefficients derived from annual inventories to those of a change detection between the first and last survey only. In all cases, the data only follow a negative power law for events greater than~1 × 10 -3 m 3 , i.e. with no rollover for the smallest events . This can be attributed to censoring by under-sampling and other biases, such as the relatively high threshold that was set for the minimum detectable change (0.10m) during data processing as well as differences in the way that direct cloud-to-cloud comparison methods such as M3C2 identify and treat insignificant change. All parameters quoted hereafter are for power laws fitted to the uncensored data only.

Rockfall magnitude, frequency and erosion
Over 58000 rockfalls were captured along 20.5km of cliffs between August 2014 and March 2017 (32 months), in three separate, approximately equal intervals. The area monitored constitutes~805739m 2 of cliff face, with an average cliff height of~40m (Table 2). Rockfalls ranged in volume from <0.0001 ± 0.00m 3 to~15498.05 ± 552.36m 3 , with a mean rockfall volume of 2.15 ± 0.24m 3 (Table 2). Along the coastline, rock yield totalled 124843m 3 from August 2014 to March 2017, representing an average rate of recession of 0.06myr -1 when derived by spatially averaging the observed rockfall volume over the monitored area. This is the same order of magnitude as erosion rates derived from previous terrestrial monitoring of rockfalls between Boulby and Staithes (Table S4, Supporting Information). Across the inventories, the average meshing error for the total eroded volume is approximately ±10.73%, which is relatively low despite the conservative approach used to calculate error margins (Abellán et al., 2014).
The magnitude-frequency distributions for rockfalls captured over the three monitoring periods were modelled using negative power law scaling relationships (Figure 3a), where the magnitude of the exponent β ranges from 1.54 (2016-2017) to 1.69 (2014)(2015). These values fall inside the 1.00-2.00 range commonly found for non-cumulative plots of rockfall magnitude-frequency density (Table S3). As demonstrated by Barlow et al. (2012), the superimposition and coalescence of rockfalls has the effect of lowering the power law coefficient of the single change detection over a period (2014-2017; three years) relative to that derived from more frequent sampling over the same period (annually; Figure 3b). This explains the decrease in the overall number of rockfalls observed (from 58 032 to 25969), and a corresponding increase in individual rockfall volumes. Rockfall inventories derived from the annual change detections were therefore used for further analysis and for deriving all of the power law scaling coefficients used here. We recognize that our distributions (captured annually), as opposed to previous local scale monitoring (captured monthly), could be prone to under-sampling and other biases, such as the threshold that was set for the minimum detectable change (0.10m) during data processing, as well as differences in the way that cloud-to-cloud comparison methods identify and treat insignificant change. Previous research has demonstrated the relative stability and reliability of power laws fit to rockfall volume probability distributions at monitoring intervals >12hours (Williams et al., 2019). We therefore consider that any bias in the data collected most likely relates to rockfalls that were not captured (e.g. as a result of the level of detection or resolution).
The mean coastal erosion rate varied between years, increasing by an order of magnitude from 0.02myr -1 (2014-2015) to 0.10myr -1 (2016-2017; Table 2). This change is partly driven by an increase in the density or rate of rockfalls year-on-year, although the majority of this increase can be accounted for by the occurrence of eight large (>1000m 3 ), full-scale cliff collapses that occurred along the coast, comprising half of the total volumetric flux observed during 2016-2017. The mean rockfall volume more than trebled in this time, from 0.99 ± 0.04m 3 (2014-2015) to 3.31 ± 0.28m 3 (2016-2017), while the median (~0.01m 3 ) and mode (~0.2 × 10 -2 m 3 ) remained relatively constant. The dominance of larger events is also reflected in the magnitude-frequency distribution; the exponent, β, of which decreases in magnitude through time (Figure 3a). Year-on-year, events of all magnitudes increased in frequency, marking an overall increase in the rate of rockfall activity along the coastline. Rock yield averaged 2.32m 3 per linear metre of coastline per year, doubling that previously recorded along seven sites between Staithes and Boulby (~710m of coastline) over seven years (Rosser et al., 2013). This totalled 124843m 3 (equivalent to the volume of a~50m cube) during the monitoring period, despite a low mean erosion rate of 0.06myr -1 ( Table 2).
Rates of cliff retreat along the North Yorkshire coast are highly variable within a single year (Figure 4a-c), with retreat rates ranging from as much as 1.12 × 10 -5 to 1.63myr -1 (2016)(2017). Erosion rates at a number of cliff sections also sharply increased over the monitoring period, most notably at Port Mulgrave, which retreated at an average rate of <0.002 myr -1 during 2014-2015 but, due to a number of large rockfalls and landslips (31943.47m 3 in total), this rate increased to 0.04 myr -1 (2015-2016) and then 0.28myr -1 (2016-2017). Counter to this variability, in many places there are consistent patterns in the spatial variation of erosion derived from binned (100m lengths of coast) rockfall over time (Figure 4). A gradual gradient in erosion rates is common across contiguous 100m sections of the coast, indicating that the nature of cliff erosion and rockfall captured here represents a response to controls that vary systematically over scales >100m. This is in contrast to a stochastic distribution, which would suggest that erosion rates vary independently, and hence randomly, between immediately adjacent 100m cliff sections.
The~3.5km stretch of cliffs between Boulby and Cowbar Nab retreated at an average rate of 0.02myr -1 during 2014-2015, increasing to a rate of 0.07myr -1 during 2016-2017. The highest rates of retreat along the coastline are observed here, with local erosion rates regularly exceeding 0.05myr -1 and even reaching 1.47myr -1 in the event of a full-scale cliff collapse (15498.05 ± 552.36m 3 , 2016-2017). The highest frequencies of both small (≤0.1m 3 ) and large (>0.1m 3 ) rockfalls also occurred along this stretch of cliffs (Figure 4), with over 11000 of the recorded 58032 rockfalls located along this section. This is partly attributable to the fact that the cliffs between Staithes and Boulby are the highest along the coastline, reaching up to 150m towards Boulby, and so have a greater surface area available to release rockfalls. North facing areas such as Boulby show the highest rates of retreat, most likely 2778 J. BENJAMIN ET AL.
due to their exposure to easterly and northerly North Sea storm waves in comparison to sheltered embayments such as Runswick Bay, which eroded at an average rate of 0.005m yr -1 between 2014 and 2017. By contrast, erosion rates remained low at Kettleness throughout the monitoring period, where the cliffs eroded at a mean rate of~9.70 × 10 -3 myr -1 .

Spatial variation in scaling parameters
Although the length of time over which rockfall frequency estimates are made is known to exert a profound influence on β (Barlow et al., 2012;Williams et al., 2019), the length scale, L, over which power laws can be applied remains poorly constrained. It is also apparent that many previously reported observations do not indicate whether power laws derived on a site-specific basis (e.g. L <10 2 m) can adequately explain behaviour over larger spatial scales (e.g. L >10 3 m). Given the emerging need and capability to address this scale of investigation (Kennedy et al., 2017), we considered the effect of using a variable length scale, or here the area of coastal cliff monitored, to calculate β. In order to ensure that the monitored frequency density of large events was not influenced by a sampling bias associated with short-term monitoring of a stochastic process (an infrequent large event captured, or missed, by chance within a short time window), we considered the size of the longest rockfall axis across the three inventories ( Figure 3a). We therefore set L to 100m, which also approximates the scale of site-specific monitoring (for example, Lim et al., 2010), increasing in 100m increments up to a maximum window length of 24 km, which is equal to the length of coastline monitored (L max ). For each value of L, we estimated β using a sliding window of length L, repeating for 240 iterations, I. The sliding distance, S, of the window is inversely proportional to its length, such that: The relationship between the β value and the extent or horizontal length scale of monitoring, L, was modelled using a  2779 CHARACTERISTICS OF REGIONAL SCALE ROCKFALL INVENTORIES two-term power series model (Figure 5a, b). Over each monitoring period, an inflection occurs around L ≈ 2.5km (Figure 5a), implying that a magnitude-frequency distribution that is physically meaningful for modelling longer-term or wider-scale cliff erosion is captured only at length scales that exceed this distance (equivalent to~1 × 10 5 m 2 assuming an average cliff height of 40m), when surveyed at approximately annual intervals. Where rates of rockfall are low, or where the largest possible event is larger than that observed here, this length scale is likely to increase further.
Values of β derived from previous terrestrial monitoring of rockfalls at a number of sites along the North Yorkshire coast range between~0.80 and 2.40 (Figure 5c). These sites were monitored at total length scales ≤600m and at approximately equal monitoring intervals (~30days). Rates of erosion derived from previous monitoring vary between 0.004myr -1 and 0.128myr -1 , capturing well the mean rate of retreat along the coastline but poorly representing the extremes of the distribution (Figure 6a, b). Given that our data indicate that values of β converge when the spatial extent of monitoring is increased beyond L = 2.5km, we also considered the influence of whether or not monitoring in a continuous section that exceeds this length scale, or whether an extent of multiple segments that total the same length derives the same statistical results. Neighbourhood differences in the erosion rates shown in Figure 5 indicate that, for years with less variability in erosion rates (2014)(2015)(2015)(2016), there is a relatively small difference in erosion rate between each 100m bin and its neighbours (Figure 6c).

Rockfall shape
We explored the self-similarity of rockfalls by considering their 3D shape. This provides evidence for changes in detachment mechanism with rockfall scale, and hence the applicability of a single power law. Rockfalls were classified in relation to the aspect ratios of the principal length axes (a, b, and c) to give one of three end-member shapes, either equant (blocky), platy (slabs) or elongate (rods), which themselves are divided into 10 sub-categories (Sneed and Folk, 1958). We observe a clear scaling effect across the 58032 rockfalls captured along the coast between 2014 and 2017 (Figure 7), which, importantly, is a trend that is consistent year-on-year irrespective of overall rockfall or erosion rate ( Figure S4, Supporting Information). Rockfalls occurring at the tails of the distribution (<10 -3 m 3 and >10 1 m 3 ) present approximately equal proportions of equant, platy and elongate shapes, with a tendency towards very platy and elongate forms, particularly where V R >10 2 m 3 . However, there is a clear peak in the distribution for rockfalls of volume 0.002m 3 to 0.064m 3 (n = 30766), the majority of which (>70%) are blocky in shape. This equates to a cube of dimensions between~0.13m and~0.40m, which is of a similar order to discontinuity-defined structural control on block size and shape at this scale (see de Vilder et al., 2017). This observation also suggests that the larger rockfalls are not just larger versions of small events, and therefore that different detachment mechanisms may operate at different scales of rockfall.

Discussion
The ability to precisely quantify and subsequently understand rockfall behaviour is critical for successfully modelling the present and future dynamics of actively failing rock slopes. The reliability and efficiency of rockfall hazard protection measures depends on the outcome of these modelling practices . However, quantifying rockfall activity has proven problematic, particularly at the regional scale, with a range of approaches currently used to measure the retreat, area or volume of changes in rock slopes (Abellán et al., 2014). These approaches often resort to gridding the data into erosional cells in order to reduce processing time and complexity (see for example, Leyland et al., 2017). We have presented a novel workflow, from data collection through to analysis, that is applicable to object-based change detection on extensive near-vertical surfaces, such as cliffs, riverbanks, canyon walls and transport corridors. We used this workflow to detect and characterize regional-scale rockfall activity from airborne LiDAR data, giving an inventory of >58000 rockfalls along 20.5km of coastal cliffs over 32months. We discuss the insights that this dataset generates, and the implications of our findings for monitoring and modelling rockfall activity.

Implications for monitoring rockfall activity
Our data show that a full 3D treatment of new forms of high-resolution airborne LiDAR provides a robust means to monitor rockfall activity and the resulting cliff retreat continuously (in space), in 3D, and over large spatial scales (>10 3 m). When mounted on a helicopter, continuous swaths from this type of airborne LiDAR can be used to collect data along narrow corridors even with steep, near-vertical slopes, presenting a considerable advantage over mobile terrestrial LiDAR when scanning areas that are often limited by range, view direction, and can feature extensive occlusion (Lato et al., 2009;Dunham et al., 2017). Point cloud data obtained from LiDAR gives precision on slope angle, aspect, and, on bare rock faces, lineaments and other structural features, while multi-return and intensity data can be used during segmentation and classification routines (Axelsson, 1999;Sithole and Vosselman, 2004;Vosselman et al., 2005) as a means for biomass filtering and assessment of slope stability (see Chen, 2013Chen, , 2015, for a The workflow developed and used to process airborne LiDAR data here is semi-automatic and, uniquely, provides a means of analysing regional-scale variations in rockfall activity along rock slopes with a non-linear planform geometry. The approach outlined here gives a 3D watertight mesh, centre of gravity, principal axes, volume, and volumetric uncertainty for each rockfall. Information such as rockfall shape, volume, source area location and cliff surface geometry are required by 3D models used for rockfall susceptibility mapping (Leine et al., 2014). These models are advantageous over their two-dimensional (2D) equivalents as they allow for calculating both the onset (or triggering) probability as well as modelling trajectories (Crosta and Agliardi, 2003;Lambert et al., 2013;Macciotta et al., 2015;Fanos et al., 2016), meaning that their outputs can be used for optimizing maintenance and remediation works, as well as for designing countermeasures. Although the effect of rockfall size on runout length is well-known (for example, Statham, 1976;Evans and Hungr, 1993;Okura et al., 2000;Copons et al., 2009), the importance of rockfall shape for determining runout distance and the depositional patterns of debris has only recently been quantified (for example, Sellmeier, 2015). The 3D nature of the rockfall inventory derived here would permit real observed rockfall detachments to be introduced into rockfall modelling on a regional scale, allowing for a true representation of the natural variability of rock shapes and sizes with relation to their geological provenance. This gives the potential to link cliff geology and structure to rockfall shape and geometry on a regional scale.
Our results show that the probability distribution of rockfall shape as a function of volume is broadly predictable, and, crucially, that rockfall shape is not scale invariant. While a much greater proportion of very small (<10 -3 m 3 ) and very large (>10 1 m 3 ) rockfalls tend towards very platy and elongate forms, the majority of rockfalls (>70%) of volume 0.002m 3 to 0.064 m 3 (n = 30766) are blocky or 'equant' in shape. This pattern is consistent and occurs year-on-year, marking a transition from rockfall as a structurally-defined process to rockfall as either small scale consequences of incremental weathering (V R < 10 -3 m 3 ), or fracturing-related mass-movements that break through rock bridges to generate larger, predominantly cliff face-parallel rockfall (V R >10 2 m 3 ). The relationship between rockfall shape and volume that we observe along the North Yorkshire coast holds important consequences for modelling future cliff erosion, both along cliffed shorelines and in non-coastal settings, as such relationships could be used to place approximate bounds on the dimensions of modelled rockfall, and the step-back events that result (for example, Dong and Guzzetti, 2005;Young et al., 2011), using rockfall magnitude-frequency relationships.
A key finding of the analysis presented here is that localized estimates of rockfall activity captured over small extents (<10 3 m 2 ) do not generate stable magnitude-frequency distributions, and so cannot be upscaled for the purpose of modelling wider-scale or longer-term cliff evolution. We show that monitoring at length scales <2.5km (equivalent to a minimum area of~1 × 10 5 m 2 assuming an average cliff height of 40m) has a significant effect on the frequency estimates of the largest events, potentially giving rise to considerably higher (where, by chance, a large event is captured) or lower (where no large event occurs) frequencies than is actually the case. Here, the surface area of the largest recorded event is~7.5 × 10 3 m 2 , equating to 7.5% of this 'minimum' area. Along the North Yorkshire coast, this figure represents~12% of the total cliff length or area monitored. This extent, which is a function of the probability of being able to capture a statistically valid sample of the largest possible events, is likely to vary between settings based on differences in weathering and other environmental conditions (for example, precipitation, temperature, and the frequency of triggering events), triggering mechanisms, and lithological characteristics, each of which are thought to influence the power law scaling of rockfalls (Barlow et al., 2012).
The relatively small difference in erosion rate between each 100m bin of coastline and its neighbours (Figure 6c) indicates that there is likely to be more structure (or less variation) in the erosional signal if contiguous cliff sections are monitored. Our data therefore show that, in order to increase the likelihood of capturing a stable magnitude-frequency relationship and therefore a complete distribution of rockfall activity, monitoring of the cliffs under examination here should be undertaken at multiple sites totalling at least 2.5km in length. This is at odds with previous site-specific monitoring undertaken here and in other settings. We argue that this would act to overcome local (~10 2 m) structure in the data and allows an assessment of the general behaviour of the wider coastline. The length scale that we have found for the North Yorkshire coast, 2.5km, is likely to vary between settings based on, for example, variability in rock strength and structure, and ideally should be constrained elsewhere. Assuming that the patterns of erosion observed form part of a longer-term cycle of cliff failure and profile-form adjustment that is not fully captured by the relatively short duration of monitoring undertaken here (twoyears and sevenmonths), then neighbouring bins are more likely to be at a similar stage of this process than more distal sites. Previous research has demonstrated that the propagation of instability and failure along these cliffs operates at 10 1 year timescales. For example, Rosser et al. (2013) estimated an average resurfacing time of 28.1years along~710m of cliffs at sites A-G ( Figure 6) based on extrapolating the spatial footprints of failures derived from monthly monitoring over a seven-year period. The assumption of similarities between neighbouring bins can be justified given the longer wavelength variations of rock mass structure along the coastline (~10 3 m), which are likely to moderate the erosional effects of shorter wavelength variations in wave loading (<10 2 m). If this cyclical nature of cliff erosion and retreat does exist, and it is characterized by time-dependent failure processes (for example, incremental oversteepening leading to large-scale failure), then a stable magnitude-frequency relationship can only be observed over a longer total duration of monitoring, or over a more widely-distributed area of monitoring.
Moreover, sampling in a distributed manner is more likely to capture the erosional response of the cliffs to a wider variety of controls and, by inference, at different stages of the longer-term failure cycle. However, it is important to note that this space-for-time substitution might yet break down and under-sampling could remain if environmental extremes promoting rockfall occurrence are not adequately represented in a given monitoring period.
Implications for modelling rockfall activity and coastal cliff erosion Given the increasing tendency towards collecting large inventories on potentially hazardous geomorphic processes (Korup et al., 2012), our findings attest to the importance of collecting a dataset that is both temporally and spatially complete. Although large pools of data were previously thought to enable a statistically robust analysis of magnitudes, frequencies and the formulation of exceedance probabilities for hazard appraisals (Korup et al., 2012), we argue that this is only the case if the quantitative input for these is complete over a large enough spatial scale relative to the scale of the events experienced. This is pertinent for research using terrestrial LiDAR to monitor rockfall activity, as such approaches typically operate at short length scales (see Abellán et al., 2014, for a review). More widely, our analysis has considerable implications for modelling the evolution of non-coastal slopes, where magnitude-frequency scaling is often used to inform hazard assessment and mitigation (Abellán et al., 2011;Dewez et al., 2013), and where accurately defining rockfall recurrence intervals is essential (for example, Budetta and Nappi, 2013;Budetta et al., 2016;Moos et al., 2017).
The literature on coastal cliff behaviour and evolution often categorizes cliffed coastlines into perceivably stable 'hard' rock cliffs and actively retreating 'soft' rock cliffs, with the implication of this distinction being that soft rock coasts are considered more vulnerable to instability and rapid change (Allison, 1989;Sherman and Gares, 2002). Along the North Yorkshire coast, erosion rates ranged from as much as 1.12 × 10 -5 -1.63myr -1 over the monitoring period of this study. This reflects wider trends across the British Isles, where rates of landward retreat range from <0.001myr -1 in what are apparently the most resistant rocks, to >10myr -1 where cliffs are composed of soft, conformable glacial tills (Brooks and Spencer, 2012). The susceptibility of soft rock coastlines to erosion and retreat has meant that they are often prioritized in coastal management schemes (Lee and Clark, 2002), leaving the binary distinction implying some level of uniformity of behaviour or erosion within each sub-category. The results of our analysis demonstrate that local (10 2 m) erosion rates along stretches of perceivably 'hard' rock coastline can reach approximately the same order of magnitude as the highest rates of erosion occurring on soft rock coastlines, even over the annual timescales considered here. Along the North Yorkshire coast, almost half (~53 722m 3 ) of the total (~124843m 3 ) volume eroded by rockfalls between August 2014 and March 2017 accrued in only 12 large (>1000m 3 ) cliff collapses, eight of which occurred between April 2016 and March 2017, resulting in an instantaneous step-back of the coastline by up to 6m in places. Assuming that the longest axis of these events is cliff-parallel (representing rockfall width) and the shortest axis is cliff-normal (representing rockfall depth), then these events alone caused an average step-back of~1.92m over 972m of the coastline, equating to approximately 4.8% of the total cliff length monitored over a period of twoyears and sevenmonths. 2783 CHARACTERISTICS OF REGIONAL SCALE ROCKFALL INVENTORIES Present models of cliff retreat fail to capture the timing and scale of episodic events, and, at present, little is known of how the long-term rates of erosion derived by these models arise from the accumulation of individual, instantaneous events. Rising global sea-levels in conjunction with projected changes in winds, tides, precipitation, storm events, and wave climate are expected to accelerate coastal cliff retreat and threaten coastal populations in many areas (Sunamura, 1988;Bray and Hooke, 1997;Dickson et al., 2007;Nicholls et al., 2007;Trenhaile, 2010Trenhaile, , 2014, resulting in a pressing need to understand and model the erosional response of hard rock coastlines like the North Yorkshire coast to these processes (Trenhaile, 2011). We see the approach presented here as being fundamental to this, as it shifts focus to a regional scale and allows variability in drivers to be considered at scales previously not possible.

Conclusions
The results presented here provide unique insight into rockfall dynamics and how they vary on a regional scale. Using high-resolution, multi-temporal airborne LiDAR data we have explored regional-scale variations in rockfall magnitude and frequency, to what extent these relations are sensitive to the spatial scale of monitoring, and the implications of our findings for modelling cliff evolution, both in coastal and non-coastal settings. Our main conclusions are summarized as follows: 1. High-resolution airborne oblique LiDAR provides a robust means to monitor rockfall activity continuously over large spatial scales (>10 3 m). The workflow presented here is semi-automatic, providing a 3D mesh, principal axes and volumetric uncertainty for each rockfall, and, uniquely, it takes into account spatial variations in both cliff profileand plan-form. Given that rockfall shapes, volumes, source area locations and cliff surface geometry are known to influence rockfall trajectories, our research provides a point of reference for future investigations concerned with regional scale rockfall hazard modelling. 2. Rockfall magnitude-frequency relationships are highly sensitive to the spatial scale of monitoring, such that monitoring at length scales <2.5km significantly increases the frequency estimates of the largest events observed on this coastline. This has profound implications for research using approaches that measure rockfall scars to generate inventories, such as structure-from-motion or terrestrial LiDAR, both in coastal and non-coastal environments, as any scaling relationships derived may be subject to significant bias as a function of spatial monitoring extent. 3. Rockfall shape is not scale invariant, and the relationship between rockfall volume and block shape is persistent year-on-year, giving the potential to predict rockfall future shape distributions using magnitude-frequency distributions. Variations in mean rockfall shape with volume also imply a systemic shift in the underlying mechanisms of detachment with scale, questioning the validity of applying a single probabilistic model to the full range of observed rockfall volumes at these sites. 4. Local (10 2 m) rates of coastal cliff erosion can vary over six orders of magnitude along a 20.5km stretch of hard rock cliffs. This erosion is marked by the widespread occurrence of episodic step-back events, 12 of which accumulated an average step-back of~1.92m over nearly 5% of the cliff length monitored. Our findings dispel the concept that hard rock coastlines are relatively stable and highlights the importance of understanding and modelling the erosional response of hard rock coastlines under a changing climate.
Our findings have demonstrated the importance of sea cliff retreat as an episodic process, where sudden, large (>1000 m 3 ) rockfall events punctuate periods of relative stability. How the erosional work done by episodic, large-scale events accumulates into a long-term rate of erosion remains to be seen, with present models of cliff retreat failing to capture the timing, scale or drivers of these events. Uniquely, the methodologies and data presented here mark a step-change in our ability to understand the competing effects of different processes in determining the magnitude and frequency of rockfall activity, meaning that we are therefore better placed to investigate relationships between process and form/erosion at critical (regional) scales.

Supporting Information
Additional supporting information may be found online in the Supporting Information section at the end of the article. Figure S1. Map of the North Yorkshire coast, showing (a) ground control points used in the four airborne LiDAR surveys undertaken. Red points denote sites that were repeated in all four surveys, while blue points denote those that were used in the first survey only, (b) orthophoto tiles, which cover the full extent of the surveyed area, and (c) the 38 blocks that the data were divided into for processing. Map produced using shapefiles from Ordnance Survey © Crown Copyright and Database Right 2017. Ordnance Survey (Digimap Licence). Figure S2. Principles of the Power Crust algorithm, depicted on a two-dimensional object for simplicity. For a given group of points, S, the Power Crust algorithm is able to extract a simplified skeletal shape, or the medial axis, which is then used to produce a surface representation of the points. The medial axis transform (MAT) represents a solid by the set of maximal balls completely contained in its interior (a). The MAT is approximated by a subset of Voronoi vertices of S, called poles, which lie near the medial axis (b). The balls surrounding the poles are known as polar balls (c), the radius of which determines the weighting of each pole. An inverse transform is approximated by using a power diagram of the weighted poles. This acts like a weighted Voronoi diagram by dividing space into polyhedral cells (d). These are then divided into interior and exterior faces, where the boundary of separation of these subsets forms the output surface, or power crust (e). The power crust is therefore a watertight boundary of the three-dimensional solid described by the approximate MAT (or power shape), and eliminates the need for the polygonalization, hole-filling, or manifold extraction post-processing steps required in other surface reconstruction algorithms. Diagram adapted from Amenta et al. (2001). Figure S3. Three three-dimensional triangular surface meshes generated using the Power Crust algorithm. The minimum, average and maximum possible mesh sizes for the given rockfall are shown. The algorithm uses each set of rockfall points and a tolerance, defined between 0 and 1, as inputs. The tolerance is used to determine the inner and outer poles of the power diagram, which defines the boundary of the power crust ( Figure S2d). In the majority of cases, a higher tolerance yields a more robust fit, although this varies from mesh to mesh. To ensure the most robust fit in every case, each rockfall was meshed at nine different tolerances and the mesh closest to the average volume of the nine resulting meshes was chosen. The lower and upper bounds of the calculated rockfall volumes were then determined for each rockfall by using the smallest and largest possible meshes, respectively. Figure S4. Rockfall shape monitored along the North Yorkshire coast, UK, from (a) 2014-2015, (b) 2015-2016, and (c) 2016-2017. The results are plotted as a stacked bar graph, with colours corresponding to the 10 sub-categories shown in Figure  7b, defined by Sneed and Folk (1958). Histograms show the frequency of rockfalls in each volume bin. Table S1 Summary statistics of the ground control data for each of the four airborne LiDAR surveys undertaken. Table S2 Summary of the properties recorded for each rockfall in the inventory. Table S3 Absolute β/ρ values derived from previous terrestrial monitoring of rockfalls, sorted by β. Values in bold are those reported in their respective studies. Table S4 Erosion rates derived from previous terrestrial monitoring of rockfalls between Boulby and Staithes. Data are used in Figure 6(a). Table S5 Absolute β values derived from previous terrestrial monitoring of rockfalls along the North Yorkshire coast. Data are used in Figure 5(c). 2787 CHARACTERISTICS OF REGIONAL SCALE ROCKFALL INVENTORIES