Assessing the impacts of land use–land cover changes on direct surface runoff: a remote sensing approach in Khulna City Open Access

The massive urban expansions that have occurred over the past few decades have altered land use-land cover (LULC) change through different anthropogenic activities on landscape (Vojtek & Vojteková 2019; Moniruzzaman et al. 2020). As a part of the urbanization process, large agricultural areas and other non-urban land are regularly converted into impervious land, resulting in land use shifts that instantaneously disrupt the natural hydrological processes (Hu et al. 2020). For urban planning, water resource management, and early flood warning in major cities, quantitative measurement of the effects of urbanization on surface runoff is critical (Vojtek & Vojteková 2019; Hu et al. 2020; Moniruzzaman et al. 2020).

By 2050, 64% of the developing nations and 86% of the developed nations are expected to become fully urbanized (United Nations Department of Economic and Social Affairs 2014; Hu et al. 2020). Asia, like the other continents, has witnessed significant LULC changes in the twenty-first century which has caused a land use shift of 6 million hectares per year at the expense of forest cover from 2000–2005 compared with 1990–2000 (Lambin et al. 2003; Lindquist et al. 2012). Bangladesh has ranked 4th among the Asian countries in terms of urban growth with a growth rate of 4.8% between 2000 and 2011 (World Bank, 2007). As the 3rd largest metropolitan city of Bangladesh, Khulna has experienced a drastic growth rate of 3.8% which has changed the urban texture of the city with unforeseen LULC changes over the past two decades (Moniruzzaman et al. 2018). This has caused serious disruption to natural drainage patterns and damage to the flood retention zones (Rahman et al. 2009) that account for unprecedented alteration to hydrological processes like increased surface runoff (Rahman et al. 2009; Hu et al. 2020), reduced runoff response time and changed hydrological regimes (Hu et al. 2020).

Classification of LULC is one of the very fundamental tasks prior to any decision making regarding LULC change detection (Weng 2001; Miller et al. 2014; Roy et al. 2014; Addae & Oppelt 2019). Maximum likelihood classification (MLC) and Artificial neural network (ANN) have been considered among the frequently used methods for LULC classification (Candade & Dixon 2004; Bahari et al. 2014). However, these methods have their own limitations; ANN accounts for overfitting and local minima problems, while MLC requires large amounts of training samples and assumes the data are normally distributed which is rarely met for unbalanced real world data (Mountrakis et al. 2011). Support vector machine (SVM) is one of the more reliable classification techniques that has been developed in recent years. SVM is distinguished by an effective hyperplane searching strategy that uses a small training area and hence takes less time to process. This method avoids the problem of overfitting and makes no assumptions about the data type. Despite being non-parametric, the technique can develop effective decision boundaries and hence reduce misclassification (Vapnik 1995; Bahari et al. 2014). SVM is a binary classifier that works by determining the best hyperplane and correctly classifying data points into two groups. The classifier can be expanded from binary to multiclass using a variety of strategies, including one versus all and one against one (Sap & Kohram 2008). SVM can categorize data in both linear and nonlinear ways. For nonlinear data, the kernel function is utilized (Sap & Kohram 2008; Huang et al. 2010; Mountrakis et al. 2011; Karimi et al. 2019). Earlier studies done by Foody and Mathur showed that when comparing numerous classification approaches, such as discriminant analysis, decision trees, feed forward neural networks, and SVM classification; SVM produced the best results (Foody & Mathur 2004). Similarly, Shi and Yang found that SVM outperforms MLC in terms of quantitative accuracy (Shi & Yang 2015). Candade and Dixon also examined the result of SVM with ANN and claimed that SVM demonstrated greater performance even with a small number of training samples (Candade & Dixon 2004).

Again, several well established models are in practice for the quantification of the rainfall–runoff relationship viz. variable infiltration capacity (VIC) hydrological model, water flow and balance simulation model (WaSiM), soil and water assessment tool (SWAT), etc. (Beven & Kirkby 1979; Thakur et al. 2017; Moniruzzaman et al. 2020). Almost all of these models demand extensive data, rigorous calibration, and site-specific applicability. The collection, management, and analysis of these data are costly, time-intensive, and challenging as well (Beven & Kirkby 1979; Thakur et al. 2017). Hence, Soil Conservation Service-Curve Number (SCS-CN) method is predominantly practiced for catchments all over the world for the estimation of direct runoff potential from rainfall as it is simple, requires less data, and has static assumptions (Hu et al. 2020; Moniruzzaman et al. 2020; Psomiadis et al. 2020). As SCS-CN does not require rigorous calibration like the other models, it can be applied to ungauged a smaller watershed where measured runoff data are unavailable (Zare et al. 2016). With its high storage, superior data handling, data processing, data analysis, and manipulation capabilities, GIS and remote sensing techniques can enrich this model's input parameters by providing data for inaccessible locations and high spatio-temporal variability (Gajbhiye 2015). In 2003, Nayak and Jaiswal discovered that utilizing GIS and CN, there was a high correlation between measured and estimated runoff depth (Nayak & Jaiswal 2003). This technique is easy to understand and implement because it relies on a limited number of readily available parameters that are responsive to changing land use patterns and can be applied to a wide range of basin types and sizes (Nayak & Jaiswal 2003; Bansode & Patil 2017; Ara & Zakwan 2018; Psomiadis et al. 2020).

In this study, the Khulna City Corporation (KCC) area was selected as the research case. LULC classification of satellite images for the years 2005, 2010, 2015 and 2020 was performed using the SVM algorithm. Geospatial techniques (GIS and remote sensing) and the SCS-CN model have been incorporated to explore the association between land use change and surface-runoff potential. The major objectives of this study were to (a) detect the spatio-temporal LULC changes from 2005 to 2020, (b) estimate the potential surface runoff and its distribution over the study space due to the changes in land use, (c) evaluate the impact of land use change on surface runoff potential.

Khulna is the 3rd largest city of Bangladesh and the country's second most important port of entry (BBS 2011). KCC has a total area of 45.65 square km. and a population of 15,00,689 people (KCC 2021), with an annual population growth rate of 0.68 percent (BBS 2011). It is split into 31 wards (Figure 1) defined by a geographical boundary between 22°46′ and 22°58′ north latitude, and 89°28′ and 89°37′ east longitudes (Moniruzzaman et al. 2018). The city is surrounded by Bhairab and Mayur rivers on the north, and the Rupsa River at the core of the city (Roy et al. 2005). Rupsa is a river highly prone to flooding during the rainy season and high tides. Khulna's climate is humid during summer, and sunny during winter, with an average annual temperature of 26.3 °C (79.3°F) and monthly average temperatures ranging from 12.4 °C (54.3°F) in January to 34.3 °C (93.7°F) in May (Moniruzzaman et al. 2018). Khulna's average annual total rainfall is 1,630 mm which is increasing by 4.960 mm each year over a course of 55 years, from 1960 to 2015 (Mondal et al. 2017). The KCC area is not subjected to direct flooding from Rupsha-Bhairab River but the low-lying areas of the western and southern parts of KCC are prone to be flooded during monsoon period due to heavy rain and tidal flooding (Zannat & Islam 2015).

Temporal satellite images from WorldView-2 and Quickbird-2 for four different years viz. 2005, 2010, 2015, and 2020 were used in this study. A digital elevation model (DEM) of the study area with a spatial resolution of 1 m was also collected from Khulna Development Authority (KDA) for the year 2015. Table 1 lists the specifications of satellite data used in this study to detect the LULC and their changes over the study periods.

Some additional datasets were also used for demographic (population) study, rainfall analysis, and soil map preparation for various hydrological soil groups (Table 2).

The collected DEM was used to perform watershed delineation and identify different catchments for the study area using hydro-processing methods in a GIS (ArcGIS 10.5) environment.

Past studies revealed that enormous LULC changes have occurred due to urbanization, which plays a significant role in influencing direct runoff estimate and infiltration processes, and are thus regarded as a fundamental input for the SCS-CN technique (Dunjó et al. 2004; Moniruzzaman et al. 2018, 2020). Cloud-free time-series satellite images from Worldview-2 and QuickBird-2 were chosen for the preparation of LULC layers for 2005, 2010, 2015, and 2020. From Table 1 it is evident that instead of year 2020 the satellite data were collected at the very end of year 2019 to ensure cloud-free quality image and this was considered as the image for 2020. Necessary atmospheric corrections and image enhancements were done as part of the image preprocessing works. A machine learning algorithm i.e. SVM-based supervised classification (Traoré et al. 2019) technique was applied to classify the images, for different years into four land cover categories e.g. vacant land, urban area, vegetation, and water body. The details of the land use classes are given in Table 3.

An ArcGIS SVM toolbox that connects SVM with ArcGIS through a freely available support vector machine classification library ‘LIBSVM’, developed by Chih-Chung Chang and Chih-Jen Lin at National Taiwan University, was used to classify the satellite images into different LULC classes (C. C. Chang & Lin, 2011; Wehmann, 2015). LIBSVM was created to enhance the use of SVM in other scientific fields (Chang & Lin 2011).

The training samples were collected based on the visual inspection of natural color composite (RGB combination 3-2-1) of satellite images. VL=vacant land, UA=urban area, VG=vegetation, and WB=water body; were used as classification categories. Training pixels were sampled at random for each class in each image by drawing several polygons all over the images and then partitioned into training and testing datasets at random (Table 4).

Each image was classified using the training data via a SVM classifier that used a radial basis function (RBF) kernel with 1-vs-1 multiclass strategy. Here, RBF kernel was chosen for classification since this kernel accounts for better performance for remotely sensed data compared to other kernels (linear, polynomial, sigmoid etc.) (Huang et al. 2010; Mountrakis et al. 2011; Shi & Yang 2015; Cardoso-Fernandes et al. 2020). In order to train the SVM model for classification; a regularization parameter, C and a bandwidth parameter of the RBF kernel, must be chosen. C indicates how confident the classifier should be on the training data's accuracy, with higher values suggesting more confidence. Whereas, determines the smoothness of the function approximated with the kernel. Greater degree of smoothing is indicated by higher values which basically refers to the decision surface generated during SVM training and has nothing to do with the output classification's spatial properties (Hsu et al. 2003). A five-fold cross-validated grid search was utilized to select the value of these hyperparameters. The SVM regularization parameter C was varied in 10 stages in log-2 space from −5 to 15, whereas the kernel bandwidth parameter was varied in 10 steps in log-2 space from −15 to 3. For each of the images over the study period (2005, 2010, 2015 and 2020), the parameter set with the highest training accuracy was chosen (Table 5). Then, using the predict function of the toolbox final classified LULC maps have been produced.

The accuracy assessment was performed by taking 100 randomly selected reference points for each class in ArcGIS and ERDAS IMAGINE platform. High-resolution areal images from Google Earth were employed for manual ground truthing. Based on the true and predicted classes of the test data a confusion matrix has been formulated which provides insight for each class about how many samples are correctly classified, misclassified, and for which class the samples were incorrectly classified.

For validation purpose user's accuracy (UA), producer's accuracy (PA), overall accuracy (OA) and Kappa statistics were estimated. Producer's and user's accuracy were found out to evaluate the classification at class level. The samples that were successfully identified in each class are represented by the matrix diagonal, and the OA is determined by dividing the sum of the correctly classified instances by the total number of pixels (Story & Congalton 1986). The PA (PA) is calculated by dividing the number of right pixels in a given class by the total number of pixels in the test set for that class, and it represents the likelihood of correctly classifying a training sample (Congalton 1991). The UA, on the other hand, can be calculated by dividing the number of correct pixels in a given class by the total number of pixels classified in that class. It demonstrates the map's reliability by indicating the likelihood that a pixel that is categorized as belonging to a particular class truly belongs to that class on the ground (Story & Congalton 1986; Congalton 1991). The Kappa statistic is a measure of agreement between the output of the classifier and the reference data (Cohen 1960). Over the course of the study, receiver operating characteristic (ROC) curves were created for each LULC class, which plotted the true-positive rate against the false-positive rate (Cardoso-Fernandes et al. 2020). The area under curve (AUC) score was also computed to compare among the accuracies of different LULC classes.

Furthermore, this study employed a transfer matrix to see the LULC changes in the process of urbanization for different time intervals.

KCC is located in the southwest of Bangladesh. Its north and south sides are located on young Holocene-Recent Alluvium of the Ganges deltaic plain and Ganges tidal plain, respectively. Up to a depth of 300 meters, the region comprises of mostly coarse to very fine sand, silt, and silty clay, and has peaty soil and calcareous and non-calcareous soil at the top (Roy et al. 2005).

For digitization and representation of soil texture types of the study area, the source data were collected from the Soil Resources Development Institute, Bangladesh. According to observation, soil textures in the study area are prominently peaty soil, clay, clay loam, loam, silt, silty clay loam, etc. among which peat, silty clay loam type soil are mostly found around the study unit. Again, the HSG serves as a key input for the SCS-CN model and according to the United States Department of Agriculture (2009), the HSG layer was formulated by reclassifying soils based on their characteristics and potential for surface runoff (U.S. Department of Agriculture 2009). The overall study area exhibits three types of HSGs, e.g. A, B, and C, which cover 62.33%, 32.64%, and 5.03% shares of the region respectively (Figure 2). Table 6 depicts the behavior of these HSGs towards infiltration rate and runoff potential of a given rainfall event. The northern part of the study area mostly exhibits HSG A, which is of low potential for surface runoff. However, the southern part of KCC area consists of HSG C, which has moderate to high runoff potential. A very small portion of HSG D (of high runoff potential) is found on the north-eastern part along the river running beside KCC area.

Based on the types of LULC and varying HSG a set of CN values has been recommended by Cronshey in 1986 (Cronshey 1986). The values for curve number are computed through the CN grid by incorporating LULC classes, HSG, and the standard CN-II value (Table 3) lookup tables (Cronshey 1986; U.S. Department of Agriculture 2009; Subramanya 2013) in the ArcGIS platform using the HEC Geo-HMS extension (Figure 3).

The moisture content of soil at the start of a rainfall–runoff event is referred to as the AMC and it is widely recognized for governing early abstraction and infiltration (Subramanya 2013). There are three types of AMC viz. AMC I, AMC II, and AMC III, which stand for dry, average, and saturated soil condition, respectively, in terms of soil moisture content (Silveira et al. 2000) (Table 7).

In this study, both of the AMC and CN were considered to be in moderate condition. Hence, AMC II- and CN-II-based runoff values were estimated for this research.

where S denotes the potential retention of the area concerned and P is the maximum daily rainfall with the selected -year return period. The rainfall P was calculated using the Bangladesh Meteorological Department's maximum daily rainfall data of years 1987–2017.

It is to be mentioned that model validation could not be performed with the observed runoff data since they were unavailable for the study area. This would not be much of a concern since this study is predominantly focused on relative changes in runoff depth over the study period (Perry & Nawaz 2008; Zare et al. 2016). Therefore, a comparison of these study findings with similar studies done before was outlined as validation for runoff assessment.

where and represent the surface runoff depth (mm) in the initial and final land use scenarios of each stage respectively, denotes the rainfall depth (mm), and are the absolute amount of runoff change, denotes the relative degree of change. If and have a positive value, the land use change at this stage is indicated to lead towards an increased runoff.

Finally, correlation analysis was conducted to assess the association between the change of surface runoff and land use driving factors. Spearman Rank Correlation analysis was used in this study because the data did not satisfy the assumption of normal distribution (Hu et al. 2020). The stronger the correlation co-efficient, the greater the influence on the change of surface runoff, it may be deduced.

Figure 4 depicts the hypsometry of the KCC area. In terms of elevation from the mean sea level (MSL), the highest share was recorded in the range of 1 m to 3 m (59.58% of the total KCC area), which was followed by the elevation intervals of 3 m to 4 m (21.08% of the KCC area), mostly found in areas alongside the Rupsha River. The next highest share and the lowest share were recorded for the hypsometric intervals of −1 m to 1 m (18.41% of the KCC area) and 4 m to 7.61 m (0.93% of the KCC area) respectively.

Data of 30 years (1987–2017) were collected from Bangladesh Meteorological Department, Khulna Station as the source data for the rainfall analysis. The highest daily maximum rainfall was recorded to be 271 mm/day on August 22, 2016. Figure 5 shows the distribution of annual daily maximum rainfall during 1987–2017. Again, the annual maximum daily rainfall for 30 years has been synthesized for trend analysis and the return period was calculated using the Weibull's formula.

Lastly, a logarithmic trendline was fitted between annual daily maximum rainfall and return period to project the rainfall amount of a rainfall event of 100 years (Figure 6). It is estimated as 360 mm/day, which was further used in this study for simulating different runoff scenarios. The value for R² of the trendline was found to be 0.96, which indicates that the trendline fits the data pretty well. This study considered a rainfall event with a return period of 100 years to see the runoff response of an extreme rainfall event given a significant amount of LULC change over the study region.

To investigate classification accuracies, four measures viz. PA, UA, OA, and Kappa were calculated using the confusion matrix method (Kamusoko et al. 2014; Hu et al. 2020).

The classified LULC namely vacant land, urban area, vegetation and water body had overall accuracies of 82.45% (2005), 85.08% (2010), 84.89% (2015), and 86.50% (2020) with Kappa coefficients of 0.80, 0.82, 0.81, and 0.84 respectively (Table 8).

All the Kappa values are close to 0.80, indicating that the LULC classes of the actual and classifier's output images are approximately 80% similar. Again, the ROC curves and AUC scores greater than 0.80 suggest satisfactory goodness of fit for all of the LULC classes for 2005, 2010, 2015 and 2020 (Figure 7). This result also affirms a substantial agreement between the reference data and simulated classified LULC data over the study period (Mukherjee et al. 2009; Nath et al. 2020).

The overall changes of four LULC types viz. vacant land, urban area, vegetation, and water body were analyzed based on satellite imagery for 2005, 2010, 2015, and 2020. Figure 8 illustrates the chronological LULC changes and Table 9 shows the share of LULC change in areas and percentage for the chosen years. In addition, Figure 9 depicts relative LULC changes for different years of the study area. The vacant land showed an overall reduction of 3.67% from 2005 to 2020. The highest share of vacant land was observed to be 23.63% in the year 2015, whereas it was only 7.22% in 2010. Although the change of vacant land was boosted by 16.41% during the period 2010–2015, it showed a declining trend of 11.89% and 8.19% for the years 2005–2010 and 2015–2020, respectively.

The urban area experienced significant growth of 7.96% and 1.87% in the years 2005–2010 and 2015–2020, respectively. However, for the years 2010 and 2015 it almost remained unchanged (14.67 km² and 14.66 km², respectively). Overall, during the years 2005 and 2020, the urban area share was 24.42% and 32.38%, respectively, indicating a growth of 9.82% in the study period. The impact of rapid urbanization, physical development, rural-to-urban migration, etc. made a potential contribution to the growth of urban areas of KCC during 2005–2020 (Rahman et al. 2009; Moniruzzaman et al. 2018; Sarkar et al. 2021).

Vegetation cover among the LULC classes had the largest share in all of the study years except 2015. It is to be mentioned that vegetation cover for the years 2005 and 2010 made up approximately 50% of the total area. In addition, the share of vegetation cover reduced until 2015 and became approximately equal with the share of urban area in 2015 and 2020. Moreover, in the years 2005–2010, 2010–2015, and 2005–2020, it was decreased by 1.59%, 18.8%, and 13.24%, respectively, but was increased by 7.15% in the period 2015–2020.

Conversely, the share of water bodies had an increasing trend until 2015 and escalated significantly from 2005 to 2010 and it has almost doubled. During the years 2005–2020, the share of water bodies increased by 7.11%, and this trend continued for the years 2005–2010 and 2010–2015, rising by 5.53% and 2.41%, respectively. However, between 2015 and 2020, it remained almost the same. Various government attempts to rehabilitate water bodies have resulted in an increase in the overall share of water bodies throughout the study period (Zannat et al. 2020).

The direct surface runoff for the whole study unit was quantified using the SCS-CN technique based on a rainfall event with a 100-year return period (360 mm/day) and different hydrological soil cover combinations. The resultant surface runoff depth ranged from 192 mm/day to 360 mm/day (Figure 10) based on the retention capacity of different surfaces, which are classified into Low (192–200 mm/day), Moderate (200–300 mm/day) and High (300–360 mm/day). Table 10 depicts that areas with low runoff potential were almost the same both in 2005 and 2010 (32.33% and 32.39%, respectively).

These areas predominantly belong to vegetation land cover class with higher potential of rainwater interception, a superior infiltration rate, and a soil texture (peat, sand etc.) that falls into HSG-A. The area shares were decreased in years 2015 and 2020 (17.11% and 24.29%, respectively) and lost area shares became areas with moderate to high runoff potential.

Urban areas are particularly vulnerable (Figures 8 and 10) to surface runoff (as well as possible flooding), as seen by high runoff depth (300–360 mm/day interval) and the areal extent of surface runoff with high potential was found to be 38.37%, 47.86%, 53.99%, and 55.68% in years 2005, 2010, 2015, and 2020, respectively (Table 10). This signifies a gradual increase in surface runoff depth across KCC area over the study period. The areas with high runoff potential correspond more or less to urban area LULC class with primary soil types of HSG C and D (silty clay loam, clay, etc.) with a lower rainwater infiltration rate. Furthermore, the moderate runoff depth interval of 200–300 mm/day has experienced an overall decline in area share by 9.27% during the period 2005–2020 and the lost areas have been replaced by areas with high runoff potential.

The study encompasses an assessment of the variation in urban rainfall–runoff due to LULC changes in KCC from 2005 to 2020. The findings of Mondal and colleaues annotate an increase in rainfall by 4.960 mm/year over a course of 55 years, from 1960 to 2015 (Mondal et al. 2017). The upward trend in precipitation results in inundation of the city, causing distress and sudden halts to urban life (Molla 2019; Roy 2019). The findings of Sarkar and colleagues affirmed that haphazard urban development, encroachment of canals, and the dysfunctionality of the existing sluice gates in KCC propel severe inundation during intense rainfall incidences (Sarkar et al. 2021). The convictions of (Vojtek & Vojteková (2019), Roy et al. (2020), Hounkpè et al. (2019), and Azizi et al. (2021) assert that urban inundations are notably associated with the LULC changes of an area. The results of this study are in line with the mentioned inferences on the correlation between LULC changes and inundation trends. Several national electronic media correspondents have captured the severity of such flooding incidences in Khulna City through digital lenses (Figure 11).

Table 11 illustrates the increment of surface runoff over different periods for a rainfall event of a 100-year return period. The overall study area was divided into 60 catchments and the average runoff depth of the catchments was estimated in this regard. During these study periods, the increment in surface runoff was observed as 7.13 mm, 14.71 mm, and −4.83 mm, respectively, for years 2005–2010, 2010–2015, and 2015–2020. The increment for surface runoff coefficient was observed as 0.02, 0.04, and −0.01 for three different time intervals respectively. Again, the relative degree of runoff change at three stages were 2.46%, 4.96%, and −1.55% respectively. During the overall study period (2005–2020), KCC area experienced an increase of 17.00 mm in average surface runoff depth which justifies the continual upward trend of runoff depth for a given rainfall event.

The association between surface runoff change and land use change was investigated using Spearman rank-order correlation analysis. Table 12 shows the correlation coefficients between and the different land use change rates for 60 catchments at various intervals of the study periods for a rainfall event of 100-year return period.

This result depicts that changes in runoff depth were positively correlated with changes in urban areas and negatively correlated to changes in vegetation cover for all of the study stages. The degree of correlation between and land use variables dropped in the order of VG (−0.746) and UA (0.346) between 2005 and 2010, all of which were statistically significant at a significance level of 1%. VG (−0.821) was the most important factor relating to from 2010 to 2015, followed by UA (0.620) and VL (0.536) at 1% level of significance. Between 2015 and 2020, there were strong associations between and VG (−0.838), VL (0.556), and UA (0.223) when the level of significance was 1%. Throughout the study period, the major driving forces for surface runoff change have been identified as urban area and vegetation cover change.

This study used high resolution WorldView-2 and QuickBird-2 satellite images, which are more likely to produce better LULC classification and hence better runoff estimation of the study area as compared to studies done previously, using Landsat images. Figure 10 justifies the findings of runoff assessment of this study portraying enormous sufferings of the people due to stormwater induced small flood events and waterlogging in and around different parts of Khulna City. The findings of this study are in line with other study findings which are already discussed above in the validation section. Again, this study attempts to explore the major land use drivers responsible for causing runoff variations over the geographical study space. Urban area and vegetation cover are found significantly responsible to generate surface runoff which are similar to studies done by Li et al. (2019) and Hu et al. (2020). According to Zannat & Islam (2015) among many factors increased impervious land in urban areas and decreased vegetation cover cause water congestion during incessant rainfall in Khulna City and the duration of this water logging is often 2–3 days, which also supports the study findings. The prime focus of this study was to see the spatio-temporal variability of direct surface runoff potential and their drivers using CN-II-based runoff values rather than a full-scale hydrological or flood models. Surface runoff calibration can be done separately as separate research, with observed river flow data being used to calibrate CN-II values. Nevertheless, the findings of this study would help the city authorities to plan and manage stormwater related crisis more efficiently by adopting necessary development control measures.

The study region has experienced an increase in urbanization in the past few years. In recent years, urban planning and city administration have worked together to decrease the catastrophic risk that urbanization brings with it. Quantitative study on urbanization's impact on surface runoff is vital for city planning and preventing urban flooding under a given rainfall. This research was able to approximate runoff fluctuations in KCC's main urban region using GIS and remote sensing technologies as well as the SCS-CN model. The findings of this study are summarized below:

• (1)

  From 2005 to 2020, the urban area has experienced a significant growth of 9.82% and a decline of 13.24% in the share of vegetation cover.

• (2)

  Throughout the study period (2005–2020) the relative degree of average runoff depth has increased by 5.44% (17.00 mm) on a given day with a 100-year rainfall event and this trend was more or less consistent for three different time intervals shown in Table 11.

• (3)

  The rise in runoff depth was found positively correlated (p-value<0.01) to the growth of urban impervious areas and negatively correlated (p-value<0.01) to the changes in vegetation land cover class for all three study intervals (2005–2010, 2010–2015, and 2015–2020). Hence, from 2005 to 2020, impervious urbanized area and vegetation cover land use in KCC region were the predominant drivers of surface runoff change.

All relevant data are included in the paper or its Supplementary Information.