# ESA CCI Soil Moisture Assimilation in SWAT for Improved Hydrological Simulation in Upper Huai River Basin

Authors: Yongwei Liu; Wen Wang; Yuanbo Liu
Journal: Advances in Meteorology
 (2018)

---

## Abstract

The assimilation of satellite soil moisture (SM) products with coarse resolution is promising in improving rainfall-runoff modeling, but it is largely impacted by the data assimilation (DA) strategy. This study performs the assimilation of a satellite soil moisture product from the European Space Agency (ESA) Climate Change Initiative (CCI) in a physically based semidistributed hydrological model (SWAT) in the upper Huai River basin in China, with the objective to improve its rainfall-runoff simulation. In this assimilation, the ensemble Kalman filter (EnKF) is adopted with full consideration of the model and observation error, the rescaling technique for satellite SM, and the regional applicability of the hydrological model. The results show that the ESA CCI SM assimilation generally improves the streamflow simulation of the study catchment. It is more effective for low-flow simulation, while for very high-flow/large-flood modeling, the DA performance shows uncertainty. The less-effective performance on large-flood simulation lies in the relatively low dependence of rainfall-runoff generation on the antecedent SM as during which the SM is nearly saturated and the runoff is largely dominated by precipitation. Besides, the efficiency of DA is deteriorated by the dense forest coverage and the complex topography conditions of the basin. Overall, the ESA CCI SM assimilation improves the streamflow simulation of the SWAT model in particular for low flow. This study provides an encouragement for the application of the ESA CCI SM in water management, especially over low-flow periods.


---

## Body

## 1. Introduction

Soil moisture (SM) significantly impacts the rainfall-runoff process as it dominates the partitioning of precipitation into infiltration, runoff, and evaporation. In recent years, a large body of studies have been implemented to explore the approaches to improving rainfall-runoff modeling via enhancing SM estimation [1–10].One promising approach to improving SM estimation in turn improving rainfall-runoff modeling is to integrate the observed SM into the hydrological modeling process using data assimilation (DA) techniques [11–16]. In general, the SM data for integration can be obtained from field measurements and satellite observations. The in situ measurements are insufficient in the availability and spatial representativeness due to the high spatial heterogeneity of SM. Major researches on in situ SM assimilation focus on discussing the DA approaches and exploring the potential of SM assimilation in improving the hydrological process [13, 17, 18]. However, the satellite observations are capable of capturing the spatial distribution and temporal dynamics of SM on large scales. Despite the fact that the satellite remote sensing (RS) can only detect the surface SM information with a few centimeters (∼5 cm), it could represent the fastest response of SM dynamics to meteorological conditions [19]. A large number of studies have been implemented to assimilate the RS SM in the land surface model for the purpose of obtaining a more accurate and reliable profile SM data set on a regional or global scale [20–26]. Nevertheless, the assimilation of coarse-scale RS SM in the hydrological model targeted at improving the rainfall-runoff process is implemented in relatively few studies [10, 27–31]. Currently, there is still no consensus on the improvement of streamflow modeling through satellite soil moisture assimilation [4, 7]. For instance, almost no improvement of streamflow simulation was obtained by Brocca et al. [4] in the assimilation of the surface ASCAT SM retrievals, while up to 10–30% improvements were achieved in such other studies as Massari et al. [32], Lopez et al. [7], and Loizu et al. [10]. The large discrepancies of the DA performance in previous studies are likely due to the fact that it is influenced by various factors in the DA framework setup, such as the quantification of the model and observation error, the mismatch between the observed and simulated SM, the data quality and rescaling technique for RS SM, and the model physical mechanism and its regional applicability. To date, the added value of satellite soil moisture data in hydrological modeling is still underexplored [5, 32]. The performance of RS soil moisture assimilation in streamflow modeling presents certain specificity on the satellite data itself, the hydrological model, and the different configuration schemes in the DA framework setup. Therefore, specific studies on satellite soil moisture assimilation with comprehensive consideration of the DA implementation strategies are essential for exploring the significance of satellite soil moisture in hydrological modeling.In this paper, a case study for satellite soil moisture assimilation is implemented in the upper Huai River basin in China, with full consideration of the factors in the DA framework including the quantification of the model and observation error, the rescaling technique for RS SM, and the regional applicability of the hydrological model. This data assimilation is performed in a physically based semidistributed hydrological model (SWAT) based on a robust sequential data assimilation approach (the ensemble Kalman filter (EnKF)). A multisatellite-merged soil moisture data set from the European Space Agency (ESA) Climate Change Initiative (CCI) is adopted as the assimilation data source. The main objective of this study is to explore the potential of coarse-scale RS soil moisture in improving runoff modeling and to provide recommendations on the assimilation strategy.


## 2. Study Area and Data Used

### 2.1. Study Area

The study catchment is located in the upstream basin of the Huaibin hydrologic station in the Huai River basin, China (Figure1). The watershed covers about 16000 km2. The whole watershed is located in the transition region between the northern subtropical zone and the warm temperate zone. Its annual average rainfall is around 900 mm, 50%–80% of which falls during June to September. Here, the annual average temperature is about 15°C. The major land cover is the agriculture land (AGRC 32.5%, RICE 35%) and forest (FRST 23.6%) (Figure 2).Figure 1
The upper Huai River basin: location (the inset map), elevation, digital river network, location of the meteorological and hydrological stations (ST), and the subbasin delineation in SWAT model building.Figure 2
Land use/land cover (LU/LC) classes and soil distribution for the upper Huai River basin.


### 2.2. Data for SWAT Model

SWAT model building requires meteorological and underlying surface data. The meteorological data mainly include precipitation, maximum and minimum temperature, solar radiation, wind speed, and relative humidity. The precipitation data are drawn from 106 local rainfall stations within the catchment (Figure1). The other five meteorological data come from the three meteorological gauges (Xinyang, Gushi, and Guangshui) (Figure 1). The underlying surface data are the digital elevation (DEM), land cover, and soil category data. The DEM data are downloaded from the Shuttle Radar Topography Mission with a spatial resolution of 90 m (http://datamirror.csdb.cn/index.jsp). The land use/land cover (LU/LC) data are collected from the Chinese Cold and Arid Regions Science Data Center (http://westdc.westgis.ac.cn/) with a spatial resolution of 1 km (Figure 2). The soil data are resampled from a soil map at a scale of 1 : 100000 from the Soil Handbook of Henan Province. The soil for the whole catchment is divided into seven categories (Figure 2). The soil texture and its corresponding United States Department of Agriculture (USDA) classification for each category are shown in Table 1. Besides, there are six hydrologic stations (Dapoling, Changtaiguan, Zhuganfu, Xixian, Huangchuan, and Huaibin) (Figure 1) with daily streamflow measurements of 1992–2008 (the data quality issue exists for the years 2000 and 2001) in this basin.Table 1
Soil classification and its area proportions in the upper Huai River basin.
Soil category
Clay
Silt
Sand
Rock
USDA soil texture
Soil layer (mm)
Area proportion (%)
(%)
(%)
(%)
(%)
Cugutu
7.05
34
35.9
23.1
Sandy loam
0-50-150
7.35
Huanghetu
23.43
65.1
11.5
0
Silt loam
0-50-150-400-800
21.52
Huangzongrang
17.03
39.4
43.6
0
Loam
0-50-150-400-800
13.7
Huichaotu
12.86
51.8
35.3
0
Silt loam
0-50-150-400-800
4.72
Shajiangheitu
20.32
65.5
14.2
0
Silt loam
0-50-150-400-800
9.81
Shizhitu
9.16
44.1
46.7
0
Loam
0-50-150
7.8
Shuidaotu
16.46
71
12.5
0
Silt loam
0-50-150-400-800
34.62


### 2.3. ESA CCI Soil Moisture Data

The ESA CCI soil moisture data are a merged multisatellite surface soil moisture product developed in the Climate Change Initiative (CCI) by the European Space Agency (ESA). It combines the soil moisture retrievals from four microwave radiometers (SMMR, SSM/I, TMI, and AMSR-E) and two scatterometers (AMI and ASCAT) into a 0.25° global daily data set over 30 years from 1978. The data integration relies on their respective sensitivity to vegetation density and uses a Noah GLDAS-1 surface soil moisture product [33] as a climatology reference [34]. The ESA CCI SM consists of active, passive, or combined products. The active/passive products are the integration of the scatterometer/radiometer-based SM retrievals, respectively, while the combined product is the fusion of both the active and passive products. In this study, the combined product (ESA CCI SM v03.2) is adopted for soil moisture assimilation.


## 2.1. Study Area

The study catchment is located in the upstream basin of the Huaibin hydrologic station in the Huai River basin, China (Figure1). The watershed covers about 16000 km2. The whole watershed is located in the transition region between the northern subtropical zone and the warm temperate zone. Its annual average rainfall is around 900 mm, 50%–80% of which falls during June to September. Here, the annual average temperature is about 15°C. The major land cover is the agriculture land (AGRC 32.5%, RICE 35%) and forest (FRST 23.6%) (Figure 2).Figure 1
The upper Huai River basin: location (the inset map), elevation, digital river network, location of the meteorological and hydrological stations (ST), and the subbasin delineation in SWAT model building.Figure 2
Land use/land cover (LU/LC) classes and soil distribution for the upper Huai River basin.


## 2.2. Data for SWAT Model

SWAT model building requires meteorological and underlying surface data. The meteorological data mainly include precipitation, maximum and minimum temperature, solar radiation, wind speed, and relative humidity. The precipitation data are drawn from 106 local rainfall stations within the catchment (Figure1). The other five meteorological data come from the three meteorological gauges (Xinyang, Gushi, and Guangshui) (Figure 1). The underlying surface data are the digital elevation (DEM), land cover, and soil category data. The DEM data are downloaded from the Shuttle Radar Topography Mission with a spatial resolution of 90 m (http://datamirror.csdb.cn/index.jsp). The land use/land cover (LU/LC) data are collected from the Chinese Cold and Arid Regions Science Data Center (http://westdc.westgis.ac.cn/) with a spatial resolution of 1 km (Figure 2). The soil data are resampled from a soil map at a scale of 1 : 100000 from the Soil Handbook of Henan Province. The soil for the whole catchment is divided into seven categories (Figure 2). The soil texture and its corresponding United States Department of Agriculture (USDA) classification for each category are shown in Table 1. Besides, there are six hydrologic stations (Dapoling, Changtaiguan, Zhuganfu, Xixian, Huangchuan, and Huaibin) (Figure 1) with daily streamflow measurements of 1992–2008 (the data quality issue exists for the years 2000 and 2001) in this basin.Table 1
Soil classification and its area proportions in the upper Huai River basin.
Soil category
Clay
Silt
Sand
Rock
USDA soil texture
Soil layer (mm)
Area proportion (%)
(%)
(%)
(%)
(%)
Cugutu
7.05
34
35.9
23.1
Sandy loam
0-50-150
7.35
Huanghetu
23.43
65.1
11.5
0
Silt loam
0-50-150-400-800
21.52
Huangzongrang
17.03
39.4
43.6
0
Loam
0-50-150-400-800
13.7
Huichaotu
12.86
51.8
35.3
0
Silt loam
0-50-150-400-800
4.72
Shajiangheitu
20.32
65.5
14.2
0
Silt loam
0-50-150-400-800
9.81
Shizhitu
9.16
44.1
46.7
0
Loam
0-50-150
7.8
Shuidaotu
16.46
71
12.5
0
Silt loam
0-50-150-400-800
34.62


## 2.3. ESA CCI Soil Moisture Data

The ESA CCI soil moisture data are a merged multisatellite surface soil moisture product developed in the Climate Change Initiative (CCI) by the European Space Agency (ESA). It combines the soil moisture retrievals from four microwave radiometers (SMMR, SSM/I, TMI, and AMSR-E) and two scatterometers (AMI and ASCAT) into a 0.25° global daily data set over 30 years from 1978. The data integration relies on their respective sensitivity to vegetation density and uses a Noah GLDAS-1 surface soil moisture product [33] as a climatology reference [34]. The ESA CCI SM consists of active, passive, or combined products. The active/passive products are the integration of the scatterometer/radiometer-based SM retrievals, respectively, while the combined product is the fusion of both the active and passive products. In this study, the combined product (ESA CCI SM v03.2) is adopted for soil moisture assimilation.


## 3. Methodology

### 3.1. Soil and Water Assessment Tool (SWAT)

The SWAT is a physically based semidistributed watershed model, which has been widely used in rainfall-runoff modeling over recent years [35, 36]. In hydrological modeling, the catchment is firstly delineated into several subbasins according to its topography. Then, each subbasin is further divided into several hydrological response units (HRUs) based on the land use, soil, and slope. HRUs are basic calculation units for the land phase of the hydrologic cycle, on which the processes for surface runoff, lateral flow, and ground water are generated accompanied by evapotranspiration and soil water routing.Soil moisture lies in the center of the hydrologic cycle and makes different impacts on the above process. The initial profile soil water content influences surface runoff generation through the curve number in the SCS method [37]. After surface runoff generation, the water infiltrated to the soil profile is redistributed based on a storage routing technique with the soil water field capacity as the threshold. The water balance for each soil layer can be expressed as follows:(1)SWly′=SWly+Δwperc,ly−Qlat,ly−Ea,ly,where SWly and SWly′ are the soil water content (mm) at the start and end of the day, Δwperc,ly is the net percolation from the overlying layer (i.e., the layer ly+1), Qlat,ly is the lateral flow generated from the layer ly, and Ea,ly is the evapotranspiration drawn from the layer ly.The evapotranspiration from soil mainly includes two parts: soil evaporation and plant uptake/transpiration. In the SWAT, the potential evapotranspiration is firstly calculated using the Penman–Monteith equation [38]. Based on the potential evapotranspiration, the leaf area index, and the aboveground biomass and residue conditions, both the demand for transpiration/plant uptake and the demand for soil evaporation are determined. Then, the soil evaporation demand and the plant uptake demand for each soil layer are estimated using a depth distribution function. Finally, relying on the soil evaporation demand and the plant uptake demand with the available soil water as a constraint, the actual soil water evaporation and plant uptake are determined. In the processes mentioned above, the actual soil water extraction of a given layer is not allowed to be compensated by the extraction from other layers. However, the soil water deficiency can be made up by adjusting the soil compensation (esco) and plant compensation (epco) factors via changing the depth distribution of the soil evaporation demand and the plant water uptake demand. Besides, the calculation for the soil water percolation (wperc,ly) and lateral flow (Qlat,ly) is [39] omitted here.In the water routing phase, the SWAT adopts a storage feature to calculate the surface runoff and lateral flow generated from each HRU to the main channel, while it applies a linear reservoir similar technique to account for the ground water to the main channel. The channel water routing is performed using a variable storage routing method [40].


### 3.2. The Ensemble Kalman Filter (EnKF) for Soil Moisture Assimilation

The EnKF is a sequential DA approach evolved from the standard Kalman filter [41]. It is based on an ensemble of model states produced by adding the Monte Carlo noise to model forcing and states and/or parameters to approximate the model state error covariance matrix for the purpose of optimally merging the model predictions with observations.The state ensemble forecast at timet can be expressed as follows:(2)Xtf=FXt−1u,ut,δ+wtwt∼N0,σs2,where Xtf is the forecasted state ensemble at time t and Xt−1u is the updated state ensemble at t−1. In this study, it is constructed by the profile SM with up to four layers (Table 1) for all HRUs of the study basin (the HRUs delineation is detailed in Section 3.1). ut represents the model forcing inputs. In this study, it mainly includes the observed precipitation P and temperature T at each site. The precipitation error is assumed to be independent both in time and in space; that is, both the autocorrelation between time steps at each rainfall station and the error correlation among different stations are ignored. The perturbation (ηp) to precipitation is assumed to be a lognormal multiplicative distribution with mean 1 and covariance σp2 (3). The perturbation to temperature (ηT) is assumed to be an additive normal distribution with mean 0 and covariance σT2 (4). Besides, δ represents the model parameter with a perturbation (ηpar) of normal multiplicative distribution of mean 1 and covariance σpar2 (5). wt is the stochastic perturbation to the forecasted SM, being assumed to be an additive normal distribution with mean 0 and covariance σs2:(3)Pe=P⋅ηp,ηp∼lnN1,σp2,(4)Te=T+ηT,ηT∼N0,σT2,(5)δe=δ⋅ηpar,ηpar∼N1,σpar2.The state update for soil moisture can be obtained by(6)Xtu=Xtf+KtZt−HXtf,where Zt is the observation ensemble at time t. It is constructed by the RS SM for all grids covering the basin and being stochastically perturbed by an additive normal distribution with mean 0 and covariance σR2. H is the observation operator, being used to map the model states to the observations. It is constructed by the area proportions of HRUs in RS grids as the SWAT model-simulated SM is on the HRU level, while the observed SM is on RS grids. Kt is the Kalman gain, which is calculated based on the forecast and observation error covariance:(7)Kt=Pms,tPs,t+Rs,t−1,where Pms,t is the cross-error covariance between the predicted SM (Xtf) and the measurement prediction HXtf at time t, Ps,t is the error covariance of the measurement prediction at time t, and Rs,t is the error covariance of the RS SM at time t.


### 3.3. Implementation of the ESA CCI SM Assimilation in SWAT

To set up the ESA CCI SM data assimilation in the SWAT model, the SWAT model for the study catchment should be firstly built up. The model applicability is evaluated over the model calibration and validation processes. Then, the model error for the validated SWAT model is estimated based on the in situ streamflow observations at the Huaibin hydrologic station. Besides, the observation error estimation for the ESA CCI SM is implemented before the bias correction on the basis of the model-simulated soil moisture. A flow chart for implementing ESA CCI SM assimilation in SWAT is shown in Figure3.Figure 3
Flow chart for implementing ESA CCI SM assimilation in SWAT.


#### 3.3.1. SWAT Calibration

The SWAT model for the upper Huai River basin is built up based on the meteorological forcing and land surface data (Section2.2). This catchment is partitioned into 37 subbasins (Figure 1) and 146 hydrological response units (HRUs). To simplify and improve the model calibration process, a sensitivity analysis of the model parameters to the hydrologic modeling is implemented using the Latin-hypercube and one-factor-at-a-time method [42, 43]. Thereafter, the model is calibrated and validated using the daily runoff records over 1992–1999 and 2002–2008 at the interior and outlet hydrologic stations (i.e., Dapoling, Changtaiguan, Zhuganfu, Xixian, Huangchuan, and Huaibin in Figure 1), respectively. The parameter optimization is achieved by a combination of the autocalibration using the Sequential Uncertainty Fitting (SUFI2) [44] with the Nash–Sutcliffe coefficient of efficiency (detailed in Section 3.4) as the objective criteria and the manual fine-running method.


#### 3.3.2. Model Error Estimation

The determination of the model error is significant for the performance of DA as the model predictions and observations are merged based on the relative weight between the model and observation error (as in (6) and (7)). In this study, the model error is mainly contributed by the model input error for precipitation and temperature, the model parameter error for parameters sensitive to SM simulation (e.g., the available soil water capacity), and the model state error for simulated SM. The various errors mentioned above are characterized by additive/multiplicative normal/lognormal distribution specified in Section 3.2. The assumed distribution is only controlled by the standard deviation (SD), that is, σp, σT, σpar, and σs in (3), (4), (5), and (2), respectively. Therefore, the model error estimation is to quantify the model error parameters σp, σT, σpar, and σs.The quantification ofσp, σT, σpar, and σs is performed by analyzing the statistical characteristics of the simulated streamflow ensemble driven by model error perturbations on the basis of the observed streamflow at in situ sites. If the ensemble spread after perturbation is too large, the overfitting of observation exists in DA. Otherwise, the observed information cannot be fully utilized in DA. Therefore, the two ensemble verification measures (8) and (9) should be satisfied [32]. That is, if the ensemble spread sp is large enough, the temporal mean of the ensemble skill sk¯ should be similar to the temporal average of the ensemble spread sp¯:(8)sp¯sk¯≅1.And the observation should be indistinguishable from a member of the ensemble (N is the ensemble size):(9)skmse=N+12N,where(10)sp¯=1T∑kT1N∑i=1NQki−Q¯k2,sk¯=1T∑k=1TQ¯k−Qobs,k2,sk=1T∑k=1TQ¯k−Qobs,k,mse=1T∑k=1T1N∑i=1NQki−Qobs,k21/2,where Qki is the model simulated streamflow of the ensemble member i at time k, Q¯k is the ensemble mean of the model simulated streamflow, Qobs,k is the observed streamflow at time k, and T is the total time step. Different sp¯/sk¯ and sk/mse can be obtained with different σp, σT, σpar, and σs, that is, fσp,σT,σpar,σs=sp¯/sk¯ and gσp,σT,σpar,σs=sk/mse. Therefore, the optimal estimation of σp, σT, σpar, and σs can be realized by searching for the minimum value of the following function:(11)F=fσp,σT,σpar,σs−12−gσp,σT,σpar,σs−N+1/2N2.In the soil moisture assimilation, the temperature error (σT) is not very sensitive to DA performance. Hence, σT is set to be 1°C as referenced from the study of Chen et al. [12]. Besides, the predicted/simulated SM error σs is set to be 0.01 m3/m3 to avoid rapid changes of soil water content between continuous time steps [45, 46]. Finally, the precipitation (σp) and model parameter error (σpar) are determined by searching for the minimum F based on the streamflow measurements at the catchment outlet (Huaibin station) over 2002–2004 with an ensemble size of 200.


#### 3.3.3. ESA CCI SM Bias Correction

Remote sensing (RS) retrieval of SM often has systematic bias to the in situ observed and model-simulated soil moisture due to their large differences in spatial resolution and detection depth. The model-simulated SM can generally meet the water balance of the basin/region. In order to keep the basin’s water balance in DA, the systematic bias in RS SM needs to be corrected before DA [47]. In this study, the bias correction for the ESA CCI SM uses the cumulative distribution function (CDF) approach [48], where the probability of the RS SM and the simulated SM is assumed to be the same. The spatial matching between them uses the area-weighted average method to aggregate the simulated SM from HRUs to RS grids. Here, the rescaling is performed over the complete model validation period of 2002–2008, considering that CDF estimation typically requires a long record of observed and model simulated data [6].


#### 3.3.4. ESA CCI SM Error Estimation

Rational quantification on the uncertainty of RS SM is important for its optimal application. RS SM still has considerable uncertainty although its accuracy and reliability have been largely improved in recent years [49, 50]. In this study, referencing from previous researches [5, 51, 52], the error for ESA CCI SM is assumed to be an additive Gaussian distribution with the standard deviation (SD) of σR. Here, the estimation of σR is obtained from the equation referring to the study of Lievens et al. [30]:(12)σR=a0+b0sm_uncertainty+c0frc,where sm_uncertainty is an indicator of the data uncertainty for the ESA CCI SM [53, 54], which is not fully considered as the representativeness error (e.g., the error caused by vegetation or different layer depths). The representativeness error is accounted by the parameter a0, which represents the minimum retrieval error for ESA CCI SM. frc is the fraction of the ESA CCI SM grid cell covered by the forest. The calculation of frc is based on the land cover data collected from the Chinese Cold and Arid Regions Science Data Center (http://westdc.westgis.ac.cn/) (Figure 2). a0, b0, and c0 are given parameters, and b0, c0 ∈ (0, 1). Referencing from the study of Lievens et al. [30], a0, b0, and c0 are given as 0.02, 0.5, and 0.02, respectively, in this study. It should be noted that when the ESA CCI SM is high orderly rescaled (Section 3.3.3), the observation error parameter σR needs to be rescaled according to(13)σR∗=σsimσobsσR,where σR∗ is the standard deviation (SD) of the rescaled ESA CCI SM observation error, and σsim and σobs are the SD of the simulated SM error and the ESA CCI SM error, respectively.


### 3.4. Evaluation Metrics

The relative error (RE), the root mean square error (RMSE), the Nash–Sutcliffe coefficient of efficiency (NSE), and Pearson’s correlation coefficient (R) are used to measure the coincidence level of the simulated streamflow to the field observations. Meanwhile, the effectiveness criterion (EFF) [55] and the normalized error reduction index (NER) are used to directly assess the performance of soil moisture assimilation.RE describes the deviation rate (%) of the predicted streamflow to its field measurements. It can be expressed by(14)RE=∑i=1nQisim−Qiobs∑i=1nQiobs⋅100%,where n is the total time step and Qisim and Qiobs are the simulated and observed streamflow at time i.NSE is expressed by(15)NSE=1−∑i=1nQisim−Qiobs2∑i=1nQiobs−Q¯obs2,where Q¯obs indicates the mean value of the measured streamflow for the whole period. This NSE expression puts more importance on high flow. In order to give more weight to low flow, a modified version of the Nash–Sutcliffe coefficient of efficiency is adopted. It is actually a calculation of the NSE in a logarithmic form of the variable (NSElog):(16)NSElog=1−∑i=1nlogQisim−logQiobs2∑i=1nlogQiobs−logQ¯obs2.EFF reflects the data assimilation effects by comparing the sum of square error between the streamflow under assimilated and nonassimilated cases. It can be expressed as(17)EFF%=100⋅1−∑i=1nQiEnKF−Qiobs2∑i=1nQiEnOL−Qiobs2,where QiEnKF and QiEnOL are the predicted streamflow under assimilated and nonassimilated cases at time i.NER is expressed by the following [55]:(18)NER%=100⋅1.0−RMSEEnKFRMSEEnOL,where RMSEEnKF and RMSEEnOL are the root mean square errors of the variable in EnKF and EnOL (detailed in Section 4), respectively. The expression for RMSE and R can be found in the study of Liu et al. [56]. The EFF and NER > 0 means that the DA performance is effective. The larger EFF and NER values indicate the better performance of DA.


## 3.1. Soil and Water Assessment Tool (SWAT)

The SWAT is a physically based semidistributed watershed model, which has been widely used in rainfall-runoff modeling over recent years [35, 36]. In hydrological modeling, the catchment is firstly delineated into several subbasins according to its topography. Then, each subbasin is further divided into several hydrological response units (HRUs) based on the land use, soil, and slope. HRUs are basic calculation units for the land phase of the hydrologic cycle, on which the processes for surface runoff, lateral flow, and ground water are generated accompanied by evapotranspiration and soil water routing.Soil moisture lies in the center of the hydrologic cycle and makes different impacts on the above process. The initial profile soil water content influences surface runoff generation through the curve number in the SCS method [37]. After surface runoff generation, the water infiltrated to the soil profile is redistributed based on a storage routing technique with the soil water field capacity as the threshold. The water balance for each soil layer can be expressed as follows:(1)SWly′=SWly+Δwperc,ly−Qlat,ly−Ea,ly,where SWly and SWly′ are the soil water content (mm) at the start and end of the day, Δwperc,ly is the net percolation from the overlying layer (i.e., the layer ly+1), Qlat,ly is the lateral flow generated from the layer ly, and Ea,ly is the evapotranspiration drawn from the layer ly.The evapotranspiration from soil mainly includes two parts: soil evaporation and plant uptake/transpiration. In the SWAT, the potential evapotranspiration is firstly calculated using the Penman–Monteith equation [38]. Based on the potential evapotranspiration, the leaf area index, and the aboveground biomass and residue conditions, both the demand for transpiration/plant uptake and the demand for soil evaporation are determined. Then, the soil evaporation demand and the plant uptake demand for each soil layer are estimated using a depth distribution function. Finally, relying on the soil evaporation demand and the plant uptake demand with the available soil water as a constraint, the actual soil water evaporation and plant uptake are determined. In the processes mentioned above, the actual soil water extraction of a given layer is not allowed to be compensated by the extraction from other layers. However, the soil water deficiency can be made up by adjusting the soil compensation (esco) and plant compensation (epco) factors via changing the depth distribution of the soil evaporation demand and the plant water uptake demand. Besides, the calculation for the soil water percolation (wperc,ly) and lateral flow (Qlat,ly) is [39] omitted here.In the water routing phase, the SWAT adopts a storage feature to calculate the surface runoff and lateral flow generated from each HRU to the main channel, while it applies a linear reservoir similar technique to account for the ground water to the main channel. The channel water routing is performed using a variable storage routing method [40].


## 3.2. The Ensemble Kalman Filter (EnKF) for Soil Moisture Assimilation

The EnKF is a sequential DA approach evolved from the standard Kalman filter [41]. It is based on an ensemble of model states produced by adding the Monte Carlo noise to model forcing and states and/or parameters to approximate the model state error covariance matrix for the purpose of optimally merging the model predictions with observations.The state ensemble forecast at timet can be expressed as follows:(2)Xtf=FXt−1u,ut,δ+wtwt∼N0,σs2,where Xtf is the forecasted state ensemble at time t and Xt−1u is the updated state ensemble at t−1. In this study, it is constructed by the profile SM with up to four layers (Table 1) for all HRUs of the study basin (the HRUs delineation is detailed in Section 3.1). ut represents the model forcing inputs. In this study, it mainly includes the observed precipitation P and temperature T at each site. The precipitation error is assumed to be independent both in time and in space; that is, both the autocorrelation between time steps at each rainfall station and the error correlation among different stations are ignored. The perturbation (ηp) to precipitation is assumed to be a lognormal multiplicative distribution with mean 1 and covariance σp2 (3). The perturbation to temperature (ηT) is assumed to be an additive normal distribution with mean 0 and covariance σT2 (4). Besides, δ represents the model parameter with a perturbation (ηpar) of normal multiplicative distribution of mean 1 and covariance σpar2 (5). wt is the stochastic perturbation to the forecasted SM, being assumed to be an additive normal distribution with mean 0 and covariance σs2:(3)Pe=P⋅ηp,ηp∼lnN1,σp2,(4)Te=T+ηT,ηT∼N0,σT2,(5)δe=δ⋅ηpar,ηpar∼N1,σpar2.The state update for soil moisture can be obtained by(6)Xtu=Xtf+KtZt−HXtf,where Zt is the observation ensemble at time t. It is constructed by the RS SM for all grids covering the basin and being stochastically perturbed by an additive normal distribution with mean 0 and covariance σR2. H is the observation operator, being used to map the model states to the observations. It is constructed by the area proportions of HRUs in RS grids as the SWAT model-simulated SM is on the HRU level, while the observed SM is on RS grids. Kt is the Kalman gain, which is calculated based on the forecast and observation error covariance:(7)Kt=Pms,tPs,t+Rs,t−1,where Pms,t is the cross-error covariance between the predicted SM (Xtf) and the measurement prediction HXtf at time t, Ps,t is the error covariance of the measurement prediction at time t, and Rs,t is the error covariance of the RS SM at time t.


## 3.3. Implementation of the ESA CCI SM Assimilation in SWAT

To set up the ESA CCI SM data assimilation in the SWAT model, the SWAT model for the study catchment should be firstly built up. The model applicability is evaluated over the model calibration and validation processes. Then, the model error for the validated SWAT model is estimated based on the in situ streamflow observations at the Huaibin hydrologic station. Besides, the observation error estimation for the ESA CCI SM is implemented before the bias correction on the basis of the model-simulated soil moisture. A flow chart for implementing ESA CCI SM assimilation in SWAT is shown in Figure3.Figure 3
Flow chart for implementing ESA CCI SM assimilation in SWAT.


### 3.3.1. SWAT Calibration

The SWAT model for the upper Huai River basin is built up based on the meteorological forcing and land surface data (Section2.2). This catchment is partitioned into 37 subbasins (Figure 1) and 146 hydrological response units (HRUs). To simplify and improve the model calibration process, a sensitivity analysis of the model parameters to the hydrologic modeling is implemented using the Latin-hypercube and one-factor-at-a-time method [42, 43]. Thereafter, the model is calibrated and validated using the daily runoff records over 1992–1999 and 2002–2008 at the interior and outlet hydrologic stations (i.e., Dapoling, Changtaiguan, Zhuganfu, Xixian, Huangchuan, and Huaibin in Figure 1), respectively. The parameter optimization is achieved by a combination of the autocalibration using the Sequential Uncertainty Fitting (SUFI2) [44] with the Nash–Sutcliffe coefficient of efficiency (detailed in Section 3.4) as the objective criteria and the manual fine-running method.


### 3.3.2. Model Error Estimation

The determination of the model error is significant for the performance of DA as the model predictions and observations are merged based on the relative weight between the model and observation error (as in (6) and (7)). In this study, the model error is mainly contributed by the model input error for precipitation and temperature, the model parameter error for parameters sensitive to SM simulation (e.g., the available soil water capacity), and the model state error for simulated SM. The various errors mentioned above are characterized by additive/multiplicative normal/lognormal distribution specified in Section 3.2. The assumed distribution is only controlled by the standard deviation (SD), that is, σp, σT, σpar, and σs in (3), (4), (5), and (2), respectively. Therefore, the model error estimation is to quantify the model error parameters σp, σT, σpar, and σs.The quantification ofσp, σT, σpar, and σs is performed by analyzing the statistical characteristics of the simulated streamflow ensemble driven by model error perturbations on the basis of the observed streamflow at in situ sites. If the ensemble spread after perturbation is too large, the overfitting of observation exists in DA. Otherwise, the observed information cannot be fully utilized in DA. Therefore, the two ensemble verification measures (8) and (9) should be satisfied [32]. That is, if the ensemble spread sp is large enough, the temporal mean of the ensemble skill sk¯ should be similar to the temporal average of the ensemble spread sp¯:(8)sp¯sk¯≅1.And the observation should be indistinguishable from a member of the ensemble (N is the ensemble size):(9)skmse=N+12N,where(10)sp¯=1T∑kT1N∑i=1NQki−Q¯k2,sk¯=1T∑k=1TQ¯k−Qobs,k2,sk=1T∑k=1TQ¯k−Qobs,k,mse=1T∑k=1T1N∑i=1NQki−Qobs,k21/2,where Qki is the model simulated streamflow of the ensemble member i at time k, Q¯k is the ensemble mean of the model simulated streamflow, Qobs,k is the observed streamflow at time k, and T is the total time step. Different sp¯/sk¯ and sk/mse can be obtained with different σp, σT, σpar, and σs, that is, fσp,σT,σpar,σs=sp¯/sk¯ and gσp,σT,σpar,σs=sk/mse. Therefore, the optimal estimation of σp, σT, σpar, and σs can be realized by searching for the minimum value of the following function:(11)F=fσp,σT,σpar,σs−12−gσp,σT,σpar,σs−N+1/2N2.In the soil moisture assimilation, the temperature error (σT) is not very sensitive to DA performance. Hence, σT is set to be 1°C as referenced from the study of Chen et al. [12]. Besides, the predicted/simulated SM error σs is set to be 0.01 m3/m3 to avoid rapid changes of soil water content between continuous time steps [45, 46]. Finally, the precipitation (σp) and model parameter error (σpar) are determined by searching for the minimum F based on the streamflow measurements at the catchment outlet (Huaibin station) over 2002–2004 with an ensemble size of 200.


### 3.3.3. ESA CCI SM Bias Correction

Remote sensing (RS) retrieval of SM often has systematic bias to the in situ observed and model-simulated soil moisture due to their large differences in spatial resolution and detection depth. The model-simulated SM can generally meet the water balance of the basin/region. In order to keep the basin’s water balance in DA, the systematic bias in RS SM needs to be corrected before DA [47]. In this study, the bias correction for the ESA CCI SM uses the cumulative distribution function (CDF) approach [48], where the probability of the RS SM and the simulated SM is assumed to be the same. The spatial matching between them uses the area-weighted average method to aggregate the simulated SM from HRUs to RS grids. Here, the rescaling is performed over the complete model validation period of 2002–2008, considering that CDF estimation typically requires a long record of observed and model simulated data [6].


### 3.3.4. ESA CCI SM Error Estimation

Rational quantification on the uncertainty of RS SM is important for its optimal application. RS SM still has considerable uncertainty although its accuracy and reliability have been largely improved in recent years [49, 50]. In this study, referencing from previous researches [5, 51, 52], the error for ESA CCI SM is assumed to be an additive Gaussian distribution with the standard deviation (SD) of σR. Here, the estimation of σR is obtained from the equation referring to the study of Lievens et al. [30]:(12)σR=a0+b0sm_uncertainty+c0frc,where sm_uncertainty is an indicator of the data uncertainty for the ESA CCI SM [53, 54], which is not fully considered as the representativeness error (e.g., the error caused by vegetation or different layer depths). The representativeness error is accounted by the parameter a0, which represents the minimum retrieval error for ESA CCI SM. frc is the fraction of the ESA CCI SM grid cell covered by the forest. The calculation of frc is based on the land cover data collected from the Chinese Cold and Arid Regions Science Data Center (http://westdc.westgis.ac.cn/) (Figure 2). a0, b0, and c0 are given parameters, and b0, c0 ∈ (0, 1). Referencing from the study of Lievens et al. [30], a0, b0, and c0 are given as 0.02, 0.5, and 0.02, respectively, in this study. It should be noted that when the ESA CCI SM is high orderly rescaled (Section 3.3.3), the observation error parameter σR needs to be rescaled according to(13)σR∗=σsimσobsσR,where σR∗ is the standard deviation (SD) of the rescaled ESA CCI SM observation error, and σsim and σobs are the SD of the simulated SM error and the ESA CCI SM error, respectively.


## 3.3.1. SWAT Calibration

The SWAT model for the upper Huai River basin is built up based on the meteorological forcing and land surface data (Section2.2). This catchment is partitioned into 37 subbasins (Figure 1) and 146 hydrological response units (HRUs). To simplify and improve the model calibration process, a sensitivity analysis of the model parameters to the hydrologic modeling is implemented using the Latin-hypercube and one-factor-at-a-time method [42, 43]. Thereafter, the model is calibrated and validated using the daily runoff records over 1992–1999 and 2002–2008 at the interior and outlet hydrologic stations (i.e., Dapoling, Changtaiguan, Zhuganfu, Xixian, Huangchuan, and Huaibin in Figure 1), respectively. The parameter optimization is achieved by a combination of the autocalibration using the Sequential Uncertainty Fitting (SUFI2) [44] with the Nash–Sutcliffe coefficient of efficiency (detailed in Section 3.4) as the objective criteria and the manual fine-running method.


## 3.3.2. Model Error Estimation

The determination of the model error is significant for the performance of DA as the model predictions and observations are merged based on the relative weight between the model and observation error (as in (6) and (7)). In this study, the model error is mainly contributed by the model input error for precipitation and temperature, the model parameter error for parameters sensitive to SM simulation (e.g., the available soil water capacity), and the model state error for simulated SM. The various errors mentioned above are characterized by additive/multiplicative normal/lognormal distribution specified in Section 3.2. The assumed distribution is only controlled by the standard deviation (SD), that is, σp, σT, σpar, and σs in (3), (4), (5), and (2), respectively. Therefore, the model error estimation is to quantify the model error parameters σp, σT, σpar, and σs.The quantification ofσp, σT, σpar, and σs is performed by analyzing the statistical characteristics of the simulated streamflow ensemble driven by model error perturbations on the basis of the observed streamflow at in situ sites. If the ensemble spread after perturbation is too large, the overfitting of observation exists in DA. Otherwise, the observed information cannot be fully utilized in DA. Therefore, the two ensemble verification measures (8) and (9) should be satisfied [32]. That is, if the ensemble spread sp is large enough, the temporal mean of the ensemble skill sk¯ should be similar to the temporal average of the ensemble spread sp¯:(8)sp¯sk¯≅1.And the observation should be indistinguishable from a member of the ensemble (N is the ensemble size):(9)skmse=N+12N,where(10)sp¯=1T∑kT1N∑i=1NQki−Q¯k2,sk¯=1T∑k=1TQ¯k−Qobs,k2,sk=1T∑k=1TQ¯k−Qobs,k,mse=1T∑k=1T1N∑i=1NQki−Qobs,k21/2,where Qki is the model simulated streamflow of the ensemble member i at time k, Q¯k is the ensemble mean of the model simulated streamflow, Qobs,k is the observed streamflow at time k, and T is the total time step. Different sp¯/sk¯ and sk/mse can be obtained with different σp, σT, σpar, and σs, that is, fσp,σT,σpar,σs=sp¯/sk¯ and gσp,σT,σpar,σs=sk/mse. Therefore, the optimal estimation of σp, σT, σpar, and σs can be realized by searching for the minimum value of the following function:(11)F=fσp,σT,σpar,σs−12−gσp,σT,σpar,σs−N+1/2N2.In the soil moisture assimilation, the temperature error (σT) is not very sensitive to DA performance. Hence, σT is set to be 1°C as referenced from the study of Chen et al. [12]. Besides, the predicted/simulated SM error σs is set to be 0.01 m3/m3 to avoid rapid changes of soil water content between continuous time steps [45, 46]. Finally, the precipitation (σp) and model parameter error (σpar) are determined by searching for the minimum F based on the streamflow measurements at the catchment outlet (Huaibin station) over 2002–2004 with an ensemble size of 200.


## 3.3.3. ESA CCI SM Bias Correction

Remote sensing (RS) retrieval of SM often has systematic bias to the in situ observed and model-simulated soil moisture due to their large differences in spatial resolution and detection depth. The model-simulated SM can generally meet the water balance of the basin/region. In order to keep the basin’s water balance in DA, the systematic bias in RS SM needs to be corrected before DA [47]. In this study, the bias correction for the ESA CCI SM uses the cumulative distribution function (CDF) approach [48], where the probability of the RS SM and the simulated SM is assumed to be the same. The spatial matching between them uses the area-weighted average method to aggregate the simulated SM from HRUs to RS grids. Here, the rescaling is performed over the complete model validation period of 2002–2008, considering that CDF estimation typically requires a long record of observed and model simulated data [6].


## 3.3.4. ESA CCI SM Error Estimation

Rational quantification on the uncertainty of RS SM is important for its optimal application. RS SM still has considerable uncertainty although its accuracy and reliability have been largely improved in recent years [49, 50]. In this study, referencing from previous researches [5, 51, 52], the error for ESA CCI SM is assumed to be an additive Gaussian distribution with the standard deviation (SD) of σR. Here, the estimation of σR is obtained from the equation referring to the study of Lievens et al. [30]:(12)σR=a0+b0sm_uncertainty+c0frc,where sm_uncertainty is an indicator of the data uncertainty for the ESA CCI SM [53, 54], which is not fully considered as the representativeness error (e.g., the error caused by vegetation or different layer depths). The representativeness error is accounted by the parameter a0, which represents the minimum retrieval error for ESA CCI SM. frc is the fraction of the ESA CCI SM grid cell covered by the forest. The calculation of frc is based on the land cover data collected from the Chinese Cold and Arid Regions Science Data Center (http://westdc.westgis.ac.cn/) (Figure 2). a0, b0, and c0 are given parameters, and b0, c0 ∈ (0, 1). Referencing from the study of Lievens et al. [30], a0, b0, and c0 are given as 0.02, 0.5, and 0.02, respectively, in this study. It should be noted that when the ESA CCI SM is high orderly rescaled (Section 3.3.3), the observation error parameter σR needs to be rescaled according to(13)σR∗=σsimσobsσR,where σR∗ is the standard deviation (SD) of the rescaled ESA CCI SM observation error, and σsim and σobs are the SD of the simulated SM error and the ESA CCI SM error, respectively.


## 3.4. Evaluation Metrics

The relative error (RE), the root mean square error (RMSE), the Nash–Sutcliffe coefficient of efficiency (NSE), and Pearson’s correlation coefficient (R) are used to measure the coincidence level of the simulated streamflow to the field observations. Meanwhile, the effectiveness criterion (EFF) [55] and the normalized error reduction index (NER) are used to directly assess the performance of soil moisture assimilation.RE describes the deviation rate (%) of the predicted streamflow to its field measurements. It can be expressed by(14)RE=∑i=1nQisim−Qiobs∑i=1nQiobs⋅100%,where n is the total time step and Qisim and Qiobs are the simulated and observed streamflow at time i.NSE is expressed by(15)NSE=1−∑i=1nQisim−Qiobs2∑i=1nQiobs−Q¯obs2,where Q¯obs indicates the mean value of the measured streamflow for the whole period. This NSE expression puts more importance on high flow. In order to give more weight to low flow, a modified version of the Nash–Sutcliffe coefficient of efficiency is adopted. It is actually a calculation of the NSE in a logarithmic form of the variable (NSElog):(16)NSElog=1−∑i=1nlogQisim−logQiobs2∑i=1nlogQiobs−logQ¯obs2.EFF reflects the data assimilation effects by comparing the sum of square error between the streamflow under assimilated and nonassimilated cases. It can be expressed as(17)EFF%=100⋅1−∑i=1nQiEnKF−Qiobs2∑i=1nQiEnOL−Qiobs2,where QiEnKF and QiEnOL are the predicted streamflow under assimilated and nonassimilated cases at time i.NER is expressed by the following [55]:(18)NER%=100⋅1.0−RMSEEnKFRMSEEnOL,where RMSEEnKF and RMSEEnOL are the root mean square errors of the variable in EnKF and EnOL (detailed in Section 4), respectively. The expression for RMSE and R can be found in the study of Liu et al. [56]. The EFF and NER > 0 means that the DA performance is effective. The larger EFF and NER values indicate the better performance of DA.


## 4. Results

The efficiency of ESA CCI SM assimilation is highly dependent on the quality of model calibration, model, and observation error estimation. Hence, the ESA CCI SM assimilation effects on streamflow simulation accompanied by the results for model calibration and validation and for model and observation error estimation are analyzed. To illustrate the efficiency of DA, the ensemble open-loop (EnOL) cases and the EnKF cases are compared. The EnOL is an ensemble running of the SWAT model with perturbations on model inputs, model parameters, and model states without the integration of observed SM, while the EnKF is an ensemble running of the SWAT model with the same perturbation to EnOL, but with the integration of ESA CCI SM during the model propagation process.


### 4.1. Model Calibration and Validation

Table2 presents the SWAT model parameters being calibrated, which are obtained from the parameter sensitivity analysis detailed in Section 3.3.1. Figure 4 plots the simulated and observed daily series of runoff at the catchment outlet during the calibration (1992–1999) and validation (2002–2008) periods. The hydrograph of the simulated streamflow is highly consistent with that of the observed streamflow for both the calibration and validation stages, although slight underestimation exists in flood peak modeling over some periods. The statistics (Table 3) for the simulated streamflow at the catchment outlet (Huaibin) suggest that it agrees well with the measured runoff as RE < 5%, NSE > 0.8, and R > 0.9. In addition, the statistics for the other five hydrological sites (Dapoling, Changtaiguan, Zhuganfu, Xixian, and Huangchuan) also indicate that the SWAT model has fairly good applicability in the upper Huai River basin. In the calibration stage, for all six stations, RE < 15%, NSE falls between 0.65∼0.81, and R > 0.83. In the validation state, RE < 15%, NSE falls between 0.69∼0.9 (except for Huangchuan), and R > 0.83. The unexpected NSEs at Huangchuan are caused by the serious disturbances of the human activities (principally the reservoir impacts) on runoff over the years 2004, 2006, and 2008.Table 2
SWAT model parameters being calibrated.
Sensitivity sequence
Parameters
Description
1
CN2
SCS curve number for moisture conditon II
2
surlag
Surface runoff lag coefficient
3
α
gw
Baseflow recession constant
4
K
ch
Effective hydraulic conductivity
5
esco
Soil evaporation compensation constant
6
aqshthr,q
Threshold depth of water in the shallow aquifer required for the return flow to occur
7
SOL_AWC
Available soil water capacity
8
N
Manning’sn value for the main channel
9
canmx
Maximum amount of water that can be trapped in the canopy when that canopy is fully developed
10
epco
Plant uptake compensation factor
11
β
rev
Revap coefficient
12
δ
gw
Delay time for aquifer rechargeFigure 4
The observed (Obs) and simulated (Sim) daily runoff at the outlet (Huaibin) of the upper Huai River basin over the calibration (1992–1999) and validation (2002–2008) periods.Table 3
Statistical comparison of the observed and simulated daily runoff over the calibration (1992–1999) and validation (2002–2008) periods.
Hydrologic station
Calibration (1992–1999)
Validation (2002–2008)
RE (%)
RMSE (m3/s)
NSE
R
RE (%)
RMSE (m3/s)
NSE
R
Dapoling
12.97
28.06
0.66
0.84
15.14
35.49
0.81
0.91
Changtaiguan
−0.74
47.59
0.69
0.85
6.21
63.27
0.78
0.89
Zhuganfu
−7.74
38.53
0.64
0.83
−13.29
51.29
0.69
0.83
Xixian
−1.65
115.84
0.75
0.88
2.24
173.05
0.75
0.88
Huangchuan
0.13
39.61
0.67
0.83
9.89
68.32
0.42
0.73
Huaibin
3.56
133.74
0.81
0.91
1.74
154.13
0.9
0.95


### 4.2. Model Error Estimation for Precipitation and Model Parameter

Figure5 shows the objective function F (11) with varying standard deviation (SD) of the lognormal multiplicative perturbation on precipitation (σp) from 0.05 to 0.5 along with the varying SD of the normal multiplicative perturbation on parameter from 0.1 to 0.5. The objective function F reaches its minimum value when σp approximates 0.35, which suggests that, in this case, the simulated outlet streamflow is the best matching to the observed streamflow from its ensemble statistics. However, it can be seen that σpar is not that sensitive to the objective function, and the allotropism or nonuniqueness issue exits in its optimal parameter estimation. Considering the good performance of the SWAT model in the study basin (Section 4.1), the small values of σpar (<0.25) are more credible. Besides, in consideration of the robustness of the EnKF method [45], 0.25 for σpar is adopted. Note that, in this estimation of σp and σpar, the observed streamflow at the subbasin outlet is regarded as the truth, that is, the observation error is ignored. In this case, σp and σpar are likely to be slightly overestimated.Figure 5
Contour plot of the objective functionF (11) with respect to the standard error deviation for precipitation (σp) and for model parameter (σpar).


### 4.3. ESA CCI SM Error Estimation

Figure6 shows the standard deviation (SD) of the observation error σR∗ (13) for the ESA CCI SM at each grid (34 grids in total) within the catchment. The location for 1–34 grids is present in Figure 7. In general, σR∗ falls between 0.03 and 0.05 m3/m3 for all grids, which is consistent with the accuracy of the ESA CCI SM on average (0.04∼0.05 m3/m3) [57]. At each grid, σR∗ presents certain ranges (∼0.01 m3/m3), which is related to the soil moisture dynamic with time changes over 2003–2006. This also indicates the necessity for considering the temporal characteristics of the observation error for RS SM. In addition, σR∗ shows a considerable difference among different grids. The high σR∗ mainly appears on grids with the dense forest coverage, for example, the grids 10, 11, 25, 26, 33, and 34, in particular for the grid 33. The dense forest obscures the emitted radiance of the soil surface, which results in large uncertainty to the surface SM retrieval. Besides, major dense forests are distributed over the catchment with high altitudes (Figure 1), where the complex topography also impedes the accuracy and reliability of remote sensing for SM [58].Figure 6
Standard deviation (σR∗) of the observation error for the rescaled ESA CCI SM at each grid within the basin. The upper limb, the lower limb, and the red line in the box plot represent the upper quartile (h1), the lower quartile (h2), and the median value of σR∗ over 2003–2006. The dotted line stretched from the box is the range between h1 − 1.5 (h2 − h1) and h2 + 1.5 (h2 − h1).Figure 7
Forest coverage for the remote sensing grids of the ESA CCI SM over the upper Huai River basin. 1–34 is the grid coding for the ESA CCI SM.


### 4.4. ESA CCI SM Assimilation on Streamflow Simulation

Table4 statistically compares the model simulated streamflow with (EnKF) and without (EnOL) ESA CCI SM assimilation at the six hydrologic sites in the upper Huai River basin except for Huangchuan (Figure 1). The reason for Huangchuan not being taken into account is that it has data quality issue over 2004 and 2006 caused by severe human activities. Table 4 shows that the RE and RMSE are decreased and the NSE, NSElog, and R are increased at the five gauges except for Zhuganfu due to ESA CCI soil moisture assimilation. The improvement is more significant in terms of the NSElog as its increase rate is greater than NSE and R, which indicates that the RS soil moisture assimilation is more effective for low flows than high flows. Besides, EFF/NER > 0 for four sites, in particular for Dapoling, Xixian, and Huaibin (where NER > 5% and EFF > 10%), which suggests the good performance of the assimilation. The noneffective performance of ESA CCI SM assimilation on runoff simulation of Changtaiguan and Zhuganfu is probably related to their large proportions of the dense forest and complex topography coverage upstream (Figures 1 and 7). Both dense forest coverage and complex topographical conditions reduce the data quality of RS SM retrievals, thus impeding its performance in DA.Table 4
Statistical comparison of the estimated streamflow in EnOL and EnKF cases based on the observed streamflow over 2003–2006.
Hydrologic station
EnOL
EnKF
EFF (%)
NER (%)
RE (%)
RMSE (m3/s)
NSE
NSElog
R
RE (%)
RMSE (m3/s)
NSE
NSElog
R
Dapoling
15.15
35.28
0.82
0.05
0.92
6.65
33
0.84
0.1
0.92
12.51
6.47
Changtaiguan
7.43
61.71
0.78
0.62
0.89
0.09
60.89
0.79
0.68
0.89
2.65
1.33
Zhuganfu
−11.4
43.87
0.67
0.04
0.83
−18.9
44.45
0.66
0.16
0.82
−2.67
−1.3
Xixian
4.1
165
0.74
0.44
0.88
−3.61
156.37
0.77
0.57
0.88
10.19
5.23
Huaibin
5.74
149.9
0.88
0.6
0.94
0.02
140.2
0.9
0.71
0.95
12.47
6.44Figure8 compares the daily series of the model simulated streamflow at the catchment outlet during 2003–2006 with (EnKF) and without (EnOL) ESA CCI SM assimilation on the basis of the observed runoff (Obs). It can be seen that ESA CCI SM assimilation improves the streamflow modeling over low-flow periods. The predicted runoff with ESA CCI SM assimilation (the red line) is closer to the observed runoff (the black line), and the NSElog increases from 0.6 to 0.71. However, the impact of soil moisture assimilation on streamflow modeling over very high-flow/large-flood periods presents certain uncertainty. For instance, data assimilation improves the streamflow simulation over the periods August 2, 2004–August 12, 2004, and July 3, 2006–July 23, 2006, while it deteriorates the streamflow modeling over the period June 29, 2003–July 9, 2003. These results are consistent with those of the previous researches, for example, the study of Alvarez-Garreton et al. [5] and the study of Massari et al. [32]. The uncertain performance of soil moisture assimilation on large-flood simulation mainly lies in the relatively low dependence of runoff generation on antecedent soil moisture because during large-flood periods, the soil moisture is nearly saturated and the runoff is largely controlled by precipitation inputs.Figure 8
The model-simulated streamflow at the catchment outlet (Huaibin) with (EnKF) and without (EnOL) the ESA CCI SM data assimilation over 2003–2006. The upper plot is for the whole period with a base 10 logarithm coordinate. The four plots below are for the flood season from 2003 to 2006. The red and blue lines are the EnKF and EnOL ensemble members, and the red and blue bold lines represent the mean value of EnKF and EnOL. Obs is the abbreviation of observation.
(a)
(b)


## 4.1. Model Calibration and Validation

Table2 presents the SWAT model parameters being calibrated, which are obtained from the parameter sensitivity analysis detailed in Section 3.3.1. Figure 4 plots the simulated and observed daily series of runoff at the catchment outlet during the calibration (1992–1999) and validation (2002–2008) periods. The hydrograph of the simulated streamflow is highly consistent with that of the observed streamflow for both the calibration and validation stages, although slight underestimation exists in flood peak modeling over some periods. The statistics (Table 3) for the simulated streamflow at the catchment outlet (Huaibin) suggest that it agrees well with the measured runoff as RE < 5%, NSE > 0.8, and R > 0.9. In addition, the statistics for the other five hydrological sites (Dapoling, Changtaiguan, Zhuganfu, Xixian, and Huangchuan) also indicate that the SWAT model has fairly good applicability in the upper Huai River basin. In the calibration stage, for all six stations, RE < 15%, NSE falls between 0.65∼0.81, and R > 0.83. In the validation state, RE < 15%, NSE falls between 0.69∼0.9 (except for Huangchuan), and R > 0.83. The unexpected NSEs at Huangchuan are caused by the serious disturbances of the human activities (principally the reservoir impacts) on runoff over the years 2004, 2006, and 2008.Table 2
SWAT model parameters being calibrated.
Sensitivity sequence
Parameters
Description
1
CN2
SCS curve number for moisture conditon II
2
surlag
Surface runoff lag coefficient
3
α
gw
Baseflow recession constant
4
K
ch
Effective hydraulic conductivity
5
esco
Soil evaporation compensation constant
6
aqshthr,q
Threshold depth of water in the shallow aquifer required for the return flow to occur
7
SOL_AWC
Available soil water capacity
8
N
Manning’sn value for the main channel
9
canmx
Maximum amount of water that can be trapped in the canopy when that canopy is fully developed
10
epco
Plant uptake compensation factor
11
β
rev
Revap coefficient
12
δ
gw
Delay time for aquifer rechargeFigure 4
The observed (Obs) and simulated (Sim) daily runoff at the outlet (Huaibin) of the upper Huai River basin over the calibration (1992–1999) and validation (2002–2008) periods.Table 3
Statistical comparison of the observed and simulated daily runoff over the calibration (1992–1999) and validation (2002–2008) periods.
Hydrologic station
Calibration (1992–1999)
Validation (2002–2008)
RE (%)
RMSE (m3/s)
NSE
R
RE (%)
RMSE (m3/s)
NSE
R
Dapoling
12.97
28.06
0.66
0.84
15.14
35.49
0.81
0.91
Changtaiguan
−0.74
47.59
0.69
0.85
6.21
63.27
0.78
0.89
Zhuganfu
−7.74
38.53
0.64
0.83
−13.29
51.29
0.69
0.83
Xixian
−1.65
115.84
0.75
0.88
2.24
173.05
0.75
0.88
Huangchuan
0.13
39.61
0.67
0.83
9.89
68.32
0.42
0.73
Huaibin
3.56
133.74
0.81
0.91
1.74
154.13
0.9
0.95


## 4.2. Model Error Estimation for Precipitation and Model Parameter

Figure5 shows the objective function F (11) with varying standard deviation (SD) of the lognormal multiplicative perturbation on precipitation (σp) from 0.05 to 0.5 along with the varying SD of the normal multiplicative perturbation on parameter from 0.1 to 0.5. The objective function F reaches its minimum value when σp approximates 0.35, which suggests that, in this case, the simulated outlet streamflow is the best matching to the observed streamflow from its ensemble statistics. However, it can be seen that σpar is not that sensitive to the objective function, and the allotropism or nonuniqueness issue exits in its optimal parameter estimation. Considering the good performance of the SWAT model in the study basin (Section 4.1), the small values of σpar (<0.25) are more credible. Besides, in consideration of the robustness of the EnKF method [45], 0.25 for σpar is adopted. Note that, in this estimation of σp and σpar, the observed streamflow at the subbasin outlet is regarded as the truth, that is, the observation error is ignored. In this case, σp and σpar are likely to be slightly overestimated.Figure 5
Contour plot of the objective functionF (11) with respect to the standard error deviation for precipitation (σp) and for model parameter (σpar).


## 4.3. ESA CCI SM Error Estimation

Figure6 shows the standard deviation (SD) of the observation error σR∗ (13) for the ESA CCI SM at each grid (34 grids in total) within the catchment. The location for 1–34 grids is present in Figure 7. In general, σR∗ falls between 0.03 and 0.05 m3/m3 for all grids, which is consistent with the accuracy of the ESA CCI SM on average (0.04∼0.05 m3/m3) [57]. At each grid, σR∗ presents certain ranges (∼0.01 m3/m3), which is related to the soil moisture dynamic with time changes over 2003–2006. This also indicates the necessity for considering the temporal characteristics of the observation error for RS SM. In addition, σR∗ shows a considerable difference among different grids. The high σR∗ mainly appears on grids with the dense forest coverage, for example, the grids 10, 11, 25, 26, 33, and 34, in particular for the grid 33. The dense forest obscures the emitted radiance of the soil surface, which results in large uncertainty to the surface SM retrieval. Besides, major dense forests are distributed over the catchment with high altitudes (Figure 1), where the complex topography also impedes the accuracy and reliability of remote sensing for SM [58].Figure 6
Standard deviation (σR∗) of the observation error for the rescaled ESA CCI SM at each grid within the basin. The upper limb, the lower limb, and the red line in the box plot represent the upper quartile (h1), the lower quartile (h2), and the median value of σR∗ over 2003–2006. The dotted line stretched from the box is the range between h1 − 1.5 (h2 − h1) and h2 + 1.5 (h2 − h1).Figure 7
Forest coverage for the remote sensing grids of the ESA CCI SM over the upper Huai River basin. 1–34 is the grid coding for the ESA CCI SM.


## 4.4. ESA CCI SM Assimilation on Streamflow Simulation

Table4 statistically compares the model simulated streamflow with (EnKF) and without (EnOL) ESA CCI SM assimilation at the six hydrologic sites in the upper Huai River basin except for Huangchuan (Figure 1). The reason for Huangchuan not being taken into account is that it has data quality issue over 2004 and 2006 caused by severe human activities. Table 4 shows that the RE and RMSE are decreased and the NSE, NSElog, and R are increased at the five gauges except for Zhuganfu due to ESA CCI soil moisture assimilation. The improvement is more significant in terms of the NSElog as its increase rate is greater than NSE and R, which indicates that the RS soil moisture assimilation is more effective for low flows than high flows. Besides, EFF/NER > 0 for four sites, in particular for Dapoling, Xixian, and Huaibin (where NER > 5% and EFF > 10%), which suggests the good performance of the assimilation. The noneffective performance of ESA CCI SM assimilation on runoff simulation of Changtaiguan and Zhuganfu is probably related to their large proportions of the dense forest and complex topography coverage upstream (Figures 1 and 7). Both dense forest coverage and complex topographical conditions reduce the data quality of RS SM retrievals, thus impeding its performance in DA.Table 4
Statistical comparison of the estimated streamflow in EnOL and EnKF cases based on the observed streamflow over 2003–2006.
Hydrologic station
EnOL
EnKF
EFF (%)
NER (%)
RE (%)
RMSE (m3/s)
NSE
NSElog
R
RE (%)
RMSE (m3/s)
NSE
NSElog
R
Dapoling
15.15
35.28
0.82
0.05
0.92
6.65
33
0.84
0.1
0.92
12.51
6.47
Changtaiguan
7.43
61.71
0.78
0.62
0.89
0.09
60.89
0.79
0.68
0.89
2.65
1.33
Zhuganfu
−11.4
43.87
0.67
0.04
0.83
−18.9
44.45
0.66
0.16
0.82
−2.67
−1.3
Xixian
4.1
165
0.74
0.44
0.88
−3.61
156.37
0.77
0.57
0.88
10.19
5.23
Huaibin
5.74
149.9
0.88
0.6
0.94
0.02
140.2
0.9
0.71
0.95
12.47
6.44Figure8 compares the daily series of the model simulated streamflow at the catchment outlet during 2003–2006 with (EnKF) and without (EnOL) ESA CCI SM assimilation on the basis of the observed runoff (Obs). It can be seen that ESA CCI SM assimilation improves the streamflow modeling over low-flow periods. The predicted runoff with ESA CCI SM assimilation (the red line) is closer to the observed runoff (the black line), and the NSElog increases from 0.6 to 0.71. However, the impact of soil moisture assimilation on streamflow modeling over very high-flow/large-flood periods presents certain uncertainty. For instance, data assimilation improves the streamflow simulation over the periods August 2, 2004–August 12, 2004, and July 3, 2006–July 23, 2006, while it deteriorates the streamflow modeling over the period June 29, 2003–July 9, 2003. These results are consistent with those of the previous researches, for example, the study of Alvarez-Garreton et al. [5] and the study of Massari et al. [32]. The uncertain performance of soil moisture assimilation on large-flood simulation mainly lies in the relatively low dependence of runoff generation on antecedent soil moisture because during large-flood periods, the soil moisture is nearly saturated and the runoff is largely controlled by precipitation inputs.Figure 8
The model-simulated streamflow at the catchment outlet (Huaibin) with (EnKF) and without (EnOL) the ESA CCI SM data assimilation over 2003–2006. The upper plot is for the whole period with a base 10 logarithm coordinate. The four plots below are for the flood season from 2003 to 2006. The red and blue lines are the EnKF and EnOL ensemble members, and the red and blue bold lines represent the mean value of EnKF and EnOL. Obs is the abbreviation of observation.
(a)
(b)


## 5. Discussion

In general, our results indicate that the ESA CCI soil moisture assimilation in SWAT performs well in runoff modeling of the whole basin. Streamflow improvements over five in situ sites (except for Zhuganfu) are shown after proper configurations of the model and observation error. However, the improvements are not significant, which can be attributed to the following factors: First, the model error is estimated based on analyzing the ensemble characteristics of the streamflow simulations driven by the model error perturbations in reference to the ground-based runoff observations, during which the observation error for the streamflow is ignored. It might lead to an overestimated model error. In observation error estimation for satellite soil moisture, subjectivity does exist in parameter assignment of the estimation equation although the temporal and spatial variability has been taken into account. These two factors are likely to deteriorate the model and observation error estimation, which eventually degrade the DA performance. Second, the runoff improvements are obtained by updating the profile SM using the satellite SM products, which highly relies on the physical vertical coupling-based model. SWAT soil layers have limited vertical coupling [12, 18] as it does not allow actual soil water compensation from other soil layers in the storage routing technique, and the soil water deficiency is only made up by adjusting the depth distribution of soil evaporation demand. The exponential filter [59] used to derive the profile SM indicator from the surface SM observations is a common solution to the inconsistency of the shallow (surface) RS detection and the runoff root zone control mechanism [5, 32]. However, this approach is more applicable to the hydrological model with a single soil layer setup. For the multilayer setup model (e.g., SWAT), a more promising approach to the physical coupling issue would be adopting Richard’s equation because it is more representative of the real-world water movement of soil water. Finally, the runoff improvements show large discrepancies over different hydrological sites with different geographical locations, which suggests that the land surface conditions considerably influence the DA performance (similar to the results from the study of Massari et al. [32]). The dense forests and complex topographical conditions reduce the data quality of microwave soil moisture retrievals, thus deteriorating the efficiency of satellite soil moisture assimilation.


## 6. Conclusions

The ESA CCI soil moisture assimilation in the SWAT model using the ensemble Kalman filter (EnKF) with the objective to improve its rainfall-runoff simulation is undertaken in the upper Huai River basin in China. In the assimilation framework, the bias correction for the ESA CCI SM is based on the model-simulated SM using the cumulative distribution function method. The model error is estimated by analyzing the statistical character of the simulated streamflow ensemble under various model error perturbations based on the in situ observed runoff at the catchment outlet. In observation error assessment for the ESA CCI SM, both the spatial heterogeneity and the temporal variability are considered. The observation error is obtained from a linear combination of the minimum retrieval error, the uncertainty indicator in retrieval, and the forest proportions on SM grids. Besides, the SWAT model applicability for the study catchment is assessed before DA as it has significant impacts on the performance of satellite SM assimilation.The SWAT model has good applicability to the study catchment as the model-simulated daily runoff series are highly consistent with those of the in situ measurements, and the evaluation statistics for all six hydrologic stations are satisfying. In general, the ESA CCI SM assimilation improves the streamflow modeling of the study basin. The DA is more effective for the improvement of low-flow simulation, while for very high-flow/large-flood modeling, the DA performance presents uncertainty. Besides, the DA efficiency is likely to be deteriorated by the dense forest coverage and the complex topographical conditions as it shows large discrepancy over the stations with large and small proportions of the mountainous region at the upstream. Overall, the coarse-scale ESA CCI SM assimilation could improve the streamflow modeling of a physically based semidistributed model, especially for low flow. This study provides an encouragement for the application of the ESA CCI SM in water management over dry seasons or low-flow periods.



---

*Source: 7301314-2018-06-07.xml*