Subgrid Parameterization of the Soil Moisture Storage Capacity for a Distributed Rainfall-runoff Model

Spatial variability plays an important role in nonlinear hydrologic processes. Due to the limitation of computational efficiency and data resolution, subgrid variability is usually assumed to be uniform for most grid-based rainfall-runoff models, which leads to the scale-dependence of model performances. In this paper, the scale effect on the Grid-Xinanjiang model was examined. The bias of the estimation of precipitation, runoff, evapotranspiration and soil moisture at the different grid scales, along with the scale-dependence of the effective parameters, highlights the importance of well representing the subgrid variability. This paper presents a subgrid parameterization method to incorporate the subgrid variability of the soil storage capacity, which is a key variable that controls runoff generation and partitioning in the Grid-Xinanjiang model. In light of the similar spatial pattern and physical basis, the soil storage capacity is correlated with the topographic index, whose spatial distribution can more readily be measured. A beta distribution is introduced to represent the spatial distribution of the soil storage capacity within the grid. The results derived from the Yanduhe Basin show that the proposed subgrid parameterization method can effectively correct the watershed soil storage capacity curve. Compared to the original Grid-Xinanjiang model, the model performances are quite consistent at the different grid scales when the subgrid variability is incorporated. This subgrid parameterization method reduces the recalibration necessity when the Digital Elevation Model (DEM) resolution is changed. Moreover, it improves the potential for the application of the distributed model in the ungauged basin.


Introduction
For the last few decades, the development of numerous distributed rainfall-runoff models enables the spatial variations to be represented by a network of grid elements.Advances in geographic information systems (GIS), remote sensing (RS) and computational technology have also offered the potential to build complex distributed hydrologic models and improve the accuracy of hydrologic prediction in time and space [1][2][3].
However, many recent studies have suggested that distributed modelling approaches may not always provide improved simulations at the outlet compared to lumped conceptual models [4].One of the underlying reasons is that most models do not represent the spatial variability well [5,6].Limited by the resolution of available data and the computational efficiency, most grid-based distributed models do not take into account the subgrid variability of model input, parameters and model state [7][8][9].With the assumption of the uniform subgrid, the high frequency information of hydrologic variables and parameters will be lost as the large sampling dimensions act as a filter [10,11].
Because of the nonlinearity of hydrologic processes, model simulations of hydrologic responses are inherently sensitive to the spatial variability of input forcing and watershed characteristics [12].There is a growing awareness of the sensitivity of small scale variability in hydrologic systems modelling.For example, Wood et al. [13] introduced the concept of the Representative Elementary Area (REA).The REA is defined as a critical threshold at which the implicit continuum assumption can be used without knowledge of the patterns when building the catchment modelling.Arora et al. [14] compared the performances of land surface simulations with and without incorporating the subgrid variability of precipitation intensity and soil moisture.The results indicated that the inclusion of subgrid variability results in significant changes of magnitude, time and frequency of surface runoff generation and partitioning.Ghan et al. [15] assessed the relative influence of the subgrid variability of meteorology, vegetation characteristics, soil properties on surface runoff with a land surface model.They found that neglecting subgrid variability leads to the underestimation of runoff and overestimation of evaporation.Vázquez et al. [16] illustrated the scale-dependence of model effective parameters.They noted that the effective parameters need a recalibration when the grid resolution is changed.The scale-dependence of the model performances and effective parameters makes it a challenge to determine the scale at which the spatial variability can be replaced by grid-averaged characteristics.In order to reduce the necessity of recalibration, it is required to provide a well representation of subgrid variability.
The Grid-Xinanjiang model is an improved version of the Xinanjiang model, which is widely used for flood forecasting in humid and semi-humid area of China [17,18].The soil storage capacity is a key variable in the Grid-Xinanjiang model, which controls runoff generation and partitioning.Due to the relative short spatial structure, the soil storage capacity is more sensitive to the uniform grid assumption.The bias of estimation of soil storage capacity will introduce the uncertainty of hydrologic simulation, especially for application in the ungauged basin.Therefore, a good representation of subgrid variability of soil storage capacity is necessary for a scale-invariant hydrologic response in the Grid-Xinanjiang model.In this paper, based on the analysis of the model performances of the original Grid-Xinanjiang model at different grid scales, a subgrid parameterization method is presented to account for the subgrid variability of the soil storage capacity.
The rest of paper is organized as follows: in the next section, we provide a description of the study area; a brief description of the Grid-Xinanjiang model and the subgrid parameterization method are presented in the model description section; the following section includes a set of numerical experiments conducted in the Yanduhe Basin to compare the model performance with and without incorporating the subgrid variability; and the final section provides conclusions and perspectives.

Study Area
The Yanduhe Basin is located in central China with a catchment area of about 601 km 2 .The Yanduhe River, the main river of the Yanduhe Basin, originating from south of Shennongjia Mountain, flowing into the Yangtze River at 31°14′ N, 110°18′ E, is 60.6 km in length, with a mean slope of 9.5‰.More than 70% of the area is covered by vegetation.The watershed climate is humid, and the average annual precipitation is approximately 1300-1700 mm.Dominated by a monsoon climate, most rainfall occurs during the wet season between April and September.
As shown in Figure 1, the hourly rainfall is monitored at five rainfall gauging stations located in Duizi, Xiagu, Banqiao, Songziyuan and Yanduhe.Hourly streamflow and evaporation data are recorded at the outlet of the watershed.The morphology of the basin is described by the Digital Elevation Model (DEM) with 30 m resolution from the National Aeronautics and Space Administration (NASA).The elevation of the Yanduhe Basin ranges from 130 to 3031 m.

The Grid-Xinanjiang Model
DEMs are the primary computational elements for rainfall-runoff modeling in the Grid-Xinanjiang model.In each grid cell, canopy interception, direct channel precipitation, evapotranspiration, runoff generation and partitioning are simulated.A simple storage based equation is adopted to calculate the cumulative interception during rainfall events.The precipitation that falls on the channel can be treated as the direct runoff without loss.The evapotranspiration is simulated as a function of moisture content of the vertical soil profile.The soil profile is divided into three layers (upper, lower and deeper layers) in each grid cell to account for the uneven vertical distribution of the soil moisture content.The Grid-Xinanjiang model employs the saturation excess mechanism, and adopts the concept of soil storage capacity to partition the rainfall into runoff and infiltration.Depending on the free water capacity, the runoff is further subdivided into surface runoff, interflow and groundwater runoff.
For the daily simulation, the surface runoff directly routes to the channel as the fast component.The interflow and groundwater runoff represent the slow component that routes to the outlet of the corresponding subcatchment by two linear reservoirs with different lag times.The subcatchments are connected by the channel network.The outflow of each subcatchment is routed to the watershed outlet by the multiple-reach Muskingum method.More details are available in Yao et al. [17] and Liu et al. [19].

Subgrid Scale Parameterization of the Soil Storage Capacity
In the original Xinanjiang model, the spatially uneven distribution of the soil storage capacity is depicted by a power function [20]: where WM is the soil storage capacity, WMM is maximum watershed soil storage capacity, f/F is fraction of basin with the soil storage capacity less than WM, B is exponent of the curve.A similar method is also adopted in the ARNO [21] and the Variable Infiltration Capacity (VIC) models [22].The power function provides an analytical solution of runoff generation but lacks definite physical interpretation.Following the work of Williams et al. [23] and Anderson et al. [24], Yao et al. [25] utilized the soil texture and the land cover attributes to estimate the spatial distribution of the soil storage capacity.However, limited by resolution of soil property, the small-scale variability information is not readily available.Therefore, the uniform grid assumption is usually adopted.However, the high frequency information will be smoothed as the uniform grid acts as a filter.Especially for the soil storage capacity, due to the relative short spatial correlation structure, the spatial distribution of the soil storage capacity is more easily affected by the smoothing effect [26][27][28], which may further lead to the scale-dependence of the model performance.Therefore, it is necessary to take into account the subgrid variability of the soil storage capacity in the Grid-Xinanjiang model.
By analyzing the spatial organization of soil moisture and different terrain attributes (slope, mean curvature, topographic index, specific area, etc.), Western et al. [29] demonstrated that the topographic index (TI) can better explain the spatial variability at all scales from 10 m up to the catchment scales.
Shi et al. [30] also found a similar spatial pattern between the soil storage capacity and the topographic index.On one hand, topography is an important control of the soil moisture distribution [31,32].The areas that have large topographic indices usually correspond to the riparian areas, where it is easier to reach saturation because of the relatively small water deficit, and vice versa.On the other hand, the high resolutions of topographic indices are more readily obtained than the soil properties.Shi et al. [30] suggested the topographic index as an auxiliary variable to correlate with the soil storage capacity as following: where TI min is the minimum topographic index, a and b are two shape parameters.The topographic index TI is defined as ln(a/tan β), in which a is the cumulative area per unit length of contour line and β is the slope.As discussed by Pradhan et al. [33], a distinct shift of the value of the topographic index occurs as the DEM resolution becomes coarser.As a function of the topographic index, the soil storage capacity curve also changes correspondingly.Therefore, it is necessary to correct the grid-averaged soil storage capacity.By integrating the soil storage capacity in each grid, we can estimate the total watershed soil storage capacity at each grid scale.The deviation of the total watershed soil storage capacity between the finest grid and other grids can be used to correct the mean value of the soil storage capacity within the grid cell: where WMc is the soil storage capacity in each grid after correction, WMMt is the total watershed soil storage capacity at the finest resolution and WMMc is the total watershed soil storage capacity at the modelling scale.As suggested by Li and Avissar [34], here we use a beta distribution to describe the subgrid variability of the soil storage capacity.Compared to the other commonly used probability distribution (gamma, lognormal, etc.), the bounded character of beta distribution is more suitable for the description of soil property.It is defined as following: where Γ() is the gamma function, and α and β are two shape parameters.The parameters of the beta distribution can be estimated by the mean and variance values of the soil storage capacity within the grid.We can obtain the mean value of the soil storage capacity of the grid cell by Equations ( 2) and (3).The variance indicates the fluctuation range of the soil storage capacity within grid which can be derived by spatial aggregation from the finest resolution.Due to the uneven distribution of WM within the grid, part of the effective rainfall (Pe) is partitioned into the runoff (R), and the other part is stored in the soil (ΔW).As shown in Figure 2, the runoff R in a grid can be expressed as:

Model Validation
The proposed subgrid parameterization method is validated against the observed data from 1981 to 1986 using a daily time step.The model parameters are calibrated at 30 m resolution with the observations from 1981 to 1984.The rainfall data at five gauging stations is interpolated to each grid using the inverse distance weighting method.Adopting the seeding technology [35], the Yanduhe Basin is subdivided into 58 subcatchments with a threshold of 5 km 2 .We ran a drainage experiment presented by Vivoni [36] to obtain a reasonable spatial distribution of the initial soil moisture.The watershed is assumed to drain from a full saturation status without rainfall and evapotranspiration, until the watershed soil moisture content reaches the condition that we choose.The annual runoff deviation (ARD) and the Nash-Sutcliffe coefficient (NSC) are used as two criteria for the model evaluation.The ARD provides a general sense of the water balance, and the NSC is used to assess the goodness of fit between the simulated and the observed discharge.The ARD and NSC are calculated as follows: ( ) where sim R is the annual simulated runoff; obs R is the annual observed runoff; t sim Q is the simulated discharge at time step t ; t obs Q is the observed discharge at time step t ; and obs Q is the mean value of the observed discharge, n is the total number of values in the time period.
The simulated results at the Yanduhe Station are summarized in Table 1.The averaged ARD and NSC are 9% and 0.81 respectively in the calibration periods.In the validation period, the ARD and NSC also reach 10% and 0.76. Figure 3 shows an illustration of the simulated hydrograph against the observed data at the Yanduhe Station from April to October in 1986.The results demonstrate consistency between the simulation and observation.Further analysis is also based on the observed data in 1986.

The Scale-Dependence of the Grid-Xinanjiang Model with a Uniform Grid
The scale-dependence of model performance is the reason that we propose the subgrid parameterization method.Before evaluating the reasonability of the subgrid parameterization, we are also interested in, to what extent, the hydrologic responses of the original Grid-Xinanjiang model depend on the grid scale.Figure 4 depicts the probability density functions of the soil storage capacity from four different grid scales (50 m, 100 m, 500 m and 1000 m) in the Yanduhe Basin.A distinct shift is observed as the grid scale increases.The mean value of the watershed soil storage capacity increases from 30 to 57 mm when the grid scale aggregates from 50 to 1000 m.In addition, the peak of the density function also moves towards a high value.As a key variable in the Grid-Xinanjiang model, the bias indicates that for the coarse grid, more rainfall can be stored in the soil, and the runoff process will be more sensitive to the higher rainfall intensity.The hydrologic responses of the original Grid-Xinanjiang model were evaluated with the observed data from April to October in 1986 in the Yanduhe Basin. Figure 5 shows the processes of the average monthly precipitation, evapotranspiration, runoff and soil moisture at the different resolutions.For comparison, the soil moisture content θ is normalized to the relative soil moisture: θ * = (θ − θr)/(θs − θr).θs and θr are the saturation and the residual soil moisture, respectively, the values of which relate to the soil texture [37].The calibrated parameters at the 30 m resolution are directly adopted without recalibration at each grid scale.As can be seen in Figure 5, all these processes show an obvious scale-dependence.The spatially averaged precipitation reduces with the increasing grid scale.The bias of the precipitation mainly concentrates before March and after August, when the amount of rainfall is relatively low.The maximum margin of the monthly rainfall volume between 50 and 1000 m grids reaches the peak of 8% in September.For the relatively large rainfall volume, such as during June and July, the bias is only approximately 2%.This can be interpreted by the different spatial characteristics of rainfall types.As discussed by Ciach and Krajewski [38], rainfall fields are less variable for higher rainfall intensity.It is obvious that the smoothing effect will be more significant for the storm with higher variability.Compared to the precipitation, the runoff is more sensitive to the change of the grid scale.The bias of the estimation of precipitation, combined with the bias of totally watershed soil storage capacity, results in the inconsistencies of runoff at the different grid scales.For the coarse grid, the precipitation is underestimated.In addition, large soil storage capacity also allows more rainfall stored in the soil, rather than converting into runoff.The averaged bias of runoff between 50 and 1000 m resolution is approximately 28%.The bias of runoff reaches the peak around August.The runoff derived from 50 m is almost three times that from 1000 m.Since more rainfall infiltrates into the soil, the relative soil moisture increases with the grid scale.Therefore, as a function of moisture content of soil profile, the evapotranspiration is correspondingly higher for the coarse grid.The evapotranspiration and the relative soil moisture processes exhibit marked differences between March and August.This is because the evapotranspiration is more intensive during this period.The averaged bias of evapotranspiration and relative soil moisture caused by the spatial aggregation from 50 to 1000 m are 13% and 14%, respectively.Due to the decline of the measured evaporation, only slight differences occur in September and October.To evaluate the scale effect on the effective parameters, the model parameters require a recalibration at the different grid scales.To avoid the subjectivity of parameter estimation, automated calibration methods are needed.There has been a great deal of work on the research of automated calibration methods [39][40][41].In this study, the eleven model parameters are auto-calibrated based on the SCE-UA (Shuffled Complex Evolution method developed at the University of Arizona) algorithm, which has been found to be consistent, effective, and efficient in searching the globally optimum parameters [42,43].The initial parameters are referred to in Yao et al. [17].The Nash-Sutcliffe coefficient is used as an objective function.The iteration of SCE-UA algorithm will stop when the objective function cannot improve 0.1% over five iterations or the number of iteration exceeds 10,000.Table 2 lists the value of the effective parameters at four different grid resolutions.It is clear that the model parameters vary with grid scale in varying degrees.The ratio of the potential evapotranspiration to pen evaporation (K), the maximum watershed soil storage capacity (WMM) and the recession constant of the interflow storage (Ci) are relatively sensitive.The ratio of the potential evapotranspiration to pen evaporation decreases from 1.0 to 0.93 when the grid resolution increases from 50 to 1000 m.The maximum watershed soil storage capacity also decrease from 124 to 112.No significant changes were observed for the recession constant of groundwater storage (Cg) and two parameters in the Muskingum routing method (X,V).

Evaluation of the Subgrid Parameterization
As discussed above, the bias of the soil storage capacity at the different grid scales is an important reason for the scale-dependence of the model performances and the effective parameters.Therefore, a good representation of soil moisture capacity at the subgrid scale is necessary for a scale-invariant hydrologic response in the Grid-Xinanjiang model.Following the subgrid parameterization described in Section 3, Table 3 compares the annual water budget with and without the subgrid variability.The total annual evapotranspiration (Ep), runoff (R), surface runoff (RS), interflow (RI) and groundwater runoff (RG) at four different resolutions are calculated with the same parameters.As seen in Table 3, with the uniform grid, the annual evapotranspiration increases with the grid scale.The bias of the annual evapotranspiration between 50 and 1000 m resolution can reach 48 mm.The total runoff is significantly underestimated for the coarse grid.The annual total runoff at the 50 m resolution is 852 mm, and 763 mm at the 1000 m resolution.Note that the differences of interflow and groundwater runoff are quite small, and the bias of runoff mainly originates from surface runoff.This is because the runoff is partitioned into different components by the free water capacity, which is relatively small.The impact of the change of the free water capacity on the partitioning of interflow and groundwater runoff is limited.The runoff in excess of the free water capacity converts into surface runoff.Thus, the decrease in the total runoff from fine to coarse grid firstly leads to a change of the surface runoff.However, when the subgrid variability of the soil storage capacity is incorporated, the bias of the annual evapotranspiration between 50 and 1000 m significantly reduces to 3 mm.For the 1000 m resolution, the bias of evapotranspiration with and without the subgrid parameterization are considerably large, accounting for almost 22% of the annual evapotranspiration.For the non-uniform subgrid, the runoff increases by 23 mm and 101 mm for the 50 m and 1000 m resolutions, respectively.Compared to the uniform subgrid, both interflow and groundwater runoff increase marginally after incorporating the subgrid variability, and the differences at different resolutions are almost negligible.In Figure 6, the Nash-Sutcliffe coefficient is used as the criteria to evaluate the model performances of two subgrid schemes at the different grid scale.Because there is no measured data available for the actual evapotranspiration and soil moisture, the simulated results at the 30 m resolution are taken as the true values.In Figure 6, the dash-dot line is obtained from the original Grid-Xinanjiang model, and the solid line represents the model performance after incorporating the subgrid variability.With the uniform grid, the NSC of runoff process ranges between 0.93 and 0.98.The difference is more apparent for evapotranspiration and soil moisture processes.With the same parameters, the NSC of the evapotranspiration and the relative soil moisture decline to 0.88 and 0.82 at 1000 m from 0.98 and 0.97 at 50 m.However, when the subgrid variability is taken into account, the model performances at the different resolutions are quite consistent.The improvement of the model performance is significant, especially at the larger grid scales.The increases of the NSC for the three processes at 1000 m resolution reach 4%, 5% and 12%, respectively.It can be concluded that the subgrid parameterization can partly eliminate the scale-dependence of the model.Not only the total water budget, but also the performance at each time step is quite consistent at each resolution when the subgrid variability is taken into account, but it needs to be pointed out that the proposed subgrid parameterization only depicts the subgrid variability of soil storage capacity, and the variability of precipitation, evapotranspiration or flow routing velocity at the subgrid scale may still cause the inconsistency of model performances at the different grid scales.To further investigate the effect of subgrid parameterization, Figure 7 plots the simulated discharge at the outlet from the original Grid-Xinanjiang model versus the discharge after incorporating the subgrid variability at 50 m, 100 m, 500 m and 1000 m grid resolutions.It seems that the model tends to underestimate the discharge if the subgrid variability of soil storage capacity is neglected.The gap between the two schemes increases with grid scale.For the 50 m grid resolution, the discharge is closely fit to the 1:1 line.While for the 1000 m resolution, the discharge increases by an average of 6% after incorporating the subgrid variability.The bias of the discharge with and without subgrid variability can be as high as 40 m 3 /s.The large bias between the two schemes mainly occurs when the discharge is less than 100 m 3 /s.Even for small discharge, the gap is also evident.However, no significant change is observed for large discharge.It is because, for the saturation excess mechanism employed by the Grid-Xinanjiang model, large discharge is usually associated with relatively high soil content and intense rainfall.Once the watershed is fully saturated, all of the effective rainfall will directly convert into the runoff without loss.In this case, the effect of the spatial variability on runoff will no longer exist.It demonstrates that the subgrid variability is more necessary for the relatively dry period.This conclusion not only works for the Grid-Xinanjiang model, but also for the other saturation excess dominated hydrologic model.

Conclusions
Spatial variability is known to increase with scale.With the assumption of a uniform grid, the spatial averaging operation leads to the scale-dependence of model performance.In this study, we take the Grid-Xinanjiang model as an example to investigate the effect of subgrid variability on the hydrologic response.The processes of precipitation, runoff, evapotranspiration and soil moisture at four different grid scales (50 m, 100 m, 500 m and 1000 m) are compared.The calibrated parameters based on 30 m resolution are used for each simulation without recalibration.The results derived from the Yanduhe Basin show that the model performances vary with grid scales.As a key variable in the Grid-Xinanjiang model, the total watershed soil storage capacity increases from 30 to 57 mm when the grid resolution spatially aggregates from 50 to 1000 m.Depending on the different rainfall structure, the precipitation is underestimated for coarse grid.The bias of precipitation process and watershed soil storage capacity leads to the inconsistency of the hydrologic responses across scales.Calibrated by the SCE-UA method against the observed data, the model effective parameters also vary with grid scale in varying degrees.
Therefore, for a scale-invariant hydrologic response, the subgrid variability is required to be taken into account.In light of the sensitivity of the Grid-Xinanjiang model to the soil storage capacity, this study presents a subgrid parameterization method to incorporate the subgrid variability of the soil storage capacity.The topographic index is adopted as an auxiliary variable to correlate the soil storage capacity.On one hand, the topographic index and the soil storage capacity have similar spatial distributions and physical basis.On the other hand, the high resolution topographic index is more readily obtained than the soil property information.A beta distribution is introduced to represent the spatial distribution of the soil storage capacity within grid cell.Compared to the original Grid-Xinanjiang model, the simulated runoff, evapotranspiration and soil moisture processes are quite consistent at different grid scales when the subgrid variability is incorporated.The inclusion of the subgrid variability results in a significant increase of runoff and a decrease of evapotranspiration.For the annual water budget, the differences of runoff and evapotranspiration with and without incorporating the subgrid variability can reach 101 mm and 93 mm respectively at the 1000 m resolution.Due to the saturation excess mechanism employed by the Grid-Xinanjiang model, the correction of hydrologic responses resulting from subgrid parameterization is more obvious for relatively small rainfall event or dry conditions.The subgrid parameterization method benefits to reduce the necessity of recalibration when DEM resolution is changed, but it also needs to be noted that the proposed parameterization method only partly eliminates the scale-dependence of the runoff generation.Further research deserves more attention on the scale effect on flow routing.

Figure 1 .
Figure 1.Location and gauging stations of the study area.
where A is the soil storage capacity corresponding to the soil content and W 0 ; IWM is the cumulative density function of WM.

Figure 2 .
Figure 2. Schematic representation of runoff generation within a grid.

Figure 3 .
Figure 3.An illustration of the model performance of the Grid-Xinanjiang model with the subgrid parameterization in the Yanduhe Basin in 1986.Red dot, observed discharge; black line, simulated discharge; blue bar, precipitation.

Figure 4 .
Figure 4. Effect of the grid resolution on the density function of the soil storage capacity.

Figure 5 .
Figure 5.The monthly precipitation, evapotranspiration, runoff and relative soil moisture processes.

Figure 6 .
Figure 6.Model performances with and without subgrid parameterization.The solid line represents the NSC from the Grid-Xinanjiang model with uniform grid, and the dash-dot line represents the NSC after incorporating the subgrid variability.

Table 1 .
Performance of the Grid-Xinanjiang model for the calibration and validation periods.

Table 2 .
Model parameters and calibrated values at 4 different grid resolutions.

Table 3 .
Performance of the Grid-Xinanjiang model with and without the subgrid parameterization.