Modeling Surface Runoff and Evapotranspiration using SWAT and BEACH for a Tropical Watershed in North Vietnam, Compared to MODIS Products

Accurate estimation of surface runoff (Q) and evapotranspiration (ET) is a challenging task but an important research topic because both Q and ET play vital roles in the study of the hydrological cycle, of climate change, water resources, flood management and so on. In this paper we will present the modeling method to estimate the daily Q and ET for a medium-sized watershed in the tropical region of the North of Vietnam using the Soil and Water Assessment Tool (SWAT) and Bridging Event and Continuous Hydrological (BEACH) models. The models were calibrated and validated for the river discharge for SWAT and evaporation (E) for BEACH in a 12 year period from 2001 to 2012. The simulated ETs by the models were compared with the satellite-based ET of MODIS products. Our simulation results show that the SWAT and BEACH models are capable of satisfactorily reproducing (with the NSE > 0.62 and R 2 > 0.78) the stream-gauged river discharge and the observed E, respectively. Daily ET varied from 0.3 to 14 mm day -1 and was highest from May to August and lowest from December to March. Although the monthly and yearly MODIS ETs were slightly higher than those of SWAT and BEACH, a strong relationship between them was found with a standard deviation ranging from three to 40 mm. A light decrease of ET values in the 12 years can be seen in the result analyses; however a longer simulation time might be needed to ensure this trend.

in streams to the sea, all of which is dependent on the amount of rainfall, rainfall intensity and infiltration capacity [14]. The evapotranspiration is considered as the total water loss to the atmosphere [15] by evaporation both from the vegetative and non-vegetative surfaces and transpiration from plants [11,16]. Since the in situ measurements of surface runoff and ET for large areas are time consuming, costly and extremely difficult [3], indirect estimations of Q and ET using the modeling approach and satellite-based products are necessary [17], particularly for sparse data-available areas such as in northern Vietnam.
Although many studies have been conducted using the modeling method at watershed scale for estimating Q, to name some of them [18][19][20][21] employing the SWAT, KINEROS2 and the HBV and VIC models, little attention has been paid to tropical regions using the SWAT model such as [22,23] and for regions of Vietnam [24]. Similarly with studies on surface runoff, there is a large amount of literature available on research into ET including direct or in situ measurement [25,26] and indirect methods (referring to modeling [3,16] and using satellite technology [27,28]), from watershed/catchment scales [29,30] to the global scale [31][32][33] and others. Surprisingly, despite many scientific works having been conducted on ET elsewhere in the world as mentioned above, few or even no previous studies on that topic were found for the case study in Yen Bai province or in North Vietnam.
As mentioned at the very beginning, the importance of surface runoff and ET information and the lack of Q and ET literature on the study region have motivated the study objective of modeling and estimating that information. In addition, data scarcity has been reduced thanks to more ground-based hydrologic and meteorological stations being established in the region (for model calibration and validation) and in addition satellite data (MODIS, Landsat etc.) have become more available to make up for the lack of temporal and spatial data in ungauged regions [34] and be used for model validation [35] as well. The interesting point is that both large-scale in situ Q and ET measurements are difficult, expensive and time consuming [13] but the modelling approach could be appropriate to deal with this challenge. Finally, the study results might be helpful for research into other fields such as Q information for flood, ground/surface water management and ET for cropping irrigation, drought detection and management etc., particularly rapid change in climate, water and the increasing population have all become a great concern both for the environment and society [31].
In summary five main outcomes of expertise could be identified to cope with this study topic: (i) calibrating and validating daily river discharge and evaporation for the SWAT and BEACH model respectively, (ii) estimating daily ET by the two models, (iii) comparing SWAT, BEACH discharge and with observed data, (iv) comparing monthly ETs and analyzing their trend in the 2001-2012 period and, (v) mapping SWAT and MODIS annual accumulative ET. Our results elucidate the abilities of the two models to produce surface runoff and evaporation closely to field measured data, the strong relationships between SWAT, BEACH and MODIS ET and a slightly decreasing trend of ET in the simulating time.

Study Site
The study watershed is Nam Kim located in western Yen Bai province ( Figure 1) has its center coordinate as 104°07'51.3"E and 21°49'10.7"N. The mean elevation is 1571 meters above sea level and the area of the watershed is about 268 km 2 . The main land use land covers (LULC) are forest, agricultural (rice) and shrub lands with contributions of 40, 33 and 23% of the watershed area respectively. The yearly mean air temperature varies from 20 to 24 o C [36] and annual mean precipitation is 1695 mm (2001-2012 statistic from the rain gauge presented on the map). The dominant soils are Humic Ferralsols, Humic Acrisols and Humic Alisols and the cropping system is mainly rice and cassava harvested twice a year in May and October. There are three gauges established for measurements of rainfall, discharge and other meteorological data.

The Soil and Water Assessment Tool (SWAT)
SWAT is a physically based, continuous model and is computationally efficient, enabling users to study long-term impacts and is a long-term yield model [37]. The model simulates surface runoff volumes and peak runoff rates for each hydrologic response unit (HRU) using a modification of the SCS curve number method [38] or the Green & Ampt infiltration method [39]. The evaporation from soils and plants are computed separately by the method described in [40]. Wherein, the potential evapotranspiration (PET) is the ET rate of a large area completely and uniformly covered with growing vegetation and an unlimited supply of soil water. In addition, the micro-climatic processes such as advection or heat-storage should not have effects on PET. The model uses Hargreaves [41], Priestley-Taylor [42] and Penman-Monteith [43] approaches for the estimation of PET.

The Bridging Event and Continuous Hydrological (BEACH)
The BEACH is a simple two-layer soil water balance model developed by Sheikh et al. [44] using freely available GIS and programming language, PCRaster. It is a spatially distributed daily basic hydrological model and has the ability to predict antecedent soil moisture, actual E and ET as its primary objective. The model also estimates the infiltration and runoff using either the Bucket or CN approach (for more details see Sheikh et al., 2009). Although the BEACH model offers two options of calculating runoff, the study results of [44] showed that the CN approach performed better than the Bucket. That is reason for selecting the CN approach for simulating Q by the BEACH in this study.
In this model, the ET is separated into evaporation and transpiration and calculated based on the Penman-Monteith (FAO56) method [45]. Based on an assumption of a hypothetical grass cover with an assumed height of 0.12m, a surface resistance of 70sm -1 and an albedo of 0.23, the FAO56 computes a reference evapotranspiration rate (ET 0 ). The PET of the other surface is calculated by multiplying the ET 0 with a crop-specific coefficient (K c ).
The BEACH was developed for just one year (2004) hydrologic assessments, but for this study objective of long-term hydrological modeling (12 years) we recoded the code and made the model adapt to the requirement with an assumption that the crop harvesting seasons have not changed in the computational time.

MODIS Evapotranspiration
We used the MODIS ET of the improved algorithm [32] from the old algorithm [31] calculated as the sum of the E from moist soil, the transpiration from vegetation during the daytime, intercepted canopy precipitation and ET in the nighttime. Some other improvements are the vegetation cover fraction (F c ), the soil heat flux (G), the wet surface fraction (F wet ), evaporation from the wet canopy surface (E wet_C ), plant transpiration (E trans ) etc. (more details see Mu et al. [32]).
The total daily ET is calculated as in following relation: (1) Where is the actual daily soil evaporation.
The monthly ET is summarized from daily ETs.

SWAT Input
SWAT requires four groups of input data including topography, LULC, soils and climatic data.
 The topographic dataset commonly used is a Digital Elevation Model (DEM) which is employed to generate drainage network, flow directions, flow accumulation etc. for the topographical model's parameterization. We used a Triangular Irregular Network (TIN) extracted from Yen Bai geodatabase (provided by the Vietnam Resources and Environment Corporation, Hanoi Vietnam) for generating the 30x30 m DEM for the study area ( Figure 1).
 The Land use land cover (LULC) map ( Figure 2) of nine classes (including cloud) was produced from Landsat scenes acquired in 2009 with 30x30 m spatial resolution, seven bands of reflections. The map was classified using the supervised method in ENVI 4.7 software. The ground truth data was derived from the Yen Bai geodatabase. The accuracy of the classification was examined by calculating the user, producer accuracy and Kappa statistics for categorized classes ( Table 1). The data shown within the table indicates the precision of the classification is at a good level with an overall accuracy of 78.6 % and the Kappa coefficient of 0.73.  Climatic data required by the model includes daily rainfall, minimum and maximum temperature, wind speed, relative humidity, solar radiation and dew point temperature. This data can be read from measured data set or generated by a weather generator algorithm. We collected this data from the gauges (

BEACH Inputs
Similar to SWAT, BEACH requires four main groups of input data: (i) topography, (ii) lands use, (iii) soil physical properties and (iv) climatic data. For the land use information, we need extra crop-specific information that could be obtained from literature of FAO paper no. 56 [45]. The overview of input for BEACH was indicated in the Table 2.

Model Calibration, Validation and Simulation
 SWAT: Before model calibration, the watershed was modelled by dividing it into 20 HRUs ( Figure  1) and parameterized with the inputs of topography, LULC and soils. The parameters were summarized in Table 3. Climatic parameterization was separately implemented and written to climatic datasets.  ) and Nash-Sutcliffe efficiency (NSE) [46]. We used a simple averaging method (to average ET from all crops appearing in the HRU 3) to calculate the mean BEACH ET for HRU 3 and then compared it with SWAT and MODIS ET (used zonal statistics) for this HRU.

MODIS Evapotranspiration
Monthly MODIS ET datasets with 1-km spatial resolution generated based on the MOD16 algorithm in [32] from 2001 to 2012 was used to extract actual ET. These ET was verified by comparing with measured data of both globally and locally at 46 eddy flux towers presented in the study of Mu et al. [32]. We developed a simple model in the ModelBuilder of ArcMap-ArcInfo to extract the monthly and yearly ET using ArcGIS functions as in the Figure 3.

SWAT Calibration and Validation for Surface Runoff
Overestimates  Although there were still some over/underestimates of river discharge at some peaks, the result of model calibration and validation revealed the important capability of SWAT to produce the Q matching closely to measured data in a medium-size tropical watershed (SWAT was developed for arid/semi-arid river catchments). In addition, this point was supported by some similar previous studies of [19,22,23] even for a case study in Vietnam [24].

BEACH Calibration and Validation for Evaporation
The BEACH started calibrating for a 2001-2004 period for four types of crops (mentioned in Table 4) and the before-calibrated simulated E (averaged from rice and shrub E) was much lower compared to observed E ( Besides the climatic aspects such as humidity, wind, radiation and precipitation, the cropping irrigation, crop growing [30] and land use changes [4] have important effects on E and ET. For this case study, some differences between simulated and observed E might come from changes in crop and land use because we kept the cropping sowing dates and land use unchanged over the simulation period. To improve the model performance, an updated crop type and time might be needed, however for long-term observation it is very time consuming.

Comparing SWAT and BEACH Daily ET (of HRU3 from 2001 to 2012)
The time series daily actual evapotranspiration values calculated by SWAT (for HRUs) and BEACH (averaged from the four crop's ET) and measured rainfall were presented and compared to each other by the line graph ( Figure 6). The graph indicated not only the variations of ET (0.3 to 14 mm day -1 ) and rainfall but also the yearly routine of ET. The ET values appeared very high before rainy seasons (May to August) and lowest in winter (December to March). The changing in ET is basically related to stand density (number of trees per hectare) and also incoming solar radiation (referring to season) [49]. A better understanding of the ET routine and trends might help to predict the ET in the future but a longer time of simulation is needed (instead of 12 years). There were two periods May to September 2003 and 2004 with large differences in ETs estimated by SWAT and BEACH. Unfortunately, without the in situ measured ETs for comparison we did not know whether the SWAT ET was more accurate than the BEACH ET or the reverse. However, based on the good relationship between BEACH E and observed E (discussed in section 4b) the BEACH ET (lower value) would be more reliable in these computational times.

Figure 6: Line Graph of Time Series SWAT and BEACH ET (Calculated for HRU 3) and Measured Rainfall
We analyzed the relationship between the two estimates of ET represented in the Figure 7 with the moderate correlation (R 2 of 0.73). However, both models predicted ETs closely together when the rates were low (< 4 mm). These rates consisted of more than 90% of the estimates. The interesting correlation between measured daily rainfall and SWAT ET was also indicated on the figure. Basically, the wet days (with precipitation > 0 mm day -1 ) had ET rates from one to four mm day -1 . On the other hand, all the ET values are higher than six mm day -1 generated in dry day (no rainfall). The ET of 0.3 to one mm day -1 and rainfall less than 10 mm day -1 were calculated for the winter time. This relationship could provide helpful information for cropping irrigation [35], the drought season [50], particularly when we could forecast the ET values [51], and for improved understanding of the interactions between land surfaces and the atmosphere [31,52]. . The SWAT appeared to predict the discharge closer to observed information than the BEACH with NSE of 0.67 and 0.61 for the models respectively. Generally, both models slightly underestimated the discharge. However, close relationships between simulated and field measured discharge can be seen in the figure in both cases with R 2 of 0.81 for the SWAT and 0.83 for the BEACH.
Although the two models use the same CN method [53] to simulate the discharge value, many reasons could be explained for the differences of Q values shown in the Figure 8 such as the difference in model structure, in inputs etc. Many previous studies found that the modeled runoff is sensitive to CN values [19,54,55] we altered the CN value of each land use of each class (of crops for the case of he BEACH) in oder to calculate to discharge close to the observed data.

Monthly SWAT and BEACH vs. MODIS ET of HRU 3 from 2001-2012
The variations of monthly ET calculated from MODIS scenes, from SWAT and BEACH models and their standard deviations from 2001 to 2012 are presented in the bar graph of Figure 9. Although the ET from all estimates was highest from May to October and lowest from December to February, the MODIS ET varied gradually over the twelve years. Both the SWAT and BEACH estimated the ET in 2003 and 2004 much lower than the MODIS ET. The MODIS ET was estimated higher than the BEACH ET and the SWAT ET was lowest in general. The standard deviations of the three ET sources calculated for each month fluctuated from three to 40 mm and the average value was 17. 5 mm. The MODIS ET was estimated using the improved ET algorithm at global scale (global parameters) [32] and compared to the watershed scale in this study we assert that the agreement between the three estimated ETs was at a good level.
The monthly ET varied over different seasons and the most driving factors on ET found in the study of [11] were relative humidity and wind speed at two meters high. Both the SWAT and BEACH models thoroughly take into account these two factors but the algorithm for MODIS ET does not include wind speed and topographic influences (elevation, slope and aspect) [33].

ET Trend Analyses
The slightly downward trend of ET in the 2001-2012 periods has been shown in the Figure 10 for all cases of estimations illustrated by the negative slope line of the equations. However, in the first four years the ET decreased gradually but in the rest of the years the rates were more stable. While all data sources represented the decrease of monthly ET, the 12 years of simulation might be too short to come to a conclusion that the annual ET in the study region is currently declining and will in the future. Therefore, longer term assessment could be needed (40 years or more) and the reasons for this trend were also out of this paper's scope. . In general, the MODIS maps had higher ET rates than SWAT maps. The yearly ET calculated for the HRUs also varied over time in all the maps. Some similarities in the maps can be seen in the HRUs such as HRU 7, 8, 10, 13, 14 and 20 (except the year 2004) with lower ET and 1, 5, 11 and 19 with higher ET. Differences between SWAT and MODIS distributed ET appeared in the figure as well, for example the HRU 12, 16 and 14 with lower ET in MODIS maps and higher ET in the SWAT ones. These dissimilarities were thought to be normal because the results were estimated from different approaches and scales. However, the agreements were dominant.
The changes of ET patterns might be related to changes in land use land cover for example the decrease in canopy interception causes a decrease in ET and percolation and increase in runoff [4,56]. In addition, as discussed in the section 4f that SWAT used distributed soil and topographic inputs, therefore the SWAT ET might be more precise than MODIS due to the fact that global satellite products generally contain noises [34].

Summary and Conclusions
Daily discharge and evaporation were calculated and validated by SWAT for the 268 km 2 watershed in the tropical region in North Vietnam with acceptable agreement between simulated and observed data verified by the NSE and R 2 values. This might reveal a possible application of the two hydrological models for tropical regions. The good correlations between discharges produced by the models were also represented in the study results as well.
Although the daily estimated ETs by the SWAT and BEACH models were not validated, they matched well with each other and the monthly ETs were compared with the published MODIS product. Despite the SWAT and BEACH ETs being slightly lower than the MODIS ET, basically a close correlation between them can be seen in the study results (analyzed using standard deviations) and also all the monthly ETs showed the slight downward trend in the simulation time (2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012).
The results of zonal statistics applied for the yearly MODIS and SWAT ETs were mapped, providing interesting information of temporal and distributed ET patterns in the watershed. Both differences and similarities could be found in the map but the correspondence between them was dominant. We conclude that the MODIS ET was very helpful for verifying the smaller scale of ET estimation by the two models.