Hydrological Modeling and Runoff Mitigation in an Ungauged Basin of Central Vietnam Using SWAT Model

The A-Luoi district in Thua Thien Hue province of Vietnam is under extreme pressure from natural and anthropogenic factors. The area is ungauged and suffering from data scarcity. To evaluate the water resources availability and water management, we used Soil and Water Assessment Tools (SWAT). A multi-approach technique was used to calibrate the hydrological model. The model was calibrated in three time scales: daily, monthly and yearly by river discharge, actual evapotranspiration (ETa) and crop yield, respectively. The model was calibrated with Nash-Sutcliffe and R2 coefficients greater than 0.7, in daily and monthly scales, respectively. In the yearly scale, the crop yield inside the model was calibrated and validated with Root Mean Square Error (RMSE) less than 2.4 ton/ha. The water resource components were mapped temporally and spatially. The outcomes showed that the highest mean monthly surface runoff, 323 to 369 mm, between September and November, resulted in extreme soil erosion and sedimentation. The monthly average of actual evapotranspiration was the highest in May and lowest in December. Furthermore, installing “Best Management Practices” (BMPs) reduced surface runoff in agricultural lands. However, using event-based hydrological and hydraulic models in the prediction and simulation of flooding events is recommended in further studies.


Introduction
Hydrological modeling of water cycle in areas with extreme events and natural hazards (e.g., flooding, droughts) is imperative for sustainable management of soil and water resources.Understanding water resources availability would help stakeholders and policymakers to plan and develop an area.There are various hydrological models that can estimate water resource availability (e.g., lumped models, physical distributed and semi-distributed models, empirical models, statistical models).Among these, semi-distributed hydrological models can simulate water balance spatially based on various soils, land uses, topography and climate conditions.One of these hydrological models is SWAT [1] which is tested in various world climates from arid and semi-arid regions [2] to humid and tropical areas [3].Moreover, it is able to simulate water resources in large scales to regional scales.For example, Schuol et al. [4] estimated blue and green water availability in the African continent.Phuong et al. [5] evaluated the surface runoff and soil erosion in a regional area in Vietnam.One of the concerns in hydrological modeling is the uncertainty of parameters; therefore, Qiao et al. [6] used a multi-objective function to reduce this uncertainty in the SWAT model.Surface runoff, soil erosion and flooding are the main issues in humid and tropical areas due to heavy rainfall [7].Various researchers have studied climate variabilities and human activities on soil erosion, water quality and hydrological process [8,9].
Vietnam's tropical region is very prone to soil erosion and surface runoff.Approximately 40% of natural areas are at risk of erosion [10].In recent years, most of the forested areas, particularly in Central Vietnam, were converted into agricultural lands, which resulted in surface runoff and soil erosion [11].Additionally, the hydrological behavior of watersheds in Vietnam is changing due to anthropogenic factors.The evidence of this comes from rising of soil erosion, surface runoff and floods in Central Vietnam [5].There is no comprehensive framework to mitigate the environmental hazards in this area; therefore, we created a hydrological model of Central Vietnam with a case study in the area of A-Luoi district in Hue province.
The A-Luoi district in Central Vietnam is under high socio-economic pressure, which resulted in various natural hazards such as flooding and soil erosion.Lack of data is one of the main obstacles for policymakers in this area.The main goal of this research is to set up a hydrological model in the A-Luoi area as a representative case study of Central Vietnam and to simulate the water balance components (e.g., water yield, soil moisture, actual evapotranspiration) in order to understand the status of water resources in this area.Assessments of surface runoff, and soil erosion are the further goals of this study.We calibrated the hydrological model using multiple calibration approaches by satellite images and regionalization technique.Finally, we presented the best management practices adapted to A-Luoi area to reduce surface runoff and soil erosion.

Study Area
The study area is a part of A-Luoi district, which is located in the center of Vietnam.The A-Luoi district consists of three watersheds.The first watershed in the west of the A-Luoi, where the main river, called the "Xe Khong River," flows toward Laos.This watershed has an area of approximately 542 km 2 .The second watershed is situated in the middle of A-Luoi and the main river, Song Sia River, flows into the "Tam Giang lagoon".The watershed has an area of approximately 410 km 2 .The third watershed in the east of A-Luoi district has an area of approximately 506 km 2 .The river in this watershed, called "Song Huu Trach River," flows into the "Song Huong lagoon" (Figure 1).We selected the second watershed as our study area, while it is the main part of A-Luoi district.
The annual rainfall in this watershed varies from 2228 mm to 5495 mm (long term climate data).According to the Soil Taxonomic classification [12], Ferralic Acrisols is the dominant soil type in this area.There are Humic Acrisols and Lithic Leptosols types in the northwestern and southeastern part of A-Luoi as well.The watershed has diverse elevations from 19 to 1777 meters above sea level from the north to the west of the basin.Figure 1 shows the study area in A-Luoi district including climate and hydrometric stations around the district.

Data Acquisition
There are four climate stations, three stations outside the boundary (Khe Sanh, Hue and Nam dong stations) with the monthly data available and one station inside the boundary with the daily data from 2005 to 2013 (Aluoi station) (Figure 1).We used daily climate data (e.g., daily rainfall, maximum and minimum temperature, solar radiation, and relative humidity) for the model inputs, and monthly data for the Weather Generator (WGN) parameters.There are no hydrometric stations inside or at the outlet of the A-Luoi district.The closest hydrometric station to A-Luoi, called "Thuong Nhat," is located east of the study area that used for parameterization in A-Luoi watershed.The FAO (Food and Agriculture Organization of the United Nations)-soil map [13] was improved based on local soil data and was also used in this study.This soil map contains information for two soil depths: in topsoil (0-30 cm) and subsoil (30-100 cm).Land use statistics in 2010 was obtained from Department of Natural Resources and Environment in Thua Thien Hue Province.The land use The Nam Dong watershed is a nearby watershed to A-Luoi with measured discharge data, which used for regionalization approach.

Data Acquisition
There are four climate stations, three stations outside the boundary (Khe Sanh, Hue and Nam dong stations) with the monthly data available and one station inside the boundary with the daily data from 2005 to 2013 (Aluoi station) (Figure 1).We used daily climate data (e.g., daily rainfall, maximum and minimum temperature, solar radiation, and relative humidity) for the model inputs, and monthly data for the Weather Generator (WGN) parameters.There are no hydrometric stations inside or at the outlet of the A-Luoi district.The closest hydrometric station to A-Luoi, called "Thuong Nhat," is located east of the study area that used for parameterization in A-Luoi watershed.The FAO (Food and Agriculture Organization of the United Nations)-soil map [13] was improved based on local soil data and was also used in this study.This soil map contains information for two soil depths: in topsoil (0-30 cm) and subsoil (30-100 cm).Land use statistics in 2010 was obtained from Department of Natural Resources and Environment in Thua Thien Hue Province.The land use statistics contains different classes such as protected forest area, evergreen forests, bare lands, residential area, paddy rice lands and other agricultural lands.Forest is dominant land use in the study area and more than 90% of the area coved by evergreen forests.Rice is dominant agriculture crops and cultivated two times per year: Summer-Autumn and Winter-Spring periods (Figure 2).The surface water is the sources of irrigation/inundation for rice.Digital Elevation model (DEM) was derived from Shuttle Radar Topography Mission (SRTM) with 30m resolution.Crop schedule and yield data [14] were implemented in the model in order to simulate crop growth in SWAT.Moderate Resolution Imaging Spectroradiometer (MODIS) time series data (i.e., MOD16 product, including potential and actual evapotranspiration) were used to extract monthly evapotranspiration.These data were further used to calibrate the hydrological model.statistics contains different classes such as protected forest area, evergreen forests, bare lands, residential area, paddy rice lands and other agricultural lands.Forest is dominant land use in the study area and more than 90% of the area coved by evergreen forests.Rice is dominant agriculture crops and cultivated two times per year: Summer-Autumn and Winter-Spring periods (Figure 2).The surface water is the sources of irrigation/inundation for rice.Digital Elevation model (DEM) was derived from Shuttle Radar Topography Mission (SRTM) with 30m resolution.Crop schedule and yield data [14] were implemented in the model in order to simulate crop growth in SWAT.Moderate Resolution Imaging Spectroradiometer (MODIS) time series data (i.e., MOD16 product, including potential and actual evapotranspiration) were used to extract monthly evapotranspiration.These data were further used to calibrate the hydrological model.

Model Configuration
SWAT is able to predict the impact of land management practices on soil and water.Additionally, SWAT can simulate water quality and quantity as well as crop growth.The main components of the SWAT model consist of weather, hydrology, plant growth, nutrient, pesticide, bacteria and land management [1].In this study, we focused only on the hydrologic, plant growth and land management components of the SWAT model.In SWAT, a watershed is partitioned into sub-basins using DEM, which are then further subdivided into HRUs (Hydrological Response Units) with homogenous soil, land use and slope characteristics.HRUs are the basis of water balance calculation.The water in each HRU can be stored in four levels: snow, soil profile (0-2 m), shallow aquifer (0-2 m) and deep aquifer (2-20 m).Rafiei Emam et al. [15] revealed that a single HRU could not show the characteristics of a sub-basin.Therefore, we used multiple HRUs for sub-basin analysis in this research.In SWAT, a water balance calculation is simulated in two separate components: land phase and routing phase.The land phase controls the amount of water, sediment, nutrient and pesticide loadings to the main channel in each sub-basin, while the routing phase defines the movement of water, sediment and nutrient loadings from the channel networks to the outlet.The water balance is calculated in the land phase using the following Equation (1): where is the soil water content, is the initial soil water content, is the amount of precipitation, indicates surface runoff, shows actual evapotranspiration, is the amount of percolation and shows the return flow.In this study, surface runoff was simulated using a Soil Conservation Service Curve Number (SCS-CN).The runoff curve number is an empirical parameter used for predicting direct runoff or infiltration from rainfall excess.The curve number method was developed by the USDA Natural Resources Conservation Service [16].SCS is widely used and is an efficient method for determining the approximate amount of direct runoff from a rainfall event in a particular area.It is based on the area's hydrologic soil group, land use, treatment and hydrologic condition.Percolation was

Model Configuration
SWAT is able to predict the impact of land management practices on soil and water.Additionally, SWAT can simulate water quality and quantity as well as crop growth.The main components of the SWAT model consist of weather, hydrology, plant growth, nutrient, pesticide, bacteria and land management [1].In this study, we focused only on the hydrologic, plant growth and land management components of the SWAT model.In SWAT, a watershed is partitioned into sub-basins using DEM, which are then further subdivided into HRUs (Hydrological Response Units) with homogenous soil, land use and slope characteristics.HRUs are the basis of water balance calculation.The water in each HRU can be stored in four levels: snow, soil profile (0-2 m), shallow aquifer (0-2 m) and deep aquifer (2-20 m).Rafiei Emam et al. [15] revealed that a single HRU could not show the characteristics of a sub-basin.Therefore, we used multiple HRUs for sub-basin analysis in this research.In SWAT, a water balance calculation is simulated in two separate components: land phase and routing phase.The land phase controls the amount of water, sediment, nutrient and pesticide loadings to the main channel in each sub-basin, while the routing phase defines the movement of water, sediment and nutrient loadings from the channel networks to the outlet.The water balance is calculated in the land phase using the following Equation (1): where SW t is the soil water content, SW 0 is the initial soil water content, R day is the amount of precipitation, Q sur f indicates surface runoff, ET a shows actual evapotranspiration, W seep is the amount of percolation and Q gw shows the return flow.In this study, surface runoff was simulated using a Soil Conservation Service Curve Number (SCS-CN).The runoff curve number is an empirical parameter used for predicting direct runoff or infiltration from rainfall excess.The curve number method was developed by the USDA Natural Resources Conservation Service [16].SCS is widely used and is an efficient method for determining the approximate amount of direct runoff from a rainfall event in a particular area.It is based on the area's hydrologic soil group, land use, treatment and hydrologic condition.Percolation was simulated with a layered storage routing technique combined with a crack flow model; and potential evapotranspiration (PET) was simulated using the Penman-Monteith method.The actual evapotranspiration was predicted based on PET in addition to soil and crop factors based on the methodology developed by Ritchie [17].The daily value of the Leaf Area Index (LAI) was used to divide the PET into potential soil evaporation and potential plant transpiration.LAI and root development were predicted by the crop growth component of SWAT.In SWAT, soil erosion is calculated with MUSLE algorithm (Modified Universal Soil Loss Equation) [18], which is a modified version of USLE (Universal Soil Loss Equation) [19].USLE is a function of rainfall energy, while MUSLE is a function of runoff.
To better understand the hydrological behavior, which leads to a better calibration of the model, we used an automated digital filter program [20] to separate discharge into base flow and surface runoff.Using this approach, the low-frequency base flow component was separated from high-frequency runoff components through a low-pass filter [21].These data further were used to calibrate the hydrological model.Meaurio et al. [22] used the same method to calculate a base flow constant.More details are presented by Arnold and Allen [23].In this research, watershed was delineated into 20 sub-basins using an SRTM map and stream network in an ArcGIS environment.The area of watershed is approximately 41,000 ha.Some 54 HRUs were created by integrating land use, soil and slope maps.Because A-Luoi is a mountainous area, we estimated precipitation (PLAPS) and temperature (TLAPS) difference by elevation.PLAPS is precipitation lapse rate (mm H 2 O/km), and TLAPS refers to temperature lapse rate ( • C/km).For this aim, all sub-basins with a difference in elevation greater than 100 m were selected and PLAPS and TLAPS were applied.The PLAPS and TLAPS were estimated at 243.2 mm and −5.8 • C per km using climate and elevation data.Rice is the dominant agriculture crop in A-Luoi; therefore, the crop schedule for rice was applied to the model.

Calibration, Validation and Uncertainty Analysis
The A-Luoi district is ungauged; therefore, to calibrate the models we used different approaches.The model was calibrated and validated by the daily river discharge, yearly crop yield and monthly evapotranspiration during January 2006 and December 2010, with the three-year warm up period based on data in 2005, repeated three times.Figure 3 shows the work plan of the calibration process.
simulated with a layered storage routing technique combined with a crack flow model; and potential evapotranspiration (PET) was simulated using the Penman-Monteith method.The actual evapotranspiration was predicted based on PET in addition to soil and crop factors based on the methodology developed by Ritchie [17].The daily value of the Leaf Area Index (LAI) was used to divide the PET into potential soil evaporation and potential plant transpiration.LAI and root development were predicted by the crop growth component of SWAT.In SWAT, soil erosion is calculated with MUSLE algorithm (Modified Universal Soil Loss Equation) [18], which is a modified version of USLE (Universal Soil Loss Equation) [19].USLE is a function of rainfall energy, while MUSLE is a function of runoff.
To better understand the hydrological behavior, which leads to a better calibration of the model, we used an automated digital filter program [20] to separate discharge into base flow and surface runoff.Using this approach, the low-frequency base flow component was separated from highfrequency runoff components through a low-pass filter [21].These data further were used to calibrate the hydrological model.Meaurio et al. [22] used the same method to calculate a base flow constant.More details are presented by Arnold and Allen [23].In this research, watershed was delineated into 20 sub-basins using an SRTM map and stream network in an ArcGIS environment.The area of watershed is approximately 41,000 ha.Some 54 HRUs were created by integrating land use, soil and slope maps.Because A-Luoi is a mountainous area, we estimated precipitation (PLAPS) and temperature (TLAPS) difference by elevation.PLAPS is precipitation lapse rate (mm H2O/km), and TLAPS refers to temperature lapse rate (°C/km).For this aim, all sub-basins with a difference in elevation greater than 100 m were selected and PLAPS and TLAPS were applied.The PLAPS and TLAPS were estimated at 243.2 mm and −5.8 °C per km using climate and elevation data.Rice is the dominant agriculture crop in A-Luoi; therefore, the crop schedule for rice was applied to the model.

Calibration, Validation and Uncertainty Analysis
The A-Luoi district is ungauged; therefore, to calibrate the models we used different approaches.The model was calibrated and validated by the daily river discharge, yearly crop yield and monthly evapotranspiration during January 2006 and December 2010, with the three-year warm up period based on data in 2005, repeated three times.Figure 3 shows the work plan of the calibration process.

Regionalization Approach
The hydrological model's parameters for ungauged basins (i.e., without any river discharge measurements) can be predicted using the regionalization method [24,25].Usually watersheds with similar characteristics show similar hydrological treatments and therefore the hydrological parameters can be transferred from the same watersheds.This transformation is called regionalization.A number of basin characteristics (e.g., soil, topography, land use, precipitation, and temperature) were presented in Samaniego and Bardossy [26].There are different methods for this transformation including regression methods [27], spatial proximity [28] and physical similarity [29].Parajka et al. [30] used geostatistical methods (i.e., Kriging) to transfer model parameters.Sun et al. [31] used a new approach to calibrate hydrological parameters in ungauged basins.They used the river water level obtained from satellite radar altimetry observations at the outlet of the basin instead of river discharge.In another work, Sun et al. [32] used the river flow width obtained through satellite observation as a surrogate for river discharge records to calibrate a hydrological model.
Because there is not a river discharge station at the outlet of our study area and we cannot calibrate the relevant parameters directly, we used regionalization approach to transfer the relevant parameters from the nearby catchment to our study area.To do that, Thuong Nhat catchment, 45 km away from the center of A-Luoi, with the observed discharge data was selected.Data is available at daily scale from 2005 to 2010 at Nhat station in this catchment.The characteristics of the two basins are similar with difference annual precipitation of 2.6%.More physical characteristics are presented in Table 1.We set up a new model in this basin and the model was calibrated based on the river discharge data.The calibrated parameters were further transferred into our original model of A-Luoi basin.The SUFI-2 algorithm (Sequential Uncertainty Fitting, ver.2) [33] in the SWAT-CUP package was used to calibrate and analyze uncertainty.Nash-Sutcliffe coefficient (Equation ( 2)) was computed as the objective function for river discharge in a gauged basin: where Q obs and Q sim are the measured and simulated data, respectively, n is the total number of data records and Q obs is the mean measured data.
In SUFI-2 all sources of uncertainty (e.g., uncertainty in input data, conceptual model and parameters) are represented in model output.The degree of uncertainty is measured by two factors: P-factor, which represents the percentage of measured data bracketed by 95% prediction uncertainty (95PPU), and R-factor, which is the ratio of average width of 95PPU to standard deviation of measured data.The model attempts to bracket most of the data (P-factor = 1) within a narrow band (R-factor = 0).However, often a balance must be reached between these two factors.

Calibration Based on Crop Yield
We calibrated a crop yield because calibration using a crop yield gives more confidence on actual evapotranspiration, and hence give the precise estimation of other water components such as, soil moisture, surface runoff and percolation [15].The SUFI-2 algorithm was used to calibrate the crop parameters and uncertainty analysis.A mean square error was used as the objective function.There are several parameters that affect crop yield.These parameters include the heat unit (HU), harvest index (HI), and bio target (BIO-TARG).The HU is the amount of heat that brings plants to maturity, the HI is the fraction of above-ground dry biomass to yield crops, and the BIO-TARG is the potential plant yield.We optimized these parameters in order to calibrate crop yield in our hydrological model.
Plant growth in the model is usually inhibited by temperature, water, nitrogen and phosphorous stress factors.Rafiei Emam et al. [34] evaluated the water and temperature stresses of winter wheat using a climate change model.Temperature is a very important factor for plant growth; plants can grow only if the mean temperature exceeds the base temperature of plant.The accumulation of daily mean air temperatures above the plant's optimum temperature over the plant-growth period is called a heat unit.The base (or optimum) temperature varies for different plants, for instance, paddy rice has base temperature of 10 • C.

Calibration Based on ETa Using MODIS Products
Because agriculture covers only a small part of A-Luoi and most of the lands in the study area are covered by forest, we used ETa products of MODIS in order to calibrate the SWAT model.There have been several attempts to use satellite observation of evapotranspiration to calibrate a hydrological model.Immerzeel and Droogers, [35] used MODIS satellite images to calibrate the hydrological model using evapotranspiration data.Cheema et al. [36] calibrated the SWAT model in the irrigated Indus basin based on spatially distributed ET derived from remotely sensed data.
Nash-Sutcliff efficiency (NSE) and coefficient of determination were used to assess the SWAT model based on ETa calibration.NSE was used as the objective function to compare the monthly observed and simulated ETa.
We used MOD16 products [37] of MODIS from 2006 to 2013 on a monthly scale.We downloaded 90 images (6 images are missing) in Central Vietnam (h = 28 and v = 7).We further extracted PET and ETa of MOD16 products in ArcGIS environment in our study area and summarized data in each sub-basin.We analyzed the reliability of the data and finally these data were used to calibrate the SWAT hydrologic model.

Regionalization Outcome
The results of calibration and validation of Nam Dong watershed, a nearby watershed to A-Luoi, showed a good performance of the model with high NS (0.74, 0.51) and R 2 (0.71, 0.88) both for calibration and validation, respectively.The sensitivity analysis shows that the curve number (CN) is the most sensitive parameter followed by Base Flow Alpha Factor for Bank Storage (Alpha-BNK) and effective hydraulic conductivity in the main channel (CH-K2).The importance of these parameters has been reported by some other researchers.Rafiei Emam et al. [2,38] mentioned that CN is the most sensitive parameter in mountainous area.Luo et al. [39] revealed that the parameters of CH-K2 and Alpha-BNK play important roles in the calibration of the SWAT model in Tropical areas.
Comparing observed with simulated discharge shows that there is a good agreement between the simulated daily stream flow and the observed data.Additionally, the scatter plot between observed and simulated data for both calibration and validation show the pretty high performance of the model (Figure 4).Therefore, the calibrated parameters were transferred into the A-Luoi model for further processing.

Parameterization for Crop Yield
The sub-basins with paddy rice cultivations (i.e., sub-basins located south of the A-Luoi area) were selected for the parameterization of the model based on crop yield.The calibration based on crop makes the model results more confident [15].There is a direct relationship between crop yield and evapotranspiration; therefore, the calibration of the model based on yield improves the simulation of other water components.
The sensitivity analysis showed that all crop parameters (e.g., HEAT-UNIT, BIO-TARG, HI) were sensitive to crop yield. Figure 5 shows the results of calibration and validation of the paddy rice yield.As illustrated, observed yields are inside or quite close to simulated yield band, indicating good results.The yield varies between 26.5 and 30 ton/ha, with the highest yield found in the year 2010.The RMSE for yield was 2.1 kg•ha −1 and 2.39 kg•ha −1 in the calibration and validation period, respectively, which shows pretty good accuracy of calibration.Akhavan et al. [40] revealed the high range of RMSE from 80 to 4220 kg•ha −1 for calibration and validation of wheat yield in Hamedan, Iran.They mentioned that the large RMSE could be due to the lack of data concerning management practices, which can be accounted in our study area as well.The P-factor is quite good, with a value larger than 0.8 for calibration and 0.89 for the validation period.The large R-factor (i.e., R-factor = 3) represents a large uncertainty.The reason for this amount of uncertainty is probably due to low data availability and insufficient accounting for agriculture and industrial water use in the model.The other reason is due to Nitrogen (N) stress, because the N is removed from the farmlands quickly by rising rainfall and surface runoff.

Parameterization for Crop Yield
The sub-basins with paddy rice cultivations (i.e., sub-basins located south of the A-Luoi area) were selected for the parameterization of the model based on crop yield.The calibration based on crop makes the model results more confident [15].There is a direct relationship between crop yield and evapotranspiration; therefore, the calibration of the model based on yield improves the simulation of other water components.
The sensitivity analysis showed that all crop parameters (e.g., HEAT-UNIT, BIO-TARG, HI) were sensitive to crop yield. Figure 5 shows the results of calibration and validation of the paddy rice yield.As illustrated, observed yields are inside or quite close to simulated yield band, indicating good results.The yield varies between 26.5 and 30 ton/ha, with the highest yield found in the year 2010.The RMSE for yield was 2.1 kg•ha −1 and 2.39 kg•ha −1 in the calibration and validation period, respectively, which shows pretty good accuracy of calibration.Akhavan et al. [40] revealed the high range of RMSE from 80 to 4220 kg•ha −1 for calibration and validation of wheat yield in Hamedan, Iran.They mentioned that the large RMSE could be due to the lack of data concerning management practices, which can be accounted in our study area as well.The P-factor is quite good, with a value larger than 0.8 for calibration and 0.89 for the validation period.The large R-factor (i.e., R-factor = 3) represents a large uncertainty.The reason for this amount of uncertainty is probably due to low data availability and insufficient accounting for agriculture and industrial water use in the model.The other reason is due to Nitrogen (N) stress, because the N is removed from the farmlands quickly by rising rainfall and surface runoff.

Reliability of MODIS Data
To calibrate the model based on the evapotranspiration, first we evaluated the reliability of the MOD16 of MODIS time series data.Toward this aim, the PET of MOD16 was extracted monthly from January 2006 to December 2013.The average PET-MODIS in the basin was plotted against PET calculated by measuring data in the SWAT model using the Penman-Monteith method (Figure 6).As Figure 6 shows, there is a high correlation between PET of MOD16 and PET calculated by SWAT.The dynamics of time series of SWAT and PET-MODIS was mapped in Figure 7.The dynamics of PET of SWAT and MODIS are very similar to each other, but there is a bias between these data.The amount of PET in winter, when there is less precipitation, is more reliable than the other months.The peak of PET-MODIS in summer, May to August, is higher than measured data (approximately 10 to 70 percent more than measured one), which makes it unreliable in these months.Therefore, we calculated the bias between two data time series.Toward this aim, we divided the year into two periods: first from October to March and second from April to September.Then we analyzed data based on rational analysis and found two constants for these two periods.To correct the bias, for the first period, the PET data was multiplied by 0.45, and for second period, the PET data was multiplied by 0.70.Then, the corrected PET of MODIS (PETmc) multiplied to the ratio of PET/ETa ( ) of MODIS in order to bias correction of ETa.

=
(3) where PETm is the potential evapotranspiration of MODIS, ETam is the actual evapotranspiration extracted from MODIS, ETamc is the bias corrected of ETa-MODIS, and PETmc is the bias corrected PET-MODIS results.

Reliability of MODIS Data
To calibrate the model based on the evapotranspiration, first we evaluated the reliability of the MOD16 of MODIS time series data.Toward this aim, the PET of MOD16 was extracted monthly from January 2006 to December 2013.The average PET-MODIS in the basin was plotted against PET calculated by measuring data in the SWAT model using the Penman-Monteith method (Figure 6).As Figure 6 shows, there is a high correlation between PET of MOD16 and PET calculated by SWAT.The dynamics of time series of SWAT and PET-MODIS was mapped in Figure 7.The dynamics of PET of SWAT and MODIS are very similar to each other, but there is a bias between these data.The amount of PET in winter, when there is less precipitation, is more reliable than the other months.The peak of PET-MODIS in summer, May to August, is higher than measured data (approximately 10 to 70 percent more than measured one), which makes it unreliable in these months.Therefore, we calculated the bias between two data time series.Toward this aim, we divided the year into two periods: first from October to March and second from April to September.Then we analyzed data based on rational analysis and found two constants for these two periods.To correct the bias, for the first period, the PET data was multiplied by 0.45, and for second period, the PET data was multiplied by 0.70.Then, the corrected PET of MODIS (PETm c ) multiplied to the ratio of PET/ETa (ρ) of MODIS in order to bias correction of ETa.
where PETm is the potential evapotranspiration of MODIS, ETam is the actual evapotranspiration extracted from MODIS, ETam c is the bias corrected of ETa-MODIS, and PETm c is the bias corrected PET-MODIS results.

Reliability of MODIS Data
To calibrate the model based on the evapotranspiration, first we evaluated the reliability of the MOD16 of MODIS time series data.Toward this aim, the PET of MOD16 was extracted monthly from January 2006 to December 2013.The average PET-MODIS in the basin was plotted against PET calculated by measuring data in the SWAT model using the Penman-Monteith method (Figure 6).As Figure 6 shows, there is a high correlation between PET of MOD16 and PET calculated by SWAT.The dynamics of time series of SWAT and PET-MODIS was mapped in Figure 7.The dynamics of PET of SWAT and MODIS are very similar to each other, but there is a bias between these data.The amount of PET in winter, when there is less precipitation, is more reliable than the other months.The peak of PET-MODIS in summer, May to August, is higher than measured data (approximately 10 to 70 percent more than measured one), which makes it unreliable in these months.Therefore, we calculated the bias between two data time series.Toward this aim, we divided the year into two periods: first from October to March and second from April to September.Then we analyzed data based on rational analysis and found two constants for these two periods.To correct the bias, for the first period, the PET data was multiplied by 0.45, and for second period, the PET data was multiplied by 0.70.Then, the corrected PET of MODIS (PETmc) multiplied to the ratio of PET/ETa ( ) of MODIS in order to bias correction of ETa.

=
(3) where PETm is the potential evapotranspiration of MODIS, ETam is the actual evapotranspiration extracted from MODIS, ETamc is the bias corrected of ETa-MODIS, and PETmc is the bias corrected PET-MODIS results.

Procedure of Model Calibration Using ET-MODIS
In this step, the SWAT Model was calibrated in the sub-basins with dominant forest cover.To calibrate the SWAT model by ETa data, first we kept the calibrated parameters for river discharge and crop yield and then we chose the parameters that affected ETa (Table 2).Five parameters were used for sensitivity analysis using a one-at-a-time method [33].Comparing PET extracted from MODIS-MOD16 with one measured using the Penman-Monthieh method.For bias correction, MODIS data were multiplied by 0.45 for the period of October-March and 0.70 for April-September.

Procedure of Model Calibration Using ET-MODIS
In this step, the SWAT Model was calibrated in the sub-basins with dominant forest cover.
To calibrate the SWAT model by ETa data, first we kept the calibrated parameters for river discharge and crop yield and then we chose the parameters that affected ETa (Table 2).Five parameters were used for sensitivity analysis using a one-at-a-time method [33].The results of the global sensitivity analysis are shown by p-value and t-test, which show the measure of sensitivity and the significance of sensitivity, respectively.The results show that BLAI (t-test = [2.5],p-value = 0.05) is the most sensitive parameter for ET calibration followed by DLAI (t-test = [1.5],p-value = 0.3) and ESCO (t-test = [0.6],p-value = 0.5).BLAI represent the maximum potential leaf area index of trees.BLAI is one of the six parameters used to determine leaf area development of plants during the growing season.DLAI is the "fraction of the growing season when the leaf area begins to decline" [41].ESCO depends on soil characteristics and controls the soil evaporation.
The results of the model calibration were satisfactory.Figure 8a shows the results of the ETa calibration in the forest areaof watershed.As Figure 8 depicts, the observed data are captured by the predicted band (95PPU).In all of the sub-basins, the temporal changes of ETa were simulated quite well, with a Nash-Sutcliffe (NS) and coefficient of determination (R 2 ) higher than 0.79.Moriasi et al. [42] and Santhi et al. [43] reported that if the NS and R 2 are greater than 0.5, then the performance of the model is satisfactory.The results also show that the simulated ETa has variation in different seasons.Petković et al. [44] mentioned that relative humidity and a wind speed of two meters are two factors that play a role in the variation of ET.The same results were also presented by Nguyen and Kappas [3].As Figure 8a,b shows, the lowest evapotranspiration is from November to February, when the average temperatures are lower than over the other months.
* r parameter value is multiplied by 1 + given value, v parameter value is replaced by a value from the given range.
The results of the global sensitivity analysis are shown by p-value and t-test, which show the measure of sensitivity and the significance of sensitivity, respectively.The results show that BLAI (ttest = [2.5],p-value = 0.05) is the most sensitive parameter for ET calibration followed by DLAI (t-test = [1.5],p-value = 0.3) and ESCO (t-test = [0.6],p-value = 0.5).BLAI represent the maximum potential leaf area index of trees.BLAI is one of the six parameters used to determine leaf area development of plants during the growing season.DLAI is the "fraction of the growing season when the leaf area begins to decline" [41].ESCO depends on soil characteristics and controls the soil evaporation.
The results of the model calibration were satisfactory.Figure 8a shows the results of the ETa calibration in the forest areaof watershed.As Figure 8 depicts, the observed data are captured by the predicted band (95PPU).In all of the sub-basins, the temporal changes of ETa were simulated quite well, with a Nash-Sutcliffe (NS) and coefficient of determination (R 2 ) higher than 0.79.Moriasi et al. [42] and Santhi et al. [43] reported that if the NS and R 2 are greater than 0.5, then the performance of the model is satisfactory.The results also show that the simulated ETa has variation in different seasons.Petković et al. [44] mentioned that relative humidity and a wind speed of two meters are two factors that play a role in the variation of ET.The same results were also presented by Nguyen and Kappas [3].As Figure 8a,b shows, the lowest evapotranspiration is from November to February, when the average temperatures are lower than over the other months.

Quantification of Water Components
After the model was calibrated by river discharge, crop yield and ET, we computed the water resource availability in the watershed.Figure 9 shows the water resources availability including the monthly (for the time period 2006 to 2013) water yield (i.e., the net amount of water that leaves the sub-basin and contributes to streamflow in the reach during the time period), green water resources (i.e., soil water content) and green water flow (i.e., actual evapotranspiration) in the A-Luoi watershed of Vietnam.

Quantification of Water Components
After the model was calibrated by river discharge, crop yield and ET, we computed the water resource availability in the watershed.Figure 9 shows the water resources availability including the monthly (for the time period 2006 to 2013) water yield (i.e., the net amount of water that leaves the sub-basin and contributes to streamflow in the reach during the time period), green water resources (i.e., soil water content) and green water flow (i.e., actual evapotranspiration) in the A-Luoi watershed of Vietnam.There is a significant spatial variation in the hydrological components across the watershed.The distributed maps showed that in the north of basin, where the rainfall and water yield are small, and actual evapotranspiration (ETa) is high.The ETa shows the amount of water that is consumed by plants.The annual average ETa (M95PPU) during the period of 2006 to 2013 varied from 719 mm/year to 838 mm/year in different parts of the watershed.The lowest amount of ETa was found in the southern area, while the northern area had the highest ETa (Figure 9a).The soil water content map (Figure 9b) shows the area where rain-fed Rubber and other perennial crops have a better chance of high production due to soil moisture.The southern area has the highest soil moisture, while the eastern area has the lowest soil moisture.This is probably due to soil characteristics in different parts of the watershed.Water yields (Figure 9c) varied from 1679 mm/year to 2222 mm/year in different sub-basins from 2006 to 2013.Water yield is the net amount of water that leaves the sub-basin and contributes to streamflow in the reach during the time period.
The annual average rainfall from 2006 to 2013 was 4049 mm/year.We calculated the monthly average water components in the entire basin.The results of average monthly evapotranspiration showed that ETa is highest in May, due to the high demand of water by forests and plants in this There is a significant spatial variation in the hydrological components across the watershed.The distributed maps showed that in the north of basin, where the rainfall and water yield are small, and actual evapotranspiration (ETa) is high.The ETa shows the amount of water that is consumed by plants.The annual average ETa (M95PPU) during the period of 2006 to 2013 varied from 719 mm/year to 838 mm/year in different parts of the watershed.The lowest amount of ETa was found in the southern area, while the northern area had the highest ETa (Figure 9a).The soil water content map (Figure 9b) shows the area where rain-fed Rubber and other perennial crops have a better chance of high production due to soil moisture.The southern area has the highest soil moisture, while the eastern area has the lowest soil moisture.This is probably due to soil characteristics in different parts of the watershed.Water yields (Figure 9c) varied from 1679 mm/year to 2222 mm/year in different sub-basins from 2006 to 2013.Water yield is the net amount of water that leaves the sub-basin and contributes to streamflow in the reach during the time period.
The annual average rainfall from 2006 to 2013 was 4049 mm/year.We calculated the monthly average water components in the entire basin.The results of average monthly evapotranspiration showed that ETa is highest in May, due to the high demand of water by forests and plants in this month.The prediction uncertainty in summer is higher than in winter due to the high consumption of water by plants during these months (Figure 10a).Paddy rice lands affect the ETa considerably.The total ratio of ETa to rainfall for eight years (2006-2013) was calculated at 36%, showing that the ETa is always less than rainfall in tropical areas.Nevertheless, this ratio is lowest in October and highest in January.
month.The prediction uncertainty in summer is higher than in winter due to the high consumption of water by plants during these months (Figure 10a).Paddy rice lands affect the ETa considerably.The total ratio of ETa to rainfall for eight years (2006-2013) was calculated at 36%, showing that the ETa is always less than rainfall in tropical areas.Nevertheless, this ratio is lowest in October and highest in January.

Surface Runoff
The result of the SWAT hydrological model shows that the surface runoff is high in upstream areas of basins.Vegetation cover, soil type and slope gradient [45] are the important factors affecting surface runoff.Land use change and deforestation lead to increasing surface runoff and hence soil erosion.The results showed that the spatial average of surface runoff in the watershed varied from 915 mm/year to 1072 mm/year from 2006 to 2013.The area upstream has the highest amount of surface runoff probably due to its high slope.
The monthly average 95PPU of surface runoff was extracted in order to understand the temporal distribution and uncertainty of surface runoff (Figure 10b).This figure reveals that October has the highest surface runoff, with an average of approximately 368 mm/month.The surface runoff during September to November has the highest uncertainty, probably due to a conceptual model and parameter uncertainties.However, all uncertainties are depicted in the 95PPU band, shown in gray in the figure.As the Figure 10b shows, there is a strong relationship between surface runoff and rainfall.In February, the rainfall is the lowest and the uncertainty surface runoff prediction is also low, while in October, the month with the highest rainfall, the uncertainty is high.Surface runoff mainly occurred in areas with a slope higher than 25%.This slope class includes the farmlands, which have a surface runoff larger than forest land use.Nevertheless, the surface runoff in the paddy rice lands is low because the slopes in this area usually break down by terracing, therefore, the time concentration is increased.The soil erosion and sedimentation followed the surface runoff behavior [5].Therefore, with an increase in surface runoff, the soil erosion and hence sedimentation is increased.The same results were presented by Phuong et al. [5].The results revealed that surface runoff is highest when the amount of rainfall reaches the maximum amount during the year.

Surface Runoff
The result of the SWAT hydrological model shows that the surface runoff is high in upstream areas of basins.Vegetation cover, soil type and slope gradient [45] are the important factors affecting surface runoff.Land use change and deforestation lead to increasing surface runoff and hence soil erosion.The results showed that the spatial average of surface runoff in the watershed varied from 915 mm/year to 1072 mm/year from 2006 to 2013.The area upstream has the highest amount of surface runoff probably due to its high slope.
The monthly average 95PPU of surface runoff was extracted in order to understand the temporal distribution and uncertainty of surface runoff (Figure 10b).This figure reveals that October has the highest surface runoff, with an average of approximately 368 mm/month.The surface runoff during September to November has the highest uncertainty, probably due to a conceptual model and parameter uncertainties.However, all uncertainties are depicted in the 95PPU band, shown in gray in the figure.As the Figure 10b shows, there is a strong relationship between surface runoff and rainfall.In February, the rainfall is the lowest and the uncertainty surface runoff prediction is also low, while in October, the month with the highest rainfall, the uncertainty is high.Surface runoff mainly occurred in areas with a slope higher than 25%.This slope class includes the farmlands, which have a surface runoff larger than forest land use.Nevertheless, the surface runoff in the paddy rice lands is low because the slopes in this area usually break down by terracing, therefore, the time concentration is increased.The soil erosion and sedimentation followed the surface runoff behavior [5].Therefore, with an increase in surface runoff, the soil erosion and hence sedimentation is increased.The same Our results showed that the ETa product of MODIS can be used for the calibration of the hydrological model in case of data scarcity.However, it is important to note that bias correction should apply to this product before using it in the calibration process.
Based on the results, surface runoff was high, probably due to the geometry of the basin and its soil properties.The surface runoff is high, which leads to rising soil erosion and sedimentation.Our results show that using BMPs leads to decreased surface runoff and hence soil erosion in agricultural lands.Therefore, we suggest implementing terracing/contouring systems in the steep farmlands of the A-Luoi area.

Figure 1 .
Figure 1.Study area showing stream networks, climate and hydrometric stations and the A-luoi district border.The A-Luoi district consists of three watersheds, which are shown in this figure.The Nam Dong watershed is a nearby watershed to A-Luoi with measured discharge data, which used for regionalization approach.

Figure 1 .
Figure 1.Study area showing stream networks, climate and hydrometric stations and the A-luoi district border.The A-Luoi district consists of three watersheds, which are shown in this figure.The Nam Dong watershed is a nearby watershed to A-Luoi with measured discharge data, which used for regionalization approach.

Figure 2 .
Figure 2. Paddy rice schedule plan in A-Luoi district.

Figure 2 .
Figure 2. Paddy rice schedule plan in A-Luoi district.

Figure 4 .
Figure 4. Scatter plot of daily river discharge for calibration (a) from 1 January, 2006 to 31 December, 2008 and (b) the validation period from 1 January, 2009 to 31 December, 2010.

Figure 4 .
Figure 4. Scatter plot of daily river discharge for calibration (a) from 1 January, 2006 to 31 December, 2008 and (b) the validation period from 1 January, 2009 to 31 December, 2010.

Figure 5 .
Figure 5. Calibration (a) and Validation (b) of model for crop yield (rice) based on the average of yield in the selected sub-basins.

Figure 6 .
Figure 6.Correlation between potential evapotranspiration (PET) of (MODIS products and SWAT.

Figure 5 .
Figure 5. Calibration (a) and Validation (b) of model for crop yield (rice) based on the average of yield in the selected sub-basins.

Figure 5 .
Figure 5. Calibration (a) and Validation (b) of model for crop yield (rice) based on the average of yield in the selected sub-basins.

Figure 6 .
Figure 6.Correlation between potential evapotranspiration (PET) of (MODIS products and SWAT.

Figure 6 .
Figure 6.Correlation between potential evapotranspiration (PET) of (MODIS products and SWAT.

Figure 7 .
Figure 7. Comparing PET extracted from MODIS-MOD16 with one measured using the Penman-Monthieh method.For bias correction, MODIS data were multiplied by 0.45 for the period of October-March and 0.70 for April-September.
Figure 7.Comparing PET extracted from MODIS-MOD16 with one measured using the Penman-Monthieh method.For bias correction, MODIS data were multiplied by 0.45 for the period of October-March and 0.70 for April-September.

Figure 8 .
Figure 8. Results of calibration (a) and validation (b) of ETa using MODIS data.

Figure 8 .
Figure 8. Results of calibration (a) and validation (b) of ETa using MODIS data.

Figure 10 .
Figure 10.Average monthly (a) evapotranspiration and (b) surface runoff against rainfall.The gray bands show 95% prediction uncertainty (95PPU) and the red stars show the median of iteration (during 2006-2013).

Figure 10 .
Figure 10.Average monthly (a) evapotranspiration and (b) surface runoff against rainfall.The gray bands show 95% prediction uncertainty (95PPU) and the red stars show the median of iteration (during 2006-2013).

Table 1 .
Properties of the donor and target catchments.
* Dominant landuse and soil type are mentioned; ** L/W: Length/Width of catchment.

Table 2 .
Soil and Water Assessment Tools (SWAT) parameters and their initial and final range for river discharge, evapotranpiration (ETa) and crop yield calibration.

Table 2 .
Soil and Water Assessment Tools (SWAT) parameters and their initial and final range for river discharge, evapotranpiration (ETa) and crop yield calibration.r parameter value is multiplied by 1 + given value, v parameter value is replaced by a value from the given range. *