Watershed Modelling of the Mindanao River Basin in the Philippines Using the SWAT for Water Resource Management

This study aims to simulate the watershed of the Mindanao River Basin (MRB) to enhance water resource management for potential hydropower applications to meet the power demand in Mindanao with an average growth of 3.8% annually. The soil and water assessment tool (SWAT) model was used with inputs for geospatial datasets and weather records at four meteorological stations from DOST-PAGASA. To overcome the lack of precipitation data in the MRB, the precipitation records were investigated by comparing the records with the global gridded precipitation datasets from the NCDC-CPC and the GPCC. Then, the SWAT simulated discharges with the three precipitation data were calibrated with river discharge records at three stations in the Nituan, Libungan and Pulangi rivers. Due to limited records for the river discharges, the model results were, then, validated using the proxy basin principle along the same rivers in the Nituan, Libungan, and Pulangi areas. The R values from the validation are 0.61, 0.50 and 0.33, respectively, with the DOSTPAGASA precipitation; 0.64, 0.46 and 0.40, respectively, with the NCDC-CPC precipitation; and 0.57, 0.48 and 0.21, respectively, with the GPCC precipitation. The relatively low model performances in Libungan and Pulangi rivers are mainly due to the lack of datasets on the dam and water withdrawal in the MRB. Therefore, this study also addresses the issue of data quality for precipitation and data scarcity for river discharge, dam, and water withdrawal for water resource management in the MRB and show how to overcome the data quality and scarcity.


Introduction
Among the developing countries, the Philippines faces a considerable challenge regarding development due to the continuous increase in electricity demands, with an annual average rate of increase of 4.3% [1]. The power demand of the Mindanao island group in the Philippines has increased by 3.8% annually over the past decades [2]. In April 2017, the maximum power peak demand in Mindanao reached approximately 1,696 MW [3]. However, the Mindanao water resources contributed 38%, or 1,947 GWh, of the gross power generation from hydropower in June 2017 [3]. Regardless of the current contribution of water resources to renewable energy, the power demand continues to outpace the supply. Thus, to address this emerging problem, assessment for a potential source of sustainable renewable energy is needed. The purpose of this study is to enhance water resource management for hydropower application in Mindanao to improve the electrification situation and support the implementation of the Renewable Energy Act of the Philippines. Thus, this study carefully investigates the precipitation inputs in watershed modelling by comparing the observational data from the DOST-PAGASA with the global gridded precipitation datasets from the NCDC-CPC and GPCC. In addition, water management can include agricultural water footprint analysis by considering increasingly complex indicators, such as multiple climate variables, soil characteristics, and crop properties, to simulate the water cycle [13]. Therefore, weather-related events must be understood in terms of their possible effects on the watershed since hydropower generation is driven by river discharges. Moreover, SWAT considers other climate variables such as temperature, humidity and solar radiation, while soil characteristics and crop properties are some of the inputs to the hydrologic response unit (HRU) of SWAT [14].
Therefore, the objectives of this study are to evaluate the observed rainfall data with gridded global precipitation datasets to overcome the present data scarcity and to simulate the river discharges of the MRB for water resource management and future hydropower development. The literature review on watershed modelling with SWAT is presented in Section 2. The study area, material and method are presented in Section 3, 4 and 5, respectively. The results and discussion are described in Section 6 and 7, respectively, followed by conclusion in Section 8.

Watershed Modeling with SWAT in the Available Literature
The SWAT model was developed in the early 1990s by Jeff Arnold [14], and it has been recognized worldwide as an effective tool in water resource management for assessing the impact of the climate on water supplies and nonpoint sources of pollution in watersheds [15,16]. Moreover, SWAT is a scientific tool used to evaluate streamflow, agricultural chemicals and sediment yield in a large basin [17]. For instance, SWAT was applied to a semi-arid climate in India for rainfall-runoff modelling of river basins [18]. SWAT was also applied to short-term climate data for the assessment of potential hydropower in Assam, India [19]. Correspondingly, the climate change impact on hydropower safety in Dak Nong, Vietnam, was carried out using SWAT [20].
Similarly, hydrological modelling of the Hoa Binh reservoir in Vietnam was conducted to optimize the utilization of flood control and hydropower generation. The results revealed a significant reduction in the peak flood downstream during the rainy season and a stable reservoir level during the dry season [21]. However, the hydrological model in the upper Mekong Basin identified a significant variation from the normal seasonal characteristic of river discharges since the hydropower began operation [22]. Moreover, water balance analysis in SWAT was used to quantify agricultural water demand for the sub-arid Mediterranean watershed [23]. SWAT modelling was carried out in a snowy area of Istanbul, in both Asia and Europe, to evaluate the water budgets of water resources in the context of uncertainties due to climate change and population growth in urbanized areas [24]. Land use change was evaluated using SWAT by simulating the streamflow of Murchison Bay in Uganda to further estimate sediment yield and nutrient loss for water resource management [25]. The groundwater analysis of the Taleghan Dam in Iran was also analysed by simulating the runoff river simulation using SWAT [26]. In addition, SWAT was used to demonstrate the importance of precipitation inputs as the main cause of uncertainties during the simulation of the Adige River basin in Italy, using multiple types of precipitation inputs [27]. Researchers found that monthly simulation produced better results than daily simulation in the ungauged Tonle Sap Basin in Cambodia [28]. Furthermore, SWAT was introduced to simulate runoff of the Mekong River to evaluate the hydrological application of tools in large basins [29].
In some parts of the Philippines, SWAT is utilized in different applications, such as for assessment of potential hydropower in the Visayas [30], Misamis Occidental [5], and the Agusan River basin [31], for predicting runoff in an ungauged watershed in Mabacan [32], and for simulating sediment yield in the Layawan watershed in Mindanao, to investigate land use change [33].
However, most of the applications in developing countries face a lack of precipitation and river discharge data, resulting in reliability issues in the validation process [34]. Therefore, we attempted to address how to overcome the data scarcity in the selected study area. Thus, this work used three types of precipitation datasets. Two are gridded datasets with a resolution of 0.5˚ latitude and 0.5˚ longitude, the NCDC-CPC dataset, and 1˚ latitude and 1˚ longitude, the GPCC dataset, as presented in Table 1. Hence, precipitation datasets were assigned to four stations to proportionally represent large areas of the MRB. Moreover, the calibration was conducted for 3 rivers: the Nituan, Libungan and Pulangi rivers. Then, validation of the calibrated models was conducted through the proxy basin to facilitate data scarcity. Finally, the SWAT interface was implemented in ArcGIS to carry out a watershed model of the MRB despite the limited hydrological datasets available for validation.

The Philippines and Mindanao
The Philippines is situated in Southeast Asia and includes three major zones: Luzon, Visayas and Mindanao, as shown in Figure 1(a). Mindanao is located in the southern Philippines and is the second-largest group of islands, after Luzon [35]. The Philippines had a population of 100,981,437 during the 2015 census [8] and includes a total of 7,107 islands with a land area of 300,000 km 2 (Figure 1(b)) [36]. Moreover, the Philippines comprises 17 regions, 80 provinces, 143 cities, 1,491 municipalities, and 42,028 barangays (villages) [36]. Mindanao has a land area of 120,812 km 2 and had a population of 24,135,775 during the 2015 census, as shown in Figure 1(c) [8], and it is subdivided into six administrative regions that further split into 27 provinces, 35 cities, and 422 municipalities [36].
The Philippines has 421 principal river basins and 18 major basins according to the National Water Regulatory Board (NWRB) [30]. Additionally, the Philippines has four types of climate, as defined by the spatial distribution of monthly rainfall, and experiences an average of 20 typhoons annually [7]. From 1990 to 2006, 520 disasters were induced by seven major natural hazards in this country, affecting 19,298,190 families (approximately 95 million people), who were repeatedly hit by natural hazards such as tropical cyclones, floods and landslides within the same period [37]. Considering these characteristics, water resource management is very challenging in the Philippines because of seasonal weather changes. Modelling is an alternative approach to account for the weather factors that influence the watershed.

Mindanao River Basin (MRB)
The MRB is the second-largest basin in the country [38], with a total area of 21,503 km 2 [39]. It lies in four regions covering 72 municipalities and 1,732 villages in eight provinces, namely, Maguindanao, Lanao del Sur, Bukidnon, Sultan Kudarat, Davao del Sur, Davao del Norte, North Cotabato, and South Cotabato, as shown in Fig. 2 [39]. Due to the dependency on rain throughout the year, the MRB climate was classified as Types 3 and 4 under the modified Corona Climate Classification System of the PAGASA [4]. Moreover, the MRB includes major rivers, such as the Mindanao River and the Tamontaka River, which enters the sea of the lowest part of Cotabato City in Maguindanao [2]. The Pulangi River originates from Bukidnon Province. The Ambal-Simuay River has its waterhead in Lanao del Sur, and the Ala River navigates the Ala Valley in the south [2].
Moreover, the MRB is located at the coordinates of 124º47'35.71' longitude and 7º12'17.06" latitude [4]. MRB has three vast marshes, namely, the Ligawasan, Ebpanan, and Libungan marshes, located within the central and lower parts of the basin. Thus, this large water resource in the MRB will be a potential asset to enhance hydropower development in support of the economic growth of the nearby regions. Therefore, the main reason to choose this study area is to maximize the application of the potential water resources for sustainable hydropower development in Mindanao.

Datasets
This study mainly used the available datasets from certain government agencies in the Philippines, as presented in Table 1. These datasets were requested from the corresponding listed agencies, but the available records are limited. The precipitation datasets from global gridded models were obtained online from the corresponding websites.

Digital Elevation Model
The synthetic aperture radar digital elevation model (SAR-DEM) with a 10-m resolution was obtained from the University of the Philippines Training Center for Applied Geodesy and Photogrammetry (UP-TCAGP) [30,31]. These datasets were collected from point cloud data at a rate of 300 to 400 km 2 per day at every sensor by the use of airborne light detection and ranging (LiDAR) technology and appended with SAR in some areas of concern [40]. This DEM was mostly used in the related previous studies because of its high resolution and accessibility [5,38,40]. The DEM was projected with the Universal Transverse Mercator (UTM) Zone 51 projection and World Geodetic System (WGS) 1984 as the horizontal datum, as displayed in Figure 2(a).

Administrative Boundaries
The Philippines administrative boundaries were obtained from a global administrative map and compared with other shapefiles from the Philippines GIS organization. The administrative shapefile was then projected with UTM 51N and WGS 1984, and then the MRB areas were overlaid onto the provincial boundaries, as shown in Figure 2(b), and municipal boundaries, as shown in Figure 2(c). The MRB lies in 72 municipalities of 6 provinces within 4 regions of Mindanao.

Land Use and land Cover
Land use and land cover data from Landsat 8 of 2010 with a 30-m resolution were obtained from the National Mapping and Resource Information Authority (NAMRIA) and validated on the ground in 2015 by the agency. The shapefile of this dataset also used the same projection as the DEM. This dataset is among the main components of the model structure. Thus, the HRU was determined using this dataset, and the reclassification results are presented in Figure 2(d). The HRU results reported the following figures: the total area of the watershed is 2,041,449.74 ha, including the comprising agriculture area (52.65%), bushland (23.78%), open forest (7.84%), closed forest (5.71%), marshland (4.01%), grassland (4.37%), water (1.12%), and built-up area (0.53%). The large agricultural land of the MRB with an area of 1,074,869.37 ha comprises perennial crops and annual crops.

Weather Records
The weather records for 22 years, as shown in Figure 3, were obtained from DOST-PAGASA. Three datasets, temperature, humidity and wind speed, are available at four DOST-PAGASA stations in Cotabato, General Santos, Davao and Malaybalay, for the period from 1995 to 2017, but solar radiation is available at only the General Santos station for 2016-2017. The weather dataset is a main input for the SWAT model [14,17]. Thus, these weather records were applied to simulate the MRB watershed model. The four stations are shown in Figure 1(d) and were used simultaneously during the model simulations.

Precipitation
The precipitation is a sensitive input in SWAT modelling because of the direct effects on the streamflow output [27,41]. However, the study area has only 2 weather stations, the Cotabato and Malaybalay stations, located within the domain of the MRB. Two other weather stations, General Santos and Davao, are located outside the MRB. These four stations are not enough to represent the large MRB. Therefore, to address this concern, the datasets from the global gridded precipitation model from NCDC and GPCC were compared with the observational datasets from DOST-PAGASA [12]. Thus, the precipitation records from the abovementioned four DOST-PAGASA stations were then used in a comparison with datasets from the closest points of two global gridded precipitation models to evaluate the quality of the rainfall dataset. The multiple precipitation types were used to address the concern on lack of access to quality inputs in a developing country [34].
The NCDC-CPC precipitation is described as a global daily spatial coverage with a resolution of 0.5° latitude and 0.5° longitude covering 1979 to 2017 [42], while the GPCC precipitation [43], a global daily land surface precipitation with a resolution of 1° latitude and 1° longitude, covers the period from 1982 to 2016. The comparison results of the precipitations are summarized in Table 2.
Moreover, the three precipitation datasets were applied for simulating the watershed of MRB to improve the results of simulated discharges. The MRB was simulated from 2000 to 2017 and assigned 3 years of warm-up. Then, datasets from 2005 to 2010 were used to evaluate the precipitation responses against the simulated discharges during calibration and validation of the models based on the available period of river discharge records, as shown in Figure 4.

Method
As described, SWAT was designed for agricultural, non-point source pollution and runoff river flow research. However, it has many features; for example, it can model stream flow by validating the simulated discharge from measured discharges. The hydrological cycle in SWAT is controlled by the water balance equation presented in Eq.
(1). Thus, this water balance equation drives the physics of SWAT, allowing it to model the watershed of a certain basin [24,44]. Where

Analysis Procedure
This study applied the following procedure to model the MRB watershed with SWAT. Each step of the procedure was clearly stipulated to ensure the process during the model simulation. Moreover, the analysis procedure is summarized in Figure 5 to provide a clear overview of the methods being applied with the SWAT model.

Model flow Input Output
Data preparation First, input weather datasets were prepared to match required input formats. The geospatial datasets such as the DEM, land use and land cover were processed and unmasked within the basin boundaries. Then, the watershed is delineated by using the processed inputs of the geospatial datasets. The processed land use and land cover were used to generate the HRUs by ArcSWAT. The processed weather dataset was loaded and utilized for the entire modelling process. Then, initial values of model parameters are set up which can be adjusted after calibrating simulation results in SWAT-CUP interface.
Then, calibrations were carried out to improve model performance by adjusting the parameters based on the sensitivity analysis from the SWAT-CUP (SUFI2) outputs. The watershed of the MRB was delineated into 107 subbasins, as shown in Figure 6. Then, the calibration period was carried out from 2005 to 2010, as reflected in Figure 4. Due to a limited number of river gauges in the study area, the only rivers with a record were the Nituan River, from 2005 to 2010 the Libungan River, from 2006 to 2008, and the downstream Pulangi River, from 2009 to 2010. Therefore, these rivers were utilized for calibration: sub-basin 28 was assigned to the Nituan River (Figure 7), subbasin 40 was assigned to the Libungan River (Figure 8), and sub-basin 45 was assigned to the downstream of the Pulangi River (Figure 9). Then, validations were carried out using the proxy basin principle [45] due to the relatively short records of river discharges. For instance, the calibrated/fitted parameters of River A (the Nituan River) were applied to River B (the Libungan River and Pulangi River) to validate the simulated discharges against the observed discharges. Since the record data of stream flow is not enough to split into two equal parts therefore the same datasets of the Nituan, Libungan, and Pulangi were also used in the proxy basin validation.

SWAT Calibration
The SWAT-CUP program was established to support the SWAT tool to minimize concerns about uncertainties [46] [47]. Moreover, a sensitivity analysis was included inside the SWAT-CUP to tune the parameters according to the recommended results from a number of iterations [48]. Hence, the range of parameters was also enumerated in the absolute values section of SWAT-CUP [47]. Furthermore, SWAT-CUP also incorporates the statistical formulas for the Nash-Sutcliffe coefficient (NSE), coefficient of determination (R 2 ) and percent bias (PBIAS), as presented in Appendix A, to evaluate the model performance. Thus, this study used the sequential uncertainty fitting version 2 (SUFI2) in the SWAT-CUP program for the calibration of the model performances for the Nituan, Libungan and Pulangi rivers.
The calibration was executed with 500 simulations in every iteration, and we conducted 5 iterations per recommendation in the previous studies [46,47,49]. Model calibration does not guarantee the improvement in the model performance. However, it helps modellers evaluate uncertainties by elucidating the sensitivity of parameters that can be adjusted with the observational datasets to improve the statistical results and model fitness [46]. Here, the calibration process was executed cautiously, and the 11 sensitivity parameters shown in Table 3 were investigated, but it was observed that the model performances always depend on the quality of the weather inputs, especially the precipitation, and model structure.
In addition, there are two ways to adjust the parameters: by manual calibration in the ArcSWAT itself or by the use of the SWAT-CUP tool, as previously described. Regardless of these options, the purpose of the calibration process is to improve the model fitness and statistical indicators to ensure the quality of the model outputs. Thus, this study used the general evaluation criteria recommended for the watershed to evaluate the model performances [49,50], as presented in Table 4.

Calibration
The calibration results were measured through the statistical indices shown in Table 5  A negative NSE means that the model performance is unsatisfactory and is characterized by extreme values [50]. A negative statistical performance indicates that the observed average streamflow is better than the simulated streamflow. The simulated discharges at the Libungan and Pulangi rivers are underestimated due to a lack of information on the reservoir management of dams and the water withdrawal from agricultural irrigation in the study area. The system of water storage, release and distribution has a very significant effect on the behaviour of river discharge in a watershed [51]. Then, this has an excessive impact on model performance, as shown in Figures 8 and 9 In addition, sub-basins 40 and 45 have a wetland (marshland) component, and the elevation difference between upstream and downstream is high. Therefore, the terrain of the MRB has a large heterogeneous component that may not be able to be accounted for during the modelling process. The SWAT application is very challenging in a largescale model and in wetlands. The wetlands normally absorb the surface and subsurface water between the inlet and outlet at several points, even if the inlet is well-defined [52]. Although SWAT already employs the basic concept of wetlands (marshlands), its ability to emulate the riparian wetland-river interaction is still under-studied [53]. Furthermore, most of the statistical model results do not show a significant difference between them. Results of this study are similar to the previous studies conducted in the Philippines due to lack of accessibility of datasets [34]. For instance, the SWAT modelling results in Pagsanjan-Lumban Basin in the Philippines show an R 2 and NSE of 0.42 and 0.22, respectively, because of the lack of information on dams and paddy areas (Marshland) [51]. The SWAT water balance simulates the seasonal average discharges of 6.53 m 3 /s for DOST-PAGASA, 5.26 m 3 /s for NCDC-CPC, and 4.73 m 3 /s for GPCC in Nituan River. The seasonal average simulated discharge in Libungan River are 5.08 m 3 /s for DOST-PAGASA, 3.97 m 3 /s for NCDC-CPC, and 4.37 m 3 /s for GPCC. The seasonal average simulated discharges in Pulangi River are 261.69 m 3 /s for DOST-PAGASA, 215.63 m 3 /s for NCDC-CPC, and 250.45 m 3 /s for GPCC. The observed seasonal average discharges are 6.20 m 3 /s in Nituan, 11.37 m 3 /s in Libungan, and 470.64 m 3 /s in Pulangi. Among the simulated discharge models, the DOST-PAGASA model has closer values with the observed river discharges. Therefore, the DOST-PAGASA simulated discharges in the river mouth of the Mindanao river basin which is located in sub-basin 42 has an average simulated discharge of 502.03 m3/s, with the peak and lowest discharges of 1,239 m 3 /s and 2.29 m 3 /s, respectively.

Validation: Proxy Basin
Access to enough data is a common problem in a developing country, and the Philippines is into exception [34]. Due to the challenges of limited access to river discharge records and a limited number of river gauges located in the study area, model validation was used to create a proxy basin. This principle is not new; it was introduced under the hierarchical scheme for systematic testing of hydrological simulation [45]. It was explained that streamflow at Basin C is to be selected as ungagged, and two basins within the region will be selected as gauged rivers, for example, Basins A and B. Then, the models will be calibrated in one of these basins and validated in another basin within the region. For instance, the model is calibrated in Basin A and validated in Basin B, or vice versa. In addition, the proxy basin is useful if the available streamflow record in a basin is not insufficient for equal split-sample and only if two validation results are acceptable and identical [45]. Moreover, the proxy basin was also used for model development of ungagged basins and at regional scales of watershed models [54]. The validation of the calibrated model was also carried out with SUFI2 of SWAT-CUP at the same rivers and in the same period as mentioned in the calibration section. However, the conventional way to split the available dataset is not applicable due to the insufficient length of the data record for river discharge and inconsistency of the duration, as shown in Figure 10. Therefore, from the proxy basin, the resulting parameters from the Nituan River were applied to the Libungan and Pulangi rivers, and vice versa, in this study to conduct the validation of the calibrated model and address the issues on data scarcity.  Table 4. In contrast, the NSE values for all the sub-basins remain unsatisfactory with negative values, as depicted in Table 6. The simulated discharges at Libungan and Pulangi remain underestimated. In addition, the p-factor and r-factor of all the models did not change significantly from calibration to validation, even though proxy basins were applied. With these results, the calibration and validation satisfied only the Nituan River model. The Libungan and Pulangi rivers are the subject of more in-depth studies. Thus, the model performance will improve only if the dam management and irrigation water withdrawal will be accounted for in the model. Although datasets for dams and irrigation are not available; however, these results may be useful for understanding the role of dams and irrigation systems, especially downstream of rivers. In addition, the purpose of validation is to evaluate the fitness of the model after tuning the parameters during the calibration process. The proxy basin validation was carried out to evaluate the application of the SWAT model in a large basin with limited hydrological datasets. Therefore, regional modelling of large-scale basins with limited datasets does not accurately account for all the heterogeneous characteristics of the watershed.

Effects of Precipitation on Simulated Discharge
Since river discharge is basically characterized by precipitation patterns, the observed precipitation inputs were carefully examined by evaluating the global gridded datasets and ground-truth dataset discussed in the previous precipitation section. The results of the correlation between the global gridded datasets and the observational dataset are 0.46 for the NCDC-CPC data and 0.78 for the GPCC data at General Santos, 0.90 for the NCDC-CPC data and 0.83 for the GPCC at Cotabato, 0.95 for the NCDC-CPC data and 0.63 for the GPCC data at Davao, and 0.92 for the NCDC-CPC for the GPCC and 0.52 for the GPCC for the GPCC at the Malaybalay station.
These strong correlations of the precipitation at the 3 stations indicate that the precipitation datasets have similar intensity and seasonal precipitation patterns in Mindanao. Thus, both global gridded precipitation models, either NCDC-CPC or GPCC, and DOST-PAGASA produce a trend of simulated discharges. However, among them, the precipitation model from DOST-PAGASA produced higher simulated discharges, as observed in the peak and base flow in Figure 10. In addition, the characteristics of simulated discharges physically influence the precipitation behaviours, as presented in Figures 7 to 9, respectively.
The three precipitation models proved that a high intensity of rainfall usually occurs in the months of June, October, and November, as reflected in the peak flow of the discharges in Figure 10. Moreover, the simulated peak discharge is slightly higher than the observed peak discharge in June for the Nituan and Libungan rivers, unlike that at the Pulangi River, which has a lower simulated peak discharge in June and November. Therefore, this situation might be affected by the dam system upstream of the Pulangi River. Thus, precipitation patterns for the Nituan River are valuable information for the nearest basin and for the upstream sub-basin to create a scenario for monthly seasonal discharge. However, for the Libungan and Pulangi rivers, the precipitation pattern is not beneficial for hydropower analysis because the simulated model is underestimated, and the area of interest is in wetlands with low elevations. Then, this terrain is not suitable for hydropower development; the ground is soft, and the elevation difference is not enough to generate power. Therefore, for the purpose of this study, utilization of the Nituan River will be suggested for the planning and development of hydropower in the study area.

Quality of Observation
As stated in the previous sections, the number of river gauges and river discharge records are very limited in MRB, and there are only 3 gauges with inconsistent records. Hence, the methods and types of instruments used for data collection were not mentioned in the source. This also affects the model performances; for instance, s shorter period of recorded data is most likely to produce low statistical index results, as observed for the Libungan and Pulangi rivers. Moreover, the ideal calibration and validation require enough river gauges data to split a dataset into two equal periods.
The quality of the river discharge has substantial effects on model calibration and validation. Aside from the length of records and the number of gauges available in the study domain, the method of collecting and processing the data are very serious factors to be considered. As explained previously, these datasets are part of a pilot project to support the technical capability of responsible agencies in certain regions. In short, a transitional process for transferring technology know-how from the service provider to the agency might lead to transitional development. Thus, the datasets duration is inconsistent among the stations, and the number of gauges is very limited despite the large size of this basin, with hundreds of streams and rivers. Therefore, as an alternative, using the regional watershed model and proxy basin process is an efficient way to implement model validation in the MRB.

Conclusion
The main purpose of this study was to contribute to the improvement in water resource management for hydropower applications. Watershed modelling in the MRB in Philippines was carried out using SWAT. Since precipitation is a critical input that has a direct influence on river discharges, measured precipitation dataset within the the MRB is desirable as many as possible. However, only 2 DOST-PAGASA weather stations are available in MRB while another 2 stations are located close to MRB, despite the large area. Therefore, global precipitation gridded datasets from NCDC and GPCC were utilized to investigate the quality of precipitation in the MRB. Each precipitation dataset was individually used in SWAT simulations. Then, the simulated discharge was calibrated at 3 river gauges in the Nituan, Pulangi, and Libungan rivers using SWAT-CUP (SUF1). Moreover, due to limited short records from the 3 river gauges, the proxy basin process was applied for the validation of calibrated models of the same rivers. The models of the Nituan River provide better results compared with those of the Libungan and Pulangi rivers, even though a calibration was executed, and a proxy basin was also applied.
Lack of enough and qualified data is a common problem in a developing country. Due to the challenges of limited access to river discharge records and a limited number of river gauges located in the study area, model validation based on a proxy basin principle was applied in this study; for instance, the model was calibrated in Basin A and validated in nearby Basin B, or vice versa, to overcome the data scarcity in the study area. The proxy basin was also used for model development of ungagged basins at regional scales of watershed models. Therefore, this approach of calibration and validation of watershed modelling based on the proxy basin principal with various precipitation inputs demonstrates a method of watershed modelling for regions with insufficient precipitation and discharge data in developing countries.
The underestimated results of the Libungan and Pulangi rivers are mainly due to a lack of information on dams and irrigation water withdrawal. The study area has a vast wetland (marshland) and is characterized by a high elevation difference between the upstream and downstream, contributing to the uncertainty in the models and weak performance of the models in wetlands with limited datasets. Thus, these findings will be applicable to only the upstream side of the MRB for identifying potential sites for hydropower development.
In addition, comparison results for the precipitation inputs may be useful references to improve the meteorological measurements by adding additional weather stations in the regions. The trend of simulated discharges against precipitation inputs at 3 stations demonstrates the monthly seasonal characteristics of a watershed. Thus, this information can be used to determine which month might have a high potential for hydropower generation. Finally, this study specifically discussed the importance of meteorological agencies improving data collection and the application of the collected data in large areas; additionally, this study discussed the very significant role of river discharge, dam management and agriculture water withdrawal data for watershed analysis.

Appendix A: Equations of the Statistical Indices Used to Evaluate the Model Performance
The index of agreement between two variables is computed by the following formula [55,56]: Where R 2 is the coefficient of determinants is the measured discharge Ǭ m is the mean of the measured discharge is the simulated discharge Ǭ is the average of the simulated discharge The root mean square error (RMSE) is used to measure the absolute fitness between the observed and the modelled results: Where RMSE is the root mean square error of the samples P i is the predicted values in a sample is the observed values in a sample n is the number of observations The Nash-Sutcliffe coefficient is calculated as follows [57,58]: Where NSE is the Nash-Sutcliffe coefficient is the measured discharge Ǭ is the mean of the measured discharge is the simulated discharge PBIAS measures the model fitness in terms of average tendency, and a small value of PBIAS indicates a better model fitness. PBIAS is calculated by [47,57]: Where is the measured discharge is the simulated discharge n is the number of observations RSR is a standardization of the RMSE to measure how well the model results fit with the observed values, and a lower value of RSR indicates a better model fitness. RSR is calculated by the following equation [47,57]: Where is the measured discharge Ǭ is the mean of measured discharge is the simulated discharge n is the number of observations The Kling-Gupta efficiency (KGE) is used to examine the decomposition of the Nash-Sutcliffe efficiency, with a value close to 1 indicating a better performance. KGE is calculated by the following equation [47,59]: Where r is a linear regression coefficient between the simulated and the observed data α is a ratio of standard deviation between the simulated and measured data (α= ) β is a ratio of the means between simulated and measured data ( = ).
The r-factor is the thickness of the 95% predicted uncertainties (95PPU) [60]: Where ,97.5% is the upper boundary of the 95PPU (95% predicted uncertainties) ,2.5% is the lower boundary of the 95PPU (95% predicted uncertainties) is the observed standard deviation.