A Flexible Framework HydroInformatic Modeling System — HIMS

It is important to simulate streamflow with hydrological models suitable for the particular study areas, as the hydrological characteristics of water cycling processes are distinctively different due to spatial heterogeneity at the watershed scale. However, most existing hydrological models cannot be customized to simulate water cycling processes of different areas due to their fixed structures and modes. This study developed a HydroInformatic Modeling System (HIMS) model with a flexible structure which had multiple equations available to describe each of the key hydrological processes. The performance of the HIMS model was evaluated with the recommended structure for semi-arid areas by comparisons with two datasets of observed streamflow: the first one of 53 Australian watersheds, the second one of the Lhasa River basin in China. Based on the first dataset, the most appropriate watersheds were identified for the HIMS model utilization with areas of 400–600 km2 and annual precipitation of 800–1200 mm. Based on the second dataset, the model performance was statistically satisfied with Nash-Sutcliffe Efficient (NSE) greater than 0.87 and Water Error (WE) within ±20% on the streamflow simulation at hourly, daily, and monthly time steps. In addition, the water balance was mostly closed with respect to precipitation, streamflow, actual evapotranspiration (ET), and soil moisture change at the annual time steps in both the periods of calibration and validation. Therefore, the HIMS model was reliable in estimating streamflow and simulating the water cycling processes for the structure of semi-arid areas. The simulated streamflow of HIMS was compared with those of the Variable Infiltration Capacity model (VIC) and Soil and Water Assessment Tool (SWAT) models and we found that the HIMS model performed better than the SWAT model, and had similar results to the VIC model with combined runoff generation mechanisms.


Introduction
Hydrological models are important tools to simulate water cycling processes in studying water resource assessment, management, and utilization.Based on the introduction of the unit hydrograph concept [1] and the recognition of the water cycling processes including infiltration, evaporation, streamflow concentration, lumped hydrological models have been developed including the Stanford model [2], Soil Conservation Service model (SCS) [3], and Xinanjiang Model [4].Since the concept was introduced of the distributed hydrological model by Freeze and Harlan [5], semi-distributed and distributed models have been developed and widely applied such as the SWAT model [6], the Hydrologic System Program-Fortran (HSPF) [7], and the VIC model [8], etc.
Spatial heterogeneity is taken into consideration by the distributed hydrological model in the model structures and methods describing the hydrological processes [9,10].Hydrological models discretize the watersheds into sub-watersheds or grid cells [6,8], or hydrological response units (HRUs) [6], and representative elementary areas (REAs) [11] based on the datasets of digital elevation model (DEM), land cover, and soil.They are defined by the distributed hydrological models as the basic spatial modeling unit [9,10] with homogenous soil properties, land covers, and hillslopes, on which the hydrological processes are simulated with the models.Spatial heterogeneity within the unit is presented with the infiltration/rainfall excess estimation process.The Xinanjiang model [4] simulates the variable infiltration capacity within the grid cell with the probability distribution function and then incorporates the effects of soil properties on the hydrological conductivity.The soil conservation service curve number (SCS-CN) method determines the standardized rainfall-runoff relationship accounting for the spatial heterogeneity of the soil, land use, and topography [12].Runoff generation mechanisms of intensity excess in the arid areas and saturation excess in the humid areas are combined dynamically with the time compression analysis (TCA) method to improve the runoff predicting accuracy in semi-arid areas [13,14].
Spatial heterogeneity among units is presented by the hydrological models with the methods of overland flow concentration and channel routing through which overland flow concentrates in the channels and reaches the watershed outlet [9].The kinematic wave model simulates the overland flow movement along hillslopes through transferring the rainfall excess into hydrographs and accounting for the slopes, aspects, and hillslope roughness [15].The lumped Geomorphologic Instantaneous Unit Hydrograph (GIUH) model accounts for the influence of the basin geomorphic structure on the watershed outlet hydrograph by taking topographic characteristics into consideration based on the Horton-Strahler theory [16,17].The distributed channel routing model of the Saint-Venant (S.V.) equation estimates the flow rate, velocity, and depth at cross sections along the channel accounting for the distance along the watercourse, the cross section area, and the watercourse bottom slope based on the laws of mass and momentum conservation [15].Numerical schemes have been developed to obtain solutions of the S.V. equations such as the finite element, finite difference method, and the Preissman implicit scheme which work under the specific hydrological conditions of the basins [18].
To explicitly describe the water cycling processes of specific study areas with hydrological models, the model structure needs to be flexible enough such that alternative models/additional processes can easily be implemented.It is hard to customize the hydrological model according to the hydrological characteristics of the study areas with the existing models due to their fixed structures such as the VIC model, SWAT model, HSPF, etc.In this paper, a hydrological model HIMS is developed with modularized structures and multiple choices describing the key hydrological processes.The HIMS model with flexible structures can be customized by the users according to the distinctive hydrological characteristics of the study areas.With the recommended structure of the HIMS model for the semi-arid area, the most appropriate watersheds were identified for the ranges of size and annual precipitation for the observed streamflow of 53 Australian watersheds.The performance of the HIMS model was evaluated with the observed streamflow datasets of the Lhasa River watersheds in China, and the simulated streamflow with the HIMS model was compared with those simulated with the VIC and SWAT models.

Materials and Methods
The HIMS model is a hydrological model with a grid cell or sub-watershed as its basic unit of water balance.Each grid/sub-watershed includes hillslopes and one channel, grids/sub-watersheds are linked with stream networks of different orders based on the spatial topological relationships as a distributed hydrological model.The HIMS model has multiple choices of equations within each modular of the hydrological processes including canopy interception, snow melt, evapotranspiration (ET), runoff generation/infiltration, sub-surface flow, overland flow routing, and channel flow routing.The flexible structure enables the HIMS model to be customized to the specific study areas.
Figure 1 shows the customization of the HIMS model with the key processes based on the omission of other processes.The runoff generation modular contains the SCS-CN, the Green-Ampt equation, the LCM developed by Liu et al. [19,20], and the Time-Variant Gain Model by Xia et al. [21].Though both the SCS and LCM indicate reciprocal linear relationship between precipitation losses and total precipitation volume [22], the LCM equation is much more advanced than the SCS-CN method when it comes to the utilization at un-gauged areas.The key parameters of SCS are hard to obtain, while the parameters of LCM are easily obtained with field experiments.Therefore, LCM is much more suitable for estimating runoff of the un-gauged study areas [22].Hillslope flow routing contains the kinematic wave equation and the unit hydrograph, which are used to calculate the flow velocity, concentration time from hill-slopes to the channel, and/or the flow concentration process of the sub-watershed/grid cell, including peak flow, peak lag time, and duration of the grid cell.Channel routing is used to obtain the water quantity and velocity of a cross section of a channel or watershed in a short time.This modular contains the Muskingum method and the kinematic wave equation with the newly developed solving method by Li [23,24].This new solving method enables the kinematic wave equation to be applied to watersheds with dry channels or complex stream networks or discretized boundary conditions.If a reservoir exists between channels, then the most widely used linear reservoir method can be selected from the reservoir routing modular to simulate the regulating process of streamflow.
Figure 1 shows the customization of the HIMS model with the key processes based on the omission of other processes.The runoff generation modular contains the SCS-CN, the Green-Ampt equation, the LCM developed by Liu et al. [19,20], and the Time-Variant Gain Model by Xia et al. [21].Though both the SCS and LCM indicate reciprocal linear relationship between precipitation losses and total precipitation volume [22], the LCM equation is much more advanced than the SCS-CN method when it comes to the utilization at un-gauged areas.The key parameters of SCS are hard to obtain, while the parameters of LCM are easily obtained with field experiments.Therefore, LCM is much more suitable for estimating runoff of the un-gauged study areas [22].Hillslope flow routing contains the kinematic wave equation and the unit hydrograph, which are used to calculate the flow velocity, concentration time from hill-slopes to the channel, and/or the flow concentration process of the sub-watershed/grid cell, including peak flow, peak lag time, and duration of the grid cell.Channel routing is used to obtain the water quantity and velocity of a cross section of a channel or watershed in a short time.This modular contains the Muskingum method and the kinematic wave equation with the newly developed solving method by Li [23,24].This new solving method enables the kinematic wave equation to be applied to watersheds with dry channels or complex stream networks or discretized boundary conditions.If a reservoir exists between channels, then the most widely used linear reservoir method can be selected from the reservoir routing modular to simulate the regulating process of streamflow.

Structure of the HIMS Model and Key Equations
In each grid cell/sub-basin, the hydrological cycling processes are described in Figure 2 and the equations are listed in Table 1 including precipitation, evapotranspiration (ET), snowmelt, interception, infiltration, surface flow, subsurface flow, and baseflow.

Structure of the HIMS Model and Key Equations
In each grid cell/sub-basin, the hydrological cycling processes are described in Figure 2 and the equations are listed in Table 1 including precipitation, evapotranspiration (ET), snowmelt, interception, infiltration, surface flow, subsurface flow, and baseflow.
Precipitation falls on the ground as rainfall or snow after the interception process with involvement of the Degree-Day equation [25] taking the snow water equivalent (SWE) into consideration.Potential Evapotranspiration (ET) is simulated with the Hargreaves equation [26] accounting for the impacts of the average, maximum, and minimum air temperatures.The actual ET is calculated for each vegetation type with the single crop coefficient method of FAO-56 [27,28] by calculating evaporation and transpiration separately and accounting for the impacts of leaf area index (LAI) and soil moisture of the root zone.For each sub-watershed, the actual ET is the sum of each vegetation type within the sub-watershed.Runoff is generated by HIMS accounting for two runoff generation mechanisms: intensity excess and saturation excess with a soil structure of three layers-root zone, transition zone, and saturated layer.The mechanism selection on runoff generation depends on the comparison result of soil water content between the transmitted zone and the saturated soil layer.If the water content of the transmitted zone and root zone is less than that of the saturated soil layer, the intensity excess mechanism is adopted, otherwise the saturated excess mechanism is selected.
The streamflow of hillslopes is composed of the surface flow, lateral flow, and the partial baseflow indicated by Figure 2 and calculated with Equations 6-9 in Table 1.The streamflow of hillslopes within the sub-watershed is concentrated into one channel with the overland flow routing models like Equation 10 in Table 1.The streamflow of sub-watersheds is concentrated to the watershed outlet with channel routing models like Equation 11 in Table 1.The HIMS ensures that the remaining part of the baseflow from hillslopes goes into the next sub-watershed directly through the ground water movement, as a way of describing the sub-watersheds water communication with underground water.If reservoirs exist, reservoir routing is simulated with equations as Equation 12 in Table 1.The HIMS model is suitable for utilization at areas with baseflow dominated the hydrographs, such as the areas with Karst landform, or areas with temporal seasonal rivers, or the head areas of sub-watersheds with low undercut channels.As irrigation is taken into consideration by the root zone water storage calculation, HIMS is also suitable for the agriculture dominated watersheds.
Water 2018, 10, x FOR PEER REVIEW 4 of 17 Precipitation falls on the ground as rainfall or snow after the interception process with involvement of the Degree-Day equation [25] taking the snow water equivalent (SWE) into consideration.Potential Evapotranspiration (ET) is simulated with the Hargreaves equation [26] accounting for the impacts of the average, maximum, and minimum air temperatures.The actual ET is calculated for each vegetation type with the single crop coefficient method of FAO-56 [27,28] by calculating evaporation and transpiration separately and accounting for the impacts of leaf area index (LAI) and soil moisture of the root zone.For each sub-watershed, the actual ET is the sum of each vegetation type within the sub-watershed.Runoff is generated by HIMS accounting for two runoff generation mechanisms: intensity excess and saturation excess with a soil structure of three layersroot zone, transition zone, and saturated layer.The mechanism selection on runoff generation depends on the comparison result of soil water content between the transmitted zone and the saturated soil layer.If the water content of the transmitted zone and root zone is less than that of the saturated soil layer, the intensity excess mechanism is adopted, otherwise the saturated excess mechanism is selected.
The streamflow of hillslopes is composed of the surface flow, lateral flow, and the partial baseflow indicated by Figure 2 and calculated with Equations 6-9 in Table 1.The streamflow of hillslopes within the sub-watershed is concentrated into one channel with the overland flow routing models like Equation 10 in Table 1.The streamflow of sub-watersheds is concentrated to the watershed outlet with channel routing models like Equation 11 in Table 1.The HIMS ensures that the remaining part of the baseflow from hillslopes goes into the next sub-watershed directly through the ground water movement, as a way of describing the sub-watersheds water communication with underground water.If reservoirs exist, reservoir routing is simulated with equations as Equation 12in Table 1.The HIMS model is suitable for utilization at areas with baseflow dominated the hydrographs, such as the areas with Karst landform, or areas with temporal seasonal rivers, or the head areas of sub-watersheds with low undercut channels.As irrigation is taken into consideration by the root zone water storage calculation, HIMS is also suitable for the agriculture dominated watersheds.
Canopy interception depends on minimum of precipitation canopy interception ability and interception storage lack where canopy interception Canopy interception depends on minimum of precipitation canopy interception ability and interception storage lack where canopy interception ability is relative to leaf area index Actual ET is influenced by soil, climate, and crop type Conceptual model, where R and r can be chosen for various land types (other infiltration equations also exist in HIMS as choice candidates) Runoff is generated according to the soil water content, two scenarios of Horton and Dunne flows are taken into consideration.
Based on water balance equation d c , vegetation coverage degree, reflecting the spatial distribution of vegetation; P(t), rainfall (mm); W cu (t), water for shortage the canopy (mm); K c , canopy interception coefficient, related to vegetation types; LAI(t), the vegetation's leaf area index; D f , the factor for degree-day model, ranging from 3 × 10 −4 -6 × 10 −4 (m • C −1 h −1 ); T(t), the air temperature ( • C); T b , the initial temperature when snow begins to melt ( • C); a and b, parameters, arranging from 0.0023 to 0.0032 and from 0.5 to 0.6, respectively; RA max , the greatest solar radiation if possible (MJm −2 d −1 ); T, T max, and T min , daily average, maximum and minimum temperature, respectively ( • C); L, latent heat of vaporization (MJkg −1 ); Ws(t), non-saturated soil moisture content (mm); SM w , maximum soil storage capacity (mm) of root zone; SW, soil water content of root zone (mm); R and r, parameters related to land use and soil moisture; f, infiltration capacity (mm); FW, soil water content of transition zone (mm); FW m , maximum soil water content of transition zone (mm); L a , coefficient of sub-surface runoff; REC, recharge into groundwater (mm); V s , average water flow velocity in the slope section (m s −1 ); K b , coefficient of baseflow; GW s , ground water storage-volume (mm); h, hydraulic radius (m); S i , hydraulic gradient for the hillside slope; m, average roughness of the slope; x and y, indexes; V B , flood wave velocity along the way (m s −1 ); S, average gradient; P n , net rainfall (mm); k 1 , ratio of average flow velocity to the flood wave velocity; K B , ratio of flood wave velocity to mean velocity, related to the shape of the river section; K and m, parameters in the power relationship between peak flow and net rainfall; av and am, parameters in the power relationship between the average velocity and peak flow; belta1, convergence index; Q in1 and Q in2 , flow into river at the beginning and end of a time step (m 3 s −1 ); Q out1 and Q out2 , flow out of river at the beginning and end of a time step (m 3 s −1 ); C 1 , C 2 and C 3 , parameters; I r , flow reservoir (m 3 s −1 ); Q r , flow released from the reservoir (m 3 s −1 ).

Model Calibration and Evaluation
It is necessary to conduct hydrological analysis of the study areas before model utilization, in order to construct the model structure according to the hydrological characteristics based on the observed dataset.In the semi-arid areas, it is recommended to use the LCM as the runoff generation equation and the Muskingum equation as the channel routing model based on the literature review [22][23][24] for setting the HIMS model structure.It makes senses as the HIMS model has the structure of the combined runoff generation mechanisms and LCM is an empirical model developed on the base of the observed data from the semi-arid areas including the Lhasa River basin [21].
The HIMS model was evaluated with the recommended structure for the semi-arid areas in this paper by comparisons with observed streamflow: (1) a dataset of 53 watersheds of southeastern Australia was used to identify the appropriate ranges of watershed size and annual precipitation quantity for the HIMS model.(2) The other dataset of the Lhasa River watershed in China was used to test the HIMS model performance in estimating flood events at hourly time steps, the streamflow of time series at daily and monthly time steps, and maintaining water balance at annual time step.
(3) The simulated streamflow of the Lhasa River watershed by the HIMS model was compared with those by the VIC and SWAT models at daily and monthly time steps.Parameters adjusted in the calibration process are listed in Table 2 together with definitions and the changing ranges.They were identified by Liu et al. [29] with the GLUE method as highly sensitive parameters in influencing the estimation of actual ET, overland runoff, lateral flow, groundwater recharge, baseflow, sub-watershed flow, and channel watershed flow.Statistical parameters used to evaluate the model performance in estimating streamflow are the Nash-Sutcliffe efficiency (NSE) and the water error (WE).If the calculated WE is within the range of ±20%, the simulated results are regarded as being satisfactory.The commonly used parameter of standard deviation of measured data (RSR) and percentage of bias (PBIAS) referred to in Moriasi et al. [30] are used for model evaluation.The model performance is satisfactory when RSR is smaller than 0.6, and PBIAS is in the range of ±25% for the streamflow simulation.
Water 2018, 10, 962 where Y obs i is the ith observation, Y sim i is the ith simulated value, Y mean is the mean of the observed data, n is the number of the data.4) and elevation from 3550m to 7126 m.This watershed is in the continental climate zone with average air temperature of 7 °C, highest air temperature of 12 °C, and lowest air temperature of -2.5 °C.Annual averaged precipitation is around 500 mm with the wet

HIMS Evaluation with Observed Streamflow of the 53 Australia Watersheds
The HIMS model is run with the recommended structure of the semi-arid area for each of the 53 watersheds with the calibrated model parameters listed in Table 2.The streamflow is well simulated for the 53 watersheds after calibration indicated in Table 4 with the averaged NSE above 0.5 and WE

HIMS Evaluation with Observed Streamflow of the 53 Australia Watersheds
The HIMS model is run with the recommended structure of the semi-arid area for each of the 53 watersheds with the calibrated model parameters listed in Table 2.The streamflow is well simulated for the 53 watersheds after calibration indicated in Table 4 with the averaged NSE above 0.5 and WE less than 10% in both calibration and validation periods of the 53 watersheds.The statistical metrics are shown in Table 4 for watershed 401210 which is selected randomly from the 53 watersheds together with the adjusted parameters in Table 2.The daily streamflow comparison is shown in Figure 5 between the simulation and observation on the streamflow.Both the statistical metrics and the comparison of the time series indicate that the streamflow of watershed 401210 is well simulated by the HIMS model with the structure for the semi-arid areas.are shown in Table 4 for watershed 401210 which is selected randomly from the 53 watersheds together with the adjusted parameters in Table 2.The daily streamflow comparison is shown in Figure 5 between the simulation and observation on the streamflow.Both the statistical metrics and the comparison of the time series indicate that the streamflow of watershed 401210 is well simulated by the HIMS model with the structure for the semi-arid areas.For each of the 53 watersheds, NSE and WE are obtained through streamflow comparisons between simulation and observation in the calibration and validation periods.The relationships between NSE/WE and watershed sizes are explored in Figure 6a, Figure 6b respectively in order to identify the appropriate range of watershed size for the HIMS model utilization with ensured simulation accuracy.The watersheds of 400-600 km 2 are identified as the most appropriate ones for the utilization of the HIMS model with the present structure, as the streamflow of these watersheds is well simulated with the NSEs reaching the peaks and the WEs touching the bottoms (Figure 6a,b, Table 5).For the watersheds larger than 600 km 2 , the model performance is not good with high WE, as one set of HIMS parameters cannot represent the spatial variability of the large watershed on the stream networks, landscape, soil types, and land uses.This in turn increases the complexity of the hydrological processes and the uncertainty of the model prediction.For the watersheds smaller than 400 km 2 , the HIMS model performs decently instead of well due to the strong nonlinear relationship of rainfall-runoff which is different from the large area watersheds which have a simpler rainfallrunoff relationship by integrating the small-scale complexity.
The relationships between NSE/WE and precipitation quantity are explored in Figure 6c, Figure 6d respectively in order to clarify the ranges of annual precipitation quantity which are appropriate for the utilization of the HIMS model with a semi-arid structure.The watersheds with annual precipitation quantity exceeding 800 mm are more appropriate for the HIMS model as it has better predicting performance for these watersheds with greater NSE and lower WE than the watersheds with precipitation out of this range (Figure 6c,d).For the arid watersheds with precipitation less than 800 mm, the poor model performance of HIMS is possibly attributed to the high nonlinearity and heterogeneity of the rainfall-runoff process in these regions [29,32,33].The model performance is greater with the increase of size and precipitation quantity within the ranges of the watershed area and annual precipitation (Table 5, Figure 6).For each of the 53 watersheds, NSE and WE are obtained through streamflow comparisons between simulation and observation in the calibration and validation periods.The relationships between NSE/WE and watershed sizes are explored in Figure 6a,b respectively in order to identify the appropriate range of watershed size for the HIMS model utilization with ensured simulation accuracy.The watersheds of 400-600 km 2 are identified as the most appropriate ones for the utilization of the HIMS model with the present structure, as the streamflow of these watersheds is well simulated with the NSEs reaching the peaks and the WEs touching the bottoms (Figure 6a,b, Table 5).For the watersheds larger than 600 km 2 , the model performance is not good with high WE, as one set of HIMS parameters cannot represent the spatial variability of the large watershed on the stream networks, landscape, soil types, and land uses.This in turn increases the complexity of the hydrological processes and the uncertainty of the model prediction.For the watersheds smaller than 400 km 2 , the HIMS model performs decently instead of well due to the strong nonlinear relationship of rainfall-runoff which is different from the large area watersheds which have a simpler rainfall-runoff relationship by integrating the small-scale complexity.
The relationships between NSE/WE and precipitation quantity are explored in Figure 6c,d respectively in order to clarify the ranges of annual precipitation quantity which are appropriate for the utilization of the HIMS model with a semi-arid structure.The watersheds with annual precipitation quantity exceeding 800 mm are more appropriate for the HIMS model as it has better predicting performance for these watersheds with greater NSE and lower WE than the watersheds Water 2018, 10, 962 10 of 17 with precipitation out of this range (Figure 6c,d).For the arid watersheds with precipitation less than 800 mm, the poor model performance of HIMS is possibly attributed to the high nonlinearity and heterogeneity of the rainfall-runoff process in these regions [29,32,33].The model performance is greater with the increase of size and precipitation quantity within the ranges of the watershed area and annual precipitation (Table 5, Figure 6).

HIMS Evaluation with the Observed Streamflow of the Lhasa River Watershed in China
The performance of the HIMS model with the semi-arid structure is evaluated on streamflow estimation from two perspectives.The first one is on the simulated streamflow comparisons between the simulation and observations at hourly, daily, monthly time steps of the Lhasa River basin.The second one is on the water balance analysis with precipitation, actual ET, streamflow, and soil moisture changes at the annual time step.simulated streamflow at hourly time step is compared with the observed flood events between 2000 and 2012 with the peak flow, flood amount, and lag time listed for the entire watershed (Table 6).They are compared with the measurements in corresponding time periods with the simulation error within ±20%, and NSE above 0.8 for most of the flood simulation events.The simulated flood amounts are close to that of the measurements (Figure 7).The hydrograph shows floods with sharp peaks around July 2003 and gentle ones around September 2001, as the Lhasa River watershed locates in the arid and semi-arid area, where runoff is generated with the two mechanisms existing simultaneously in the watershed.Surface flow is the major component of the floods generated with the intensity excess mechanism, therefore, the times for flow concentrations are short with high velocity and the shape of the hydrograph is sharp and symmetrical [26].For floods generated with the saturation excess mechanism, baseflow is the major component of the streamflow, therefore, the channel routing time is relatively long and the shape of the hydrograph is gentle.The HIMS model performs better in simulating streamflow than the TOPMODEL and Xinanjiang model

HIMS Evaluation with the Observed Streamflow of the Lhasa River Watershed in China
The performance of the HIMS model with the semi-arid structure is evaluated on streamflow estimation from two perspectives.The first one is on the simulated streamflow comparisons between the simulation and observations at hourly, daily, monthly time steps of the Lhasa River basin.The second one is on the water balance analysis with precipitation, actual ET, streamflow, and soil moisture changes at the annual time step.
The simulated streamflow at hourly time step is compared with the observed flood events between 2000 and 2012 with the peak flow, flood amount, and lag time listed for the entire watershed (Table 6).They are compared with the measurements in corresponding time periods with the simulation error within ±20%, and NSE above 0.8 for most of the flood simulation events.The simulated flood amounts are close to that of the measurements (Figure 7).The hydrograph shows floods with sharp peaks around July 2003 and gentle ones around September 2001, as the Lhasa River watershed locates in the arid and semi-arid area, where runoff is generated with the two mechanisms existing simultaneously Water 2018, 10, 962 11 of 17 in the watershed.Surface flow is the major component of the floods generated with the intensity excess mechanism, therefore, the times for flow concentrations are short with high velocity and the shape of the hydrograph is sharp and symmetrical [26].For floods generated with the saturation excess mechanism, baseflow is the major component of the streamflow, therefore, the channel routing time is relatively long and the shape of the hydrograph is gentle.The HIMS model performs better in simulating streamflow than the TOPMODEL and Xinanjiang model of Peng et al. [34] in the overlaid study period from 2000 to 2001 with the correlation coefficient (R 2 ) 0.71 and 0.75, and WE −5% and −3%. of Peng et al. [34] in the overlaid study period from 2000 to 2001 with the correlation coefficient (R 2 ) 0.71 and 0.75, and WE −5% and −3%.8 on the simulated streamflow against observational streamflow at daily and monthly time steps in the periods of model calibration and validation.The peaks of the observed daily and monthly streamflow are not well reached by the simulated values as in reality the streamflow quantity is the target in the calibrating process instead of the peaks of the streamflow at the daily and monthly time steps based on the model function.The evaluating parameters are statistically satisfied with NSEs exceeding 0.80 and WEs within ±10% at the daily and monthly time steps in the calibration period and WE within ±20% in the validation period (Table 7).The simulated monthly streamflow is a summary of the daily streamflow for each month, therefore, it is reasonable that the NSE is higher and the WE is lower for the monthly streamflow than that at the daily step.Another two statistical parameters are calculated as well for model performance evaluation: RSR and PBIAS recommended by Moriasi et al. [30].RSR is less than 0.6 and PBIAS is within the range of ±25%, therefore, the model performance in the streamflow estimation is satisfied at both daily and monthly time steps.
Water 2018, 10, x FOR PEER REVIEW 12 of 17 evaluating parameters are statistically satisfied with NSEs exceeding 0.80 and WEs within ±10% at the daily and monthly time steps in the calibration period and WE within ±20% in the validation period (Table 7).The simulated monthly streamflow is a summary of the daily streamflow for each month, therefore, it is reasonable that the NSE is higher and the WE is lower for the monthly streamflow than that at the daily step.Another two statistical parameters are calculated as well for model performance evaluation: RSR and PBIAS recommended by Moriasi et al. [30].RSR is less than 0.6 and PBIAS is within the range of ±25%, therefore, the model performance in the streamflow estimation is satisfied at both daily and monthly time steps.To make sure the outputs of the HIMS model besides the streamflow are reasonably simulated, water balance analysis is conducted with elements including the actual ET and soil moisture changes of the entire Lhasa River watershed (Table 8).The streamflow of the entire watershed is the sum of the surface flow, lateral flow, and baseflow.The unit of the streamflow is transferred from cubic meter per second into millimeter through dividing by the watershed area for the convenience of indicating the water balance situation of each year.In the periods of calibration and validation, precipitation is almost equal to the sum of streamflow, actual ET, and soil moisture change.The streamflow fluctuates at the same step with precipitation, which is the dominant influencing factor of the streamflow generation and concentration processes.The soil moisture changes in the opposite direction to that of the actual ET.Therefore, the water balance is in the closed status on four hydrological elements: precipitation, actual ET, streamflow, and soil moisture change, in both the calibration and validation periods.To make sure the outputs of the HIMS model besides the streamflow are reasonably simulated, water balance analysis is conducted with elements including the actual ET and soil moisture changes of the entire Lhasa River watershed (Table 8).The streamflow of the entire watershed is the sum of the surface flow, lateral flow, and baseflow.The unit of the streamflow is transferred from cubic meter per second into millimeter through dividing by the watershed area for the convenience of indicating the water balance situation of each year.In the periods of calibration and validation, precipitation is almost equal to the sum of streamflow, actual ET, and soil moisture change.The streamflow fluctuates at the same step with precipitation, which is the dominant influencing factor of the streamflow generation and concentration processes.The soil moisture changes in the opposite direction to that of the actual ET.Therefore, the water balance is in the closed status on four hydrological elements: precipitation, actual ET, streamflow, and soil moisture change, in both the calibration and validation periods.

Simulated Streamflow Comparison between HIMS, VIC, and SWAT Models
The simulated streamflow of the Lhasa River basin by the HIMS model is compared in this study with that of the SWAT published by Dawaciren [35] and VIC models by Liu et al. [36] with three evaluating statistical parameters at both daily and monthly time steps.Both the HIMS and VIC models perform better than the SWAT model with NSE and R 2 greater than 0.85 at daily and monthly time steps (Table 9).It should be attributed to the runoff generation equation of the HIMS and VIC model that they performed better than the SCS-CN method in the SWAT model.However, the HIMS model has a similar performance to the VIC model with NSE, R 2 , and WE in both calibration and validation periods at daily and monthly time steps.The similarity of the simulation results between the HIMS and VIC model is attributed to the combined runoff generation mechanisms of the two models.The VIC model is the improved version which dynamically simulates the combined runoff generation mechanisms of the saturation excess and intensity excess with the time compression analysis method.The HIMS model has structurally combined the two mechanisms with three soil layers.The peaks of the observed streamflow are not reached and the bottom is underestimated by the HIMS simulated streamflow at both daily and monthly time steps (Figure 8).It also occurs in the streamflow simulation of the SWAT model and VIC model [35,36].This is probably due to the non-linear relationship of rainfall-runoff, which has not been clearly presented with equations in the academic field, therefore, the simulated streamflow error occurs in the three hydrological models.In addition, the "peak" and "bottom" are not reached by the simulated streamflow as the evaluating model needs to simulate the large proportion of the simulated rainfall-runoff events of a long time series containing the wet and dry periods, while the extreme events cannot be simulated properly by the model with one set of model parameters.Therefore, the HIMS with the recommended structure performs better than the SWAT model and has the similar performance to the VIC model of the improved version, even though all of the three hydrological models do not simulate the streamflow of "peak" and "bottom" point very well.

Discussion
It is important to have the hydrological model customized according to the hydrological characteristics of the study areas.The HIMS model provides multiple choices describing each of the water cycling processes, accounts for the combined runoff generation mechanisms, and adds a runoff generation model like LCM [19,20] as well as a newly explored solving method of channel routing model [23,24].Even though the existing models like the MIKE-SHE [37][38][39] and HEC-HMS [40,41] modularize the water cycling processes and technologies like the open modeling interface (OpenMI) technology [42][43][44][45] provide an environment of integrating the hydrological modules, the advantages of HIMS are not covered.Besides the characteristic of flexible structure, the HIMS model is suitable for un-gauged areas and areas of various channel routing scenarios with the friendly user interface and theoretically improved runoff generation and concentration models.
With the recommended structure of the HIMS model for semi-arid areas, the most appropriate area range is from 400 to 600 km 2 , and the appropriate annual precipitation is greater than 800 mm with ensured model performance on NSE and WE.This finding is close to that obtained by Parajka et al. [46,47] whose runoff-hydrograph predictions tend to be more accurate in humid than in arid catchments and more accurate in large than in small catchments based on the meta-analysis of 34 reports of 3874 catchments with over 20 hydrological models and detailed analysis of some individual basins among the 34 reports.If the area and precipitation are out of these ranges, the performance of the HIMS model cannot achieve high predicting accuracy, as the HIMS model with one set of parameters cannot represent the spatial variability of large areas and the strong non-linear rainfall-runoff relationship of the small and arid watersheds.
The HIMS model has been proved to be reliable in estimating streamflow at hourly, daily, monthly time steps by comparisons between simulation and observation and analysis of the water balance in the Lhasa River basin.The performance of the HIMS model was evaluated by comparing the simulated streamflow with the streamflow simulated by the widely used models VIC and SWAT.Though HIMS and VIC have a better performance than the SWAT model due to the combined runoff generation mechanism, HIMS is structurally advanced due to the modular of each hydrological process with multiple choices.All three models do not perform well in simulating the "peaks" and "bottoms" imbedded in the streamflow of the long time series at daily time step.More attention was focused on common precipitation events in the calibrating process to ensure decent statistical evaluating results, instead of the extreme storm events.However, HIMS needs to be improved in the non-point source pollution fields by adding in the water quality modules relative to the SWAT model with the modules of N/P/pesticides.

Conclusions
It is important to simulate the streamflow with the customized hydrological model according to the distinctive hydrological characteristics of the study areas, due to the influence of heterogeneity on the water cycling processes.However, most existing hydrological models cannot be customized to the study areas due to their fixed structures and modes.To solve this problem, this study presents the HIMS model with flexible structures and evaluates the HIMS model with the recommended structure for semi-arid areas.
The HIMS model is characterized by the modules of each hydrological process and the multiple sub-models in each modular.It structurally combines the two runoff generation mechanisms of infiltration excess and saturation excess with three soil layers.Particularly, it is suitable for un-gauged areas due to the addition of the LCM model in the runoff generation modular, and areas with various channel scenarios due to the new solving method of the Kinematic wave equation in the channel routing modular.The HIMS model has proved to be reliable with the recommended structure for semi-arid areas based on evaluation with two observed streamflow datasets.With the first dataset of 53 Australian watersheds, the appropriate watersheds were identified for the utilization of the HIMS model with area of 400-600 km 2 and annual precipitation over 800 mm.With the second dataset of the Lhasa River basin in China, the model performs decently in streamflow simulation based on the satisfactory statistical results with NSE higher than 0.8 and WE within ±20% on comparisons between simulation and observation on the streamflow at hourly, daily, and monthly time steps.In addition, the HIMS model is reliable in simulating the water cycling process, as the water balance is mostly closed at the annual time step on precipitation, streamflow, actual ET, and soil moisture change.Finally, HIMS has proved to be more advanced relative to the widely used VIC and SWAT models by comparisons on the simulated streamflow, as HIMS has a better performance than the SWAT model and a flexible structure which the VIC model does not have.
As an integrated hydrological modeling platform, HIMS was developed for water related researches and water resource managers to estimate flood prediction and streamflow estimation.The HIMS model was combined with graphical software and visualized with Geographical Information System (GIS), which realizes the interactions between the modeling simulation process and result visualization.The code of the HIMS is available on the official website of the Chinese Academy of Sciences.HIMS will become a tool for estimating nutrient losses with high accuracy by imbedding the modules of soil erosion and nutrient losses to estimate the spatial and temporal distribution of soil loss and loss of N/P and pesticides of watersheds.Additionally, experiments on basic theory are also needed to be conducted, as the HIMS model does not take into account the impact of the initial water content of soil layers on the storm events for flood prediction.This has been found to be non-negligible by analysis of the observational datasets used for this analysis.

Figure 2 .
Figure 2. Structure of the HIMS model.

Figure 2 .
Figure 2. Structure of the HIMS model.

Dataset 1
is used to identify the most appropriate size of sub-watershed for the HIMS model with the recommended structure for the semi-arid area.It is from 53 watersheds of the Murray Darling Basin (MDB) located in southern Australia.Figure3shows the spatial distribution of the 53 watersheds with the locations of weather stations and hydrological stations.The 53 watersheds have areas from 52 to 1251 m 2 (Table3), elevations from 40.45 to 1064.26 m, and forest cover percentages from 0 to 100%.They are located at 36-39 • S, 141-150 • E and belong to the humid zone, semi-humid zone, and semi-arid zone of climate type with annual precipitation of 514-1612 mm, air temperature of 8.19-14.69• C, and streamflow of 41-540 mm.The observed dataset comes from Commonwealth Scientific and Industrial Research Organization (CSIRO) land and water, Australia, and includes daily precipitation, maximum air temperature, minimum air temperature, daily streamflow from 1972 to 1990, and the soil data.DEM and land use maps are downloaded from the seamless United States Geological Survey (USGS) data server with resolution of 30 m × 30 m. Table 3 shows the watersheds distribution on watershed area and annual precipitation.The observed data of streamflow is used for model warming up from 1972 to 1973, calibration from 1974 to 1980, and validation from 1981 to 1990.The hydrological model setting for the 53 watersheds is the LCM model for runoff generation and the Kinematic wave equation for channel routing.Water 2018, 10, x FOR PEER REVIEW 7 of 17 2.3.Observational Datasets Dataset 1 is used to identify the most appropriate size of sub-watershed for the HIMS model with the recommended structure for the semi-arid area.It is from 53 watersheds of the Murray Darling Basin (MDB) located in southern Australia.Figure 3 shows the spatial distribution of the 53 watersheds with the locations of weather stations and hydrological stations.The 53 watersheds have areas from 52 to 1251 m 2 (Table 3), elevations from 40.45 to 1064.26 m, and forest cover percentages from 0 to 100%.They are located at 36-39° S, 141-150° E and belong to the humid zone, semi-humid zone, and semi-arid zone of climate type with annual precipitation of 514-1612 mm, air temperature of 8.19-14.69°C, and streamflow of 41-540 mm.The observed dataset comes from Commonwealth Scientific and Industrial Research Organization (CSIRO) land and water, Australia, and includes daily precipitation, maximum air temperature, minimum air temperature, daily streamflow from 1972 to 1990, and the soil data.DEM and land use maps are downloaded from the seamless United States Geological Survey (USGS) data server with resolution of 30 m × 30 m. Table 3 shows the watersheds distribution on watershed area and annual precipitation.The observed data of streamflow is used for model warming up from 1972 to 1973, calibration from 1974 to 1980, and validation from 1981 to 1990.The hydrological model setting for the 53 watersheds is the LCM model for runoff generation and the Kinematic wave equation for channel routing.

Figure 3
Figure 3 Spatial distribution of 53 watersheds in the southwestern area of Australia.

Figure 3 .
Figure 3. Spatial distribution of 53 watersheds in the southwestern area of Australia.

Dataset 2
is of the Lhasa River watershed (90 • 05 -93 • 20 E, 29 • 20 -31• 15 N) located in south Tibet with area of 31,109 km 2 (Figure4) and elevation from 3550m to 7126 m.This watershed is in the continental climate zone with average air temperature of 7 • C, highest air temperature of 12 • C, and lowest air temperature of -2.5 • C. Annual averaged precipitation is around 500 mm with the wet season from June to September accounting for 80% of the annual precipitation.The watershed is covered by the grassland, unused land, and woodland about 68%, 23%, and 9% correspondingly.The major soil types are Felty soils, Frigid calcic soils, and Brown cold calcic soils.There are three weather stations within the watershed and seven outside the watershed which are used to intercept the weather data for the entire watershed with the inverse distance weighted (IDW) method.The climate datasets of the 10 stations are obtained from the information center of the Chinese meteorological bureau, and include precipitation, maximum, and minimum air temperature from 1973 to 2008 at daily time steps.DEM data comes from the Global land cover facility, University of Maryland with resolution of 90 m × 90 m.The land cover map is from the Data Center of Environmental Sciences of Chinese Academies with resolution of 1 km × 1 km.The soil dataset is from the same organization with the same resolution.The Lhasa River watershed is divided into 15 sub-watersheds with average area of 2173 km 2 based on the utilization of the hydrological modular of ArcGIS 10.0.The dataset of the observed streamflow comes from the Lhasa River station which is a controlled hydrological station for the entire watershed, as the streamflow from the sub-watersheds concentrates into channels and reaches the watershed outlet after going through this station (Figure4).Model calibration and validation are done with the observed streamflow from 1976 to 1995 and 1996 to 2008 respectively.The Lhasa River watershed has one hydrological station (Lhasa station) and three meteorological stations (Lhasa, Damxung, Maizhokunggar) located in the watershed.Therefore, the hydrological model setting is as follows: the runoff is estimated with the LCM model which has good performance in the ungauged areas and the channel routing estimation is done with the Muskingum model recommended by Bai et al.[31] as an appropriate method for the arid and semi-arid areas.Water 2018, 10, x FOR PEER REVIEW 8 of 17 season from June to September accounting for 80% of the annual precipitation.The watershed is covered by the grassland, unused land, and woodland about 68%, 23%, and 9% correspondingly.The major soil types are Felty soils, Frigid calcic soils, and Brown cold calcic soils.There are three weather stations within the watershed and seven outside the watershed which are used to intercept the weather data for the entire watershed with the inverse distance weighted (IDW) method.The climate datasets of the 10 stations are obtained from the information center of the Chinese meteorological bureau, and include precipitation, maximum, and minimum air temperature from 1973 to 2008 at daily time steps.DEM data comes from the Global land cover facility, University of Maryland with resolution of 90 m × 90 m.The land cover map is from the Data Center of Environmental Sciences of Chinese Academies with resolution of 1 km × 1 km.The soil dataset is from the same organization with the same resolution.The Lhasa River watershed is divided into 15 sub-watersheds with average area of 2173 km 2 based on the utilization of the hydrological modular of ArcGIS 10.0.The dataset of the observed streamflow comes from the Lhasa River station which is a controlled hydrological station for the entire watershed, as the streamflow from the sub-watersheds concentrates into channels and reaches the watershed outlet after going through this station (Figure 4).Model calibration and validation are done with the observed streamflow from 1976 to 1995 and 1996 to 2008 respectively.The Lhasa River watershed has one hydrological station (Lhasa station) and three meteorological stations (Lhasa, Damxung, Maizhokunggar) located in the watershed.Therefore, the hydrological model setting is as follows: the runoff is estimated with the LCM model which has good performance in the ungauged areas and the channel routing estimation is done with the Muskingum model recommended by Bai et al. [31] as an appropriate method for the arid and semi-arid areas.

Figure 4 .
Figure 4. Lhasa River watershed with the geographical locations of the hydrological and meteorological stations in the study area.

Figure 4 .
Figure 4. Lhasa River watershed with the geographical locations of the hydrological and meteorological stations in the study area.

Figure 5 .
Figure 5. Daily streamflow comparison between simulation and observation for the watershed 401210.

Figure 5 .
Figure 5. Daily streamflow comparison between simulation and observation for the watershed 401210.

17 Figure 6 .
Figure 6.Changes of NSE and WE with the watershed sizes (top) and annual precipitation (bottom).

Figure 6 .
Figure 6.Changes of NSE and WE with the watershed sizes (top) and annual precipitation (bottom).

Figure 7
Figure 7 Comparisons between simulated and measured streamflow at hourly time step.(a) Strom events simulation in 2000; (b) Strom events simulation in 2001; (c) Strom events simulation in 2003; (d) Strom events simulation in 2010.Streamflow comparisons are shown in Figure 8 on the simulated streamflow against observational streamflow at daily and monthly time steps in the periods of model calibration and validation.The peaks of the observed daily and monthly streamflow are not well reached by the simulated values as in reality the streamflow quantity is the target in the calibrating process instead of the peaks of the streamflow at the daily and monthly time steps based on the model function.The

Figure 8 .
Figure 8. Comparisons between simulated and measured streamflow at daily and monthly time steps.

Figure 8 .
Figure 8. Comparisons between simulated and measured streamflow at daily and monthly time steps.

Table 1 .
Main equations of the customization HIMS model.

Table 1 .
Main equations of the customization HIMS model.

Table 2 .
Parameters of HIMS for adjustment and values in the calibration period.

Table 3 .
Watersheds distribution on watershed area and annual precipitation.

Table 3 .
Watersheds distribution on watershed area and annual precipitation.
NOTE: -denotes no watershed locates in that range.

Table 4 .
The average statistical parameter values for 53 watersheds in calibration and validation.
Note: std* denotes the standard deviation.

Table 4 .
The average statistical parameter values for 53 watersheds in calibration and validation.
Note: std* denotes the standard deviation.

Table 5 .
Model performance with NSE on streamflow prediction for watersheds with various sizes and precipitation quantities.

Table 5 .
Model performance with NSE on streamflow prediction for watersheds with various sizes and precipitation quantities.

Table 6 .
Comparisons between the measurements and simulation on flood events from 2000 to 2012.

Table 6 .
Comparisons between the measurements and simulation on flood events from 2000 to 2012.

Table 7 .
Evaluation results of precipitation of the time series.

Table 7 .
Evaluation results of precipitation of the time series.

Table 8 .
Water balance annual from 1996 to 2008 together with that of calibration and validation periods.

Table 9 .
Statistical comparison with the observed streamflow and the streamflow simulated with the HIMS, VIC, SWAT model at monthly time step.