Assessing the Impact of Deforestation on Decadal Runoff Estimates in Non-Homogeneous Catchments of Peninsula Malaysia

Abstract

This study calibrated the Soil Conservation Service Curve Number (SCS-CN) model to predict decadal runoff in Peninsula Malaysia and found a correlation between the reduction of forest area, urbanization, and an increase in runoff volume. The conventional SCS-CN runoff model was found to commit a type II error in this study and must be pre-justified with statistics and calibrated before being adopted for any runoff prediction. Between 1970 and 2000, deforestation in Peninsula Malaysia caused a decline in forested land by 25.5%, resulting in a substantial rise in excess runoff by 10.2%. The inter-decadal mean runoff differences were more pronounced in forested and rural catchments (lower CN classes) compared to urban areas. The study also found that the CN value is a sensitive parameter, and changing it by ±10% can significantly impact the average runoff estimate by 40%. Therefore, SCS practitioners are advised not to adjust the CN value for better runoff modeling results. Additionally, NASA’s Giovanni system was used to generate 20 years of monthly rainfall data from 2001–2020 for trend analysis and short-term rainfall forecasting. However, there was no significant uptrend in rainfall within the period studied, and occurrences of flood and landslide incidents were likely attributed to land-use changes in Peninsula Malaysia.

1. Introduction

Humans have significantly altered the earth’s surface over time by transforming natural areas into industrial or agricultural regions. In the case of Malaysia, as a developing country, there has been extensive land-use and land-cover change since the 1970s. The rapid pace of development in Peninsula Malaysia has resulted in changes to the landscape and vegetation. These activities have had an uneven impact on the drainage basins, leading to changes in runoff patterns. Human activities, especially deforestation, have a greater impact on vegetation cover and the environment at local, regional, and global levels. Land-use and water cycle patterns are also influenced by anthropogenic activities, particularly deforestation, which has cleared the way for urban development. As a result, it is crucial to explore and estimate the impact of development activities on runoff patterns in rural catchments in Malaysia [1,2,3].

Forests are an essential element in sustaining our supply of water [4,5,6,7,8]. However, to make room for development, many forested areas have been cleared [9,10,11,12,13,14,15]. Between 2002 and 2020, Malaysia lost 27,000 km² (roughly 2/3 of Switzerland or nearly 39 Singapores) [16] of humid primary forest area at the clearing rate of 1421 km²/year within 19 years. This had resulted in a 17% reduction of humid primary forest in Malaysia within the given time period [17]. Compared to the deforestation of an area of 13,000 km² from 1978 to 1994 [18], the forest-clearing rate in Malaysia has doubled in recent decades. Deforestation challenged the quality and quantity of water [19] and caused significant hydrological changes which included an increase in runoff [20,21].

Interception losses in tropical and subtropical rainforests have been shown to vary from 6 to 42% of rainfall [22,23,24]. Meanwhile in Malaysia and Indonesia, forest interception loss was reported from 12.7% to 21% [25,26]. Forest removal will convert the interception losses into a contributing factor to increase surface runoff. Past research studies have investigated rainforest conversion and its hydrological impact [27,28,29,30,31,32,33,34]. Zhang et al. [35] reported that a reduction of forest cover from 22% to 10% caused changes in river discharge in China. Sharma et al. [36] concluded that under projected land-use scenarios, runoff would increase when forest areas were converted into agricultural land in the central Himalayas.

The conversion of natural land has increased the surface runoff while rapid urbanization and industrialization have been cited as the main causes of major flooding in Malaysia [37,38]. Many researchers quantified the effects on hydrological parameters as due to changes in vegetation cover or forest logging activities [39,40,41,42,43,44], while urbanization was also identified to be a key factor in landscape alteration, causing an increase in runoff. River discharge can also be affected by human interventions. The landscape alteration caused faster runoff from storms, increased peak flows and changed hydrologic cycles [45,46,47]. Sahin and Hall [48] reported that a 10% reduction in canopy cover resulted in a 20 to 25 mm increase in annual water yield while Bosh and Heewlett [39] found that a 10% removal increased the water yield by an average of 40 mm. Although flooding is a response to the complex hydrological system, the dynamic changes are caused by anthropogenic factors to mother nature.

The Soil Conservation Services (SCS) of the United States has developed a widely accepted methodology, called the Curve Number (CN), to estimate the direct runoff from rainfall in hydrological studies. This method is based on the concept that the amount of runoff produced by a rainfall event depends on various factors, such as land use, soil type, and the antecedent moisture condition of the soil. The SCS-CN method uses a CN value that represents the combined effect of these factors on runoff generation, ranging from 0 to 100. This value is used in an equation to estimate the amount of direct runoff from rainfall. The SCS-CN method is commonly used in watershed management, flood forecasting, and erosion control planning, and has been widely adopted by government agencies worldwide, integrated into software, and taught in hydrology courses.

However, studies around the world in recent decades reported runoff prediction accuracy problems with the conventional SCS-CN model [49,50]. Previous studies have utilized the SCS Curve Number (CN_0.2) to demonstrate variations in runoff response due to agricultural land-use and seasonal changes [51]. This paper applied the SCS-CN model calibration methodology with inferential statistics that was developed by authors in a previous study [50] and demonstrated the extended application to model decadal rainfall-runoff conditions in Peninsula Malaysia. The correlation between deforestation and urbanization on runoff increment in Peninsula Malaysia was also established.

2. Materials and Methods

2.1. Study Site and Data Collection

Peninsula Malaysia, with a land area of 132,265 km², is slightly larger than England (130,395 km²). It is bordered by Thailand to the north and Singapore across the strait of Johor to the south. This study utilized the calibrated SCS lump rainfall-runoff model in conjunction with the most recent rainfall-runoff dataset published by the Malaysian federal agency. The dataset, which can be located in the appendix of the Department of Irrigation and Drainage’s Hydrological Procedures no. 27 (DID HP 27), documented 227 storm events across 41 distinct catchments between 1970 and January 2000 in Peninsula Malaysia [52]. The smallest recorded storm event had a rainfall depth of 19 mm, with a measurable runoff depth of 4.8 mm, while the largest event measured 420 mm in rainfall depth and 258 mm in runoff depth.

In this study, the DID HP 27 dataset was divided into three periods to investigate the potential relationship between decadal runoff conditions, urban population growth in Peninsula Malaysia, and deforestation data provided by the Forestry Department Peninsula Malaysia [53,54,55,56,57,58,59]. This study grouped 58 recorded storm events from 1970 to 1979 as the M70 dataset, 81 events from 1980 to 1989 as the M80 dataset, and 88 events from 1990 to January 2000 as the M90 dataset for runoff analyses and comparison.

2.2. Calibration of SCS-CN Model

The SCS initially developed the CN rainfall-runoff model (Equation (1)) for the federal flood control program in 1954, and it has since become the basic model for runoff estimation.

$$Q=\frac{{\left(P-I_a\right)}^2}{P-I_a+S}$$

(1)

Q = Runoff depth (mm)

P = Rainfall depth (mm)

I_a = Rainfall initial abstraction amount (mm)

S = Catchment maximum water retention potential (mm)

where I_a = λS. If P < I_a, Q = 0.

The SCS also put forth the hypothesis that I_a = λS = 0.2S, where λ represents the initial abstraction ratio coefficient, which was proposed as a constant value of 0.2. The justification for Equation (1) was based on daily rainfall and runoff data, rather than event measurements, and its only official documentation source can be found in the National Engineering Handbook, Section 4 (NEH-4) [60,61]. By substituting I_a = 0.2S, Equation (1) is simplified into Equation (2):

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

(2)

DID HP 27 dataset was used to assess the SCS-CN rainfall-runoff model in authors’ previous study [50] with inferential statistics and concluded that the existing SCS-CN model is not even statistically significant at alpha = 0.05 level for runoff predictions in Peninsula Malaysia. Therefore, it must be calibrated according to the local rainfall-runoff dataset and the key parameter λ must be derived to formulate a statistically significant rainfall-runoff model for Peninsula Malaysia. All statistical analyses were performed using the IBM SPSS software version 26.0 in this study [62].

The SCS model was calibrated according to each decadal rainfall-runoff data batch and the calibrated runoff predictive models were formulated with the newly derived λ value to represent the aforementioned decadal rainfall-runoff conditions. Inter-decadal runoff difference can then be mapped with the runoff difference between newly formulated decadal rainfall-runoff models while non-parametric statistics were used for runoff trend analyses.

This study also re-assessed the validity of Equation (1) on the entire M70, M80, and M90 datasets through the reverse derivation of the λ value. Equation (2) will be discarded if the inferential statistics of the P-Q dataset reject the validity of λ = 0.2. Equation (1) will then be calibrated according to the P-Q data pairs of M70, M80, and M90. Non-parametric inferential statistics, bootstrapping, and Bias Corrected and Accelerated (BCa) procedure (2000 samples with replacement) were conducted in SPSS on each batch of derived λ and S values in order to search for optimum values within the confidence intervals at alpha 0.01 level. The λ confidence intervals’ span can be used to assess the null hypothesis H₀ as shown below:

Null Hypothesis (H₀): Equation (2) (λ = 0.2) is valid to model runoff estimates with DID HP 27 dataset.

The SCS-CN model calibration of this study consists of the following steps:

• Rearrange Equation (1) into: $S=\frac{{(P-I_a)}^2}{Q}+I_a-P$
• For each P-Q event pair, calculate the corresponding S value with the above equation under the SCS constraint that I_a < P value and calculate the λ value with λ = I_a/S.
• Conduct bootstrap, BCa (at α = 0.01 level) inferential statistical analyses (2000 samples with replacement) for the calculated λ and S datasets separately for each decadal model.
• Generate 99% confidence intervals for λ and S datasets for each decadal model.
• Test the null hypothesis (H₀) by referring to the λ confidence intervals’ span and its standard deviation for each decadal model. If the λ = 0.2 value exists within the λ confidence interval, use Equation (2) to model rainfall-runoff. Otherwise, move to step 6.
• Find the optimum λ and S values from BCa confidence intervals and calculate I_a for each decadal model using supervised optimization technique by minimizing the overall model bias (BIAS) near to the value of zero.
• Formulate the calibrated SCS model by substituting I_a and S into Equation (1).
• According to a group of researchers [63], when λ value other than 0.2 is detected, its corresponding S value (denoted by S_λ) must be correlated to the S_0.2 values for CN calculation. As such, correlate S_λ and S_0.2 with the S general formula which was derived by a past researcher [50]: $S_{\lambda }=\frac{\left[P-\frac{\left(\lambda -1\right)Q}{2\lambda }\right]-\sqrt{PQ-P^2+{\left[P-\frac{\left(\lambda -1\right)Q}{2\lambda }\right]}^2}}{\lambda }$
• Substitute optimum λ and S_λ into Equation (1) to formulate the decadal model.
• Lastly, substitute $S_{0.2}=\frac{\mathrm{25,400}}{{\mathrm{C}\mathrm{N}}_{0.2}}-254$ into each decadal model to express Q in term of P and CN_0.2.

Note: Appendix B shows step 9 and 10 using an example.

2.3. IMERG Satellite Rainfall Trend Analysis

The Integrated Multi-Satellite Retrievals for Global Precipitation Measurement (IMERG) [64] is a unified satellite product that has been launched by NASA and JAXA. This product combines, calibrates, and integrates satellite microwave precipitation estimates with microwave-calibrated infrared satellite estimates, precipitation gauge analyses, and other precipitation estimators, in order to provide precipitation measurements at fine time and space scales across the globe. There are three types of IMERG systems: Early (IMERG-E), Late (IMERG-L), and Final (IMERG-F) runs, with latency times of 6 h, 18 h, and 3 months, respectively. All IMERG runs are available at half-hourly, daily, and monthly temporal resolutions, and have global coverage at a spatial resolution of 0.1°.

In this study, the IMERG-F monthly version 6 product at a spatial resolution of 0.1° was employed, as it has been bias-corrected using precipitation gauges from the Global Precipitation Climatology Centre (GPCC) and is considered more accurate for scientific research compared with other IMERG products. The data can be accessed through NASA’s Giovanni system [65]. Monthly rainfall data from 2001 to 2020 (Appendix A) was obtained from the Giovanni system to demonstrate the annual rainfall distribution across Malaysia.

Trend analysis was conducted on the monthly rainfall data in Malaysia. Before conducting the analysis, a normality test was performed to determine which central of tendency (mean or median) was to be used. If a dataset has less than 2000 samples, the Shapiro-Wilk test is recommended instead of the Kolmogorov-Smirnov test. This study used a dataset with less than 2000 samples; therefore, the Shapiro-Wilk test was employed. If the p-value is greater than 0.05, the dataset is considered to be normally distributed, and therefore, mean is used to measure the central tendency of the dataset [50]. The non-parametric inferential statistics using the BCa bootstrapping method was conducted on 2000 random samples (with replacement) to calculate the 99% confidence intervals of the mean or median for each interval.

Rainfall time series forecasting models were created using monthly rainfall data in Malaysia from January 2001 to December 2020 (N = 240) with the Expert Modeler in SPSS. The rainfall amount was forecasted for the next 24 months, from January 2021 to December 2022 in Malaysia, to determine if there will be a significant change in rainfall pattern in the near future that could affect runoff prediction.

3. Results and Discussion

3.1. The Optimum λ and S of Decadal Models

Table 1. Inferential statistics results of derived λ for M70 decadal dataset at α = 0.01.

λ	Statistics	Bootstrap, BCa 99%
		Bias	Std. Error	Confidence Intervals
				Lower	Upper
Skewness	3.815
Kurtosis	19.768
Mean	0.098	0.0003	0.014	0.069	0.142
Median	0.065	0.0006	0.00294	0.049	0.089

,

Table 2. Inferential statistics results of derived λ for M80 decadal dataset at α = 0.01.

λ	Statistics	Bootstrap, BCa 99%
		Bias	Std. Error	Confidence Intervals
				Lower	Upper
Skewness	3.217
Kurtosis	13.028
Mean	0.095	0.0002	0.014	0.066	0.135
Median	0.047	0.0022	0.007	0.034	0.064

and

Table 3. Inferential statistics results of derived λ for M90 decadal dataset at α = 0.01.

λ	Statistics	Bootstrap, BCa 99%
		Bias	Std. Error	Confidence Intervals
				Lower	Upper
Skewness	5.393
Kurtosis	34.674
Mean	0.076	0.00005	0.013	0.051	0.115
Median	0.042	0.00183	0.007	0.031	0.063

display the BCa 99% confidence intervals of λ for the M70, M80, and M90 decadal datasets. The BCa 99% confidence intervals of both the mean and median for each decadal dataset do not include the value of 0.2, leading to the rejection of the null hypothesis (H₀) at the alpha = 0.01 level. Equation (2) was found to be statistically invalid and, thus, cannot be utilized to model runoff conditions in Peninsula Malaysia for the M70, M80, and M90 decadal datasets. The rejection of H₀ necessitates the search for a new, optimal value of λ to develop a new rainfall-runoff prediction model.

The λ dataset is skewed and tested to be non-normally distributed in SPSS for all three decadal groups and therefore, the search for an optimal representative λ value using a supervised optimization technique will be concentrated on the range of median confidence intervals. These intervals are [0.049, 0.089] for the M70 dataset (Table 1), [0.034, 0.064] for the M80 dataset (Table 2), and [0.031, 0.063] for the M90 dataset (Table 3).

The BCa 99% confidence intervals of S_λ for the M70, M80, and M90 decadal datasets are presented in

Table 4. Inferential statistics results of derived S_λ for M70 decadal dataset at α = 0.01.

Sλ	Statistics	Bootstrap, BCa 99%
		Bias	Std. Error	Confidence Intervals
				Lower	Upper
Skewness	1.298
Kurtosis	0.975
Mean	151.592	−0.482	14.954	117.083	187.008
Median	123.615	−0.166	11.367	91.255	157.815

,

Table 5. Inferential statistics results of derived S_λ for M80 decadal dataset at α = 0.01.

Sλ	Statistics	Bootstrap, BCa 99%
		Bias	Std. Error	Confidence Intervals
				Lower	Upper
Skewness	1.794
Kurtosis	5.201
Mean	180.994	−0.339	14.954	141.892	231.088
Median	147.950	−3.921	19.112	113.890	183.670

and

Table 6. Inferential statistics results of derived S_λ for M90 decadal dataset at α = 0.01.

Sλ	Statistics	Bootstrap, BCa 99%
		Bias	Std. Error	Confidence Intervals
				Lower	Upper
Skewness	1.132
Kurtosis	1.407
Mean	161.827	0.221	11.934	131.989	192.939
Median	142.610	−2.566	21.298	95.236	189.506

. The normality of the S_λ dataset was tested using SPSS for all three decadal groups, and found to be normally distributed. Therefore, the optimal S_λ value was searched for within the range of the mean confidence intervals. There intervals are [117.083, 187.008] for M70 dataset (Table 4), [141.892, 231.088] for M80 dataset (Table 5), and [131.989, 192.939] for M90 dataset (Table 6).

The optimal λ and S_λ values for the M70, M80, and M90 decadal datasets using a supervised optimization technique are presented in

Table 7. Optimal λ, S_λ and I_a for M70, M80 and M90 decadal datasets.

Dataset	Optimal λ	Optimal Sλ	Ia = λSλ
M70	0.049	160 mm	7.904 mm
M80	0.034	190 mm	6.431 mm
M90	0.031	160 mm	4.956 mm

. The product of the optimal λ and S_λ values gives the representative initial abstraction value for each dataset, which can be calculated as I_a = λS_λ.

3.2. The Decadal Rainfall-Runoff Models

The decadal rainfall-runoff models for Peninsula Malaysia are presented in

Table 8. Decadal rainfall-runoff models for M70, M80 and M90 decadal datasets.

Dataset	Runoff Predictive Model	Nash-Sutcliffe Index	Equation Number
M70	$Q=\frac{{\left(P-7.904\right)}^2}{P-152.096}$	0.958	(3)
M80	$Q=\frac{{\left(P-6.431\right)}^2}{P-183.569}$	0.910	(4)
M90	$Q=\frac{{\left(P-4.956\right)}^2}{P-155.044}$	0.907	(5)

Q = runoff depth (mm), P = rainfall depth (mm).

by substituting the respective I_a and S_λ values from Table 7 into Equation (1). Equations (3)–(5) were then formulated to model the decadal rainfall-runoff conditions in Peninsula Malaysia. To further analyse decadal runoff trend across multiple rainfall depths and CN_0.2 scenarios in Peninsula Malaysia, Equations (3)–(5) need to be re-expressed in terms of CN_0.2 to benefit SCS practitioners as they are more familiar with the use of curve number [50].

To calculate S_λ and S_0.2 for each decadal dataset, the general S_λ formula (step 8 in methodology Section 2.2) can be used with the optimum λ values. Through SPSS, this study identified statistically significant power function correlation between S_λ and S_0.2 for the M70, M80, and M90 decadal datasets, which is consistent with previous research findings [66,67,68]. The final equations are listed in

Table 9. Correlation equations between S_λ and S_0.2 for M70, M80 and M90 decadal datasets.

Dataset	Correlation Equation	Adjusted R-Squared	Standard Error of Estimate	Equation Number
M70	$S_{0.049}=1.184{S_{0.2}}^{1.081}$	0.939	0.134	(6)
M80	$S_{0.034}=1.107{S_{0.2}}^{1.094}$	0.910	0.201	(7)
M90	$S_{0.031}=1.179{S_{0.2}}^{1.069}$	0.907	0.165	(8)

All correlations are significant at p < 0.001.

.

Equations (6)–(8) played a crucial role in converting S_λ to S_0.2, enabling SCS practitioners to use a rainfall-runoff model with CN_0.2, which they are more familiar with. Furthermore, by establishing a correlation between the newly derived S_λ and S_0.2, Equations (3)–(5) were modified to be expressed in CN_0.2 terms, facilitating decadal trend analyses with CN_0.2.

Equations (3)–(5) can be expressed in CN_0.2 by substituting S_λ in Equation (1) with Equations (6)–(8), as well as the SCS-CN formula (Step 10 in methodology Section 2.2). By doing so, the decadal runoff predictive models can be re-expressed as shown in Appendix B. The resulting alternate representations for decadal runoff predictive models in Peninsula Malaysia are presented in

Table 10. Alternate form of decadal rainfall-runoff models for M70, M80 and M90 decadal datasets.

Dataset	Runoff Predictive Model	Equation Number
M70	$Q_{0.049}=\frac{{\left[P-23.077{\left(\frac{100}{{\mathrm{C}\mathrm{N}}_{0.2}}-1\right)}^{1.081}\right]}^2}{\left[P+447.876{\left(\frac{100}{{\mathrm{C}\mathrm{N}}_{0.2}}-1\right)}^{1.081}\right]}$	(9)
M80	$Q_{0.034}=\frac{{\left[P-15.992{\left(\frac{100}{{\mathrm{C}\mathrm{N}}_{0.2}}-1\right)}^{1.094}\right]}^2}{\left[P+456.589{\left(\frac{100}{{\mathrm{C}\mathrm{N}}_{0.2}}-1\right)}^{1.094}\right]}$	(10)
M90	$Q_{0.031}=\frac{{\left[P-13.618{\left(\frac{100}{{\mathrm{C}\mathrm{N}}_{0.2}}-1\right)}^{1.069}\right]}^2}{\left[P+425.963{\left(\frac{100}{{\mathrm{C}\mathrm{N}}_{0.2}}-1\right)}^{1.069}\right]}$	(11)

in term of CN_0.2.

3.3. The Decadal Runoff Trend Analyses

The decadal runoff models (Equations (9)–(11)) enable the quantification of runoff conditions for various decades under different rainfall depths (P) and CN_0.2 scenarios, facilitating the analysis of runoff changes. The DID HP 27 dataset contains the lowest and highest recorded rainfall depths, ranging from 20 mm to 430 mm across CN_0.2 classes from 46 to 94. Runoff difference tables can then be calculated between any two decadal models.

This study evaluated the inter-decadal runoff differences between M70 and M80, M80 and M90, and M70 and M90 in Peninsula Malaysia. The runoff amount of the earlier decade was subtracted from the latter to determine the inter-decadal runoff difference. For instance, the inter-decadal runoff difference between M70 (Equation (9)) and M80 (Equation (10)) was calculated by subtracting the runoff amount of M70 from that of M80. A positive inter-decadal runoff difference amount indicated a larger runoff amount in M80 compared to M70 and vice versa. Statistical analyses were performed to determine significant runoff trends between different decades. This study also correlated decadal runoff changes with deforestation and urbanization data in Peninsula Malaysia.

Non-parametric Kendall’s Tau b and Spearman’s Rho statistics were used to evaluate the inter-decadal runoff trend in SPSS. Both statistics showed a significant positive correlation (2-tailed) at alpha = 0.01 for all inter-decadal periods, rainfall depths, and CN_0.2 classes mentioned above. This positive correlation indicates an upward trend in inter-decadal runoff, which can be visually represented in Figure 1, Figure 2 and Figure 3. To assess the magnitude of this upward trend in each inter-decadal scenario and CN_0.2 class (ranging from 46 to 94) according to rainfall depths from 20 mm to 430 mm, Sen slopes and its collective inferential statistics were calculated. The Sen slopes and inferential statistics of all CN_0.2 classes were then analysed collectively for each inter-decadal scenario at the alpha = 0.01 level, and the results are tabulated in

Table 11. Inferential statistics of Sen Slopes for inter decadal runoff difference between M80 and M90.

Sen Slopes M80 to M90	Statistics	Bootstrap, BCa 99%
		Bias	Std. Error	Confidence Intervals
				Lower	Upper
Skewness	−0.022
Kurtosis	−1.512
Mean	0.0121	−0.00002	0.00225	0.00701	0.01720
Median	0.0127	−0.00029	0.00422	0.00439	0.02119
Std. Deviation	0.0087	−0.00039	0.00103	0.00588	0.01032
Range	0.0247

,

Table 12. Inferential statistics of Sen Slopes for inter decadal runoff difference between M70 and M80.

Sen Slopes M70 to M80	Statistics	Bootstrap, BCa 99%
		Bias	Std. Error	Confidence Intervals
				Lower	Upper
Skewness	0.272
Kurtosis	−0.940
Mean	0.0051	−0.00001	0.00081	0.00313	0.00713
Median	0.0048	0.00008	0.00126	0.00230	0.00824
Std. Deviation	0.0031	−0.00015	0.00046	0.00178	0.00396
Range	0.1000

and

Table 13. Inferential statistics of Sen Slopes for inter decadal runoff difference between M70 and M90.

Sen Slopes M70 to M90	Statistics	Bootstrap, BCa 99%
		Bias	Std. Error	Confidence Intervals
				Lower	Upper
Skewness	0.065
Kurtosis	−1.473
Mean	0.0178	−0.00003	0.00322	0.01049	0.02506
Median	0.0175	0.00015	0.00589	0.00712	0.03125
Std. Deviation	0.0124	−0.00057	0.00149	0.00823	0.01479
Range	0.0353

.

Calculated Sen slope values were tested as normally distributed in SPSS. Therefore, the mean Sen slope value was used to represent each inter-decadal runoff scenario. On average, the collective Sen slope value was 0.0121 (p = 0.01, 99% confidence interval from 0.00701 to 0.01720) to indicate the runoff incremental trend between M80 and M90. The Sen slope between M70 and M80 was 0.0051 (p = 0.01, 99% confidence interval from 0.00313 to 0.00713), while between M70 and M90 it was 0.0178 (p = 0.01, 99% confidence interval from 0.01049 to 0.02506). The Sen slope values also estimated the percentage of rainfall depth that becomes incremental runoff. For instance, the average expected runoff increment from a rainfall depth of 100 mm across CN_0.2 classes from 46 to 94 can be estimated to be 1.21 mm between M80 and M90 (i.e., 0.0121 × 100 mm).

The study conducted a repetition of all statistical analyses with a CN_0.2 range of 46 to 70 to assess the runoff changes across lower CN_0.2 classes, which correspond to rural and forested catchments. This was done to obtain a more accurate estimate of the inter-decadal runoff increment conditions in these areas. The results showed that the runoff incremental trend between M80 and M90 of CN_0.2 (46 to 70) had a Sen slope value of 0.0190 (p = 0.01, 99% confidence interval from 0.01595 to 0.02216). The Sen slope value between M70 and M80 was 0.0075 (p = 0.01, 99% confidence interval from 0.00606 to 0.00898), and between M70 and M90, it was 0.0276 (p = 0.01, 99% confidence interval from 0.02262 to 0.03239). For example, the expected runoff increment from rainfall of 100 mm across CN_0.2 classes from 46 to 70 was estimated to be 1.90 mm between M80 and M90. The study found that runoff increments were significant (p = 0.01) between all inter-decadal scenarios and were more apparent in forested and rural areas (highlighted area in Figure 4).

Positive inter-decadal runoff difference in Peninsula Malaysia is depicted in Figure 1, Figure 2 and Figure 3. High rainfall depths and low CN_0.2 groups, which are associated with forested catchments, are particularly affected. These study outcomes are in line with previous studies [67,68,69,70]. Inter-decadal runoff differences are more pronounced under high rainfall depths. The mean runoff of different decades across different CN_0.2 classes was calculated and compiled, as shown in Figure 5. M90 had the highest runoff, while M70 had the lowest. Greater percentage changes in mean runoff were observed in lower CN_0.2 classes (forested catchments) compared to higher CN_0.2 classes (urban area). The largest mean runoff incremental percentage was 12.6% (6.6 mm) from M70 to M90 at CN_0.2 = 46, while the smallest change was 0.1% (0.1 mm) from M80 to M90 at CN_0.2 = 94.

3.4. The Impact of CN_0.2 Variation on Runoff

According to [71], due to variation in hydrological conditions, CN_0.2 value is often calibrated to match observed runoff dataset in modelling practice. Researchers observed that a variation of ±10% in CN_0.2 might lead to ±50% runoff variation [72] while [73] it was reported that even 1% increase in CN_0.2 with rainfall depth of 175 mm had caused 2.03% increase in runoff. References [73,74] concluded that CN_0.2 variations will have a larger impact on runoff than other parameters in Equation (1).

CN_0.2 tweaking becomes a convenient way to calibrate and validate hydrological models. However, other studies reported that CN_0.2 value of a catchment was unstable and decreased when rainfall increased [74,75,76]. The error and sensitivity analysis results by some researchers stated that CN_0.2 variations would induce a larger impact on runoff calculation with inherent error rather than rainfall depth variations [72,74,77]. CN_0.2 tweaking might achieve or improve temporal hydrological modelling accuracy through the trial-and-error technique, but the practicality of the end result was often uncertain and lack of statistical justification.

This study modelled the impact of CN_0.2 variation with the DID HP 27 dataset. According to [78], the practical CN values were likely to be within the range of 40 to 98. The optimum best collective CN_0.2 was 71 for the entire DID HP 27 dataset, thus, CN_0.2 variation up to 40% was chosen to cover the range of CN_0.2 from 43 to 99 and rainfall from 20 mm to 430 mm. CN_0.2 upscaling induced larger runoff change than downscaling while both effects were largely felt at rainfall depths below 100 mm. On average, runoff would reduce by 37% when CN_0.2 was downscaled up to 40% between 20 mm and 430 mm. On the other hand, the average runoff increased by 306% when CN_0.2 was upscaled up to 40%. The average runoff for both scenarios was almost identical when rainfall depths were limited to higher rainfall depths (100 mm to 430 mm). The average runoff reduced by 34% when CN_0.2 was downscaled up to 40%, while average runoff increased by 35% when CN_0.2 was upscaled to the same range. Varying the CN_0.2 value by ±10% resulted in an average runoff change of 40%, which is consistent with the findings reported in [72]. Similarly, upscaling the CN_0.2 value by 1% with a rainfall depth of 175 mm caused a 2% increase in runoff, which matches the range reported by [73]. Sen slope analyses showed that in both CN_0.2 upscaling and downscaling scenario, runoff reduction and incremental rates reduced toward the high rainfall depths but increased according to the CN_0.2 variation percentage. Lower rainfall depths (20 to 100 mm) had higher runoff variation percentages than higher rainfall depths (100 to 430 mm), as reported by previous studies [67,68,69,70]. Figure 6 and Figure 7 present the overview of the impact of CN_0.2 variation on runoff with equations to estimate the percentage change in runoff.

3.5. The Impact of Deforestation and Urbanization on Runoff in Peninsula Malaysia

According to statistics from the Malaysia Department of Forestry, Peninsula Malaysia went through extensive deforestation from the 1970s. Forest area decreased at fast rates in the 1980s and started to stabilise in the 1990s. Forest area had reduced by 21.1% from M70 to M80 and 25.5% from M70 to M90 (Figure 8). Figure 9 was created to show the relationship between these decadal forest area reduction rate and its corresponding mean inter-decadal incremental runoff difference (Q_v%) across different CN_0.2 classes in Peninsula Malaysia. During the period between M70 and M90 in Peninsula Malaysia, the mean excess (incremental) runoff volume difference for CN_0.2 classes ranging from 46 to 70 was calculated to be 6.8 mm, equivalent to 6.8 million litres per square kilometre. This corresponds to a 10.2% increase in excess runoff, and it occurred simultaneously with a 25.5% decrease in forest area. These findings provide insights into the hydrological impacts of deforestation on non-homogeneous catchments. In general, inter-decadal mean runoff differences were more pronounced in forested and rural catchments (lower CN_0.2 classes) than urban areas. Inter-decadal runoff difference between M70 and M90 is significantly greater than runoff difference between M70 and M80 (Figure 9).

According to the published data and figures from the Department of Statistics Malaysia [79,80,81,82,83,84,85,86,87], the urban population in Peninsula Malaysia had been increasing rapidly (Table 14). In comparison to the forest area statistics from the Department of Forestry [53,54,55,56,57,58,59], an inverse correlation was identified in SPSS as:

FA = 5.533 + (4.922/Urb-pop)

(12)

R²_adj = 0.964, SE = 0.175, p < 0.012

FA = Forest area (Million hectare)

Urb-pop = Urban population in Peninsula Malaysia (Millions)

The inverse correlation between urban population and the forest area in Peninsula Malaysia implies that urban development has significant correlations with deforestation (Figure 8). On the other hand, the deforestation has a direct impact on runoff amount as shown in Figure 9.

Table 14. Decadal urban population and forest area in Peninsula Malaysia [12,80,81,82,83,84,85,86,87].

Year	Urban Population (Millions)	Forest Area (Millions Hectare)
1970	2.03	8.01
1980	4.81	6.35
1990	7.97	6.27
2000	12.26	5.97

Table 14. Decadal urban population and forest area in Peninsula Malaysia [12,80,81,82,83,84,85,86,87].

Year	Urban Population (Millions)	Forest Area (Millions Hectare)
1970	2.03	8.01
1980	4.81	6.35
1990	7.97	6.27
2000	12.26	5.97

3.6. Decadal λ and I_a

In recent decades, urbanization in Peninsula Malaysia has caused uneven land development, leading to non-homogeneity in catchments. This study aimed to address this issue by calibrating the SCS-CN model using rainfall-runoff data from different decades to develop decadal models. The models demonstrated a strong ability to estimate runoff amounts, achieving a Nash-Sutcliffe Index ranging from 0.907 to 0.958 (Table 8), even in non-homogeneous catchments. These findings suggest that recalibrating the SCS-CN models based on regional and decadal specific rainfall-runoff conditions could be an effective approach for estimating runoff in non-homogeneous catchments.

In a previous study by the authors, the optimum λ value for the entire DID HP 27 dataset was identified as 0.051 to model overall runoff conditions. However, in this study, different optimum λ values were derived for the decadal datasets of M70, M80, and M90, which also led to changes in the corresponding initial abstraction (I_a) values (Table 7). Over time, the optimal λ and I_a values for each decade were found to decrease, indicating changes in land cover resulting from deforestation and urbanization that impact runoff conditions in rural catchments. The decreasing trend in λ leads to a corresponding increase in runoff over time in Peninsula Malaysia.

SCS practitioners commonly calibrate CN_0.2 with one batch of runoff data and validate the final results against another batch to determine the optimum CN_0.2 value for modelling a combined dataset. However, this study highlights concerns with this practice due to land-use and land-cover changes in Peninsula Malaysia, which directly affect catchment runoff conditions over time.

There is a statistically significant upward trend in runoff (at alpha = 0.01) across all CN_0.2 classes from M70 to M90 due to changes in land use. Therefore, SCS practitioners must be cautious and aware that blindly accepting the λ value as 0.2 is not advisable, and it is strongly recommended to derive a regional-specific λ value. Although an optimum λ value of 0.051 was used in a previous study [50] to model the entire dataset with a Nash-Sutcliffe value of 0.92, it differed significantly from the optimum λ values of different decades. Hence, runoff predictive models formulated with different optimum λ values will yield differences in runoff predictions.

3.7. Rainfall Trend Analyses

The BCa bootstrapping method with a 99% confidence interval was used to analyse the monthly rainfall trend in Peninsula Malaysia (Figure 10), which revealed that there has been no significant trend in the monthly rainfall over the past 20 years. The forecasted rainfall from 2021 to 2022 was consistent with the current trend. The results were supported by the model generated by Expert Modeler (Figure 11), which indicated that there would be no significant alteration in the rainfall trend in the near future.

Flood occurrences are strongly influenced by changes in land use, including deforestation, agricultural activities, and urbanization. The conversion of natural land cover to urban and other developed land use can significantly alter the hydrological cycle, resulting in increased surface runoff and reduced infiltration. This alteration of the landscape can lead to changes in the frequency, magnitude, and timing of floods, as well as increased erosion and sedimentation in rivers and streams [37].

The increase in surface runoff resulting from the conversion of natural land for human use has caused major flooding in Malaysia [88]. Floods are the most destructive natural disasters in the country, with an estimated 85 river basins, mainly in Peninsula Malaysia, prone to recurrent flooding. Based on a study in 2014, approximately 9% of the total area of Malaysia, covering 29,800 km², is vulnerable to flood disaster, affecting almost 4.82 million people, equivalent to 22% of the total population [89]. In 2014, the states of Johor, Kelantan, Pahang, Perak, and Terengganu in Peninsula Malaysia, which were severely affected by floods, also recorded high rates of forest loss [38].

Recent floods and landslides in Malaysia have been attributed to deforestation, which results in the release of sediment and weakened soil. Trees help to prevent sediment runoffs and hold water, making them an essential factor in maintaining a stable environment. The excessive clear-cutting of trees for oil palm plantations has been identified as the primary cause of mudslides in recent times, with poor construction standards also contributing to the problem. Therefore, it can be inferred that deforestation and poor construction standards are the key factors responsible for these events, rather than El Niño or global warming, as suggested by studies [90,91,92,93,94].

4. Conclusions

This study adopted the CN hydrological model calibration technique developed by the authors in a previous study [50] and applied it in decadal runoff prediction study in Peninsula Malaysia. The correlation between forest area reduction, urbanization, and runoff volume increment was established. Highlights of the study are summarized as below:

• The use of the conventional SCS-CN runoff model will commit type II error in this study to predict runoff conditions of different study periods. It must be pre-justified with statistics and calibrated prior to adoption for any runoff prediction. It is also not recommended to conduct calibration and validation on the entire DID HP 27 dataset of this study as each demarcated decade was represented with its unique and statistically significant runoff predictive model. Calibration and validation methodology based on the conventional SCS-CN runoff model fail to quantify accurate runoff conditions spanning across different time periods with significant land-use and cover change.
• CN adjustment practice to formulate a hydrological model can have a large inherent error as small adjustments on the curve number can lead to large variation in the runoff. Given sufficient sample size, SCS-CN runoff model should be calibrated and formulated according to its unique optimum λ values to represent rainfall-runoff conditions of different time periods. In this study, when CN value was varied ± 10%, the average runoff changed by 40%. This study found a significant increase in runoff across all CN_0.2 classes in Peninsula Malaysia due to changes in land use, emphasizing the importance of deriving a regional-specific λ value and cautioning that different optimum λ values for different decades will yield differences in runoff predictions.
• This study emphasizes the significance of accounting for regional and decadal-specific rainfall-runoff conditions to estimate runoff in non-homogeneous catchments effectively. The calibrated SCS-CN model using data from different decades showed a remarkable ability to accurately estimate runoff amounts, even in non-homogeneous catchments. The models achieved a strong ability to estimate runoff amounts, attaining a Nash-Sutcliffe Index ranging from 0.907 to 0.958, even in non-homogeneous catchments.
• Calibrated SCS decadal (lump) runoff models show significant decadal runoff uptrend which coincides with the overall deforestation rate in Peninsula Malaysia. The presented methodology may become more apparent with regional specific deforestation rate and its corresponding rainfall-runoff dataset. The reduction of forest area by 25.5% in Peninsula Malaysia between 1970 and 2000 was found to be directly proportional to an increase in excess runoff volume of 10.2%. In general, inter-decadal mean runoff differences were more pronounced in forested and rural catchments (lower CN classes) than urban areas.
• NASA’s Giovanni system was used to generate 20 years of annual rainfall maps while monthly rainfall data (2001 to 2020) was also extracted for trend analysis and short-term forecast. This study found no significant uptrend in the rainfall within the period, and the occurrence of flood and landslide incidents can likely be attributed to land-use changes in Peninsula Malaysia.