Mapping Landform and Landslide Susceptibility Using Remote Sensing, GIS and Field Observation in the Southern Cross Road, Malang Regency, East Java, Indonesia

: There has been an increasing trend of land area being brought under human’s use over time. This situation has led the community to carry out land-use development activities in landslide hazard-prone areas. The use of land can have a positive impact by increasing economic conditions, but it can have negative impacts on the environment. Therefore, this study aimed to identify the landslide hazard, focusing on the development of a landform map to reduce the risk of landslide disaster in JLS, Malang Regency. The integration of remote sensing and geographic information systems, as well as ﬁeld observation, were used to create a landform map and a landslide susceptibility map. Using the geomorphological approach as a basic concept in landform mapping, the morphology, morphogenesis, and morphoarrangement conditions were obtained from the remote sensing data, GIS, and ﬁeld observation, while morphochronological information was obtained from a geological map. The landslide susceptibility map was prepared using 11 landslide conditioning factors by employing the index of entropy method. Thirty-nine landform units were successfully mapped into four landslide susceptibility classes. The results showed that the study area is dominated by a high level of landslide susceptibility with a majority of moderate to strongly eroded hill morphology. It also reafﬁrms that landform mapping is a reliable method by which to investigate landslide susceptibility in JLS, Malang Regency.


Introduction
Indonesia is a country that has a high risk of landslides [1][2][3]. Of the 9383 landslide disasters that have occurred in Indonesia, 1483 were in a single year in 2019 [4]. They are spread across almost all provinces in Indonesia, including the East Java Province in general and the Malang Regency specifically [5]. More than 80% of the southern region in the Malang Regency is categorized as being at a high risk of landslide occurrence [6]. Landslides can have an impact on environmental damage and cause both physical and economic losses [7,8].
The occurrences of a landsalide is triggered by many factors, such as topography, climate, vegetation, land use, earthquakes, and others [2,9]. Factors affecting landslides can be divided into intrinsic factors and extrinsic factors. An intrinsic factor is the main factor originating from the condition of the land itself, while an extrinsic factor is a trigger factor from outside which can increase the potential for landslides [10]. Parameters included in intrinsic factors are topographic conditions, soil material, and geology. Increased steep slope conditions can lead to the more frequent occurrence of landslides due to low soil stability [11]. Meanwhile, extrinsic parameters are the result of factors such as human activity, e.g. the development of residential areas on steep slopes and the construction of roads with slope cuts and improper slope loading, rain, and earthquake activity [3].
There has been an increasing trend of land areas being brought under human's use over time. This situation has led the community to carry out land-use development activities in landslide hazard-prone areas, such as hilly areas and steep slopes. The southern part of the Malang Regency of Indonesia is dominated by hills with fairly steep slopes, where the Southern Cross Road (JLS) has been built with the aim of increasing the accessibility of people between regions, especially for tourism promotion. With this increased accessibility, it is likely that the use of new areas in hilly terrain will increase too. However, the region's morphological condition, dominated by steep slopes, means that unless development is planned properly the area possesses a very high risk of landslide occurrences triggered by extrinsic factors. In order to reduce the risk, a comprehensive understanding, in particular of the physical aspect of the area, is urgently needed.
Traditionally, one of the first steps in landslide disaster management is to create a landslide susceptibility map [12] through which an estimate of the level of potential landslides can be investigated [13]. Thus, a landslide susceptibility map serves as a reference for disaster management [14]. Information on landform conditions are very useful for preparing the landslide susceptibility maps [15] by identifying, mapping, and analyzing the landslide risk areas [16][17][18]. Since the study of landform involves a detailed investigation of the morphology, flow patterns, and processes related to landslide susceptibility, a landslide susceptibility map is useful in determining landslide-prone areas [19]. The basic concept used in landform study consists of understanding four elements of analysis, namely: morphology, morphochronology, morphogenesis, and morphoarrangement [20,21]. In conducting a landform analysis morphological identification has an important role, because for each individual slope the process that occurs is different, which in turn produces different materials [22]. In addition, these differences of process and material will have different effects on the frequency of landslide occurrences.
Several studies have been developed and carried out to create a landslide susceptibility and hazards assessment [23]. Scholars such as Gökceoglu and Aksoy have proposed the landslide susceptibility mapping using deterministic stability analyses and image processing technique for use in the Mengen region, Turkey [24]. This approach studied the mechanisms of the landslides by two-dimensional stability analysis from field observation and the parameters controlling such slide development. Furthermore, Clerici et al. applied the procedure for landslide susceptibility zonation by the conditional analysis method, which is conducted in the Parma river basin, Italian Northern Apennines [25]. This approach uses a multivariate method that simultaneously considers all the factors contributing to instability by means of a geographic information system (GIS). Landslide susceptibility is determined by calculating landslide density according to the combination of different instability factors. Another scholar used the artificial neural network for determination and application of the weights for landslide susceptibility mapping, which is applied in Yongin, Korea, using a GIS as the basic analysis tool for spatial data management and manipulation [26]. Landslide location, slope, curvature, soil texture, soil drainage, effective thickness, wood type, and wood diameter were used for analyzing landslide susceptibility. In addition, Rawat and Joshi proposed the landslide susceptibility index (LSI) method for use in the Igo river basin, Eastern Himalaya, India [27], and Spinetti et al. proposed the use of the LSI method in the Sorrentina Peninsula, Southern Italy [28]. However, using a landform map as a basis for the creation of a landslide susceptibility map is rare, particularly in Indonesia.
In order to create a landform map, the integration of remote sensing and a GIS plays an important role [29][30][31][32][33]. While the remote sensing data can be used to produce land cover information, the digital elevation model (DEM) can be used for hill shade and slope analysis, which is then reduced to a topographic map to obtain detailed information about the morphological conditions of the study area [21]. The GIS has a role in the presentation and storage of spatial data, field observation, and others data in the form of attributes, therefore it can store the information needed in making landform maps [20]. Since the identification of landslide hazards based on a landform map is rare in Indonesia, this study aimed to identify landslide susceptibility through the development of a landform map, in order that a land susceptibility map can be an important input in the planning of land use development, to reduce landslide disaster risk in JLS, Malang Regency.

Methodology
This study used a GIS and remote sensing technology to create both landform and landslide susceptibility maps. The formulation of landform maps is based on: (1) morphological aspects for slope information, (2) morphoprocess for erosion and deposition information, (3) morphochronology for rock lithology information, and (4) morphoarrangement for the detailed slope information. The compilation of landform maps is used as one of the landslide conditioning factors and it then models the distribution of landslide susceptibility. The model used is the bivariate statistical analysis method, which is the index of entropy. The model was applied and developed to analyze the spatial relationship between landslide conditioning factors and the distribution of landslides in the study area. The procedures followed for landform and landslide susceptibility mapping in this study are presented in Figure 1.

Study Area
This research was conducted in JLS, Malang Regency, East Java of Indonesia, located between 8.269 • to 8.447 • S latitudes and 112.362 • to 112.785 • E longitudes with a total area of 437.95 km 2 , and elevation ranging from 14 to 748 masl. The study area included five districts, namely Donomulyo, Bantur, Gedangan, Sumbermanjing Wetan, and Dampit. The physiographic conditions of the study area are very diverse and consist dominantly of bumpy and hilly morphological conditions prone to landslide risk [34], and hence they serve as good site for this study ( Figure 2).

Data Availability
The remote sensing data used in this study were Sentinel-2B (10 m resolution) and ALOS PALSAR DEM (12.5 m resolution) acquired from the Alaska Satellite Facility website (https:/asf.alaska.edu). The Sentinel-2B that has been launched by the European Space Agency (ESA) is utilized for monitoring land use/land cover change [35]. The product characteristics are provided in Table 1 [36]. The data acquired on 22 April 2020 was processed to level-1C through cloud masking and geometric corrections. ALOS satellite, which is equipped with a PALSAR radar sensor, is useful for monitoring the landslide susceptibility [37,38]. We used ALOS PALSAR DEM for topographic and morphometric analysis as an input for determining landslide susceptibility ( Table 2). In addition, the information on lithology at the study area was obtained from the geological map of the Turen Sheet and Blitar in 1992, with the scale of 1:100,000 from the Center for Geological Research and Development. Table 2 describes the data characteristics and their sources from which various intrinsic and extrinsic factors were derived for creating landform and landslide susceptibility maps.

Data Analysis
In this study, landform and landslide susceptibility maps were created using the On-Screen Image Interpretation (OSII) method and the index of entropy model as suggested by [43].

Data Extraction Using DEM
The elevation data was derived through the neighborhood technique to determine topographic and morphometric features, such as slope, topographic position index, stream power index, and hill shading effects [44]. This technique can produce the detailed topographic position index (TPI) attributes. TPI is an index value that determines the difference between the elevation of the central point z 0 in a grid (DEM data pixel) and the average elevation of z, and can be calculated using the following Equations [45].
Positive values on the TPI attribute indicate that the central point is in a higher position than the average elevation conditions in the vicinity. The radius (R) value for the neighborhood technique can be set with the provisions, the higher value of R will produce more general landform units, while the small value is indicated by the presence of minor landscapes, such as valleys and ridges. The TPI attribute is useful to determine the dominant morphology. The combination of the radius size of the neighborhood technique is used to extract complex morphological information [45,46]. Furthermore, the morphoprocess identification for landform mapping requires morphometric attributes, which is derived from DEM data. We used the stream power index (SPI) attribute for this process due to its ability to detect the potential of erosion processes. Greater SPI values indicate the erosion hazards due to a higher amount of overland water runoff from the upper slopes [44]. The SPI product was created by using the Raster calculation tool available in the ArcGIS Pro v2.5 software. SPI can be calculated as shown in the following Equation.
where, SPI is the stream power index. A s is the specific catchment area (m 2 /m) and β is the slope of the degree unit.

On-Screen Image Interpretation (OSII)
The OSII method is used to generate landform aspects, such as morphological data (slope and hill shade), morphochronology (geological maps), morphoprocess (stream power index), and morphoarrangement (topographic position index). The hybrid approach was implemented in making landform maps by superimposing all landscape features to identify unique landforms [43]. The reason behind using this method is to follow the mapping scales of the maps that were used in the analysis, including the spatial resolution of remote sensing data. The scale of the map can be determined through mathematical calculations based on the selection of optimal spatial resolution [47]. The mapping scale of 1: 50,000 requires at least remote sensing data with a spatial resolution of 10-20 m.

Landslide Susceptibility Mapping
Landslide events in the study area represent the interaction of several conditioning factors, such as slope aspect (slope orientation), compound topographic index, elevation (m), landform unit, land use, normalized difference vegetation index (NDVI), plan curvature (100/m), profile curvature (100/m), slope ( • ), stream density (km/km 2 ), and stream distance (m) [48,49] as presented in Figure 3. There is no universal standard or rules to select the conditioning factors, rather the selection is based on the condition of the study area itself because of these are specific to the different areas [50]. All conditioning factors were prepared from several sources as presented above in Table 2. ALOS PALSAR DEM data was used to extract the topographic and morphometric attributes, such as aspect, elevation, plan curvature, profile curvature, and slope. Environmental factors, such as NDVI and land use are obtained by classifying the sentinel-2B remote sensing data, and hydrologic factors, which include stream density and stream distance were prepared by processing vector data of river networks into raster data obtained from the Geospatial Information Agency [42].
The mapping of landslide susceptibility provides the basis for predicting landslide events in the form of landslide susceptibility zoning. There are several methods and techniques for landslide susceptibility zoning, such as qualitative vs quantitative, and direct vs indirect [51]. This study used the quantitative method based on bivariate statistics for the prediction of landslide susceptibility in JLS, Malang Regency. The index of entropy model used in deriving land susceptibility accommodates the calculation and measurement of instability, imbalance, interference, and uncertainty within a system. After that, the system accepts entropy values having a one-to-one relationship based on Boltzmann's principle of the level of interference of the system itself [52]. The Shannon-refined Boltzmann's principle specifically uses an entropy model related to information theory.
The index of entropy model has been widely used to determine the natural hazard weighting index and is also implemented for environmental modeling, such as the prediction of landslide susceptibility, drought level index, and groundwater quality assessment [53][54][55]. Entropy information on landslides refers to the degree of influence on various landslide conditioning factors. Additional value is given to several important factors into the index system, and hence the entropy value can be converted to an objectively weighted value in the index system. Calculation of entropy values for the prediction of landslide susceptibility were done using the following Equations [53,56].
H jmax = log 2 S j , S j = number o f classes (7) where, a and b are the percentage of landslide events and domains (classes of each variable), respectively, P ij is the probability density, H j and H jmax are entropy values, I j is the coefficient of information, and W j is the weight value generated for the variable as a whole. The landslide susceptibility index is calculated by adding the weight value of the variable, which has been reclassified for the second time. The values of landslide susceptibility index using the index of entropy model is presented below Equation ( where, y is the landslide susceptibility index value, aspect rec is the second reclassification result and so on the next variable, W j is the weight value for each variable that has been calculated based on the entropy value information. The results of the susceptibility index have a continuous interval data type. The division of susceptibility index classes refers to the natural-breaks method of classification [53,57,58].

Results and Discussion
The data analysis produced 39 landform units in the study area. Compared to other landforms in the East Java province of the country, the landforms in the study area were found to be more diverse in characteristic, as the study area is dominated by Karst landforms. Figure 4 presents the spatial distribution of the landform factors and landform units and Table 3 presents the area coverage of each landform unit.

Morphological and Morphoarrangement Conditions in the Study Area
The morphological conditions play an important role in determining landform boundaries, because the morphological configuration is able to present different processes. As described in the above section, the information on morphological conditions (hill shade and slope) were extracted from ALOS PALSAR DEM data, and the morphoarrangement condition from the TPI. Thus, it was possible to determine the boundary difference of morphology and morphoarrangement. The morphological conditions in this study were divided into three classes: hilly, moderately hilly, and flat, while the morphoarrangement was arranged into five classes, namely: peak of hills, upper slopes of hills, middle slopes of hills, lower slopes of hills, and colluvium and alluvium foot slopes.

Morphochronological Condition in the Study Area
The landform condition, especially the morphochronological condition, is closely related to the lithology conditions. In this study, information on the lithology conditions was done through manual delineation using the geological map, which showed six lithological units in JLS, Malang Regency (Table 4). In the study area, Karst materials are the dominant lithological units and most rock formations are of tertiary Miocene age. It affects the level of rock resistance to erosion, and the rock resistance is an important parameter for determining landslide susceptibility [59].

Morphogenesis Condition
Morphogenesis condition represents the dominant process that occurred in the study area. This information was obtained from the SPI OSII attribute. The index value of SPI is related to the strength of water flow, which eventually depends on the slope and catchment area. The greater the index value, the greater the potential for strong erosion, and vice versa. Although, different processes produce different landform conditions, the process that occurs is usually closely related to lithology and the slope gradient. Areas with steep slopes are associated with the erosion process, while gentle slope areas are associated with the deposition process. Rocks that have a high level of resistance will have a low erosion tendency [60]. The rock resistance to erosion can be known from the flow density. Greater flow density indicates that rock conditions are prone to erosion [59].
The results of morphoprocess identification at the study site revealed four erosion categories, namely significantly eroded, moderately eroded, slightly eroded, and deposited. The study area is dominated by significantly eroded areas (62.52%), and moderately (30.36%) eroded areas as can be explained by the hilly and undulating topographic characteristics. The process of strong erosion occurs mostly in the upper steep slopes of the hills. Meanwhile, the deposition process mostly occurs on the foot slopes because the slope is gentler and flat, covering a small proportion of the study area.

Analysis of Landslide Susceptibility
The zoning of landslide susceptibility was done based on the previous landslide inventory data, as these data are very useful for reconstructing landslide activity, building databases, and assessing landslide susceptibility [61]. The prediction of landslide susceptibility also requires historical landslide data using entropy models and is calculated based on the value of the probability density P ij on the occurrence of landslides in each variable class.
The results of the calculation of the weight valuesn Table 5 show the importance of each slope instability class. It is indicate that the higher of weight value, the higher contribution of landslide occurrences. The largest weight value W j is in the variable landform unit (0.003061), followed by slope (0.00164) and land use (0.000918). A class that has no zero occurrence indicates no effect on the landslide itself. The small weight value indicates that the proportion between the number of landslide pixels and each class of variables is small. The most important landslide conditioning factors for modeling landslide susceptibility, especially using the index of entropy model, in this study are landform units, slope, and land use. Topographic conditions with steep slopes have the potential for very high landslide susceptibility [3]. From the data analysis, the resulting landslide susceptibility in this study was divided into four categories; low, medium, high, and very high. It was found that the high susceptibility class has the largest proportion covering 51.75% of the total study area. About 13.59% of the area falls into the very high susceptibility class and about one-third of the area falls combinedly in the moderate and low susceptibility classes. Spatially, the areas that are very prone to landslides are widely spread in the eastern part of JLS ( Figure 5), mainly due to the morphological factors dominated by hills with a steep slope.
In terms of lithology, the study area is dominated by limestone composed of a variety of carbonate minerals. The mineral has high dissolving characteristics and thus makes the pores of the constituent rocks easy to crack. The continuing process of Karstification caused the limestone resistance to become weaker [62], resulting in the increased occurrences of landslides [63]. In addition, the result of the study showed that landslides often occur in the corridor of a road, where translational landslides are common due to the weakening of rock resistance as the road construction cuts across this landform ( Figure 6). The pictures, as selected example sites in the study area, depict most of the landslide occurrences due to a cut-slope failure during the road construction, posing a highly hazardous zone for landslides.
Human activity, such as deforestation and land conversion, are extrinsic factors in high risk areas that are frequently the reason for the occurrence of landslides. Another significant reason for the occurrence of landslides is the expansion of settlements in the study area not only near the river but also in the hilly areas with high susceptibility. The process of urbanization and road construction activities have changed the topographic structure, causing slope deformation and instability [64] making the Malang Regency, especially along JLS, highly susceptible to landslides. However, with better-planned comprehensive JLS Malang Regency infrastructure development, it is expected that the flow of logistics distribution, improvement of accessibility, agropolitan production centers, and industrial zones can develop properly with less disturbance in high-risk areas [65]. In addition, the tourism sector in the Malang Regency, where there is a total of 44 beaches scattered along the south coast, is also very dependent on the development of this access [66,67].
The analysis of landslide susceptibility is expected to help minimize the impact and improve the accessibility of the community and the development of tourist areas by providing information to avoid the high-risk areas, as the consideration of disaster management is crucial.

Conclusions
Landform mapping in the JLS Malang Regency was carried out based on an understanding of the concepts of morphology, morphochronology, morphogenesis, and morphoarrangement. Geospatial technology in the form of remote sensing and GISs were used to conduct all analyses. The hybrid approach and OSII methods were used as they are able to produce more realistic landform units by being able to integrate the data and maps of various scales and resolutions. A zoning map of landslide susceptibility used bivariate statistical methods, namely the index of entropy model.
The most influential factor in the occurrence of landslides in the study area is the landform unit, slope, and land use. Thirty-nine landform units were successfully mapped with four landslide susceptibility classes. The study area is dominated by a high level of landslide susceptibility with a majority of moderate to strongly eroded hill morphology. However, the landslide conditioning factor will be different and vary depending on the geographical conditions of the region itself, hence the weight values in this study may not be suitable for other regions although general methodology can be adopted for landslide susceptibility mapping [68]. The results of the landslide susceptibility assessment model index of entropy can be used as a reference to minimize disaster risk and optimize regional development along the southern cross road (JLS) in South Malang.