Projected Curve Number (CN) Changes and Surface Runoff Response to Corrected Land Cover Dynamics in the Citarum-Majalaya Catchment

2.1 Study area

The study location is a catchment area of the Citarum-Majalaya stream gauging station (around 204.52 km²) and is part of the upstream Citarum sub-watershed (around 269.43 km²). The catchment is geographically located between 7°3'1.96"-7°14'33.1"S and 107°37'50.7"-107°48'40.8"E, and administratively is part of the Bandung District, West Java Province, Indonesia. The distribution of the catchment area is in 6 subdistricts, namely the majority are in Kertasari Subdistrict (33.5%), Pacet Subdistrict (30.8%), Ibun Subdistrict (26.6%), and a small portion is in Majalaya Subdistrict, Pangalengan Subdistrict, and Paseh Subdistrict. The location of the catchment in Bandung District can be seen in Figure 1.

image005.png

Figure 1. The study area and the rainfall station distribution

2.2 Land cover

The land cover map of the study location was created and corrected using the ArcGIS 10.7.1 software application for 2000 to 2021. The basic data in shapefile format from 2000-2021 was obtained from the Ministry of Environment and Forestry of the Republic of Indonesia. Land cover map data is corrected before being used further for land cover change analysis. Corrections are carried out for validation and to avoid errors in interpretation, so that land cover data is obtained with good accuracy. The land cover maps of the study area before correction, for the years 2000 and 2021, are presented in Figure 2.

ping_mu_jie_tu_2025-05-20_101108.png

Figure 2. The land cover maps of the study area before correction

Land use/cover and conservation structures (water and soil) greatly influence surface runoff [22], and in this research will be more on discussing the effect of land cover changes on runoff. Analysis of land cover changes in this location was carried out for 22 years (2000-2021). The land cover of this location was classified based on SNI 7645-2010 concerning land cover classification [23], namely forest, shrubs, plantation, settlements (built environment), open land, water bodies, dryland farming, and ricefields. Land cover influences identifying flood-prone zones and directly impacts infiltration rates [24]. Impervious surfaces in the settlement (built environment) can increase runoff and reduce water infiltration [25]. This analysis is conducted to obtain information on the pattern of change for each type of land cover, which becomes the basis for projecting changes and estimating CN values.

This study will also analyze the variation in land cover between 2000-2021, which is divided into several phases, namely 2000, 2006, 2012, 2018, and 2021. Variations in land cover changes are quantified using the dynamic degree of single land cover k and the dynamic degree of comprehensive land use ${{L}_{c}}$. The K reflects the level of annual variation in the area of a type of land between the initial and final phases, and ${{L}_{c}}$ reflects the comprehensive annual variation level of the area of all land types. The positive K and ${{L}_{c}}$ value indicate that the area of land cover type increases, and if the value is negative, then the area of land cover type decreases. The greater the absolute value of K and ${{L}_{c}}$, the more significant the variation. The K and ${{L}_{c}}$ value can be calculated using the Turner method [4, 26], with the formula as in Eqs. (1) and (2).

$K=\frac{{{U}_{b}}-{{U}_{a}}}{{{U}_{a}}}\times \frac{1}{T}\times 100%~$                    (1)

${{L}_{c}}=\frac{\mathop{\sum }_{i=1}^{n}\left| {{U}_{bi}}-{{U}_{ai}} \right|}{2\mathop{\sum }_{i=1}^{n}{{U}_{ai}}}\times \frac{1}{T}\times 100%$                           (2)

where, U_a and U_b are the area of land type in the initial and final phases; T is the length of the sample period (unit: year); n is the number of land types in the study area.

2.3 Rainfall

Rainfall analysis in this study utilizes precipitation data collected from three rainfall stations located around the study area, namely the Cibeureum Station (area of influence: 87.8 km²), the Cipaku-Paseh Station (area of influence: 84.4 km²), and the Situ Cisanti Station (area of influence: 31.6 km²). The spatial distribution of these stations, the availability of data, and their respective areas of influence are illustrated in Figure 1 and summarized in Table 1. The rainfall data were obtained from the Hydrology and Water Environment Center, Directorate General of Water Resources, Ministry of Public Works and Housing, Republic of Indonesia. Due to the spatial and temporal variability of rainfall, it is necessary to express the data in the form of regional rainfall, which in this study is calculated using the Thiessen polygon method (Figure 3). For analytical purposes, the regional daily rainfall data were subsequently converted into monthly values, with the average rainfall for the period 2000–2021 presented in Figure 4.

Table 1. Rainfall data availability

image008.jpg

Figure 3. The rainfall station and the Thiessen Polygon

image009.png

Figure 4. The monthly average regional rainfall

2.4 Hydrologic soil group

The Hydrologic Soil Group (HSG) map of the study area is derived from the broader HSG map of the Citarum Watershed. The map data, provided in shapefile format, were obtained from the Oak Ridge National Laboratory Distributed Active Archive Center (ORNL DAAC), NASA [27]. The HSG map specific to the study area (Figure 5) was generated through a Geographic Information System (GIS) process, involving an intersection between the study area boundary (Citarum-Majalaya catchment) and the HSG map of the Citarum Watershed. Subsequently, an overlay analysis was performed between the HSG map and the land cover map of the study area to generate Hydrologic Response Unit (HRU) data. These HRU data serve as the basis for calculating the composite Curve Number (CN) values within the study area.

image010.jpg

Figure 5. The HSG map of the study location

HSG classification based on soil texture class, along with runoff potential, is presented in Table 2 [27]. The four standard classes, namely A, B, C, and D, correspond to soils with low, moderately low, moderately high, and high runoff potential, respectively. Wet soils have high runoff potential (regardless of texture) due to the presence of a groundwater table within 60 cm of the surface [27]. Wet soils have a dual HSG (e.g., HSG A/D): they can act as HSG A or HSG D depending on drainage conditions or groundwater levels. This usually occurs when the surface soil has high infiltration characteristics (HSG A), but the subsoil has very low permeability (HSG D), which can cause water to pool or drain slowly.

Table 2. HSG classification based on soil texture class

Source: Study [27].

2.5 Methods

This research calculates the CN value and its projections using two approaches. The first is the assumption that there is no effort (do nothing) related to land cover control and planning, and the second is the land cover change scenario based on the spatial pattern plan. The projected CN value for the first approach (do nothing) uses a trend model from the observed CN value and uses a scenario land cover change trend model to calculate the scenario CN.

The novelty of this research is the analysis and projection of CN values at the study location using corrected land cover and land cover scenarios, as well as calculating surface runoff predictions using forecast rainfall data using the ANN method. The output of this research is surface runoff information based on scenario CN values and CN do nothing values. The output from these two approaches is then compared and analyzed. The results of this research can be used as recommendations for regional officials in formulating policies related to the preparation of spatial plans and natural resource management in their region. The research stage diagram is presented in Figure 6.

ping_mu_jie_tu_2025-05-20_101851.png

Figure 6. Flow chart of research stages

2.5.1 Analysis of changes in curve number 

Curve number (CN) value describes the soil's ability to infiltrate and turn excess rain into surface runoff. The CN value of a watershed or catchment depends on the soil type and land cover. For catchments that consist of several soil types and land covers, composite CN can be used to estimate direct runoff using the following calculation method [28, 29].

$CN=\frac{\sum {{A}_{i}}C{{N}_{i}}}{\sum {{A}_{i}}}$                       (3)

where, CN is the composite CN value, ${{A}_{i}}$ is the subcatchment area i, and $C{{N}_{i}}$ is the CN value of catchment area i. CN values range from 1 (minimum) to 100 (maximum) depending on the hydrological soil group HSG), land use/cover type, and previous soil moisture conditions [30]. Changes in the CN value of the study location were analysed from 2000 to 2021. Examples of overlay maps of land cover in 2021 and HSG are presented in Figure 7. The CN values for several HRU conditions are presented in Table 3 [31-33].

image012.jpg

Figure 7. The HRU maps

Table 3. CN value for the hydrologic response unit

2.5.2 Projection of CN values

The analysis of changes in CN values in this research uses the trend method. This method is an analysis technique used to identify data patterns within a certain period and can be used to make predictions for the future. The model generally consists of a linear trend and a nonlinear trend. One of the advantages of the non-linear trend method is that it is excellent for predicting long-term data, and the forecast results are close to actual values [34]. The trend pattern of the CN value data series was detected over 22 years (2000-2021). The trend model for CN values can be expressed in a general equation as in Eq. (4), where t is time (year) and ${{a}_{i}}$ are model coefficients [35]. If n=1, then the trend has a linear pattern; if n=2 and n=3, then it has a quadratic and cubic polynomial pattern, respectively.

$CN=\underset{i=0}{\overset{n}{\mathop \sum }}\,{{a}_{i}}{{t}^{i}}$                     (4)

$R=\frac{n\sum xy-\sum x\sum y}{\sqrt{\left( n\sum {{x}^{2}}-{{\left( \sum x \right)}^{2}} \right)\left( n\sum {{y}^{2}}-{{\left( \sum y \right)}^{2}} \right)}}$                        (5)

The first step in trend analysis is to visualize data by making graphs or plots of data. This is done to detect data patterns over time as a basis for selecting an appropriate trend model. The best (most appropriate) trend model is determined based on the R-squared (R²) statistical criteria value. The best model is the one with R² closest to 1. Then, based on the selected trend model, the CN value for the future is projected. The R calculation method is presented in Eq. (5). The R value measures the strength of the relationship between the projected CN and the actual CN, and R² measures how well the model fits the actual CN data.

The projection of CN values in this research is based on two approaches. The first approach is the assumption that there will be no efforts regarding land cover change, and the second is a scenario based on spatial planning. The second approach is created based on spatial pattern plans by the Regional Regulation of Bandung Regency No. 1 of 2024 concerning Regional Spatial Planning of Bandung Regency 2024-2044. Maps and data of land cover plans are presented in Figure 8 and Table 4. Projections of CN values using the second approach were carried out to understand the impact of land cover changes on surface runoff.

image013.png

Figure 8. Spatial pattern plan of Bandung Regency based on the spatial plan 2024-2044

Table 4. Land cover of Bandung Regency according to the spatial pattern plan in 2044

Evaluation of model performance is carried out to determine the model's reliability by taking into account errors from predicted data to actual data. Model performance is measured using the statistical value criteria. These values are Mean Absolute Percentage Error (MAPE), Nash-Sutcliffe Efficiency (NSE), and Root Mean Square Error (RMSE). The MAPE value can be calculated using Eq. (6).

$MAPE=\underset{i=1}{\overset{n}{\mathop \sum }}\,\left| \frac{{{y}_{i}}-{{{\hat{y}}}_{i}}}{{{y}_{i}}} \right|\times ~100%$                  (6)

RMSE measures the standard deviation of a model's prediction error; a smaller value means the model performs better [36]. The RMSE value is calculated using Eq. (7). NSE is used as an index of suitability between calculated and observed data. The model fits very well if the NSE value is equal to 1, and it is poor if it approaches $-\infty$ [37]. The NSE value is calculated using Eq. (8).

$RMSE=\sqrt{\frac{1}{n}\underset{i=1}{\overset{n}{\mathop \sum }}\,{{\left( \frac{{{y}_{i}}-{{{\hat{y}}}_{i}}}{{{Y}_{i}}} \right)}^{2}}}$                 (7)

$NSE=1-\frac{\mathop{\sum }_{i=1}^{n}{{\left( {{Y}_{i}}-{{{\hat{y}}}_{i}} \right)}^{2}}}{\mathop{\sum }_{i=1}^{n}{{\left( {{Y}_{i}}-\bar{y} \right)}^{2}}}$                  (8)

where, ${{Y}_{i}}$ = observation value; ${{\hat{y}}_{i}}$= model value.

The statistical test used in validating the CN value of the projection results is the student's t distribution test (t-test) to test the suitability of the model [26, 38], and using the MAPE statistical value criteria (4) in testing the model performance. Some of the t-test formulas used are shown in Eqs. (9) and (10). A t-test is also employed to evaluate the differences in projected Curve Number (CN) values obtained from two distinct approaches.

$t=\frac{\left| {{{\bar{X}}}_{1}}-{{{\bar{X}}}_{2}} \right|}{\sigma \sqrt{\frac{1}{{{N}_{1}}}+\frac{1}{{{N}_{2}}}}}$                 (9)

$\sigma =\sqrt{\frac{{{N}_{1}}S_{1}^{2}+{{N}_{2}}S_{2}^{2}}{{{N}_{1}}+{{N}_{2}}-2}}$                   (10)

where, t: the student’s distribution (t-test); ${{\bar{X}}_{i}}:$ mean CN of the i-th sample; $\sigma $: standard deviation of the population; ${{N}_{i}}:$ size of the i-th sample, and $S_{i}^{2}:$ variance of the i-th sample. The test hypothesis is two-way (H1: Actual CN ≠ Projected CN) with the following rejection criteria: Reject H0 if “${{t}_{count}}$ < $-{{t}_{{}^{\alpha }/{}_{2}}}\left( table \right)$ ” or “${{t}_{count}}>{{t}_{{}^{\alpha }/{}_{2}}}\left( table \right)''$.

2.5.3 Forecasting monthly rainfall

Rainfall forecasting is carried out to predict monthly rainfall data in the Citarum-Majalaya catchment for 23 years (2022-2044). Then this rainfall data will be used to predict surface runoff discharge during that period. The rainfall forecasting process uses the Artificial Neural Network Backpropagation algorithm (ANN-BP). It is an information processing system with certain performance characteristics in common with biological neural networks [39]. It has advantages in terms of a small error rate and quite good generalization capabilities by using a training process according to weights to forecast future events [40]. This method was first introduced by Rumelhart et al. [41]. Backpropagation works is by adjusting the weight values backwards based on the error value in the learning process [42]. The outline of the stages of ANN analysis with the backpropagation algorithm is: (1) analyze the number of hidden layer neurons; (2) determine the number of output layer neurons; (3) initialization of weights and initial bias; (4) performance of network training on training data; (5) performance of network testing using test data; and (6) model validation.

2.5.4 Surface runoff

Accurate surface runoff estimation is one of the most fundamental and important points for planning and managing water resources and assessing the environmental quality of soil and water [43]. Several methods are available to determine surface runoff, namely the rational, Green-Ampt, and SCS-CN methods [44]. In this study, surface runoff was calculated using the SCS-CN method equation (Eq. (11)) because it is easy to apply and considers all catchment factors that influence runoff [45], and is suitable for small watersheds with an area of ≤250 km² [46]. This method is a popular rainfall-runoff model used to estimate losses and runoff directly from rainfall events [47].

$R=\frac{{{\left( P-0.2S \right)}^{2}}}{P+0.8S}$                            (11)

where, P is rainfall; S is a retention parameter that describes the maximum potential for rainfall to be infiltrated or become runoff. The S value is calculated by Eq. (12).

$S=25.4~\times ~\left( \frac{1000}{CN}-10 \right)$                           (12)