Calibration of Spatially Distributed Hydrological Processes and Model Parameters in SWAT Using Remote Sensing Data and an Auto-Calibration Procedure : A Case Study in a Vietnamese River Basin

In this paper, evapotranspiration (ET) and leaf area index (LAI) were used to calibrate the SWAT model, whereas remotely sensed precipitation and other climatic parameters were used as forcing data for the 6300 km2 Day Basin, a tributary of the Red River in Vietnam. The efficacy of the Sequential Uncertainty Fitting (SUFI-2) parameter sensitivity and optimization model was tested with area specific remote sensing input parameters for every Hydrological Response Units (HRU), rather than with measurements of river flow representing a large set of HRUs, i.e., a bulk calibration. Simulated monthly ET correlations with remote sensing estimates showed an R2 = 0.71, Nash–Sutcliffe Efficiency NSE = 0.65, and Kling Gupta Efficiency KGE = 0.80 while monthly LAI showed correlations of R2 = 0.59, NSE = 0.57 and KGE = 0.83 over a five-year validation period. Accumulated modelled ET over the 5-year calibration period amounted to 5713 mm compared to 6015 mm of remotely sensed ET, yielding a difference of 302 mm (5.3%). The monthly flow at two flow measurement stations were adequately estimated (R2 = 0.78 and 0.55, NSE = 0.71 and 0.63, KGE = 0.59 and 0.75 for Phu Ly and Ninh Binh, respectively). This outcome demonstrates the capability of SWAT model to obtain spatial and accurate simulation of eco-hydrological processes, also when rivers are ungauged and the water withdrawal system is complex.


Introduction
Managing river basins and environmental systems in a sustainable manner is receiving growing attention from national water resources institutes, the United Nations, non-governmental-reporting ecosystem services to governments and the United Nations. The novelty of this study is the application of SUFI-2 in the optimization of 15 input parameters of the soil vegetation processes using observations of actual water management processes, such as irrigation and conservation of water in wetlands in a local context. Such level of detail and reflection of real world interferences of mankind on the natural hydrological cycle can never be achieved from flow measurements, and opens better opportunities for simulation of local eco-hydrological processes that occur locally in ungauged basins, which can never be interpreted from bulk flow measurements, if there are any.

Study Area
The Day Basin is located between 19 • 55 to 21 • 10 N and 105 • 20 to 106 • 25 E. The Day Basin is a sub-basin of the transboundary Red River basin (See Figure 1). The total area of the basin is nearly 6300 km 2 . The highest elevation is 1256 m in the western part of the basin. The Day Basin comprises several river tributaries, among which the largest is the Day River, with a total length of approximately 250 km. The Day Basin has a high biodiversity, with abundant flora and fauna in the forested hills, freshwater aquatics, and wetland. The land use is also diversified, although agricultural land use is dominant (64%).
Water 2018, 10, x 3 of 21 opportunities for simulation of local eco-hydrological processes that occur locally in ungauged basins, which can never be interpreted from bulk flow measurements, if there are any.

Study Area
The Day Basin is located between 19°55′ to 21°10′ N and 105°20′ to 106°25′ E. The Day Basin is a sub-basin of the transboundary Red River basin (See Figure 1). The total area of the basin is nearly 6300 km 2 . The highest elevation is 1256 m in the western part of the basin. The Day Basin comprises several river tributaries, among which the largest is the Day River, with a total length of approximately 250 km. The Day Basin has a high biodiversity, with abundant flora and fauna in the forested hills, freshwater aquatics, and wetland. The land use is also diversified, although agricultural land use is dominant (64%).
The Day Basin encompasses the capital city of Hanoi (population in 2015: 7.5 million inhabitants) in the northeast and several major economic centers located downstream, such as Nam Dinh (population: 1.8 million) and Ninh Binh (population: 0.9 million). Both the Red River and Day Basin have been exposed to various hydrological research activities before, e.g., [8,[50][51][52] on hydrological and water quality modelling. Using the rating curve suggested in [50], discharge at two locations (i.e., Ninh Binh and Phu Ly) was reconstructed from the year 2000 up to 2013. The annual total precipitation is around 1700 mm per year, and reference evapotranspiration (ETo) is approximately 1100 mm per year. The climate in the Day Basin has a monsoonal character. The wet season lasts from May to September, and dry season from October to April. Precipitation can reach up to 450 mm per month in some parts of the basin, and as low as a few mm during January and February ( Figure 2). Precipitation is measured at nine stations across the basin, and is available up to 2013. These measurements will be used to validate the open access precipitation product based on satellite measurements.
The water withdrawals for irrigation in the Day Basin are rather difficult to assess because various pumping stations lift water from the Red River, and also many inlets divert water from the river gravitationally. This diffusive and unmetered inter-basin water withdrawal complicates the computation of the irrigation hydrology and the water accounts related to that. The irrigation supplies in SWAT will therefore be adjusted to reproduce an ET value that matches with ET estimates from satellites. The Day Basin encompasses the capital city of Hanoi (population in 2015: 7.5 million inhabitants) in the northeast and several major economic centers located downstream, such as Nam Dinh (population: 1.8 million) and Ninh Binh (population: 0.9 million). Both the Red River and Day Basin have been exposed to various hydrological research activities before, e.g., [8,[50][51][52] on hydrological and water quality modelling. Using the rating curve suggested in [50], discharge at two locations (i.e., Ninh Binh and Phu Ly) was reconstructed from the year 2000 up to 2013.
The annual total precipitation is around 1700 mm per year, and reference evapotranspiration (ET o ) is approximately 1100 mm per year. The climate in the Day Basin has a monsoonal character. The wet season lasts from May to September, and dry season from October to April. Precipitation can reach up to 450 mm per month in some parts of the basin, and as low as a few mm during January and February ( Figure 2). Precipitation is measured at nine stations across the basin, and is available up to 2013. These measurements will be used to validate the open access precipitation product based on satellite measurements.
The water withdrawals for irrigation in the Day Basin are rather difficult to assess because various pumping stations lift water from the Red River, and also many inlets divert water from the river gravitationally. This diffusive and unmetered inter-basin water withdrawal complicates the computation of the irrigation hydrology and the water accounts related to that. The irrigation supplies in SWAT will therefore be adjusted to reproduce an ET value that matches with ET estimates from satellites.

Soil and Water Assessment Tool (SWAT)
The Soil and Water Assessment Tool (SWAT) [10,53] has been set up for the Day Basin to compute flows, fluxes, and stocks. A total amount of 119 sub-basins and 7909 HRU has been included for ensuring sufficient detail. HRU is a modeling unit that exists of a unique combination of land slope, land use, and soil type [54]. SWAT simulates eco-hydrological processes, i.e., surface runoff, groundwater recharge, baseflow, water stocks, erosion, plant production, and water quality. The production of food, feed, and timber, and the sequestration of carbon can be inferred from the biomass production [54]. SWAT estimates the fate and transport of nutrients, sediment, pesticides, and bacteria in both land and water phases [55]. This mathematical framework provides a great basis for the determination of various ecosystem services and SDG indicators. The soil water balance is conceptualized in SWAT using Equation (1), as described in [53]: in which t is the time (days), SWt is the final soil water content at day t (mm H2O); SWo is the initial soil water content, P is the amount of precipitation, Qsurf is the amount of surface runoff, Ea is the amount of actual evapotranspiration, wseep is the amount of percolation entering the vadose zone from the soil profile, and Qgw is the volume of streamflow originating from groundwater, all measured in mm H2O on day i. The reference evapotranspiration (ET0) is computed with Global Land Data Assimilation System (GLDAS) meteorological input data. SWAT does not allow reading layers of ET0 directly, and therefore, meteorological records need to be prescribed. The chosen method for reference ET in SWAT is the Penman-Monteith method, as described in [53]. Compared with other equations in SWAT, i.e., Hargreaves [56] and Priestley-Taylor [57], the Penman-Monteith method fully takes advantage of remote sensing and global data source, such as GLDAS.

Model Calibration Using SUFI-2
The calibration of a semi-distributed and physical based model, such as SWAT, requires various plant and soil model parameters to be optimized, to ensure a rigorous representation of a basin's processes, e.g., streamflow, ET, ecological change, etc. The calibration task can become difficult and almost infeasible in large-scale applications [54]. A number of auto-calibration and uncertainty analysis tools for SWAT were developed to support solving this problem, and are currently available to assist the optimization process. This study is based on the SUFI-2 model that is part of the SWAT-CUP supporting

Soil and Water Assessment Tool (SWAT)
The Soil and Water Assessment Tool (SWAT) [10,53] has been set up for the Day Basin to compute flows, fluxes, and stocks. A total amount of 119 sub-basins and 7909 HRU has been included for ensuring sufficient detail. HRU is a modeling unit that exists of a unique combination of land slope, land use, and soil type [54]. SWAT simulates eco-hydrological processes, i.e., surface runoff, groundwater recharge, baseflow, water stocks, erosion, plant production, and water quality. The production of food, feed, and timber, and the sequestration of carbon can be inferred from the biomass production [54]. SWAT estimates the fate and transport of nutrients, sediment, pesticides, and bacteria in both land and water phases [55]. This mathematical framework provides a great basis for the determination of various ecosystem services and SDG indicators. The soil water balance is conceptualized in SWAT using Equation (1), as described in [53]: (1) in which t is the time (days), SW t is the final soil water content at day t (mm H 2 O); SW o is the initial soil water content, P is the amount of precipitation, Q surf is the amount of surface runoff, E a is the amount of actual evapotranspiration, w seep is the amount of percolation entering the vadose zone from the soil profile, and Q gw is the volume of streamflow originating from groundwater, all measured in mm H 2 O on day i. The reference evapotranspiration (ET 0 ) is computed with Global Land Data Assimilation System (GLDAS) meteorological input data. SWAT does not allow reading layers of ET 0 directly, and therefore, meteorological records need to be prescribed. The chosen method for reference ET in SWAT is the Penman-Monteith method, as described in [53]. Compared with other equations in SWAT, i.e., Hargreaves [56] and Priestley-Taylor [57], the Penman-Monteith method fully takes advantage of remote sensing and global data source, such as GLDAS.

Model Calibration Using SUFI-2
The calibration of a semi-distributed and physical based model, such as SWAT, requires various plant and soil model parameters to be optimized, to ensure a rigorous representation of a basin's processes, e.g., streamflow, ET, ecological change, etc. The calibration task can become difficult and almost infeasible in large-scale applications [54]. A number of auto-calibration and uncertainty analysis tools for SWAT were developed to support solving this problem, and are currently available to assist the optimization process. This study is based on the SUFI-2 model that is part of the SWAT-CUP supporting software package [58]. SWAT-Calibration and Uncertainty Program (SWAT-CUP) [24,58] is an auto-calibration and uncertainty analysis module program based on the SWAT engine that can deal with a range of input parameters. The optimization algorithm of SWAT-CUP allows model parameters to be predefined and optimized throughout the auto calibration process or manually adjusted iteratively between calibration batches. The employment of SUFI-2 is suitable for both new and advanced users of hydrological models, even though a good understanding of hydrologic processes and of parameter sensitivity is recommended in general terms [54]. In order to assess the efficacy of using remote sensing data and to compare with traditional discharge station-based for calibration, SWAT-CUP was run with two settings: (a) entirely based on remote sensing data (ET and LAI) and (b) entirely based on river discharge measurements. Among various evaluation coefficients allowed in SUFI-2, Nash-Sutcliffe (NSE) and Kling Gupta Efficiency (KGE) were chosen. For this particular study, a total number of 15 plant and soil parameters were pre-selected according to their sensitivity to the evolution of ET and LAI. The selection of these parameters was based on detailed reviews and analyses on SWAT parameters carried out before by various authors (e.g., [16,42,43,54,59]).

Physiographical Maps
The Digital Elevation Model (DEM) was downloaded from the Shuttle Radar Topography Mission (SRTM) with a resolution of 1 arc-second or 30 m (Figure 1b). The DEM is used to calculate slope, slope lengths, and to extract the stream network, solar angles, and air temperature corrections. The land use map is downloaded from Globcover [60]. Globcover is developed by the European Space Agency (ESA) and University of Louvain, with a spatial resolution of 300 m × 300 m. The satellite input data used for the classification was the MEdium Resolution Imaging Spectrometer (MERIS) sensor on the Environmental Satellite (ENVISAT) during 2009. While this dataset exists for several years, it captures the time span of the SWAT analysis very well. The major land cover in the Day Basin is agricultural land (64%), followed by forests (24%), and mixed mosaic (12%). Three thousand hectares (76%) of agricultural land is irrigated. The soil map used in this study originates from the International Soil Reference Information Centre (ISRIC) [61] and Food and Agricultural Organization (FAO) Digital Soil Map of the World [62]. The SoilGrids database (ISRIC, Wageningen, The Netherlands) has a spatial resolution of 1 km × 1 km, and is produced during 2014. The physical properties included in the dataset are (i) soil organic carbon (g/kg), (ii) pH index (H 2 O solution) (%), (iii) sand, silt and clay content (kg/kg), (iv) coarse fragments (volumetric) (%), (v) bulk density (kg/m 3 ), (vi) cation-exchange capacity (fine earth fraction) (cmol+/kg), and (vii) depth to bedrock (cm). A new soil map was created in this study by combining the two ISRIC and FAO soil maps with the aims to both (i) increase the spatial representation, and (ii) maintain the soil classification and soil properties from the FAO database. This task was accomplished by using standard unsupervised classification procedures (see Figure 2). Based on the distribution of land slope, soil type, and land use classes, the basin is divided into 119 sub-basins and 7909 HRUs.

Precipitation
Satellite precipitation data offers an attractive alternative to supplement in situ precipitation measurements in hydrological modelling, particularly in poorly gauged basins [29,63] Precipitation Estimation from Remotely Sensed Imagery Using Artificial Neural Networks (PERSIANN, The University of Arizona, Tucson, AZ, USA) at a resolution of 0.25 • × 0.25 • for a subtropical watershed in China, concluded that TRMM 3B42 had the best performances, and was deemed to be reliable for hydrological applications, while PERSIANN had the worst performance [64]. In a similar study for Australia and Southeast Asia, TRMM and CMORPH outperformed, and an ensemble precipitation product was suggested as a reduction of system-specific and random errors [65]. TRMM data, in situ measurements and other atmospheric and climatology models were assimilated in [66] to create an ensemble precipitation product Climate Hazards Group InfraRed Precipitation with Station data (CHIRPS, University of California Santa Barbara -Climate Hazards Group, Santa Barbara, CA, USA) with a superior resolution at 0.05 • × 0.05 • . Precipitation from CHIRPS performed well statistically for flood and drought monitoring, particularly for meteorological complex regions [67].
The current study combines TRMM7.0 and CHIRPS2.0 rainfall products. The absolute precipitation data are taken from TRMM, and the spatial patterns from CHIRPS. The refined TRMM dataset with a resolution of 0.05 • × 0.05 • , so obtained, has been used as input data for SWAT. The combined precipitation product, so obtained, showed a good performance when compared to rain gauge measurements ( Figure 3). Percent bias (PBIAS), NSE, Mean Absolute Error (MAE) denoted percent bias, Nash-Sutcliffe efficiency, and mean absolute error, respectively. The newly created precipitation product significantly improved the performance of the two original datasets in term of bias correction (PBIAS = −0.6) when averaging the errors from TRMM and CHIRPS. Nash-Sutcliffe efficiency slightly improved when comparing the ensemble precipitation to the CHIRPS (0.75 compared to 0.74), even though the MAE was marginally larger (45.98 compared to 44.31).
Water 2018, 10, x 6 of 21 precipitation product was suggested as a reduction of system-specific and random errors [65]. TRMM data, in situ measurements and other atmospheric and climatology models were assimilated in [66] to create an ensemble precipitation product Climate Hazards Group InfraRed Precipitation with Station data (CHIRPS, University of California Santa Barbara -Climate Hazards Group, Santa Barbara, California, U.S.) with a superior resolution at 0.05° × 0.05°. Precipitation from CHIRPS performed well statistically for flood and drought monitoring, particularly for meteorological complex regions [67]. The current study combines TRMM7.0 and CHIRPS2.0 rainfall products. The absolute precipitation data are taken from TRMM, and the spatial patterns from CHIRPS. The refined TRMM dataset with a resolution of 0.05° × 0.05°, so obtained, has been used as input data for SWAT. The combined precipitation product, so obtained, showed a good performance when compared to rain gauge measurements ( Figure 3). Percent bias (PBIAS), NSE, Mean Absolute Error (MAE) denoted percent bias, Nash-Sutcliffe efficiency, and mean absolute error, respectively. The newly created precipitation product significantly improved the performance of the two original datasets in term of bias correction (PBIAS = −0.6) when averaging the errors from TRMM and CHIRPS. Nash-Sutcliffe efficiency slightly improved when comparing the ensemble precipitation to the CHIRPS (0.75 compared to 0.74), even though the MAE was marginally larger (45.98 compared to 44.31).

Meteorology
The meteorological dataset includes daily estimates of solar radiation, wind speed, air temperature (maximum, minimum) and relative humidity. The dataset was derived from GLDAS (NASA Goddard Earth Sciences Data and Information Services Center (GES DISC), Greenbelt, Maryland, U.S.). GLDAS simulates meteorological data with a numerical weather prediction model having a cell size of 0.25 degrees. The NOAH Land Surface Model (NOAH, National Center for Atmospheric Research (NCAR), Boulder, Colorado, U.S.), coupled to an atmospheric boundary layer model, assimilates satellite and in situ measurements to produce various land surface states and fluxes [68]. GLDAS meteorological data was downloaded as 3-hour intervals, and accumulated into 1-day time step as required by SWAT model to calculate evapotranspiration, plant growth etc.

Actual Evapotranspiration
The most common global scale ET dataset, developed with energy balance models using remote sensing data as input, a Penman-Monteith type of equation, is MOD16 (University of Montana, Missoula, Montana, U.S.) [69]. Several years of data can be downloaded if one is registered. MOD16 is based on a simplified stomatal conductance model governed by LAI, vapor pressure deficit, and air temperature. Soil evaporation is limited by a complementary relationship hypothesis which defines land-atmospheric interactions from vapor pressure deficit and relative humidity [69]. Various ET

Meteorology
The meteorological dataset includes daily estimates of solar radiation, wind speed, air temperature (maximum, minimum) and relative humidity. The dataset was derived from GLDAS (NASA Goddard Earth Sciences Data and Information Services Center (GES DISC), Greenbelt, MD, USA). GLDAS simulates meteorological data with a numerical weather prediction model having a cell size of 0.25 degrees. The NOAH Land Surface Model (NOAH, National Center for Atmospheric Research (NCAR), Boulder, CO, USA), coupled to an atmospheric boundary layer model, assimilates satellite and in situ measurements to produce various land surface states and fluxes [68]. GLDAS meteorological data was downloaded as 3-hour intervals, and accumulated into 1-day time step as required by SWAT model to calculate evapotranspiration, plant growth etc.

Actual Evapotranspiration
The most common global scale ET dataset, developed with energy balance models using remote sensing data as input, a Penman-Monteith type of equation, is MOD16 (University of Montana, Missoula, MT, USA) [69]. Several years of data can be downloaded if one is registered. MOD16 is based on a simplified stomatal conductance model governed by LAI, vapor pressure deficit, and air temperature. Soil evaporation is limited by a complementary relationship hypothesis which defines land-atmospheric interactions from vapor pressure deficit and relative humidity [69]. Various ET comparison studies using MOD16 data have been undertaken, see for instance [32,39,70]. MOD16 ET was validated [71] using flux towers in South Africa, and found that ET was systematically underestimated by 7.5 to 26.3 mm per month, something confirmed by [72,73]. On the contrary, studies conducted in Asia showed that MOD16 consistently overestimates ET for forested land cover [38].
Another example of a global energy balance model is SEBS (University of Twente, Enschede, The Netherlands) [74,75] that has a quasi-open accessibility to acquire the data. SEBS applies an analytical solution of surface roughness for heat transfer, and it was used to create a global scale dataset that is quasi open access. It limits the sensible heat flux estimates with upper and lower boundaries of the surface resistance to evaporation. The upper boundary is determined by latent heat flux equal to zero; the lower limit by potential evapotranspiration with a minimum bulk surface resistance. CMRSET (Commonwealth Scientific and Industrial Research Organisation (CSIRO), Canberra, Australia) [76] calculates actual ET from the Priestley and Taylor reference ET for water unlimited land surfaces [77] and a crop factor (K c ) based on an enhanced vegetation index (EVI). A global vegetation moisture index (GVMI) is used to account for non-optimal moisture conditions. This method is generally empirical, and aims to develop an ET dataset that GVMI is independent from land cover classification.
Another energy balance method is the Operational Simplified Surface Energy Balance (SSEBop, USGS EROS Center, Sioux Falls, SD, USA) [35] that employs a relationship between ET 0 and a land surface temperature-based scalar (ETrf) to express land wetness. An operational version of SSEB was proposed [78,79] by assimilating air temperature to account for the topographical and latitudinal heterogeneity impact on surface temperature. The SSEBop model defines the temperature scalar using hot and cold reference values for any pixel. The cold reference value is estimated as an empirically established fraction of the daily maximum air temperature; the hot reference value is obtained by adding a vertical temperature difference (dT) to the cold reference value [35].
Though progressing significantly in the past years and gaining maturity enough for calibration of hydrological models, remote sensing-based ET products still induce various uncertainties, e.g., [72,73,80]. These uncertainties come in various forms, such as from temperature and net radiation [35], and the aerodynamic component [81]. An uncertainty of up to 4 to 18% can be linked to sources of net radiation [35]. The largest uncertainty originates, however, from the impact of soil moisture on the regulation of the ET process. While thermal infrared measurements are excellent indicators of evaporative cooling [82], they are sensitive to cloud cover. Solutions based on microwave measurements are therefore suggested simultaneously for a long time [83,84]. Consequently, a suite of ET algorithms have been built during the last 25 years, and every algorithm will have its own model formulation and accuracy. Only physically-based algorithms with thorough validation in several environmental conditions should be considered for the calibration of hydrological models.
Because existing global scale ET products have different predictive capabilities, and there is no reliable ground truth dataset available in Day Basin to select any one of them, an ensemble ET product has been created in addition to the original datasets, on the basis of a simple linear average value for the Day Basin. The ensemble ET product used in this study is based on the combination of SEBS (5 km × 5 km), CMRSET (5 km × 5 km), SSEBop (1 km × 1 km), and MOD16 data (1 km × 1 km), and has a spatial resolution of 1 km × 1 km grid. A finer ET map is deemed necessary to assess water balances at HRUs spatial level. The same downscaling procedure as described in CMRSET using the Enhanced Vegetation Index (EVI), and Global Vegetation Moisture Index (GVMI) was applied. The Residual Moisture Index (RMI) was used to describe the impact of vegetation moisture content on the crop coefficient.
The selection criteria for a certain ET product were based on the hydrological consistency between the annual totals of precipitation (ΣP) and discharge (ΣQ) for a time span of 10 consecutive years (2003 to 2012). The stream flow data (ΣQ) from the 2 stations (Phu Ly and Ninh Binh) were used (see Figure 4). In the case of the Day Basin, there is an unknown withdrawal from the Red River, which makes a direct comparison weaker. The analysis of the Day Basin demonstrates that SEBS produces the highest ET values, and MOD16 the lowest. CMRSET and MOD16 performed similarly for annual and seasonal periods, and both were lower than the ensemble ET. Of all five ET datasets (4 individual and 1 ensemble product), the ensemble ET, SSEBop, and CMRSET delivered similar annual ET rates, averaged for the drainage area at Phu Ly (ΣET of 1073 mm, 1041 mm, and 1007 mm, respectively) and Ninh Binh (ΣET of 1103 mm, 1044 mm, and 1031 mm, respectively). This is because the ensemble ET compensated for the difference between higher end and lower end ET products. Of all four ET datasets, SSEBop gave the most similar results compared to the ensemble ET. Albeit minor, differences were spotted during the dry period. For the dry season, ET from MOD16 is comparable with the ensemble ET (ΣET of 361 mm and 359 mm for Phu Ly, 376 mm and 369 mm for Ninh Binh, respectively), while SSEBop tends to give lower ET. The performance of a certain actual evapotranspiration (ETact) algorithm is dependent on factors such as land use land cover type, climate, and the presence of mountains [8], meaning that the accuracy of ETact predictions will vary across the basin. The seasonal performance of the ensemble ET values mismatch, due to storage changes in the soil water balance and the regulating role of lakes and reservoirs on river discharge. During the dry season, ET was much higher than precipitation hence, the displayed axis scale differed from the yearly average and the seasonal wet period. The interim conclusion is that the ensemble ET product generated from linear average SEBS, CMRSET, SSEBop, and MOD16 provided accurate and most stable results for the Day Basin. Accordingly, the ensemble ET data was tested further before being used in the optimization process.  [8], meaning that the accuracy of ETact predictions will vary across the basin. The seasonal performance of the ensemble ET values mismatch, due to storage changes in the soil water balance and the regulating role of lakes and reservoirs on river discharge. During the dry season, ET was much higher than precipitation hence, the displayed axis scale differed from the yearly average and the seasonal wet period. The interim conclusion is that the ensemble ET product generated from linear average SEBS, CMRSET, SSEBop, and MOD16 provided accurate and most stable results for the Day Basin. Accordingly, the ensemble ET data was tested further before being used in the optimization process.

Crop Coefficient
Crop coefficient (Kc) was defined in [85] for unlimited soil water conditions, as the ratio of actual evapotranspiration to reference evapotranspiration (ETo). The Kc value for paddy rice is suggested in various studies such as by [86][87][88]. Kc for rice varies between 1.02 and 1.23 [89] throughout various growth stages. The cropping pattern in the Day Basin consists of mainly paddy rice (2 seasons from February to April/May, and from May/June to September) and some other crops (vegetables, September/October to January). The ensemble ET data is assumed to represent paddy rice as the dominant crop. The Kc derived from the ensemble ET and reference ET was tested to see if it falls within the range found in the literature.
In order to validate the accuracy of actual ET derived from satellites, reference evapotranspiration (ETo) was calculated using the FAO56 Penman-Monteith equation [85]. GLDAS is a good example of global standardized datasets, in this case, being climate related. Figure 5 illustrates the average monthly

Crop Coefficient
Crop coefficient (K c ) was defined in [85] for unlimited soil water conditions, as the ratio of actual evapotranspiration to reference evapotranspiration (ET o ). The K c value for paddy rice is suggested in various studies such as by [86][87][88]. K c for rice varies between 1.02 and 1.23 [89] throughout various growth stages. The cropping pattern in the Day Basin consists of mainly paddy rice (2 seasons from February to April/May, and from May/June to September) and some other crops (vegetables, September/October to January). The ensemble ET data is assumed to represent paddy rice as the dominant crop. The K c derived from the ensemble ET and reference ET was tested to see if it falls within the range found in the literature.
In order to validate the accuracy of actual ET derived from satellites, reference evapotranspiration (ET o ) was calculated using the FAO56 Penman-Monteith equation [85]. GLDAS is a good example of global standardized datasets, in this case, being climate related. Figure 5 illustrates the average monthly K c in the Day Basin. The maximum K c values are experienced in the wet season (June to July). K c is much lower in the dry months of January to February, because the crop intensity is strongly reduced by the lack of rainfall, and rice becomes fully dependent on irrigation water supply.
Water 2018, 10, x 9 of 21 Kc in the Day Basin. The maximum Kc values are experienced in the wet season (June to July). Kc is much lower in the dry months of January to February, because the crop intensity is strongly reduced by the lack of rainfall, and rice becomes fully dependent on irrigation water supply.  [88] with mid-season Kc reported to be around 1.23. The Kc during end-season (April to May for dry season rice and September for wet season rice) varies around 0.90 to 1.10, similar to crop coefficients estimated by [88,89]. Hence, the Kc values are falling within the acceptable range, indicating that the ensemble ET performed well for the Day Basin and can be used for SWAT model calibration.  The K c values are plotted with a 3-month Simple Moving Average (SMA) filter ( Figure 6). The K c values during different growing stages typically vary between 0.70 and 1.00 (initial stage), and 0.90 to 1.20 (mid-season stage). The differences in K c are due to the variety of cultivated paddy rice, as well as irrigation management. The K c in June to July and October to November 2008 was exceptionally high, since this was a very wet year. The mid-season K c reached 0.81 to 1.00, and 0.96 to 1.21, in March and June respectively. These values are similar to those of [88] with mid-season K c reported to be around 1.23. The K c during end-season (April to May for dry season rice and September for wet season rice) varies around 0.90 to 1.10, similar to crop coefficients estimated by [88,89]. Hence, the K c values are falling within the acceptable range, indicating that the ensemble ET performed well for the Day Basin and can be used for SWAT model calibration.  The Kc values are plotted with a 3-month Simple Moving Average (SMA) filter ( Figure 6). The Kc values during different growing stages typically vary between 0.70 and 1.00 (initial stage), and 0.90 to 1.20 (mid-season stage). The differences in Kc are due to the variety of cultivated paddy rice, as well as irrigation management. The Kc in June to July and October to November 2008 was exceptionally high, since this was a very wet year. The mid-season Kc reached 0.81 to 1.00, and 0.96 to 1.21, in March and June respectively. These values are similar to those of [88] with mid-season Kc reported to be around 1.23. The Kc during end-season (April to May for dry season rice and September for wet season rice) varies around 0.90 to 1.10, similar to crop coefficients estimated by [88,89]. Hence, the Kc values are falling within the acceptable range, indicating that the ensemble ET performed well for the Day Basin and can be used for SWAT model calibration.

Leaf Area Index
Leaf Area Index (LAI) is defined as the area of green leaf per ground area. It is an important variable for eco-hydrological modelling and quantifying ecosystem services [90]. LAI influences the evapotranspiration rate and its partitioning into transpiration (T), soil evaporation (E), and interception (I). At the same time, LAI determines the amount of Absorbed Photosynthetically Active Radiation (APAR), which determines the energy level for photosynthesis. MOD15-LAI data (NASA EOSDIS Land Processes DAAC, USGS Earth Resources Observation and Science (EROS) Center, Sioux Falls, SD, USA) with an 8-day temporal resolution was downloaded from open access databases. Monthly average LAI values with a spatial resolution of 1 km have been reconstructed for the period of 2005-2011.
The SWAT model estimates LAI values assuming a certain upper boundary function of growth that is corrected for stress factors (temperature, water, and nutrients). These stress factors can vary greatly and have an empirical character. The calibration of the model including LAI measurements is therefore crucial. The empirical LAI parameters to be prescribed in SWAT use an internal crop development database [53].

SWAT Parameters to Be Optimized
A summary of various sets of parameters used in model calibration of different processes, i.e., surface runoff, snow, plant growth etc. can be found in [54]. A similar list of parameters was also suggested by [16,42,59]. Calibration of ET is less common, because this dataset is usually not available. More recent studies, including [8,43,45], demonstrated that SWAT can be calibrated against spatial ET data as well. Key parameters for the calibration of ET and their value range were derived from the SWAT user's manual [53]. Since the purpose of this study is to calibrate the model for ET and LAI, the model parameters to be optimized were divided into two groups: ET and LAI. This grouping indicates which parameters are most sensitive to ET and LAI, and thus, suitable for optimizing the model performance for these two processes. The list is long because one parameter might control more than one process, e.g., the available water capacity of soil layers affects the generation of surface runoff, but simultaneously determines the amount of water that is evaporated (ET) and percolated to the underground. SUFI-2 automatically optimizes the selected parameters within their predefined range, hence this aspect of parameter selection is following the default guidelines of the tool.
Six parameters ESCO, EPCO, REVAPMN, SOL_K, SOL_AWC, and SOL_BD have been selected to optimize the ET simulations. The baseflow recession constant ALPHA_BF and SCS runoff curve number (CN2) were also included because of their influences on the surface-subsurface hydrological processes, and thus, on the water availability for evapotranspiration. Seven other parameters are identified that influence the leaf area index development, see [53]: BLAI, ALAI_MIN, DLAI, LAIMX1, LAIMX2, FRGRW1, and FRGRW2. SUFI-2 thus optimizes 15 model parameters for different sub-basins and HRUs. Since SWAT-CUP is unable to read spatial data, time series of ET and LAI was derived and used as input into the model. While ET process is more dependent on the local climate and land cover properties, LAI is strictly driven by the plant type, and hence, the calibration should be done at sub-basin, land cover, and HRU levels. The utilization of high resolution spatial datasets as input into HRU analysis in SWAT (i.e., land use, DEM, soil maps) ensures fair distribution of HRUs, and hence, spatially distributed ET and LAI across the watershed.
In the case of conventional calibration using discharge, two river discharge points (i.e., Phu Ly and Ninh Binh) were used (500 simulations for both calibration and validation) and the acquired result was used to compare against the output from calibration technique using remotely sensed data.
In the first step of calibration technique using remote sensing data, average ET and LAI for 119 sub-basins and all land-cover groups was extracted from ensemble ET and Moderate Resolution Imaging Spectroradiometer (MODIS) and used as observed data in SWAT-CUP. Derived time series of ET and LAI from each sub-basin and land cover type was used to calibrate 15 parameters in each sub-basin and land cover group individually (500 simulations each run for calibration and validation). This first estimation of ET and LAI ensures a good agreement in term of average ET, LAI, and subsequently, water balance. In the second batch of calibration, a number of HRUs from each sub-basin representing all land use classes and ensuring a total coverage of at least 50% the sub-basin's area were selected in the calibration (500 simulations each run for calibration and validation). The calibration result yields a set of soil and plant parameters which varies from sub-basins and land cover groups, and that reflects the local eco-hydrological processes, and is therefore highly valuable. This unprecedented spatial variability from SUFI-2 optimization cannot be obtained from soil surveys and soil maps, which makes this investigation interesting. The classical SUFI-2 calibration will give the same parameter value to all sub-basins and HRUs encompassed in a certain drainage area, and hence, the new system of calibration provides much more insight into the local characteristic of the soil and vegetation system.

Results and Discussions
The The monthly simulation for the entire Day Basin, as one bulk system, was presented in Figure 7. The simulated values related satisfactorily to the observed ensemble ET values. The ET peak value from remote sensing did not exceed 140 mm/month, while the modelled ET was as high as 160 mm/month. For remote sensing-based calibration, simulated ET in SWAT is underestimated by 5 to 10 percent in the irrigated agricultural land downstream, while the forest land showed good correlation. NSE yielded from 0.61 (calibration) to 0.65 (validation). R 2 for calibration and validation was 0.71, while Kling Gupta Efficiency (KGE) ranged from 0.80 to 0.83, respectively. The high values from KGE were due to the fact that low values were well represented for bias and correlation, as compared to NSE's reputation on overemphasizing peak values. By contrast, traditional discharge-based calibration performed statistically more inferior with NSE, which were only at −0.01 and −0.20, and KGE at 0.46 and 0.47, while R 2 was of significant values 0.77 and 0.76, for calibration and validation, respectively. Discharge-based calibration results also encountered a consistent underestimation of ET temporally as compared to remote sensing-based method. Another observation in both cases is that the lower ET values simulated during winter were always lower than the observed values from remote sensing, even though this difference was much lower when using ET and LAI as calibration inputs. This could be related to the dry period in which SWAT computes water stress due to a lack of soil moisture. This could suggest that the storage capacity of the soil in reality is higher, or it could also be related to lower vertical water fluxes between the top soil, sub-soil, and the unconfined shallow aquifer, or the sensitivity of vegetation to soil moisture. A very low reference ET 0 during winter could also be an explanation. SUFI-2 ensured, however, that the spatial patterns match rather well. Some local differences occur unavoidably due to the inaccuracy of the various mathematical expressions used to compute a complex hydrological process such as ET. Remote sensing ET values reflect more the real world conditions as they are based on observations [33][34][35]42]. The agreement between SWAT and remote sensing data was expressed by means of the correlation coefficient and the bias. Using spatial dataset in the calibration yields resulted in much higher agreement, and reveals that both the ET formulations in SWAT as well as the SUFI-2 optimization techniques are adequate. The same conclusion was drawn earlier by other researchers that validate SWAT on the basis of ET data series. Further to the evapotranspiration simulations of SWAT, the vegetation response to water can be evaluated by means of a comparison of the simulated LAI (see Figure 8) in two cases: (a) remote sensing-driven, and (b) discharge-driven calibration. The remote sensing-based calibration yielded acceptable NSE ranging from 0.50 for calibration and 0.57 for validation. R 2 for calibration and validation was 0.51 and 0.59, while Kling Gupta Efficiency (KGE) ranged from 0.59 to 0.75. Conventional discharge-based calibration performed poorly, in term of statistics, with NSE of only 0.35 and 0.43, KGE of −0.34 and 0.03, while R 2 was higher at 0.67 and 0.77 for calibration and validation, respectively. LAI values in discharge-based calibration were consistently underestimated compared with observation from MODIS. With regard to remote sensing-based calibration, the timing of the green cover development seems acceptable. The peak LAI values during Spring 2005 and 2006 are not accurate. This may be due to constancy of the LAI related calibration parameters for all the different simulation years. It would be better to make the maximum LAI parameter variable for each year, to enable it to better respond to years with weather anomalies. A good description of LAI evolution will improve the timing of rice emergence, which in turn, affects the irrigation and transpiration processes. Further to the evapotranspiration simulations of SWAT, the vegetation response to water can be evaluated by means of a comparison of the simulated LAI (see Figure 8) in two cases: (a) remote sensing-driven, and (b) discharge-driven calibration. The remote sensing-based calibration yielded acceptable NSE ranging from 0.50 for calibration and 0.57 for validation. R 2 for calibration and validation was 0.51 and 0.59, while Kling Gupta Efficiency (KGE) ranged from 0.59 to 0.75. Conventional discharge-based calibration performed poorly, in term of statistics, with NSE of only 0.35 and 0.43, KGE of −0.34 and 0.03, while R 2 was higher at 0.67 and 0.77 for calibration and validation, respectively. LAI values in discharge-based calibration were consistently underestimated compared with observation from MODIS. With regard to remote sensing-based calibration, the timing of the green cover development seems acceptable. The peak LAI values during Spring 2005 and 2006 are not accurate. This may be due to constancy of the LAI related calibration parameters for all the different simulation years. It would be better to make the maximum LAI parameter variable for each year, to enable it to better respond to years with weather anomalies. A good description of LAI evolution will improve the timing of rice emergence, which in turn, affects the irrigation and transpiration processes.  The total simulated ET is 5855 mm, compared with 5727 mm for the ensemble ET during the calibration, hence a difference of 128 mm (2.2%). The model performs consistently during the validation period with 5713 mm (simulated ET), as compared to 6015 mm (ensemble ET), leading to a difference of 302 mm (5.3%). In very general terms, the set of ET related equations in SWAT has a limited capacity to mimic the complex processes of soil evaporation, plant interception, and plant transpiration that occur in reality, due to the dynamic meteorological and hydrological processes. The results indicate that SUFI-2 succeeded in generating a close to reality ET in both calibration and validation of ET, both spatially and temporally.   The total simulated ET is 5855 mm, compared with 5727 mm for the ensemble ET during the calibration, hence a difference of 128 mm (2.2%). The model performs consistently during the validation period with 5713 mm (simulated ET), as compared to 6015 mm (ensemble ET), leading to a difference of 302 mm (5.3%). In very general terms, the set of ET related equations in SWAT has a limited capacity to mimic the complex processes of soil evaporation, plant interception, and plant transpiration that occur in reality, due to the dynamic meteorological and hydrological processes. The results indicate that SUFI-2 succeeded in generating a close to reality ET in both calibration and validation of ET, both spatially and temporally.  The total simulated ET is 5855 mm, compared with 5727 mm for the ensemble ET during the calibration, hence a difference of 128 mm (2.2%). The model performs consistently during the validation period with 5713 mm (simulated ET), as compared to 6015 mm (ensemble ET), leading to a difference of 302 mm (5.3%). In very general terms, the set of ET related equations in SWAT has a limited capacity to mimic the complex processes of soil evaporation, plant interception, and plant transpiration that occur in reality, due to the dynamic meteorological and hydrological processes. The results indicate that SUFI-2 succeeded in generating a close to reality ET in both calibration and validation of ET, both spatially and temporally.    To verify simulated streamflow from the remote sensing-driven calibration, the relationship between monthly river flow simulations and flow measurements at the Phu Ly and Ninh Binh stations was assessed, as shown in Figure 11. The discharges in several cross sections was constructed using the rating curve of measured flows and water levels during restricted periods, measured using acoustic Doppler current profilers [51].
The agreement between simulated river flow and station discharge measurements was expressed through R 2 , ranging from 0.71 (Phu Ly) to 0.78 (Ninh Binh), and Nash-Sutcliffe from 0.55 to 0.63, respectively. Considering that river discharge was computed from measured water level [51], the agreement between simulated river flow and station discharge measurements is good. Measured river flow data were not used in the calibration process, which shows that a good simulation of ET facilitates the prediction of stream flow. Note that the flow in these rivers is far from being natural, due to all the irrigation practices that are occurring. Measured river flow data were not used in the calibration process, which shows that a good simulation of ET will make it possible to calculate flows directly (streamflow and ET being the two largest components of the water balance) without having to optimize flow in the To verify simulated streamflow from the remote sensing-driven calibration, the relationship between monthly river flow simulations and flow measurements at the Phu Ly and Ninh Binh stations was assessed, as shown in Figure 11. The discharges in several cross sections was constructed using the rating curve of measured flows and water levels during restricted periods, measured using acoustic Doppler current profilers [51].
The agreement between simulated river flow and station discharge measurements was expressed through R 2 , ranging from 0.71 (Phu Ly) to 0.78 (Ninh Binh), and Nash-Sutcliffe from 0.55 to 0.63, respectively. Considering that river discharge was computed from measured water level [51], the agreement between simulated river flow and station discharge measurements is good. Measured river flow data were not used in the calibration process, which shows that a good simulation of ET facilitates the prediction of stream flow. Note that the flow in these rivers is far from being natural, due to all the irrigation practices that are occurring. Measured river flow data were not used in the calibration process, which shows that a good simulation of ET will make it possible to calculate flows directly (streamflow and ET being the two largest components of the water balance) without having to optimize flow in the calibration process. Hence, this shows that flow data can be replaced by ET to optimize SWAT model performance, and that the ET is reflecting actual anthropogenic processes of the water cycle.
Water 2018, 10, x 15 of 21 calibration process. Hence, this shows that flow data can be replaced by ET to optimize SWAT model performance, and that the ET is reflecting actual anthropogenic processes of the water cycle. To make the fit between remotely sensed precipitation and ET feasible, extra water has to be supplied to the land use class irrigated land in certain months. The total amount of irrigation water supply is computed to be 1.934 billion m 3 /y, and this amount of water is thus withdrawn from the Red River through various ungauged inlet points. This number can now be estimated with more precision because ET is known [16,38,43]. The main water intake for irrigation is between January and March (during the Winter-Spring paddy rice) with a maximum amount of 98 mm during February. During the Summer-Autumn paddy rice season, the water intake is concentrated in July with an amount of 113 mm per month. The total storage change ∆S for the entire Day Basin indicates the difference between all inflow and outflow terms. ∆S during 2003 to 2013 for unsaturated and saturated zones are 11.7 and 9.6 mm/y, respectively, due to the fact that there is water locally stored in lakes, streams, and also moving within the saturated layers to the deep aquifer. This is a realistic number, and confirms that the water balance of the total system is rather accurately simulated.
The total water balance for the basin in average rainfall, a dry (2007), and a wet (2008) year, is shown in Table 1. While the average annual rainfall was 1710 mm/y, it reached 2121 mm/y in 2008 and 1496 mm/y in 2007. The ET in 2007 and 2008, was 985 and 1010 mm/y, respectively, compared to 958 mm/y on average. This quasi-constancy of ET is noticed more often in other studies. The regulating role of soil water storage in the vadose zone, and the lower evaporative demand in wet years and higher demand during dry years, are some major causing factors for restriction of ET dynamics. Due to the impact of rainfall variation to overall water budget, the change in soil water content (∆SW) for 2007 was indeed negative (∆SW= −27 mm/y) while for 2008, it was positive (ΔSW = +80 mm/y). The final calibrated range of 15 parameters is expressed in Table 2. To make the fit between remotely sensed precipitation and ET feasible, extra water has to be supplied to the land use class irrigated land in certain months. The total amount of irrigation water supply is computed to be 1.934 billion m 3 /y, and this amount of water is thus withdrawn from the Red River through various ungauged inlet points. This number can now be estimated with more precision because ET is known [16,38,43]. The main water intake for irrigation is between January and March (during the Winter-Spring paddy rice) with a maximum amount of 98 mm during February. During the Summer-Autumn paddy rice season, the water intake is concentrated in July with an amount of 113 mm per month. The total storage change ∆S for the entire Day Basin indicates the difference between all inflow and outflow terms. ∆S during 2003 to 2013 for unsaturated and saturated zones are 11.7 and 9.6 mm/y, respectively, due to the fact that there is water locally stored in lakes, streams, and also moving within the saturated layers to the deep aquifer. This is a realistic number, and confirms that the water balance of the total system is rather accurately simulated.
The total water balance for the basin in average rainfall, a dry (2007), and a wet (2008) year, is shown in Table 1. While the average annual rainfall was 1710 mm/y, it reached 2121 mm/y in 2008 and 1496 mm/y in 2007. The ET in 2007 and 2008, was 985 and 1010 mm/y, respectively, compared to 958 mm/y on average. This quasi-constancy of ET is noticed more often in other studies. The regulating role of soil water storage in the vadose zone, and the lower evaporative demand in wet years and higher demand during dry years, are some major causing factors for restriction of ET dynamics. Due to the impact of rainfall variation to overall water budget, the change in soil water content (∆SW) for 2007 was indeed negative (∆SW= −27 mm/y) while for 2008, it was positive (∆SW = +80 mm/y). The final calibrated range of 15 parameters is expressed in Table 2.

Conclusions
The availability of spatially distributed precipitation, ET and LAI gridded data from open access-or partially open access-earth observation data platforms, makes it feasible to calibrate soil and vegetation process parameters of eco-hydrological models, also when rivers are ungauged and the water distribution system is complex. In this study, four individual ET models were averaged linearly to match the simulations of ET from SWAT. No in situ measurements were available to verify the performance of individual ET models. From inspection of the water balance, a simple linear average value gave best results as compared to streamflow and K c crop coefficients, described in the international literature. Nevertheless, it is required to undertake more studies on the ensemble ET product to provide further progressive insights on averaging of individual estimates.
Secondly, this paper demonstrates the capabilities of SWAT and the auto-calibration SUFI-2 to render bio-physical processes in data scarce basins in a distributed manner. A total of 15 essential bio-physical parameters of the unsaturated zone and exchange processes seem to be adequately calibrated in SWAT for each of the 7907 HRUs. Otherwise, there would not have been such a good agreement with the spatio-temporal variability of the remote sensing parameters, as in the case of conventional discharge-driven calibration. Furthermore, this level of spatial detail cannot be obtained from soil maps. The hydrological formulations in SWAT are thus adequate for simulating eco-hydrological processes. In the near future, remote sensing data on soil moisture, net primary production, and water quality will become available as well, which will undoubtedly further enrich the options to calibrate additional SWAT model parameters.
The approach proposed in this study using SWAT-CUP and SUFI-2 will improve the facilitation and standardization of calibration process for basins with scant field data. By optimizing evapotranspiration and photosynthesis for Hydrological Response Unit, swift estimates of surface runoff, erosion, groundwater recharge, baseflow, storage changes, withdrawals, and carbon assimilation can be determined and used to quantify ecosystems services. The availability of the system parameters will allow future predictions of the basin water cycle in response to external factors, such as climate and land-use changes, and computing scenarios for green growth, i.e., conservation plans, reduction of greenhouse gas emissions, etc.