Physically-Based One-Dimensional Distributed Rainfall-Runoff Model Using the Finite Volume Method and Grid Network Flow Analysis

This work develops a grid based rainfall-runoff model (GRM), which is a physically based and spatially distributed model. Surface flow was analyzed using a kinematic wave model with the governing equations discretized using the finite volume method (FVM). This paper suggests a grid network flow analysis technique using variable rainfall intensity according to the flow directions to analyze one-dimensional flows between the grids. The model was evaluated by applying it to the Wuicheon watershed, a tributary of the Nakdonggang (Riv.), in Korea. The results showed that the grid-based, one-dimensional kinematic wave model adopted the FVM and the grid network flow analysis technique well. The simulation results showed good agreement with the observed hydrographs and the initial soil saturation ratio was most sensitive to the modeling results.


INTRODUCTION
Since Freeze and Harlan (1969) reported the potential of physically-based distributed models, these models have been developed for precise analysis of the effects of rapid land use and land development changes on watershed hydrological cycles (Abbott et al. 1986).As suggested previously (Beven 1989;Beven 1992;Stuart and Stocks 1993;Pullar and Springer 2000), in addition to the technical advances in computer performance and GIS (Geographic Information System), several studies have focused on the development and application of physically-based distributed models (Jain and Singh 2005;Todini 2011).
Many distributed models have been proposed to simulate watershed storm runoff (Singh and Woolhiser 2002).A range of model structures with a combination of various governing equations and numerical methods have been used to solve these equations (Refsgaard 1997;Du et al. 2007;Vansteenkiste et al. 2014).Some of them use the finite difference method (FDM) and the finite element method (FEM).Wang et al. (2002) developed a distributed watershed model to calculate the hillslope discharge using diffusion wave equations and the FDM in overland flow.James and Kim (1990) solved the diffusion wave equations in overland flow and the Saint-Venant equations in channel flow using the FDM.Vieux et al. (1990) developed a distributed model using kinematic wave equations with the FEM in space and the explicit scheme for time, and this model was developed into the Vflo TM model (Vieux 2004).Goodrich et al. (1991) divided the surface into two-dimensional TIN (Triangular Irregular Network) elements and analyzed the overland flow using the kinematic wave equations and the FEM.The FDM was used to calculate the channel flow.Recently, Du et al. (2007) suggested a simplified distributed model using the kinematic wave equations and the FDM for simulating storm runoff.Amaguchi et al. (2012) proposed a distributed model using vector-based catchment delineation to analyze the urban storm runoff with the FDM considering the effects of buildings, urban streets and sewer systems.
The finite volume method (FVM) is generally used to solve the diffusion wave equations or dynamic wave equations for watershed runoff modeling.Zhao et al. (1994) applied the FVM to solve shallow water flow equations.Jain et al. (2004) and Jain and Singh (2005) suggested using distributed models for calculating the runoff in grid based watersheds using the two-dimensional diffusion wave equations and the FVM.The RSM (Regional Simulation Model) (SFWMD 2005) uses the shallow water flow equations and the FVM for calculating runoff in a TINs based watershed.In contrast to these studies, the present study tried to solve the kinematic wave equations using the FVM for grid-based, watershed-runoff modeling.
In general, a two-dimensional model can provide detailed hydraulic information on the two-dimensional flow properties, but it requires more computational time than the one-dimensional model (Ahmad and Simonovic 1999;Leandro et al. 2011;Vacondio et al. 2014).In the case of watershed runoff simulations, previous studies (Du et al. 2007;Looper and Vieux 2012;Amaguchi et al. 2012) reported that the one-dimensional model can reproduce the observed flood hydrograph well with a shorter calculation time than the twodimensional model.Therefore, the kinematic wave model can be used to successfully simulate storm runoff and onedimensional analysis can reduce the computational time.
The aim of this study was to develop a grid based, onedimensional distributed rainfall-runoff model, GRM (Grid based Rainfall-runoff Model), using the kinematic wave equations and the FVM.A grid network flow analysis technique was suggested for a one-dimensional flow calculation in a grid-based distributed model.The model was evaluated with 4 rainfall-runoff events, and the sensitivities of the model parameters were analyzed.

Governing Equations
The GRM is a physically-based distributed model for simulating rainfall-runoff events.This is a grid-based distributed model that utilizes raster type distributed topographic data and time series hydrologic data as the input data.Figure 1 shows the runoff processes of the hydrological components in the model.Surface flow was divided into overland flow and channel flow.Infiltration excess flow and saturation excess flow from rainfall contribute to overland flow, infiltrated rainfall contributes to overland flow by return flow, and subsurface flow and base flow go into the channels.
Figure 2 presents a schematic diagram of the runoff in a grid where overland flow and channel flow occur simultaneously.The kinematic wave model was applied for runoff analysis.Infiltration processes were calculated using the Green-Ampt model (Green and Ampt 1911).In subsurface flow, the kinematic wave subsurface flow model proposed by Beven (1981) was applied.Continuity equations of overland flow and channel flow are as follows.Continuity equation for overland flow: Continuity equation for channel flow: where q (q = uh) is the flow per unit width of the control volume (L 2 T -1 ), u is the overland flow velocity in the x direction (L T -1 ), h is the flow depth (L), r is the rainfall intensity (L T -1 ), f is the infiltration rate (L T -1 ), q r is the return flow into the overland flow (L 2 T -1 ), A is the channel cross sectional area (L 2 ), Q is the discharge in the channel (L 3 T -1 ), q L is the lateral flow from overland flow (L 2 T -1 ), q ss is the subsurface flow (L 2 T -1 ), q b is the base flow (L 2 T -1 ), and ∆y is the width of the control volume (L).In Eq. ( 1) the effective rainfall and return flows that occur due to partial saturation on the surface are applied to the lateral inflow components (Dunne and Black 1970;Beven and Kirkby 1979).Equation (2) includes the rainfall occurring in the channels, and the overland flow, subsurface flow and base flow that go into channels are included as lateral inflow components.Infiltration can be calculated using the Green-Ampt equations (Chow et al. 1988).( ) where = -, h is the porosity, e i is the effective porosity, and K is the hydraulic conductivity (L T -1 ).
Subsurface flow (q ss ) is calculated by applying the kinematic wave model proposed by Beven (1981).In Eq. ( 5), saturated soil depth (D s ) is calculated by cumulative infiltration using Eqs.( 3) and ( 4), which are the sink term in the overland flow continuity equation [Eq.( 1)].The calculated subsurface flow is applied to the source term in the channel flow continuity equation [Eq.( 2)].
where D s is the saturated soil depth (L), and S a is the slope angle (Degree).S a uses the land surface slope of each grid cell calculated using the D8 method (O'Callaghan and Mark 1984).If the current control volume is saturated, the subsurface flows from the upstream control volumes contribute to return flow (q r ) in the current control volume using Eq. ( 6).
( ) q q y r s s i i where (q ss ) i is the subsurface flow from the upstream i th control volume neighboring the current control volume that is saturated.Percolation from soil layer A to soil layer B is calculated using Eq. ( 7) and base flow (q b ) in soil layer B is calculated using Eq. ( 8), which is based on Darcy's law (Freeze and Cherry 1979).
where p is the percolation during t D (L), K Bv is the vertical hydraulic conductivity in soil layer B (L T -1 ), K Bh is the horizontal hydraulic conductivity in soil layer B (L T -1 ), and D B is water depth in soil layer B (L).The hydraulic conductivity in the soil layer A [K in Eqs. ( 3) and ( 4)] is applied to set K Bv and K Bh , and K can be obtained from Chow et al. (1988) for each applied soil texture data set.D B is calculated using the percolation depth using Eq. ( 7).
The base flow interchanges with soil layer B and the channel are calculated using Eqs.( 9) and ( 10).If the water depth of soil layer B is higher than that of the channel, the base flow moves to the channel through the channel bed, and if the water depth of soil layer B is lower than that of the channel, the base flow moves to soil layer B through the difference in water depth.
Where h B is the water depth of the soil layer B, h ch is the water depth of the channel, and b is the channel bed width.

Numerical Solution
Figure 3 shows the runoff calculation process including the modeling components in the governing equations.The effective rainfall, subsurface flow and percolation are calculated in all cells.The return flow and overland flow are calculated for the cells with the 'overland flow cell' attribute, and the base flow and channel flow are calculated for the cells with the 'channel flow cell' attribute.In addition, for the cells with the 'channel and overland flow cell' attribute, all of the hydrological components are calculated.
The GRM discretizes the governing equations using the FVM and the Newton-Raphson method is used to calculate the convergent solutions for nonlinear terms in Eqs. ( 1) and (2) while calculating the flow area (h, A) and discharge (q, Q), and infiltration.Subscripts i, p, w, and e denote the number for the control volume, the center of the control volume, the control volume surface in the upstream direction (-x direction) from where inflows into the control volume occur, and the control volume surface in the downstream direction (+x direction) to which outflows occur, respectively (Patankar 1980).
GRM discretizes the governing equations by integrating Eqs. ( 1) and (2) in terms of the x-direction and time with the implicit scheme.The discretization equation for overland flow is expressed as Eq. ( 11).
where S i is the source term in overland flow ( ) S r f q y , and b i is the channel width for the control volume, CV i .The discretization equation for channel flow can be written as Eq. ( 12).
where S i is the source term in the channel flow ( ) S r y q q q i i i L i s si b D = + + + .

Grid Network Flow Analysis
When one-dimensional kinematic wave equations are analyzed using the FVM, the flows between the grids are divided according to the inflow directions from the upstream grids and outflow directions to the downstream grids into the orthogonal flows, diagonal flows and orthogonal and diagonal flow mixtures.In this study the lengths of the control volumes in them x direction ( x i D ) for the individual flows were established and the rainfall amount in a control volume was calculated based on the rainfall amount preservation conditions.

Orthogonal Flow
In orthogonal flows in Fig. 4a in all the distances from the control volume center to the control volume surface have a value equal to 1/2 of the grid size and ∆x i are calculated as follows.
x w e x i i i where x D is the grid size in the x direction (L).

Diagonal Flow
Figure 4b shows the diagonal flow scheme.In the di-agonal flow all distances from the control volume center to the control volume surface have a value equal to 1/2 of the length of the diagonal line for the grid ( x 2 D ).Therefore, x i D is calculated using Eq. ( 14).When x 2 D is used as the flow distance ( x i D ) in the diagonal direction, the control volume area ( x y i # D D ) becomes larger than that in the orthogonal flow case because the widths ( y D ) of all grids are not changed.Therefore, even if the rainfall intensities are equal to each other, the rainfall amount for each control volume can be changed by the flow direction.Equation ( 15) should be established because the rainfall intensity (r) used in the source term calculations is distributed evenly along the control volume length in the x direction ( x i D ) and the rainfall amount occurring in a grid should be constant regardless of the flow directions.Therefore, for grids where diagonal flow occurs, the rainfall [r * (L T -1 )] calculated from Eq. ( 16) is applied to the source term calculations.
x y r x y r where y D is grid size in the y direction (L).

Mixture of Orthogonal Flows and Diagonal Flows
If two or more upstream grids are adjacent to a grid, the inflow paths from the upstream grids may be mixtures of the orthogonal direction and the diagonal direction, as shown in Fig. 4c.In this case, the control volume length in the x direction ( x i D ) can vary with the number of grids where runoff occurs among adjacent upstream grids and their flow directions.The average values calculated using the flow directions and the number of upstream grids [Eq.( 17)] are used to apply x i D , which changes during the calculation process, as the same value in a control volume.In this case the rainfall (r * ) calculated by the rainfall conservation condition in one control volume is expressed as Eq. ( 18).
where N 0 is the number of adjacent upstream grids, ( ) w i n d is δw i in the n th direction among the adjacent upstream grids.

Study Area
The study area was the Wuicheon watershed, which is located upstream of the Museong stream gauge on the tributary of Nakdonggang (Riv.) in Korea.Most of the Wuicheon watershed consists of agricultural areas and mountains.The watershed area is approximately 472 km 2 .Figure 5 illustrates the Wuicheon watershed and hydrological observatories.

Topographic Input Data
Hydrological topographic data and thematic maps constructed in a raster layer were used as the main input data.If high resolution watershed data is used it could reduce the scale problem of spatially-distributed data but the application of high resolution data requires longer model runtime.In this study the spatial resolution of a grid layer was 200 × 200 m, considering the model runtime for the study area of 472 km 2 .
Table 1 lists the spatial data applied in this paper.Topographic data on the watershed was created using DEMs of 200 × 200 m, and used to set the target cells for runoff analysis, flow direction information and the slope on each cell.Data on soil and land cover was constructed using the soil and land cover maps provided by the National Academy of Agricultural Science and the Ministry of Environment.The maps were converted to 200 × 200 m grid layers.The soil map was applied to set the Green-Ampt model parameters for each soil texture and the soil depth parameters.The land cover map was used to set the impervious ratio and roughness coefficient for each land cover type.
Table 2 lists the areal ratio for each of the soil texture and infiltration parameters in the Green-Ampt model (Chow et al. 1988) in the Wuicheon watershed.Table 3 shows the area ratio and depth for each soil depth class (National Institute of Agricultural Science and Technology 1992).The soil texture map and soil depth map were extracted from the detailed soil map attributes in Korea.Table 4 presents   ratio, roughness coefficient, and impervious ratio for each land cover class in the Wuicheon watershed.In Table 4, the 'Impervious ratio' was calculated using the ratio of the effective rainfall to the total rainfall.In the water area, the model assumes that the rainfall is not infiltrated into soil directly and the total rainfall is equal to the effective rainfall, which contributes to direct runoff.Therefore, the model sets the impervious ratio of the water area to 1.

Rainfall Events
This study used the rainfall and flow data observed from 11 rain gauges (Miseong, Museong, Gomae, Hyoryeong, Daeyul, Sanseong, Hwasan, Hwasu, Seoksan, Seobu, and Uiheung) and one stream gauge (Museong).Four rainfall events, observed in 2007 and 2008 (Table 5), were selected, which have one hour time intervals.The GRM uses the distributed rainfall grid layers or mean areal precipitations in a text file as the input rainfall data.In this study, the distributed grid rainfall layers from Krigging method using the 11 rain gauges in the Wuicheon watershed were applied to the runoff simulation.
'Event 1' was applied to calibration and the others were used for validation.In the general calibration processes, two or more rainfall events are used for model calibration.On the other hand, in this study, the parameters calibrated by 'Event 1' could be applied properly to the other events, and no further calibrations using other rainfall events were needed.

Model Calibration
The model was calibrated using 'Event 1'.Table 6 and Fig. 6 present the estimated parameters and hydrographs before and after calibration.The results before calibration were calculated using the model in which the parameters were set as the default values shown in Table 6.The simulated results using the default parameters showed higher values than the observed hydrograph in the rising period and peak runoff rate.Therefore, the initial soil saturation ratio for 'Event 1' may be lower than 1, the default value of the model.Therefore, the initial soil saturation ratio parameter used for adjusting the rising period and peak runoff rate was estimated to be 0.93, which is smaller than the default value.At the time the runoff simulation was started the moisture content [θ in Eqs.(3) and ( 4)] was set using the initial soil saturation ratio parameter.The initial soil saturation ratio was set with a real positive value < 1, and applied to all cells in a watershed.Different initial soil moisture quantities were then assigned to each cell in a watershed because the cells have different soil depths and porosities to each other.To adjust the time to the peak runoff rate, the minimum slope of the channel bed was estimated to be 0.005, which is lower than the default.In Table 7, 'After cal.' shows the results of the runoff simulation using the calibrated model.
The simulation results were evaluated using the root mean square error (RMSE), Nash-Sutcliffe model efficiencies (ME), and correlation coefficient (CC). [ where Q obs and X are the observed flow (L 3 T -1 ), Q sim and Y are the simulated flow (L 3 T -1 ), Q obs is the average of observed flow (L 3 T -1 ), and N is the number of data.
In Fig. 6 the model calibration results showed good agreement with the observed hydrographs and the runoff reactions to rainfall were appropriate.In addition, the results from the calibrated model shown in Table 7 represent the observed hydrographs with 3.36, 1.6%, and ±0.0 hr relative errors in the total runoff volume, peak runoff rate and time to peak runoff rate, respectively.The RMSE, CC, and ME showed good statistical values.Therefore, in the validation for three rainfall events with the exception of 'Event 1', all of the parameters other than the initial soil saturation ratio, which are the same values as those shown in Table 6, were used.The initial soil saturation ratio was a temporally variable parameter that must be estimated for each event in the validations (Refsgaard and Storm 1996;Jain and Singh 2005;Sahoo et al. 2006, Tavakoli andde Smedt 2013).

Model Validation
Figure 7 and Table 8 show the simulation results for the 'Event 2', 'Event 3', and 'Event 4' storm events using the calibrated model.The model used in the validation had parameters calibrated using 'Event 1' ('Applied value' in   Table 6) except for the initial soil saturation ratio parameter.Because the initial soil saturation ratio is changeable according to storm events, it must be estimated properly for each storm event.Figure 7 shows the hydrographs and scatter plots of the simulated and observed data for the storm events.The simulated discharges represented the observed data very well.In Table 8, the total runoff volume, peak runoff rate and time to peak runoff rate showed relative errors of 0.3 -2.3%, 1.1 -1.7%, and ±0.0 to +1.0 hr, respectively, compared to the corresponding observed values.In addition, the RMSE, CC, and ME showed good agreement with the observed data.Therefore, the calibrated model could simulate the hydrologic events selected for the validation properly.These results show that the calibrated model has temporal transferability and can simulate storm events in different time periods.

SENSITIVITY ANALYSIS
The sensitivities of the parameters estimated by the user including the initial soil saturation ratio, channel roughness coefficients and minimum slopes of channels to the total runoff volume, peak runoff rate and time to peak runoff rate were analyzed and the individual parameter effects on the model calibration were assessed.The channel roughness coefficient of 0.045 and channel minimum slope of 0.005 estimated for model calibration and the value of initial soil satu-ration ratio of each event presented in Tables 7 and 8 were established as reference values for sensitivity analysis.The reference values for the initial soil saturation ratio, channel roughness coefficient and minimum slope of channel were changed to ±10% at ±2% intervals and the values were applied to runoff simulations, and the changes in the peak runoff, peak time and total runoff volume were reviewed.
Figure 8 shows the rates of changes in peak runoff, peak time and total runoff volume, resulting from changes in the initial soil saturation ratio, channel roughness coefficients and minimum channel slopes.The initial soil saturation ratio is an important parameter that affects the level of possible infiltration into the soil.The peak runoff, peak time and total runoff volume reacted sensitively, suggesting that the estimated initial soil saturation ratio has significant effects on the model calibration.The channel roughness coefficients were less sensitive to the peak runoff and total  runoff volume compared to the initial soil saturation ratio while showing similar sensitivity to peak time to that shown by the minimum channel slope.Therefore, the estimation of channel roughness coefficients can be applied to fit the peak runoff and peak time and also to the fit the hydrograph forms in detail.The minimum slopes of the channels showed similar sensitivity to the peak runoff rate and peak time to that shown by the channel roughness coefficients.On the other hand, the sensitivity of the minimum slopes of the channels to the total runoff volume was very low.Therefore, an estimation of the minimum slopes of the channels can be applied to fit the peak runoff rate and peak time rather than the total runoff volume and.

DISCUSSIONS AND CONCLUSION
The GRM, which is a physically based distributed rainfall-runoff model, was developed in this study using the FVM.The model simulates the rainfall, infiltration, overland flow, channel flow, subsurface flow, and base flow hydrological components.The kinematic wave model was used to simulate the overland flow, channel flow and subsurface flow, and the Green-Ampt model was used to calculate infiltration.Subsurface flow was calculated using a kinematic wave model and base flow was analyzed based on Darcy's law.The FVM was applied to discretizing the governing equations and the Newton-Raphson method was applied to calculate the nonlinear terms.
The grid network flow analysis technique was suggested to calculate the one-dimensional flows in the twodimensional grid domain.The width of the control volume was not changed, but the length of the control volume in the x direction was changed in accordance with the flow relations between the grids.To conserve the amount of rainfall occurring in each control volume, the rainfall to be applied to the source terms of the discretization equations was recalculated in accordance with the changing lengths in the x direction based on the mass conservation of discretization equations concept in the FVM.
When applying the model to the Wuicheon watershed, the physical parameters established based on the DEM, soil and land cover data were applied without corrections, and the observed hydrographs were well reproduced by correcting only the initial soil saturation ratio and minimum channel slope among the parameters estimated by the users.This means that the data and physical parameters applied to the runoff simulations reflect the characteristics of the watershed appropriately and that the GRM analyzes appropriately the rainfall-runoff phenomena in the watershed physically.
Sensitivity analysis revealed the initial soil saturation ratio to have relatively high sensitivity to the total runoff volume, peak runoff rate and peak time.The channel roughness coefficient showed a lower sensitivity to the peak runoff rate and total runoff volume compared to the initial soil saturation ratio and showed similar sensitivity to the peak time as the minimum channel slope.The minimum channel slope showed sensitivity to the peak runoff rate and peak time similar to that shown by the channel roughness coefficient.Therefore, the initial soil saturation ratio had a very large effect on the model calibration.Moreover, the estimated channel roughness coefficient had slightly greater effects on the fit of the total runoff volume and the peak runoff rate than the peak time and the minimum channel slope showed low sensitivity.

Table 2 .
Areal ratio for each soil texture and parameters for the Green-Ampt model in the Wuicheon watershed.

Table 3 .
Areal ratio and depth for each soil depth class in the Wuicheon watershed.

Table 4 .
Areal ratio, roughness coefficient, and impervious ratio for each land cover in the Wuicheon watershed.

Table 6 .
Estimated values for the GRM parameters for model calibration.

Table 7 .
Model results for the calibration.

Table 8 .
Model results for the validation.