A Comprehensive Approach to Develop a Hydrological Model for the Simulation of All the Important Hydrological Components: The Case of the Three-River Headwater Region, China

Abstract

The objective of the study was to configure the Hydrological Modeling System (HEC-HMS) in such a way that it could simulate all-important hydrological components (e.g., streamflow, soil moisture, snowmelt water, terrestrial water storage, baseflow, surface flow, and evapotranspiration) in the Three-River Headwater Region. However, the problem we faced was unsatisfactory simulations of these hydrological components, except streamflow. The main reason we found was the auto-calibration method of HEC-HMS because it generated irrational parameters, especially with the inclusion of Temperature Index Method and Soil Moisture Accounting (an advanced and complex loss method). Similar problems have been reported by different previous studies. To overcome these problems, we designed a comprehensive approach to estimate initial parameters and to calibrate the model manually in such a way that the model could simulate all the important hydrological components satisfactorily.

1. Introduction

Understanding how water is used and cycled through a watershed provides the foundation for understanding and describing the interactions of landscapes and water. The most basic tool for understanding these interactions is the water cycle [1]. The water cycle is one of the most important cycles (e.g., carbon cycle, the nitrogen cycle, and other biogeochemical cycles) of Earth’s climate system [2,3] and comprises different hydrological components such as precipitation, evaporation, transpiration, soil water, and streamflow [1]. Therefore, for precise water resources management and planning and sustainable development for a region, it is essential to study all the important components (e.g., precipitation, evaporation, transpiration, soil water, baseflow, streamflow, snowmelt, and terrestrial water storage) of the water cycle [4,5]. However, the observed data of these components, except precipitation and streamflow, are hardly available because they are difficult to measure for longer periods and larger watersheds [6,7]. For the past two decades, hydrological models have become the most powerful tools to study these hydrological components [8,9]. However, their performance is mainly dependent on the calibration, which is the most important and central process in the development of a model [10,11]. Generally, calibration methods are divided into two categories: manual and automatic [12]. In manual calibration, parameters are adjusted manually to achieve an objective function (e.g., the minimization of root mean square error). However, this is time-consuming, subjective, and tedious, especially when the number of parameters increases, and requires expertise in the field of hydrology and modeling [9,12,13].

Recently, auto-calibration has become mainstream for solving the issues of optimizing the model’s parameters [9,14]. These computer-aided methods automatically search a set of optimal parameters to some predefined criteria (e.g., percent error in peak) [15]. Additionally, these methods are considered objective, speedy, and relatively easy to implement [9,16]. Although intensive advancement has been made in auto-calibration, most applications have been restricted to research and not operational modeling [17], and still many problems have been reported by researchers, especially related to the reliability of calibrated parameters [18,19,20,21]. Similar problems were faced in this study when we tried to calibrate the Hydrological Engineering Center–Hydrological Modeling System (HEC-HMS) automatically using built-in optimizing techniques in the Three-River Headwater Region (TRHR). According to Bashar and Zaki [22] and Dariane, Javadianzadeh and James [14], the optimization methods (i.e., Nelder–Mead and Univariate Gradient) used in HEC-HMS are either inefficient or do not work rationally, especially with the inclusion of Soil Moisture Accounting (SMA) and Temperature Index Method (TIM). Therefore, many studies, such as, Bhuiyan, McNairn, Powers and Merzouki [18], Chea and Oeurng [23], Fleming and Neary [24], Gebre [19], Gyawali and David [25], Koch and Bene [20], and Neary, et al. [26] suggested using both the automated and manual (preferably) methods during calibration of HEC-HMS when using SMA and TIM, because the parameters obtained by the automated methods are mostly irrational, which yields misleading results [21,22].

As mentioned above, many studies suggest using manual calibration; however, we could not find any comprehensive approach that can be used to configure HEC-HMS in such a way that it can be used to generate all the important hydrological components satisfactorily. Moreover, these studies focused only on the simulation of streamflow. Therefore, in this study, we designed a comprehensive approach to configure HEC-HMS for the simulation of all the important hydrological components (e.g., streamflow, baseflow, soil moisture, snowmelt water, evapotranspiration, and terrestrial water storage) satisfactorily. The focus of the approach was to provide a way to estimate some logical initial values for the process parameters, explore the sensitive parameters, and determine the relationships between parameters and hydrograph components, e.g., which parameter affects which hydrograph component (e.g., low flow, high flow, rising limb, and falling limb). This not only can reduce the time of calibration but also enhance the understanding of the hydrological processes and provide some rational parameters and, consequently, some reasonable results. Additionally, this approach can be applied to similar kinds of hydrological models (e.g., SWAT, TOPMODEL, and HBV) when it is necessary to simulate all the important hydrological components of the hydrological cycle of a region.

2. Materials and Methods

2.1. Study Area

The Three-River Headwater Region (TRHR) is the source region of the Lancang, Yellow, and Yangtze Rivers located in the Qinghai Tibetan Plateau (QTP), China (Figure 1). The source region of the Yangtze River (SYAR) above the Zhimenda streamflow gauge, the source region of the Yellow River (SYER) above Tangnaihai, and the source region of the Lancang River (SLR) above Xiangda contribute areas of 54% (158,027 km²), 40% (117,755 km²), and 6% (16,919 km²), respectively, to the TRHR [27], which covers a total area of 292,700 km². It extends between 32°–36° N and 89°–103° E, with an elevation ranging from 2600 to 6600 m AMSL (Figure 1). The TRHR is considered the water tower of China and is an important nature reserve situated in QTP [28]. The region has a very complex topography, with highly elevated mountains as the basic structure of the landscape [29]. The region has a harsh and dry climate, with temperature ranging from −5.6 °C to 7.8 °C and precipitation stretching between 262 mm and 772 mm [28,30].

2.2. Data Description

The observed daily streamflow data of 5 hydrometric stations (i.e., Zhimenda, Xiangda, Jimai, Maqu, and Tangnaihai) were obtained from the Hydrology and Water Resources Survey Bureau of Qinghai province for the period 1980–2015 (Figure 1). Meteorological data (i.e., precipitation, maximum temperature, minimum temperature, wind speed, solar radiation, and relative humidity) of 23 climate stations (Figure 1) were obtained from the Qinghai Meteorological Bureau (QMB) for the period 1980–2015. All data (e.g., climate, discharge, elevation, soil, land cover, soil moisture, and water storage) required for this study were obtained from different sources, and data name, temporal resolution, spatial resolution, available period, and sources are described in

Table 1. Spatial and temporal data used in the present study.

SR	Data Type	Spatial/Temporal Resolution	Source	Availability
1	Streamflow data	Daily	Hydrology and Water Resources Survey Bureau of Qinghai province	1980–2015
2	Climate data	Daily	Qinghai Meteorological Bureau (QMB)	1980–2015
3	DEM	90 m	NASA’s Shuttle Radar Topography Mission (SRTM), Version 004 [31] (http://srtm.csi.cgiar.org) (accessed on 20 August 2019).	Updated 2008
4	Land Use Land Cover	1 km	Global Land Cover Characteristics [32] (https://earthexplorer.usgs.gov/) (accessed on 20 August 2019).	1993
5	Soil characteristics	1 km	Harmonized World Soil Database Version 1.2 (http://www.fao.org/soils-portal/) (accessed on 20 August 2019) [33]	Update 2013
6	Snow Water Equivalent (SWE)/snow depth	25 km/daily, monthly	Advanced Microwave Scanning Radiometer-Earth Observing System (AMSR-E), Version 2 [34]	2002/6–2011/10
		1°/daily	Canadian Sea Ice and Snow Evolution (CanSISE) [35]	1981–2010
		0.25°/daily	Environmental and Ecological Science Data Center for West China (WESTDC) [36]	1979–2019
7	Soil Moisture Content	0.25°/daily	European Space Agency (ESA) Climate Change Initiative Soil Moisture product (ESA-CCI-SM_v4.7),	1978–2019
		0.25°/daily	China Soil Moisture Dataset from Microwave Data Assimilation (ITP-LDAS) [37]	2002–2011
		0.5°′0.625°/diurnal	Modern-Era Retrospective analysis for Research and Applications (MERRA-2) [38]	1980–2020
8	GRACE data	300 km/monthly	Geo-forschungs-Zentrum Potsdam (GFZP), University of Texas-Center for Space Research (UT-CSR), and UT-CSR Mascons [39],	2002–2020
9	Leaf Area Index	0.25°/monthly	Global Monthly Mean Leaf Area Index Climatology, 1981–2015 [40]	1981–2015
10	Evapotranspiration	4 km/monthly	TERRACLIMATE [41]	1958–2019
		0.1°/monthly	Terrestrial evapotranspiration dataset across China, version 1.5 [42]	1982–2017
		500 m/8-daily	MOD16A2 Version 6 [43]	2000–2020

.

2.3. Setup of HEC-HMS

The Hydrological Modeling System (HEC-HMS), developed at the Hydrologic Engineering Centre (HEC) [44], has been applied around the world for flood modeling [45], water resource assessment [46], climate change impacts assessment [47], urban flooding [48], and streamflow simulation [49]. HEC-HMS encompasses seven methods to calculate precipitation losses (e.g., deficit and constant loss and SMA), seven transformation methods (e.g., Soil Conservation Services and Clark Unit Hydrograph), five baseflow estimation methods (e.g., Linear Reservoir and Recession methods), and six routing methods (e.g., Muskingum). There are two methods (i.e., Simple and Dynamic Canopy) to estimate interception and transpiration through plants, and a method to estimate water stored in surface depressions [44].

In this study, Thiessen Polygon, Penman–Monteith, TIM, SMA, Clark Unit Hydrograph, Muskingum Channel Routing, Linear Reservoir Baseflow Method, Dynamic Canopy, and Simple Surface Storage were included in the setup of HEC-HMS. The above-mentioned methods are elaborated in detail in the Supplementary Materials.

Thiessen Polygon was applied to interpolate daily climate data over each sub-basin (created in Section 2.4.1). This method is useful when gauges are scarce in a basin as compared to the basin area, as in Meenu, et al. [50]. For snowmelt modeling, each sub-basin was divided into different elevation bands (ranging from 3 to 5) for proper application of TIM.

2.4. Calibration and Validation

The auto-calibration does not provide optimal results and rational parameters in HEC-HMS, especially with TIM and SMA [21,22], and when we need to simulate other hydrological components such as soil moisture and snow water equivalent. Therefore, for the generation of all-important hydrological components, we adopted the manual calibration process in this study. For this purpose, we designed a comprehensive approach for the development of HEC-HMS, which can reduce time and easily manage many parameters. The whole procedure for the development of HEC-HMS is shown in Figure 2, which shows the required input data, tools used, model parameters, methods, and possible outputs, and is elaborated in the following steps.

2.4.1. Estimation of Physical Parameters

One of the tedious and time-consuming but important steps in the development of a complex hydrological model is to prepare input data [51], which can generally be divided into time-series data (e.g., climatic and hydrometric) and model parameters. The model parameters are of two kinds: physical parameters and process parameters [9,52]. Physical parameters represent the physical characteristics of watersheds [52], which are obtained by the watershed delineation that is considered the first step in hydrological modeling [53]. The delineation divides a basin into small sub-basins or grids, and their topographic characteristics such as sub-basin area, river length, sub-basin slope, and river slope are calculated. In this study, the TRHR was divided into 82 sub-basins (44 for SYAR, 29 for SYER, and 9 for SLR) during the delineation process of 90 m-SRTM-DEM [31], and physical characteristics were estimated for each sub-basin, as shown in Figure 3.

2.4.2. Estimation of Process Parameter

Generally, process parameters cannot be measured directly from the watershed, and these parameters need to be estimated indirectly through calibration [52]. To start the calibration process, some initial values of these parameters are required [9]. Sometimes, these initial values are assumed in cases when there is a simple model setup, with a very limited number of process parameters, and when we do not need a high level of satisfaction and must simulate only streamflow. However, this still requires a high level of expertise in the field. On the other hand, in the case of a very complex model setup including many parameters, it is essential to estimate some logical initial values of these parameters. This will not only help to simulate satisfactory hydrological components but also reduce the time of calibration. Additionally, this will provide rational and logical calibrated parameters. Therefore, in this study, we estimated these parameters in a systematic way.

The SMA method requires 5 initial conditions (i.e., canopy storage, surface storage, soil storage, groundwater layer 1 (GW1) storage, GW2 storage) and 13 process parameters (i.e., maximum infiltration, soil storage, tension storage, soil percolation, GW1 storage, GW1 percolation, GW1 coefficient, GW2 storage, GW2 percolation, GW2 coefficient, impervious surface area, maximum canopy storage, maximum surface depression storage), as described in Figure 2 [44]. The initial canopy and surface storages were taken as 0%, assuming no precipitation occurrence for a long time, as in Ahbari, et al. [54], and the rest were set as 30%. The initial values can be taken from 0% to 100%, with 0% showing almost no water available at the start of the simulation [25]. In the case of taking 0% as the initial conditions, a longer period (e.g., 5–10 years) will be required as a model spin-up period. In this study, we used 3 years of daily data as a spin-up (warm-up) period to reach the model at a steady state.

In contrast to the initial conditions, we need some reasonable and logical values for process parameters before starting the calibration process. These are estimated by using soil data and LCLU for each sub-basin on overage basis [2,22,25]. In this study, we used HWSD (Table 1) to extract soil properties such as depth and texture for the region. The soil texture is mainly used to identify the soil properties [24] because the variability in soil moisture patterns is most closely related to the variability in soil texture as compared to other soil properties [55]. Sandy loam and loam soil are dominant soil types in the SYER, covering an area of about 61% and 34% of the total area, respectively (Table 2). Porosity, field capacity, and saturated hydraulic conductivity (Table 2) were obtained from the literature, such as Clapp and Hornberger [56], J. Rawls, et al. [57], and Saxton and Rawls [58], using the soil texture of the region. The infiltration rates mentioned by FAO [59] against the texture of soil were assigned as initial values for maximum infiltration, and hydraulic conductivity as initial values for soil, GW1, and GW2 percolations, as in Fleming and Neary [24]. Soil and tension storages were estimated by multiplying the soil depth with porosity and field capacity, respectively, as in Ahbari, Stour, Agoumi and Serhir [54] and Samady [60]. The storage coefficients and storage depths of GW1 and GW2 were estimated using recession analysis on the observed streamflow, described in Fleming and Neary [24] and Samady [60]. The following steps were used to estimate the storage coefficients and storage depths of GW1 and GW2:

• First, search and separate single isolated storm event from streamflow time series.
• Plot a streamflow hydrograph for a single isolated storm event on a semi-logarithmic graph, as shown in Figure 4.
• Separate the baseflow and interflow from the streamflow hydrograph. We used the graphical method [61] for the separation of hydrograph components (Figure 4), as in Ahbari, Stour, Agoumi and Serhir [54] and Samady [60].
• Estimate the recession constant (k) for both the recession curves of baseflow and interflow. We used regression analysis on a semi-logarithmic graph. Baseflow is contributed by the deep groundwater layer (GW2) and interflow from the shallow groundwater layer (GW1) [24].
• Estimate the storage coefficient (S_c), storage capacity (S_t), and storage depth (S_t/area) with the following equations for both GW1 and GW2:

$$S_c=\frac{-1}{lnK}$$

(1)

$$S_t=S_c\times Q_t=\frac{-1}{lnK}{Q}_t$$

(2)

$$Q_t=Q_o{K}^t$$

(3)

where Q_o represents the initial baseflow (t = 0) and Q_t is the baseflow at time t (day). Equation (2) is the integration of the recession equation (Equation (3)) to estimate the storage capacity of a watershed. Similarly, estimate S_c and S_t for different storms in different seasons, where no rain occurs for some time after the peak flow [24]. The main step in the recessional analysis is to determine the recession constant. Typically, it ranges from 0.7 to 0.94 for interflow and 0.93 to 0.995 for baseflow [62]. However, a high value (e.g., >0.9) indicates the dominance of baseflow in a watershed [63]. In this study, we determined recession constants ranging between 0.92 and 0.98 for baseflow and between 0.6 and 0.88 for interflow. This procedure was applied to 5 gauging stations (Tangnaihai, Maqu, Jimai, Zhimenda, and Xiangda) and used for each sub-basin because data are available only on these stations. However, the exact values were estimated after calibration and validation processes.

Table 2. Soil properties in the source region of the Yellow River.

Soil Texture	Area (km2)	% of the Total Area	Soil Depth (mm)	Porosity (cm3/cm3)	Field Capacity (cm3/cm3)	KS (mm/h)	Infiltration (mm/h)
Clay	2481	2	1000	0.49	0.41	0.6	2
Clay loam	1926	2	1000	0.48	0.36	2.3	5
Loam	40,205	34	654	0.46	0.28	13.2	10
Sand	630	1	1000	0.4	0.1	210	25
Sandy loam	72,088	61	300	0.44	0.18	25.9	20
Silt loam	670	1	1000	0.49	0.31	6.8	7

Note(s): Sources: FAO (infiltration) [57,58], HWSD [33].

Table 2. Soil properties in the source region of the Yellow River.

Soil Texture	Area (km2)	% of the Total Area	Soil Depth (mm)	Porosity (cm3/cm3)	Field Capacity (cm3/cm3)	KS (mm/h)	Infiltration (mm/h)
Clay	2481	2	1000	0.49	0.41	0.6	2
Clay loam	1926	2	1000	0.48	0.36	2.3	5
Loam	40,205	34	654	0.46	0.28	13.2	10
Sand	630	1	1000	0.4	0.1	210	25
Sandy loam	72,088	61	300	0.44	0.18	25.9	20
Silt loam	670	1	1000	0.49	0.31	6.8	7

Note(s): Sources: FAO (infiltration) [57,58], HWSD [33].

Figure 4. Hydrograph components to estimate recession constant for GW1 and GW2 layers.

Figure 4. Hydrograph components to estimate recession constant for GW1 and GW2 layers.

The impervious surface area was obtained from the NASA Socioeconomic Data and Applications Center [64], which was only 0.09% of the total area of the SYER, and even less (0.03%) has been mapped by Zhang, et al. [65]. Max canopy (

Table 3. Maximum canopy storage in the headwater of the Yellow River, described in Ahbari et al. (2018), Samady (2017), and Verbeiren et al. (2016).

Code	Description	% of the Total Area	Canopy Storage (mm)
1	Urban and Built-Up Land	0.00	0.5
2	Dryland Cropland and Pasture	0.47	1.5
5	Cropland/Grassland Mosaic	0.01	2.0
6	Cropland/Woodland Mosaic	0.02	2.0
7	Grassland	88.61	2.0
8	Shrubland	1.99	2.5
9	Mixed Shrubland/Grassland	1.66	2.2
10	Savanna	1.05	2.0
11	Deciduous Broadleaf Forest	0.15	3.0
12	Deciduous Needleleaf Forest	0.13	2.0
14	Evergreen Needleleaf Forest	0.01	2.0
15	Mixed Forest	0.81	3.0
16	Water Bodies	1.51	0.0
17	Herbaceous Wetland	0.02	1.0
18	Wooded Wetland	0.00	1.0
19	Glacier	0.12	0.0
21	Wooded Tundra	3.43	2.0

) and surface-depression storages (

Table 4. Maximum surface storage according to percent slope [24].

Surface	Slope (%)	Max Surface Storage (mm)
Paved impervious area	NA	3.2–6.6
Steep, smooth slopes	>30	1.0
Moderate to gentle slopes	5–30	12.7–6.4
Flat, furrowed land	0–5	50.8

) were obtained from the literature, such as Ahbari, Stour, Agoumi and Serhir [54], Samady [60], and Verbeiren, et al. [66], using the vegetation types and basin slope in the region. Vegetation types were extracted using Global Land Cover Characteristics (Table 1), with 88% of grassland, and the slope was extracted using SRTM-DEM (Table 1), ranging from 0% to 185% (17% mean) in the SYER (Figure 3).

The LRBM requires 3 parameters (i.e., initial baseflow, groundwater coefficient, and reservoir) for each GW layer. Initial baseflow can be easily estimated from the observed streamflow, which can be low flows of streamflow, as used in this study. The groundwater coefficient, which is a response time of a basin [21], is the same as described above for SMA. Reservoirs are used to attenuate the streamflow; the greater the number of reservoirs, the more attenuation in the streamflow [44]. There are two parameters for the CUH method: time of concentration (TC) and basin storage coefficient. TC was estimated using a built-in function in HEC-GeoHMS, named the TR55 method, in which a basin (sub-basin) is divided into three flows such as sheet flow over the plane surface, shallow flow in small streams, and open channel flow. This requires much data about land and channel characteristics (e.g., land slope, land cover, width, depth, velocity, and roughness coefficient) for precise estimation of TC, which does not seem possible for a large and remote area. Nonetheless, we used Google Earth to measure the channel characteristics (e.g., cross-sectional area) at different locations to make an informed guess about TC. A detailed description of the TR55 method is given in [67]. The storage coefficient represents the temporary storage of excess precipitation in the basin as it drains to the outlet point [68], which can be computed as the flow at the inflection point on the falling limb of the hydrograph divided by the time derivative of flow [69]. Muskingum routing requires two parameters to estimate: K and X [44]. K is the travel time of a flood wave through a reach. Therefore, if inflow and outflow hydrographs are available, it can be estimated as the interval between similar points on inflow and outflow hydrographs (e.g., the centroid areas of two hydrographs) [23]. X is a dimensionless weighting factor, which ranges between 0 and 0.5 [44,70].

2.4.3. Sensitivity Analysis

After preparing input data and before calibration, another key process in manual calibration is to explore the sensitive parameters. Sensitivity analysis (SA) helps to recognize the parameters that have strong effects on the model outputs and hence influence the model response. SA also supports analyzing the interaction between parameters, its preferable range, and its spatial variability, which in turn influence the model outcomes. Consequently, SA accelerate the calibration process in hydrological modeling [71]

For this purpose, we used the local sensitivity analysis to determine the sensitive parameters in the region. According to [72], there are two types of sensitivity analysis: local and global. In the local, each parameter is changed one by one while keeping all other parameters unchanged, and its impact is measured on outputs. On the other hand, in the global, a model is run with the initial values of parameters, the changes in parameters and output are measured, and the percentage effect of each parameter on the output is calculated to determine the sensitive parameters. Sensitivity analyses were carried out by changing each parameter by 10% each time (between −40% and +40%) to observe the changes in simulated discharge. Two indicators, peak flow (PF) and total flow volume (TFV), were used to observe the effect of each parameter on discharge. All the parameters were ranked from top to bottom (the most sensitive to the least sensitive) and shown in Figure 5. In the case of PF, the top four sensitive parameters were maximum infiltration, storage coefficient, soil percolation, and surface storage; in the case of TFV, they were the percolations and storage coefficients of GW1 and GW2. This showed that maximum infiltration was most sensitive in the case of PF and percolation of GW2 in the case of TFV in the SYER. Instead of observing the parameter effects on PF and TFV, we also detected the parameter’s effects on other hydrograph components such as the rising limb, falling limb, and low flows. For this purpose, each parameter was increased by 10% (one by one) and its impact was observed on the hydrographs, which are shown in Figure 6. This gave a clear picture of which parameters affected which components of the hydrograph and made it easy to calibrate the model accurately.

In this study, HEC-HMS was calibrated for 2006–2010 and validated for the forward period (2011–2015) and the backward period (1981–2005) at Jimai, Maqu, Tangnaihai, Zhimenda, and Xiangda hydrometric stations. The other hydrological components such as soil moisture, baseflow, and TWS were validated for the periods that were available. For the evaluation of calibration and validation, four commonly used indicators (i.e., coefficient of determination (R²), Nash–Sutcliffe efficiency (E), normalized root mean square error (NRMSE), and percent volume deviation (PVD)) were used along with graphical plots, as in García, et al. [73], Mahmood, et al. [74], Meenu, Rehana and Mujumdar [50], and Verma, et al. [75].

$$R^2=\frac{{{\displaystyle \sum }}^{ }\left(Q_{obs }-\overline{Q_{obs}}\right)\times \left(Q_{sim }-\overline{Q_{sim}}\right)}{\sqrt{{{\displaystyle \sum }}^{ }{\left(Q_{obs }-\overline{Q_{obs}}\right)}^2\times {\left(Q_{sim }-\overline{Q_{sim}}\right)}^2}}$$

(4)

$$PVD (\%)=100\times \frac{Q_{sim}-Q_{obs}}{Q_{obs}}$$

(5)

$$E=1-\frac{{{\displaystyle \sum }}^{ }{\left(Q_{sim }-Q_{obs}\right)}^2}{{{\displaystyle \sum }}^{ }{\left(Q_{obs }-\overline{Q_{obs}}\right)}^2}$$

(6)

$$NRMSE=\frac{RMSE}{\overline{Q_{obs}}}=\frac{\sqrt{ \left(\frac{1}{n}{{\displaystyle \sum }}_{i=1}^n{(Q_{sim}-Q_{obs})}^2\right) }}{\overline{Q_{obs}}}$$

(7)

2.5. Baseflow Separation

Baseflow is one of the most important parameters of a hydrological cycle, contributing a large part to streamflow, but it is difficult to measure [6]. Many baseflow separation methods (graphical and digital filters) have been developed to estimate baseflow from the streamflow, which are described in [62,76,77]. Since observed baseflow is not available in the TRHR, and we need to evaluate baseflow simulated by HEC-HMS, we applied a frequently used digital filter known as the recursive digital filter described by Nathan and McMahon [62], as below:

$$f_k=\alpha f_{k-1}+\frac{\left(1+\alpha \right)}{2}\left[y_k-y_{k-1}\right]$$

(8)

where $f_k$ is quick flow at the kth time, α is a filter parameter, and $y_k$ is streamflow. Nathan and McMahon [62] recommended an α value of 0.925 for daily flow and 0.995 for hourly flow. In the present study, daily streamflow was used in the TRHR.

2.6. Terrestrial Water Storage

Terrestrial water storage (TWS) change is crucial for global as well as regional hydrological cycles and water resources management [78]. It is also essential for understanding a wide range of hydrological, climatological, and ecological processes. TWS mostly includes surface water, snow and ice, soil moisture, and groundwater storage. However, it is difficult to measure TWS directly [7]. In the present study, we used a hydrological model (i.e., HEC-HMS) to estimate the TWS in the TRHR. In the present study, TWS included canopy storage ($S_{canopy})$, surface storage ($S_{surface})$, soil storage ($S_{soil})$, groundwater storages from GW1 and GW2 ($S_{ground1},S_{ground2})$, and water stored in snow and ice as SWE ($S_{swe})$, as below:

$$TWS=S_{canopy}+S_{surface}+S_{soil}+S_{ground1}+S_{ground2}+S_{swe}$$

(9)

Since the in situ observations of TWS are not available in the TRHR, remote sensing data constitute a good source to evaluate TWS obtained by HEC-HMS. For this purpose, we obtained GRACE data of three solutions, i.e., GFZ, TU-CSR, and TU-CSR-Mascon for the period 2002–2015. TWS from GRACE is in the form of anomalies relative to 2004–2009. Therefore, we also calculated anomalies of TWS relative to 2000–2009 for a proper comparison.

2.7. Actual Evapotranspiration Estimation

Evapotranspiration is one of the most important components of the hydrological cycle and one of the most important physical processes in natural ecosystems. It links water and energy between the soil, land surface, and the atmosphere in the climate system [79,80]. The estimation of actual evapotranspiration (AET) is of vital importance in hydrological modeling, water resources management and planning, and related ecosystem services, as well as in addressing global climate change [79,81]. Evapotranspiration is also important for predicting plant productivity and species richness [80]. In the present study, the water balance equation was used to estimate AET, as given below:

$$AET=P-Q_{sim}-\Delta S$$

(10)

where P, Q_sim, and ΔS represent precipitation, simulated flow, and changes in water storage, respectively, in the basin. ΔS is generally assumed negligible for analysis in a longer period (≥10 years) [82] because it is difficult to measure. Nonetheless, in the present study, ΔS was estimated from TWS, as below:

$$\Delta S=TW{S}_t-TW{S}_{t-1}$$

(11)

3. Results

3.1. Calibration and Validation with Streamflow

Table 5. Performance evaluation of HEC-HMS in the TRHR on a daily basis.

	Jimai	Maqu	Tangnaihai	Zhimenda	Xiangda
Calibration
E	0.66	0.85	0.83	0.74	0.69
R2	0.73	0.88	0.89	0.83	0.73
PVD (%)	4.45	3.60	6.24	8.96	−4.29
NRMSE	0.46	0.30	0.30	0.51	0.43
Validation-1
E	0.69	0.83	0.90	0.77	0.74
R2	0.77	0.83	0.90	0.82	0.78
PVD	12.8	3.29	−1.32	8.68	6.66
NRMSE	0.39	0.35	0.27	0.49	0.49
Validation-2
E	0.62	0.82	0.82	0.79	0.61
R2	0.71	0.82	0.82	0.81	0.65
PVD	21.07	−4.38	−0.78	1.36	−1.13
NRMSE	0.52	0.42	0.40	0.48	0.55

describes the results of calibration and validation obtained from daily simulated and observed streamflow at five hydrometric stations. The E values ranged from 0.66 to 0.85 during calibration, 0.69 to 0.90 for validation-1, and 0.62 to 0.82 for validation-2. The results of R² were almost similar or slightly better than the E results during both calibration and validation. The PVD values varied between −4.5% and 9.0% for calibration, −1.5% and 13% for the first validation period, and −4.5% and 21% for the second validation period. On the other hand, NRMSE varied from 0.30 to 0.51 during calibration, 0.27 to 0.49 during validation-1, and 0.40 to 0.52 during validation-2. According to Van Liew and Garbrecht [83], E values greater than 0.75 are greatly appreciated, and the values between 0.36 and 0.75 are referred to as satisfactory depending on the objectives of the study. The model showed acceptable results on all gauges, though the numbers of climate stations are sparse, especially in the SYAR. Similar results have been reported by Meng, et al. [84], Yang, et al. [85], and Zhang, et al. [86] in the SYER. For example, E values ranging from 0.78 to 0.91 for calibration and 0.70 to 0.89 for validation at Tangnaihai have been reported by Meng, Su, Yang, Tong and Hao [84], Yang, Wang, Yu, Krysanova, Chen, Schwartz and Sudicky [85], and Zhang, Su, Hao, Xu, Yu, Wang and Tong [86]. Furthermore, E values ranging from 0.67 to 0.81 for calibration and 0.61 to 0.78 for validation at Maqu and from 0.51 to 0.76 for calibration and 0.49 to 0.61 for validation at Jimai have been reported by Meng, Su, Yang, Tong and Hao [84]. Guo, et al. [87], Wang, et al. [88], and Zhang, et al. [89] also showed similar results of calibration and validation for the SYAR and SLR, ranging from 0.72 to 0.87.

The graphical presentation of results is another important visual tool for the evaluation of model performance subjectively and qualitatively, showing how well simulated data follow the variations (e.g., peak and low flows) in observed data. Therefore, the simulated streamflow was plotted against the observed discharge to visualize the performance of the developed model at five gauges in the region, as shown in Figure 7. At Jimai, the peaks were not well followed by the model during both calibration and validation as compared to low flows. For example, during calibration, the model overestimated the peaks in 2007, 2008, and 2010, and during validation-2, the model not only overestimated the peaks in some years but also overestimated low flows. This might be due to the lack of climate data, the existence of snow/ice, and complex topography. At Maqu and Tangnaihai, both peaks and low flows were better simulated by the model except in some years. For example, the model overestimated peaks in 2011 during validation-1 at Maqu and in 1995 and 1996 during validation-2 at Maqu and Tangnaihai stations. At Zhimenda, the model overestimated slightly in the case of peaks and low flows during calibration and in the case of low flow during validation-1. The model also followed well during the second validation period, except for some years such as 1982 and 2001. At Xiangda, the peaks were also not well followed by modeled streamflow, especially during both the validation periods. Nonetheless, low flows were simulated well by the model. It was also noted that the rising and recession limbs of observed streamflow were well followed by the simulated streamflow on all gauges. From the above discussion and comparison with other studies, the model showed the capability of simulating streamflow.

3.2. Validation of Other Hydrological Components

In this study, soil moisture content, baseflow, terrestrial water storage anomalies, snow water equivalent, and AET were evaluated with reanalysis and remote sensing data (RRSD) obtained from freely available sources because in situ data are not available in the region. Recently, the RRSD has been considered a suitable alternative in the scientific community because it provides promising, convenient, and relatively easy-to-use input data for modeling in remote regions, such as the TRHR [90,91]. For the evaluation, three statistical indicators, i.e., correlation coefficient (R), percent volume deviation (PVD), and root mean square error (RMSE) were calculated on a monthly basis (

Table 6. Evaluation of HEC-HMS using different hydrologic components in the TRHR.

	Soil Moisture Content			Snow Water Equivalent			Baseflow	Terrestrial Water Storage Changes			Actual Evapotranspiration
	ESA-CCI-SM	MERRA-2	ITP-LDAS	CanSISE	AMSR-E	WESTDC	DRF	GFZ	UT-CSR	UT-CSR_MASCON	TERRACLIMATE	MOD16A2	TEDAC
Analysis Period 2001–2010				2001–2010			1981–2015	2002–2015			2000–2010
The Source of the Lancang River (SLR)
R	0.72	0.60	0.44	0.51	0.75	0.72	0.95	0.72	0.70	0.76	0.89	0.79	0.87
RMSE	0.031	0.027	0.033	6.5	12.5	6.0	29.0	20.0	22.1	17.5	15.3	23.9	24.3
PVD (%)	4.5	−6.6	9.0	64.1	−51.7	111.8	18.1	11.0	−5.3	34.3	0.1	−38.6	−36.3
The Source of the Yellow River (SYER)
R	0.72	0.74	0.56	0.82	0.73	0.81	0.96	0.70	0.75	0.79	0.85	0.80	0.83
RMSE	0.029	0.020	0.022	3.6	6.3	2.5	100.8	21.7	20.4	17.7	21.9	24.5	22.2
PVD (%)	6.7	5.7	−2.4	−38.4	−59.8	22.4	7.4	17.4	9.1	18.6	−13.5	−33.2	−25.4
The Source of the Yangtze River (SYAR)
R	0.80	0.84	0.49	0.24	0.33	0.32	0.96	0.64	0.66	0.55	0.86	0.73	0.81
RMSE	0.034	0.015	0.035	9.7	8.2	9.9	113.2	19.2	19.2	28.7	17.3	23.6	24.6
PVD (%)	11.5	−0.9	−10.2	318.4	−11.3	410.4	20.1	16.3	5.1	−23.8	−3.5	−42.2	−37.6

) because the RRSD data are mostly available in a monthly format.

3.2.1. Soil Moisture Content

Soil moisture content (SMC) simulated by HEC-HMS was compared with the data obtained from MERRA-2, ITP-LDAS, and ESA-CCI-SM products (described in Table 1) for the period 2001–2010 over the TRHR. On the whole, the R values ranged from 0.44 (ITP-LDAS) to 0.84 (MERRA-2), the RMSE values from 0.015 (MERRA-2) to 0.035 (ITP-LDAS), and the PVD values from −10.2 (ITP-LDAS) to 11.5 (ESA-CCI-MS) in the TRHR (Table 6). In the case of the SLR, the highest correlation and lowest deviation (PVD) were obtained with ESA-CCI-MS, followed by MERRA-2. However, RMSE was the lowest for MERRA-2, followed by ESA-CCI-MS (Table 6). Figure 8a shows that the model’s SMC followed all three datasets well, but overestimated from the ITP-LDAS data in the SLR. In the case of the SYER, the highest R and the lowest RMSE values were calculated with MERRA-2, followed by ESA-CCI-MS. However, PVD values were the lowest for ITP-LDAS, followed by MERRA-2 (Table 6). Figure 8b shows that the model also captured the annual cycle of all three datasets well. In the case of the SYAR, better results were obtained for MERRA-2, followed by ESA-CCI-MS (Table 6). However, the model overestimated ESA-CCI-MS and underestimated IPT-LDAS (Figure 8c). We also compared SMC with previous studies such as Deng, et al. [92], Yuan, et al. [93], Zeng, et al. [94], and Zhang, et al. [95] that have used observed SMC. According to these studies, observed SMC at Maqu and Maduo varies between 0.07 and 0.38 (mm/mm) and 0.01 and 0.18 (mm/mm), respectively, during an annual cycle, with the maximum values in the wet season. Nonetheless, the mean annual values of the present study at Maqu and Maduo were comparable with the observations.

3.2.2. Baseflow

For the evaluation of baseflow simulated by HEC-HMS, we used baseflow separated by the recursive digital filter (RDF) [62] for the period 1981–2015. Table 6 shows that the simulated baseflow has a high correlation (>0.94) with RDF in all three river basins. The values of RMSE were 29 m³/s, 100 m³/s, and 113 m³/s in the SLR, SYER, and SYAR, respectively. PDV was above 15% in the SLR and SYAR and below 10% in the SYER. Figure 9 shows that the model overestimated in all months except in July in the SYER. In the low flow months (November–April), the snowfall period, baseflow must be the same as the streamflow, but baseflow from RDF was lower than the streamflow. Nevertheless, the baseflow simulated by HEC-HMS was almost the same as the streamflow in the low flow months. Therefore, we can conclude that baseflow simulated by HEC-HMS was more rational than RDF. Baseflow simulated by the model contributed 77–78% to streamflow, while RDF contributed 65–73% to streamflow on an annual basis. Some studies, such as Chen, et al. [96], Liu, et al. [97], Lu, et al. [98], and Qian, et al. [99], have simulated or separated baseflow in the TRSR, but all showed different results. Chen, Zheng, Chen and Liu [96] and Qian, Wan, Wang, Lv and Liang [99] determined a baseflow of 65–78% of streamflow using different digital filters in the SYER and SYAR; Liu, Yao, Wang and Yu [97] explored only about 30% of baseflow using Global Land Data Assimilation System (GLDAS) products; and Lu, Wang, Shao, Yu, Hao, Xing, Yong and Li [98] simulated 50–54% of baseflow using the VIC model in the SYAR and SYER. During the recession analysis (Section 2.4.2, we estimated high values of recession constant, confirming baseflow dominance in the region.

3.2.3. Terrestrial Water Storage

In the present study, terrestrial water storage (HEC-TWS) estimated by HEC-HMS was evaluated with GRACE data, as in Yuan, Ji, Wang, Liang, Yang, Ye, Su and Wen [93]. Since GRACE-TWS is in the form of anomalies relative to the mean of 2004–2009, we also calculated anomalies of HEC-TWS in the region, the same as GRACE. For the evaluation, we used three solutions of GRACE (i.e., GFZ, UT-CSR, and UT-CSR-Mascon) for 2002–2015. The HEC-TWS anomalies had high correlations with three solutions, ranging from 0.70 to 0.76 in the SLR, 0.70 to 0.79 in the SYER, and 0.55 to 0.66 in the SYAR. RMSE ranged between 17 mm and 29 mm in all three basins for three solutions, and PVD ranged between −24% and 35%. According to all three indicators, the HEC-TWS anomalies were best captured by UT-CRS, followed by the GFZ. Figure 10 shows that the annual patterns of GRACE-TWS anomalies were well followed by HEC-TWS in all three headwaters. However, the model underestimated during the dry months (November–March) and overestimated during the wet period, especially in June–August. For the comparison of TWS with the previous studies, we could not find any study simulating TWS using a hydrological model, especially using a conceptual model, for a longer period than GRACE-TWS, except Yuan, Ji, Wang, Liang, Yang, Ye, Su and Wen [93]. They simulated for the period 1979–2014 using a land surface model, but they also validated their results with GRACE-TWS. However, most studies such as Huang, et al. [100], Jing, et al. [101], and Xu, et al. [102] used GRACE-TWS either for the evaluation of Interim Reanalysis Data and GLDAS data or the assessment of other hydrological components such as groundwater and evapotranspiration.

3.2.4. Snow Water Equivalent

For the evaluation of SWE simulated by HEC-HMS (HEC-SWE), we obtained SWE data from three datasets, i.e., CanSISE, AMSR-E, and WESTDC (Table 1), for the period 2001–2010. The R values ranged from 0.51 (CanSISE) to 0.75 (AMSR-E) in the SLR, from 0.73 (AMRS-E) to 0.82 (CanSISE) in the SYER, and from 0.24 (CanSISE) to 0.33 (AMRS-E) in the SYAR (Table 6). PVD showed that the model underestimated (11–60%) AMSR-E and overestimated (about 22–410%) from CanSISE and WESTDC in all three basins, most of the time. Although correlations were high in the SLR and SYER, PVD and RMSE were also high. Figure 11 clearly shows that the model overestimated from CanSISE and WESTDC but underestimated from AMRS-E in all three basins, except in the SYER, where the model also underestimated from CanSISE. Dai, et al. [103] and Yang, et al. [104] also explored that the estimated SWE by observed snow depth was underestimated by AMRS-E-SWE, especially over grassland areas, in QTP and the Tibetan Plateau, respectively. Significant differences between in situ and remote sensing products have also been reported over the Contiguous United States [105].

Figure 11 displays that CanSISE and WESTDC were comparable with each other but underestimated from AMSR-E in all three basins. On the whole, we can say that HEC-SWE is not comparable with the remote sensing products, but we cannot say that it is wrong because these products are also not comparable with each other and have some uncertainties. For the evaluation of HEC-SWE, we estimated accumulated snowfall during the snowfall period (October–March) in the SLR and plotted along with remote sensing products and HEC-SWE, as shown in Figure 12. The pattern of HEC-SWE was well in line with the pattern of snowfall in the SLR. For example, HEC-SWE started to accumulate from mid-October and reached a peak in March, very similar to the accumulated snowfall. This indicates that our results are rational because SWE should be maximum in March because the temperature remains below the freezing point from October to March in the SLR (Figure 12). However, the remote sensing products showed peaks for SWE in November or December. Since uncertainties are high in the remote sensing datasets, we need field observations for an accurate evaluation of HEC-SWE.

3.2.5. Actual Evapotranspiration

For the evaluation of AET estimated by HEC-HMS (HEC-AET), three products, MOD16A2, TERACLIMATE, and TEDAC, were used from 2000 to 2010. Table 6 shows high correlations with all three products, ranging from 0.73 (MOD16A2) to 0.89 (TERACLIMATE) in all three basins. The highest correlations were found with TERACLIMATE, followed by TEDAC in all three basins. The PVD values were in the satisfactory range for TERACLIMATE in all three basins, ranging from −13.5% to 0.1%. However, the PVD values for MOD16A2 and TEDAC were not satisfactory in all three basins, ranging from −42% to −25%. The PVD values showed that the model underestimated all three products in all three basins, except TERACLIMATE in the SLR.

Figure 13 shows the mean monthly AET estimated by HEC-HMS and three products for the period 2000–2010. The monthly precipitation was also plotted with AETs, which helped to evaluate the estimated AET. The patterns of TERACLIMATE and TEDAC were followed well by HEC-AET in all three basins. However, HEC-AET was underestimated by both products in most months, especially during June–August. MOD16A2 highly overestimated estimated AET during cold and dry months in all three basins, which does not appear rational. We can also evaluate HEC-AET with precipitation for a long period annual mean. It was noted that the mean annual HEC-AET, TERACLIMATE, TEDAC, and MOD16A2 were 60%, 60%, 87%, and 90% of precipitation, respectively, in the SLR; 77%, 79%, 122%, and 131% of precipitation in the SYAR, respectively; and 72%, 83%, 97%, and 108% of precipitation in the SYER, respectively. The AET values of TEDAC and MOD16A2 do not look rational because AET from these products contained a high percentage (>90%) of precipitation. Previous studies such as Li, et al. [106], Sato, et al. [107], Xue, et al. [108], and Zhang, Su, Yang, Hao and Tong [89] have also reported an AET of about 70% of the precipitation, similar to this study, and an even smaller percentage (60% of precipitation) was reported by Bei, et al. [109]. From the above discussion, it can be concluded that HEC-AET and TERACLIMATE appear rational compared to TEDAC and MOD16A2 in the TRHR.

3.3. Uncertainties and Limitations

In general, there are four main sources of uncertainties in hydrological modeling [110], i.e., input data, output data, model structural, and parametric uncertainties. One uncertainty related to data in the region is that the data required to simulate hydrological components are not enough in this region of complex topography. For example, there are only 6–10 climate stations in and around the SYAR (158, 027 km²). This means that the gauge density (about 13,000 km²/gauge) is so sparse compared to the density suggested by the World Meteorological Organization (250–1000 km²/gauge) in similar kinds of harsh conditions [111]. Although reanalysis and remote sensing datasets (RRSD) are suitable alternatives, these also contain uncertainties, as we observed in Section 3.2.4 For example, we observed that SWE data from RRSD products showed large differences from each other (Figure 11). This can be due to the different tools and techniques used in the retrieval of these datasets. Despite all advances in the development of models, they still cannot entirely represent the hydrological system of a region because of a lack of enough understanding of the hydrological cycle. According to [112], the estimation of parameters during the calibration process can be another source of uncertainty in a hydrological model because it is impossible to determine the exact parameters of a study area.

There are also some limitations in HEC-HMS used in the present study. For example, it has only one option (Temperature Index Method to estimate the snowmelt water, which mainly depends on temperature. However, snowmelt is a complex method, which requires a complete energy balance along with the water balance approach. Additionally, some basic parameters such as PX, base temperature, and wet-melt rate in TIM are used for the whole basin. However, these should be used on sub-basin level, even on an elevation-band level in each sub-basin, because these can vary spatially in a region, especially a region having complex mountainous terrain such as the TRHR. In HEC-HMS, instead of precipitation and temperature, the other climatic variables such as wind speed, relative humidity, and sunshine hours are not spatially distributed properly using some sophisticated techniques such as inverse distance weighting or Thiessen Polygon methods. For example, each sub-basin is assigned by only one climate station located in or around the sub-basin.

In this study, we used point-based climate station data, which are very small in number as compared to the area of the region. In addition, the model was calibrated and validated for only streamflow gauges. In the future, we can use gridded satellite climate data such as the Tropical Rainfall Measuring Mission (TRMM). In addition, the model results can be compared with other hydrological model (e.g., SWAT) using the same mythology.

4. Conclusions

In this study, we wanted to configure HEC-HMS so that it could simulate all the important hydrological components (e.g., streamflow, baseflow, soil moisture content, terrestrial water storage, and snowmelt water) in the Three-River Headwater Region (TRHR). However, the main problem is the inefficient built-in auto-calibration methods in HEC-HMS, especially when we include Soil Moisture Accounting (SMA) and Temperature Index Method (TIM). Secondly, no studies have been reported having a comprehensive approach that can be used to configure HEC-HMS in such a way that all the important hydrological components (e.g., streamflow, baseflow, soil moisture content, terrestrial water storage, snowmelt water, and evapotranspiration) can be simulated. Therefore, in this study, we designed a comprehensive approach to configure HEC-HMS so that the model can generate all the important components satisfactorily. The approach was based on providing a systematic way to determine some rational initial values for the process parameters, sensitive parameters, and relations between parameters and hydrograph components (e.g., low flow, peaks, and rising and falling limbs). This will reduce the time of calibration and enhance the deep understanding of the hydrological processes.

The results showed that HEC-HMS has the capability to simulate almost all the important components satisfactorily if we configure the model properly and rationally. This study will be a guideline for the modelers and users of HEC-HMS and similar kinds of other hydrological models.