EXPLORING CLIMATE CHANGE EFFECTS ON WATERSHED SEDIMENT YIELD AND LAND COVER-BASED MITIGATION MEASURES USING SWAT MODEL, RS AND GIS: CASE OF CAGAYAN RIVER BASIN, PHILIPPINES

The impact of climate change in the Philippines was examined in the country’s largest basin—the Cagayan River Basin—by predicting its sediment yield for a long period of time. This was done by integrating the Soil and Water Assessment Tool (SWAT) model, Remote Sensing (RS) and Geographic Information System (GIS). A set of Landsat imageries were processed to include an atmospheric correction and a filling procedure for cloud and cloud-shadow infested pixels was used to maximize each downloaded scene for a subsequent land cover classification using Maximum Likelihood classifier. The Shuttle Radar Topography Mission (SRTM)-DEM was used for the digital elevation model (DEM) requirement of the model while ArcGISTM provided the platform for the ArcSWAT extension, for storing data and displaying spatial data. The impact of climate change was assessed by varying air surface temperature and amount of precipitation as predicted in the Intergovernmental Panel on Climate Change (IPCC) scenarios. A Nash-Sutcliff efficiency (NSE) > 0.4 and coefficient of determination (R) > 0.5 for both the calibration and validation of the model showed that SWAT model can realistically simulate the hydrological processes in the study area. The model was then utilized for land cover change and climate change analyses and their influence on sediment yield. Results showed a significant relationship exists among the changes in the climate regime, land cover distributions and sediment yield. Finally, the study suggested land cover distribution that can potentially mitigate the serious negative effects of climate change to a regional watershed’s sediment yield.


INTRODUCTION
The Philippines cannot escape the negative impacts of climate change. The country was tagged as a climate hotspot and vulnerable to some of the worst manifestations of climate change (Jabines & Inventor, 2007). As with other developing countries in Asia, the Philippines is highly subject to natural hazards as exemplified by the 2006 landslide and the havoc wreaked by typhoons Frank, Ondoy andPedring in 2008, 2009 and 2011 respectively. The country is also prone to various hydro-meteorological and geological hazards because of its geographic and geologic setting, threatening the country by the passage of tropical cyclones and occurrences of extreme or prolonged rainfall, strong earthquakes, volcanic eruptions and tsunamis and these hazards will be aggravated and the impact of geological events can be worsened by global warming (Solidum, 2011). Furthermore, climate change threatens the country by increasing the intensity and frequency of storms and droughts. CAD-PAGASA (2004) reported that the country is likely to be adversely affected by climate change since its economy is heavily dependent on agriculture and natural resources. Given these scenarios, it is timely that research pertaining to the impact of climate change to the country be quantitatively assessed.
The study particularly explored the influence of land cover on sediment yield and suggested land cover conversions that can potentially mitigate the serious negative effects of climate change to the sediment yield of a large basin. The study is significant for a proposed watershed management in the country that will incorporate the possible impacts of climate change on sediment yield.

Geographical and Political boundaries
The Cagayan River Basin (CRB) is the largest river basin in the Philippines. It is located in the northeastern portion of the island of Luzon and between 15 0 52'N-18 0 23'N latitudes and 120 0 51'E-122 0 19'E longitudes ( Figure 1). CRB has a drainage area of approximately 27,700 km 2 covering the provinces of Regions 2, Cordillera Autonomous Region (CAR) and small parts of Region 3 (RBCO, 2007).

Climate, Topography and Physiography
CRB falls under Type III climate zone which is characterized by no pronounced maximum rain period and a short dry period (BRS-DPWH, 2002). According to PAGASA (2009), the northern part of the basin has an average annual rainfall of 1,000 mm and 3,000 mm in the southern mountains. The mean annual temperature and average relative humidity are 23.6-26.0 0 C and 75-85%, respectively (DPWH & JICA, 2001).
The area is relatively flat plain but is broken by low rising ridges and hummocks in some places (BRS-DPWH, 2002). Approximately 50% of the area is relatively flat with slope that varies from 0-17%. About 33% of the area has slopes between 17-42% while the rest are with slopes greater than 42% based on a slope map derived from the SRTM-DEM. It is also surrounded by three mountain ranges: Sierra Madre, Cordillera Central and Caraballo-Maparang in the east, west and south. respectively (DPWH & JICA, 2001

MODEL DESCRIPTION
The Soil and Water Assessment Tool (SWAT) model is a river basin scale model specifically used in predicting the effect of land management practices, over long periods of time, on variables such as flow and sediment to areas of varying soils, land use and management conditions. It is physically based, uses readily available inputs, computationally efficient and enables users to study long-term impacts (Neitsch et. al, 2005). The study used the SWAT 2005 version via the ArcSWAT interface for ArcGIS™.
SWAT is a continuous time model and is not designated to simulate detailed single-event flood routing (Neitsch et. al, 2005) and operates on a daily time step (Hao et. al, 2003). To predict surface runoff yield, the model uses a modified version of the SCS CN method (USDA-SCS, 1972): where Q and R are the daily surface runoff and daily rainfall, respectively, both in mm H 2 O. S is a retention parameter which varies spatially under various soil, land use, management and slope conditions, and temporally to respond to changes in soil water content (Hao et. al, 2003). The retention parameter is related to the curve number (CN) and defined as: The model estimates erosion and sediment yield from each subbasin using the Modified Universal Soil Loss Equation (MUSLE) (Williams, 1995): where sed is the sediment yield on a given day in metric tons, Q surf is the surface runoff volume in mmH 2 0/ha, q peak is the peak runoff rate in m 3 /s, area hru is the area of the HRU in ha, K is the USLE soil erodibility factor, C is the USLE cover and management factor, P is the USLE support practice factor, LS is the USLE topographic factor and CFRG is the course fragment factor. For a detailed description of these variables, the reader may refer to the theoretical documentation of SWAT 2005 (Neitsch et. al, 2005).

METHODOLOGY
The general procedure used in the study can be seen in Figure  1. Several datasets are inputted to the model after data preparation. The basin is automatically delineated via the SWAT model interface (ArcSWAT) using the input DEM while subbasins and finer subdivisions in the basin called the hydrologic response units (HRU) are defined by setting threshold limits for land use/land cover, soil type and slope class. Available flow and sediment data were used to calibrate and validate the model. The calibrated model was then rerun for five scenarios, the results of which are compared and became the basis of the author's final analysis.

Digital Elevation Model:
The Shuttle Radar Topography Mission Digital Elevation Model (SRTM-DEM) was used for the DEM requirement of SWAT. SRTM data are products of processed raw radar signals spaced at different intervals at the Jet Propulsion Laboratory (JPL) (USGS, 2011). The DEM used was a 3 arc-second medium resolution elevation data (approx. 90 m) resampled using cubic convolution interpolation and downloaded at the EarthExplorer website (URL: http://edcsns17.cr.usgs.gov/NewEarthExplorer). The DEM was used to generate percent slope values, to do automatic watershed delineation and in defining stream networks and gage outlets.

Land Cover/Land
Use Map was generated from two sets of Landsat 7 TM and ETM+. CRB covers three Landsat scenes within the coverage of path 116, rows 47, 48 and 49. These satellite imageries passed through a processing scheme as shown in Figure 2. All downloaded Landsat ETM+ images were Level 1T products in Geographic Tagged Image-File Format (GeoTIFF). These were geo-referenced to include terrain correction that corrected parallax error from local topographic relief with a digital elevation model (Helmer & Ruefenacht, 2007).
A simple atmospheric correction called the Dark Object Subtraction (DOS) technique was applied to the set of final images (reference images) covering the CRB before they were processed for cloud and cloud shadow masking and filling method developed by Martinuzzi et. al (2006). Image classification was done in per scene basis. The five general land cover classes used in the study are forest, vegetation, water bodies, bare soil and built-up areas. There were four supervised classifiers (Maximum likelihood, Parallelepiped, Minimum Distance and Mahalanobis Distance) and two unsupervised classifiers (ISODATA and K-means) that were tested to classify the reference images. Maximum likelihood was the final classifier used for image classification since it produced the highest over-all accuracy and a kappa coefficient () nearest to a value of 1 (Table 1). The overall accuracy should be at least 85% as prescribed by Anderson et. al (1976) while the value for  is preferably close to 1.00 for this represents a situation where the classification is perfectly superior to random assignment of classes (Gao, 2009). The resulting classified image was post-processed using Majority Analysis in ENVI™ to eliminate "salt and pepper" effects.
The of LULC map using the more recent image set for the Cagayan River Basin is shown in Figure 3. Additionally, Table  2 shows the user-defined LULC classes and their corresponding SWAT LULC classes.

Weather Stations.
Three main weather stations located in Cagayan and Isabela provinces were used in the study. These stations have daily and monthly rainfall, humidity and temperature (minimum and maximum) data obtained from the Philippine Atmospheric, Geophysical and Astronomical Services Administration (PAGASA). Three additional rainfall stations located in Kalinga, Abra and Benguet from PAGASA and two additional precipitation data from stations in the Aurora province were also used from the Weather Underground (http://www.wunderground.com). The eight weather stations are shown in Figure 1. The weather data were obtained in text format and reformatted for SWAT input.

Climate Change Data.
The A1B and A2 climate change scenarios are utilized in the study. Data for these two scenarios are extracted from PAGASA's run of the Providing Regional Climates for Impact Studies (PRECIS) model which generated projected changes in seasonal mean temperature ( 0 C) and rainfall (%) (PAGASA, 2010). These projected incremental changes in rainfall and temperature are inputted to the model by editing the RFINC(mon) and TMPINC(mon), respectively, in the .SUB files of the SWAT model.

Sensitivity Analysis.
SWAT parameters will have to undergo sensitivity analysis first before model calibration to help in identifying and ranking parameters that have significant impact on specific model outputs such as streamflow and sediment yield (Saltelli et. al, 2000). The most sensitive parameters for flow are the Baseflow alpha factor (Alpha_BF), Initial SCS runoff curve number for moisture condition II (CN2) and the threshold depth of water in the shallow aquifer required for return flow to occur (Gwqmn). For sediments, the most sensitive parameters are the linear parameter for calculating the maximum amount of sediment that can be reentrained during channel sediment routing (Spcon), Exponent parameter for calculating sediment reentrained in channel sediment routing (Spexp) and the USLE equation support practice factor (Usle_P). These parameters are given priority in manual adjustments since minor adjustments in their values can translate to significant change in simulated values.

Model Performance Evaluation
The model was evaluated using four quantitative statistics as recommended and used by Moriasi et. al (2007) and Duan et. al (2009). These statistics are the Nash-Sutcliffe efficiency (NSE), percent bias (PBIAS), ratio of the root mean square error to the standard deviation of measured data (RSR) and the coefficient of determination (R 2 ). NSE indicates how well the plot of observed versus simulated values fits the 1:1 line (Alansi et. al, 2009), PBIAS measures the average tendency of the simulated data to be larger or smaller than their observed counterparts (Gupta et. al, 1999), RSR is the ratio of the Root Mean Square Error and the standard deviation of measured data (RMSE) (Moriasi et. al, 2007) and R 2 is an indicator of relationship strength between the observed and simulated values (Alansi et. al, 2009). In general, model simulation can be judged as satisfactory if NSE>0.40 and R2>0.5 (Duan et. al, 2009), and if RSR≤0.70, PBIAS ±25% for streamflow and PBIAS ±55% for sediment (Moriasi et. al, 2007). Tables 4-5 provides a summary of model performance during model calibration and validation in different time steps.

Base Scenario
The soil loss rate map derived from sediment yield for the base scenario is shown in Figure 3. The minimum and maximum values for this scenario are 0.45 ton ha -1 yr -1 and 12.05 ton ha -1 yr -1 which correspond to subbasins 28 and 21, respectively. The simulated maximum value is beyond the upper limit of tolerable soil loss of 11.2 t ha -1 yr -1 according to Hudson (1995) as cited by Alibuyog et. al (2009). Subbasin 28 is characterized by a relatively flat terrain with the whole area having slope less than or equal to 17% and predominantly an agricultural area (AGRR) while most of the subbasin 21 area have steep slopes (>17%) and with almost equal distribution of RNGE, AGRR and FRST areas.

Model Simulation incorporating Climate Change
Climate change data are incorporated in the model by inputting the projected seasonal change in rainfall and temperature for each subbasin. After manipulating the subbasin parameters for climate change analysis (RFINC and TMPINC), the calibrated model was rerun for two scenarios (A1B and A2) each under two time slices centered at year 2020 and 2050. These two scenarios have been the focus of climate change model intercomparison studies according to IPCC (2007). The average total generated sediment yield of the whole basin for this run is shown in Figure 6a.

Model Simulation incorporating Land Use/Land Cover Change
The projected land cover change rates derived from subtracting the sets of Landsat images are inputted to the model by modifying the individual HRU files (with file extension *.hru) for each subbasin. The calibrated model is rerun with the modified HRU files and the average sediment yield for the simulation period is evaluated. The result of this run is also shown in Figure 6a.

Model Simulation incorporating Climate Change and Land Use/ Land Cover Change
Using all the modified files in Sections 5.2 and 5.3, a rerun of

Model Simulation incorporating land cover-based mitigation measures
The study used riparian reforestation and afforestation of hilly and mountainous areas as two land cover-based mitigation measures. Important SWAT input files (e.g., *.mgt, *.sol, *.hru) were modified to reflect these land cover changes. These mitigation measures were applied to the model under the four previously mentioned scenarios in Sections 5.1 to 5.4. Figure 7 shows the spatial distribution of HRUs where these mitigation measures are applied while Figure 6b shows their effect on the generated sediment yield of the basin.

The Effects of LULC and Climate Changes on the Sediment Yield of the Cagayan River Basin
The Cagayan River Basin was found to have a total average annual sediment yield of 114.76 ton ha -1 yr -1 under the base scenario that is, if the present land use/land cover (LULC) distribution of the basin were to use in model simulation. This simulated sediment yield agrees with the reported observed range of average erosion rate for Regions CAR, I and II according to FAO (1998) as cited by Asio et. al (2009). The projected changes in LULC and climatic parameters have produced an increase in sediment yield (+4.5% to +28.8%) compared to the base scenario except for the climate change scenario A2 2020.
The said climate change scenario has produced a decrease (-2.1%) in the basin's sediment yield possibly due to the projected lesser rainfall and less increase in temperature compared to its counterpart for the year 2050 (A2 2050) and scenario A1B. It should be noted, however, that the same scenario will eventually produce an increase in sediment yield if coupled with LULC change (A220LC in Figure 6a).
In general, it has been demonstrated that sediment yield will increase if the current LULC change rate and projected change in the climatic parameters in the Cagayan River basin would persist and climate change scenario A1B would bring a higher increase in sediment yield compared to the A2 scenario.

Effects of Applying Proposed Mitigation Measures
The land cover-based mitigation measures applied to both climate change scenarios with LULC change have produced lesser sediment yields compared to the base scenario ( Figure  6b) by about -26.3% to -45.32%. It is interesting to note that the highest generated sediment yield (157.21 t ha -1 yr -1 under A1B50LC scenario) was decreased to 109. 44 t ha -1 yr -1 which was lower than the current sediment yield of the basin (114.76 t ha -1 yr -1 ) under the base scenario. This demonstrated how riparian reforestation and afforestation of hilly and mountainous areas have successfully mitigated the ill-effects of climate change to the sediment yield of the Cagayan river basin even if coupled with LULC change.

CONCLUSIONS
The study has demonstrated how the integration of SWAT model, Remote Sensing and GIS can be a powerful tool in simulating watershed variables such as the sediment yield of a large river basin. Moreover, the study has validated the applicability of the model in simulating flow and sediment discharge dynamics of the Cagayan river basin based on the satisfactory values of the statistical measures of model efficiency. Lastly, climate change data were successfully utilized to quantify the impact of changes in the climate regime to the sediment yield of the study basin.  Figure 6. Simulated sediment yield of the basin under different scenarios: (a) BS=Base Scenario, LC=land use/land cover change scenario, codes for other scenarios: A1B=climate change scenario A1B, A2= climate change scenario A2, 20=corresponds to year 2020, 50=corresponds to year 2050 (e.g., A1B20LC= climate change scenario A1B at year 2020 with land use/land cover change scenario); (b) Comparison of sediment yields generated before and after applying mitigation measures.