ANALYZING THE IMPACTS OF LAND COVER CHANGE TO THE HYDROLOGIC AND HYDRAULIC BEHAVIOURS OF THE PHILIPPINES' THIRD LARGEST RIVER BASIN

Changes in land cover can have negative impacts on the hydrological and hydraulic processes in river basins and watersheds such as increase in surface runoff and peak flows, and greater incidence, risk and vulnerability of flooding. In this study, the impacts of landcover changes to the hydrologic and hydraulic behaviours of the Agusan River Basin (ARB), the third largest river basin in the Philippines, was analysed using an integrated approach involving Remote Sensing (RS), Geographic Information System (GIS), and hydrologic and hydraulic models. Different land-cover classes in the ARB for the years 1995 and 2017 were mapped using Landsat 5 TM and Landsat 8 OLI images. Using a post-classification change detection approach, changes in land-cover were then determined. The impacts of these changes in land-cover to the to the basin discharge were then estimated using a calibrated hydrologic model based on the Hydrologic Engineering Center Hydrologic Modeling System (HEC-HMS) under different extreme rainfall conditions. The impact of the changes in land-cover to flood depth and extent was also determined using a hydraulic model based on the HEC-RAS (River Analysis System). Land cover classification results revealed that the ARB is 67.7% forest in 1995 but have decreased to 62.8% in 2017. Agricultural areas in the basin were also found to have increased from 12.2% to 15.5% in the same period. Other notable land cover changes detected include the increase in built-up lands and range lands, and decrease in barren lands. HEC HMS and HEC RAS model simulation results showed that there was an increase in discharge, flood depth, and flood extents between 1995 and 2017, implying that that the detected changes in land cover have negative impacts to hydrologic and hydraulic behaviours of the ARB. * Corresponding author


INTRODUCTION
Many studies have shown that changes in land cover due to anthropogenic activities have negative impacts on the hydrological and hydraulic processes in river basins and watersheds (Matheussen et al., 2009, Zope et al., 2016, Koneti et al., 2018. Human-induced modifications that alter the land cover characteristics such as deforestation, urbanization, and croplands expansions can change how a river basin responds to rainfall (Isik et al., 2008), and they have been found to increase surface runoff and peak flows (Koneti et al., 2018, Guzha et al., 2018. In most cases, these changes led to greater incidence, risk and vulnerability of flooding (Bronstert et al., 2002, Gao et al., 2013, Tarigan, 2016, Zope et al., 2016, Liu and Shi, 2017, which could bring catastrophic damages to human lives and properties (Tran et al., 2010). Hence, quantifying anthropogenic activities, and assessing and understanding its impacts to the hydrologic and hydraulic behaviours of river basins have become very important nowadays not only for efficient water resources planning and management (Koneti et al., 2018), but also for effective flood mitigation and disaster management (Petchprayoon et al., 2010, Zope et al., 2016. Studying the impacts of land cover change on flood behaviour is considered a complex and time-consuming process because the factors that determine river flow and flood intensity, e.g., land cover, vary both spatially and temporally (Petchprayoon et al., 2010). However, such problems have been overcome through the application of Remote Sensing (RS) and Geographic Information System (GIS) technologies together with hydrologic and hydraulic models. RS" capability to provide the most fundamental information describing the nature and extent of land cover especially over large areas in an accurate and timely manner, together with GIS as an efficient tool for spatial data analysis, have been proven to be effective in characterizing the patterns, intensity and dynamics of land cover change (Petchprayoon et al., 2010, Santillan et al., 2011, Thakkar et al., 2017. On the other hand, hydrologic models (also called rainfall-runoff models) can be used to simulate the hydrology of the river basins using rainfall and land-cover information as important inputs (Nie et al., 2011, Santillan et al., 2011. Using this model, a river basin"s response to rainfall in the form of surface runoff and peak flows can be estimated based on different landcover conditions, and the model outputs (e.g., flow hydrographs) are compared to determine the impact of land cover change (Santillan et al., 2011). The flow hydrographs generated by hydrologic models for certain land cover conditions can also be used as inputs into hydraulic models to dynamically simulate the movement of water along the rivers and in the floodplains, thereby generating spatially-distributed flood depths and extents (Santillan et al., 2016). The outputs of the hydraulic model for different land cover conditions can then be subjected to spatial analysis using GIS to ascertain any changes in flood depth and extent, where the detected changes, if there are any, may be attributed to changes in land cover (Zope et al., 2016).
The Agusan River Basin (ARB) is the third largest river basin in the Philippines (Figure 1), with a drainage area of approximately 10,921 km 2 , and a main tributary (Agusan River) measuring 350 km long (ADB, 2008). The basin is of great importance due to its socio-economic contributions that include abundant water resources, logging, agriculture, and mining. An essential feature of the ARB is the Agusan Marsh, an extensive floodplain in the basin where rivers, creeks and tributaries converge and drain northward to the Agusan River and into Butuan Bay (Primavera and Tumanda Jr., 2008). The marsh has long been considered to acts as storage for rain water and reduces the downstream flow of flood water into population centers situated within the Agusan River (Ramsar, 1999), particularly Butuan City, a highlypopulated and urbanized city located in the downstream of the basin. However, recent flooding incidents caused by extreme rainfall brought about by tropical storms, low pressure systems, and tail end of a cold front that resulted to fatalities, agricultural and infrastructure damages around the marsh and in populated centers (NDRRMC, 2014, NDRRMC, 2017 seem to suggest that the marsh is no longer functioning as expected. In fact, high siltation in the marsh due to the deforestation and other anthropogenic activities, particularly mining, in the watersheds upstream of the marsh is a problem that remains to be solved through the years (Ramsar, 1999). Land-use conversions for agricultural plantation development, particularly rice, oil palm and banana plantation developments, have also encroached some of the river basin"s natural biodiversity and structure (Varela et al., 2013).
With the increasing presence of anthropogenic activities and natural changes occurring within the ARB, it is important to quantify how such activities have modified the basin"s land-cover. However, there is a current lack of studies and information on the spatial extent and dynamics of land-cover in the ARB through the years. Therefore it cannot be ascertained how such changes have impacted the hydrologic and hydraulic behaviours of the basin. Undertaking such studies is of tremendous importance for purposes of assessing the current status of basin, and for understanding the impacts of natural changes and humaninduced modifications to the incidence of flooding.
The primary aim of this study is to analyse of the impacts of land cover changes to the hydrologic and hydraulic behaviours of the ARB using an integrated approach involving RS, GIS, and hydrologic and hydraulic models. Using the years 1995-2017 as the period of analysis, the study specifically aims to (i) map the different land cover classes in the ARB for 1995 and 2017 using Landsat images, (ii) detect the changes in land-cover between these years, (iii) simulate the impacts of the changes in landcover to the basin runoff and discharges using a calibrated hydrologic model based on HEC HMS under different extreme rainfall conditions, and (iv) simulate the impact of the changes in land-cover to flood depth and extent using a hydraulic model based on the HEC-RAS (River Analysis System) under different extreme rainfall conditions. The year 1995 was considered the baseline data for the analysis as it can represent the state of the ARB prior to the implementation of many government laws and programs such as the Mining Act of 1995 and the National Greening Program in 2011.

The Study Area
The analysis is particularly focused on how much change has occurred in the portions of the ARB that are upstream of Butuan City (specifically in the upstream of "Dankias Station" as shown in Figure 1), and how much discharge is generated in the upstream that flows in the Agusan River towards Butuan City and Butuan Bay for the two land cover conditions and for different extreme rainfall conditions. For this study, the boundary of the ARB was delineated from a 10-m DEM using watershed delineation tools available in ArcGIS 10.2 software. This DEM, which was provided by the DREAM Program of the University of the Philippines-Diliman, was generated from RADARSAT-2 data gathered during 2012-2013 (Ramirez, 2014). The delineated boundary of the ARB with its outlet located in Butuan Bay has an area of 12,079 km 2 which is larger than the reported drainage area of 10,921 km 2 . On the other hand, the area of the ARB upstream of Dankias Station is approximately 11,523 km 2 .

Landsat images:
The land-cover maps for 1995 and 2017 with 30-m spatial resolutions were generated using a total of 23 Landsat images under path 112 and rows 54 and 55 of the Worldwide Reference System (WRS) where the ARB is entirely contained. Processing such number of images was necessary in order address the problem of cloud cover in all available images of the study area, thereby allowing to substitute or composite areas having cloud cover (Toure et al., 2018). All images were downloaded as Standard Terrain Correction L1TP products from the United States Geological Survey (USGS) Earth Explorer website (http://earthexplorer.usgs.gov/). Of these, 13 images acquired by Landsat 5 Thematic Mapper (TM) were utilized for deriving the 1995 land-cover map, while 10 images acquired by Landsat 8 Operational Land Imager (OLI) were utilized to generate the 2017 land cover map.
For the 1995 land cover mapping, two Landsat 5 TM images (in rows 54 and 55) acquired on April 12, 1995 have the least cloud cover and were utilized as the main source for land-cover mapping of the entire ARB. Missing land cover information in these images were substituted using information derived from Landsat 5 TM images acquired on May 30, October 5, October 12, 1995. Some missing data portions that cannot be substituted using row 54 and 55 1995 images were supplemented using images acquired on November 3, 1994 and September 21, 1996. In a similar manner, the main source of land cover information for 2017 was derived from two Landsat 8 OLI images acquired on August 30 (row 54) and July 13 (row 55), 2017. Missing data was supplemented using images acquired on March 23, April 8, April 24, May 26, June 27, July 13, July 29, and August 30, all in the year 2017.
The geometric accuracies and co-registrations of the images were checked using at least 9 ground control points for each image and using the path 112, row 54 image acquired on April 12, 1995 as the base image. All images were found to be consistently coregistered to the base image, with global Root Mean Square Errors (RMSEs) less than half of a pixel. All the images were then subjected to radiometric calibration (to generate at-sensor radiance images) and atmospheric correction (using dark object subtraction to generate surface reflectance images). Normalized Difference Vegetation Index (NDVI) images were also derived from the reflectance images. Finally, for each imagery date, the reflectance images, NDVI together with a 30-m resolution Shuttle Radar Topography Mission (SRTM) Digital Elevation Model (DEM) were layer-stacked into a single file that will be used as input for image classification. Overall, a total of 23 layer-stacked datasets were produced. Clouds and cloud shadows were manually digitized from images, and from it a mask band was prepared to exclude all clouds and cloud shadows during the classification stage. All image processing and subsequent analyses were performed using ENVI 5 software.

Image classification:
The Maximum Likelihood (ML) classifier was used to independently classify each of the layerstacked image datasets. The goal of the classification is (i.) to produce land cover maps for assessment of land cover change in the ARB and as sources for land cover-related parameterization of the HEC HMS and HEC RAS models. A classification scheme consisting of seven (7) level 1 classes (agriculture, barren land, built-up land, forest, rangeland, water, and wetland) based on the land cover and land use classification system of Anderson et al (1976) was initially used. Then, each level 1 class were expanded to into several sub-classes based on their presence in the ARB except for barren land, built-up land, rangeland, and forest (Table 1). This resulted to a total of 16 classes for Level 2 classification.
A minimum of 1000 pixels per class were selected for classifier training except for the "Sago palm" class where only 333 pixels were selected for the 1995 classifications, and 504 pixels for the 2017 classifications, due to their limited distribution in the ARB. These Regions of Interests (ROI) were collected through visual interpretation of the images which was made possible by the familiarity of the interpreter with the study area, and with the help of high resolution satellite images in Google Earth. Various RGB band combinations were displayed to facilitate image interpretation and accurate collection of training ROIs for each class. Prior to classification, histograms of training ROIs for each class were checked for normality to ensure that the basic assumption of the ML is satisfied.

Level 1 Class
Description (after Anderson et al., 1976)  Only portions of the ARB that are cloud free were included for the classification. Each classification result was visually examined, and post-classification refinements were applied to reduce classification errors caused by the similarities in spectral responses of certain classes (Yuan et al., 2005) such as barren land (barren areas, exposed river beds and beaches) and built-up land (i.e., barren land pixels that were classified as built-up were changed as barren land). For each covered year, the individual results of the post-classification refinements were then mosaicked to create a Level 1 (7 classes), and Level 2 (16 classes) land cover maps. The Level 1 land cover maps were generated by combining the appropriate classes in the Level 2 land cover maps.

Classification Accuracy Assessment:
The accuracies of the 1995 and 2017 land cover maps were determined using the procedures suggested by Congalton (1991) in assessing the accuracy of land-cover maps with more than 12 classes. An independent sample of 100 pixels was randomly selected for each class in the Level 2 land cover maps (total of 1600 pixels). The actual classes of these pixels were then determined by overlaying in high resolution Google Earth images, and through visual interpretation of the Landsat images. The comparisons between the actual and classified Level 2 land-cover classes for each random point were then summarized using confusion or error matrices. The values of Overall Classification Accuracy, Producer"s Accuracy and User"s Accuracy were then computed from the error matrices.

Land cover change detection:
The 1995 and 2017 landcover maps were then subjected to post-classification comparison change detection analysis to identify the location and extent of land-cover change in the ARB. This approach can provide ""from-to"" change information and the kind of transformations that have occurred can be easily calculated and mapped (Yuan et al., 2005).

Hydrologic Modeling using HEC HMS
A hydrologic model based on the Soil Conservation Service -Curve Number (SCS-CN) approach was used to integrate the land-cover information derived from the analysis of the Landsat images into a hydrologic model as a means of estimating the impacts of the differences in land cover conditions to the hydrological processes of the watershed (Santillan et al., 2011), specifically in the generation of runoff under different extreme rainfall scenarios. The goal of hydrologic modeling is to estimate the amount of runoff generated in the headwaters of Agusan River, especially in the portions upstream of Dankias Station, and to determine the amount of discharge that flows into the portion of the river that traverses Butuan City. The SCS-CN model computes runoff (Q) from precipitation or rainfall (P) at the sub-basin level using the following equations ( In the above equations, Ia (initial abstraction) includes short-term losses due to evaporation, interception, surface detention, and infiltration while S stands for potential maximum retention which characterizes the sub-basin"s direct runoff potential (Ponce and Hawkins, 1996). S is directly related to land cover and soil infiltration through the CN or Curve Number which ranges from 0-100, depending on the sub-basin"s antecedent moisture condition (Santillan et al., 2011). In a sub-basin, CN values can vary depending on the number of combinations of land-cover and hydrologic soil group (HSG). Such combination is called a hydrologic soil-cover complex (HSSC), and for each HSSC, a CN value is assigned (NRCS, 2004). For hydrologic modelling purposes, the CN value of a heterogeneous sub-basin can be estimated through area-based weighted averaging of its CN values. . These values were based on Rainfall-Intensity Duration Frequency (RIDF) generated by the Philippine Atmospheric Geophysical and Astronomical Service Administration (PAGASA) based historical rainfall data recorded at Butuan City. The results of the simulations for each year were then analysed and checked for differences, and correlated with the changes in land-cover in the ARB. The simulated discharge hydrographs at Dankias Stations were also used as input into a two-dimensional (2D) HEC RAS model of Butuan City.

Hydraulic Modeling using HEC RAS
The 2D

Land-cover classification and accuracy
The Level 1 land cover classification results for 1995 and 2017 ( Figure 2) have overall accuracies of 93.81% and 93.94%, respectively (Table 3). User"s and Producer"s Accuracies for agriculture, forest, water, and wetlands classes for both classifications were consistently high, ranging from 90% to 99%. Inconsistent PA and UA values were found for barren land, builtup land, and rangeland classes. In the 1995 classification, the PA and UA of built-up land class are high at 94.9% and 93%, respectively; for the 2017 classification, the PA is still high at 93.48% but the UA decreased to 86% implying higher commission errors than that of 1995 classification. For the barren land class, the 1995 PA and UA values significantly increased from 85% to more than 91%. PA and UA values for the rangeland class were relatively low for both classifications except for the 91% UA in the 1995 classification. The accuracy of change detection that can be expected from the two land cover maps was at 88.13% (Table 3) which was obtained by multiplying the individual classification accuracies (Yuan et al.,2005). Relatively higher change detection accuracy rates are found for agriculture and wetland. High commission and omission errors in change detection ranging from 20% to 25% are found for barren land, built-up land, and rangeland. All these imply that some of the detected changes are not actual changes but misclassification errors.

Land-cover change in the Whole ARB
The individual class areas and change statistics for 1995 and 2017 are summarized in Table 4. The 1995 classification shows that 67.7% of the ARB is covered by forest while agricultural areas are at 12.3%; barren land, rangeland and wetland follow at 7.7%, 6.6%, and 4.7%, respectively. Water and built-up land occupies the least area of the ARB at 0.8% and 0.2%, respectively. An almost similar pattern of land cover dominance can be observed for its 2017 condition except that areas classified as rangeland has more area covered (12.2%) than barren land (2.8%), making it the third most dominant land cover in 2017 next to forest (62.8%) and agriculture (15.5%).

Land-cover change in the ARB Upstream of Dankias Station
Since the focus of this study is on the analysis of land cover change impacts to the hydrology and hydraulics in the downstream portion of the ARB, particularly in Butuan City, it would be more relevant to present and discuss the results of the land cover change detection and analysis that were conducted for the portions of the ARB that are upstream of Dankias Station. The individual class areas and change statistics for 1995 and 2017 for these portions are summarized in Table 5. The land cover change matrix is shown in Table 6 where the unchanged pixels are located along the diagonal of the matrix. Just like in the whole ARB, forest and agriculture are the dominant land cover classes in this portion of the ARB for both years. From 1995 to 2017, agricultural areas increased to 1,663 km 2 , with the amount of change equivalent to 3.2% of the total area under consideration. Major contributors in this increase are the agricultural conversions of forest (448 km 2 ), rangeland (172 km 2 ), and rangeland (138 km 2 ). On other hand, barren land decreased to 331 km 2 in 2017, majority of which was due to conversion to forest, agriculture and rangeland. The detected change in agriculture can be said to reflect the actual conversion of lands for agricultural purposes due to high classification accuracies for this class in the 1995 and 2017 classifications. However, the same cannot be said for the barren lands because some of the detected decrease in area can be attributed to the low classification accuracies in 1995 classification.
The built-up land in 1995 significantly increased by 170 km 2 in 2017. It can be said that built-up land area statistics for 1995 is reliable due to higher classification accuracies compared to that of 2017. For both classifications, there is also high confidence that all built-up land were detected. However, overestimation in built-up land in 2017 is undeniable due to the 14% commission errors. While increase in built-up land through the years in the ARB is a normal process due to increasing population and economic development, the majority of the detected changes in built-up land can be more explained by misclassification errors than true changes such as in the detected changes of built-up to agriculture and forest.
Forest cover in this portion of the ARB decreased its area by 588 km 2 in 2017, which is equivalent to 4.9% of this portion of the ARB. This reduction in forest cover appears have been due, in majority, to conversion of 1995 forest cover into rangeland (939 km 2 ), agriculture (448 km 2 ), and barren land (175 km 2 ). While large tract of lands have been deforested, 229 km 2 of agricultural land, 381 km 2 of barren land, and 281 km 2 of rangeland have been converted into forest. These detected changes may be reflective of the National Greening Program which the Philippine government initiated from 2011 to 2016.
Overall (by summing the non-diagonal values of the land cover change matrix), changes in this portion of the ARB accounts to 3,676 km 2 (31.9%) while unchanged areas was computed to be 7,847 km 2 (68.1%).

Flood Simulation Results
The flood maps and statistics generated by the HEC RAS models for Butuan City under three extreme rainfall scenarios (Figure 4) revealed flood depths and extents for 2017 that are generally larger than in 1995. This can be attributed to larger volume of discharge generated in 2017 than in 1995. However, the differences appear to be more pronounced for 5-year rainfall scenario (29.3%), and of little difference for the 25-year and 100-year rainfall scenarios (Table 8). This may imply that the impact of land cover changes in the ARB to flood depth and extent can be of little difference as the rainfall scenarios become more extreme. However, many other factors still need to be considered for this result to be conclusive.   Table 8. Area of simulated flood extents for Butuan City (in km 2 ).

CONCLUSIONS
In this study, the impacts of land cover change to the hydrologic and hydraulic behaviours of the ARB was analysed using an integrated approach involving RS, GIS, and hydrologic and hydraulic models. From the analysis of Landsat images, land cover maps of the ARB for the year 1995 and 2017 were generated with 93% overall classification accuracies. These maps revealed that ARB is 67.7% forest in 1995 but have decreased to 62.8% in 2017. Agricultural areas in the basin were also found to have increased from 12.2% to 15.5% in the same period. Other notable land cover changes detected include the increase in builtup lands and range lands, and decrease in barren lands. The land cover maps were used inputs into HEC HMS and HEC RAS models simulate the impacts of the different land cover conditions to total discharge, flood depth and extents. Simulation results showed that there was an increase in discharge, flood depth, and flood extents between 1995 and 2017, implying that that the detected changes in land cover have negative impacts to hydrologic and hydraulic behaviours of the ARB.