Interactive comment on “ Automated landform classification in a rockfall prone area , Gunung Kelir , Java ”

Comments by I S Evans, for Earth Surface Dynamics, on Samodra et al. – Automated landform classification in a rockfall prone area, Gunung Kelir, Java. G. Samodra, G. Chen, J. Sartohadi, D. S. Hadmoko, and K. Kasama GENERAL COMMENTS: This paper classifies slopes in a small area of Java into 7 units of the 9-unit slope model (catena), using unsupervised classification. An inventory of 521 rockfall deposits (ranging from 0.0018 to 3627 m3) is related to four of these units, and power functions are fitted to each subset. Some of the writing is repetitive or vague. Re-writing of the paper is needed, more


Introduction
In attempts to study and understand landform, people have tried to map and document landform features since long time ago. Summerfield (1991) explained that the first attempt of human to document landform has been started from the age of Herodotus 25 (5th century BC) and Aristotle (384-322 BC). It was described in the simple way. In 20 the genesis focused on geology and geomorphology structure, process, and the development phase of landform, whereas Penck implicitly classified landform with the consideration of the forms or shapes of terrain features, the processes involved in their formation, and the grouping of associated forms into typical assemblages (Beckinsale and Chorley, 1991). Later on, Lobeck (1939) classified landform into constructional and destructional landform. Demek and Embleton (1978) made a more systematic landform classification for geomorphological mapping purpose. The landform classification was based on morphology, relief and genesis, age (morphochronological point of view), morphoassociations and morphoregions. Meanwhile, Verstapen (1983) and Zuidam (1983) classified landform based on the genesis as follows: structural, volcanic, fluvial, 15 marine, glacial, eolian, solutional, organic, and denudational.
The latter classification has been widely used in Indonesia. It would be well performed to classify landform only based on the genesis especially on the small-scale map. However, it is necessary to add the other landform information in order to map geomorphological features in the medium to large scale. In a later development of Introduction phometry as a preliminary tool for risk assessment, the spatial planning manager can make a balance between minimizing risk and promoting some development priorities. Risk can be defined as "the expected number of lives lost, persons injured, damage to property and disruption of economic activity due to a particular damaging phenomenon for a given area and reference period" (Varnes, 1984). The definition is orig- 5 inally used to describe landslide risk. Later, the terminology is used for all types of mass movement including rockfall. The word "rockfall" is often distinguished from more general landslide phenomena due its typical material, size and failure mechanism. It is defined as rock fragments (Hungr and Evans, 1988) with size from a few dm 3 to 10 4 m (Levy et al., 2011) started by the detachment of blocks from their original posi-10 tion (Crosta and Agliardi, 2003) and followed by free falling, bouncing, rolling or sliding (Peila et al., 2007). Rockfall risk can be expressed by the simple product of temporal probability, spatial probability, reach probability, vulnerability and value of the element at risk (Fell et al., 2005;Westen et al., 2005;Agliardi et al., 2009) as follows: 15 where P (L) j km is the temporal probability (exceedance probability) of rockfall in the magnitude scenario (i.e. boulder volume) class j and crossing landform k for different period m; P (T |L) i j is the probability of the rockfall in the volume class j reaching the element at risk i ; P (I|T ) i is the temporal spatial probability of the element at risk i , V i j is the vulnerability of the element at risk i to the magnitude class j and E i is the economic 20 value of the element at risk i . Based on the Eq. (1), the magnitude and exceedance probability of rockfall are diverse in time and places. The 9-unit slope model (Dalrymple et al., 1968)  expose the spatial distribution of potentially high damage of elements at risk affected by rockfall. Thus, selection of preventive mitigation measure type, structural protection location, and structural protection dimension should be supported by rockfall risk assessment based on landform analysis. Traditionally, landform analysis (delineation and classification procedure) is based 5 on the stereoscopic technique of aerial photo and field investigation. This method is very common in Indonesia. It has been applied for soil mapping, land evaluation analysis, land suitability analysis, spatial planning, and so on. There is also mentioned in Indonesia's National Standard document of Geomorphological Mapping that the technical requirement for geomorphological mapping is an interpretation of remote sensing 10 data combined with field measurement (SNI, 2002). The standard landform classification in Indonesia is based on the ITC Classification System (Zuidam, 1983). However, the traditional method of landform classification requires simultaneous consideration and synthesis of multiple different criteria (MacMillan and Shary, 2009) and the quality depends on the skill of interpreter. The development of landform classification has been 15 applied mostly in soil landscape studies. Thus, we try to automated classify landform based on the 9-unit slope model which is appropriate to rockfall analysis. Even though, the 9-unit slope model is significant for pedogeomorphic process response (Conacher and Dalrymple, 1977), it is also capable to explain rockfall deposition. 20 Gunung Kelir is located in Yogyakarta Province, Indonesia. It lies in the upper part of Menoreh Dome that is located in the central part of Java Island (Fig. 1). The area is dominated by Tertiary Miocene Jonggrangan Formation that consists of calcareous sandstone and limestone. Bedded limestone and coralline limestone which form isolated conical hills may also be found in the highest area surrounding the study area. well explained by Bemmelen (1949). The evolution of Kulon Progo Dome was started by the rising up of geosynclines of South-Java in Eocene Period. It made the magma of Gadjah Volcano consisted of basaltic piroxene andesites reached up to the surface. Then, it was followed by the activity of Idjo Volcano in the south with more acid magma consisted of hornblende-augite andesites and dasites intrusion. After the strong de- Gunung (Mountain) Kelir, of Javanese origin, literally means a curtain that is used to perform wayang (Javanese traditional shadow puppet). Its toponym describes a 100- 15 200 m high escarpment that has a maximum slope nearly 80 • . The complex of Gunung Kelir consists several generic landforms which are prone to rockfall. Its slope gradient varies between 0 • and 80 • , meanwhile mean of slope gradient is 23.14 • with the standard deviation 13.05 • . Altitude ranges from 297.75 to 837.5 m. There are 152 buildings exposed as elements at risk on the lower slope of the escarpment.

Data and methods
Rockfall risk analysis requires assessment of susceptibility and identification of an element at risk. To portray the susceptible area, geomorphologist opinion is commonly used to classify landform through interpretation of aerial photos and field survey. However, subjectivity of investigator hinders application of this method. Therefore, unsu- 25 pervised landform classification of 9-unit slope model is applied in the present study.
The main objective of this study is to provide automated landform classification partic- ularly for rockfall analysis. To achieve the primary objective, several works ( Fig. 2) are conducted in this study: (1) fieldwork, (2) DTM preprocessing, (3) DTM processing, (4) rockfall modeling, (5) landform classification based on fuzzy k-means, and (6) rockfall volume statistic. Fieldwork was intended to identify rockfall boulders and elements at risk. A field in-5 ventory of fallen rockfall boulders of different size has been done to obtain the spatial distribution and dimension of rockfall deposition. The dimension and potential rockfall source were determined to simulate rockfall trajectory, velocity, and energy. The buildings on the lower slope of the escarpment were also plotted in order to obtain the spatial distribution of elements at risk. Finally, DGPS profiling was conducted to improve the 10 performance of DTM. The objective of DTM preprocessing was to improve the quality of DTM-derived products. We applied DTM preprocessing proposed by Hengl et al. (2004) including reduction of paddy terraces, reduction of outliers, incorporation of water bodies, and reduction of errors by error propagation. DTM was produced by interpolation from a 1 : 25 000 15 Topographical Map with contour interval 12.5 and elevation data from the DGPS profiling. DTM processing generated several morphometric and hydrological variables such as slope, plan curvature, SPI (Stream Power Index) and SCI (Shape Complexity Index) (Fig. 3). DTM-derived products were processed in ILWIS software with several available scripts in Hengl et al. (2009). 20 The other morphometric variables were rockfall velocity and energy. There were processed by RockFall Analyst as an extension of ArcGIS (Lan et al., 2007). It included modeling of rockfall trajectory by kinematic algorithm and raster neighbourhood analysis to determine velocity and energy of rockfall. Rockfall velocity and energy analysis needed information about slope geometry and other parameters such as mass, initial 25 velocity, coefficient of restitution (Table 1), friction angle and minimum velocity offset. Slope parameter was derived from corrected DTM. The other parameters were derived from secondary data and field data. For example, coefficient restitution and friction angle were derived from literature review based on landuse map and geological map, ESURFD 2, 2014 Automated landform classification in a rockfall prone area G. Samodra et al. whereas mass was determined from the dimension of boulders derived from field data measurement. The landform elements were derived, as the 9-unit slope model, by using the unsupervised fuzzy k-means classification (Burrough et al., 2000) as where µ is the membership of i th object to the cth cluster, d is the distance function which is used to measure the similarity or dissimilarity between two individual observations, q is the amount of fuzziness or overlap (q = 1.5). Unsupervised k-means classification was written and applied in ILWIS script with an additional class center for each morphometric variable (Table 2). Modified 9-unit slope model was applied by excluding 10 alluvial toe slope and seepage slope into a final landform classification. Observed volume of rockfall and cumulative distribution in four generic landforms were plotted in a log-log chart. Hungr et al. (1999) and Dussauge et al. (2003) investigated the frequency volume distribution of rockfall volume obey a negative power-law scaling. They found that the rockfall volumes follow a power law distribution with rel- 15 atively similar exponent value. The observed cumulative volume distribution was adjusted by a power law distribution as follows: where N R is the number of events greater than V R , V R is the rockfall volume and b is a constant parameter (cumulative power-law scaling exponent). Linear regression was 20 adopted to estimate b value in rockfall volume statistic in this paper.

Landform classification
Geomorphometry defined as a quantitative landform analysis (Pike et al., 2008) was initially applied for the assessment and mitigation of natural hazard (Pike, 1988). Dijke and Westen (1990), for example, introduced rockfall hazard assessment based on 5 geomorphological analysis. Later, Iwahashi et al. (2001) analyzed slope movements based on landform analysis. Both utilized DTMs derived from interpolation of 1 : 25 000 scale contour map to analyze geomorphological hazard. Nowadays, the interpolation of contour map is still powerful to create medium scale mapping when better resolution DTMs are not available. However, reduction of error in interpolation of contour map is 10 needed to obtain plausible geomorphological feature. Reduction of paddy terraces, reduction of outliers, incorporation of water bodies, and reduction of errors by error propagation were applied in this study to improve the performance of DTM. The result shows that paddy terraces still exist where the sampling of elevation data are absent. In addition, "flattening" topography can also be found on 15 slopes less than 2 %. Remaining paddy terraces mostly occur in the transportational middle slope and flattening phenomenon mostly occurs in the interfluves. Both errors influence the plausibility of slope, (Fig. 3a) but those do not much influence the final classification of landform elements.
The consideration to decide upon the number of final classification of landform ele-20 ments and morphometric variable to be employed in the automated landform classification is essential. Final classification of landform elements should represent appropriate semantic description related to rockfall processes. Modified 9-slope model (Dalrymple et al., 1968) i.e. interfluves, convex creep slope, fall face, transportational middle slope, colluvial foot slope, slope, and channel bed was used to represent conceptual entities 25 of rockfall deposition in each slope segment. Convex creeps slope represents a potential rockfall source. Selecting morphometric variables should also consider rockfall processes, besides morphology shapes of the landscape. It should reflect the movement and deposition of rockfall boulders. Morphometric variable usage is a priori knowledge of geomorphologist in term of utilization of spatial thinking toward the process of recognizing the 10 generic landform in relation to rockfall processes. The experience and former knowledge are involved during the selection of morphometric variables.
Morphometric variables derivation through DTM processing was divided into two parts, i.e. morphometric variables derived from RockFall analyst (velocity, energy) and from ILWIS script (slope, plan curvature, shape complexity index, stream power in-15 dex). Rockfall velocity and energy are second derivative of DTM (Lan et al., 2007). The first derivative DTM i.e. slope angle and aspect angle were employed to compute the rockfall trajectory. Then, rockfall trajectory was used to model the rockfall velocity and rockfall energy by using neighborhood and geostatistical analysis. The highest velocity occurs in the transportational middle slope. Velocity gradually increases in the fall 20 face and decreases in the colluvial footslope. Since the energy is also calculated from rockfall velocity, the spatial distribution pattern of energy is rather similar to rockfall velocity. Both velocity and energy of rockfall influence the area of fall face, transportational middle slope and colluvial footslope. The first change of a pixel into zero velocity and energy of its neighborhood operation is determined as the end of boulder movements. 25 It means that the rockfall boulders are deposited on this site.
Plan curvature and stream power index influence the pattern of the convex creep slope and the channel bed. It forms water divide and stream channel. Shape complexity index, sliced using an equal interval 25 m, was measured as the complexity of outline of 2-D object. It predominantly influences the spatial distribution of the interfluves which has low value around 1. Low value of shape complexity index represents how simple and compact a feature is. Its effect on the other landforms is not apparent because the value of the shape complexity index in the lower slope is relatively homogeneous i.e. 4-5.

5
Selecting morphometric variables is important in automated landform mapping. The generic landform result will depend on how well morphometric variables are selected to perform automated landform classification. It represents how well morphometric variable can describe the specific process working on a landform element. Its spatial dependency influences the application of automated landform mapping in different places 10 and different geomorphological process. Thus, understanding the specific genetic in an area will give benefit in selecting morphometric variables. Some morphometric variables such as velocity and energy will not be useful to be used in fluvial, marine and aeolian landform. In a fluvial landform, the scenario of flood inundation map may replace it. 15 The final classification result (Fig. 4) was draped over DTM. The volume statistic rockfall deposit was employed to evaluate the coincidence between landform classification and rockfall frequency-magnitude. Since landform classification considers surface form and process, we argue that landform classification in a rockfall prone area exhibit scale-specifity (Evans, 2003). The magnitude (volume) and frequency of boulder 20 deposits may have a specific scale related to each generic landform.

Volume statistic of rockfall based on landform
The 521 rockfall deposits obtained from a geomorphological inventory range in size from 18 × 10 −4 to 3.6 × 10 3 (Fig. 5). Volume statistic of rockfall was observed based on the main landforms corresponding to rockfall deposition, i.e. fall face, transportational 25 middle slope, colluvial foot slope and lower slope. The cumulative volume distribution of Fig. 6 and Table 3  Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | b = 0.58, 0.73, 0.68, 0.64 respectively. The power law distribution is well fitted to explain 55 %, 30 %, 36 % and 38 % boulder deposits population in the fall face, transportational middle slope, colluvial foot slope and lower slope respectively. The model is not well fitted to explain small size rockfall deposit due to rollover phenomenon. It needs many reports obtained from complete historical rockfall data in many places with different 5 physical characteristic to obtain "general" agreement in assessing rockfall distribution. In many places, such complete historical data are absent. However, several authors agreed that volume distribution of rockfall follows a power law distribution. However, there is still lack knowledge of the b value due to the absence of complete inventory data. Dussauge et al. (2003) argue that the variation of b 10 value is due to the variability of cliff dimension, area scale, lithology, geometrical and mechanical parameters of rockfall. Hungr et al. (1999) proposed that the jointed rock (b = 0.65-0.70) has higher b value than massive rock (b = 0.40-0.43). Gunung Kelir is subvertical cliff dominated by calcareous limestone. It has b = 0.58-0.73 for volume larger than 2 m 3 and 10 m 3 . It shows that this study may also confirm the b value for 15 jointed rock proposed by Hungr et al. (1999). Lithology and surface material play important role for rockfall volume distribution. It influences the bouncing velocity during the impact between rockfall boulder and surface material. Soft rock tends to reduce the energy and decrease the velocity of rockfall. Landform also influences the value of scaling exponent. Fall face has the smallest b 20 value compared to another landform. It also indicates that lower frequency of smaller event is more dominant in fall face. Whereas, higher frequency of greater event is more dominant in transportational middle slope and colluvial foot slope. It shows the gradation pattern of rockfall deposition around generic landforms. This may correspond to the morphometric condition. The shape and characteristic of surface i.e. morphometric 25 variables determine how a rockfall was deposited in a generic landform. Formerly, we considered that the distribution pattern along the x-axis and y-axis was influenced by the number of measurements in the rockfall boulder datasets in each landform. The distribution pattern of fall face seems similar to lower slope meanwhile transportational for fall face and lower slope; and < 80 m 3 for transportational middle slope and colluvial foot slope. As stated by Brunetti et al. (2008), we also consider that the distribution pattern is not influenced by the number of measurements in the dataset. All landforms exhibit a rollover of the frequencies rockfall boulders. It is similar to the 5 rollover identified by (Dussauge et al., 2003;Hungr et al., 1999;Guthrie et al., 2004). Roll over occurs in the volume size around 3 m 3 for fall face, 11 m 3 for lower slope, and 6 m 3 for colluvial foot slope and transportational middle slope. Since the rockfall process is more related to the deposition zone rather than failure zone, a "rollover" to frequencies should be addressed to the process during the impact between rockfall 10 boulder and surface. The rockfall boulders were deposited on the site when the local surface can decrease volume and energy into zero velocity. It can be influenced by soft surface condition and or obstacle which can interfere the movement of a boulder. The gradation pattern of rockfall deposition may be addressed to scale-specifity (Evans, 2003). The volume of the individual rockfall deposit in the fall face spans 5 15 orders magnitude. The landforms which have higher order magnitude are lower slope, colluvial slope and transportational middle slope respectively. Careful attention should be addressed to the maximum individual boulder deposited in the lower slope. Figure 3 shows that the volume rockfall deposit in the lower slope spans 6 orders magnitude. However, it indicates a long missing gap between the largest boulder and the second 20 largest boulder. The local surface parameter may influence this problem. We consider that maximum order magnitude of lower slope is rather similar with colluvial foot slope (around 400 m 3 ). The likelihood for the deposition of greater rockfall volume can be defined. The higher magnitude of rockfall is more likely to be deposited on the transportational middle slope rather than colluvial foot slope, transportational middle slope 2,2014 Automated landform classification in a rockfall prone area G. Samodra et al.

Implication for preliminary rockfall risk analysis
In the past, many people used to consider that the natural hazard should be approached on the domain of engineering. However, both structural and nonstructural mitigation should be included in natural hazard mitigation comprising geomorphological, geographical, and geological approach (Oya, 2001). Furthermore, natural hazard 5 could be found spatially through unique characteristic of an area. Specific geomorphology features may pose a specific hazard. Fall face, transportational middle slope, colluvial foot slope and lower slope exhibit scale specifity. The most susceptible places, in order, for rockfall hazard in Gunung Kelir area are fall face, transportational middle slope, colluvial footslope and lower slope respectively.

10
Automated landform analysis and rockfall statistic can pose the likelihood of rockfall magnitude in a specific landform. Each generic landform indicates the susceptibility degree to rockfall events. The magnitude-frequency relation of rockfall can be calculated to estimate the annual frequency of rockfall events in each generic landform. It can be defined with reference to specific event magnitude class in a specific generic 15 landform. Similar with Eq. (3), Dussauge et al. (2003) and Malamud et al. (2004) show that magnitude-frequency distribution of rockfall events in given volume class j followed power law distribution and can be described as: where N(V ) is the cumulative annual frequency of rockfall events exceeding a given 20 volume V , N 0 is the total annual number of rockfall events, and b is the power law exponent. Adopting a Poisson model for the temporal occurrence of rockfall, the probability of experiencing n rockfall during time t (adapted from landslide) is given by (Crovelli, 2000): where P [N(t) ≥ 1] is the probability of experiencing ≥ 1 rockfalls during time t, λ is the estimated average rate of occurrence of landslides which corresponds to 1/µ, with µ the estimated mean recurrence interval between successive failure events. The variable µ can be derived from the cumulative frequencies for each considered volume class N(V ) in Eq. (4) (Hungr et al., 1999). Thus, P (L) j km (Eq. 1), the temporal probabil-5 ity (exceedance probability) of rockfall in the magnitude scenario for a different period, can be obtained in a specific generic landform. Preliminary rockfall analysis can be delivered by evaluating elements at risk located in the susceptible place for rockfall hazard. There are 3 buildings located on the transportational middle slope and 14 buildings located on the colluvial footslope. This is 10 useful information on which to base prioritization action for countermeasures policy and design. Geomorphologic analysis should be taken into account to locate structural measures (e.g. barriers, embankments, rock sheds) in suitable location. It will improve cost efficiency to optimize budget and design. The information of building located on the landform classified as high hazard can also be an input to the prioritization 15 of an evacuation procedure. Therefore, the prioritization of mitigation action based on geomorphometric analysis can meet the technical suitability and the effectiveness of selected mitigation options.

Conclusions
Geomorphometry application can be an alternative tool to minimize the subjectivity of 20 Indonesia's standard landform classification applied in disaster risk reduction. Application of unsupervised landform classification in Gunung Kelir poses reasonable result for preliminary rockfall risk assessment. Rockfall protection through structural measures and landuse planning should take into account landform analysis.
However, the original classification of 9-unit slope model should be modified if it 25 is applied in different places. It should consider the genetic working on the specific landforms. slope, fall face, transportational middle slope, colluvial footslope, slope and channel bed (Fig. 2) is different with the original classification of the 9-unit slope model. Alluvial toe slope and seepage slope were excluded from the final landform classification in the present study. Channel wall was also modified as lower slope. Since the study area is located in the upper part of Kulon Progo Dome, the depositional process of alluvium 5 does not work in that such area. Seepage slope was merged with interfluves slope because interfluves slope and seepage slope landform classification is more related to pedogeomorphic process rather than gravitational process. The considerations to merge and exclude some landforms were based on the experience and the judgement of researchers. The proposed methodology applied in the rockfall prone area should be 10 tested in different areas which have similar genetic. Further studies should also explain the effects of scale and spatial dependency on the landform classification. 2,2014 Automated landform classification in a rockfall prone area G. Samodra et al.   Fall Face 53 18 × 10 −4 -1.0 × 10 2 2-1.0 × 10 2 0.58 0.98 Transportational Middle Slope 211 39 × 10 −4 -3.6 × 10 3 11-3.6 × 10 3 0.73 0.99 Colluvial Foot Slope 199 37 × 10 −4 -4.8 × 10 2 10.5-4.8 × 10 2 0.68 0.99 Lower Slope 58 21 × 10 −4 -3.6 × 10 3 11-3.6 × 10 3 0.64 0.97