The Pinacanauan de Ilagan River, with a catchment area of approximately 3000 km², serves as a drain or outlet for the southeast region of the Cagayan catchment. The Pinacanauan de Ilagan River continues its stream as the Ilagan River and branches out into two channels—the Abuan and Bintacan catchments. Particularly, the Abuan and Bintacan rivers are located at the northeast portion of the Pinacanauan de Ilagan River. The Abuan catchment upstream of the Abuan bridge has an area of 490 km², while the upper Bintacan catchment has an area of about 111 km² (Dingle et al., 2019). Figure 1 shows the location of the study area and the discharge monitoring stations that were used in this study. A water level monitoring station located mid-stream was also used for the calibration of the hydraulic model using the river’s rating curve.

The Abuan River Basin, located about 10 km upstream of the confluence of the Cagayan and Pinacanauan de Ilagan Rivers, has an upper watershed consisting of the Sierra Madre Mountain Range and a lower watershed containing foothills and floodplains at the merging point of the Abuan River and Pinacanauan de Ilagan River. For the whole Cagayan River basin which includes the Ilagan basin, the rainy season usually occurs from May to November, with maximum monthly rainfall typically present in the months of July or August due to the southwest monsoon, and in November due to typhoons. Near Ilagan, the average annual precipitation is about 1489.5 mm in the years 1976–1985 (Japan International Cooperation Agency, 2002) and flood events usually occur in the basin due to abundant rainfall brought about by these typhoons.

This study utilized numerical modeling to determine flood depths and inundated areas for storm events and the possible effects of climate and land cover change on the flood hazard of the area. The data needed for the study were categorized as either spatial data, meteorologic data, or hydraulic data. Table 1 summarizes the input data needed to set-up the model and how they were utilized to meet the study objectives.

The topography used as input in the hydraulic model was the Interferometric Synthetic Aperture Radar (IFSAR) digital elevation model of the National Mapping and Resource Information Authority (NAMRIA), the local mapping agency of the country. This was used to represent the geospatial conditions of the study area particularly to characterize the flow and inundation of the flood within the study site boundary given its 5-m-by-5-meter fine resolution as compared to other satellite-derived DEMS which are usually coarser. On the other hand, the land cover map was used to classify the land use of different areas within the region that serves as guide in assigning Manning’s roughness coefficient (n) corresponding for each land use during the calibration and scenario simulation stages. The daily rainfall from the Philippine Atmospheric Geophysical and Astronomical Services Administration (PAGASA) served as input to simulate the rainfall-runoff to inundation process for chosen storm events while the river discharge data from the Department of Public Works and Highways (DPWH) were used to fine-tune the model before carrying-out the climate change scenario analysis.

The numerical modeling tool HEC-RAS, created by the US Army Corps of Engineers, was used in this study because of its capability to simulate flood development in either one or two dimensions. Moreover, it has a considerably vast amount of documentation, tutorials, and manuals to aid the user in understanding and navigating through the software. Besides the significance of its functionalities for the study, this is also selected among other numerical modeling software because it is widely used in river hydraulic analysis in both research and practical applications in the Philippines. Shown in Figure 2 is the general flow of the study.

The base geometric model of the river was first developed using the provided terrain and land cover maps of the area. Boundary conditions that act as forcings such as precipitation and river discharge were added into the model and unsteady flow analysis was performed. To ensure acceptable model performance and results, calibration and validation was done before carrying out scenario analysis on the effects of climate and land cover change to the flood hazard of the area.

Calibration of the hydraulic model was done to ensure a realistic model behavior (Mostafaie et al., 2018). The main process of calibration was done by varying the Manning’s roughness coefficients incorporated into the area’s land use classification. The simulated discharges were compared against observed data for a historical storm event that has affected the study area namely Tropical Storm Mekkhala (local name “Amang”). During this stage, the roughness coefficients of the base model were further varied until the simulated flows are comparable with the selected observed data based on statistical parameters. For model validation, the resulting flows of the calibrated model were compared against observed data from two other storm events namely Typhoon Rammasun (local name “Glenda”) and Typhoon Kalmaegi (local name “Luis”) which were different storm events from the ones used during calibration. This evaluation of the model performance considering other storm events is necessary to ensure that the model will yield satisfactory results for all other simulations including scenario analyses. Summarized in Table 2 are the storm events used during model calibration and validation and their simulation periods.

To quantify the sufficiency of the model’s performance, a series of efficiency and performance criteria were used such as the Nash-Sutcliffe efficiency (NSE), ratio of the root mean square error (RMSE) to the standard deviation of measured data (RSR), and percent bias (PBIAS). This method of calibration and evaluation was used in previous studies on coupled watershed and flood modeling such as that from (Puno et al., 2018), and it is also a part of the guidelines set for the quantification of accuracy in watershed simulations by (D. N. Moriasi et al., 2007).

The Nash-Sutcliffe efficiency (NSE) determines how well the observed versus the simulated data fits the 1:1 line. The difference between the observed and simulated data are considered acceptable if the NSE is between 0 and 1. Otherwise, the performance of the simulation will be considered unacceptable. More specifically, the acceptable value for NSE should be greater than 0.65 (D. N. Moriasi et al., 2007). The NSE value is calculated using Equation 1. N S E = 1 − ∑ i = 1 n Q i o b s − Q i s i m 2 ∑ i = 1 n Q i o b s − Q m e a n 2 (1)

where: NSE: Nash-Sutcliffe efficiency value of the simulation.

Q i o b s : ith observed depth of the river

Q i s i m : ith simulated depth of the river

Q m e a n : mean value of the observed depths of the river

n: total number of observations.

Root Means Square Error-observations standard deviation ratio (RSR) standardizes the RMSE using the standard deviation of the observations and includes error index statistics. For this parameter, lower RSR and RMSE values indicate better model simulation performance where RSR values of less than 0.60 are considered acceptable (D. N. Moriasi et al., 2007). The RSR is calculated using Equation 2. R S R = R M S S T D E V o b s = ∑ i = 1 n Q i o b s − Q i s i m 2 ∑ i = 1 n Q i o b s − Q m e a n 2 (2)

where: RSR: RMSE-observations standard deviation ratio.

RMSE: Root Mean Square Error.

S T D E V o b s : observations standard deviation

Lastly, the percent bias (PBIAS) evaluates the mean tendency of predicted simulated results to overestimate or underestimate the observed data. A PBIAS of 0 would mean an accurate model. A positive value pertains to underestimation while a negative means overestimation. For this parameter, a percent bias of less than 15% is deemed acceptable (D. N. Moriasi et al., 2007). The percent bias of the simulation is calculated using Equation 3. P B I A S = ∑ i = 1 n Q i o b s − Q i s i m x 100 ∑ i = 1 n Q i o b s (3)

where: PBIAS: Percent bias.

After satisfying the calibration and validation criteria, the calibrated model representing the Pinacanuaman de Ilagan river hydraulics was then used in generating a flood hazard map during Typhoon Haima (local name “Lawin”) with a simulation period of October 11-27, 2016, as the baseline condition.

Using the calibrated river model, the effects of land cover change on the flood hazard parameters such as flood depth and extent in the area were examined using land cover from NAMRIA. Two land cover maps (2015 and 2020) were obtained and used as inputs in the hydraulic model to compare flood depths based on simulation results. Since the baseline condition will use Typhoon Haima that occurred in 2016, the 2015 land cover will serve as the baseline input for this scenario during the analysis of results.

The effect of climate change was also examined by applying the percent change in precipitation for the area based on the CLIRAM tool of PAGASA. The climate projection of the tool is derived from two prominent scenarios outlined by the Intergovernmental Panel on Climate Change (IPCC), the RCP 4.5 (moderate greenhouse gas emissions) and the RCP 8.5 (high-level greenhouse gas emissions). In assessing the corresponding effects of the change in precipitation for the most critical condition, only the projected scenario for the upper limit of the High-Emission Scenario (RCP8.5) was selected. In this study, it was assumed that the upper limit of the High-Emission Scenario (RCP8.5) scenario already simulate the worst flooding circumstance that will warrant risk management. Shown in Table 3 are the projected changes in rainfall in the mid-21st century for Isabela, Philippines.

The percent change in precipitation was applied to the observed rainfall during Typhoon Haima to simulate the scenario of a projected storm with a higher intensity driven by climate change. The resulting flood hazard maps were then compared to the baseline condition for analysis.

The hydraulic model of Pinacanauan de Ilagan River was calibrated to a medium-intensity storm event namely Tropical Storm Mekkhala. Channel roughness coefficients of the Abuan, Bintacan, Ilagan, and Cagayan Rivers were adjusted until the model simulated discharges were comparable to observed discharges and acceptable according to the statistical performance criteria. The river model was then tested for two other storm events (Typhoon Rammasun and Typhoon Kalmaegi) to further validate model results using the calibrated roughness coefficients of the channels. Figure 3 shows the visual comparison between the observed flow from the calibration point at midstream and the model simulated discharges for the calibration and validation runs.

Modeled and observed discharges were then used to calculate the Nash-Sutcliff Efficiency (NSE), Percent Bias (PBIAS), and the RMSE-observations standard deviation ratio (RSR) to further assess model performance. These criteria were satisfied with final channel roughness values of 0.035 for Abuan river, 0.025 for Bintacan, 0.045 for Ilagan, and 0.025 for Cagayan river. Table 4 shows the summary of the statistical parameters that were used to assess the model’s performance for this study.

The calibrated manning’s roughness for the river channels produced model discharges that are identical to the observed discharges as shown in Figure 3. Successful model calibration is further shown by the calculated values of the statistical parameters shown in Table 4. Using the modeled and observed discharges, the NSE of the model was determined to be 0.86 and according to (D. N. Moriasi et al., 2007), the model is already acceptable when the NSE is greater than 0.65. In addition, PBIAS values of less than 15% and RSR values of less than 0.60 also indicate satisfactory model performance that were also exhibited by the model with PBIAS and RSR values of 2.26% and 0.37, respectively after calibration.

For model validation, the same channel roughness coefficients were used for the two other storm events and the same statistical parameters were used to evaluate the sufficiency of the model. The model was able to produce an NSE of 0.70, PBIAS of 6.03%, and RSR of 0.55 that were all indications of satisfactory model performance for simulation considering Typhoon Rammasun, while validation using Typhoon Kalmaegi also produced acceptable results having an NSE, PBIAS, and RSR values of 0.80, 6.03%, and 0.45 respectively. Overall, results show that the calibrated river model has an acceptable performance, with an underestimation of about 2%–6%. Having satisfied the statistical criteria for model performance in modeling the discharges which signifies no systematic errors in simulation (Shahid et al., 2018), the Pinacanauan de Ilagan River model was then used to analyze the flood hazard considering a historical storm (baseline) and the two other scenario analyses.

A high-intensity storm event, Typhoon Haima, was used as baseline condition for the hydraulic simulation. The choice of the event was based on the assumption that high-intensity storms may represent extreme cases of flooding specifically for areas that are located on floodplains or near the banks of the river. Depicted in Figure 4 is the results of hydraulic simulation for the baseline typhoon event.

Based on simulation results, the maximum precipitation that has fallen in the area during the simulation period was 24 mm on 15 October 2016. At around the same day, both the discharges and water level elevations started rising to values of 166.83 m³/s and 27.94 m, respectively. The hydrograph at the monitoring station mid-stream reached its peak discharge on the 20th of October where this date also corresponds to the highest water level elevation that was simulated which was 35.74 m. These peak values of discharge and water level occurred a day after the second highest rainfall during the simulation was recorded which is 14.5 mm.

The rainfall-runoff process in the upstream catchments has introduced flows going into the main river in addition to the precipitation that was transformed to streamflow which have resulted in a peak flow of 8250.99 m³/s as shown. The calibrated roughness of the main channel was determined to be 0.045 which could also have affected the flow. Furthermore, comparing the value of the peak discharge during Typhoon Haima to the other storm events used for calibration and validation shown in Figure 3, it can be seen that it is indeed a storm of higher intensity and is fit for the objective of determining the resulting flood hazard in the study area.

Considering the same storm event and the 2015 land cover as input to represent the surface roughness that affect the flow of water in the channel and floodplains as baseline, the inundated areas within the study site boundaries were analyzed. Three observation points within the study site boundary were selected and the flood depths were determined based on hydraulic simulation results. A summary of the selected observation points and description is provided in Table 5.

Figure 5 shows the flood inundation within the study area. The residential area (O1) which is adjacent to the river shows flood rising to a depth of about 1.25 m during the storm event. Aside from being close to the main tributary, this residential area is also located near the confluence of a smaller channel and the main river that could have added to the overflow volume reaching the banks. Meanwhile, the second observation point (O2) located at midstream shows that flood depths reach about 1.49 m. This point was chosen to be analyzed in the study since it is located near a bridge that serves as a route for evacuation in the event of natural disasters such typhoons and floods. Higher flood depths are observed at midstream as compared to O1 which can be attributed to the additional flood plain area southeast of this point as can be seen in the flood hazard map. The most upstream observation point (O3) was chosen since it is located in an agricultural area which is one of the primary livelihood sources in the area. Simulation results show that flood depths reach to about 1.87 m at this point which can be attributed to it being located in a low elevation area as shown in Figure 1. Moreover, this point is near the river confluence thus making it more susceptible to flooding once the river exceeds its capacity especially during typhoon events.

To determine the changes in cover within the study area and its probable effects on the flood hazard, the 2015 and 2020 Land Cover Maps from NAMRIA are shown in Figure 6 and a summary for each cover is provided in Table 6. A simplistic approach of merging some covers such as Annual Crop and Perennial Crop into a single representative cover that is Agricultural was done for better visualization on the map and simpler analysis. From the figure shown, it can be observed that the concentration of built-up areas located near the river and downstream have decreased from the year 2015–2020. In terms of location, built-up areas can also be seen to have spread with greater proximity from the river in 2020 as compared to the historical 2015 land cover.

The forest cover upstream particularly on the side of the Bintacan River has shown to have increased by 15.48% from 2015 to 2020 while the grasslands and brushlands as well as agricultural areas have decreased by 17.65% and 19.72% respectively during the same period. The 2020 land cover also shows an increase in barren lands which are located primarily along the river. This increase in barren lands can be viewed as a result of sediments being transported by the natural river processes through time. Moreover, small creeks can also be seen in the land cover of the year 2020 that are not present in the 2015 available map. Because the change in land cover can cause an increase in the flood risk of an area (Ma et al., 2024), the resulting flood hazard in terms of flood depth considering these changes in cover in the study site was quantified in this study.

The calibrated model was again used to simulate the Typhoon Haima storm event now using the 2020 land cover to determine the flooding in the area. Table 7 contains summary of the flood depths at the observation points comparing the two different land covers used in the study.

From the results of using the 2020 land cover, the flood depth at the residential area downstream reached 1.30 m as compared to 1.25 m depth of flood when using the 2015 cover. The 4% increase in depth can be attributed to the presence of barren land downstream as shown in Figure 6. On the other hand, a 2.01% increase in maximum flood depth was observed midstream at the second observation point where the flood depth increased from 1.49 to 1.52 m. Based on the land cover maps, there is also an increase in barren land near this observation point. Moreover, it can be noticed that some channels became narrower and thus could have also decreased the volume capacity. The largest increase in flood depth was determined in the upstream observation point where the depth increased from 1.87 to 2.11 m and a change in cover from agricultural to grasslands can be seen on the available maps.

The effects of climate change on the flood inundation within the study area were also examined by multiplying the projection factor from the Climate Risk Analysis Matrix (CLIRAM) of the local weather agency to the time-series rainfall of Typhoon Haima while maintaining the 2015 land cover as input. Depicted in Figure 7 is the flood hazard map of the site considering the effects of climate change.

An increase of 0.23 m in the depth of flood for the residential area downstream was calculated while the midstream point increased its depth by 0.11 m considering the effect of climate change. Flood inundated areas near the upstream observation point can be seen to have increased compared to the baseline condition. Flood can also be seen to have propagated to more areas given the relatively level topography near this observation point. Aside from more areas being inundated, calculations show that maximum depths have reached 2.24 m which is about a 19.79% greater compared to the baseline peak flood depth. In general, the flood depths at the three observation points increased when factoring in an increase in precipitation over the area since more rainfall could be transformed to runoff and eventually ended up in the streams. Furthermore, since the capacity of the channels remain the same while the volume of water going into the streams is greater as a result of increased precipitation, higher flood depths are expected especially in areas near the river and at lower lying elevations.

Quantifying and describing flood hazards considering the effects of climate change and urbanization is vital for proactive disaster risk management strategies. This study utilized numerical modeling to calculate the flood depths and identify inundated areas resulting from a historical storm event and perform scenario analyses considering the changes in precipitation and land cover within the study area. Performing roughness calibration for the channels has produced more comparable discharges when compared to the observed values. The graphical comparisons in 3.1 showed good agreement between simulated and measured discharges while calculated statistical parameters are all within acceptability according to the study of Moriasi et al. (2015). In relation to the main objective of this study, the roughness parameter calibration is significant in predicting inundation extent and flow as highlighted by (Pappenberger et al., 2005).

After achieving acceptable results based on statistical criteria, modeled peak discharges are still underestimated which could have been caused by the use of riverbed elevations from the digital topography rather than using this DEM in conjunction to actual river bathymetry which is a method that has proven to significantly improve model performance (Muthusamy et al., 2021). But aside from this underestimation of the peak flows, the model was able to capture the general trend in streamflow for the channels in the study site that is important in describing the flood hazard in the area. In addition, the final roughness values for the channels within the study area ranged from 0.025 to 0.045 that is within the range of natural channels with stones and some weeds based on (Chow, 1985). This may also be related to the barren lands and grasslands along channels depicted in the land covers in 3.4 while this can still be validated further by integrating details from actual field surveys.

Results of the baseline hydraulic simulation revealed that the topography of an area plays a critical role in conveying streamflow resulting from a rainfall event. The timing of the peak discharge calculated at mid-stream came sometime after the peak precipitation fell into the area which could have been affected by the travel time of the stream flow coming from the upstream tributaries of Abuan, Bintacan, and Ilagan. This can also be described from the idea of higher-to-lower energy direction of free surface flow of open channels. In addition, the roughness in the main channel could also have caused the delay in the arrival of the flood wave (Asselman et al., 2022). This observation of river flow and water level is however limited to the calibration point located at midstream and incorporating additional observation points along the stream can help describe more the routing behavior of the river.

The generated flood map shown in 3.3 for the storm event also shows that existing stream networks can also aid in predicting areas that are likely to be inundated. The areas where the tributaries from the Abuan, Bintacan, and Ilagan Rivers form a confluence is where prominent swelling can be observed. This is primarily due to the flows coming from these different rivers combining in the main channel, but considering the same capacity, leads to overflow to the banks. In fact, according to the study of Osawa (2024), high flood frequencies are often experienced by areas within 1 km distance from a river confluence. Aside from this, the topography also significantly affects what areas are to be inundated as studies have also pointed out that topography-related factors highly influence flood occurrence (Zhao et al., 2024). Flood inundation maps from models were traditionally validated by comparing it to observed inundation areas. There had already been attempts on the use of satellite videos (Masafu and Williams, 2024) and crowd-sourced data (Songchon et al., 2023) to validate 2D flood simulation results. With reference to the scope of the study, since there is no available satellite image during the period of storm events simulated, the generated flood maps were not validated to check the flood extents spatially which can either underestimate or overestimate the total inundated areas of the study site. However, these maps can be further validated to improve their accuracy but are still deemed acceptable after completing satisfactory streamflow calibration.

Findings from the flood maps generated considering the change in cover show the effects of the land cover-based surface roughness on the resulting flood depths in the area. The presence of barren lands along the channel as possibly an effect of the natural sediment transport at segments of the river often tend to decrease the volume capacity and thus results to river overflow and eventually flooding in the area. In addition, the change in land cover also has affected channels by making them narrower as shown in 3.4 which can have an effect on the river capacity. In relation to the change in land cover, the roughness of grasslands is lower compared to agriculture cover (Asante and Bandaragoda, 2007; Huang, 2009) that may lead to more water flowing from the surface into the stream resulting to overflow. Furthermore, flood peaks and stream flows also are affected by the change in land cover (Gao et al., 2018; Sanyal et al., 2014) that directly contributes to induced flood risks.

Using the local weather agency’s climate projection tool, it is estimated that there will be an increase in precipitation in the study area. This increase in precipitation is expected to exacerbate the flooding in the area and has shown that climate change has more pronounced effects on the flood risk of the study site as compared to the effect of land cover change. Although this study utilized only the RCP 8.5 (High emission) scenario to represent climate change effects, the resulting flood depths are idealized to be the most conservative that can aid future flood mitigation initiatives as also shown in the results of (Sarchani and Tsanis, 2024) where both RCP 4.5 and RCP 8.5 were compared in terms of flood depth and extent.

The model was able to calculate flood depths within its domain, however the flood inundation results are based on the primary influence of topography in the study area. This would imply that any development or alteration in the topography after its acquisition are not reflected. Furthermore, the digital topography (IFSAR) that was used as input cannot capture accurately the riverbed elevation since it usually reflects water surface and thus could affect the accuracy of the model calculated flood depths. Model results could further be improved by combining the digital topographic data with river bathymetry from field surveys as demonstrated by the study of Karim et al. (2013). The layout and impacts of existing drainage systems were also not factored into the hydraulic model that was used in this study. Moreover, the applied rainfall in the area was taken conservatively in hydraulic modeling and was not subjected to losses as it would have been accounted in hydrologic modeling. This may lead to overestimation of the runoff volume since hydrologic losses such as canopy interception, surface depressions, and losses due to infiltration play an important role on the rainfall-runoff process. The spatial accuracy of the model-generated flood maps can further be improved by comparing them to observed flood extents during the storm events simulated through available images from the local disaster management office or on-the-ground interviews from residents within the study areas.

The results of this study can help policymakers in updating the disaster risk reduction and management strategies of local government units. By factoring in the effects of their future developments as well as the effects of climate change, the risks associated with natural hazards can be reduced by limiting the exposure and managing vulnerabilities found within the community. In addition, the methodology that was used in this study can also serve as a framework in understanding behaviors of local systems that are subject to changing landscapes and varying hydraulic and meteorologic forcings for better management and planning.