A new modeling approach for simulating microtopography-dominated, discontinuous overland flow on infiltrating surfaces

Advances in Water Resources 78 (2015) 80–93

Contents lists available at ScienceDirect

Advances in Water Resources journal homepage: www.elsevier.com/locate/advwatres

A new modeling approach for simulating microtopography-dominated, discontinuous overland flow on infiltrating surfaces Jun Yang, Xuefeng Chu ⇑ Department of Civil and Environmental Engineering, North Dakota State University, PO Box 6050, Fargo, ND 58108, United States

a r t i c l e

i n f o

Article history: Received 13 May 2014 Received in revised form 30 January 2015 Accepted 5 February 2015 Available online 14 February 2015 Keywords: Overland flow Microtopography Puddle dynamics Depression Infiltration Digital elevation model

a b s t r a c t Realistic modeling of discontinuous overland flow on irregular topographic surfaces has been proven to be a challenge. This study is aimed to develop a new modeling framework to simulate the discontinuous puddle-to-puddle (P2P) overland flow dynamics for infiltrating surfaces with various microtopographic characteristics. In the P2P model, puddles were integrated in a well-delineated, cascaded drainage system to facilitate explicit simulation of their dynamic behaviors and interactions. Overland flow and infiltration were respectively simulated by using the diffusion wave model and a modified Green–Ampt model for the DEM-derived flow drainage network that consisted of a series of puddle-based units (PBUs). The P2P model was tested by using a series of data from laboratory overland flow experiments for various microtopography, soil, and rainfall conditions. The modeling results indicated that the hierarchical relationships and microtopographic properties of puddles significantly affected their connectivity, filling–spilling dynamics, and the associated threshold flow. Surface microtopography and rainfall characteristics also exhibited strong influences on the spatio-temporal distributions of infiltration rates, runoff fluxes, and unsaturated flow. The model tests demonstrated its applicability in simulating microtopography-dominated overland flow on infiltrating surfaces. Ó 2015 Elsevier Ltd. All rights reserved.

1. Introduction The importance of surface microtopography in surface runoff, infiltration, and other hydrologic processes has been emphasized [17,25,26,24,16,13,42,12,7]. The spatial variability of surface microtopography influences overland flow generation [32], delays the initiation of surface runoff [16], and enhances the retention of runoff water [1]. It has been found that surface microtopography controls spatial and temporal variations of overland flow depth and velocity [47,27,18,19]. In addition, greater depths of water ponded in microtopographic depressions tend to increase infiltration [18,37] and unsaturated flow [30]. The spatial variability in surface microtopography also affects surface–subsurface exchange and runoff generation for riparian wetlands [20]. A series of modeling studies have demonstrated the importance of quantifying the actual surface topography and the related hydraulic effects in hydrologic models [47,41,18]. Due to the existence of depressions on topographic surfaces, overland flow often features with discontinuous characteristic and thresholdcontrolled puddle filling and spilling behaviors [12]. Such a topography-dominated, threshold-driven overland flow process varies ⇑ Corresponding author. Tel.: +1 701 231 9758; fax: +1 701 231 6185. E-mail address: [email protected] (X. Chu). http://dx.doi.org/10.1016/j.advwatres.2015.02.004 0309-1708/Ó 2015 Elsevier Ltd. All rights reserved.

spatially and temporally and is referred to as puddle-to-puddle (P2P) dynamics [12]. While the impacts of surface topography on runoff and infiltration processes have been well understood, limited effort has been made to quantitatively characterize such impacts [26,24,16]. In the recent decade, research efforts have been made to conceptually model the hydrologic role of surface microtopography and quantify hydrologic connectivity and the associated dynamic variability in contributing areas. Darboux et al. [15] developed a conditional-walker method to simulate depression filling and water redistribution on rough surfaces, and further evaluate the influence of surface microtopography on discharge. Antoine et al. [4] applied a model similar to the one developed by Darboux et al. [15] and further proposed a functional connectivity indicator (relative surface connection function) that linked surface connection to the filling of surface depression storage. Shaw et al. [38] developed a conceptual model to simulate the fill-spill of depressions and applied it to quantify runoff contributing areas of the wetlands in the Prairie Pothole Region with assumptions of impervious surfaces and uniformly-distributed water input. The aforementioned conceptual models did not take into account the spatio-temporally varying rainfall and infiltration processes. Instead, a steady and uniform rainfall distribution was assumed. In addition, instantaneous water transfer over a

J. Yang, X. Chu / Advances in Water Resources 78 (2015) 80–93

topographic surface was commonly assumed in such conceptual models (e.g., [15,14,4,5,38], which may not be valid for a larger scale hydrologic system. To consider the runoff water transferring time and depression storage, Antoine et al. [3] integrated the connectivity information of subgrids (relative surface connectivity function) into the hydrograph by implementing two corrective procedures (weighted-source and weighted-surface). A uniform slope and 1D sheet flow were assumed in their modeling. Various physically-based distributed models have been widely used to simulate overland flow processes. In such models, the Saint–Venant equations or the simplified forms (diffusion wave and kinematic wave equations) are numerically solved. However, these models are often limited to the modeling of continuous overland flow. Particularly, they may suffer from stability and convergence problems due to the highly non-linear nature of the governing equations and abrupt slopes caused by topographic variations [44,39,23,41]. To avoid the numerical oscillations, many hydrologic models rely on smoothing the digital elevation models (DEMs) by removing depressions [41,29,40]. In addition, such numerical models are computationally intensive and time consuming. For instance, Antoine et al. [3] found that the computation time of a numerical diffusion wave model for a small plot (1  1 m2 and 2-mm DEM resolution) was three orders of magnitude slower (100 days) than that of their simplified conceptual model (10 min). Importantly, those numerical models may not be able to account for the important hydrologic roles of surface depressions and realistically simulate the P2P dynamics and the resulting discontinuous overland flow (e.g., [18,19,40]). Chu et al. [12] developed a conceptual P2P modeling framework, which focused on characterization of surface microtopography and modeling of the dynamic P2P filling, spilling, merging, and splitting processes by assuming instantaneous water transfer. The conceptual model was also applied to quantify the spatio-temporal variability in hydrologic connectivity [46]. This study focuses on improving the process-based overland flow modeling techniques. It is aimed to develop a new hydrotopographic modeling framework to characterize the dynamic P2P filling–spilling–merging–splitting overland flow processes and the associated threshold behaviors under the influence of surface microtopography. The specific objective is to develop a quasi-3D physically-based model that couples a hierarchical overland flow model with a subsurface model. The former is a diffusion wave model for simulating areal surface runoff, while the latter is a modified Green–Ampt model [9] for simulating point (cell based) infiltration and unsaturated flow to account for the infiltration-excess overland flow (Hortonian overland flow) mechanism. Furthermore, testing of the P2P model is performed by using real observed data from overland flow experiments for various microtopography, soil, and rainfall conditions. This study builds on the conceptual P2P modeling framework [12]. However, the P2P model presented herein simulates the physical overland flow processes, and accounts for infiltration and unsaturated flow. To the best of our knowledge, this is the first attempt to simulate microtopographydominated, discontinuous overland flow over infiltrating surfaces in conjunction with the modeling of threshold-controlled, P2P dynamics under varying microtopography, soil, and rainfall conditions.

2. Development of the P2P overland flow model 2.1. Overview of the P2P overland flow modeling framework The P2P modeling framework includes both cell-to-cell (C2C) and P2P flow routing for a set of puddle-based units (PBUs) in a well-delineated, cascaded drainage system [12]. The flowchart of

81

the P2P overland flow model developed in this study is shown in Fig. 1. The model input data include basic simulation parameters, surface DEM, soil hydraulic property parameters, meteorologic data, and puddle delineation results from the PD program [13] (Fig. 1). In the P2P model, a series of nested loops are implemented for water routing, including time loop, basin loop, PBU loop, C2C loop, and P2P loop (Fig. 1). A topographic surface may consist of a number of basins, depending on the surface microtopographic characteristics and the number of outlets at the boundaries of the surface. Within a time step, water routing is performed for all basins (Fig. 1). Overland flow in a basin is essentially controlled by a cascaded P2P drainage system (Fig. 2). This system consists of a series of PBUs determined by surface depressions, which break the continuity and connectivity of the topographic surface (Fig. 2). The PBUs have unique hydrotopographic characteristics and exhibit strong spatial and temporal variability (Fig. 2). Based on the puddle delineation results from the PD program [13], the P2P model tracks the PBUs of each basin and detects their upstream–downstream relationships. PBU is a basic simulation unit within the basin routing loop (Fig. 1). A PBU consists of a number of contributing cells (DEM grids) and a highest-level puddle, which may include a group of lower-level embedded puddles [12] (Fig. 2). The threshold of the highest-level puddle connects this PBU to its downstream PBU (Fig. 2). Overland flow is routed for all PBUs by following their sequences in the welldelineated, cascaded P2P flow drainage system [12] (Fig. 1). Two routing procedures (i.e., C2C and P2P) are implemented for all PBUs (Figs. 1 and 2) [12]. The C2C water routing transfers water from upstream to downstream cells, eventually to the water-ponded cells (i.e., puddle cells), while the P2P water routing simulates the dynamic filling and depleting of puddles and their interactions (i.e., spilling, merging, and splitting) (Fig. 2). The simulation for a basin ends when flow routing is completed for all PBUs in the basin, and the modeling continues until all basins are simulated (Fig. 1). In the P2P model, the water sources of any overland cell include lateral inflow from its upstream cell(s) and rainfall input while the sink terms consist of lateral outflow to its downstream cell, infiltration, and evaporation (Fig. 2). The modified Green–Ampt model [9] is used to simulate infiltration-excess overland flow. Non-uniform and unsteady rainfall and evaporation are considered in the P2P model by specifying a number of rainfall and evaporation zones. The evaporation of ponded water in puddles is estimated based on Allen et al. [2]. 2.2. C2C overland flow routing The PBUs on a topographic surface may have varying sizes and irregular geometric shapes. To facilitate the flow routing on such irregular spatial domains, a DEM-based drainage network is identified for each PBU based on the D8 method [36] (Fig. 3). The drainage network of a PBU is determined by tracking all contributing cells using their flow directions backward from the puddle center or outlet of the PBU. The sequence of the cells in the derived flow drainage network essentially records the contributing relationships from one cell to another cell, or from multiple cells to one cell (Fig. 3). By taking advantage of these contributing relationships of cells in the flow drainage network, the C2C overland flow routing is conducted by simulating water movement from the most upstream cells to water-covered cells in the puddle(s) of the PBU (Fig. 3). Note that the water-ponded puddle cells are excluded in the C2C routing, and that the C2C simulation domain may change with the rising/falling in the ponded water levels of the puddles as more/fewer cells are covered by water (Fig. 3). The P2P routing for a PBU initiates after the C2C routing is completed (Fig. 1).

82

J. Yang, X. Chu / Advances in Water Resources 78 (2015) 80–93

Fig. 1. Flowchart of the P2P overland flow model.

Studies have demonstrated that the diffusion wave model is sufficiently accurate for modeling overland flow and adequate for a variety of practical applications [44,34,22,47,41,43,29,7]. In the P2P model, a 1D diffusion wave model is used for the C2C overland flow modeling for a time-varying spatial domain of each PBU. The

continuity and momentum equations of the 1D diffusion wave overland flow model can be respectively expressed as:

w

  @hðx; tÞ @½Qðx; tÞ þ ¼ w qs;in ðx; tÞ  qs;out ðx; tÞ @t @x

ð1Þ

83

J. Yang, X. Chu / Advances in Water Resources 78 (2015) 80–93

no flow boundary C2C routing for contributing cells

unsteady nonuniform Rainfall

overland flow contributing cells

C2C routing for contributing cells

P1 infiltration

A constant head boundary

P2

ET

threshold

P3 P4 A

puddle filling

puddle spilling

P5 P6

connected area2 (AC2)

AC4

AC5

AC6

B

zero-depth gradient boundary

basin outlet

PBU2

puddle-based unit1 (PBU1)

Fig. 2. Topography-dominated overland flow dynamics, cell-to-cell (C2C) and puddle-to-puddle (P2P) water routing procedures, and boundary conditions [note: P1 and P3 are higher-level merged puddles (e.g., P3 includes P4 and P5); connected area (AC) corresponds to a non-fully-filled puddle (e.g., P4 and AC4); PBU represents a puddle-based unit (e.g., PBU1 and PBU2)].

M

Center of a puddle cell (e.g., N+1)

i

Center of a contributing cell (e.g., i) Center of an outside dummy cell (e.g., M)

N+1

Contributing cell (cc) of a PBU Water ponded puddle cell

Puddle p

Puddle threshold cell Cell in a downstream PBU that receives water spilled from its upstream PBU No flow boundary (i.e., PBU boundary) Constant head boundary (i.e., puddle boundary) Flow direction

x (grid size) Fig. 3. Flow drainage network, contributing cells, water-ponded cells, no flow boundary, and constant head boundary of a puddle-based unit (PBU).

and

@hðx; tÞ ¼ S0 ðxÞ  Sf ðx; tÞ @x

The MacCormack method [31] is used to solve the diffusion wave model [Eqs. (1) and (2)]. Specifically, the predictor procedure

ð2Þ

where h = water depth [L]; x = distance along the flow direction [L]; w = water width [L]; t = time [T]; Q(x, t) = volume flow rate [L3/T]; qs,in (x, t) = inflow rate of source terms [L/T]; qs,out (x, t) = outflow rate of sink terms [L/T]; S0 = surface slope; and Sf = friction slope. The Manning’s equation is applied to calculate flow velocity and account for flow resistance. Rainfall is the only water source in the P2P model while the sink terms include infiltration and evaporation of the ponded water. In addition, the cells that directly connect to the threshold(s) of an upstream puddle (or PBU) may receive the spilled water when the upstream puddle is fully filled (Fig. 2). Thus, the inflow and outflow rates from the source and sink terms are respectively given by:

qs;in ðx; tÞ ¼ rðx; tÞ þ hspðx; tÞ

ð3Þ

and

qs;out ðx; tÞ ¼ eðx; tÞ þ f ðx; tÞ

ð4Þ

where r = rainfall intensity [L/T]; e = evaporation rate [L/T]; f = infiltration rate [L/T]; and hsp = flow rate of the water spilled from the upstream puddle or PBU [L/T].

is implemented to calculate the predicted water depth hði; k þ 1Þ by using the forward finite difference approximations for both time and space. Following the MacCormack method, the predictor finite difference form of Eq. (1) can be expressed as:

Dx

hði; k þ 1Þ  hði; kÞ Q ði þ 1; kÞ  Q ði; kÞ ¼ Dt Dx   þ Dx qs;in ði; k þ 1=2Þ  qs;out ði; k þ 1=2Þ ð5Þ

where Dx = discretized cell size [L]; k = time step; i = cell number; and qs;in ði; k þ 1=2Þ and qs;out ði; k þ 1=2Þ = inflow and outflow rates of source and sink terms for cell i at time k + 1/2, respectively. Solving Eq. (5) yields the predicted overland flow depth hði; k þ 1Þ:

Dt½Qði þ 1; kÞ  Q ði; kÞ hði; k þ 1Þ ¼ hði; kÞ  D x2   þ Dt qs;in ði; k þ 1=2Þ  qs;out ði; k þ 1=2Þ

ð6Þ

Q (x, t) is calculated by using the Manning’s equation:

Qðx; tÞ ¼

5 1 Dx ½hðx; tÞ3 ½Sf ðx; tÞ2 nðxÞ

ð7Þ

84

J. Yang, X. Chu / Advances in Water Resources 78 (2015) 80–93

where n(x) = Manning’s roughness coefficient. The friction slope Sf is calculated using the momentum equation [Eq. (2)]:

@hðx; tÞ Sf ðx; tÞ ¼ S0 ðxÞ  @x

ð8Þ

Q (i, k) and Q (i + 1, k) in Eq. (5) are given by incorporating the backward finite difference form of the momentum equation [Eq. (2)] into Eq. (7):

Q ði;kÞ ¼

 1 5 Dx hði;kÞ  hði  1;kÞ 2 ½hði;kÞ3 S0 ðiÞ  ¼ Q up-in ði;kÞ nðiÞ DLi ;i1

ð9Þ

where DL = distance between the centers of two cells [L], which is pffiffiffi equal to Dx if flow is along x or y direction, and 2Dx if flow is along a diagonal direction (Fig. 3); Qup-in(i, k) = upstream inflow of cell i at time step k [L3/T]; and Qout(i, k) = outflow of cell i at time step k [L3/ T]. Substituting Eqs. (9) and (10) into Eq. (6) gives:

ð11Þ

According to the DEM-based D8 flow drainage network, cell i may receive water from a number of upstream cells (Fig. 3). The total inflow of cell i, Qin(i, k), can be calculated as the summation of the outflow Qout from its upstream contributing cells as follows:

Q in ði; kÞ ¼

nncc Xi

Q out ðj; kÞ

ð12Þ

j¼1

where nncci = number of upstream cells that contribute water to cell i. Thus, based on the D8 flow drainage network, Eq. (11) can be further expressed as:

Dt½Q in ði; kÞ  Q out ði; kÞ hði; k þ 1Þ ¼ hði; kÞ þ Dx2   þ Dt qs;in ði; k þ 1=2Þ  qs;out ði; k þ 1=2Þ

ð13Þ

For the corrector procedure, a backward finite difference method is implemented for both time and space. As shown below, the time increment used in the corrector procedure is Dt/2.

Dt Qði; k þ 1Þ  Q ði  1; k þ 1Þ 2 Dx2   Dt þ q ði; k þ 3=4Þ  qs;out ði; k þ 3=4Þ 2 s;in

hði; k þ 1Þ ¼ hði; k þ 1=2Þ 

ð14Þ

in which qs;in ði; k þ 3=4Þ and qs;out ði; k þ 3=4Þ are the source and sink terms

at

time

point

k + 3/4,

respectively.

Q ði; k þ 1Þ

and

Q ði  1; k þ 1Þ are calculated by incorporating the forward finite difference approximation of the momentum equation [Eq. (2)] as follows:

" #12 i53 Dx h hði þ 1; k þ 1Þ  hði; k þ 1Þ Q ði; k þ 1Þ ¼ hði; k þ 1Þ S0 ðiÞ  nðiÞ DLi ;iþ1 ¼ Q out ði; k þ 1Þ ð15Þ " i53 Dx h Q ði  1;k þ 1Þ ¼ hði  1;k þ 1Þ S0 ði  1Þ nði  1Þ #12 hði;k þ 1Þ  hði  1; k þ 1Þ  ¼ Q upin ði;k þ 1Þ DLi ;i1

Dt Q out ði; k þ 1Þ  Q upin ði; k þ 1Þ 2 Dx2  Dt  þ qs;in ði; k þ 3=4Þ  qs;out ði; k þ 3=4Þ 2

hði; k þ 1Þ ¼ hði; k þ 1=2Þ 

ð17Þ

Similar to the predictor procedure, the total inflow of a cell can be calculated based on the D8 flow drainage network (Fig. 3):

Q in ði; k þ 1Þ ¼

nncc Xi

Q out ðj; k þ 1Þ

ð18Þ

j¼1

 1 5 Dx hði þ 1; kÞ  hði; kÞ 2 ½hði þ 1; kÞ3 S0 ði þ 1Þ  Q ði þ 1; kÞ ¼ nði þ 1Þ DLi ;iþ1 ¼ Q out ði; kÞ ð10Þ

Dt½Q upin ði; kÞ  Q out ði; kÞ hði; k þ 1Þ ¼ hði; kÞ þ Dx2   þ Dt qs;in ði; k þ 1=2Þ  qs;out ði; k þ 1=2Þ

Thus, Eq. (14) can be rewritten as:

hði; k þ 1=2Þ is given by:

hði; k þ 1=2Þ ¼

i 1h hði; kÞ þ hði; k þ 1Þ 2

ð19Þ

Hence, the final expression for the water depth at time step k + 1 for the corrector procedure is given by:

( 1 Q ði; k þ 1Þ  Q out ði; k þ 1Þ hði; k þ 1Þ ¼ hði; kÞ þ hði; k þ 1Þ þ Dt in 2 Dx2 )   ð20Þ þDt qs;in ði; k þ 3=4Þ  qs;out ði; k þ 3=4Þ

2.3. Initial and boundary conditions In the C2C modeling, all contributing cells in PBUs are assumed to be dry initially:

hðx; tÞjt¼0 ¼ 0

x 2 Xcc

ð21Þ

where Xcc = contributing cells in the spatial domain for the diffusion wave modeling (Fig. 3). Based on Eqs. (13) and (20) and the initial condition Eq. (21), the water depths at the first time step for the predictor and corrector procedures can be respectively expressed as:

  hði; 1Þ ¼ Dt qs;in ði; 1=2Þ  qs;out ði; 1=2Þ hði; 1Þ ¼

ð22Þ

( )   1 Q ði; 1Þ  Q out ði; 1Þ hði; 1Þ þ Dt in þ D t q ði; 3=4Þ  q ði; 3=4Þ s;in s;out D x2 2

ð23Þ Three types of boundary conditions are considered in the P2P overland flow model to handle the interactions between the C2C modeling domain and the outside surroundings, including the water-covered puddle cells. These boundary conditions include: (1) no flow boundary for the most upstream cells of all PBUs (Fig. 3), (2) constant head boundary for the C2C routing cells that directly connect to the water-ponded puddle cells (Fig. 3), and (3) zero-depth gradient boundary for the outlet cells. The most upstream boundary cells of the drainage network of a PBU have a closed boundary (Fig. 3). The inflow of the boundary cells from the outside surrounding is zero:

Qðx; tÞ ¼ 0

x 2 XNFB

ð24Þ

where XNFB = no flow boundary domain (Fig. 3). Eq. (24) implies that a boundary cell N (N 2 XNFB ) receives zero flow from its outside upstream cell M (i.e., outflow of dummy cell M) at time step k:

Q out ðM; kÞ ¼ 0

ð25Þ

The constant head boundary is applied for water transfer from the C2C contributing cells to the water-ponded puddle cells (Fig. 3). A constant head boundary cell N is linked to its downstream puddle cell N + 1. The water elevation of puddle cell N + 1 is equal to the uniform water level of its corresponding puddle:

ð16Þ

hðx; tÞ þ zðxÞ ¼ zhpðp; tÞ

x 2 XCHB ðp; tÞ

ð26Þ

J. Yang, X. Chu / Advances in Water Resources 78 (2015) 80–93

hðN þ 1; kÞ þ zðN þ 1Þ ¼ zhpðp; kÞ

N þ 1 2 XCHB ðp; kÞ

ð27Þ

where zhp(t) = water elevation of puddle p at time t [L]; z(x) = elevation at x [L]; XCHB(p, k) = constant head boundary domain at time step k for puddle p; and h(N + 1, k) = water depth of cell N + 1 at time step k [L] (Fig. 3). Note that the water level of a puddle or its water-ponded cells is assumed constant within one time step in the C2C water routing procedure, but it may change over time due to the dynamic P2P process. The outflow boundaries (a set of outlet cells delineated by the PD program [13] are located at the most downstream ends of the entire drainage system. The zero-depth gradient boundary is applied to these outlet cells:

@h ¼0 @x

x 2 XZDGB

hðN þ 1; kÞ ¼ hðN; kÞ

ð28Þ N 2 XZDGB

ð29Þ

where XZDGB = zero-depth gradient boundary domain; N = outlet boundary cell; and N + 1 = outside dummy cell. The aforementioned numerical scheme for solving the diffusion wave equations is stable if the Courant–Friedrichs–Lewy criterion is satisfied [28,29].

max

8 h pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffii9

Dx

;

1

i 2 Xcc

ð30Þ

where g = gravitational acceleration [L/T2]. For traditional numerical models, the computation time is long because the number of iterations for each time step is calculated as the maximum iterations for all cells within the spatial domain, and is often determined by the water-ponded puddle cells that have high flow velocity and/or water depth. In the P2P model presented herein, however, the number of iterations for the C2C routing of a PBU for each time step is determined for a limited number of cells in the C2C routing domain that excludes any puddle cells. This approach significantly improves the modeling efficiency and the computation cost. 2.4. P2P overland flow routing The P2P water routing links to the C2C modeling through the constant head boundaries and involves simulations of the puddle filling, spilling, merging, and splitting processes after water routing from contributing cells to the water-covered puddle cells of a PBU (Figs. 1 and 2). A PBU may include one first-level puddle (e.g., puddle P6 in PBU2, Fig. 2) or a group of embedded puddles at multiple levels (e.g., puddles P1–P5 in PBU1, Fig. 2). The P2P water routing consists of two major steps: (1) creating/updating a puddle routing list and (2) conducting the P2P water routing for puddles on the puddle routing list. The dynamic P2P routing for a time step may not involve all puddles of a PBU. In the P2P model, the puddle routing list for a time step is utilized to regulate the puddles that are ‘‘active’’ in the P2P processes. Whether or not a puddle is involved in the P2P routing at a time step depends on its filling condition and the relationships with its neighboring puddles. For instance, if a PBU has a single-level puddle (e.g., puddle P6 in PBU2, Fig. 2), this puddle is added to the puddle routing list. For a group of multilevel puddles, whether or not the highest-level puddle (e.g., puddle P1, Fig. 2) and its embedded puddles (e.g., puddles P2–P5 of PBU1, Fig. 2) are included in the puddle routing list depends on their filling status denoted by a puddle ponding index. Three puddle ponding index values (0, 1, and 100) are used to represent the unfilled (e.g., puddle P1, Fig. 2), partially-filled (e.g., puddle P2, Fig. 2), and fully-filled conditions (e.g., puddle P6, Fig. 2), respectively.

85

The P2P routing involves a sequence of puddle filling, spilling, merging, and splitting for all puddles on the puddle routing list. The puddle routing list is a dynamic one, which is updated by adding or removing puddles to/from the list after completing water routing due to their spilling, merging, and splitting processes (Fig. 1). The P2P routing for a puddle on the puddle routing list depends on the amount of available water in this puddle and the relationships with its neighboring puddles. First, the total volume of water available to fill the puddle (Fig. 1) is calculated based on: (1) the water ponded in the puddle at the beginning of current time step, (2) the water loss or gain from the source and sink terms, and (3) inflow from the upstream contributing cells simulated by the C2C water routing. Based on this calculated total water volume and the stage-area-volume relationship of the puddle, P2P modeling is performed for the following three conditions (Fig. 1): Condition 1: If the total available water volume of the puddle being routed is greater than zero and smaller than its maximum depression storage, a puddle filling process occurs (Fig. 1, puddle P2 in Fig. 2). Based on the stage-area-volume relationship of the puddle, the inundation area and the corresponding water-ponded puddle cells are identified and the level of ponded water of this puddle and the associated cells is determined. It is assumed that all water-ponded puddle cells have a uniform water level (e.g., puddle P2 in Fig. 2). Condition 2: Due to infiltration and evaporation of the ponded water, an individual puddle can be dry. In this case, a higher-level puddle that consists of a number of lower-level embedded puddles will split into its lower-level puddles (puddle splitting process) if the net volume of available water is negative (Fig. 1). In the P2P model, the ponding index of this higher-level puddle shifts to 0 (i.e., unfilled condition). The lower-level embedded puddles are added to the puddle routing list, and corresponding changes also are made to their ponding index values. For example, due to water losses, a splitting process is trigged for puddle P3 (Fig. 2), forming two lower-level embedded puddles P4 and P5 (Fig. 2). The puddle splitting process occurs only when the water level of a merged puddle drops below the threshold of its embedded puddles (e.g., the threshold point A of puddles P4 and P5 in Fig. 2). Condition 3: If the volume of water available for the puddle being routed is greater than its maximum depression storage, three cases that may involve both spilling and merging processes occur (Fig. 1): (1) if the puddle being routed is a single level puddle or a puddle at the highest level, the excess water spills through its threshold(s) to the downstream cell(s) (Fig. 1) (e.g., puddle P6 in Fig. 2); (2) if the puddle being routed (e.g., puddle P4 in Fig. 2) shares a common threshold with another puddle that is not fully filled (e.g., P5, Fig. 2), the excess water spills to that puddle; and (3) if the puddle being routed shares a common threshold with another puddle that is already fully filled, the two puddles merge and form a higher-level puddle. This merged puddle is then added to the puddle routing list and its ponding index is changed accordingly. After each puddle in the puddle routing list completes its water routing, the puddle routing list is updated until all puddles in the list are routed (Fig. 1). Within one time step, the water level in a puddle may rise or drop dramatically, and the water routing process generally involves puddles at different levels. If the highest level puddle in a PBU reaches the fully-filled condition (e.g., P1, Fig. 2), the excess water spills to its downstream PBU (PBU2, Fig. 2), and then water routing for next PBU initiates (Fig. 1). Following the cascaded overland flow drainage pattern, such a modeling procedure is performed for all PBUs. The simulation for a basin ends when all of the PBUs in the basin complete their flow routing (Fig. 1). Then, the modeling continues for next basin until all basins complete their P2P water routing (Fig. 1). The total discharge from

86

J. Yang, X. Chu / Advances in Water Resources 78 (2015) 80–93

the entire surface is calculated as the summation of the discharges from all basins. A surface may have some special topographic features (e.g., flats and multi-threshold puddles). In addition to water routing for regular puddles, special schemes also are developed and incorporated in the P2P model to cope with the water routing for these special topographic features. In the case of a multi-threshold puddle, the filling process is similar to that for regular puddles. However, after this puddle is fully-filled, the excess water is evenly distributed to all the thresholds of the puddle and spills to their downstream cells or puddles. A flat is treated as a ‘‘flat puddle’’ that is able to transfer runoff water, but has zero depression storage (i.e., it cannot store any water) [12]. The P2P model provides modeling details on the dynamic C2C and P2P processes, puddle filling status, ponded water distributions, discharges at basin outlets, and mass balance analyses for all time steps. 2.5. Modified Green–Ampt model for simulating infiltration and unsaturated flow Infiltration and unsaturated flow are simulated for all cells in a PBU by using a modified Green–Ampt model [9] (Fig. 1). The modified Green–Ampt model incorporated in the P2P model simulates infiltration into layered soils of arbitrary initial soil moisture distributions and rainfall excess under complex rainfall patterns. Two distinct periods, pre-ponding and post-ponding, are taken into account in the infiltration modeling. The model tracks the movement of wetting front along each soil profile, checks the ponding status, and handles the shift between ponding and non-ponding conditions. In addition, the model simulates vertical soil water drainage and redistribution for dry time periods between rainfall events [10]. 2.6. Modeling of time-varying hydrologically-connected areas Depressions break the connectivity of a topographic surface and generate a number of hydrologically localized areas, which are referred to as ‘‘connected areas’’ (ACs) (Fig. 2). At the initial stage of the rainfall–runoff process, all cells within an AC drain runoff water to the corresponding puddle (e.g., AC2 in Fig. 2). After the puddle is fully filled, it spills the excess water to its downstream cell(s) or puddle(s), which develops the connection between the two ACs, inducing expansion or evolution of ACs. Any hydrologic connection of an AC to the outlet of the surface may result in a stepwise increase in the outlet discharge. Note that an AC is different from a PBU in two aspects: (1) PBUs are a series of ‘‘static’’ topographic domains serving as basic units in the P2P overland flow routing, while ACs are time-varying contributing areas associated with unfilled or partially-filled depressions and basin outlets; and (2) PBUs include a highest-level puddle and probably a set of embedded puddles (e.g., puddles P1–P5 for PBU1 in Fig. 2), while each AC relates to one unfilled or partially-filled puddle regardless of the level of the puddle. In the P2P model, identification of initial ACs (t = 0) starts from the basin outlet(s) and the centers of all first-level puddles based on the flow directions determined by the PD program [13]. These DEM-determined ACs may change over time, depending on the surface topographic characteristics and water sources and sinks of the system. For any time steps, such time-varying ACs are determined starting from the basin outlets and unfilled or partiallyfilled puddles by adding: (1) their contributing cells and (2) any upstream ACs that are fully filled and have established connection to the AC being identified. These time-varying ACs are utilized to analyze hydrologic connectivity and the related threshold behaviors.

2.7. Testing of the P2P overland flow model The conceptual version of the P2P model has been tested using experimental data for impervious surfaces to demonstrate its capability in simulating the P2P dynamics [12] and quantifying the related hydrologic conductivity [46]. In this study, the P2P model was compared against the analytical solution to a kinematic wave overland flow model [45] for a smooth hillslope surface. The area of the smooth hillslope surface was 60.0  30.0 m2 and its slope was 5.0%. A steady and uniform rainfall with an intensity of 7.62 cm/hr and a duration of 30 min was applied in the modeling. Few existing studies combined overland flow experiments with modeling of the P2P overland flow dynamics for infiltrating soil surfaces of various depressions. In this study, 12 surfaces (S1–S6 and RR1–RR6, Fig. 4) were used to test the new P2P model. These surfaces were created in the field and laboratory settings to represent various microtopographic characteristics and soil surface roughness under these two conditions. Three laboratory overland flow experiments (Experiments 1–3) were used to verify the capability of the P2P model in simulating the microtopographyinfluenced, dynamic P2P filling–spilling–merging–splitting processes for infiltrating soil surfaces under various rainfall conditions. Table 1 shows the basic information on the three laboratory overland flow experiments. Experiments 1 and 2 were conducted using two soil surfaces with distinct microtopographic characteristics (rough and smooth, Fig. 4a and b) under the same soil (silty clay: 11.9% sand, 42.0% silt, and 46.1% clay) and initial moisture condition (0.213 cm3/cm3). The effective hydraulic conductivity and capillary head for this soil were 0.29 cm/hr and 80 cm, respectively. Similar unsteady rainfall events with the same duration (Table 1) were applied for these two experiments. A much larger soil surface (Fig. 4c), different soil type (silty clay loam: 11.1% sand, 60.1% silt, and 28.8% clay) and initial soil moisture content (0.145 cm3/cm3), and a steady rainfall event were applied for Experiment 3 (Table 1). The effective hydraulic conductivity and capillary head of this silty clay loam soil were 0.30 cm/hr and 42 cm, respectively. The initial soil moisture contents were measured for soil samples as per ASTM D-2216-10 [6]. Surfaces S1–S3 were scanned by using an instantaneous-profile laser scanner and their high-resolution DEMs were acquired. The horizontal and vertical resolutions of the original scanned data were 0.98 and 0.5 mm, respectively. Lower resolution DEMs (e.g., 10 mm) used for the modeling were then created by using the Kriging method. A four-head Norton-style rainfall simulator was utilized to generate rainfall for the three experiments. To calibrate the rainfall simulator, the area covered by the simulator was divided into 20  20 cm2 grids and rainfall intensities were measured for all grids before the experiments to account for the spatial variations of rainfall although the variations were small. The non-uniformly distributed rainfall was then used in the P2P modeling by specifying a number of rainfall zones. The rainfall intensities selected for the overland flow experiments in this study ranged from 1.15 to 5.9 cm/hr (Table 1). The selection of the rainfall intensities was primarily based on the purposes of the experiments, soil types, overland flow generation potentials, and potential adverse impacts on surface microtopography and soil hydraulic properties (e.g., raindrop effect, soil erosion, and surface sealing). The data collected from the experiments included: outlet discharges (measured directly), wetting front depths (recorded from all the transparent plexiglass sides of the soil boxes), and a set of critical times (e.g., ponding time, and puddle spilling and merging times). P2P modeling was performed by using the same experimental data (e.g., DEM, soil, and meteorologic data) and the simulation results were compared with the observed data.

87

J. Yang, X. Chu / Advances in Water Resources 78 (2015) 80–93

(a) S1

(b) S2

(d) S4

(g) RR1

(c) S3

(f) S6

(e) S5

(h) RR2

(j) RR4

(k) RR5

(i) RR3

(l) RR6

Fig. 4. DEMs of 12 topographic surfaces used for P2P overland flow modeling.

Table 1 Description of three laboratory overland flow experiments. Experiment

Surface

Surface size (m2)

Duration (min)

Rainfall intensity (cm/hr)

Initial water content (cm3/cm3)

Soil type

1 2 3

S1 S2 S3

1.0  0.6 1.0  0.6 3.0  1.2

92 92 80

Unsteady: 5.90 (25 min) – 1.15 (43 min) – 3.56 (16 min) – 5.90 (8 min) Unsteady: 6.00 (25 min) –1.17 (43 min) – 3.52 (16 min) – 6.00 (8 min) Steady: 3.23 (60 min)

0.213 0.213 0.145

Silty clay Silty clay Silty clay loam

2.8. Analyses of overland flow dynamics and infiltration for various microtopographic and rainfall conditions To better understand the combined effects of surface microtopography and rainfall characteristics on the dynamic P2P overland flow processes accounting for infiltration, the P2P model was applied to nine selected surfaces that had different microtopographic characteristics and various rainfall conditions. The surfaces included: (1) two field plots with an area of 3.5  6.0 m2 and an overall slope of 2.5% (S4 and S5, Fig. 4d and e); (2) a smooth hillslope surface (S6, Fig. 4f), which had the same area and slope as those of surfaces S4 and S5; and (3) six random roughness surfaces with an area of 0.6  2.0 m2 (RR1–RR6, Fig. 4g–l). The random

roughness values for the six RR surfaces increased from 0.34 cm for RR1 to 1.5 cm for RR6. Surfaces S4–S5 and RR1–RR6 were scanned by using the instantaneous-profile laser scanner to obtain their high-resolution DEMs. The statistical descriptions of the six random roughness surfaces were detailed by Chu et al. [11] and Chi et al. [8]. In the P2P modeling, a steady uniform rainfall event with an intensity of 2.2 cm/hr (Rain 1, Fig. 5) was firstly applied to surfaces S4–S6 and RR1–RR6 to examine the effects of surface microtopography on the P2P overland flow processes. Then, three unsteady rainfall events (Rain 2–4, Fig. 5) were applied to surfaces S4–S6 to discuss the influences of temporal variability in rainfall on the P2P dynamics. Note that the three unsteady rainfall events (Rain

88

J. Yang, X. Chu / Advances in Water Resources 78 (2015) 80–93

1), there was no ponded water on the surfaces, and all rainwater infiltrated into soil. Gradually, the rainfall intensity (5.90 cm/hr for S1) exceeded the infiltration capacity, and surface ponding occurred and outlet discharge initiated (Fig. 7a). The discharge dropped to zero due to the significant decrease in rainfall intensity (from 5.90 to 1.15 cm/hr for S1) (stage 2, Fig. 7a), and the connections between puddles for surface S1 were broken as the water levels in puddles dropped below their thresholds due to infiltra-

16 Rain1 Rain2 Rain3 Rain4

12 10 8 6 4

160

2

140 0.4

0.6

0.8

1.0

Time (hr) Fig. 5. Four rainfall events used for P2P modeling (Rain 1: steady rain with an intensity of 2.2 cm/hr; Rain 2–4: unsteady rain with the same cumulative rain as Rain 1).

Table 2 Modeling scenarios for the field surfaces S4–S6 and random roughness surfaces RR1– RR6 under Rain 1–4 (Rain 1: steady rainfall with an intensity of 2.2 cm/hr; Rain 2–4: unsteady rainfall with the same cumulative rainfall as that of Rain 1). Rainfall

Initial water content

Soil type

S4–S6 RR1–RR6

6.0  3.5 2.0  0.6

60

Rain 1–4 Rain 1

0.32

Clay loam

2–4) had the same cumulative rainfall as that of the steady rainfall event (Rain 1). Clay loam soil with an initial moisture content of 0.32 was used in these simulations (Table 2). The rainfall duration was 60 min for all the simulations (Table 2). 3. Results and discussion 3.1. P2P model vs. analytical kinematic wave model Fig. 6 shows the comparison of the hydrographs simulated by the analytical kinematic wave model [45] and the P2P model for the impervious hillslope surface. Good agreement has been achieved. As shown in Fig. 6, both hydrographs display an initial significant increase (rising limb) and reach an equilibrium state with a constant discharge of 137.15 m3/hr (equilibrium stage), which equals the rate of the rainfall input. As rainfall ceases at t = 30 min, the discharge decreases dramatically and approaches to zero (falling limb). The comparison results indicate that the P2P model is able to accurately simulate overland flow for contributing cells (i.e., C2C water routing). Fig. 6 also shows the comparison of the hydrographs from the P2P model for the impervious and infiltrating surfaces. Due to infiltration, a delayed initiation of surface runoff (rising limb), a lower steady outlet discharge, and a steeper falling limb can be observed. 3.2. P2P modeling vs. laboratory overland flow experiments Fig. 7 shows the comparison between the observed and simulated hydrographs for experiments 1–3 (Table 1). To evaluate the goodness of fit, the Nash–Sutcliffe modeling efficiency coefficient [35] was calculated for each set of the observed and simulated discharge data. The modeling efficiency values were 0.95, 0.89, and 0.99 for experiments 1–3, respectively. Four stages can be identified in the hydrographs of surfaces S1 and S2 (Fig. 7a and b) corresponding to the changes in rainfall intensity. Initially (stage

60 40

0

10

20

30

40

Time (min) P2P overland flow model for impervious surface Analytical solution for impervious surface P2P overland flow model for infiltrating surface Fig. 6. Comparisons of the hydrographs simulated by the P2P model and the analytical kinematic wave model for both impervious and infiltrating surfaces.

500 Discharge (cm 3 /min)

Duration (min)

80

0

0 stage 4

400 rainfall simulated-S1 observed-S1

300 stage 1

200

5 stage 3 10

100 stage 2

(a)

0

15 0

20

40

60

80 0

500 Discharge (cm 3/min)

Surface area (m2)

100

20

stage 4

400 simulated-S2 observed-S2 observed-S1

300 stage 1

200

5 stage 3 10

100 stage 2

(b) 15

0 0

20

40

60

80

1200 Discharge (cm 3 /min)

Surface

120

Rainfall intensity (cm/hr)

0.2

Rainfall intensity (cm/hr)

0.0

Discharge (m3 /hr)

0

0

1000

5

observed-S3 simulated-S3 rainfall

800

10

600

15

400

20

200

25 (c)

0

Rainfall intensity (cm/hr)

rainfall intensity (cm/hr)

14

30 0

20

40 60 Time (min)

80

Fig. 7. Comparison of the observed and simulated hydrographs for surfaces S1–S3.

89

J. Yang, X. Chu / Advances in Water Resources 78 (2015) 80–93

ponding time spilling time respilling time

Simulated time (min)

80 60

45o line

40 20 (a) S1 0 0

40 60 Observed time (min)

80

45o line

ponding time spilling time

80 Simulated time (min)

20

60

40

20 (b) S3 0 0

20

40

60

80

Observed time (min) Fig. 8. Comparison of the observed and simulated critical times (puddle ponding, spilling, and re-spilling times) for all selected puddles of surfaces S1 and S3.

tion. The hydrograph of experiment 1 started to rise again when the rainfall intensity changed to 3.56 cm/hr for S1 (stage 3, Fig. 7a). The discharge showed a remarkable stepwise increase when the rainfall intensity increased to 5.90 cm/hr for S1 (stage 4, Fig. 7a). A similar changing pattern and the four stages can be observed for experiment 2 (Fig. 7b). The comparisons of the hydrographs for surfaces S1 and S2 indicate that the smooth surface (S2, Fig. 4) for experiment 2 had a faster increase in hydrograph and yielded more surface runoff during the initial period of stages 1 and 3 (Fig. 7a and b). This finding also has been reported by other researchers in their laboratory and field experimental studies (e.g., [48,33,25,21]. The difference in runoff generation between the smooth and rough surfaces diminished gradually with persistent rainfall for the later period of stages 1 and 3 (Fig. 7b), which was consistent with the findings from Moore and Singer [33] and Gomez and Nearing [21]. The final observed

and simulated discharges at stage 4 were significantly higher than those at stage 1 for both surfaces (S1 and S2, Fig. 7a and b) although the two stages had the same rainfall intensity (5.90 cm/ hr). This can be attributed to the decrease in soil infiltrability at the later stage as more soil was fully saturated and the wetting front moved to deeper soil. From Fig. 7b, it also can be observed that surface S1 had a ‘‘steeper’’ increase in the hydrographs at stages 1 and 3 than surface S2 due to the threshold behaviors (spilling process) of several dominated puddles (Fig. 4a). The observed and simulated hydrographs for surface S3 exhibit a strong stepwise increasing pattern (Fig. 6c), which can be attributed to its unique microtopographic characteristics (e.g., distribution and organization of puddles and their geometric properties) (Fig. 4c) and the threshold flow associated with the puddle filling–spilling–merging dynamics. When all puddles were fully filled, the entire surface contributed runoff water to the outlet through the well-connected drainage system, which led to a significant increase in the hydrograph (Fig. 7c). These modeling details and the comparisons against the observed data demonstrate the ability of the P2P model to accurately simulate the timing and magnitude of the response of the system to various rainfall inputs and capture the microtopography-influenced P2P dynamics and overland flow generation processes. Fig. 8 shows the simulated and observed ponding, spilling, and re-spilling times (i.e., critical P2P times) for the dominated puddles on surfaces S1 and S3 during experiments 1 and 3. These puddles showed different behaviors in the P2P processes due to a number of controlling factors. Basically, the simulated critical times match the observed ones for both experiments (Fig. 8). The mass balance analyses for the three experiments and the corresponding simulations are summarized in Table 3. The average relative errors of the cumulative surface runoff (outlet discharge), cumulative infiltration, and depression storage for the three experiments were 3.85%, 1.06%, and 0.00%, respectively (Table 3). Infiltration was the dominant process, accounting for 76.87%, 74.89%, and 76.01% of the total rainfall input for experiments 1–3, respectively (based on the simulated values, Table 3). The mass balance analyses also confirmed that the smooth surface S2 generated more surface runoff (24.95% of the total rainfall) than the rough surface S1 (20.58% of the total rainfall). The difference between the simulated cumulative outlet discharges for surfaces S1 and S2 (1.40  103 cm3) was much greater than the storage loss from surface S1 (0.75  103 cm3), indicating that the rough surface S1 with depressions increased the infiltration rates. The simulated runoff water retained on the contributing area was not significant for the three experiments (<100 cm3, Table 3). By the end of the simulation period, the average ponded water depths for the puddle/depression

Table 3 Mass balance analyses for experiments 1–3 and comparison with the simulations. Observed (103 cm3)

Simulated (103 cm3)

Experiment

Mass balance terms

1 (S1)

Cumulative rainfall Cumulative runoff Cumulative infiltration Depression storage Water on contributing area

30.13 6.64 22.73 0.75 –

30.13 6.20 23.16 0.75 0.02

Relative error (%) 0.00 6.63 1.89 0.00 –

2 (S2)

Cumulative rainfall Cumulative runoff Cumulative infiltration Depression storage Water on contributing area

30.46 7.78 22.68 0.00 –

30.46 7.60 22.81 0.00 0.04

0.00 2.31 0.57 0.00 –

3 (S3)

Cumulative rainfall Cumulative runoff Cumulative infiltration Depression storage Water on contributing area

154.90 28.36 118.58 7.96 –

154.90 29.10 117.74 7.96 0.09

0.00 2.61 0.71 0.00 –

90

J. Yang, X. Chu / Advances in Water Resources 78 (2015) 80–93

(a) S4

(c) RR4

(b) S5

Fig. 9. Spatial distributions of the simulated flow depths (h) for surfaces S4, S5, and RR4.

Discharge (cm3 /min)

5000 S4 S5 S6

4000 3000 2000

(a)

1000 0 0.0

0.2

0.4

0.6

0.8

1.0

16000 RR1 RR2 RR3 RR4 RR5 RR6

Discharge (cm3 /min)

14000 12000 10000 8000 6000 4000

(b)

2000 0 0.2

0.4

0.6

0.8

1.0

Time (hr) Fig. 10. Simulated hydrographs for surfaces S4–S6 and RR1–RR6 under Rain 1 (steady uniform rainfall with an intensity of 2.2 cm/hr).

areas were 7.82 and 13.05 mm for surfaces S1 and S3, respectively. The simulated flow velocities were 2.26, 2.17, and 1.59 cm/s for surfaces S1–S3 (experiments 1–3), respectively. The modeling for experiments 1–3 demonstrated the capability of the P2P model in simulating the P2P dynamics, outlet hydrographs, infiltration, and other hydrologic processes under the influence of surface microtopography for different rainfall events.

3.3. Combined effects of rainfall and surface microtopography on rainfall partition and unsaturated flow Fig. 9 shows the spatial distributions of the simulated overland flow depths at t = 1.0 h for surfaces S4, S5, and RR4 under Rain 1. The simulated overland flow depths for surface S4 exhibited a

strong non-uniform distribution. Runoff water was mainly stored in the large depressions (Fig. 9a). Following the cascaded drainage system determined by surface microtopography, runoff water drained from upstream to downstream puddles (Fig. 9a). As shown in the delineation results for the field plot surface S5 (Fig. 4e), it was featured with many small well-connected depressions, forming a series of micro-rills (Fig. 9b). The simulated overland flow depths for surface RR4 had a scattered and random distribution (Fig. 9c). Fig. 10 shows the simulated hydrographs for surfaces S4–S6 and RR1–RR6 under Rain 1. Although surfaces S5 and S6 had different microtopographic characteristics, drainage patterns, and water distributions of overland flow depths, their hydrographs were very similar. Their cumulative runoff values were 23.51% and 23.29% of the total rainfall input, respectively (Table 4). This can be attributed to their relatively ‘‘smooth’’ microtopography, compared with the rougher surface S4. The similar hydrographs obtained in the simulations raise the question on the level of details one needs to go into when characterizing a surface, an issue out of the scope of the present study. The rough field plot surface S4 generated much less surface runoff, 14.90% of the total rainfall due to its larger depression storage (43.79  103 cm3, 9.48% of the total rainfall input, Table 4). The hydrographs of surfaces S5 and S6 (Fig. 4e and f) exhibited a smooth increasing pattern while the rough surface S4 (Fig. 4d) showed a strong stepwise increasing pattern because of the dynamic P2P processes and the associated threshold behaviors (Fig. 10a). The hydrographs of the six random roughness surfaces RR1–RR6 displayed variable increasing patterns (Fig. 10b), which can be attributed to the differences in their microtopographic characteristics (e.g., sizes of soil aggregations and maximum depression storage). The changing patterns of the simulated hydrographs for these nine surfaces (S4–S6 and RR1– RR6) confirmed that overland flow generation was strongly affected by surface microtopography, featuring a threshold-controlled, dynamic expansion of the areas that connected and contributed runoff water to the outlet(s). Fig. 11 shows the relationships of random roughness with maximum depression storage and cumulative runoff for surfaces RR1– RR6. As expected, a rougher surface with a higher random roughness index had greater maximum depression storage, and consequently generated less surface runoff (Fig. 10b). The changing pattern of a hydrograph largely depended on puddles, their organizations or relationships, and their microtopographic properties (e.g., maximum depression storage). To examine the effects of unsteady rainfall on the P2P overland flow processes and unsaturated flow (e.g., wetting front movement), a series of simulations were conducted for surfaces S4–S6

91

J. Yang, X. Chu / Advances in Water Resources 78 (2015) 80–93 Table 4 Summarized simulation results for surfaces S4–S6 and Rain 1–4.a

a

Surface

Mass balance terms

Rain 1

Rain 2

Rain 3

Rain 4

S4

Cumulative rainfall (103 cm3) Depression storage (103 cm3) Water on contributing area (103 cm3) Cumulative infiltration (103 cm3) Cumulative runoff (103 cm3) Mean depth of wetting front (cm) Coefficient of variation of wetting front depths

462.00 43.79 0.95 348.42 68.84 7.26 0.01

462.00 27.98 0.00 254.95 179.07 5.34 0.13

462.00 40.07 0.00 301.08 120.85 6.28 0.09

462.00 39.58 0.00 235.35 187.07 4.85 0.04

S5

Cumulative rainfall (103 cm3) Depression storage (103 cm3) Water on contributing area (103 cm3) Cumulative infiltration (103 cm3) Cumulative runoff (103 cm3) Mean depth of wetting front (cm) Coefficient of variation of wetting front depths

462.00 2.89 2.45 348.04 108.62 7.25 0.00

462.00 0.23 0.00 241.44 220.33 5.05 0.05

462.00 2.56 0.00 290.58 168.85 6.06 0.04

462.00 1.61 0.00 232.34 228.05 4.78 0.02

S6

Cumulative rainfall (103 cm3) Depression storage (103 cm3) Water on contributing area (103 cm3) Cumulative infiltration (103 cm3) Cumulative runoff (103 cm3) Mean depth of wetting front (cm) Coefficient of variation of wetting front depths

462.00 0.00 6.37 348.05 107.58 7.25 0.00

462.00 0.00 0.00 240.48 221.52 5.00 0.00

462.00 0.00 0.00 291.53 170.47 6.00 0.01

462.00 0.00 0.00 233.58 228.42 4.75 0.01

Rain 1: steady rainfall with an intensity of 2.2 cm/hr; and Rain 2–4: unsteady rainfall with the same cumulative rain as that of Rain 1.

6000

3

Volume (cm )

5000 4000 cumulative runoff maximum depression storage

3000 2000 1000 0 0.2

0.4

0.6

0.8

1.0

1.2

1.4

1.6

RR (cm) Fig. 11. Relationships between the simulated cumulative runoff or maximum depression storage and random roughness (RR) for surfaces RR1–RR6.

(Fig. 4) by applying Rain 2–4 (Fig. 5). The simulated hydrographs are shown in Fig. 12 and the simulation results are summarized in Table 4. It can be observed from Fig. 12 that the hydrographs

of different surfaces were significantly affected by the temporal variations of rainfall although they had the same cumulative rainfall by t = 1.0 h (Table 4). The peaks of Rain 2–4 ranged from 12.13 to 13.39 cm/hr (Fig. 5). These unsteady rain events (Rain 2–4) generated more surface runoff during the high rainfall intensity periods (Fig. 12) and yielded higher cumulative runoff than the steady rain event (Rain 1) (Table 4). The rainfall event with a peak occurring at the later stage (e.g., Rain 4) produced a higher peak flow rate (e.g., 32,118 cm3/min for surface S5, Fig. 12), and resulted in greater cumulative runoff and smaller cumulative infiltration than the rain events with an earlier peak (e.g., Rain 2 and 3) (Table 4). This can be attributed to the combined influence of the higher rainfall intensity and the decreased soil infiltrability at the later stage. For the unsteady rainfall (Rain 2–4), surface microtopography also had clear influences on the simulated hydrographs. The rough field plot surface S4 had a delayed rising limb in the hydrograph, but reached a peak flow rate similar to that of surfaces S5 and S6 (Fig. 12). More importantly, surface S4 had a greater final mean

S4 S5 S6 S4 S5 S6 S4 S5 S6 S6

Discharge (cm 3/min)

30000

20000

-

rain 2 rain 2 rain 2 rain 3 rain 3 rain 3 rain 4 rain 4 rain 4 rain 1

10000

0 0.0

0.2

0.4

0.6

0.8

1.0

Time (hr) Fig. 12. Simulated hydrographs for field surfaces S4–S6 under four rainfall events (Rain 1: steady rain with an intensity of 2.2 cm/hr; and Rain 2–4: unsteady rain with the same cumulative rainfall as Rain 1).

92

J. Yang, X. Chu / Advances in Water Resources 78 (2015) 80–93

depth of wetting front than surfaces S5 and S6 for Rain 2–4 (Table 4), indicating that increased infiltration and unsaturated flow occurred at the rough surface S4 due to its larger and deeper depressions and greater ponded water depths. The differences in the final mean wetting front depths between the rough and smooth field plot surfaces (i.e., surfaces S4 and S5) narrowed down as the peak of rain was delayed from Rain 3 to Rain 2 and Rain 4 (Table 4). In addition, as indicated by the coefficient of variation of wetting front depths (Table 4), the spatial variations in the wetting front depths decreased for the rain events with a later peak for surfaces S4 and S5. The cumulative infiltration was affected by surface microtopography and rainfall characteristics. Surface depressions were critical to the rainfall partition process. Not only did they function as water storage, but also served as a persistent water source that enhanced infiltration. Surface microtopography also influenced unsaturated flow and controlled the spatial variations in wetting front depths. Thus, it can be concluded that the combined influence of rainfall characteristics and surface microtopography altered the rainfall partition, outlet peak flow, cumulative runoff, surface depression storage, cumulative infiltration, and the spatial variations in wetting front. The findings from this study emphasized the importance of surface microtopography and rainfall features in the dynamic P2P overland flow, infiltration, and unsaturated flow processes.

overland flow and infiltration processes across scales. In particular, characterization of the P2P overland flow dynamics and quantification of the associated threshold behaviors would allow one to visualize the ‘‘real’’ puddle filling, spilling, merging, and splitting overland flow processes across a microtopographic surface, instead of a conceptualized ideal sheet flow on a uniform sloping surface, and to better understand the underlying overland flow generation mechanism. The P2P modeling system also will be useful to examine the spatio-temporally varying P2P dynamics, which can be further applied for characterizing the related threshold behaviors and analyzing hydrologic connectivity. Broader applications of the P2P modeling system may include wetland delineation and modeling, and integrated ecohydrologic modeling and assessment. A variety of microtopographic soil surfaces and different rainfall patterns have been selected and the real observed experimental data have been used to test the P2P model developed in this study and demonstrate its applicability. Its performance has been verified by comparing with the experimental data. However, the tests and applications of the current model are limited to laboratory- and field-scale soil surfaces. For large-scale watershed modeling, more hydrologic and hydraulic processes (e.g., open channel flow and surface water-groundwater interactions) need to be considered.

4. Summary and conclusions

This material is based upon work supported by the National Science Foundation under Grant No. EAR-0907588 and NSF EPSCoR Award IIA-1355466. The North Dakota Water Resources Research Institute also provided partial financial support in the form of a graduate fellowship for the first author. The authors would like to thank Jianli Zhang, Leif Sande, Daniel Bogart, Yaping Chi, Noah Habtezion, Yang Liu, and Yingjie Yang for their contributions to the related work.

The unique features of the developed P2P model include: (1) it featured a unique two-domain (C2C and P2P), hierarchical modeling structure, which improved the model performance in the simulation of overland flow under complex microtopographic conditions; (2) it physically simulated overland flow for contributing cells, which was internally coupled with the modeling of P2P dynamics by accounting for the threshold behaviors and interactions of puddles (i.e., filling, spilling, merging, and splitting); and (3) it simulated infiltration and unsaturated flow in a microtopography-dominated surface–subsurface system under the conditions of spatio-temporally varied rainfall, heterogeneous soil, and arbitrary initial soil moisture and surface ponded water distributions. The major data needed for real applications include: surface roughness and microtopographic data, soil distributions (across the surface area and along soil profiles) and soil hydraulic properties, as well as spatiotemporally variable meteorologic data (e.g., rainfall and evapotranspiration). Not only did the model provide essential information on overland flow, infiltration, surface ponding, and unsaturated flow, but also the details on the dynamic variability in hydrologic connectivity under various rainfall, microtopography, and soil conditions. This modeling study indicated that surface microtopography strongly influenced overland flow generation, featuring a threshold-controlled, dynamic expansion of the areas that contributed surface runoff to the outlet(s), which largely depended on the relationships of puddles and their microtopographic properties. Surface microtopography and rainfall characteristics also affected the quantity and distributions of infiltration and unsaturated flow, which resulted in different wetting front movement patterns. It was demonstrated that surface depressions played an important role in rainfall partition between surface runoff and infiltration. More importantly, they altered the spatio-temporal redistribution and interactions of infiltration and runoff water. Due to the influence of surface microtopography, wetting front movement exhibited high variability. This study highlighted the combined effects of microtopography and rainfall patterns on P2P overland flow dynamics, infiltration, and unsaturated flow. The new modeling methodology and the P2P model developed in this study can be used to analyze depression storage-runoffinfiltration relationships and investigate the complexity of

Acknowledgments

References [1] Abedini MJ, Dickinson WT, Rudra RP. On depressional storages: The effect of DEM spatial resolution. J Hydrol 2006;318:138–50. http://dx.doi.org/10.1016/ j.jhydrol.2005.06.010. [2] Allen RG, Pereira LS, Raes D, Smith M. Crop evapotranspiration, guidelines for computing crop water requirements, FAO irrigation and drainage, Paper 56. Rome: FAO – Food and Agriculture Organization of the United Nations; 1998. [3] Antoine M, Javaux M, Bielders C. Integrating subgrid connectivity properties of the micro-topography in distributed runoff models, at the interrill scale. J Hydrol 2011;403:213–23. http://dx.doi.org/10.1016/j.jhydrol.2011.03.027. [4] Antoine M, Javaux M, Bielders C. What indicators can capture runoff-relevant connectivity properties of the micro-topography at the plot scale? Adv Water Resour 2009;32(8):1297–310. http://dx.doi.org/10.1016/j.advwatres.2009.05. 006. [5] Appels WM, Bogaart PW, van der Zee S. Influence of spatial variations of microtopography and infiltration on surface runoff and field scale hydrological connectivity. Adv Water Resour 2011;34(2):303–13. http://dx.doi.org/ 10.1016/j.advwatres.2010.12.003. [6] ASTM. Standard test methods for laboratory determination of water (moisture) content of soil and rock by mass, ASTM D2216–10, ASTM International, West Conshohocken, PA; 2010. [7] Chen L, Sela S, Svoray T, Assouline S. The role of soil-surface sealing, microtopography, and vegetation patches in rainfall–runoff processes in semiarid areas. Water Resour Res 2013;49(9):5585–99. http://dx.doi.org/ 10.1002/wrcr.20360. [8] Chi Y, Yang J, Bogart D, Chu X. Fractal analysis of surface microtopography and its application in understanding hydrologic processes. Trans ASABE 2012;55(5):1781–92. http://dx.doi.org/10.13031/2013.42370. [9] Chu X, Marino MA. Determination of ponding condition and infiltration into layered soils under unsteady conditions. J Hydrol 2005;313:195–207. http:// dx.doi.org/10.1016/j.jhydrol.2005.03.002. [10] Chu X, Marino MA. Simulation of infiltration and surface runoff – A windowsbased hydrologic modeling system HYDROL-INF. In: Randall Graham, editor. Examining the confluence of environmental and water concerns, Proceedings of the 2006 world environmental and water resources congress, American Society of Civil Engineers; 2006. P. 1–8. doi: 10.1061/40856(200)429. [11] Chu X, Yang J, Chi Y. Quantification of soil random roughness and surface depression storage: Methods, applicability, and limitations. Trans ASABE 2012;55(5):1699–710. http://dx.doi.org/10.13031/2013.42361.

J. Yang, X. Chu / Advances in Water Resources 78 (2015) 80–93 [12] Chu X, Yang J, Chi Y, Zhang J. Dynamic puddle delineation and modeling of puddle-to-puddle filling–spilling–merging–splitting overland flow processes. Water Resour Res 2013;49(6):3825–9. http://dx.doi.org/10.1002/wrcr.20286. [13] Chu X, Zhang J, Chi Y, Yang J. An improved method for watershed delineation and computation of surface depression storage. In: Potter KW, Frevert DK, editors. Watershed management 2010: Innovations in watershed management under land use and climate change, Proceedings of the 2010 watershed management conference, American Society of Civil Engineers; 2010. p. 1113–22. doi: 10.1061/41143(394)100. [14] Darboux F, Davy P, Gascuel-Odoux C. Effect of depression storage capacity on overland-flow generation for rough horizontal surfaces: water transfer distance and scaling. Earth Surf Proc Land 2002;27(2):177–91. http:// dx.doi.org/10.1002/esp.312. [15] Darboux F, Davy P, Gascuel-Odoux C, Huang C. Evolution of soil surface roughness and flowpath connectivity in overland flow experiments. Catena 2001;46(2–3):125–39. http://dx.doi.org/10.1016/S0341-8162(01)00162-X. [16] Darboux F, Huang C. Does soil roughness increase or decrease water and particle transfer? SSSAJ 2005;69(3):748–56. http://dx.doi.org/10.2136/ sssaj2003.0311. [17] Dunne T, Zhang W, Aubry BF. Effects of rainfall, vegetation, and microtopography on infiltration and runoff. Water Resour Res 1991;27(9): 2271–85. http://dx.doi.org/10.1029/91WR01585. [18] Esteves M, Faucher X, Galle S, Vauclin M. Overland flow and infiltration modelling for small plots during unsteady rain: numerical results versus observed values. J Hydrol 2000;228(3–4):265–82. http://dx.doi.org/10.1016/ S0022-1694(00)00155-4. [19] Fiedler FR, Ramirez JA. A numerical method for simulating discontinuous shallow flow over an infiltrating surface. Int J Numer Methods Fluids 2000;32(2):219–39. http://dx.doi.org/10.1002/(SICI)1097-0363(20000130)32: 2<219::AID-FLD936>3.0.CO;2-J. [20] Frei S, Lischeid G, Fleckenstein JH. Effects of micro-topography on surfacesubsurface exchange and runoff generation in a virtual riparian – A modeling study. Adv Water Resour 2010;33:1388–401. http://dx.doi.org/10.1016/ j.advwatres.2010.07.006. [21] Gomez JA, Nearing MA. Runoff and sediment losses from rough and smooth soil surfaces in a laboratory experiment. Catena 2005;59:253–66. http:// dx.doi.org/10.1016/j.catena.2004.09.008. [22] Gonwa WS, Kavvas ML. A modified diffusion wave equation for flood propagation in trapezoidal channels. J Hydrol 1986;83:119–36. http:// dx.doi.org/10.1016/0022-1694(86)90187-3. [23] Gottardi G, Venutelli M. A control-volume finite element model for twodimensional overland flow. Adv Water Resour 1993;16:277–84. http:// dx.doi.org/10.1016/0309-1708(93)90019-C. [24] Govers G, Takken I, Helming K. Soil roughness and overland flow. Agronomie 2000;20(2):131–46. http://dx.doi.org/10.1051/agro:2000114. [25] Hairsine PB, Moran CJ, Rose CW. Recent developments regarding the influence of soil surface characteristics on overland flow and erosion. Aust J Soil Res 1992;30:249–64. http://dx.doi.org/10.1071/SR9920249. [26] Helming K, Römkens MJM, Prasad SN. Surface roughness related processes of runoff and soil loss: A flume study. SSSAJ 1998;62:243–50. http://dx.doi.org/ 10.2136/sssaj1998.03615995006200010031x. [27] Huang C, Bradford JM. Depressional storage for Markov–Gaussian surfaces. Water Resour Res 1990;26(9):2235–42. http://dx.doi.org/10.1029/ WR026i009p02235. [28] Huang J, Song CCS. Stability of dynamic flood routing schemes. J Hydraul Eng 1985;111(12):1497–505. http://dx.doi.org/10.1061/(ASCE)0733-9429(1985) 111:12(1497). [29] Jain MK, Singh VP. DEM-based modeling of surface runoff using diffusion wave equation. J Hydrol 2005;302:107–26. http://dx.doi.org/10.1016/j.hydrol.2004. 06.04.

93

[30] Liu Y, Yang J, Chu X. Infiltration and unsaturated flow under the influence of surface microtopography – Model simulations and experimental observations. In: Patterson CL, Struck SD, Murray DJ, editors. Showcasing the future, Proceedings of the 2013 ASCE world environmental and water resources congress, American Society of Civil Engineers. p. 468–75. [31] MacCormack RW. The effect of viscosity in hypervelocity impact cratering. AIAA Paper 69–354, Ohio, Cincinnati; 1969. [32] Martin Y, Valeo C, Tait M. Centimetre-scale digital representations of terrain and impacts on depression storage and runoff. Catena 2008;75:223–33. http:// dx.doi.org/10.1016/j.catena.2008.07.005. [33] Moore DC, Singer MJ. Crust formation effects on soil erosion processes. SSSAJ 1990;54:1117–23. http://dx.doi.org/10.2136/sssaj1990.0361599500 5400040033x. [34] Morris EM, Woolhiser DA. Unsteady one-dimensional flow over a plane: Partial equilibrium and recession hydrographs. Water Resour Res 1980;16(2):355–60. http://dx.doi.org/10.1029/WR016i002p00355. [35] Nash JE, Sutcliffe JV. River flow forecasting through conceptual models. Part IA discussion of principles. J Hydrol 1970;10(3):282–90. http://dx.doi.org/ 10.1016/0022-1694(70)90255-6. [36] O’Callaghan JF, Mark DM. The extraction of drainage networks from digital elevation data. Comput Vision Graphics Image Proc 1984;28:323–44. http:// dx.doi.org/10.1016/S0734-189X(84)80011-0. [37] Rossi MJ, Ares JO. Depression storage and infiltration effects on overland flow depth-velocity-friction at desert conditions: field plot results and model. Hydrol Earth Syst Sci 2012;16:3293–307. http://dx.doi.org/10.5194/hess-163293-2012. [38] Shaw DA, Pietroniro A, Martz LW. Topographic analysis for the prairie pothole region of Western Canada. Hydrol Process 2012;27(22):3105–14. http:// dx.doi.org/10.1002/hyp.9409. [39] Smith GD. Numerical solution of partial differential equations: Finite difference methods. Oxford: Clarendon press; 1985. [40] Tatard L, Planchon O, Wainwright G, Nord J, Favis-Mortlock D, Silvera N, et al. Measurement and modelling of highresolution flow-velocity data under simulated rainfall on a low-slope sandy soil. J Hydrol 2008;348:1–12. http:// dx.doi.org/10.1016/j.jhydrol.2007.07.016. [41] Tayfur G, Kavvas ML, Govindaraju RS, Storm DE. Application of St. Venant equations for two-dimensional overland flows over rough infiltrating surfaces. J Hydraul Eng 1993;119:51–63. http://dx.doi.org/10.1061/(ASCE)07339429(1993)119:1(51). [42] Thompson SE, Katul GG, Porporato A. Role of microtopography in rainfall– runoff partitioning: an analysis using idealized geometry. Water Resour Res 2010;46(W07520). http://dx.doi.org/10.1029/2009WR008835. [43] Wang M, Hjelmfelt AT. DEM based overland flow routing model. J Hydrol Eng 1998;3(1):1–8. http://dx.doi.org/10.1061/(ASCE)1084-0699(1998)3:1(1). [44] Woolhiser DA, Liggett JA. Unsteady 1-dimensional flow over a plane—rising hydrograph. Water Resour Res 1967;3(3):753–71. http://dx.doi.org/10.1029/ WR003i003p00753. [45] Woolhiser DA. Simulation of unsteady overland flow. In: Mahmood K, Yevjevich V, editors. Unsteady flow in open channels. Fort Collins, Colorado: Water Resources Publications; 1975. p. 485–508. Vol. II. [46] Yang J, Chu X. Quantification of the spatio-temporal variations in hydrologic connectivity of small-scale topographic surfaces under various rainfall conditions. J Hydrol 2013;505:65–77. http://dx.doi.org/10.1016/j.jhydrol. 2013.09.013. [47] Zhang W, Cundy TW. Modeling of two-dimensional overland flow. Water Resour Res 1989;25(9):2019–35. http://dx.doi.org/10.1029/ WR025i009p02019. [48] Zobeck TM, Onstad CA. Tillage and rainfall effects on random roughness: A review. Soil Till Res 1987;9:1–20. http://dx.doi.org/10.1016/0167-1987(87) 90047-X.