Development of an Evapotranspiration Data Assimilation Technique for Streamflow Estimates: A Case Study in a Semi-Arid Region

Streamflow estimates are substantially important as fresh water shortages increase in arid and semi-arid regions where evapotranspiration (ET) is a significant contribution to the water balance. In this regard, evapotranspiration data can be assimilated into a distributed hydrological model (SWAT, Soil and Water Assessment Tool) for improving streamflow estimates. The SWAT model has been widely used for streamflow estimations, but the applications combining SWAT and ET products were rare. Thus, this study aims to develop a SWAT-based evapotranspiration data assimilation system. In particular, SWAT is gridded at Hydrologic Response Unit (HRU) level to incorporate gridded ET products acquired from the remote sensing-based ETMonitor model. In the modeling case, Gridded SWAT (GSWAT) shows a good agreement of streamflow modeling with the original SWAT. Such a scant margin between them is due to the modeling domain mismatch caused by different HRU delineations. In the ET assimilation case, we carry out a synthetic data experiment to illustrate the state augmentation Direct Insertion (DI) method and a real data experiment for the upper Heihe River Basin. The results demonstrate the benefits of the ET assimilation for improving hydrologic processes representations. In the future, more remotely sensed data can be assimilated into the data assimilation system to provide more reliable hydrological predictions.


Introduction
Highly accurate forecasting of streamflow is critical for water resource management as fresh water shortages increase in arid and semi-arid regions [1]. As an important component of the water budget, it is also critical for sustainable crop production, food security and social sustainability [2][3][4]. With the data assimilation method and remotely sensed data, streamflow estimates can be improved a lot [5][6][7][8][9]. However, most of these studies used remotely sensed soil moisture or snow; little attention has been focused on assimilating observed evapotranspiration (ET) data. In particular, evapotranspiration is the largest mode of water loss in arid and semi-arid regions [2]. Streamflow estimates will benefit a lot from taking into account observed evapotranspiration from remote sensing in these areas.
Remote sensing-derived ET products have been used to evaluate or calibrate the hydrological model and improve streamflow forecasting [10][11][12][13]. However, such calibrations need a lot of time and labor or computing resources. Additionally, calibration is insufficient for streamflow forecasting when forcing input (i.e., precipitation) is not correct. Data assimilation as a complementary tool can be used to fill up the calibration deficiency in decreasing the discrepancy between streamflow prediction and observation. Pan et al. [14] assimilated Moderate Resolution Imaging Spectroradiometer (MODIS)-derived ET into Variable Infiltration Capacity (VIC) and Surface Energy Balance System (SEBS) models to estimate water budget terms. However, ET was only a non-state variable in the model, which could not feed the assimilation effect back to the model and update the model's state variables, so the other hydrological variables were not optimized as a whole and this method is not suitable for improving streamflow prediction.
ET-related data assimilation researches are still in the formative stage for most hydrological models [15,16]. For model state updating purposes, ET observation should be associated with the model state variable; considering ET as a state variable in a model (e.g., [15]) or associating it with other state variables are two practical methods. Remotely sensed evapotranspiration can be used for hydrological model parameter calibration or estimation, which can affect the surface runoff generation [10,12,13,16,17]. Parr et al. [18] found that hydrological predictions may be subject to the bias in past ET modeling, leading to substantial uncertainties in the modeled hydrological trend. They replaced the simulated ET components (canopy evaporation, transpiration, bare ground evaporation, canopy sublimation, and surface sublimation) of the VIC model with the bias-corrected ET calculated from a simple bias-correction algorithm and historic remotely sensed ET. The streamflow estimation was improved due to the corrected ET components. This procedure can still be considered as a model bias correction method. Understanding historic evapotranspiration trends and high-quality forcing data are prerequisites, because the mass balance is the main constrain for ET components and runoff updating in VIC [18]. Generally, monthly time scale water balance, including ET and streamflow, can be simulated well by well calibrated hydrological models. They are particularly valuable for applications where one is primarily interested in monthly, seasonal and annual streamflow volumes [19]. While the flood risk warning and management need hydrological data at high spatiotemporal resolution, utilizing fresh incoming high-resolution ET data for hydrological forecasting simultaneously is still a big challenge.
It is a major challenge to relate surface runoff modeling to the ET modeling. The reports and literatures on the relationship among precipitation, ET and surface runoff will be useful for associating daily ET and surface runoff [20]. Rawlins et al. [21] found that simulations of arctic evapotranspiration and runoff are both strongly dependent on the quality of time series data used to drive the model. The uncertainty of model input (i.e., precipitation) can be the main factor to influence the simulation accuracy of both ET and surface runoff [22]. The uncertainty of precipitation can be defined by a widely used multiplicative error model [23,24]. An estimated precipitation error can be obtained by analyzing observed and simulated ET. The error for simulated ET can be transformed to precipitation error based on the ET model. So using remotely sensed evapotranspiration to reduce input uncertainties and further to improve the hydrologic predictions is practical.
In this paper, we assess the performance of streamflow predictions with remote sensing-derived ET and state augmentation direct insertion method for a distributed hydrological model, SWAT. We have to solve two problems: (1) matching ET data with the modeling units efficiently; and (2) assimilating the ET into the model to update the streamflow estimates. In this way, a modified SWAT model and a state augmentation direct insertion assimilation method are designed to assimilate ET data.
The remainder of the paper is organized as follows. Sections 2 and 3 provide background on hydrological models. Section 4 describes the direct insertion data assimilation method. We describe case studies including the design of the synthetic and real-word data assimilation experiments in Section 5. Section 6 presents the results from two designed cases: (1) hydrological modeling scheme; and (2) evapotranspiration data assimilation scheme. Finally, conclusions and outlooks are summarized in Section 7.

Soil and Water Assessment Tool (SWAT)
The SWAT model is a semi-distributed hydrological model that includes detailed representations of hydrological processes at the watershed scale. The model calculates vertical and horizontal water Sustainability 2017, 9, 1658 3 of 21 movement through three spatial levels (HRU level, subbasin level and basin level). There are two hydrologic cycle modeling phases (land phase and routing phase) in SWAT. In Land phase, snow melt is calculated based on a linear function of the difference between the average snow pack maximum air temperature and the base or threshold temperature for snow melt. Daily surface runoff is calculated based on Soil Conservation Service curve number (SCS-CN) method. Vegetation growth and canopy intercept is calculated based on plant growth function and LAI. Potential and actual evapotranspiration (evaporation from ground and vegetation, transpiration from vegetation) can be calculated based on Penman-Monteith method (used in this study), Priestley-Taylor method or Hargreaves method. Water percolates through several soil layers and ultimately becomes aquifer recharge. Some of it moves laterally in the soil profile and contribute to streamflow. Water leaves groundwater storage primarily by discharge into rivers or lakes. In the routing phase, water routing from reaches and channels to basin outlet is calculated using Muskingum or the variable storage coefficient routing method. See Arnold, et al. [25] for a detailed description of the SWAT model. SWAT has been widely used for regional water budget estimations around the word for decades [26,27].

Gridded SWAT (GSWAT)
The modeling units in SWAT are Hydrologic Response Unit (HRU), subbasin and basin from fine to large. HRU is the combination of unique soil, land cover and topography characters. HRUs may have various sizes and shapes and no geospatial locations. The HRU attributes will cost users a lot of work to process the grid format remote sensing data to match their spatial location and irregular size. We need to match HRU with grid data efficiently. More user-friendly modeling units will save SWAT users a lot of time and labor for incorporating grid format data. We gridded SWAT through two main procedures: data gridded and model gridded. The comparison between SWAT and GSWAT is shown in Table 1. The software of GSWAT is available at http://hexiayouxi.tk.

Data Gridded
We gridded HRU in SWAT, which can be called GSWAT. Each grid cell is an individual HRU in GSWAT. Forcings and parameters for HRUs are the same for the corresponding grid cells. An HRU parameters database is generated from the SWAT database (e.g., SWAT.mdb). HRU identification number raster (hruid.txt) and subbasin identification raster (subid.txt) are generated based on the HRU and subbasin maps. Only the area and shape of HRU are changed to match with grid format data. The model will link the HRU index with the corresponding HRU parameters when reading input data. Therefore, there are no original HRU level input files (e.g., *.hru). The number of HRUs is determined by the number of cells within the watershed. The grid spatial dimension and resolution of the watershed is defined in the control file (file.cio). Forcing data are derived from the weather stations from SWAT. More details of the modification of SWAT can be seen in [28,29]. We can match gridded ET data with gridded HRUs spatially. Spatial extent and location between gridded ET and gridded HRU should be consistent.

Model Gridded
Vertical and horizontal water balance equations in SWAT are retained (Equation (1)) for GSWAT. There are two source soil water balance equation which separates evaporation and transpiration in surface and root zone [30]. The hydrologic process as simulated by SWAT and GSWAT at each HRU is based on the single source water balance equation [31]: where SW t is the soil water content above the wilting point in soil profile at the end of the day t, SW o is the initial soil water content. P j is the amount of the precipitation including rainfall and snowfall on Sustainability 2017, 9, 1658 4 of 21 day j, SURQ j is the amount of surface runoff on day j, ET a,j is the amount of actual evapotranspiration on day j, S j is the amount of water entering the vadose zone from the soil profile and G j is the amount of return flow on day j. All units are mm H 2 O. Soil water content can be expressed as volume fraction based on the thickness of each layer. The conventional SWAT model code is modified to incorporate the gridded forcing data and parameter data. The HRU related parameter variables are redefined to include the grid location dimension, for example, soil water available variable (sol_awc), whose HRU number dimension size is modified to ncolums × nrows. Ncolums and nrows are grid cells' column and row respectively. HRU id (hruid) variable is added to imply variable's spatial location. At this point, the original jth HRU's sol_awc(:,j) variable is modified to sol_awc(:,hruid(j)). Grid data reading/writing subroutines (e.g., readsub1.f, readhru1.f, pmeas2.f) are added into the code to read/write ASCII/NETCDF format data (e.g., hruid.txt).

State Augmentation Direct Insertion Data Assimilation
According to the water balance equation (Equation (1)), the water balance components are correlated with each other. During both ET and surface runoff simulations, the precipitation uncertainty plays a leading role. When rainfall occurs, water is partitioned into three main components: surface runoff, evapotranspiration and infiltration. Surface runoff is calculated based on SCS-CN conceptual method. The physically surface runoff hydrological process is not represented well [32].
Based on the data assimilation concepts [33], the observed and modeled hydrological variables relationships calculation equations can be written as follows: where P o is model input (observed) precipitation, ε p is the precipitation error, λ is the multiplicative error coefficient, ET o is the observed ET, ET m is the modeled ET. ε e is the model error for ET process, δ e is the ET measurement error, SURQ o is updated surface runoff (corresponds to observed surface runoff), SURQ m is modeled surface runoff. ε q is the model error in surface runoff process, δ q is the surface runoff measurement error. Due to the precipitation having multiplicative error most possibly [23], we define the precipitation error as a multiplicative error (see Equation (2)). The multiplicative error model extracts the systematic errors more cleanly, which is more consistent with the large variability of precipitation measurements. When the observed data has low bias and deserves to be trust (δ e = 0, δ q = 0), we can assume that the discrepancy between observed data and simulated data is due to the model uncertainty (see Equations (3) and (4)). When the model has been calibrated and validated well against long term observations, model system errors could be neglected, therefore we can assume that the model uncertainties for evapotranspiration and surface runoff are dominated by the input (precipitation) uncertainty. Additionally, precipitation, evapotranspiration and surface runoff are highly correlated with each other [34]. The error coefficient for simulated ET and surface runoff can be linearly correlated. We can also assume the model error is partially represented by simulated ET and surface runoff and precipitation error: where λ et and λ surq is the multiplicative error coefficient for simulated ET and surface runoff, α et and α q is the multiplicative error coefficient for precipitation error, representing the precipitation error contribution to the simulated ET and surface runoff. Here, we simply assume the errors sourced from precipitation error contribution equally (α et = α q ). Along with observed ET, surface runoff can be estimated with the error coefficient calculated from simulated ET. We can redeploy the water balance components based on the proportional discrepancy between observed and simulated water balance components. With observed and modeled ET data, we can define the uncertainty and calculate the correct coefficient. We already assumed that the proportional discrepancy between observed ET and simulated ET is also suitable for surface runoff variable. Then we multiply modeled surface runoff variable with the correct coefficient to acquire unbiased surface runoff. The basin states will be updated along with the hydrologic processes subsequently. This data assimilation method can be termed as state augmentation direct insertion (DI) method.

Implementation Procedures
The flowchart of the implementation procedures of GSWAT coupled with DI is shown in Figure 1. As shown in Figure 1, in each time step gridded HRU simulation, when GSWAT model detects available ET data, it will replace modeled ET value with incoming ET data at the end of the ET calculation. PET and surface runoff will be updated according to the correct coefficient derived from observed ET (ET o ) and modeled ET ET m (see Equations (3) and (4)). Finally, the updated HRU runoff will be accumulated to the subbasin level and enters routing phase. Streamflow will be simultaneously updated by the updated water storage in each reach.

Study Area and Data
The Heihe River Basin is the second largest inland river basin in China ( Figure 2). The upper Heihe River Basin is around 10,000 km 2 and located above 1700 m elevation. This watershed is located in the semi-arid area. The upper Heihe River provides 80.2% of the total river flow to the middle and lower reach. Water from glacial melt only accounted for 3.6% of the river flow. Precipitation is the main water sources in this region. The annual precipitation varies in the range 200-500 mm and 60% of the total precipitation occurred in the summer. The annual potential evapotranspiration is over 600 mm. Runoff from the upper Heihe River maintains industrial and agricultural production in the midstream and ecological system heath of the downstream [35]. There are two major tributaries originate from

Study Area and Data
The Heihe River Basin is the second largest inland river basin in China ( Figure 2). The upper Heihe River Basin is around 10,000 km 2 and located above 1700 m elevation. This watershed is located in the semi-arid area. The upper Heihe River provides 80.2% of the total river flow to the middle and lower reach. Water from glacial melt only accounted for 3.6% of the river flow. Precipitation is the main water sources in this region. The annual precipitation varies in the range 200-500 mm and 60% of the total precipitation occurred in the summer. The annual potential evapotranspiration is over 600 mm. Runoff from the upper Heihe River maintains industrial and agricultural production in the midstream and ecological system heath of the downstream [35]. There are two major tributaries originate from the Qilian Mountain in east and the west. Yingluoxia hydrological station is the outlet of the upstream watershed. The dominant land cover types are mountain forest and alpine meadow.
All input data including forcing data, parameter data for SWAT and GSWAT are all acquired from WestDC (http://westdc.westgis.ac.cn). They are resampled to 1 km for both models. Streamflow observations obtained from the Yingluoxia station are used to calibrate SWAT model parameters (e.g., CN2). Details of calibration can be seen in Zhang et al. [36]. The meteorological data (air temperature, precipitation, relative humidity, wind speed, and solar radiation) and calibrated parameters (e.g., soil texture and conductivity) are both used for SWAT and GSWAT. The SWAT model grid size is 30m based on the DEM resolution. The grid resolution for GSWAT is 1 km, which will be consistent with the following evapotranspiration data. All input data including forcing data, parameter data for SWAT and GSWAT are all acquired from WestDC (http://westdc.westgis.ac.cn). They are resampled to 1 km for both models. Streamflow observations obtained from the Yingluoxia station are used to calibrate SWAT model parameters (e.g., CN2). Details of calibration can be seen in Zhang et al. [36]. The meteorological data (air temperature, precipitation, relative humidity, wind speed, and solar radiation) and calibrated parameters (e.g., soil texture and conductivity) are both used for SWAT and GSWAT. The SWAT model grid size is 30m based on the DEM resolution. The grid resolution for GSWAT is 1 km, which will be consistent with the following evapotranspiration data.
Evapotranspiration data was produced based on remote sensing-based ETMonitor model [37] and WRF (Weather Research and Forecasting) model data by Cui and Jia [38]. With a variety of parameterization methods of the aerodynamic and bulk surface resistances, the estimation of soil evaporation and vegetation transpiration was based on the Shuttleworth-Wallace Dual-Source model. The model contains a variant form of Penman-Monteith equation that combines energy balance with plant physiology and soil water diffusion based on the assumption of aerodynamic mixing occurring at a mean canopy source height within the canopy. The rainfall interception was estimated using a reformulated remote sensing-based Gash model, which is forced by remote sensing observations of canopy structure (e.g., LAI and vegetation storage capacity) and rainfall intensity. The interception of vegetation was used to replace the interception of canopy, and a Poisson distribution function was applied to the LAI value of each pixel. A recommended method proposed by the World Meteorological Organization (WMO) was used to estimate the sublimation of snow and ice using the formula developed by Kuzmin [39]. The evaporation from water surface was estimated using the classical Penman equation [40]. Evapotranspiration data was produced based on remote sensing-based ETMonitor model [37] and WRF (Weather Research and Forecasting) model data by Cui and Jia [38]. With a variety of parameterization methods of the aerodynamic and bulk surface resistances, the estimation of soil evaporation and vegetation transpiration was based on the Shuttleworth-Wallace Dual-Source model. The model contains a variant form of Penman-Monteith equation that combines energy balance with plant physiology and soil water diffusion based on the assumption of aerodynamic mixing occurring at a mean canopy source height within the canopy. The rainfall interception was estimated using a reformulated remote sensing-based Gash model, which is forced by remote sensing observations of canopy structure (e.g., LAI and vegetation storage capacity) and rainfall intensity. The interception of vegetation was used to replace the interception of canopy, and a Poisson distribution function was applied to the LAI value of each pixel. A recommended method proposed by the World Meteorological Organization (WMO) was used to estimate the sublimation of snow and ice using the formula developed by Kuzmin [39]. The evaporation from water surface was estimated using the classical Penman equation [40].
The ETMonitor model combines three different actual ET parameterizations for water body, snow/ice surface and soil-vegetation canopy. It utilizes microwave and optical remote sensing observations including MODIS albedo product (MCD43), MODIS vegetation index product (MOD13A2), soil moisture (AMSR-E), MODIS land cover type product (MCD12Q1), snow cover (MOD10A2) and precipitation (Tropical Rainfall Measuring Mission (TRMM)) as input data to estimate the daily ETp for all sky conditions (see Equation (7)). The saturated and residual soil water contents of soil surface layer for the ETMonitor were provided by the China dataset of soil hydraulic parameters [41], which are derived from multiple Pedotranfer Functions, including the parameters in the Clapp and Hornberger Functions (FCH) and in the van Genuchten and Mualem Functions (FGM), respectively. The ETMonitor model combines three different actual ET parameterizations for water body, snow/ice surface and soil-vegetation canopy. It utilizes microwave and optical remote sensing observations including MODIS albedo product (MCD43), MODIS vegetation index product (MOD13A2), soil moisture (AMSR-E), MODIS land cover type product (MCD12Q1), snow cover (MOD10A2) and precipitation (Tropical Rainfall Measuring Mission (TRMM)) as input data to estimate the daily ET p for all sky conditions (see Equation (7)). The saturated and residual soil water contents of soil surface layer for the ETMonitor were provided by the China dataset of soil hydraulic parameters [41], which are derived from multiple Pedotranfer Functions, including the parameters in the Clapp and Hornberger Functions (FCH) and in the van Genuchten and Mualem Functions (FGM), respectively.
where ET p is the evapotranspiration product. M is ETMonitor model.
The daily actual evapotranspiration of the Heihe River Basin for the years 2009-2011 was estimated with the ETMonitor at a spatial resolution of 1 km. This product was evaluated using eddy covariance (EC) flux observations (The data were collected from the Watershed Allied Telemetry Experimental Research (WATER) and provided by the Cold and Arid Regions Sciences Data Center at Lanzhou (available at http://westdc.westgis.ac.cn/) at field scale and compared with other ET product (MOD16A2) at regional scale with smaller RMSE (low to 0.39 mm/day) and higher R 2 (up to 0.96). This ET product presents spatial and temporal consistent with other product. Merging such lots of remote sensing data and models provides us a reliable ET observation.

Hydrological Modeling Scheme
To evaluate the performance of hydrological modeling, we carried out Yingluoxia watershed experiment with SWAT and GSWAT. SWAT was used to simulate daily streamflow within the upper Heihe River Basin between 1990 and 2008. The previous SWAT calibration and validation work is based on 1990-1999 and 2000-2009 Yingluoxia hydrological station observation data [36]. Twenty-two sensitive parameters were selected to be calibrated against the streamflow observations based on the previous studies. The basin and subbasin level parameters for SWAT were automatically calibrated using Particle Swarm Optimization method with 100 iterations. We used the hydroPSO R package to conduct the automatic model calibration which is proposed by Zambrano-Bigiarini and Rojas [42], and HRU level parameters were manually calibrated (see Table 2). To exclude the influence of parameters source, we used the already calibrated parameters from SWAT for GSWAT. During the SWAT calibration period, NSE (Nash Sutcliffe Efficiency) and R 2 were 0.86 and 0.92 respectively for the simulated monthly streamflow. During the SWAT validation period, NSE and R 2 were 0.88 and 0.86 respectively. Then, the GSWAT was driven by the same forcings derived from SWAT to simulate upper Heihe River Basin variables from 1990 to 2008.

Synthetic Experiment Design
Our primary goal of synthetic data experiment is to test the data assimilation method. We designed a synthetic "twin" experiment and specifications for different cases are shown in Table 3. In the TRUE case, we first ran GSWAT to simulate upper Heihe river basin throughout the year 2000 and considered the results as true states. We conducted many years (1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)  assimilation method can really diminish model errors under different forcing patterns in the real data experiment. The ET result was selected as ET observation. The precipitation was selected as true precipitation. For data assimilation and open loop simulation, the precipitation uncertainty was defined as a multiplicative log normally distributed error (defined in Equations (8)-(11)) [43]; In the DA and NODA cases, GSWAT was driven by the perturbed precipitation (P perturbed ) to simulate basin states and fluxes with and without ET assimilated throughout year 2000. Finally, we analyzed streamflow results derived from this experiment.
where e is normally distributed random variable with 1 mean and standard deviation of 0.25. P is precipitation, and s is the variance scaling factor of precipitation with a set value of 0.5 in this study.

Real Data Experiment Design
The real data experiment design is similar to the synthetic experiment except the observation data. The observation data in the real data experiment is the ET product (see Section 5.1). The parameters derived from SWAT may be not appropriate for GSWAT, the streamflow forecasting can be further improved by an additional calibration. So we added extra calibration procedures for GSWAT modeling and ET data assimilation cases. More details of the test plan for real data experiment are shown in Table 4. Table 4. Test plan for real data experiment.

CASE Model Method
: used in the case.

Results and Discussion
As noted above, our work is centered on assimilating remote sensing-based ET data into the modified SWAT model. In particular, we first assess the modified model by comparing it with the original model; then we assess the performance of ET data assimilation with the modified SWAT model through synthetic and real data experiments. All cases are tested in the upper Heihe River Basin located in northwest China.

Hydrological Modeling Experiment
Streamflow, especially flood peak rate estimations from GSWAT, is a little higher than SWAT under the same forcing conditions (Figure 3). That is due to some parameters, e.g., time for flow from farthest point in subbasin to enter a channel whose calculation is based on HRU and subbasin area, derived from SWAT are not most suitable for GSWAT. The area of watershed delineated by GSWAT and SWAT are 9674 km 2 and 9759 km 2 respectively. Area mismatch between GSWAT and SWAT is over 85 km 2 . Another reason is that when irregular HRUs (Figure 4a) are outputted to grid cells (Figure 4b), scale issue will cause parameter mismatch [28,44]. The parameter linkages across scales from original HRU to gridded HRU will produce bias [44]. The discrepancy between GSWAT and SWAT is mainly caused by system bias (HRU delineation) [45]. It can be diminished through model calibration. Monthly and daily streamflow trend estimations from GSWAT and SWAT agree well with each other (R 2 = 0.99 and R 2 = 0.96) (Figure 3). The overestimation of peak value can be corrected with proper parameters for GSWAT at daily time scale. As high correlation exists between GSWAT and SWAT, parameters for GSWAT can be calibrated as well as SWAT against streamflow observation. Then the GSWAT model will have the ability to simulate watershed streamflow as well as SWAT. Figure 5 plots averaged annual basin hydrological variables obtained from SWAT and GSWAT from 1991 to 2008. Precipitations are the same for both models. Snow water equivalents are sublimated in both cases. GSWAT presents overestimations of surface runoff, groundwater flow and recharge, percolation and Potential ET (PET), underestimations of Lateral flow and evapotranspiration (ET) compared with SWAT. As noted above, that is due to parameter mismatching after grid HRU delineation. In the end, GSWAT produces more water yield than SWAT. For streamflow predictions, unsuitable parameters (e.g., slope length and CN2) related to HRU area and shape attributes for GSWAT can be calibrated against observations individually instead of using parameters calibrated based on SWAT.   Previously, the parameters for GSWAT are derived from calibrated SWAT. The parameters may not be appropriate for GSWAT. In this case, HRU level parameters (i.e., CN2) for GSWAT were calibrated against streamflow observations at Yingluoxia station from 1991 to 2008 individually. We calibrated the HRU level parameters by percentage manually. In validation period (2009) Table 5). With the same simple calibration procedure, calibrated GSWAT (GSWATCALI) can present a better performance than calibrated SWAT for streamflow prediction.    Previously, the parameters for GSWAT are derived from calibrated SWAT. The parameters may not be appropriate for GSWAT. In this case, HRU level parameters (i.e., CN2) for GSWAT were calibrated against streamflow observations at Yingluoxia station from 1991 to 2008 individually. We calibrated the HRU level parameters by percentage manually. In validation period (2009) Table 5). With the same simple calibration procedure, calibrated GSWAT (GSWATCALI) can present a better performance than calibrated SWAT for streamflow prediction.
We have tested another grid version of SWAT (SWATgrid) and found that there is little difference (The correlation coefficient between the two data sets is 0.998) between SWATgrid and SWAT [46]. So we just compared our model (GSWAT) with original SWAT. Additionally, SWATgrid is an interface to prepare the input files (e.g., *sub, *.hru) for running SWAT in a grid subbasin-based We have tested another grid version of SWAT (SWATgrid) and found that there is little difference (The correlation coefficient between the two data sets is 0.998) between SWATgrid and SWAT [46]. So we just compared our model (GSWAT) with original SWAT. Additionally, SWATgrid is an interface to prepare the input files (e.g., *sub, *.hru) for running SWAT in a grid subbasin-based scheme. There is a large amount of input files (the input files are determined by the number of grid cells), which will affect the application of SWATgrid in large-area and high-resolution modeling and data assimilation, which is our main focus in our future studies. So we developed our own grid-based SWAT model and also incorporated a parallel technique (OpenMP) into it to speedup our model [29].

Synthetic Experiment
The result from GSWAT with ET assimilated (DA) is closer to true streamflow than GSWAT without data assimilation (NODA) (Figure 6). Figure 6 shows that DA decreased RMSE from 83.57 m 3 /s to 28.13 m 3 /s and MAE from 36.05 m 3 /s to 10.40 m 3 /s and increased NSE from 0.41 to 0.77 and R 2 from 0.98 to 0.99. In streamflow estimation, DA performs better than NODA with lower RMSE, MAE and higher NSE, R 2 , especially for peak value (flood) estimation. It is essential to consider these forecasting accuracy when developing and implementing water quantity deployment actions in a watershed. Streamflow estimation error in DA is lower than that in NODA. The uncertainty in surface runoff modeling was successfully decreased by ET observation integration. The ET observation can correct not only ET modeling uncertainty but also surface runoff modeling uncertainty. The state augmentation direct insertion ET assimilation integration proved to be a successful method to improve the streamflow prediction for GSWAT. ET can also be used to constrain SWAT calibration [47]. The effectiveness of decreasing model errors with ET data assimilation method is similar to calibrating hydrological model with ET measurement. Both of them use additional measurement to correct water cycle components of the model based on their own accuracies. In synthetic calibration experiment, the calibration with accurate ET measurements is more resistant to errors and a robust way to produce accurate streamflow predictions [12].
From Figures 7 and 8, we can see that precipitation bias (Pbias equal to perturbed precipitation minus true precipitation), ET bias (ETbias) and surface runoff bias (SQbias) are spatially heterogeneous. The spatial patterns among ET correct coefficient, ET observation, ET bias and precipitation bias are related to each other. Surface runoff bias is decreased by ET correct coefficient especially in high bias region (Figures 7f and 8f). The most improvement for surface runoff estimation corresponds to the overlay of high ETbias and high SQbias (Figure 8b,c,f). In this synthetic data experiment case, precipitation uncertainty dominates the uncertainty in ET modeling and surface runoff modeling. ET observation is used to calculate the ET modeling bias correction coefficient, which is used to correct the bias in surface runoff generation modeling. The correction coefficient for surface runoff modeling is equal for ET modeling based on our assumption. This kind of assumption is the simplification of the relationship between uncertainties in ET and surface runoff modeling. However, the precipitation uncertainty for the ET modeling and surface runoff modeling being the same is not proper. This simplified assumption is suitable for modeling bias correction when the uncertainty relationship between ET and surface runoff modeling is within the range of linear uncertainty distribution. More work on the uncertainty relationship among precipitation, ET and surface runoff are needed in the future. Historic ET data can be used to correct the ET components estimations, which affect surface runoff generation in VIC [18]. The ET components estimations in GSWAT model have little effect on surface runoff generation. We cannot use the corrected ET components to update surface runoff generation in GSWAT. So in this study, we correct surface runoff modeling bias in GSWAT model by using state augmentation assimilation method to integrate ET observation directly. The correction coefficient calculated from modeled and observed ET for ET modeling was augmented for modeled surface runoff. The simultaneous updating of surface runoff with the method we provided and real-time ET observation is worth looking forward to. Besides the rule-based insertion method, the ensemble-based data assimilation method can also be used to deal with the correlated model errors. For example, Hartanto et al. [48] assimilated satellite-based actual evapotranspiration in a distributed hydrological model to improve the streamflow predictions based on a particle filter data assimilation framework; the hydrological states from the model runs (particles) that simulate actual ET close to the remotely sensed actual ET map are then chosen as initial conditions for the next simulation. Based on a hydrological model which expresses ET as state variables, Zou et al. [49] assimilate ET into Distributed Time Variant Gain Model (DTVGM) using ensemble Kalman filter (EnKF) to obtain more accurate soil moisture and streamflow.
Sustainability 2017, 9,1568 13 of 21 calibrating hydrological model with ET measurement. Both of them use additional measurement to correct water cycle components of the model based on their own accuracies. In synthetic calibration experiment, the calibration with accurate ET measurements is more resistant to errors and a robust way to produce accurate streamflow predictions [12]. From Figures 7 and 8, we can see that precipitation bias (Pbias equal to perturbed precipitation minus true precipitation), ET bias (ETbias) and surface runoff bias (SQbias) are spatially heterogeneous. The spatial patterns among ET correct coefficient, ET observation, ET bias and precipitation bias are related to each other. Surface runoff bias is decreased by ET correct coefficient especially in high bias region (Figures 7f and 8f). The most improvement for surface runoff estimation corresponds to the overlay of high ETbias and high SQbias (Figure 8b,c,f). In this synthetic data experiment case, precipitation uncertainty dominates the uncertainty in ET modeling and surface runoff modeling. ET observation is used to calculate the ET modeling bias correction coefficient, which is used to correct the bias in surface runoff generation modeling. The correction coefficient for surface runoff modeling is equal for ET modeling based on our assumption. This kind of assumption is the simplification of the relationship between uncertainties in ET and surface runoff modeling. However, the precipitation uncertainty for the ET modeling and surface runoff modeling being the same is not proper. This simplified assumption is suitable for modeling bias correction when the uncertainty relationship between ET and surface runoff modeling is within the range of linear uncertainty distribution. More work on the uncertainty relationship among precipitation, ET and surface runoff are needed in the future. Historic ET data can be used to correct the ET components estimations, which affect surface runoff generation in VIC [18]. The ET components estimations in GSWAT model have little effect on surface runoff generation. We cannot use the corrected ET components to update surface runoff generation in GSWAT. So in this study, we correct surface runoff modeling bias in GSWAT model by using state augmentation assimilation method to integrate ET observation directly. The correction coefficient calculated from modeled and observed ET for ET modeling was augmented for modeled surface runoff. The simultaneous updating of surface runoff with the method we provided and real-time ET observation is worth looking forward to. Besides the rule-based insertion method, the ensemble-based data assimilation method can also be used to deal with the correlated model errors. For example, Hartanto et al. [48] assimilated satellite-based actual evapotranspiration in a distributed hydrological model to improve the streamflow predictions based       Figure 9 plots daily streamflow observations at the Yingluoxia station and predictions obtained from GSWAT, GSWAT with ET assimilated (GSWATET), calibrated SWAT (SWAT), calibrated GSWAT (GSWATCALI) and calibrated GSWAT with ET assimilated (GSWATCALIET). The calibration is used to diminish some system biases (parameters biases) in models. Table 5 summarizes the statistical fit criteria for daily streamflow predictions from these cases. We get the following results and deductions:

Real Data Experiment
(1) Calibration integration: GSWATCALI performs better than GSWAT and GSWATCALIET performs better than GSWATET. For example, the RMSE for GSWATCALI is less than for the GSWAT (82 versus 97 m 3 /s). The calibrated SWAT parameters may not be suitable for GSWAT daily streamflow estimations. Daily scale streamflow estimation can be further improved by calibrating model parameters based on daily streamflow observation. With or without ET assimilated, calibration procedure can still further diminish parameter (system) errors in models. (2) ET assimilation integration: GSWATET shows clearly improved streamflow peak rates estimations compared with GSWAT; ET assimilation decreases RMSE from 97 m 3 /s to 78 m 3 /s. After being calibrated, GSWATCALIET also presents a better streamflow prediction than GSWATCALI; ET assimilation decreased RMSE from 63 m 3 /s to 57 m 3 /s. Model system bias has little effect on ET assimilation. With or without calibration procedure, ET assimilation can still further improve streamflow prediction for GSWAT.
Streamflow estimation from GSWATCALIET is closest to streamflow observation. GSWATCALIET is better than both GSWATCALI and GSWATET. Calibration and ET assimilation are equally important in streamflow prediction in this study. Calibration can diminish model system bias (parameter bias) [50] and data assimilation can correct both model input bias and model parameter bias [51]. Observation and model uncertainties cause discrepancy between streamflow prediction and observation. There is a very high peak streamflow prediction from all models and low peak streamflow observation, for example at 18 August. These are several types of errors in streamflow observation, such as instrumental error, representative error and artificial error (i.e., miscalculation and parallax). The errors in streamflow observation are mostly temporally consistent and cannot suddenly increase so much. So we deduce that the most possibility of causing such large discrepancy is that precipitation data has a large bias during those days. Calibration cannot correct this kind of bias by adjusting model parameters based on the observed streamflow, but ET assimilation can diminish the effects of these kinds of errors on streamflow estimation in GSWAT.
Sustainability 2017, 9,1568 16 of 21 representative error of weather stations. Precipitation input affects discharge estimates more than calibrating model parameters and assimilating remote sensing soil moisture and ET [54].  Figure 10 shows the monthly and annual hydrological variables (precipitation (PREC), surface runoff (SURQ), lateral flow (LATQ), groundwater flow (GWQ), percolation (PERC), soil water [55], evapotranspiration (ET), potential evapotranspiration (PET) and water yield (YIELD)) estimations from GSWAT and GSWATET. GSWATET shows underestimations of ET (observed), PET, SURQ and YIELD compared with GSWAT. As observed ET data is lower than GSWAT modeled, surface runoff in GSWATET is corrected to be lower than GSWAT modeled. Therefore, the water yield in GSWATET is also lower than that in GSWAT and closer to observed streamflow. The large corrected water volume is in summer where high rainfall rate and peak streamflow occurred ( Figure 10). In this study, the evapotranspiration is overestimated in GSWAT. GSWATET corrects this kind of error by simply replacing it with observed data and updates surface runoff with the correct coefficient calculated from modeled and observed ET. Through the routing phase modeling, the model summarizes the updated states (water storage in reach) and fluxes (evapotranspiration and surface runoff) to the outlet of the basin. Finally, the updated streamflow can be closer to the observed streamflow at the outlet of the basin.
ET assimilation updates the basin variable ET and SURQ, removes some estimation errors in GSWAT ET and SURQ modeling. In GSWAT, ET includes evaporation from rivers and lakes, bare soil, and vegetative surfaces; evaporation from within the leaves of plants (transpiration); and sublimation from snow and ice surfaces. The potential evapotranspiration is calculated based on the Penman-Monteith equation in this study. Potential soil water evaporation is estimated as a function of potential evapotranspiration and Leaf Area Index (LAI). Actual soil water evaporation is estimated by using exponential functions of soil depth and water content. There are large uncertainties (such as LAI and soil water content errors) in evapotranspiration modeling in GSWAT. The bias in modeled ET can be attributed to imperfect modeling processes for various land cover types, parameterizations and forcing input uncertainties etc. The flood peak is impacted by the magnitude and distribution of precipitation, antecedent soil moisture condition and surface runoff generation in the basin. In addition to long-term RMSE, R, MAE, streamflow results should be evaluated based on their ability to capture peak streamflow values. In 2009, observed peak streamflow does not correspond well to precipitation that indicates that the forcing data has large bias. The high-quality ET product is derived from merged multiple remote sensing data including TRMM precipitation data, WRF data and comprehensive ET model. Merged forcing data has lower bias than that from a single source [56,57]. However, there are still some discrepancies between rainfall events and streamflow observations. For SWAT and GSWAT, the daily streamflow simulation is good (NSE > 0.50) for most years. The daily streamflow simulation results from SWAT and GSWAT in 2009 are bad (NSE < 0.50) due to that the weather stations dramatically overestimate the sub-watershed rainfall during 2009 (MAE > 30 m 3 /s). For example, some of the streamflow peaks are not corresponding well to the rainfall events (see Figure 9), because the rainfall only occurred near the weather station and no rain occurred in most areas of the sub-watershed. The rainfall data representative error is too large in 2009; the NSE value will be not satisfied whenever how hard we try to calibrate the model. This kind of errors is obvious in the original SWAT, which uses one station's data to represent the whole sub-watershed when data are scarce. GSWAT model is originally developed to represent HRU level weather conditions; we can integrate more spatial representative rainfall data (i.e., radar data [52,53]) into the model to decrease this kind of errors. In this paper, we just showed the effectiveness of the ET assimilation method in decreasing rainfall volume error, which cannot handle the spatial representative error of weather stations. Precipitation input affects discharge estimates more than calibrating model parameters and assimilating remote sensing soil moisture and ET [54]. Figure 10 shows the monthly and annual hydrological variables (precipitation (PREC), surface runoff (SURQ), lateral flow (LATQ), groundwater flow (GWQ), percolation (PERC), soil water [55], evapotranspiration (ET), potential evapotranspiration (PET) and water yield (YIELD)) estimations from GSWAT and GSWATET. GSWATET shows underestimations of ET (observed), PET, SURQ and YIELD compared with GSWAT. As observed ET data is lower than GSWAT modeled, surface runoff in GSWATET is corrected to be lower than GSWAT modeled. Therefore, the water yield in GSWATET is also lower than that in GSWAT and closer to observed streamflow. The large corrected water volume is in summer where high rainfall rate and peak streamflow occurred ( Figure 10). In this study, the evapotranspiration is overestimated in GSWAT. GSWATET corrects this kind of error by simply replacing it with observed data and updates surface runoff with the correct coefficient calculated from modeled and observed ET. Through the routing phase modeling, the model summarizes the updated states (water storage in reach) and fluxes (evapotranspiration and surface runoff) to the outlet of the basin. Finally, the updated streamflow can be closer to the observed streamflow at the outlet of the basin.
ET assimilation updates the basin variable ET and SURQ, removes some estimation errors in GSWAT ET and SURQ modeling. In GSWAT, ET includes evaporation from rivers and lakes, bare soil, and vegetative surfaces; evaporation from within the leaves of plants (transpiration); and sublimation from snow and ice surfaces. The potential evapotranspiration is calculated based on the Penman-Monteith equation in this study. Potential soil water evaporation is estimated as a function of potential evapotranspiration and Leaf Area Index (LAI). Actual soil water evaporation is estimated by using exponential functions of soil depth and water content. There are large uncertainties (such as LAI and soil water content errors) in evapotranspiration modeling in GSWAT. The bias in modeled ET can be attributed to imperfect modeling processes for various land cover types, parameterizations and forcing input uncertainties etc. The flood peak is impacted by the magnitude and distribution of precipitation, antecedent soil moisture condition and surface runoff generation in the basin. In addition to long-term RMSE, R, MAE, streamflow results should be evaluated based on their ability to capture peak streamflow values. In 2009, observed peak streamflow does not correspond well to precipitation that indicates that the forcing data has large bias. The high-quality ET product is derived from merged multiple remote sensing data including TRMM precipitation data, WRF data and comprehensive ET model. Merged forcing data has lower bias than that from a single source [56,57]. This kind of ET product integrating low bias forcing data and comprehensive ET modeling process is better than observed or modeled ET alone [37]. With this kind of ET data, GSWATET can correct errors in ET modeling and other related variables modeling. In the end, GSWATET presents a superior streamflow prediction than GSWAT. Besides calibration, assimilating a reliable ET product into gridded SWAT improves the streamflow prediction in the upper Heihe River Basin. This kind of ET product integrating low bias forcing data and comprehensive ET modeling process is better than observed or modeled ET alone [37]. With this kind of ET data, GSWATET can correct errors in ET modeling and other related variables modeling. In the end, GSWATET presents a superior streamflow prediction than GSWAT. Besides calibration, assimilating a reliable ET product into gridded SWAT improves the streamflow prediction in the upper Heihe River Basin.

Summary and Outlooks
Streamflow is an important water resource for sustainable watershed management. To promote the accuracy of the streamflow estimation in a semi-arid basin, the evapotranspiration product was assimilated into the gridded SWAT (GSWAT) model using the direct insertion method. We first modified SWAT to incorporate grid format data and assessed its capability in hydrological modeling.

Summary and Outlooks
Streamflow is an important water resource for sustainable watershed management. To promote the accuracy of the streamflow estimation in a semi-arid basin, the evapotranspiration product was assimilated into the gridded SWAT (GSWAT) model using the direct insertion method. We first modified SWAT to incorporate grid format data and assessed its capability in hydrological modeling. Parameters were further calibrated to improve daily streamflow prediction for SWAT and GSWAT. We related the ET modeling uncertainty (mainly precipitation uncertainty) to surface runoff modeling and updated it with incoming ET observation at daily time scale. The observed ET along with modeled ET updated surface runoff by using correct coefficient. Correct coefficient was calculated based on the assumption of the precipitation uncertainty for surface runoff and evapotranspiration modeling is correlated. This kind of assumption can be further improved with deeper investigation. Updated fluxes enter into the reach and update the water storage, which will route to the outlet of the basin finally. A synthetic data experiment and a real data experiment were carried out for a semi-arid basin to validate the efficacy of the state augmentation direct insertion method. The results show that the state augmentation direct insertion (DI) assimilation method can achieve an acceptable streamflow estimation improvement with a reliable ET observation, which derived from the combination of multisource remote sensing data and ET model with multiple parameterizations [37].
Parameters for models (SWAT, GSWAT and GSWATET) are calibrated against daily and monthly streamflow observations from 1990 to 1999 and validated from 2000 to 2009. Spatial domain and parameter mismatch between gridded and original HRU dominate the small discrepancies between basin variables simulated by GSWAT and SWAT. Spatial resolution of GSWAT can be matched with SWAT but will cost a lot of computational time. The detailed comparison between GSWAT and SWAT has been studied and will not be discussed in greater detail here [28]. Toward perfect streamflow prediction, calibration and data assimilation are all necessary. When the calibration is not sufficient, the data assimilation method can be used as a complement tool to diminish forecasting errors sourced from model uncertainties. In this study, parameter optimization and data assimilation procedures are separate. The simultaneous optimization and data assimilation (SODA) method can be used for improving hydrological forecasting [58,59].
Not only streamflow but also other hydrological variables (e.g., soil moisture and snow water equivalent) can be updated by assimilating ET product into the model. We did not include this kind of study in this article to briefly highlight the possibility of improving streamflow estimation with ET assimilated. Different ET products generated from different methods (e.g., [60]) can also be assimilated into the model to evaluate their effects on regional water balance simulation. Besides ET products, GSWAT can also assimilate other remote sensing data such as snow state products and soil moisture products from MODIS, AMSR-E, SMOS (Soil Moisture and Ocean Salinity) and SMAP (Soil Moisture Active and Passive) etc. Remotely sensed precipitation data can also be assimilated into GSWAT to decrease the forcing input uncertainty. The use of remote sensing data for hydrological model data assimilation has important implications to the catchments with sparse or no gauging. With more and more available remote sensing data, data assimilation technology will be a promising method to merge these data to provide more accurate land surface states and fluxes products for assisting watershed management and flood/draught prediction.