Simulation of Surface Runoff and Channel Flows Using a 2D Numerical Model

1. Introduction

Numerical simulation is increasingly used for studying overland flows. Since runoff drives soil erosion and landscape evolution, the runoff models provide a foundation for modeling soil erosion, rill erosion, and related processes at the watershed scale [1, 2]. Models involving different levels of abstraction have been proposed [3, 4, 5]. Two commonly used models are the diffusion wave (DW) and kinematic wave (KW) models [6, 7, 8, 9]. The KW models set the friction slope to be equal to the bed slope and ignore the inertial terms [10]. The method has been successfully used to describe overland flows [11, 12, 13, 14]. The governing equations are highly nonlinear and do not have general analytical solutions, so one has to solve them numerically for practical cases [15]. The models based on full Saint-Venant (SV) equations have also been applied and produced better results.

Two-dimension models are generally used for cases with irregular domains. A distributed rainfall-runoff model using the KW approximation solved by an implicit finite difference scheme was developed [16], but channel flows are computed using a separate KW model. Fully two-dimensional shallow water equations are being utilized for modeling overland flows in late 1980s [17]. A two-dimensional finite difference (FD) runoff model was developed by solving 2D SV equations [18]. Shallow water equation-based 2D models [19] were used for runoff over an irregular topography of experimental scale with infiltration processes considered and in rural semiarid watersheds for overland flows generated by storms [20].

In addition to finite difference method (FDM), the two-dimensional finite element (FEM) and finite volume methods (FVM) have been used for overland flow simulations. A FEM KW model was developed by Liu et al. [21] for simulating runoff generation and concentration over an irregular bed and reproduced experimental results. Tests [15] indicated that the FVM-based 2D SV model performed better than that of FDM. Costabile et al. [22] solved the shallow water equations using the FVM and applied the resulting model to simulate a real event on a watershed of 40 km². Nunoz-Carpena et al. [23] solved the KW equation using the Petrov-Galerkin method. Venkata et al. [24] developed a Galerkin DW FEM and applied it to a small watershed. Singh et al. [25] simulated runoff processes by solving the 2D shallow water equations with a shock-capturing scheme and the FVM. Shirmeen et al. [26] showed results of a validated, FEM 2D model in predicting runoff from a flat agricultural watershed.

In order to check numerical models’ mathematical correctness and physical applicability, the developed computational models have been tested with analytical solutions, experimental, and field data. Iwagaki [27] studied runoff using analytical methods and experimental data; several specific solutions were developed based on the characteristic method. Govindaraju et al. [28] developed analytical solutions using KW and DW approximations. Comparisons of analytical solutions, numerical solutions, and experimental data were discussed. Singh [29] detailed the KW model’s analytical and numerical solutions and their wide applications. Cea [30] tested FVM using an experimental watershed with a complex shape. These overland flow models use simplified equations and need to specify pre-existing channel networks, which make it difficult to simulate soil erosion cases with hill-slope evaluation and mixed sheet-channel flow conditions.

CCHE2D is a physically based model, which treats the entire watershed including the channels and ditches as one continuous domain. One does not need to differentiate overland sheet flow and channel flow calculation areas using grid cells and 1D channel networks as is done in GSSHA [31], WASH123D [32], NIKE-SHE [33], and SHETRAN [34]. It is also not necessary to employ arbitrarily shaped sub-watersheds and 1D channel networks as is done in the CCHE1D model [35]. In these models, 2D DW equations or KW equations are solved for the overland flow using finite difference methods, and the 1D SV equation is solved in the prescribed channel networks. In contrast to these models, in CCHE2D, hydrodynamics over the entire watershed is simulated using only 2D equations discretized on an irregular quadrilateral finite element mesh, which is generated using digital elevation model (DEM) data. The simulated overland sheet flow and channel flow are seamlessly connected everywhere in the domain and the channel network is formed automatically. This method may be more applicable when sediment transport, rill erosion, or gully erosion processes in watersheds are considered.

In this study, the CCHE2D model is modified and applied to simulate watershed hydrological processes. CCHE2D is a general hydrodynamic model for unsteady, turbulent free flows, sediment transport, and pollutant transport. It has been validated and applied widely to simulations of channel flow, flooding, coastal flow, bed topographic change, and chemical contamination in aquatic environments [36, 37, 38, 39, 40].

The major objectives of the present paper are to assess the accuracy and the effectiveness of this FEM in predicting overland runoff processes, and its applicability to practical agricultural watersheds with ditches and natural stream channels. The approach of the study followed the recommendations of [41] for quality assurance that numerical models have to be verified and validated using analytical solutions, physical experimental data, and field data. The validated numerical model was used to simulate and characterize the hydrological processes of an agricultural watershed in the Mississippi River alluvial plain where farm fields are drained and separated by ditches and stream channels. A limitation was found in the interpolation method when it is applied to the water surface elevation of the sheet runoff. A numerical scheme was developed and implemented for improving the bilinear interpolation. The present study focused on watershed surface flow processes over bare soils; interception, evapotranspiration, and infiltration were not considered.

2. Mathematical model

Surface runoff due to precipitation is typically quite shallow and can be aptly represented by the 2D shallow water equations within the CCHE2D model [36, 38]. The water surface elevation of the runoff flow, η, is calculated by the continuity equation in a Cartesian coordinate system

in which h is the local water depth, t is time; R is rainfall intensity, which may vary in time and space, and u and v are depth-averaged velocity components in x and y directions, respectively. The depth-integrated 2D momentum equations for turbulent flows are as follows:

in which g is the gravitational acceleration, ρ is water density, τxx, τxy, τyx, and τyy are depth-averaged Reynolds stresses, and τbx, τby are bed shear stresses. In the overland runoff area, the Reynolds stress terms vanish, and Eqs. (2) and (3) become the shallow water equations. The Reynolds stress terms remain significant in the part of the domain with channel and concentrated flows. A special finite element method called the efficient element method is adopted in the model, in which a collocation approach is used to discretize the equations in a structured quadrilateral nonorthogonal mesh system. A partially staggered grid is used for solving these equations. A velocity correction method is used to couple the continuity equation and the momentum equations. More details about this model’s numerical methodology and techniques can be found in earlier publications [36, 38, 42].

The full Eqs. (1)–(3) are applicable for general flow conditions. In realistic cases where runoff and channel flow conditions coexist, a general flow model is necessary. Under the sheet flow condition, the advection and turbulence stress terms in the momentum equations vanish because the dominant forcing for the overland flow is the gravity and bed shear stress. The water depth is very small, and water surface slope and bed slope become almost the same:

in which b is the bed elevation. The general flow equations then become the KW equations. Under this condition, the flow is completely dominated by the bed slope. Shear stresses on the bed are evaluated in conjunction with the Manning equation as:

in which n is the Manning roughness coefficient and U=u2+v2is the total velocity magnitude.

3. Interpolation of water surface elevation

CCHE2D uses a partially staggered method: the velocities are solved at collocation points and the pressure (water surface) is solved at cell centers [36]. A bilinear interpolation method is used to interpolate the water surface elevation solution to the collocation nodes where the momentum equations are solved. The bilinear interpolation works well for general channel flow simulations because the water depth is large in comparison with the variation of bed surface and the mesh size. When overland sheet runoff is simulated; however, the water depth is very small; it is often less than the microelevation variation of bed topography represented in an element. In this case, the interpolated water surface elevation may be lower than the bed if the bed is concave down and vice versa. This is a limitation of the interpolation method. In the concave down case, dry nodes are created artificially; in the concave up case, artificial masses of water could be erroneously created. Figure 1 illustrates this problem in one dimension. The problem occurs whenever irregular bed topographies are encountered. A correction is therefore necessary to the interpolation over the surface runoff area.

A numerical scheme has been developed and implemented in CCHE2D to correct the interpolation error [43]. Figure 2 illustrates how the scheme is formulated in one dimension with an exaggerated vertical scale. Eq. (7) is the formulation to compute the correction value Δb for nonuniform meshes, and it is simplified to Eq. (8) if the mesh is uniform. It is straightforward to extend Eqs. (7) and (8) to two dimension. Water depth at the cell centers is positive, without this correction, the depth at the middle point would become negative because the interpolated water surface elevation is below the bed. This scheme is necessary and effective when cases with irregular topography are simulated

where b1, b2, and b3 are bed elevation, Δb is the interpolation correction, b2iis the linearly interpolated value, and Δx1 and Δx2 are mesh spacing (Figure 2). The interpolated water surface elevation needs to be corrected by Δb.

4. Analytical verification

Two analytical solutions were obtained by solving a one-dimensional kinematic equation analytically for rain-generated runoff by [44, 45]. The solution of sustained rains for the runoff to reach a steady state [44] and the solution for rainfall that stops before the runoff becomes steady [45], including the tailing stage solution after rainfall stops, were provided. The governing one-dimensional kinematic equation for deriving these solutions is:

in which q is the discharge of water per unit width (m²/s), k is an exponent (=5/3), and α (=5) is a coefficient (m^2−k/s). These analytical solutions were realized for a few simple cases: runoff due to steady rainfall intensity on a uniform planar area of 200 × 1 m with a slope of 1.0%. The rainfall intensity was R = 2.7 × 10⁻⁵ m/s, and the Manning’s coefficient was n = 0.02 m^−1/3s. For comparison, the same case was simulated using CCHE2D and a 10 × 100 point 2D mesh with uniform spacing. The solutions were recorded at cross sections located at 50, 100, 150, and 200 m, from the upstream end of the plane.

Figure 3 shows the comparisons of the simulated runoff and the analytical solutions for the sustained rain collected at the four cross sections. Hydrographs at each cross section indicate that equilibrium runoff (steady state) is reached before the rain stops at T = 1000 s. The runoff is always nonuniform, and the peak discharge increases in the downstream direction. At first, the flow is unsteady (rising limb), then becomes steady until T = 1000 s, and finally becomes unsteady in the falling limb. The runoff reaches equilibrium earlier at locations closer to upstream. The simulation is a little less than the analytical solution at the time approaching the peak discharge, particularly near the downstream. The solution can be improved by reducing the local mesh size effectively.

Figure 4 shows a case in which the rainfall stops before runoff reaches steady state (T = 200 s); the hydrographs, thus, have a different pattern. The peak discharge is reached at the time the rain stopped and is the same for all cross sections. The peak discharge for the lower cross sections lasts longer because the flows at the lower locations are sustained by upstream contributions. The runoff recession is earlier for upstream locations. The shape of the two sets of simulated hydrographs at all cross-section locations corresponded well with the analytical solutions.

5. Experimental validation

CCHE2D model was validated using experimental data sets collected from the literature. All of these cases were carried out on impervious overland flow planes. The only quantity measured in these experiments was the downstream runoff discharge.

5.1 Case 1

Morgali and Linsley [46] obtained two sets of experimental runoff data. Their tests were carried out over a straight turf surface of 21.95 m long with a constant slope (0.04) and width. The Manning’s coefficient, n, was found to be 0.5 m^−1/3s. The rains had two different intensities and were uniform along the slope for 1200 s (20 min). Figure 5 compares the experimental data and the numerical simulations. The analytical solution for this test condition [44] is also presented in Figure 5. It was found that these runoff experiments fit well with the analytical solution. A 110 × 10 points uniform mesh and 0.01 s time step were used for the numerical simulation. The CCHE2D numerical results showed good agreements with the analytical solution as well as the experimental results (Figure 5). The rising limb of the discharge hydrograph and the peak discharge were captured very well by the simulations. The processes of the two experiments, 1A (R = 92.96 mm/h) and 1B (R = 48.01 mm/h), look similar because the only difference in the experiments was rainfall intensity. The peak discharges of the experiments occurred at approximately 850 and 1100 s, respectively, for Case 1A and Case 1B. The numerical solutions of CCHE2D agreed well with the experimental data. The peak discharges for Case 1A and Case 1B are 5.67×10⁻⁴ and 2.93 × 10⁻⁴ m³/s, respectively. The times to peak discharge for Case 1A resulted from the analytical solution, CCHE2D and the experiment, are 760, 850 and 950 s, respectively. The differences among the three are less for the Case 1B (Figure 5).

Figure 5 also compares the simulation results of CCHE2D and the model results by Govindaraju et al. [28]; the two numerical solutions agree well for the case with the higher rain (1A), but the fit of their solutions based on the SV equations does not correspond well for the case with the smaller rainfall (1B). The results of CCHE2D also outperform the analytical solution of the DW approximation [28].

5.2 Case 2

Cea et al. [30] conducted three runoff experiments of complex topography and simulated these cases using a 2D unstructured FVM. The experimental watershed was a rectangle (2 × 2.5 m) made by three planes of stainless steel, each of them with a slope of 0.05 (Figure 6). Two dikes (1.86 and 1.01 m in length) were placed in the watershed to vary the topography. Rainfall intensity, duration, and runoff hydrographs were measured. As a result, the runoff direction, distribution, and pattern of the hydrograph were affected. The runoff was accumulated and became channel flows along intercepting lines of slopes and dikes. Since both overland flow and channel flow are involved, faithful simulation requires solving full governing Eqs. (1)–(3). The rainfall applied to each test case was different. In the first test (2A), the rainfall intensity was 317 mm/h for 45 s. In the second test (2B), the rainfall intensity was 320 mm/h for 25 s; then it was stopped for 4 s and restarted for an additional 25 s with the same intensity. In the third test (2C), rainfall intensity was 328 mm/h. The rainfall was applied for 25 s; then it was stopped for 7 s and then restarted for another 25 s.

In this study, CCHE2D was applied and the numerical results were compared with experimental data. The watershed was modeled using an irregular structured mesh with the cell size ranging from 0.034 to 0.009 m; the mesh was refined near the main channel and the outlet for improving results. The Manning’s roughness coefficient was set equal to 0.009 m^−1/3s. The simulation time was 120 s for each case. The channel flow and runoff sheet flow coexisted: the runoff from the watershed surface was accumulated in the bottom of the watershed channel with a triangle-shaped cross section formed by the side slopes. Results of cases 2A and 2C are shown in Figures 7 and 8, respectively.

Figure 7 shows the comparison between the numerical solution and experimentally observed runoff hydrograph of Case 2A. The solution of the CCHE2D model agrees very well with the experimental results. The flow discharge increased continuously once the rain started. The peak discharge occurred at the time the rainfall stopped (at 45 s). Although the rising and the falling limbs of the hydrograph were slightly overestimated, the shape of the hydrograph and the peak discharge were aptly predicted.

Figure 8 shows the comparison between the numerical and experimental runoff hydrographs of Case 2C. The shape of the hydrograph was successfully predicted. The interval between the two rainfall peaks was 7 s. The first runoff peak discharge occurred at the time the rainfall stopped, at 25 s. The runoff discharge decreased for approximately 10 s and then increased. The second runoff peak discharge occurred at approximately 57 s. The simulated processes and the observed physical processes showed a good general agreement; it also matched well with the model results of [30].

Figure 9 shows the simulation results at t = 54 s (the peak of the second rainfall) for Case 2C: (a) simulated water depth contour distribution, (b) simulated flow unit discharge pattern and (c) velocity vector distribution in the watershed. The distributions indicate how the overland sheet flow, under the influence of dikes and topography, concentrates into channels and flows out of the watershed. The flows over the slopes are sheet runoff, but complex recirculations are developed in the main channel. The water surface is no longer parallel to the bed surface. These flows cannot be represented by KW, DW, and SV models.

6. Application to a real watershed mathematical model

This section presents the application of CCHE2D to a sub-watershed of the Howden Lake watershed, an 18 km² agricultural watershed in the Mississippi River alluvial plain (Figure 10). In this region of low relief, watersheds are configured by farm fields drained by culverts, ditches, and intermittently flowing streams called bayous. During periods between runoff events, the channels contain standing water. The studied sub-watershed was upstream of a gaging station on an intermittently flowing bayou. The average annual precipitation in this region is about 1440 mm. Precipitation occurs as intense thunderstorms or low-intensity rains associated with major frontal movements. The latter type of events may stretch over several days of drizzle and sporadic showers. During growing seasons, channels experience some flow and stage fluctuation due to irrigation withdrawals and return flows.

Watershed topography was surveyed by airborne LiDAR with a 1.5 m horizontal resolution. The vertical accuracy was 15.0 cm RMSE or better. The watershed elevation ranges from approximately 43.89–48.99 m. A nearly uniform fine mesh (mesh spacing = 3.76–4.98 m) was generated for the simulation with the ditches and small streams between the plots further refined locally. Cultivated fields are connected to the streams and ditches with drainage culverts, which often convey water from one sub-watershed to another. The locations of culverts in the study watershed were identified in a field survey and incorporated in the numerical mesh.

Soil data were obtained from the Soil Survey Geographic (SSURGO) database [47]. The watershed is covered mostly by soils with high clay content, which is typical of the region [48]. Infiltration is, therefore, negligible and was not considered in the simulation. Precipitation and flow stage data were measured by field instrumentation. Because the stream instrumented with the gage station has complex conditions, it was difficult to collect reliable velocity data during a rain event. Only stage data were available. As a result, the gage station does not have a discharge-stage rating curve. Development of a rating curve using simulations and measured data for this site would be helpful for understanding the hydrologic processes in these watersheds.

Because the Howden Lake watershed is of low-relief, it was often difficult to determine the boundaries between sub-watersheds in field surveys or on topographic contour maps. For example, the runoff from a piece of field may flow in two directions into two sub-watersheds, and the location of the divide line might be identified only from the runoff flow distribution during a simulation. Normally, the outline of a watershed is a given condition for a hydrology study. In this study, the exact boundary outline was not firmly established even after field surveys. A larger area containing the studied watershed was simulated, and the watershed boundary and area were finally defined by the simulated runoff and channel flow patterns. The boundary outline of the studied watershed (Figure 11) contributing to the gage was identified by visually checking the simulated overland flow directions of CCHE2D.

In the simulations, the streams and ditches between farming plots were represented using DEM elevations like flat surface areas. No channel networks were prescribed, but the simulated surface runoff flowed logically to the ditches and to the stream channel. No other watershed analysis tools were needed. Although the study results presented later are for this identified watershed, the spatial domain of numerical simulations was several times larger (Figure 10). The northern side of the stream channel had been blocked by farmers, so the overland flow from the watershed entered the stream in the middle and flowed in a southwesterly direction (Figure 11). The water from this identified watershed pasted the gage, while runoff from the region outside this watershed was discharged from the simulation domain via other ditches and streams. The area of this watershed, including cultivated land, drainage ditches and a stream segment, was found to be 973,700 m². In this area, the topographic elevation ranges from approximately 46.77–47.49 m in one plot and from 47.27 to 48.09 m in another. The mean slope of the fields is 0.0097 and 0.0098, respectively.

Several observed storm events were selected for the model application. To reduce minor losses of water due to evaporation, soil wetting and infiltration, etc., only large rain events were considered. The rainfall event in April 2011 (Table 1) was first used for simulation. Figure 12a shows the detailed ground elevation contour of a small simulation area (dashed rectangle area in Figure 11). The elevation of this area ranges from about 46.8 to about 47.4 m. Figure 12b shows the direction vectors of the runoff near the end of the simulation. Because the water is very shallow, the flow direction is highly affected by the ground topography. Figure 12c and d shows the direction vectors and water depth distribution at the peak time of the rainfall.

Although the variation of the bed surface topography is very small, the simulation shows how the runoff is controlled by microtopography (Figure 12a). At the peak time of the rainfall, the overall water depth in this area (Figure 12d) is much deeper, and the flow directions (Figure 12c) are less affected by the local microtopographic features. The flow on the right side of the domain is still sheet runoff under the KW condition; while on the left side, the water depth is more than 0.2 m, and the flow is no longer governed by the KW condition. This model provides the outflow hydrograph as well as the temporal and spatial distribution of the water depth and flow velocity, which can be used for studying soil erosion, agro-pollutant transport, and water quality.

The gage station (Figure 11) recorded the channel water surface elevation at regular time intervals, but velocities were generally too low for accurate measurement, and therefore, water discharge was not measured. In order to better understand the watershed hydrology, a rating curve of the form:

was developed using simulated discharge, in which L is the measured water surface elevation, r and z are parameters, and L0 is the initial water surface elevation prior to a rainfall event. Eq. (11) has two unknown parameters, but there is only one relationship available for determining their values. The total volume of runoff, obtained by numerically integrating Eq. (11) in time, is equal to the rain volume, VR:

in which Li is the measured water surface elevation at the gage station. With Eq. (12) satisfied, values of r and z that best fit the shape of the discharge hydrograph computed using Eq. (11), and that of the numerical simulation, were determined for each event by trial and error.

Attempts were made to fit all simulated curves using a single set of values for r, z and a mean base stage L0, but the result showed unacceptable discrepancies. L0 varied due to antecedent precipitation, downstream discharge control, sedimentation, and water usage between events. The range of L0 for the studied events is 0.33 m (Table 1). Given the complexities of the hydraulic regime in the water body, varying from standing to moving state and with varying downstream controls, variable rating curve parameters are sensible. Event-specific rating curve parameters are not ideal but are useful in a research context.

Manning’s roughness coefficient (n) is a major factor in the determination of watershed runoff characteristics and generally reflects ground cover and management. The event on April 27–28, 2011 was used for initial calibration of Manning’s coefficient. The studied watershed is cultivated with soybeans (Glycine max L. Merr.), corn (Zea mays L.), cotton (Gossypium hirsutum L.), and rice (Oryza sativa L.). The sensitivity of the CCHE2D model in Howden Lake watershed to Manning’s n was examined using a wide range of values from 0.030 to 0.30 m^−1/3s. Smaller Manning’s n results in a higher runoff peak discharge and an earlier peak flow arrival time. A visual comparison of discharge hydrographs based on stage measurements and numerical simulation (Figure 13) indicates that n = 0.3 m^−1/3s is the most appropriate choice for the overland runoff area because the peak times of these runoff events are consistent. Considering that the depths of the sheet runoff are much smaller than the microtopographic irregularities over the fields, the calibrated n represents not only the bed resistance but also form drags due to the microbed forms, crop residue, and vegetation. This n value agrees with the recent runoff studies [25, 31, 49] in cases of overland flows, including those in the Goodwin Creek Experimental Watershed in Northern Mississippi. There are numerous trees, bushes, and weeds growing along and within the channel, thus, n = 0.16 m^−1/3s was used for the channel and kept unchanged for other rain event cases.

The total observed rainfall volume for the April 27–28, 2011 event (Figure 13) was approximately 86,000 m³ (88.32 mm). The total simulated runoff volume is about 80,600 m³ (83.78 mm), which is reasonable because the hydrograph recession limb extended past the simulation termination at 47 h. There were several small rain events that occurred before the event shown in Figure 13, so the runoff volume based on the observed water surface elevation may include recession of the earlier events.

Figure 14 compares the discharge hydrographs of several additional runoff events computed using Eq. (12) and that of the numerical simulations. The identified parameters for these events, r and z, are listed in Table 1. Events 9/2011 and 10/2013 have one major peak, while those of 11/2011 and 5/2013 each have two major peaks. The simulated hydrographs fit well with those computed using Eq. (11). The two rain peaks of the 5/2013 event were separated by about 2 h, but those of the 11/2011 event were separated by 15 h. The runoff of the 5/2013 event showed only one peak because the two rain peaks were very close, and the runoff peaks were superimposed. However, the temporal separation of the two peaks of the 11/2011 event was much longer. Therefore, the superimposed hydrologic response also displayed two peaks. These watershed responses were reproduced by the numerical simulations.

As noted above, the watershed has multiple field ditches that convey runoff into the channel (Figures 10 and 11). Ditch and channel flow were simulated together with the overland sheet runoff. Figure 15 shows the simulated flows in the channel network of the watershed. The contours represent the distribution of the unit flow discharge. The vectors in the ditches and in the stream formed a channel network indicated by the large velocity vectors; those in the runoff area are too small to be seen. The flows in the stream are turbulent when the rainfalls are large. Because no velocity data were acquired, the simulated velocity results in the channel were not validated.

7. Conclusions

The numerical model CCHE2D was used to model sheet runoff from watersheds, large and complex enough to include both overland and channel flow processes. The model was systematically verified and validated using analytical solutions and experimental data due to steady and unsteady rainfall intensity, and applied to a real world watershed. Good agreement between the analytical solutions, experimental data, and numerical simulations were obtained. For the experimental cases involving complex watershed shapes, the numerical model has the ability to simulate runoff over the slope surfaces and the channel flows.

A numerical scheme was developed to correct the bilinear interpolation of the water surface elevation from its solutions at the staggered cell centers to the collocation nodes. The scheme was necessary and effective for obtaining good sheet runoff simulation results in watersheds with irregular topography. One would have to smooth the ground topography if a model requires the interpolation of water surface solution under this condition.

The model was applied to an agricultural watershed in the Mississippi River alluvial plain. It was useful to identify the boundary of the monitored watershed and develop the rating curve at the gage station of the watershed. Several significant runoff events were selected for simulation. Each of the simulated runoff hydrographs and the rating curves agreed well with those observed in the field. The sensitivity of the model to overland sheet flow friction was studied. An increase in the bed surface friction coefficient significantly diminishes the peak of runoff discharge, delaying its time of arrival. Values of n = 0.2–0.3 m^−1/3s for overland flow were found to be adequate to best fit the numerical simulations and the observed data in the studied watershed. With a high-resolution mesh, the model can predict the complex surface runoff pattern over the agricultural land. Ditch and stream channels in the domain are a connected channel network. The model is able to simulate sheet runoff, turbulent channel flow, and their transitions seamlessly. The simulated hydrological processes for several storm events fit well to those observed at the gage station. The capability would be useful for studies related to soil erosion and agro-pollutant transport. The model is currently used for watershed applications without considering interception, evapotranspiration, and infiltration. Additional work is needed to further extend the research in these areas.

Acknowledgments

This work is supported in part by USDA Agriculture Research Service under the Research Project No. 6060-13000-025-00D (NCCHE) monitored by the USDA-ARS National Sedimentation Laboratory (NSL). Support is also in part by the Southeast Region Research Initiative (SERRI) project and the University of Mississippi (UM).