Spatial-Temporal Dynamics Land Use/Land Cover Change and Flood Hazard Mapping in the Upstream Citarum Watershed, West Java, Indonesia

The input data for this study, comprising LULC maps information for 1990, 1996, 2000, 2003, 2009 and 2019, were collected from studies conducted by Yulianto et al. (2018, 2019). The LULC information was obtained from 30-m multi-temporal Landsat imagery, which was classified using the maximum likelihood approach. There are 7 (seven) classes of LULC information, which are given as follows. Class 1, Urban/built-up area; Class 2, Primary forest; Class 3, Secondary forest and mixed garden; Class 4, Plantation; Class 5, Wet agricultural land; Class 6, Dry land farming; and Class 7, Water body. The maps and classes are presented in Fig. 3 and Table 1, respectively.

Fig. 3

Interpretation related to forest classes, secondary forests and mixed gardens, plantations and several other classes on Landsat imagery can be performed based on the recognition of spatial character. The characteristics of these objects can be identified based on the interpretation elements such as colour, shape, size, pattern, texture, location and appearance associations of the object. The training area is needed in the classification process as a class identifier on the object. The training area can reflect the characteristics of each object to be classified. Furthermore, the maximum likelihood classification method has been used to classify LULC in the study area.

Specifically for urban/built-up area, primary forest, plantation, wet agriculture land, dry land farming and water body, classes can be distinguished based on the characteristics of colour, shape, size, pattern, texture, location and appearance of the object association. Meanwhile, the class of secondary forests and mixed gardens, based on their characteristics in Indonesia, tend to have similarities and mixed with other objects. Because of this mixing with several other objects, classification of the class becomes difficult. Thus, in this study for the class, secondary forest and mixed garden are combined into one class.

The urban area expansion rate (UAER) gives the average annual rate of growth of an urban or built-up area during a given period and can be formulated based on Equation (1) (Rimal et al. 2018, Rijal et al. 2018). (1)UAER=(UAt2−UAt1)(Tt2−Tt1)×100UAER = {{(UA_{t2} - UA_{t1} )} \over {(T_{t2} - T_{t1} )}} \times 100 where: –

UAER is the rate of urban area expansion (ha a⁻¹),

A landform, as a part of geomorphology, is the smallest unit of physical feature of the earth’s surface that is produced by natural processes. Plains, hills, mountains and valleys are examples of features of the earth’s surface that can be described as landforms. It is important to study and observe the characteristics of landforms and their relationship with human life in order to be aware of the potential hazards and disasters in a region. Landforms are characterised by differences in the structure and process of geomorphology, topography and constituent material (Strahler 1957, Verstappen 1983, Ruffell, McKinley 2014). In prior decades, the classification of landform properties has been carried out by calculating the geometry manually. Recently, however, an increase and expansion in computer technology has led to new spatial analysis methods, the development of algorithms and the ease of obtaining digital elevation data, all of which have contributed to earth-oriented geomorphometrics (Tagil, Jenness 2008).

The development of computer technology has led to an acceleration in the processing of landform classification, which, in the past, was performed based on manual interpretation. Semi-automated landform-based classification is one such method that has both simplified and hastened the process of mapping and classification of landforms that was previously performed manually (e.g. Horton 1945, Coates 1958, Evans 1972). A GIS enables detailed analysis of land surfaces to be integrated by remotely sensed data, specifically digital elevation models (DEMs) (Schillaci et al. 2015). There are several methods of spatial analysis; the development of algorithms that contribute to earth-oriented geomorphometry include the curvature-based approach, fuzzy landform elements, pattern recognition, relief segmentation, object-based, morphometric feature and others. The development of these methods has thus helped to speed up the landform classification process (e.g. Irvin et al. 1997, Ehsani, Guisan et al. 1999, Burrough et al. 2000, Wilson, Gallant 2000, Weiss 2005, Ramalingan et al. 2006, Quiel 2008, Drăguţ, Eisank 2012, Jasiewicz, Stepinski 2013).

Terrain surface classification (TSC) was used as the approach for the semi-automated landform classification in this study. TSC, as proposed by Iwahashi and Pike (2007), can be used to illustrate the geometric signature, comprising a set of measures that describe the shape of the topography by distinguishing geomorphologically different landscapes. Geometric signatures are designed to create an estimated terrain-unit map that includes slope gradient, surface texture and local convexity. In this study, the semi-automated extraction of landform elements that can be derived from SRTM30 DEM data was used to create a landform-based classification using the TSC approach. SRTM30 DEM was provided by the U.S. Geological Survey (USGS). System for Automated Geoscientific Analyses (SAGA) version 6.3 software was used to automatically identify the land elements from the TSC extraction (www.saga-gis.org). According to Olaya and Conrad (2009), SAGA is full-fledged GIS software whose features have some relation with geomorphometry. It also includes a large set of geoscientific algorithms that can serve as powerful tools in the analysis of DEMs.

Ho and Umitsu (2011) conducted a research related to micro-landform classification and flood hazard assessment and integrated the methods using remotely sensed data. The development of methods for the integration of the classification of micro-landforms and flood hazard zones was based on a geomorphological approach using SRTM30 DEM and Landsat data combined with field investigation. Flood hazard zonation was created based on categories from geomorphological features that grouped the average elevation of each landform into flood hazard classes. Furthermore, the results of their study to create flood hazard zonation were validated using field surveys, topographic characteristics and the history of flood inundation maps. This study has similarities in terms of the research topics that were examined by Ho and Umitsu (2011) related to flood hazard assessment and mapping. The difference with this study lies in the proposed use of the TSC approach to create micro-landforms as a basis for creating flood hazard assessment and mapping in the study area. The geomorphological characteristics in this case are represented by micro-landform units, whereas the flood characteristics and ground elevation level/value of each landform are used to analyse and classify the flood hazards in the study area. Other supporting data, such as flood height, were obtained from field surveys, whereas historical flood inundation data were procured from Sentinel-1 data in order to ascertain flood probability information.

Multi-temporal remotely sensed data for the period 2014–2018, Sentinel-1 C-band Synthetic Aperture Radar Ground Range Detected (SAR GRD), were used to create a flood probability map capable of describing the history of floods in the study area. According to Twele et al. (2016), Sentinel-1 data with Level-1 Ground Range Detected (GRD) products consist of focused SAR data that have been detected, multi-looked and projected to ground range using an earth ellipsoid model. Several data processing steps for Sentinel-1 SAR GRD products have been carried out in this study using SNAP version 6.0.0 developed by the European Space Agency (ESA), namely, (a) pre-processing, calibration; (b) pre-processing, speckle filtering; and (c) post-processing, geometric correction.

Flood inundations can be identified using the change detection approach based on before-flood and during-/after-flood information/data available for the period 2014–2018. Implementation of the minimum error auto-thresholding algorithm was developed by Kittler and Illingworth (1986) and was used to classify flooding and non-flooding in the study areas. Field observations and interviews with residents of the surrounding area were conducted to gather information regarding the historical depth of the flood and the impact of the flood. Furthermore, flood probability can be calculated based on Equation (2). (2)Prob_total=Σx=1yProb_valuex(i,j)yProb\_total = \Sigma _{x = 1}^y {{Prob\_value_{x(i,j)} } \over y} where: –

Prob_total is the total flood probability value,

During the 26-year period of LULC observation from 1990 to 2016 (Fig. 4), LULC was dominated by dry land farming, primary forest, secondary forest and mixed garden, urban/built-up areas, wet agricultural land and plantation. Increases in LULC during this period occurred in urban/built-up areas, dry land farming and plantation. Meanwhile, decreases in LULC were recorded in respect of primary forest, secondary forest and mixed garden and wet agricultural land. For dry land farming during the period 1990–2016, the estimated change in area was an increase of 35,716 ha or 15.63%, whereas for urban/built-up areas and plantation, the increases were 11,374 ha (4.98%) and 3,097 ha (1.36%), respectively. Meanwhile, the estimated decreases in secondary forest and mixed garden, primary forest and wet agricultural farming were in the order of 32,018 ha (14.01%), 12,704 ha (5.56%) and 5,466 ha (2.39%), respectively.

Fig. 4

Furthermore, the LULC distribution during the period 1990–2016 can be divided into eight sub-watersheds within the study area (Fig. 5), namely, Cihaur, Cikapundung, Cikeruh, Ciminyak, Cirasea, Cisangkuy, Citarik and Ciwidey. In order to analyse the role of LULC changes and their effects on flooding, the dynamics of changes in LULC calculations during 1990–2016 were detailed to the sub-watershed level in the study area. Fig. 5 shows the distribution of the LULC classes within each sub-watershed. Cihaur, Cikapundung and Cikeruh are located in the northern part of Bandung, whereas Ciminyak, Cirasea, Cisangkuy, Citarik and Ciwidey are all located in the southern part of Bandung. In general, the primary forest class, as a rain catchment area, is better distributed in the southern part of Bandung than the northern part. Meanwhile, the distribution of areas in the urban/built-up area class, which have the potential to reduce the recharge area size, is concentrated in the northern part of Bandung in comparison to the southern part. During the period 1990–2016, there was an increase in the dry land farming class, which was distributed across all sub-watersheds and dominated the other classes.

Fig. 5

In addition to the primary forest that acts as a rain catchment area, the growth of urban/built-up areas can also lead to a reduction in the size of the catchment areas and thus contribute to flooding in the study area. The size of the urban/built-up areas and the rate of expansion of these areas across the eight sub-watersheds in the study area are presented in Tables 2 and 3. On the basis of Tables 2 and 3 and Fig. 6, it can be seen that of the eight sub-watersheds, three feature a high level and high rate of urban/built-up area expansion, namely, the Cikapundung, Cihaur and Cikeruh sub-watersheds located in the northern part of Bandung. The rates of urban/built-up area expansion for these 3 sub-watersheds during the period 1990–2016 were 149.9, 118.3 and 64.5 ha a⁻¹, respectively. Meanwhile, the sub-watersheds that had the rates of urban/built-up area expansion for 1990–2016 at or below 35.0 ha a⁻¹ are located in southern Bandung.

Fig. 6

The result of the semi-automated landform classification using the TSC approach in the study area is presented in Fig. 7. According to Iwahashi and Pike (2007), the parameters in TSC describe the geometric signatures used to create terrain-unit maps, which are related to slope gradient, surface texture and local convexity. In general, the micro-landform classification in the study area can be classified into 16 units. Furthermore, the Bandung Basin is used as a unit area boundary to enable a focus on the classification of potential flood areas. This is presented in Fig. 8 and combined with the information obtained using the Topographic Position Index (TPI) approach by Yulianto et al. (2019). Meanwhile, the use of TSC for micro-landform characteristics can be described in relation to flood conditions based on their characteristics, which is presented in Table 4.

Fig. 7

Fig. 8

On the basis of Table 4, TSC as a micro-landform classification approach displays characteristics in relation to flood conditions. There are 10 classes of TSC micro-landforms in the Bandung Basin, 8 of which are micro-landform classes that have the potential to be flooded based on their characteristics. Historical flood inundation data produced using multi-temporal Sentinel-1 data between 2014 and 2018 were used to analyse the probability of flood inundation in the study area, which is presented in Fig. 9. Field observations and interviews with residents of the surrounding area were carried out to gather information regarding the historical depth of flooding and the impact of the floods. Furthermore, the results obtained from mapping the TSC as a micro-landform and the probability of inundation and flood depth taken from field observations have been integrated to produce maps and a flood hazard assessment of the study area, as presented in Fig. 10.

Fig. 9

Fig. 10

The spatial distribution and analysis of LULC changes during the period 1990–2016 have shown an increase in some LULC classes that have the potential to reduce recharge areas and also contribute to flooding in the study area. The expansion of the urban/built-up area class of 447.3 ha a⁻¹ during the period 1990–2016 has thus contributed to a reduction in recharge areas, especially in the northern part of Bandung (Fig. 6), where the greatest contributions came from the sub-watershed areas of Cikapundung, Cihaur and Cikeruh. This type of expansion also leads to an increase in run-off associated with the reduction of recharge areas in the study area. In addition, there has been a change in LULC class, notably a shift towards dryland farming, which has come to dominate the other LULC classes and can also contribute to flooding in the study area. Some of the various types of plants used in dryland farming, such as potatoes, onions and others, tend to provide minimal vegetation cover and a high potential for erosion during rain that can affect the sedimentation and siltation of rivers in the study area. Reduced river capacity as a result of sedimentation tends to lead to the overflow of rivers and flooding. Meanwhile, the trend for year-on-year decreases in the coverage of both the primary forest and secondary forest and mixed garden class of LULC in the area also impact on the reduction of infiltration zones, notably in the upstream areas of the study area. According to Rosyidie (2013), one cause of flooding in the study area was the change in LULC in those upstream areas that were previously forest and have since been converted into dry land farming, comprising horticultural farms with seasonal plants such as potatoes, carrots, mustard greens, cabbage and beans that tend to require a short time to harvest. High erosion rates occur because these seasonal plants do not protect the soil from erosion and contribute to sedimentation and flooding in the lower regions. Bosch and Hewlett (1982) described the strong correlation that exists between a reduction in forested areas, or deforestation, and an increase in flooding, whereas Wan and Yang (2007) explained how the expansion of urban/built-up areas is one of the major contributors to an increase in flooding. Furthermore, Warburton et al. (2012) outlined how the contribution of LULC to flooding depends on the proportional and spatial distribution of various types of LULC within the watershed location. The results of this study (Fig. 5) show that the proportional distributions are different for each type of LULC in the sub-watersheds in the study area. The sub-watersheds of Ciwidey, Cisangkuy, Ciminyak and Cirasea are areas with the highest proportions of primary forest compared to the other sub-watersheds, which indicates that these areas remain relatively strong in terms of their catchment areas. Meanwhile, the sub-watersheds of Cikapundung, Cihaur, Cikeruh and Citarik contain more limited proportions of primary forest. This, coupled with relatively large proportions of urban/built-up areas and much in the way of dryland farming in these areas, tends to contribute significantly to flooding in the study area.

The result of the flood hazard assessment and mapping is based on the integration of TSC as a micro-landform classification approach with the flood inundation probability map and depth of flood data obtained from field observations (Fig. 10). In this study, the flood conditions and their interaction with micro-landform features are influenced by the topographic characteristics of the region in the form of the basin area, referred to collectively as the Bandung Basin or Cekungan Bandung. According to Narulita et al. (2008), the boundary of the Bandung Basin is an area based on the distribution of ancient Bandung lake sediments that morphologically formed a plain, with the surrounding area as the source of the lake sediment. This area is an ancient Bandung sedimentary plain that dried up thousands of years ago. The flooding mechanism in the study area is based on the topographic characteristics affected by the former lakes, alluvial fan plain and major river (i.e. the Citarum). The alluvial fan plain spreads to cover the northern part, and the river flows across the Bandung lake plain in meandering patterns, especially in the south of Bandung. The results of the flood hazard assessment map (Fig. 10) in the very high flood hazard class show the micro-landform of the study area in the form of a plain area with a slope of less than 9%, or gentle slope, which has high convexity (concave; for hazard 1) and low convexity (convex; for hazard 2) and is also characterised by the fine texture of sediment, the dominant proportion of which comprises finer particles such as silt and clay. The flood depth in the area is close to exceeding 5 m. Fig. 11a and 11b shows examples of locations in the very high flood hazard class located in the Baleendah district. Fig. 11c shows a location in the high flood hazard class and illustrates a micro-landform condition in the form of a plain featuring a gentle slope with coarse texture, thus indicating a greater composition of sand (compared to fine texture, dominantly clay) and which also has both high convexity (concave) and low convexity (convex). At this location in the Bojongsoang district, the depth of flood is between 2 and 4 m. Meanwhile, an example of a moderate hazard area location is given in Fig. 11d, which is located in Cibiroso, featuring a micro-landform class of open slopes, moderate slope, coarse texture, high convexity and a flood depth of 1–2 m. Hannan (2017) produced a map of potential flood hazards in the same location (in the study area) using the Topographic Wetness Index (TWI) approach. SRTM30 DEM data were also used as an input to calculate the flow accumulation and slope, which were then used to describe the wetness of the soil. According to Huang et al. (2008), Grabs et al. (2009) and Yang et al. (2015), TWI has a positive correlation with the characteristics of soil moisture, such as horizon depth, organic matter and clay content. Furthermore, TWI has been widely used to determine the wet conditions in a watershed and is able to identify areas that are saturated with water and thus have high flood potential. The results of the study by Hannan (2017) show that the TWI value in the study area ranged from −0.03 to 7.85, which represents a large potential for flood disasters; this can then be classified as a map of potential flood hazards. This is in accordance with the results of this study, whereby the micro-landform areas with TSC classes of plain, gentle slope, fine texture and high convexity, and also those with low convexity, have very high flood hazard potential in accordance with the TWI value, which indicates a high wetness level. Meanwhile, regions that are characterised as being in the low flood hazard class have a TWI value indicating a low wetness level. Nurjaman (2018) used a GIS with the Complete Mapping Analysis (CMA) approach to identify the level of flood-prone parameters in the study area examined. The weighting method used to determine the score of each parameter was empirically based on field observations and the number of flood events. The parameters used were slope, altitude, rainfall, soil texture, distance from the river, population density and land cover. The results of the study revealed that some of the sub-districts in the class that is very prone to flooding included Dayeuh Kolot, Rancaekek and Margahayu. Some of the locations included in flood-prone classes included Arjasari, Paseh and Baleendah, whereas locations that were not prone to flooding included Pengalengan, Pasir Jambu and Kertasari. The results show that the same thing related to several locations with a flood-prone level, with the results of this study related to areas that range from the very high to low flood hazard classes in the study area.

Fig. 11

This study has analysed the description of LULC change information and its contribution to flooding in the study area. However, this study has not included a detailed quantitative analysis of the response of LULC changes to flooding. The next study will conduct an analysis regarding this matter in order to determine the significant influence and contribution of LULC changes to flooding. The integration of TSC as a micro-landform classification approach, the probability map for flood inundation and the depth of flood data from field observation can be used as an alternative to developing flood mapping and hazard assessment. It would also be possible to apply these methods to other regions that have limited hydrological data and other parameters for the purpose of creating a flood hazard map as a basis for hydrological and meteorological modelling. This type of combined production related to flood hazards using various geomorphological, hydrological and meteorological approaches, in addition to several other approaches, is expected to yield more detailed and complete information. This study uses LULC information that was produced with Landsat multi-temporal data, along with SRTM30 DEM that can support scale mapping up to 1:50,000. The use of high-resolution satellite imagery and also high-resolution of digital elevation models (DEMs) will certainly produce more detailed information with respect to producing LULC information and also for the use of TSC for a micro-landform area. In addition, this study has not included issues related to the environment, which also have an impact on the causes of floods in the study area. These issues include groundwater extraction and the phenomenon of land subsidence and deposition in the river sediment. As such, this constitutes an area for potential future research by adding parameters related to the issue in order to obtain more detailed and complete information.