Probable streamflow changes and its associated risk to the water resources of Abuan watershed, Philippines caused by climate change and land use changes

One of the main impact areas of climate change (CC), and land use and land cover change (LULCC) is the hydrology of watersheds, which have negative implications to the water resources. Their impact can be indicated by changes on streamflow, which is quantifiable using process-based streamflow modelling of baseline and future scenarios. Here we include the uncertainty and associated risk of the streamflow changes for a robust impact assessment to agriculture. We created a baseline model and models of CC and LULCC “impact scenarios” that use: (1) the new climate projections until 2070 and (2) land cover scenarios worsened by forest loss, in a critical watershed in the Philippines. Simulations of peak flows by 26% and low flows by 63% from the baseline model improved after calibrating runoff, soil evaporation, and groundwater parameters. Using the calibrated model, impacts of both CC and LULCC in 2070 were indicated by water deficit (− 18.65%) from May to August and water surplus (12.79%) from November to December. Both CC and LULCC contributed almost equally to the deficit, but the surplus was more LULCC-driven. Risk from CC may affect 9.10% of the croplands equivalent to 0.31 million dollars, while both CC and LULCC doubled the croplands at risk (19.13%, 0.60 million dollars) in one cropping season. The findings warn for the inevitable cropping schedule adjustments in the coming decades, which both apply to irrigated and rainfed crops, and may have implications to crop yields. This study calls for better watershed management to mitigate the risk to crop production and even potential flood risks.


Introduction
Impacts of climate change (CC) to the environment have been evident and are forecasted to continue this century and beyond, as mankind is experiencing climate and weather abnormalities indicated by droughts, irregular seasons and extreme typhoons among others (Watson et al. 1996). One main impact area of CC is the hydrology of watersheds, where changes in precipitation and temperature directly affect the dynamics and supply of water resources, and eventually the water stakeholders who will struggle to meet their water demands (Arnell 1999).
A direct causal factor of climate change is land use and land cover change or LULCC (Dale 1997). Carbon emissions from LULCC were found 33% of the total anthropogenic emissions from 1850-1990(Houghton 1999) and 12.5% from 1990-2010 (Houghton et al. 2012). In a watershed context, LULCC are often driven by socioeconomic pressure and lapses in land use management (Overmars and Verburg 2005). Forest conversion is one of the worst forms of LULCC. Forest loss in the form of deforestation reduce the ecosystem services provided by watersheds, primarily water supply and regulation (Rawlins et al. 2017). LULCC and other factors that influence the hydrologic system of watersheds can be modelled.
Hydrologic models based on the principles of physical processes is called process-based models, first introduced in the early 1960s (Fatichi et al. 2016). This model makes use of model inputs with unique hydrologic parameters to simulate rainfall-runoff models. Using land use, soil, slope and climate data as main inputs, water resources are quantified in the form of streamflow (Q) or river flows, usually in m 3 /s. Process-based modelling has been the state-of-the-art to assess the impacts of CC and LULCC in watersheds by modifying climate (or weather) and land use (or land cover) model inputs.
The impacts of both CC and LULCC on streamflow changes can be assessed together. For instance, streamflow decline over a 50-year period was attributed mostly to climate parameters (70%) and the rest to LULCC (Beguería et al. 2003), while a recent study by Hung et al. (2020) depicted the same trend on peak discharges changes, where 83% is attributed to the CC and 17% to the LULCC. Though streamflow changes are more climatedriven, the effects of LULCC is believed to effect the seasonal streamflow (Kim et al. 2013). Scenario modelling exercises revealed that deforestation is a causal factor of drought and vice versa (Staal et al. 2020), while forests converted to bare lands have resulted into abrupt runoff during typhoon events, and hence higher flood risks (Rawlins et al. 2017, Kim et al. 2019. As climate extremes seem to be inevitable, and with the continuity of deforestation-driven LULCC, modelling and understanding their impacts to water resources of watersheds is becoming imperative. Watersheds are management units that have natural hydrologic boundaries, from upland ridges to lowland service areas. Managing watersheds should be a joint-responsibility of the stakeholders, not only by the government, but also by the public (Cruz 1997). That notion gave birth to the integrated watershed management (IWM), a concept that accounts for both ecological and socio-economic aspects of managing watersheds and its resources collaboratively (Biswas 2004). In the Philippines, IWM is implemented at major river basin levels. Each major basin has a master plan enacted by the stakeholders, which embeds short-and long-term management strategies on water resources (Almaden, 2015). That includes accounting the water resources of the basins and assuring policy-support on appropriate water uses. Special attention is given to the agriculture sector being the major stakeholder of water resources (Tabios et al. 2018); and to the energy sector due to the hydropower potential of watersheds. An example IWM is Abuan Integrated Watershed Management Project, which provided CC adaptation strategies to agriculture and livelihoods, while increasing disaster resiliency (USAID 2016). The project also acknowledged the risk to forest loss and the importance of forests in water and soil regulation. Watershed management is even integrated into land use plans of local governments for bottom-up management approach (Briones et al. 2016). This formulates an institutional framework to co-manage protected areas and watersheds, and make sure that CC adaptation and mitigation strategies exist locally (Lasco et al. 2008). IWM and efforts alike have strong demand on modelling CC-LULCC impacts (e.g. water deficit) to serve as decision-support for planning.
While the effects of CC and LULCC on water resources of Philippine watersheds can be modelled (see Principe 2012;Rawlins et al. 2017;Balderama et al. 2019), modelers and planners should assess the uncertainties not only from streamflow models and its parameters, but also from the climate projections used (Su et al. 2017). The Philippine government published newer probabilities of rainfall and temperature changes, which were not used during the earlier river basin master planning between 2011 and 2016. As of now, there is very limited information on how (and how probable) water resources would change using these new projections. For streamflow parameters with probability distribution, an uncertainty propagation is necessitated, in parallel with streamflow model calibration (Abbaspour 2013). The state-of-the-art for propagating uncertainties from streamflow models is through the twostep use of Soil and Water Assessment Tool and its Calibration and Uncertainty Program (SWAT and SWAT-CUP). SWAT is a process-based model that discretizes watersheds into smaller spatial units with unique hydrologic parameters for better water balance simulations, while SWAT-CUP is an uncertainty propagation program designed specifically for SWAT model calibration (Arnold et al. 2012). In the end, CC and LULCC streamflow models are calibrated to provide optimal streamflow simulations with confidence intervals.
Taking into account the probability of streamflow changes can provide better decision-support to the water resources stakeholders (Sivakumar 2011, Ouyang et al. 2015. Better decisions in this context can be outcomes of a robust impact assessment towards lowering climate change risks through adaptation strategies e.g. for the agriculture sector (Abbasi et al. 2020). Risk assessment has been identified as the missing step after calibrating streamflow models with uncertainties (Abbaspour et al. 2018). In this context, we aim to implement not only the best practices in streamflow modelling, but also risk assessment based on probable streamflow changes caused by CC, and both CC and LULCC. Our specific objectives are: (1) calibrate a baseline streamflow model and create models of CC and CC-LULCC (impact scenarios) using SWAT and SWAT-CUP; and (2) assess the probability of water deficit or surplus and calculate its associated risk to the agriculture sector. We selected Abuan watershed as the study area for being a critical watershed-with vast water resources for agriculture, hydropower, and domestic use; but subject to conservation and protection.
The succeeding sections starts with an overview of the study area (Sect. 1), followed by the methodology section (Sect. 1) with sub-sections about the workflow diagram; SWAT and SWAT-CUP modelling; CC and LULCC impact scenarios; and risk assessment. The results (Sect. 2) start with the outcomes from calibration and validation; sensitivity analysis of parameters; impacts indicated by waters surplus and deficit; and the risk assessment to agriculture. The discussion of the key findings (Sect. 3) revolves around the magnitude of streamflow changes and its impacts and the importance of uncertainty and risk assessment, then leading to several concluding remarks (Sect. 4).
Abuan has a size of 63,791 ha, where the largest part is inside Sierra Madre Natural Park. Its mother basin, Pinacanauan de Ilagan, is an elongated basin parallel to the Sierra Madre mountains ranges, with vast remaining forest cover and biodiversity. From the ridges of the mountains, the western part of the watershed extends to the lowlands with 373 ha built-up areas and 2,913 ha croplands (mostly corn) farmed by 2900 households (Balderama et al. 2016); while the east-bound area is oriented to the eastern seaboard of the Philippines. Since Abuan is mountainous and nearby the sea, the watershed and the rest of Isabela province is receiving an average of 2430 mm of rain annually, more or less distributed year-round (Type 4 climate), and wettest around December because of the north-east monsoon. During monsoon season and typhoon events, Abuan's lowland are high-risk to flooding; while during dry season, its croplands can also be prone to droughts (Balderama et al. 2016).

Workflow and main inputs
The workflow of this study (Fig. 2) has five main steps starting from streamflow modelling, initial assessment, calibration-validation, impact scenario building, and risk assessment. The last step is an additional step, identified as the missing piece towards comprehensive streamflow modelling (Abbaspour et al. 2018). In terms of the main inputs, weather data from local stations dated 1991-2010 and the climate projections were obtained from Philippine Atmospheric Geophysical and Astronomical Services Administration (PAGASA), the official weather office of the country. The daily river discharge data (average of 3 daily readings) of Pinacanauan de Ilagan epoch 2003-2010 was obtained from Department of Public Works and Highways, the office with extensive streamflow database from nationwide gauge stations for water infrastructure designs. Soil data was downloaded from FAO digital soil map at 9 km resolution (Batjes 2012) and paired with the pre-formatted soil database from Soil and Water Assessment Tool at Map Windows (MSWAT). Lastly, the land cover data dated 2000 was produced using a supervised classification of Landsat satellite images at 30 m resolution, having an overall accuracy of 95%. The succeeding sections sequentially follow the workflow in

Streamflow modelling using SWAT
To model streamflow, we used Soil and Water Assessment Tool (SWAT), a semi-distributed model which perform rainfall-runoff simulations among hydrological response units (HRU) or unique combinations of soil, slope, and land cover (Arnold et al. 2012). Hydrologic properties of HRUs are then aggregated into watersheds or sub-basins of Stochastic Environmental Research and Risk Assessment (2021) 35:389-404 391 the main basin and assessed through a water balance equation (Eq. 1).
where t is time in days, SWt is the final soil water content, SW is the preliminary SWC, R is rainfall, SR is surface runoff, ET is evapotranspiration, P is percolation and RF is return flow. The unit of all parameters is mm.  The number of sub-basins is mostly a user-choice during the watershed delineation step of SWAT. We delineated the sub-basins subjectively since Abuan watershed belongs to Pinacanauan de Ilagan basin, and only that larger basin has a gauge station (river discharge data or observed streamflow for model calibration). A total of 109 subbasins were delineated, discretized to 817 HRUs. Figure 3 shows the sub-basin boundary, three soil types, six land cover classes, five slope categories, and the boxplot of observed streamflow. The map also highlights the geo-location of the gauge station and the outlet of Abuan at subbasin 35. For the main inputs, majority is classified as mountain soils while the rest are clay and sandy loam. Overlapping with the mountain soils are forests as the dominant land cover of Pinacanauan (81% of total) and Abuan (92% of total). From mountains to lowlands, the whole Pinacanauan and Abuan resemble the watershed continuum form (upland-to-lowland) as shown by the slope map. The boxplot of observed streamflow depicts skewness or denser low flows.
For the weather data inputs, daily data of precipitation, wind speed, relative humidity, solar radiation, and minimum and maximum temperature were extracted from the Cagayan Valley Integrated Agricultural Research Center local weather station, located in the lowlands of the watershed. For dates with missing data, weather data from Isabela State University station was used (around 150 km apart from Abuan). In days where the two local stations have missing data, the average of global weather data from four closest stations in the watershed were used originally from NCEP or National Centers for Environmental Prediction (NCEP 2014). The SWAT model requires a ''warm-up period'' to initialize the time-series simulations. We used a 12-year warm-up period (1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002) not only to have a better model tuning, but also to match with the coverage of the observed streamflow. We run the SWAT model to produce mean monthly streamflow from 2003 to 2010. Monthly simulations are commonly used for CC impact assessments to align with the monthly projections of climate models (Teutschbein and Seibert 2010). In parallel, mean monthly observed streamflow was derived from the original daily data.

Uncertainty assessment and model calibration using SWAT-CUP
Uncertainty assessment is a must in streamflow modelling and must be explicit on the assumptions used for the uncertainty estimates (Kiang et al. 2018). Both observed and simulated streamflow have associated uncertainties. However, uncertainty assessment has commonly been implemented only for the latter because error distributions from observed data are often unavailable (also for this study). As such, we treated the observed streamflow as the ground-truth. On the other hand, the uncertainties from simulated streamflow of parametric models like SWAT are commonly assessed using stochastic methods (Abbaspour 2013).
To assess the uncertainty of our streamflow model, we used the Sequential Uncertainty Fitting (SUFI-2) of SWAT-CUP (Abbaspour 2013). This uncertainty propagation method simulates streamflow using a range of parameters (instead of a single value) simultaneously. The range gets narrower as the simulations are repeated every after iteration until the ideal fit to the observed data is attained. The simulations form an uncertainty range (95PPU), expressed as the 2.5% and 97.5% confidence interval of the cumulative distribution of the simulations. The quality of the 95PPU is measured by P-factor or the fraction of observed data inside, and the R-factor or the thickness of the range. The 95PPU is optimized through a user-defined objective function-we used Nash-Sutcliffe efficiency (NSE) by Nash and Sutcliffe (1970). NSE in Eq. 2 is commonly used as a comparative metric between the observed Q o and simulated streamflow Q s where the sum of the squared difference of the two from the simulations (i) is divided by the sum of squared difference of a simulation (i) Q o;i and the mean streamflow of all observations À Q o : Using SUFI-2, we performed the calibration from 2003 to 2007, with 5 iterations of 1000 simulations apiece. The parameter range of the last iteration was used for validation from 2008 to 2010, also with the same number of simulations. The effect of calibration is measured by NSE and coefficient of determination (R 2 ) as shown in Eq. 3where the numerator corresponds to the regression sum of squares and the denominator is the total sum of squares: Moreover, bias of the simulations was assessed using percent bias (PBIAS, Eq. 4) or the % tendency of the simulations to deviate with the observed data: There can also be a temporal bias related to seasonal flows e.g. underestimating the peaks and overestimating the baseflow or low flows (Qiu et al. 2012). We then compared the peak flows and low flows (10% and 15% of total flows, respectively) between the uncalibrated and calibrated model.
Model calibration is dependent on the choice of parameters to calibrate (Paul et al. 2017). To determine the correct choice, we first investigated the results of the uncalibrated model compared with the observed data e.g. based on low flows, peak flows, and discharge shifts (Abbaspour et al. 2015). We then selected 11 potentially sensitive parameters, which included spatial parameters like curve number (CN) and HRU-level coefficients; and nonspatial parameters like ground water, river and soil parameters (see Table 1). These parameters can be calibrated in three ways as implemented in SWAT-CUP: by replacing, adding/subtracting, and weighing. Commonly, we selected the weighing option for spatial parameters and either of the two remaining options for non-spatial parameters. The adjusted parameters of the last iteration were used for validation using the same number of simulations.
Aside from validating the calibrated parameters, we also assessed their sensitivity to the streamflow simulations. Using SWAT-CUP, a global sensitivity analysis was performed, which quantifies the relative significance of a parameter to the calibration process. Using multiple regression, global sensitivity works by comparing the average change of an objective function (e.g. NSE) when a parameter is changed while the rest are constant. Sensitivity is evaluated by T-stat equivalent to a precision measure and P-value to indicate statistical significance. The higher the T-stat and the lower the P-value, the better.

The two impact scenarios
To implement the impact assessment, SWAT models of CC and CC-LULCC were created separately using the calibrated parameters, but with modified weather data and land cover inputs. For the former, we used the new precipitation and temperature projections of the official weather department of the government (PAGASA 2018). They used Climate Information Risk Analysis Matrix (CLIRAM) downscaled model, with baseline decades between 1970 and 2000, to project emission scenarios indicated by RCP or Representative Concentration Pathways: RCP 4.5 (low) RCP 8.5 (high) both for 2050 and 2070. The CLIRAM data and projected precipitation in every emission scenario are shown in Fig. 4.
The land covers 2050 and 2070 were generated using the forest loss trend from the Global Forest Change data (Hansen et al. 2013). We first derived the net forest cover in 2018 and calculated yearly forest loss from 2000 to 2018 at 7%. Assuming continuity of trend, we estimated 12% and 26% forest loss in 2050 and 2070, respectively. The projection is reasonable given that the common forest conversion in the study are patches of slash-and-burn cultivations (Balderama et al. 2019). Using the two projected forest loss, we created pseudo deforested pixels by randomly populating points, and buffering them into 1 ha pixel each point. These pixels and the forest loss from 2000 to 2018 were labelled into open barren. The land cover maps 2000, 2050, and 2070 are shown in Fig. 5.
The impact scenarios were incorporated in the SWAT-CUP models (updated weather and land use inputs). Then, we proceeded with the streamflow simulations (1000 times) in SWAT-CUP using the calibrated parameters, this time including climate parameters namely CO2 (600-800 parts per million), and the precipitation and temperature change per emission and future year. In total, there were eight scenarios namely c50a, c50b, c70a, c70b, l50a, l50b, l70a, and l70b: where c = CC and l = LULCC; 50 = 2050 and 70 = 2070; and a = RCP 4.5 and b = RCP 8.5. We obtained the monthly simulations with 95PPU particularly for Abuan watershed (sub-basin 35).
After obtaining the CC and CC-LULC simulations for Abuan watershed, we derived mean monthly streamflows of 2050 and 2070 scenarios and measured the changes (from the baseline model) using percent change (% change). In Eq. 5, % change is computed by dividing the difference of S2 or the scenario streamflow and S1 or baseline streamflow, to S1, multiplied by 100. Then, two impact indicators of streamflow changes were derived: water surplus if there is water gain or S1 is lower than S2, and water deficit if otherwise or there is water loss.

Risk assessment to agriculture
For each impact and emission scenario, the next step was to assess the risks associated from streamflow changes. In context, risk is the potential monetary loss caused by the impact scenarios relative to the uncertainties of streamflow changes. In calculating loss (or gain) from CC and LULCC impacts, the probability of streamflow changes (95PPU) was used to determine how probable water surplus and deficit are. For this study, we calculated the risk to crop production (Risk) during the planting season from May to August. First, we adopted the crop irrigation requirement and production cost per hectare in the Philippine context (Rawlins et al. 2017;BAS 2009). We then estimated the potential irrigation water as 81.3% of the total water supply (AQUASTAT 2009), and hence the potential irrigated lands. Because of CC and LULCC, potential irrigated lands may decrease, resulting to unproductive lands. As such, the production cost for these unproductive lands cost(D) was calculated and multiplied with the probability of water deficit Pr(D), see Eq. 6. The monthly amounts were averaged and translated to million US dollars.
4 Results Figure 6 shows a hydrograph of the two streamflows (before and after calibration) relative to the observed data and precipitation, while Table 2 shows the accuracy metrics, and the peak and low flows also pre-and post-calibration. In a glance, the hydrograph reveals that the simulated streamflows were rainfall-responsive. Accuracy metrics were consistently higher after calibration and validation, except for a slightly lower R 2 (validation). Improvement was higher for NSE since it is used as the objective function in SUFI-2. The accuracy increase was reflected in the seasonal flows, where underestimated peak flows improved by 164 and 235 m 3 /s (26% and 30%), and low flows increased by 19 and 3 m 3 /s (63% and 8%) after calibration and validation, respectively.

Calibration-validation results
Despite the improvement, underestimation of the peak flows still existed.

Sensitivity of key streamflow model parameters
After calibration, the sensitivity of each streamflow parameter was measured by P-values from the global sensitivity analysis (Fig. 7). Parameters directly influencing the peak and low flows were the most sensitive (lower P-values), including the coefficients on runoff (CN2 and SURLAG), baseflow (ALPHA BF), soil evaporation (ESCO), and soil depth for baseflow (GWQN). On the other hand, insensitive parameters included hydraulic conductivity (CH K2), available water in soil (SOL AWC), and water uptake by plants (GW REVAP and EPCO) among others.

Changes in monthly streamflow as impact indicators
Monthly streamflows with uncertainty range of the two impact and eight emission scenarios were generated using all calibrated parameters (streamflow and climate). Table 3 summarizes the monthly streamflow of the baseline and impact scenarios. Figure 8 shows that result but disaggregated per emission scenario (x-axis of graph) of the two impact scenarios (purple and gold), with reference to the baseline model (dashed red line). Streamflow of the impact scenarios both deviated from the baseline streamflow in all months, on average. Water deficit started from May until October (-12.96%), while water surplus happened in November, December, and April (11.56%). LULCC contributed by -4.23% to the deficit and 8.63% to the surplus. Deficit was highest in May at -14.29%, while surplus was at peak in December 17.57%. A more striking result was the average water surplus in November and December from 1.95% (CC) to 12.79% (CC-LULCC).
Among the eight emission scenarios, water deficit was more evident in CC-LULCC scenarios, particularly in L70b    (RCP 8.5, 2070) in May and L50b (RCP 8.5, 2050) from June to October. Scenario L70b had the most water surplus in November and December. The uncertainty of emission scenarios was a function of streamflow and more evident in months with higher flows e.g. November and December. Moreover, the monthly uncertainty of streamflow change (7-8th column of Table 3) was consistently lower (higher probability) for CC than the CC-LULCC impact scenario. The uncertainty ranges had no clear seasonal trend, except for January which was the least uncertain compared to a similar month e.g. May.

Associated risk to agriculture driven by streamflow changes
Risk to crop production along with other details of its computation is shown in Table 4. The potential water supply of Abuan in one cropping season (May-August) was found to be sufficient to irrigate 16,824 ha of crop lands. Both impact scenarios caused decrease on streamflow, and thus reduction of potential production areas. Water deficit was higher at -9.46% and -18.65% as affected by CC and CC-LULCC, respectively. Because of CC, there was 1,530 ha or 9.10% croplands at risk, and the additional impact of LULCC doubled the risk at 19.13% (3,219 ha). The respective monetary value of the risks was 0.31 and 0.60 MUSD.

The effects of model calibration and the sensitive parameters
Streamflow model calibration resulted into higher accuracy and improved the simulation of peak and low streamflows than the uncalibrated model. The model improvement indicates that our calibration plan on which and how to calibrate parameters of interest is effective (Abbaspour et al. 2015). Specifically, parameters related to runoff, groundwater and soil were pre-identified, calibrated, and found sensitive to the model e.g. CN2, SURLAG, ESCO, GWQMN and ALPHA BF (see Table 1). According to the SWAT technical documentation (Neitsch et al. 2011), a higher curve number (CN2) leads to a higher surface runoff relative to the HRU; an increase in runoff lag time (SURLAG) also increases the fraction of runoff leading to the main stream in a day; and compensating for soil evaporation (ESCO) retains soil moisture. Moreover, decreasing the shallow aquifer depth (GWQMN), and increasing the shallow groundwater recharge (ALPHA BF) logically increases the low flows. The characteristics of Abuan watershed's mother basin complement these adjustments. The basin receives rain all throughout the year, without large groundwater extractions, and dominated by old-growth forests, hence should have a more regulated baseflows year-round (Rawlins et al. 2017). Slower surface runoff though may also affect the simulation of peak flows. These improvements after calibration were proven in similar watersheds in the tropics (Ndomba et al. 2008;Thampi et al. 2010), as examples of many studies that reported model limitations related to peak and low flows. The commonality of the need for calibration in this context may have been influenced by observed data scarcity in the tropics, which limits the initialization or ''warm-up years'' by SWAT models (Nyeko, 2015), and the configuration of SWAT itself, which is by default not in the tropical settings. Tropical watersheds may need to further parameterize plant and soil inputs e.g. using remote sensingassisted plant parameters (Strauch and Volk 2013) and data from local soil samples (Rawlins et al. 2017).
Despite the improvement, the calibrated peak flows were still underestimated. The complexity of quantifying peak flows is attributed to a combination, if not all, of the following factors: watershed size and aggregation of weather data (Lopez and Seibert 2016), the default CN settings of SWAT (Qiu et al. 2012), and incorrect gauge readings during typhoon events (Beschta et al. 2000).

Importance of uncertainty assessment and remaining uncertainties
The remaining errors even after calibration may propagate to further applications of the simulated streamflow. This domino effect can be mitigated if uncertainties are not only quantified, but also integrated to the said applications (Tomkins, 2014;Abbaspour et al. 2018). Uncertainty estimation using SUFI-2 comes in the form of 95% confidence interval (95PPU) after several iterations of uncertainty propagation among parameters. For this step, it is worth reiterating that we propagated future precipitation and temperate within the 10 th and 90 th percentiles instead of using only the median value. After five iterations, we obtained the desired 95PPU which are less uncertain compared to the first iteration. From that iteration, the 95PPU for any sub-basin of the study area can be derived (Abbaspour 2013). The 95PPU of Abuan watershed (sub-basin 35) was confidently used for analysis based on two reasons. First, the basin where it belongs (Pinacanauan de Ilagan basin) with the SWAT and SWAT-CUP models is mostly unregulated. There are no huge dams in the study area that could alter the calibration process if not integrated to the streamflow model (Rahman et al. 2013). Second, Pinacanauan de Ilagan is relatively a small basin (\ 3000 km 2 ) having identical hydrologic parameters with Abuan i.e., similarity on topography, land cover, soil (see Fig. 3) as well as climate and climate projections. Relatively small basins are calibrated for LULCC analysis at sub-basin levels in the Philippine landscape (Briones et al. 2016). For relatively large basins with diverse land-uses, further regionalization or similarity assessment among sub-basins may apply prior to SWAT modelling (Swain and Patra 2017). Watershed regionalization is applicable once streamflow modelling in the Philippines is implemented at major river basin scales (Araza et al. 2020).
The uncertainty range was further integrated to the risk assessment of streamflow changes to agriculture, which is the final step towards a successful and problem-oriented modelling (Abbaspour et al. 2018). That allows more Italic values indicate the average of the four months above it conservative or realistic risk estimates, which is advisable if projections from future scenarios are uncertain i.e. trends in carbon inventory of countries (Grassi et al. 2008). The remaining uncertainties from the streamflow model may come from the observed data, which is not error-free and thus needs an uncertainty estimate (Hamilton and Moore 2012). However, that is uncommon in countries with data scarcity and data sharing issues. In our case, we obtained a preprocessed observed data, but without raw data to derive a probability distribution of streamflow values for uncertainty propagation. Uncertainties may also come from the main inputs (soil, weather, and land cover inputs). We are particularly concerned with the coarse soil input used for this study, dominated by a single class. An alternative would be to use finer resolution soil data e.g. SoilGrids 250 m data, despite the longer preprocessing time to make them SWAT-formatted. We also noticed in Fig. 6 at least three peak flows that were not rainfall-responsive, which suggest that the error either comes from the precipitation data or observed data. Moreover, there can be uncertainties in LULCC inputs not only from the classification itself, but also the way SWAT uses one land cover period while producing time series outputs (Araza 2018). Lastly, there can be uncertainties from the provincial CC projections that can be substituted by spatially explicit projections once available. Nevertheless, we perceive that the uncertainties from the streamflow model itself would be the largest uncertainty source.

In-depth impact assessment
The impacts of streamflow changes were indicated by water deficit from May to October, and surplus in November, December and April. The driving factor of the deficit and surplus are the lowest projected rainfall in JJA quarter and highest in DJF quarter, respectively. These impacts were doubled after adding LULCC in the assessment-worst cases were L70b (RCP 8.5, 2070) and L50b (RCP 8.5, 2050) emission scenarios. Consequently, the projected water deficit will surely affect the crop production especially corn during May-August season. Some months with water deficit included wet months (JJA). This somehow deviates from the projections reported in Balderama et al. (2017) and Tongson et al. (2017) where water surplus was forecasted from the said months in the study area. The discrepancy is mainly because of the newer CC data used for this study e.g. with an average rainfall projection of -17.2% in JJA quarter, RCP 4.5 in 2050 (Fig. 9). The CC projections originate from the CLIRAM downscaled model, a spatially explicit climate model that accounts for every uncertainty source from global to regional scale (Daron et al. 2016). That cascaded uncertainty is common in downscaled models, thus having overly pessimistic or ''risk-adverse'' projections (Jones, 2000;Elsberry, 2002). We counteracted by setting a high number of simulations (1000) during model calibration, but could have been more e.g. 5000 to further narrow uncertainties. Moreover, more localized climate projections can be integrated to the streamflow model once the pixel-level (25 km) CLIRAM projections are publicly available.
The consequence of the timing of water deficit would be the cropping schedule changes for both irrigated and rainfed croplands. The latter is foreseen to be more affected in both wet and dry cropping season. However, slight changes in cropping calendars may affect the yield significantly (Lansigan et al. 2007). A direct solution for possible water shortages is to store the surplus water in November to December, but may need reservoirs for water storage and even flood mitigation. Indirect solutions are sustainable forest conservation and watershed management.
On the other hand, the projected water surplus, particularly in November and December, may increase the flood susceptibility of the study area. Those months showed the highest probability of streamflow change with an increase from 1.95% (CC) to 12.79% more surplus because of deforestation-driven LULCC. This projection should catch the attention of the government, given that most typhoons pass through the study area in the last quarter and there have been warnings of more intensive typhoons to come (PAGASA 2018;Tolentino et al. 2016). Flooding during typhoon events may worsen due to continuous of slashand-burn, increase in sedimentation (Balderama et al. 2019) and displacements of main streams within the study area (Dingle et al. 2019). Towards flood modelling, our SWAT model would need sub-daily data to have more accurate peak flows simulations (Yang et al. 2016

Risk to agriculture and flooding
During May-August cropping season, the impact scenarios put at most 19.13% of potential crops lands at risk: 9.10% was attributed to CC and 10.03% to LULCC. In monetary terms, the risk of CC and CC-LULCC costed around 0.31 and 0.60 MUSD. These estimates, however, are based on the potential irrigation water from Abuan watershed. In reality, most of the corn fields in the study area are rainfed (Lansigan et al. 2007), but gradually being modernized because of the agriculture modernization plan (including irrigation) in the Philippines (Inocencio et al. 2016). For comparison, assuming the actual irrigated croplands as 1/10 of the total, the croplands at risk from CC-LULCC would be lower at 8.98%. Nevertheless, rainfed crops would still be at risk especially in the May-August cropping season because of the projected average rainfall decrease at -17.2% rain during JJA quarter. The estimation of flood risks was not included mainly because of unavailable sub-daily data to accurately model peak flows for flood modelling. For future attempts to estimate flood risk, a two-step method similar to this study can be adopted: hourly streamflow simulation (using SWAT) and their probability (using SWAT-CUP), for both baseline and impact scenarios. Then, using a dynamic flood model i.e. agent-based or process-based, Monte Carlo simulations using the probable hourly streamflow and realizations of a Digital Elevation Model (DEM) will be needed to derive probable inundated areas. Lastly, a damage function based on inundation height for houses (and crops) and their land valuation will be used to calculate the risk to flooding (see Phil-WAVES 2016).

Entry points to management plans
Results of this study can be integrated from local to national plans as basis for aligning programs, projects, and activities (PPA). For instance, mainstreaming CC-LULCC effects into local plans in a participatory manner would not only allow active participation of every stakeholder, but would also integrate impacts to crop yield (Balderama et al. 2016) and localize the risk calculations e.g. farmer-level. Moreover, the results of this study can be contextualized for ecosystem services (ES) accounting. Specifically, the forest water regulation service to crop production and flooding should be given attention considering the historical and current deforestation trends in Abuan watershed (Van der Ploeg et al. 2011;Hansen et al. 2013). Moreover, results indicated that the water supply of Abuan can irrigate five times more than the actual crop lands (16,824 ha). Such water supply can be tapped for hydropower as renewable energy source. The ES context could spearhead payment schemes like PES or Payment for Ecosystem Services. PES is ideal for Abuan watershed given the ridgeto-reef topography, high conservation priority, and high risk to CC-LULCC-as shown by this study.
These recommendations are timely given that the Philippine Land Use Act is soon to be legalized, providing an opportunity to better plan and manage the Philippine watersheds. As a final thought, we call the watershed stakeholders (national and local government; non-government organizations; upland and lowland locals) to be proactive with the impacts of CC and LULCC.

Conclusions
This study aimed to associate the risks coming from probable streamflow changes caused by climate and land use change using SWAT and SWAT-CUP. We have five concluding remarks to summarize the learnings from this study. First, improvements of the simulations of the baseline model, particularly peak and low flows, were attained by knowing the right choice of parameters to calibrate and why. Second, CC and LULCC were almost equal contributors to the projected average water deficit, but water surplus was more LULCC-driven -worst emission scenarios were RCP 8.5 in 2050 and 2070 at 26% deforestation. Third, the timing of water deficit from May to August (around -18.65%) is a warning on possible cropping schedule changes, and water surplus in November to December (12.79%) may affect the production of rainfed corn lands and indicate higher flood susceptibility. Fourth, at most 1/5 potential croplands are at risk amounting 0.60 MUSD based on the probability of water deficit. As a final note, projection ranges from different climate models can be very pessimistic and may vary from one model to another, therefore their distribution should be part of the uncertainty propagation of streamflow models and the risk estimates should be conservative.

Compliance with ethical standards
Conflict of interest The authors declare no conflict of interest upon doing this study.
Availability of data and models Available upon request to authors.
Access to preprocessing and visualization scripts https://github.com/ arnanaraza/AbuanSWAT Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.