Integrating Topography and Soil Properties for Spatial Soil Moisture Storage Modeling

The understanding of the temporal and spatial dynamics of soil moisture and hydraulic property of soil is crucial to the study of hydrological and ecological processes. The purpose of this study was to derive equations that describe spatial soil water storage deficit based on topography and soil properties. This storage deficit together with the topographical index can be used to conclude the spatial distribution curve of storage capacity in a (sub-) basin for developing hydrological model. The established model was able to match spatial and temporal variations of water balance components (i.e., soil moisture content (SMC), evapotranspiration, and runoff) over the Ziluoshan basin. Explicit expression of the soil moisture storage capacity (SMSC) in the model reduced parameters, which provides a method for hydrological simulation in ungauged basins.


Introduction
The understanding of the temporal and spatial dynamics of soil moisture and hydraulic property is crucial to the study of several hydrological and ecological processes.Soil moisture is a key variable in hydrological modeling.As McNamara et al. [1] pointed out that understanding soil storage and its role in regulating catchment functions should be a priority in future observation strategies and hydrological modelling.Soil moisture and hydrological routing is computed based on the saturated and unsaturated flow equations, which can be solved by finite element or difference techniques over a three dimensional grid.These models include System Hydrological European model (SHE) [2,3] and Institute of Hydrology Distributed Model (IHDM) [4,5].While the advantages of pure, numerical simulation would seem clear, the tremendous amount of parameter evaluation required is problematic.In most cases the available data motivates the use of simple, conceptual model approaches rather than the use of a fully distributed, physically model with a large number of model parameters [6].
It is widely recognized that for saturation excess overland flow, soil storage is one of the key elements controlling runoff production in catchments in humid temperate areas.For traditionally lumped hydrological models, soil moisture is usually computed on the basis of basin average routing of unsaturated water balance, and its spatial variation is expressed in an empirical distribution curve of soil moisture storage capacity, e.g., a single parabolic curve in Xin'anjiang hydrological model [7].The model usually lacks physical expressions and distribution functions of the hydrological processes or parameters, which solely calibrated based on observed flow discharge.The semi-distributed models are capable to overcome shortage of the empirical models since they explicitly express watershed distribution of soil moisture storage using site-specific information of land surface features, e.g., soil properties and topography.For example, the topographic index [8] is widely used to analytically describe such distribution of soil moisture and runoff generation [9].It is also used for improving the empirical hydrological models with a distributed function.Such as the Xin'anjiang model, the statistical curve of spatial distribution of storage capacity was directly derived from TOPMODEL's topographic index [10][11][12] and from catchment slopes [13] or the empirical shape parameter B of the water storage capacity of the soil non-linearly was estimated in terms of the characteristic land surface slope [14].Some signal from the topography and soil type was used to explain the soil moisture variation [15].These models are benefit from a combination of reasonable computation time and physically realistic hillslope simulation.
Explicitly description of storage capacity for computation of hydrological fluxes remains challenge as it relates to topography, soil, vegetation and base rock variations.In the mountainous areas, topography controls depth to groundwater table and thus distribution of soil moisture deficit.In the hilly lowland areas, near stream saturated zones will be most extensive in locations with concave hillslope profiles and wide flat valleys (Figure 1).The storage deficit in the unsaturated zone becomes small and the runoff is easily generated because the unsaturated zone for holding soil moisture is reduced due to the high groundwater table occupation or low active depth of the unsaturated zone [12].
Water 2017, 9, 647 2 of 17 The model usually lacks physical expressions and distribution functions of the hydrological processes or parameters, which solely calibrated based on observed flow discharge.The semi-distributed models are capable to overcome shortage of the empirical models since they explicitly express watershed distribution of soil moisture storage using site-specific information of land surface features, e.g., soil properties and topography.For example, the topographic index [8] is widely used to analytically describe such distribution of soil moisture and runoff generation [9].It is also used for improving the empirical hydrological models with a distributed function.Such as the Xin'anjiang model, the statistical curve of spatial distribution of storage capacity was directly derived from TOPMODEL's topographic index [10][11][12] and from catchment slopes [13] or the empirical shape parameter B of the water storage capacity of the soil non-linearly was estimated in terms of the characteristic land surface slope [14].Some signal from the topography and soil type was used to explain the soil moisture variation [15].These models are benefit from a combination of reasonable computation time and physically realistic hillslope simulation.
Explicitly description of storage capacity for computation of hydrological fluxes remains challenge as it relates to topography, soil, vegetation and base rock variations.In the mountainous areas, topography controls depth to groundwater table and thus distribution of soil moisture deficit.In the hilly lowland areas, near stream saturated zones will be most extensive in locations with concave hillslope profiles and wide flat valleys (Figure 1).The storage deficit in the unsaturated zone becomes small and the runoff is easily generated because the unsaturated zone for holding soil moisture is reduced due to the high groundwater table occupation or low active depth of the unsaturated zone [12].The available space for moisture storage in unsaturated zone depends on not only the active depth but also the rise of the capillary fringe (Figure 1).Surface tension and soil pore capillary force lead the groundwater rise up into the soil and form a capillary fringe with a certain thickness [16], which significantly decreases soil moisture deficit.The capillary rise is associated with physical properties of soil porosity and particle diameter.For clay, where capillary is strong, a water table depth of 10 m can still be "felt" in the root zone and near surface.For sand, the water table has little role as a source if it is below the root zone [17].It is found that low runoff potential for soils having the low rise of the capillary fringe, such as deep sand, deep loess, aggregated silts, and high runoff potential for soils having the high rise of the capillary fringe, such as heavy plastic clays, and certain saline soils [18].
The main objective of this study is to derive equations for description of the spatial soil storage capacity that can make use of topography and soil property in the humid hilly watershed.The The available space for moisture storage in unsaturated zone depends on not only the active depth but also the rise of the capillary fringe (Figure 1).Surface tension and soil pore capillary force lead the groundwater rise up into the soil and form a capillary fringe with a certain thickness [16], which significantly decreases soil moisture deficit.The capillary rise is associated with physical properties of soil porosity and particle diameter.For clay, where capillary is strong, a water table depth of 10 m can still be "felt" in the root zone and near surface.For sand, the water table has little role as a source if it is below the root zone [17].It is found that low runoff potential for soils having the low rise of the capillary fringe, such as deep sand, deep loess, aggregated silts, and high runoff potential for soils having the high rise of the capillary fringe, such as heavy plastic clays, and certain saline soils [18].
The main objective of this study is to derive equations for description of the spatial soil storage capacity that can make use of topography and soil property in the humid hilly watershed.The equations were derived by combining the van Genuchten model of water retention relationship [19] Water 2017, 9, 647 3 of 15 and the TOPMODEL topographic index.They can be used to describe the site specific soil storage capacity and its statistical feature in a (sub-) basin, thus to simulate spatial distribution of hydrological variables of runoff amount, evaporation and soil water storage.The model was tested on Ziluoshan basin of Huai River watershed in the humid region in China.

Description of Ziluoshan Basin
As one of the first tributaries in the upstream of Huai River, Shaying River originated in the western mountainous area of Henan province of China (Figure 2).The watershed area above the Ziluoshan hydrological station on Shaying River is 1800 km 2 , and the mountainous area accounts for 75%.These geographical and climatic features result in extremely uneven of the annual and seasonal distributions of rainfall.The annual precipitation during 1980-1996 is 900 mm, varying from the largest in the southeast to the smallest in the north.Most of annual precipitation is concentrated in the flood seasons (June-September).Nearly 60-70% of the total precipitation occurs in the months of June and August.
DEM data of 30 m grid resolution are used to describe the spatial variations of topography [20,21].Based on the geographical map, the terrain of the catchment tilts from southwest to northeast (Figure 2).The mean elevation of topography is 820 m, varying from 284 m above the mean sea level in the downstream region to over 2122 m in the upstream mountainous areas.The topographic index in each pixel was calculated using the DTM9704 program [22,23].
Water 2017, 9, 647 3 of 17 equations were derived by combining the van Genuchten model of water retention relationship [19] and the TOPMODEL topographic index.They can be used to describe the site specific soil storage capacity and its statistical feature in a (sub-) basin, thus to simulate spatial distribution of hydrological variables of runoff amount, evaporation and soil water storage.The model was tested on Ziluoshan basin of Huai River watershed in the humid region in China.

Description of Ziluoshan Basin
As one of the first tributaries in the upstream of Huai River, Shaying River originated in the western mountainous area of Henan province of China (Figure 2).The watershed area above the Ziluoshan hydrological station on Shaying River is 1800 km 2 , and the mountainous area accounts for 75%.These geographical and climatic features result in extremely uneven of the annual and seasonal distributions of rainfall.The annual precipitation during 1980-1996 is 900 mm, varying from the largest in the southeast to the smallest in the north.Most of annual precipitation is concentrated in the flood seasons (June-September).Nearly 60-70% of the total precipitation occurs in the months of June and August.
DEM data of 30 m grid resolution are used to describe the spatial variations of topography [20,21].Based on the geographical map, the terrain of the catchment tilts from southwest to northeast (Figure 2).The mean elevation of topography is 820 m, varying from 284 m above the mean sea level in the downstream region to over 2122 m in the upstream mountainous areas.The topographic index in each pixel was calculated using the DTM9704 program [22,23].Soil in the basin is primarily shallow eluvial-slope deposits consisting of sand, loam and clay (Data Center for Resources and Environmental Sciences Chinese Academy of Sciences (RESDC)).Percentage of sand and clay within the soil depth of 0-1.0 m ranges 5.3-83.0%for sand, 1.7-31.6%for clay (Figure 3), and the remaining for loam.Value of the hydraulic parameters for the three soils refers to Tuller [24] (Table 1).The mean values of soil hydraulic parameters of VG in each pixel were estimated by proportion-weighted arithmetic mean way, and the field capacity θf was determined as Soil in the basin is primarily shallow eluvial-slope deposits consisting of sand, loam and clay (Data Center for Resources and Environmental Sciences Chinese Academy of Sciences (RESDC)).Percentage of sand and clay within the soil depth of 0-1.0 m ranges 5.3-83.0%for sand, 1.7-31.6%for clay (Figure 3), and the remaining for loam.Value of the hydraulic parameters for the three soils refers to Tuller [24] (Table 1).The mean values of soil hydraulic parameters of VG in each pixel were estimated by proportion-weighted arithmetic mean way, and the field capacity θ f was determined as a fraction of the saturated SMC (75% in this study [25]).The vegetation in the study region is primarily deciduous broadleaf forest and mixed forest.The average vegetation coverage is larger than 75% [26] and average annual Normalized Difference Vegetation Index (NDVI) is larger than 0.6 [27].a fraction of the saturated SMC (75% in this study [25]).The vegetation in the study region is primarily deciduous broadleaf forest and mixed forest.The average vegetation coverage is larger than 75% [26] and average annual Normalized Difference Vegetation Index (NDVI) is larger than 0.6 [27].The watershed was divided into a number of sub-basins (e.g., 50 sub-basins in this study) for describing spatial variations of runoff generation and river flow routing (Figure . 4).Each sub-basin includes lots of pixels in 30 m grid resolution.

Outline of the Model
The soil moisture storage model was firstly introduced at site scale, and then extended to hillslope soil moisture storage capacity (SMSC) according to TOPMODEL concept.Finally, a spatial soil moisture storage curve was concluded by all sites SMSC in a statistical way to replace the traditional statistical curve in Xin'anjiang model.The detail procedure is given blow.

The role of Soil Moisture Storage on Hydrological Fluxes
Hydrological balance on any element can be expressed as following: where W is soil moisture storage at the time interval of t−1 and t, Pt, Et and Rt is precipitation, actual evapotranspiration and runoff at the time interval of t−1 and t, respectively.In Equation ( 1), computation of both actual evapotranspiration E and runoff R depends on soil storage W (SMC in this model) and its capacity Wm (SMSC in this model), e.g., for saturation excess The watershed was divided into a number of sub-basins (e.g., 50 sub-basins in this study) for describing spatial variations of runoff generation and river flow routing (Figure 4).Each sub-basin includes lots of pixels in 30 m grid resolution.

Outline of the Model
The soil moisture storage model was firstly introduced at site scale, and then extended to hillslope soil moisture storage capacity (SMSC) according to TOPMODEL concept.Finally, a spatial soil moisture storage curve was concluded by all sites SMSC in a statistical way to replace the traditional statistical curve in Xin'anjiang model.The detail procedure is given blow.

The Role of Soil Moisture Storage on Hydrological Fluxes
Hydrological balance on any element can be expressed as following: where W t−1 and W t is soil moisture storage at the time interval of t − 1 and t, P t , E t and R t is precipitation, actual evapotranspiration and runoff at the time interval of t − 1 and t, respectively.In Equation ( 1), computation of both actual evapotranspiration E and runoff R depends on soil storage W (SMC in this model) and its capacity Wm (SMSC in this model), e.g., for saturation excess Water 2017, 9, 647 5 of 15 overland flow, runoff is generated as the soil storage reaches the capacity, which can be expressed as following: where Pe is net rainfall.For distributed modeling or physical description of spatial variation function of watershed characteristics on hydrological fluxes, expression of soil spatial variation storage and its capacity in a watershed is vital.

Site Specific Soil Moisture Storage Capacity
SMSC is usually defined as the difference between the water content at field capacity and at residual multiplied by a critical depth in unsaturated zone for moisture storage and runoff generation.In humid regions, the critical depth can be regarded as unsaturated zone thickness or the depth to free water table if there is a free water table in the soil profile.In this case, vertical water content in the critical depth can be influenced by the free water table under the capillary flux from the free water to soil moisture column (Figure 1a), which results in that soil moisture seldom reach the residual water content on the groundwater table nearby.Thus, the soil moisture profile above the water table can be determined from the balance between the pressure head gradient, which tends to draw moisture up, and gravity [28].
The soil-water pressure head distribution, ψ(z), may be modeled using Darcy's law: or where q is the steady state evaporation rate, K(ψ) is the unsaturated hydraulic conductivity function, ψ is the pressure head (m).Soil moisture with depth can be described by the van Genuchten model (hereafter VG) of water retention relationship [19,29]: where ψ is the pressure head (m), θ s and θ r is saturated and residual volume water contents, respectively, α, n are fitting parameters.Solution of Equation (3) or Equation ( 4) coupling with Equation ( 5) is a rather complicated although analytical steady state solutions can be found in the literature for some functions different with VG [30,31].These mathematical difficulties are significantly reduced if the vertical profile of soil moisture is approximated with the profile corresponding to zero-flux conditions [32].Under the assumption of zero vertical flux (quasi steady-state hydraulic assumption), the relationship between ψ and vertical distance above groundwater table (z is taken here to increase upwards) is that of hydraulic equilibrium: where D is the depth to groundwater table.Some comparisons between the vertical profile of soil moisture obtained with Equation (3) and the approximated profile provided by Equation ( 6) have been made in many literatures.It was found that difference of the soil moisture profiles between the two equations increases when the depth z increases to the nearby surface and the soil texture becomes finer.However, the difference between the Water 2017, 9, 647 6 of 15 two profiles is smaller than 5-7% for a variety of soil textures [33].For shallow water tables, the steady profile is close to hydrostatic state subject to no atmospheric forcing [34].
For the steady profile, the total profile moisture deficit will be obtained by integrating from the surface to the top of the capillary fringe.Using Equations ( 5) and ( 6), we obtain a simple mathematical expression for the unsaturated zone moisture storage S f (D), given the depth to groundwater table D: where ψ c is suction pressure for the unsaturated zone storage at field capacity.Hence, the soil moisture deficit or SMSC W m (D) defined as the difference between the moisture content at field capacity θ f and at the unsaturated zone moisture profile S f (D) in Equation ( 7) is given by: Equation ( 8) explicitly expresses the storage capacity in terms of the site specific depth to groundwater table and hydraulic parameters of soil.A large storage capacity corresponds to low groundwater table and low capillary fringe.

Hillslope SMSC
In a hilly mountain area, the topographic index is widely used to represent the influences of terrain on the spatial variations of soil wetness.In TOPMODEL, the large topographic index always being obtained in the local valley area, where happened to be the saturated zone.So a larger topographic index in a local area means less soil moisture deficit or easier runoff generation in response to rainfall input [12].Figure 1b illustrates that the local moisture deficit varies significantly along the catchment slope, with low values where the water table is near the surface (at the bottom of hills) and high values where the water table is deeper (at the top of hills) [35].
According to TOPMODEL concept, the depth to groundwater table at any location is given in terms of the watershed average depth to the water table and local topography wetness index [8]: where D i is the depth to groundwater table at the location i, D is the watershed-average depth to the water table, S zm is a scaling parameter, known as the 'effective depth' that determines the decay of hydraulic conductivity with depth, tp i is local topography wetness index (tp i = ln( a tan β ), where a is the cumulative area draining through a unit length of contour line, β is the slope of the unit area), λ is the areal average of tp i .
Equation ( 9) is written in the form: Equation ( 10) is substituted into Equation ( 8) to get: Equation ( 11) explicitly expresses soil moisture deficit at any location influenced by integration of soil properties, topography and depth to groundwater table.
The local storage capacity W m (D i ) and its basin mean value WM solely depend on the effective soil depth S zm if other parameters in the van Genuchten model in Equation (11) are measured or given according to available investigations.Beven et al. [36] described that larger value of the parameter S zm means increase of active soil depth for water storage and runoff generation.Habets and Saulnier [37] stated that the parameter S zm can be defined as one quarter of the maximum storage deficit.
In some humid area like Ziluoshan basin, streams are typically gaining streams (gaining water by drainage of baseflow from the groundwater into the stream).For the perennial rivers, the groundwater table near the surface is coincident or close to the stream water surface elevation [18,38].Thus D i and W m (D i ) approach zero in the areas nearby the stream water surface where tp i reaches its maximum value.Then, coefficient C in Equation ( 10) is: Equation ( 12) is substituted into Equation (11) to get: Thus, spatial variation of topographic index and soil properties in Table 1 are used to compute the curves of W m versus f /F for all sub-basins as shown in Figure 4.
Equation ( 11) explicitly expresses soil moisture deficit at any location influenced by integration of soil properties, topography and depth to groundwater table.
The local storage capacity Wm (Di) and its basin mean value WM solely depend on the effective soil depth Szm if other parameters in the van Genuchten model in Equation (11) are measured or given according to available investigations.Beven et al. [36] described that larger value of the parameter Szm means increase of active soil depth for water storage and runoff generation.Habets and Saulnier [37] stated that the parameter Szm can be defined as one quarter of the maximum storage deficit.
In some humid area like Ziluoshan basin, streams are typically gaining streams (gaining water by drainage of baseflow from the groundwater into the stream).For the perennial rivers, the groundwater table near the surface is coincident or close to the stream water surface elevation [18,38].Thus Di and Wm (Di) approach zero in the areas nearby the stream water surface where tpi reaches its maximum value.Then, coefficient C in Equation ( 10) is: Equation ( 12) is substituted into Equation (11) to get: Thus, spatial variation of topographic index and soil properties in Table 1 are used to compute the curves of Wm versus f/F for all sub-basins as shown in Figure 4.

The Watershed Model Development
The model can be executed at each grid based on storage capacity W m (D i ) at any local area calculated in terms of Equation (11).For runoff generation, grid runoff depth for a specific precipitation amount can be calculated based on saturation overland flow concept that runoff occurs where soil moisture reaches storage capacity (Equation ( 2)).
The site-specific storage capacity can be grouped into a number of intervals, each representing that fraction of the watershed with similar water table depth and soil moisture characteristics [39].Water balance equations (Equation ( 2)) are applied at each interval, which separate estimates of unsaturated zone soil moisture, evapotranspiration and runoff for the whole intervals [40][41][42].
Water 2017, 9, 647 8 of 15 The individual fluxes in every interval are then really weighted and combined to obtain the water balance for the entire watershed.This modeling method is efficient in computation of flow routing and parameter calibration of the model.
For a large watershed, the model can be executed at each sub-basin based on distribution curve of W m ~f /F (f /F represents proportion of W m value in a specific interval with an area f to the total (sub-) basin area F, the ratio value is between 0 and 1.) and its average value WM in a (sub-) basin.This modeling method becomes lumped structure like Xin'anjiang model, but it physically interprets the field capacity curve that allows for application in ungauged basins and consecutive routing in the large watershed.
Computations of actual evaporation, separation of runoff components and flow routing in the watershed and rivers are same as Xin'anjiang model in this study.The calibrated parameter meaning was listed in Table 2.A spatial distribution curve of free water storage capacity with a catchment average SM and an exponent of the spatial distribution curve, EX, is used to represent regulation of catchment heterogeneity on free water R.For thin soils, SM is around 10 mm, and EX is between 1.0 and 1.5 [43].The free water R is then separated into overland flow R s , subsurface flow R i , and deep layer flow (groundwater) R g .KI and KG are the outflow coefficients of the free water storage to interflow and groundwater.It is suggested that the sum KI + KG may be taken as 0.7-0.8 and the ratio of the three runoff components will be changed by altering the ratio of KG/KI [43].For application of the model in a large watershed, the whole watershed can be divided into a lot of sub-basins.These sub-basins are linked with the river channel system.The flow hydrograph at a point on the watershed from a known hydrograph of upstream or sub-basin is routed by the Muskingum method.Small tributaries of sub-basins merge to form larger stream which ultimately lead to outlet of the watershed.

Model Verification
Within the study area, daily and hourly observed precipitation data at 7 rainfall stations, pan evaporation and streamflow at the outlet of watershed were used for model testing (Figure 2).The largest and smallest streamflow discharges during 1979 and 1995 are 2080.0and 1.0 m 3 /s, respectively.Daily and hourly rainfall in each of the watershed was interpolated by the thiessen polygon method from observation data of rainfall stations.The model is calibrated and validated based on daily and hourly observed streamflow discharges at the whole watershed outlet in calibration and validation periods, respectively.The calibration period is 10 years from 1979 to 1988, and the validation period is 7 years from 1989 to 1995.Twelve and seven flood events with hourly observation Water 2017, 9, 647 9 of 15 data in calibration and validation periods, respectively, are further selected for model calibration and validation.
The calibration procedures followed the Xin'anjiang model calibration proposed by Zhao et al. [43]: (1) giving initial values of the parameters suggested; (2) calibrating the parameters of runoff generation processes to test annual water balance, e.g., KC in Table 2, by comparing totals of the simulated and observed runoff; (3) calibrating other parameters of water source separation and flow routing in Table 2 to test flow discharges by comparing the simulated and observed runoff.Model performance is evaluated with respect to the objective function of maximizing Nash-Sutcliffe efficiency coefficient (NSC) between observed and simulated runoff.Another measure of model performance used in the study is the root mean squared error (RMSE) of flow discharge.The calibrated parameters and values were listed in Table 2.The calibrated value of S zm is 42.8 mm.It is used to calculate mean value of storage capacity WM in each sub-basin and in whole watershed.The sub-basin mean WM for the 50 sub-basins ranges 107-127 mm, and whole watershed WM is 120 mm.

Result and Discussion
The calibration and validation results of flow discharges for daily and hourly data are listed in Tables 3 and 4, respectively.For daily simulation during the calibration period (Table 3), mean NSC value is 0.73, annually ranging 0.47-0.91;mean RMSE is 1.6 mm, annually ranging 0.6-4.4mm.The daily simulation results during the validation period are even better than those in the calibration period, e.g., mean NSC value is 0.77, ranging 0.67-0.87;mean RMSE is 1.0 mm, ranging 0.6-1.4mm.The model captured the flood processes (Table 4 and Figure 5) as well.During the calibration and validation periods, respectively, mean NSC value of all flood events is 0.85 and 0.83; the relative error of the simulated and observed peak discharges is 3.70% and 2.26%; mean RMSE is 44.5 and 35.6 m 3 /s.For illustrative purpose, hourly simulated and observed flood discharges are shown in Figure 5.These results demonstrated that the model is capable of reproducing both the magnitude and the dynamics of the daily and hourly flow discharges.

Sensitivity Analysis
In this study, sensitivity analysis of parameters for the effective depth Szm and soil properties was carried out to evaluate and quantify the effect of the parameter variations on model output.Curves of Wm ~ f/F for Szm values within 10-50 mm in the whole watershed were shown in Figure 6a.Sensitive analysis was executed by simulating changes of streamflow in response to these curves.For the specific soil distribution in the study watershed (Figure 3

Sensitivity Analysis
In this study, sensitivity analysis of parameters for the effective depth zm and soil properties was carried out to evaluate and quantify the effect of the parameter variations on model output.Curves of W m ~f /F for S zm values within 10-50 mm in the whole watershed were shown in Figure 6a.Sensitive analysis was executed by simulating changes of streamflow in response to these curves.For the specific soil distribution in the study watershed (Figure 3  Sensitivity of soil properties on model output was executed for the rainfall.If the soils in the whole watershed were replaced by one of the three soils (silt, sandy loam and clay), the analysis results of relationship between WM and topographic index in 50 sub-basins indicated that sub-basin mean WM decreased with the topographic index increase (Figure 7a).Because of high residual moisture content θr and capillary rise for clay, the storage capacity WM for clay was much smaller than that of sand and sandy loam.Therefore, precipitation on the clay soil generates much more runoff than that on sand and sandy loam.As an example, simulated flood discharges for the precipitation of 81.9 mm during 15-18 September 1985 for sand, sandy loam and clay are shown in Figure 7b.This amount of precipitation generates runoff of 55.8 mm for clay, 42% and 52% larger than the generated runoff of 32.4 and 26.6 mm for sandy loam and sand, respectively.For the simulated flood events, the peak discharge generated on clay is 696.0 m 3 /s, 44% and 63% larger than the peak discharges of 386.0 and 260.0 m 3 /s on sandy loam and sand, respectively.Sensitivity of soil properties on model output was executed for the rainfall.If the soils in the whole watershed were replaced by one of the three soils (silt, sandy loam and clay), the analysis results of relationship between WM and topographic index in 50 sub-basins indicated that sub-basin mean WM decreased with the topographic index increase (Figure 7a).Because of high residual moisture content θ r and capillary rise for clay, the storage capacity WM for clay was much smaller than that of sand and sandy loam.Therefore, precipitation on the clay soil generates much more runoff than that on sand and sandy loam.As an example, simulated flood discharges for the precipitation of 81.9 mm Water 2017, 9, 647 12 of 15 during 15-18 September 1985 for sand, sandy loam and clay are shown in Figure 7b.This amount of precipitation generates runoff of 55.8 mm for clay, 42% and 52% larger than the generated runoff of 32.4 and 26.6 mm for sandy loam and sand, respectively.For the simulated flood events, the peak discharge generated on clay is 696.0 m 3 /s, 44% and 63% larger than the peak discharges of 386.0 and 260.0 m 3 /s on sandy loam and sand, respectively.whole watershed were replaced by one of the three soils (silt, sandy loam and clay), the analysis results of relationship between WM and topographic index in 50 sub-basins indicated that sub-basin mean WM decreased with the topographic index increase (Figure 7a).Because of high residual moisture content θr and capillary rise for clay, the storage capacity WM for clay was much smaller than that of sand and sandy loam.Therefore, precipitation on the clay soil generates much more runoff than that on sand and sandy loam.As an example, simulated flood discharges for the precipitation of 81.9 mm during 15-18 September 1985 for sand, sandy loam and clay are shown in Figure 7b.This amount of precipitation generates runoff of 55.8 mm for clay, 42% and 52% larger than the generated runoff of 32.4 and 26.6 mm for sandy loam and sand, respectively.For the simulated flood events, the peak discharge generated on clay is 696.0 m 3 /s, 44% and 63% larger than the peak discharges of 386.0 and 260.0 m 3 /s on sandy loam and sand, respectively.

Spatial Variation of Hydrological Components
For hydrological modeling, distribution of the three components, soil moisture content, actual evapotranspiration and runoff, are vital for water resources management, environmental protection and assessment of land use and climate change impacts on hydrology.A major strength of the model is its capability to describe the spatial variations of such components.
The developed model was executed in each of sub-basins with rainfall input and W m ~f /F curves (Figure 4).The hydrological components, i.e., soil moisture content, actual evapotranspiration and runoff were calculated in each sub-basin and their spatial variations were then described.For illustrative purposes, spatial variations of mean annual precipitation and runoff during 1980-1996 simulated by the model are shown in Figure 8a,b and soil moisture content and actual evapotranspiration on 1 October 1987 after a non-rainfall period of 20 days are shown in Figure 8c,d.

Spatial Variation of Hydrological Components
For hydrological modeling, distribution of the three components, soil moisture content, actual evapotranspiration and runoff, are vital for water resources management, environmental protection and assessment of land use and climate change impacts on hydrology.A major strength of the model is its capability to describe the spatial variations of such components.
The developed model was executed in each of sub-basins with rainfall input and Wm ~ f/F curves (Figure 4).The hydrological components, i.e., soil moisture content, actual evapotranspiration and runoff were calculated in each sub-basin and their spatial variations were then described.For illustrative purposes, spatial variations of mean annual precipitation and runoff during 1980-1996 simulated by the model are shown in Figure 8a,b and soil moisture content and actual evapotranspiration on 1 October 1987 after a non-rainfall period of 20 days are shown in Figure 8c,d.
Figure 8 demonstrates that for the multi-year average, runoff distribution (Figure 8b) agrees with precipitation distribution (Figure 8a), larger in south and southwest of the high mountain area and smaller in the north and north east.Simulated results in Figure 8a and Figure 8b indicate that actual evaporation distribution generally corresponds to soil moisture content distribution, and both components are controlled by storage capacity for a specific time.As shown in Figure 8c, soil moisture content on 1 October 1987 is larger in the middle areas of watershed with a high value of WM (mostly higher than 120 mm in Figure 4b), and smaller in the north and south areas of the watershed with a low value of WM (mostly lower than 120 mm in Figure 4b).Larger soil moisture offers more water for evaporation (Figure 8d).   Figure 8 demonstrates that for the multi-year average, runoff distribution (Figure 8b) agrees with precipitation distribution (Figure 8a), larger in south and southwest of the high mountain area and smaller in the north and north east.Simulated results in Figures 8a and 8b indicate that actual evaporation distribution generally corresponds to soil moisture content distribution, and both components are controlled by storage capacity for a specific time.As shown in Figure 8c, soil moisture content on 1 October 1987 is larger in the middle areas of watershed with a high value of WM (mostly higher than 120 mm in Figure 4b), and smaller in the north and south areas of the watershed with a low value of WM (mostly lower than 120 mm in Figure 4b).Larger soil moisture offers more water for evaporation (Figure 8d).

Conclusions
For the saturation excess overland flow modeling, the critical spatial soil moisture distribution was usually represented either by a grid by grid method or a distribution curve based on the controlling component, such as soil moisture (deficit).The latter method was generally more efficient, but such curve was traditionally empirical and lacked of directly physical interpretation, which did not allow for the extended application in ungauged basins.In order to overcome the limitations, the storage capacity distribution curve in a (sub-) basin was explicitly expressed by the van Genuchten model and topographic index in TOPMODEL on the basis of digital elevation and soil map.Because the new model had capacity in describing soil moisture deficit at a grid scale or its statistical curve in a sub-basin, the model could be used to simulate temporal and spatial distributions of hydrological components, and to investigate their relation with precipitation, topography and soil variations.
The model had been proven to be a reliable and flexible tool for rainfall-runoff modelling in Ziluoshan basin.The results demonstrated that the model was not only capable of reproducing daily and hourly flow discharges, but also able to produce reasonable spatial variations of hydrological components such as runoff, soil moisture storage and actual evapotranspiration.Sensitivity analysis showed effect of the soil effective depth and soil properties on runoff.This study helped us to understand scale effects of land surface characteristics on hydrological processes.

Figure 1 .
Figure 1.Sketch of vertical profile of soil moisture deficit Wm (D) at a specific site (a) and along hillslope (b).

Figure 1 .
Figure 1.Sketch of vertical profile of soil moisture deficit Wm(D) at a specific site (a) and along hillslope (b).

Figure 3 .
Figure 3. Percentage of sand and clay within the soil depth of 0-1.0 m, (a) Percentage of sand and (b) Percentage of clay.

Figure 3 .
Figure 3. Percentage of sand and clay within the soil depth of 0-1.0 m, (a) Percentage of sand and (b) Percentage of clay.

Figure 4 .
Figure 4. Spatial distribution of storage capacity WM (basin mean value of local storage capacity).(a) Cumulative frequency (cumulative proportion of partial basin area versus total basin area) of Wm ~ f/F in terms of Equation (13) for the 50 sub-basins, (b) sub-basin mean WM.

Figure 4 .
Figure 4. Spatial distribution of storage capacity WM (basin mean value of local storage capacity).(a) Cumulative frequency (cumulative proportion of partial basin area versus total basin area) of W m ~f /F in terms of Equation (13) for the 50 sub-basins, (b) sub-basin mean WM.

Figure 5 .
Figure 5. Hourly simulated and observed runoff in the watershed for the selected flooding events.
), change of mean annual streamflow for the flood period 11 May-23 September 1985 responded to changes in Szm and WM are shown in Figure 6b.Increase of Szm significantly increased WM and thus decreased runoff.Simulated results

Figure 5 .
Figure 5. Hourly simulated and observed runoff in the watershed for the selected flooding events.
), change of mean annual streamflow for the flood period 11 May-23 September 1985 responded to changes in S zm and WM are shown in Figure 6b.Increase of S zm significantly increased WM and thus decreased runoff.Simulated results indicate that as S zm increased from 10 to 50 mm, WM increased from 12.6 to 143.8 mm, resulting in 70.0%decrease of runoff from 251.5 to 75.4 mm.Water 2017, 9, 647 12 of 17 indicate that as Szm increased from 10 to 50 mm, WM increased from 12.6 to 143.8 mm, resulting in 70.0%decrease of runoff from 251.5 to 75.4 mm.

Figure 6 .
Figure 6.Distribution of soil moisture storage capacity Wm related to Szm ('effective depth' that determines the decay of hydraulic conductivity with depth) in the whole watershed (a) and sensitivity of parameter Szm to runoff R and watershed mean storage capacity WM (basin mean value of local storage capacity) for the flood period during 11 May-23 September, 1985 (b).

Figure 6 .
Figure 6.Distribution of soil moisture storage capacity W m related to S zm ('effective depth' that determines the decay of hydraulic conductivity with depth) in the whole watershed (a) and sensitivity of parameter S zm to runoff R and watershed mean storage capacity WM (basin mean value of local storage capacity) for the flood period during 11 May-23 September, 1985 (b).

Figure 7 .
Figure 7. Relationship between sub-basin mean WM (basin mean value of local storage capacity) and topographic index for 50 sub-basins with soil types of silt, sandy loam and clay, respectively.(Note: catchment mean Szm ('effective depth') that determines the decay of hydraulic conductivity with depth) is assumed as 42.8 mm in the figure.)(a) And simulated flood discharges for sand, sand loam and clay in 1985 (b).

Figure 7 .
Figure 7. Relationship between sub-basin mean WM (basin mean value of local storage capacity) and topographic index for 50 sub-basins with soil types of silt, sandy loam and clay, respectively.(Note: catchment mean S zm ('effective depth' that determines the decay of hydraulic conductivity with depth) is assumed as 42.8 mm in the figure).(a) And simulated flood discharges for sand, sand loam and clay in 1985 (b).

Figure 8 .
Figure 8. Spatial variation of mean annual precipitation (a) and runoff (b) during 1980-1996 and spatial variation of soil moisture content (c) and actual evaporation (d) on 1 October 1987.Figure 8. Spatial variation of mean annual precipitation (a) and runoff (b) during 1980-1996 and spatial variation of soil moisture content (c) and actual evaporation (d) on 1 October 1987.

Figure 8 .
Figure 8. Spatial variation of mean annual precipitation (a) and runoff (b) during 1980-1996 and spatial variation of soil moisture content (c) and actual evaporation (d) on 1 October 1987.Figure 8. Spatial variation of mean annual precipitation (a) and runoff (b) during 1980-1996 and spatial variation of soil moisture content (c) and actual evaporation (d) on 1 October 1987.
Note: for */*, the upper and lower values represent calibrated parameter values for daily and hourly simulations in the following verification, respectively.

Table 3 .
Results of model calibration and validation using daily data.

Table 4 .
Results of model calibration and validation using hourly data.