1. Introduction
Rockfall is a persistent geological hazard that occurs all over the world. It is especially prevalent in mountainous and hilly terrain where falling rocks can harm infrastructure, housing, roadways, and people [
1,
2]. Rockfall is a form of fast-moving landslide, and the modified Varnes classification describes rockfall movement as the “detachment, fall, rolling, and bouncing of rock fragments” [
3]. Rockfall can come in the form of singular fragments or clusters of fragments, and these fragments may break apart during the fall. Large and small rockfall volumes can be dangerous when they fall from a significant height because the higher the rock falls from, the greater amount of potential energy it contains [
4]. Falling rock fragments are typically characterized separately from large masses of moving rock, such as rock avalanches, slides, and slumps, because of the dissimilar nature of their movements and interactions with substrates [
3].
An element of most approaches to evaluate the risks and avoid the consequences of rockfall hazards is to map their extents and possible areas of impact and then assess the hazard intensity in those areas [
2]. Commonly, in-person field surveys are conducted, either before or after rockfall events occur, to gather data about the nature of the hazards [
5]. However, visiting rockfall locations is not always possible, so an alternative would be to employ remote sensing methods to collect inventory data to use in mapping and assessing hazards. A rockfall hazard inventory involves identifying and recording fallen rock blocks and their known characteristics, as well as the factors that control a slope’s instability [
6]. Once an inventory is established, it can be used in the characterization of hazard intensity to produce a rockfall hazard rating map.
Data from remote methods can also be a helpful supplement to field survey data. Due to advances in remote sensing technologies, there is an abundance of data available for all parts of the world in the forms of multispectral optical imagery, digital elevation models (DEMs), synthetic aperture radar (SAR) images, and interferometric synthetic aperture radar (INSAR), to name a few, that can be used to understand and predict geological hazards [
7]. Additional remote sensing techniques that have been used extensively to collect and analyze terrain and rockfall data include light detection and ranging (LiDAR) surveys and photogrammetry methods to create 3D projections of slopes [
1]. Analyzing aerial photographs has been, and still is, the most common method used to identify and map landslide-type hazards according to [
5]. As explained in [
7], “aerial photographs capture a range of characteristics that are useful for image recognition and classification, which include the surface shape of the topography, along with the color, tone, patchiness, and texture.” Google Earth takes aerial photo mapping a step further and provides a publicly available digital platform for viewing high-resolution aerial images, and it includes elevation data for creating estimated profiles of topography. Free tools such as Google Earth allow people to learn about the hazards in their area, create their own inventories, plan mitigation, and determine locations that are safe for future development.
The motivation for this study started with the observation that there are limited characterizations and inventories of rockfall hazards for the Arequipa region of Peru, which makes the management of community development in rockfall-prone areas difficult. Specifically, the focus is on creating a method for rockfall hazard characterization, partly because rockfall is a widespread hazard across the region and partly because there is limited development of rockfall hazard rating systems (RHRSs) that use only remote methods. While many geologic and topographic factors have been used to characterize rockfall hazards, this study focuses on factors that are available remotely and analyzes how effective they would be in predicting rockfall hazards. Most RHRSs consider hazard on the scale of individual roadcuts because of the field methods employed, and others use a regional or country-wide scale, which is not as useful for individual communities. In contrast, the aim of this research is to provide a tool that can be used to understand rockfall hazards at the community scale, between these two scales, to aid researchers, local governments, communities, and non-profit organizations in the Arequipa region.
The terms “hazard” and “risk” each have distinct meanings and should be defined to outline the expectations of the final products of this study. “Hazard” refers to “a condition with the potential for causing an undesirable consequence” [
8]. In terms of rockfall, the condition would be falling rock, and the consequence would depend on parameters such as the location, volume, trajectory, and velocity of the rock [
9]. “Risk” is defined in [
8] as “a measure of the probability and severity of an adverse effect to health, property or the environment.” So, rockfall that occurs in a hard-to-access site in an uninhabited area and has no impact on life or property could have a very high hazard rating, but its risk rating would be low because the probability that the rockfall would kill or injure someone or destroy a structure is almost nonexistent. On the other hand, a small rockfall slope may have a low hazard rating due to its size, but the risk to a house located directly below the slope would be high because now there is a high probability that any given rockfall would impact the structure. This study evaluates the hazard rating of rockfall and not the risk rating. All the sites are centered around small communities, so having information on the hazard rating is useful for them to determine if any of their infrastructure or homes are potentially susceptible to rockfall impact. While this research incorporates scientific tools proven elsewhere, its novelty lies in both the scale of the work and the linking of rockfall inventories calibrated by runout modeling to improve the accuracy of hazard zonation.
2. Existing Rockfall Hazard Rating Systems
It should be emphasized that the areas over which rockfall hazard is predicted in this study span small towns in the Arequipa region, which range in size from 8 km
2 to 20 km
2. Most rockfall research studies that cover areas similar to or smaller than the areas in this study characterize the hazard based on criteria combined from both field observations and remote methods. Some researchers focus only on one site or slope at a time for which rockfall hazard is an issue. For example, both the Colorado Rockfall Rating System [
10] and the Rockfall Hazard Rating System for Oregon [
11,
12] use parameters that require field observations specific to individual rockfall slopes. Typically, such systems have a small, well-defined consequence area below the slope to consider rockfall impact to roadways and infrastructure such that the specific spatial distribution of the hazard away from the slope can be considered in an implicit, approximate manner. These systems evaluate rockfall hazard on the scale of multiple meters to a few kilometers at most. Therefore, the hazard rating for each slope can be simplified to a single number instead of mapping different hazard zones for each slope. The RHRS in [
11] created by Pierson is one of the earliest examples of a system created to characterize the hazard of individual rockfall slopes along transportation corridors from field observations. He created a chart with relevant criteria that were assigned point values for their level of hazard. The final scores could be compared so the most hazardous slopes would receive immediate attention. Since the creation of Pierson’s RHRS [
11], many other state transportation departments and research studies have used it as a base from which they tailored the criteria to their specific state or study area [
10,
12,
13,
14,
15,
16].
Other studies span entire regions similar to the scale of the Arequipa region, such as the GIS rockfall hazard modeling in Egypt by [
2], which covers an area greater than 5000 km
2. When mapping large areas, data of a lower resolution may be sufficient. Additionally, the scale of a hazard map is dependent on the intended purpose of the map, which for Omran et al. [
2] was “to examine the probable impact of such hazards on land-use planning by defining hazard patterns, securing land use, and then applying the appropriate mitigation measures”. Their final map of regional rockfall hazard was covered by designations of low, medium, or high, and they recommended that the map be used to identify potentially hazardous areas that require field investigations for more detailed evaluations of rockfall hazard. While there are substantial benefits to a regional-scale hazard map, there is a need for the development of preliminary hazard-mapping methodologies focused on the community scale as well. For example, there are studies that focus on rockfall hazard using the criterion of rockfall runout trajectory on multiple slopes above roadways and inhabited areas (instead of at a regional scale) [
6,
17]. The resulting hazard maps from these studies are helpful for direct and actionable mitigation planning, rather than requiring further modeling and field observations, as are typical next steps following the creation of a regional hazard map.
Some methods for characterizing rockfall involve using automated or semi-automated models. For example, Saroglou [
18] used GIS to model rockfall susceptibility at a national scale in Greece based on input parameters that included lithology, slope angle, rainfall intensity, and earthquake intensity. The methodology in [
18] has some crossover with other hazard rating systems, but its reliance on rainfall and earthquake data for rockfall events at a national scale is not directly applicable to community-scale characterization methodologies. Additionally, the varying roles of climate, vegetation, and other environmental variables often mean rockfall hazard assessments made for a specific region are inapplicable in other locations. Studies that depend on the frequency of rockfall are useful for regions that have an extensive rockfall inventory that includes rockfall event dates, locations, and volumes [
19]. Methods such as using mobile crowdsourcing to collect rockfall data, as described by [
20], can generate frequency and magnitude data for rockfall. For areas with limited or no existing rockfall inventories, however, it would be impractical to collect a sufficient quantity of data points to determine the overarching frequency and magnitude trends that could be used to generate hazard maps. For some areas, such as the Arequipa region, there are limited data on rainfall and earthquake intensity that can be correlated to rockfall frequency or volumes, which restricts the factors that can be included in a methodology for characterizing rockfall hazard.
A notable rockfall hazard characterization study whose methods and criteria resemble those described in this study is Stock et al.’s [
17] rockfall hazard and risk assessment for Yosemite Valley in California. They developed a system that relied on an inventory of boulders and taluses present below rockfall slopes. The process of developing this inventory included field work to record the volumes and locations of talus slopes and outlying boulders. They were also able to gather rockfall frequency information from cosmogenic exposure dating of outlying boulders. They could even determine during which season some historical rockfall events occurred. Stock et al. [
17] used STONE, a computer program that simulates three-dimensional rockfall runout trajectories, with a one-meter-resolution DEM. The one-meter DEM was a great advantage because coarser DEM resolutions have been shown to smooth topography, which may lead to inaccurate local slope angles [
21].
3. Study Area
Peru has a population of around 33.7 million people according to 2022 projections from the United Nations Population Fund [
22]. The Arequipa region, or Arequipa Department as it is sometimes referred to, in the southwest of the country only has a population of 1,382,730, according to a 2017 census [
23]. The majority of the population of the region (more than 1 million people) resides in Arequipa city, located in the east-central part of the region. Peru is divided into departments or regions, such as the Arequipa region, provinces within the regions, and districts, which are essentially towns. Each site considered for this study is a small, rural district and has a population of a few thousand, aside from Achoma, which has a population of 841 people [
23].
The Arequipa region encompasses a variety of terrain, including a coastline, rolling hills, river valleys, and high-elevation mountains, which are situated the farthest inland. To cover a variety of geologic settings, terrains, and geomorphologies that contribute to rockfall in the Arequipa region, the three sites considered are located in two different geographic areas of Arequipa. Aplao and Chaparra are located in river valleys with terrain that ranges in elevation between 600 and 2000 m, and Chivay is in the northeast mountains, at around 3600 m. These sites were chosen because of the presence of rockfall directly above community dwellings and to test the method on two different, yet prominent, geographies of Arequipa. A fourth site, Achoma, is just southeast of Chivay and was used to validate the method. Achoma was chosen because of its proximity and similarity to Chivay.
Figure 1 displays the extent of the Arequipa region and the site names and locations.
Figure 2 includes photographs taken by the research team of terrain and rockfall of each site to give the reader a sense of the terrain and geology from a different point of view.
3.1. Climate of the Arequipa Region
The entire region of Arequipa has an arid climate, but elevation plays a role in its subclassifications. The coast and highlands, including the river valleys, are considered deserts, while the mountainous region has both steppe and polar areas [
24]. In general, the Arequipa region has two main seasons—a rainy season from December to April and a dry season from May to November [
25]. The northeastern mountains where Achoma and Chivay are located experience an average annual precipitation of 550 mm [
26]. The mountains and volcanoes in the northeast reach over 6400 m, and snow lasts year-round on the highest elevation peaks, contributing snowmelt runoff to the rivers running westward to the ocean. The coastal region and the highland region, where Aplao and Chaparra are located, are arid and experience average precipitation of around 10 mm per year, which often occurs in the form of flash flooding [
25]. However, historical records have indicated annual precipitation levels up to 137.5 mm next to the coast and 412.5 mm in the highlands during very wet years [
26]. The seasonal rainfall, elevation, and snowpack in the mountains dictate the location and extent of vegetation across the region. At higher elevations (over 2700 m), native grasses, shrubs, and some trees grow in the terraced farmlands, but the amount of natural vegetation diminishes with an increase in elevation. Lower down, river valleys bloom with agricultural plots under steep hills and cliffs of bare rock and eroded soil [
27].
3.2. Geology and Geomorphology
Three quarters of Arequipa’s surface lithology consists of igneous and sedimentary rocks [
27]. Volcanic rock is mainly in the Cordillera Occidental range (high-elevation mountains), intrusive igneous rock is dispersed throughout the region, sedimentary rock is in the west-central and coastal areas, and volcanic–sedimentary rock can be found in the northwest and central-northeast part of Arequipa. The report in [
27] provides examples of ancient geomorphologic features in the Arequipa region, including plateaus, mountainous alignments, volcanic caldera, lava flows, basalt–andesitic lava fields, ignimbrite ridges, volcanic domes, stratovolcanoes, and pyroclastic flows. These features are formed from seismic activity, erosion, and volcanic processes across the region. Younger geomorphologic features in the region are due to depositional and erosional processes and include piedmonts, moraines, fluvial and glacial valleys, fluvial terraces, high plateaus, and floodplains and terraces. There are igneous and sedimentary bedrock lithologies at the three study sites and the additional validation site. However, there is minimal information available for the structural conditions of each geological unit. The influence of topography and lithology is captured within the criteria used to identify potential source zones, as explained in detail later.
3.3. Existing Hazard Inventories for the Arequipa Region
There are many geological hazards in the Arequipa Region of Peru, including landslides, subsidence, erosion, rockfall, flooding, creep, volcanic hazards, and debris flows, to name a few. One report on geological hazards in Arequipa, for example [
27], discusses an inventory and analysis of hazards from a 2013 field study that includes 2721 hazard events. The top three most abundant hazard types were debris flows, rockfall, and hillside erosion, in that order [
27]. The report is accompanied by several maps of Arequipa at a regional scale, one of which is entitled “Susceptibility to Mass Movement” and groups all forms of slope movement together, such as landslides, rockfall, debris flows, mudflows, hillside erosion, creep, and soil liquefaction. It designates susceptibility ratings as low, medium, and high across the entire Arequipa region [
27]. Rockfall is mentioned in the “Alta” or “High” category, which is the second highest category. Aplao is one of the localities listed as being within the zone of high susceptibility. Rockfall hazards are not mentioned specifically in any of the other categories, even though rockfall source zones were inventoried in areas designated as moderate and very high. This map is best used for general hazard susceptibility, as it does not indicate susceptibility to specific types of mass movement.
Other inventories such as [
28] do not include the extent of rockfall events, nor do they include inventories of source zones and runout zones for rockfall hazard. Furthermore, the occurrences listed for rockfall hazards outside of densely populated areas are sparse and mostly include recent, large events that had a significant impact, likely because these events gained enough attention or were noticeable enough to be reported. It is unknown how aware Peruvian citizens are of this online hazard inventory, or if they play a role in reporting or managing hazards. The Regional Government of Arequipa created several geological and geomorphological maps of the Arequipa region [
29]; however, these maps only identify general processes like erosion, slope movement, and weathering. There is a very limited inventory of rockfall for the Arequipa region and it does not include frequency or magnitude information for the majority of events [
27]. Since specific information is so limited, there is a need for more extensive and precise rockfall hazard inventories for the Arequipa region.
4. Methods
4.1. Available Data
The data used in this study are limited to those which are publicly available. The communities that are being mapped in this study are assumed to only have access to data and information that are publicly available, and the study’s methodology was developed with this limitation in mind. Google Earth is a reliable and publicly available software for viewing up-to-date imagery. Google Earth compiles high-resolution imagery from multiple sources with the exact resolution depending on the source. Google Earth aerial imagery was used to complete visual inventories, relying on the imagery viewer and mapping tool with scenes from June 2020 and January 2022 [
30]. Gaps in coverage were filled in with high-resolution scenes from the closest dates available, which is an automatic feature of the Google Earth imagery viewer tool. Imagery sources include Maxar Technologies with a resolution between 15 and 30 cm, and the Airbus Constellation of satellites with a resolution of 30 cm [
31,
32].
DEMs were needed for the input parameters of slope angle and slope geometry. Because of its consistent availability across the region and its public availability, a 30 by 30 m resolution ASTER Global DEM (GDEM) from April 2021 was used to determine the corresponding criteria [
33]. Publicly available data were chosen to design a method that is accessible to those who do not necessarily have the funds to purchase high-resolution data. It was important to the researchers that the communities who would be receiving the maps would also be able to access the data. The ArcGIS Pro 3.0.0 application was used to perform calculations with the DEM data and implement the hazard zone equations to generate final maps.
4.2. Rockfall Inventories
The method by which the rockfall hazards were inventoried began by using aerial imagery to locate rockfall and source zones at each of the three study sites. Next, the inventories were verified and modified based on field observations at the three sites. Finally, the inventories were revised digitally based on the observations, and inventory maps were created for each of the three sites. Complete inventories for the three study sites were created by marking fallen rocks and deposits visible below source zones in aerial imagery using the Google Earth Pro mapping application. For the validation site of Achoma, only source zones were inventoried, and their extents were verified in the field. Source zones, deposits, and runout zones were identified based on the criteria in
Table 1.
Lan et al. [
34] identified source zones in their study as slopes that exhibited sharp topographic contrast, steep slope angles, and minimal to no vegetation cover. In this study, the source zones were identified as visible rock outcroppings at the top of slopes with a slope angle of at least 20 degrees at the level of resolution available. There may be ledges within the slope that are at a higher angle and therefore more susceptible to rockfall. Due to the scale of the mapping, some source zones include short sections of a slope with no visible fractured bedrock when they exist vertically between areas of exposed bedrock or cliffs from which rock has fallen. Because of the limited DEM resolution, source zones that were less than three meters above the base of the slope were not included.
The recorded rockfalls could have occurred at any time, so there is no frequency or timing information, although the relative prevalence of rockfalls in different areas can be used to infer frequency. Objects that could not be confirmed as rocks from aerial imagery or Google Street View were indicated as “potential rocks”. The resolution of the imagery limited the size of the rocks that could be identified to boulders (approximately 1 m diameter) and deposits of smaller rocks (approximately 2 m3). Runout zones are located below the source zones and encompass the entire area where rocks and deposits have accumulated from the corresponding source zone.
In August and October 2021, the research team visited the Arequipa region to compare the rockfall inventories of the three study sites with field observations. In the field, the inventories were updated to include the observations, and photos were taken of the sites for future reference. Rockfall source zones were verified by visually confirming pathways for fallen rock below the identified source zones. The rockfall inventories were updated in Google Earth based on the observations of the field team. These are hazard
inventories because no calculations were performed to determine hazard ratings in this stage.
Figure 3 shows part of the inventory created for Aplao.
Runout distances were measured for all the fallen blocks at the study sites by drawing a line from the fallen rock to the bottom of the nearest source zone, following the most likely path the rock would take. The most likely paths were identified by looking for grooves that other rocks created in the slope, topographical features and barriers that would inhibit and redirect a rockfall path, and changes in slope angle that may redirect a rock. Statistics were recorded such as the frequency of rockfall from a certain distance from the bottom of the source zone and the percentage of fallen blocks passing certain distances in the runout zone.
4.3. Rockfall Hazard Criteria
The published literature was used to identify the parameters commonly used to indicate propensity to rockfall (
Table 2). Some of these studies adapted Pierson’s [
11] RHRS, while others tested entirely new systems, scales, and methods. The aim of this study was to create a methodology that could be implemented remotely, so parameters that could be measured remotely and were available online were prioritized. Furthermore, parameters used by a large number of studies were considered to be more important than those seldom used. Some of these included slope height (used in 75% of the systems), slope angle (70% of the systems), lithology (50% of the systems), joint persistence (40% of the systems), and joint orientation and slope roughness (also noted as “launching features” or “differential erosion”) (used in 45% of the systems).
Notwithstanding the logic above, common use does not automatically imply suitability for all sites. For example, if a site has highly variable vegetation cover, then vegetation may be a more important criterion for that site. Based on our judgment for this region, the hazard criteria chosen for this study are source zone height, source zone slope angle, rockmass joint persistence and orientation, rockfall trajectory paths, and rockfall runout distance. These criteria relate to either the “generation” (block release from slope) or “transport” (movement downhill) of rockfall and were used to determine hazard scores for both the source zones and the runout zones, respectively.
Certain criteria normally expected to influence rockfall were not included in this analysis. For example, the size of our sites is small enough that any seismic activity would have the same amount of impact across the entire site, so the presence of faults and other seismic hazard indicators were not included. Since the low resolution of the available DEM smooths the small-scale features that would dictate surface roughness, this parameter was not included. Vegetation was not included because there is a distinct lack of large vegetation, such as trees and tall bushes, that would influence the downslope path of rocks throughout the region.
Although lithology is a common criterion in RHRS, most studies determine rockfall susceptibility of formations specific to their study site. For example, Dorren and Seijmonsbergen [
21] summarized the rockfall susceptibility of various formations and correlated that to general rock types. However, the low, medium, and high ratings in their study overlapped some rock types such as limestone, sandstone, and schist, while all igneous rocks were rated as having “high” rockfall susceptibility. In Dorren and Seijmonsbergen’s [
21] study, there were not enough definitive conclusions of the rockfall susceptibility ratings of rock types to be able to apply their findings to this study. The geological map of the study area does not include structural information such as joint orientations [
27]. To include lithology in this study would require extensive field work to gather data on specific rock types located in the source zones, which is beyond the scope of this study, which is to rely solely on remote sensing data for determining inputs.
While a larger source zone area does not necessarily equate to more actual rockfall, it implies more
potential rocks that could release, so source zone height is included in this study as a surrogate for the size of a source zone. This is because the influence of source zone length is directly accounted for in the way in which the source zone hazard ratings are projected downslope (see
Section 4.4.2), and ultimately are represented spatially in the runout zone hazard ratings on the map.
Precipitation, or “water present on face”, as indicated in some studies, was not included as a criterion because the Arequipa region has a relatively uniform dry, arid to semi-arid climate, and the precipitation in all areas typically comes in the form of flash rainstorms [
26]. Irregular rainfall has the ability to trigger rockfall and it creates an environment that is more susceptible to rockfall [
18]. However, precipitation amounts vary consistently across the Arequipa region with the time of year, and the maps we produced using remote methods represent a static hazard throughout the year. Likewise, criteria that fluctuate based on the time of year such as freeze–thaw cycles, temperature, and precipitation were not considered. However, Peru generally has a wetter time of year that ranges from December to April, so it can be expected that more rockfall will likely occur during this period.
4.4. Hazard Mapping Approach Summary
To generate the hazard score for a site, rockfall hazard is split into two main processes—generation and transport—with the scoring process summarized in
Figure 4. Generation involves the formation and release of rock fragments from source zones. The rock fragments are formed when they detach along discontinuities in bedrock [
17], and then initially move by sliding, toppling, creeping, or falling from the source zones [
43]. Source zones are the most dangerous areas in a rockfall hazard zone because of the instability of the bedrock and typically steep slopes. The study in [
2] states that “when the rocks detached [
sic] from the steep slopes or cliff face, then the falling material moves abruptly with high speed downward to the slope base”, and “this type of movement may start at the source zone which is part of the high-risk zone, especially at the high elevation and high steeply sloped areas”. Note that all source zones are relatively high-hazard, and that rockfall hazard increases with an increase in slope angle [
2]. Thus, the source zone hazard rating influences the hazard score of the runout zone where the rocks are falling.
Since the generation component of rockfall hazard describes the likelihood that blocks of a given size will fall, it depends on the source zone height (which influences the total area of slope available to generate rockfall), the slope angle (which influences the number of kinematically feasible release mechanisms), and the rockmass mechanical characteristics. Because blocks that are transported below a given portion of a slope must have been generated by instability within the source zone, the hazard intensity of the source zone where rockfall originates influences the hazard of the runout zone where rockfall is transported, and it should be a factor in the hazard rating.
The transport of rockfall occurs when a rock fragment falls, bounces, slides, or rolls down a slope due to gravity after it initially releases from the source zone. The rock will lose energy when it makes contact with the slope, so lower slope angles near the base of a slope will slow down the rock, whereas steeper slopes allow for longer transport distances of rockfall fragments and allow more time for the rock to accelerate [
4]. The path a rock will follow is dependent on the terrain geometry, such that a chute or gully will have more control over the direction of a rockfall path than a smooth, open slope.
The determination of some inputs such as source zone height and slope angle can be automated using ArcGIS Pro tools, but other inputs such as structural condition and the source zone inventory must be manually input into the program. As a consequence, a semi-automated method using tools from ArcGIS Pro through the ArcPy Python package was developed to identify the boundary of the source zone (details are provided in [
44]).
Figure 5 shows a schematic of this process.
4.4.1. Source Zone Hazard Score
For the source zone, three criteria were assigned from which a score and hazard rating are determined. For this method, source zone height is the mean elevation of the bottom of the source zone (furthest downslope boundary) subtracted from the highest elevation of the source zone. The slope angle was determined from the DEM, and it is expressed in a raster format as a score per cell. A structural condition score was designated per site based on the dominant joint persistence and orientation present. For the total score, each criterion is assigned parameters that correspond to a number of points, as shown in
Table 3. The point scale is based on the RHRS created by Pierson [
11], where a higher number of points is a higher hazard, and the parameters listed are those best suited for the Arequipa region and the type and scale of hazard rating used in this study.
The following sections present justifications for the parameters chosen for each criterion and an explanation of how the scoring is implemented in ArcGIS Pro.
Source Zone Height
Almost all the RHRS methods reviewed, even the ones that rated rockfall hazard next to roadways, treated slope height as the entire height from the top of the source zone to the bottom of the slope, and they did not separately consider runout zones (i.e., a separate “transport” category for hazard assessment). This was suitable for their methods due to their smaller scale and the in-person evaluation of each slope. For an automated method that considers the spatial distribution of hazard below each slope, the location and extent of the source zone is important for determining the location of the runout zone. Importantly, the height of the source zone represents its size and the amount of exposed rockmass that could potentially produce rockfall blocks for any given bottom-of-slope position in plan-view.
Source zone height can vary considerably within a site. For example, Chivay has two distinct areas of rockfall hazard with drastically different source zone heights and extents. One is a narrow band at the edge of a mesa, while the other consists of rock outcroppings across the mountainsides that vary in size. For this reason, the source zone height is represented as the difference between the maximum source zone elevation and the average elevation of the bottom of the source zone line (indicated by the cyan line in
Figure 5). The height is calculated using Equation (1), where
hmax is the maximum elevation of the source zone,
xn represents the elevation of each point along the line that designates the bottom of the source zone,
n is the number of points along the line, and
H is the source zone height used to determine its hazard category.
After calculating the source zone height, the category of points it corresponds to for each site was determined. A single score was used for each source zone as an input to determine the overall source zone hazard score. If there were multiple source zones at a site that had different structural condition scores, then each source zone was treated as a separate rockfall zone and the scores were calculated individually. The height values were modified from Pierson’s [
11] RHRS because his system considered the height of the entire slope instead of just the source zone height, and the sites in this study are larger than the individual slopes above highways rated for many of the other studies. Furthermore, the height scores that [
11] uses represent a relationship where taller slopes mean a larger distance for transport. In this context, however, slope height is only considered with respect to its impact on the generation process, namely, as slope height increases, there is more rockmass from which potential blocks could release from a given two-dimensional cross-section. This is why the slope height scoring was changed from Pierson’s [
11] original values and scaled to designate greater slope heights in the higher-scoring categories; specifically, independent of transport factors, the hazard score should increase approximately linearly proportionally to the slope height, as the primary generation-related influence of slope height was to affect the relative size of the slope from which rockfall could be related. For the source zone heights in this study, meters were used to maintain consistency across datasets. Pierson [
11] used 25 feet for the maximum slope height for the three-point category, and 25 feet is almost 8 meters, a nice round value. The same logic applies to the other categories.
Source Zone Slope Angle
There are many ways to classify the hazard level of slope angle in an RHRS, and various studies have used sensitivity analyses, statistical analyses, and the modeling of rockfall events to relate hazard levels to slope angle ranges [
2,
3,
13,
17,
34]. The main takeaway from these studies is that there is no universal slope angle range for source zones. Stock et al. [
17] performed modeling of rockfall trajectories in Yosemite Valley and considered any slope greater than or equal to 60° in angle to be a potential source zone. Even though the authors acknowledge that rockfalls have occurred in that area on lower angle slopes, their justification is that in more than 70% of recent rockfall events in Yosemite Valley, rocks had released from slopes greater than 60°. Dorren and Seijmonsbergen [
21] chose to designate source zones as bedrock with an average slope gradient equal to or greater than 40°, which was based on the amount of rockfall activity that occurred below slopes of different mean gradients. According to a statistical analysis of rockfall in their study area, Lan et al. [
34] determined that “50% of the rockfalls occur at slope angles less than 40°, and over 84% of the rockfalls occurred at angles greater than 30°”. Because rockfall is less likely to but can release from slopes below 40°, low-angle source zones were included in the criteria, but assigned a low hazard score of three points. The slope angle of 70° was chosen to represent the upper limit of the 9-point category and the lower limit of the 27-point category to include the 60–69° slopes in the moderate category since 60° was the lower limit for Sock et al. [
17]. While multiple studies argue that the independent variable of the slope angle has the highest influence on rockfall hazard assessment, this is based on its influence on the downslope transport of a rock [
4]. For a rockfall source zone, the influence of slope angle is not nearly as potentially large on the generation of rockfall, so the criteria have been adjusted accordingly to prevent any slope from receiving the maximum of 81 points.
To score the source zone based on slope angle, the “Slope” tool in the Spatial Analyst toolbox in ArcGIS Pro was used to create a raster of the slope angle in the source zone. Then, the raster was reclassified to reflect the score categories that correspond to the slope angle ranges listed in
Table 3. When adding the scores for all three source zone criteria, the slope angle score dictates the variability in the score per cell across the source zone. The 70th-percentile hazard score is thus dependent on the highest slope angle score that makes up at least 30% of the source zone area. Although the selection of this percentile criterion is somewhat arbitrary, it is noted that slope angle is often not the determining factor in the final hazard score since the maximum score for slope angle cannot exceed 27 points. For example, an overall source zone score of 104, with 3 of those points due to slope angle, is not much lower than an overall score of 110, with 9 of the points due to slope angle. Even if the slope angle score was changed to 27 points, the source zone hazard score in this example would only increase to 128. In other words, the uncertainty in the appropriateness of the specific percentile value used does not substantially influence the overall hazard rating score.
Structural Condition
The structural condition criterion is pulled from the Colorado Rockfall Hazard Rating System modified by [
10], which was originally based on Pierson’s RHRS [
11]. The point category for the structural condition of a site is based on the persistence and orientation of the majority of the source zone rockmass.
Figure 6 includes drawings of the general joint orientations and joint continuities included in
Table 3. A single score was assigned to a site for the structural condition, which was determined by examining high-resolution aerial imagery in Google Earth. If there were distinct differences in structural condition across a site, as was the case in Chivay, then the rockfall hazard zones were separated and assigned different condition scores per zone.
4.4.2. Runout Zone Hazard Score
For the transport component of the workflow, the three inputs into the Equation are the corresponding source zone score, the amount of rockfall passing distances downslope expressed as a percentage, and a measure of the density of rockfall trajectories passing through each cell. The transport of rockfall occurs in the runout zone, so the hazard score for this portion of the area is referred to as the “runout zone score”, calculated per cell based on Equation (2), where
z is the source zone hazard score,
p is the percentage of rockfall likely to pass through a cell,
t is the “sum of rockfall trajectory paths” for a cell (see Section “Trajectory Sum”), and
r is the runout zone hazard score.
Figure 7 provides a visual interpretation of Equation (2). The source zone score is a constant input. The percent passing and trajectory path sum components vary from cell to cell, which results in a raster that contains a hazard score per cell based on the DEM’s resolution. More rockfall paths in a cell means a greater score for the “Rockfall Trajectory Path Sum”. To maintain higher hazards in gullies and basins where more rockfall paths will overlap without causing excessively high scores, the runout zone score is capped at the source zone score.
Figure 7 shows how the source hazard score, the percentage of rockfall passing, and the sum of rockfall trajectory paths are multiplied and result in a distribution of hazard ratings in the right of the image. More detail on the individual components will be provided in the following sections.
Trajectory Sum
The trajectory of a rock moving down a slope is controlled by the slope’s geometry [
34]. In their rockfall simulations, Dorren and Seijmonsbergen [
21] utilized rockfall trajectories by identifying the “fall direction” of a rock moving down a slope as “the direction of the steepest descent towards a neighboring cell,” which built on the method from [
45]. Following this general approach, in this study, the trajectories were estimated with the “Cost Path as Polyline” tool in ArcGIS Pro. The bottom of source zone line was used as the starting point for trajectories, while a “flow direction” raster that was derived from the DEM was used to force each trajectory to follow the steepest path downslope. The trajectory calculation step required a manual operation to remove the area where the trajectories should not pass, or where a rock (instead of a water particle) would not reasonably travel. This is because the “Flow Direction” tool in ArcGIS Pro is intended to be used for the hydrological analysis of water. After the trajectories were mapped, the sums of overlapping trajectories were calculated and converted into a raster, so the raster had a trajectory sum score per cell. To ensure that an absence of trajectories did not equate to zero hazard on the hazard map due to uncertainty associated with the relatively simple trajectory estimation methods, all cells with a value of zero were converted to have a value of one. Additionally, the trajectory component can cause extremely high hazard scores in very small areas (channels/chutes) at some sites, and these scores end up being outliers from the main set of hazard values. To avoid excessively large hazard scores in such areas, the runout zone hazard score had a cap, such that the highest score for a runout zone could not exceed its corresponding source zone hazard score. The trajectory sum raster was input into the source zone hazard via Equation (1).
Percent Passing of Rockfall from Bottom of Source Zone
From the developed rockfall inventories, the distance from each fallen rock to the nearest point along the bottom of the source zone was measured, while following the most likely path each rock would take. This process was performed for all rocks inventoried for the sites of Chivay, Aplao, and Chaparra. By using the inventory of runout distances to count how many rocks travel beyond certain distances away from the bottom of the source zone, the influence of factors known to influence travel distance for rocks (slope angle, slope height, and surface roughness) was embedded in the empirical runout distance criterion.
Figure 8 displays the curves generated from the runout distance data for Aplao, Chaparra, and Chivay, along with the corresponding equations for their trendlines and the R
2 values. For all three sites, the majority of the rockfall runout distances were under 100 m. The maximum rockfall runout distance for all sites was over 700 m, but very few rocks traveled that far.
The rockfall percent passing criterion was input into the runout hazard from Equation (2) as a raster where each cell is a decimal value to represent the percentage of rockfall that passes the distance that cell is from the bottom of the source zone. Most previous studies that included rockfall runout distance criteria in their RHRS did so by modeling rockfall runout using a DEM as the surface raster [
4,
17,
34,
42]. In this study, instead of simulating rockfall runout, the “Cost Distance” tool in ArcGIS Pro was used to create a raster so the value of each cell is its distance away from the bottom of the source zone. Then, the equation of the trendline for the respective site (
Figure 8) was input into the “Raster Calculator”, and the Cost Distance raster was input for
x to create a percent passing raster, which is an input for the runout zone hazard in Equation (2).
6. Discussion
The rockfall hazard characterization methodology developed in this study uses solely remote methods to characterize the hazard and develop a preliminary hazard map for different areas in the Arequipa region. While there are obvious advantages to being able to manage rockfall hazards at remote locations or those without extensive field information, there is a trade-off with accuracy and confidence in the results. For example, this study can be compared to a study with similar goals but much more site information, the rockfall hazard and risk assessment for Yosemite Valley in California by Stock et al. [
17]. They developed an inventory of boulders and taluses that existed below rockfall slopes, along with rockfall frequency information from cosmogenic exposure dating of outlying boulders to create one of the few other studies that uses rockfall inventories to inform and improve rockfall hazard mapping. In contrast, the rockfall inventory created for this study included a much larger area overall, and did not treat boulders and taluses differently. Stock et al. [
17] used a one-meter-resolution DEM to create a slope angle map and to inform STONE, a computer program that simulates three-dimensional rockfall runout trajectories. This study used flow direction derived from a 30 m DEM to map rockfall trajectories. Stock et al. [
17] had the added benefit of being able to draw on historical reports of rockfall in the Yosemite Valley, while there have been few reports on rockfall, and no rockfall modeling studies, for the Arequipa region. Stock et al. [
17] designed their method specifically for the steep rock walls in Yosemite Valley, while the criteria for this study were chosen to address the wider range of topography and geology for the Arequipa region. We suggest that these two studies represent two different end members of approaches using rockfall inventories to supplement and calibrate rockfall hazard analysis. In one case, the availability of historical records and high-resolution DEMs allows for high-confidence hazard mapping for a small area, but at the expense of requiring more detailed field and laboratory analyses. In our case, information typically available worldwide can be used to create and calibrate a system that can be applied over larger areas with little field work, but lacking resolution. The latter provides a tool for preliminary analysis and offers distinct advantages for developing countries and remote areas. The former is better suited for higher-vulnerability or concentrated features where resources can be focused on smaller areas.
Another study by Depountis et al. [
6] used a GIS-based method to perform a hazard assessment in mountainous sites in western Greece, focusing on a community and a transportation corridor. They also developed an inventory of rockfall using high-resolution aerial imagery, but they had a more robust backlog of historical imagery to supplement their inventory. They identified potential points of release for the inventoried rocks as well. They were able to use these data along with a digital surface model (DSM) to model rockfall behavior such as velocity, kinetic energy, bounce height, trajectory paths, and runout distances. They had the advantage of using a DSM with a 5 m x 5 m resolution, which is significantly higher than the DEM used in this study. Despite their higher-resolution terrain data, Depountis et al. [
6] did not use slope geometry to determine trajectory paths and instead assumed straight lines from the source point to the end point for each rock. This left out information that could be used to inform the areas that rocks pass through the most, arguably increasing the hazard in those areas, and was something that increased the accuracy of the hazard maps for this study despite the lower-resolution DEM. One of the main differences between their methodology and this study is that Depountis et al. [
6] measured the size of fallen rock blocks in the field, which was not performed in this study since the inputs needed to be remotely determined. They used the block size data to calculate the kinetic energy they produced, which were designated into intensity classes of rockfall. The intensity classifications were then multiplied by probability of occurrence grades to result in five possible hazard categories from low to very high. While the resulting hazard maps from [
6] are visually similar to the ones developed in this study, they acknowledge that further research could incorporate RHRSs or other slope stability approaches to increase the reliability of their estimations. While they did consider lithology, source zone structural conditions such as joint continuity and orientation were not a part of their assessment. The number of rocks included in their inventory also appeared to be much lower than was used in the present study, which could result in greater room for error in their calculations. However, the present study could benefit from the use of intensity classes as part of the criteria for hazard rating, but it would require data from historical events that are not currently available, as well as field observations of block size, which are not within the scope of the proposed methodology.
Several other RHRSs were described in
Section 2, including one created by Pierson in [
11] for Oregon’s Highway Division, widely considered the first of its kind for characterizing rockfall slopes in transportation corridors. Its applicability is to discrete slopes above roadways and requires extensive field observations and measurements. While the method in this study drew inspiration from the point system and several criteria in Pierson’s system, it was created for the remote characterization of entire communities that encompass multiple slopes at different aspects and angles. Pierson’s [
11] system also considers risk by including the criterion that is a measure of how often a vehicle exists under the rockfall zone. As mentioned in the introduction, risk and hazard are two different classifications and require different approaches to the scoring and identification of criteria. An advantage of this remote method is that it can characterize a larger area in less time than would be required for lengthy field investigations of individual slopes. Oftentimes, sites are far away and hard to reach, which means there are costs associated with travel and lodging while performing field investigations. A remote method like this one eliminates travel-related expenses, and the sourcing of solely publicly available data makes this method more accessible for low-income communities. However, to ensure accuracy, the maps produced from this method would benefit from field validation and the method would need to be modified based on findings in the field. For this reason, this method and the resulting hazard maps should be treated as preliminary in the development of any hazard mitigation plans or policy discussions.
Map scale makes a significant difference in the ways a given map can be applied and by whom. In contrast to methods that remotely characterize rockfall hazard at regional and national scales, such as those by [
2] and [
18], the method in this study uses criteria that identify the trajectories and runout distances of rockfall. This rating method is thus more applicable for community use since the maps show more detail of inhabited areas, including houses, roads, and farmland. Rockfall inventories provide useful data for any scale of map, including regional ones like those in [
19]. However, many countries do not create national inventories of small-scale rockfall due to the intensive time and effort they take to maintain, especially in rural areas. Community members must typically rely on local and historic knowledge passed on from those who have lived in the area longest. The novelty of this study is the creation and use of a rockfall inventory without historic information. This creates an opportunity to characterize rockfall hazard in areas with limited or no rockfall inventories.
Certain aspects of the criteria make this method more suitable for areas with minimal vegetation and low precipitation. Otherwise, there is nothing inherently site-specific about the method in this study. As long as the imagery is available, the method can reasonably be applied to communities in regions that have a similar climate to Arequipa such as western Peru, northwestern Chile, parts of eastern Argentina, southwestern Nevada, southeast California, parts of northern Mexico near the New Mexico border, central Iran, western Namibia, and parts of South Africa, to name a few. Similarly, if a high-resolution DEM or DSM is available, the method becomes more precise and easier to automate. For example, the rockfall trajectory paths would require much less manual “clean-up”, which speeds up the method. Organizations and government agencies concerned with disaster preparedness and management would find this method particularly useful when creating preliminary rockfall hazard maps of difficult-to-access communities. It would provide them with data that could be used to develop risk maps and allow them to prioritize sites that have greater exposure to rockfall. Risk maps are a logical next step to addressing hazards once they are identified. For example, simple consequence maps can be made that correspond to population density or infrastructure value. Combined with the information from the hazard maps, they would be able to create rockfall risk maps to assess what parts of a community are at greatest risk to rockfall impacts.
An overarching limitation of this rockfall hazard estimation workflow is that it is only as precise as the coarsest input data. The data for the criteria of source zone height, slope angle, and trajectory sum rely heavily on the DEM. The DEM used for this study has a resolution of 30 m, which is considered a low-resolution dataset when working with areas of 10 to 20 square kilometers. There were limited options for the DEM input, and a higher-resolution DEM would make the resulting hazard ratings more precise. On the other hand, the aerial imagery available from Google Earth Pro was of a very high resolution, and different imagery could be used from different years, which allowed the clearest imagery for a particular site to be found. This was helpful for inventory development, but there were still many rocks that were not included because they could not be clearly seen in the imagery due to their size, or an object could not be confirmed to be a rock. Developing an inventory in the field is another way to approach this aspect of the study and include more rocks of varying sizes, but that would require access to the site, extensive time in the field, and travel-related costs.
7. Conclusions
The goal of this study was to develop and test a method to use remote and widely available information to identify rockfall hazards across large and remote areas. We began with published rockfall hazard rating systems to identify parameters previously used that could potentially predict rockfall source and runout zones and could also be remotely determined. Specifically, we selected parameters that reflect the arid terrain and seasonal rainfall in the Arequipa region in southern Peru. For rockfall generation, the criteria that fit this description are slope angle, rockmass structural condition, and slope height. For rockfall transport, the criteria are rockfall runout distance, trajectory path sum, and the hazard level of the connecting source zone.
These parameters were applied to a GIS model to characterize rockfall hazard in three high-altitude Andean communities—Aplao, Chivay, and Chaparra. While other studies in the literature have tested GIS-based approaches for rockfall hazard estimation using remote methods, the work in this study to incorporate a rockfall runout distance criterion, based on rockfall inventories developed using Google Earth imagery, represents a unique approach. For example, quite often, data on the frequency and magnitude of rockfall are not available, especially in rural areas. The rockfall runout distance data collected from the rockfall inventories provide information about rockfall transport and can be used to estimate hazard levels of runout zones. The runout distance data were collected by manually identifying fallen rocks and deposits in the aerial imagery, and then measuring the distance from the rock to the nearest point at the bottom of a source zone while following a trajectory that the corresponding rock would most likely take. This method was calibrated by comparing hazard scores from three rockfall areas with different combinations of source zone criterion ratings and determining the 50th- and 70th-percentile scores as break points to designate three hazard rating ranges. From these ratings, final hazard maps were created with low-, medium-, and high-hazard zones that are simple for users to interpret along with corresponding rating descriptions of the expected rockfall activity. While there are limited options for verifying that methods based on parameters determined using remote methods can predict rockfall hazard in the Arequipa region, the density of the mapped inventory of rockfall in Aplao matches closely to the descriptions of rockfall activity in each corresponding hazard zone. The remote, semi-automated method developed for predicting rockfall hazard for the Arequipa region is an effective way to develop hazard maps in areas with limited existing data and that are difficult to access. This method could be applied to other areas in Peru that have arid climates near the coast. It could also be applied to regions and countries worldwide with a similar climate to the Arequipa region.
A hazard assessment is a starting point for researchers, governments, emergency and hazard-focused NGOs, and communities to understand and identify areas that are in danger. Once a hazard rating is made, a risk assessment is the logical next step to identify locations of high risk, inform those who live, work, or own property there, and implement mitigation. The appropriate form of rockfall mitigation is highly dependent on the site conditions and a professional in slope stabilization, such as a geotechnical engineer, should be consulted.