# Experimental and Numerical Investigations of Influence on Overland Flow and Water Infiltration by Fracture Networks in Soil

Authors: Lei Zhu; Dongjun Fan; Rong Ma; Yonggen Zhang; Yuanyuan Zha
Journal: Geofluids
 (2018)

---

## Abstract

Preferential flow is common in clay or expansive clay soils, involving water bypassing a large portion of the soil matrix. Dye tracer experiment and numerical modeling are used to simulate the surface runoff and subsurface preferential flow patterns influenced by the soil fracture network of a relatively steep hillslope system (slope angle equals to 10 degrees). The result of the experiments indicates that part of the water is infiltrated through cracks, leading to the delay of the initial runoff-yielding time and reduction of the discharge of the surface runoff. The soil water flow is mainly in the matrix when the intensity of precipitation is low. With the increasing of precipitation, soil water movement may become in the form of preferential flow through cracks. In addition, the nonuniformity of soil water infiltration and the depth of the average water infiltration increase as the precipitation intensity increases. To this end, the complete coupling model was established by using the surface-matrix-crack (SMC) model to simulate water flow within discrete fracture as well as to simulate water flow in the soil matrix based on the concept of dual permeability using the traditional Richards’ equation. In this model, the “cubic law” of fluid motion in cracks within smooth parallel plates and the two-dimensional diffusion wave approximation to Saint-Venant equations with momentum term ignored (two-dimensional shallow water equations) were used. The model divides soil water infiltration into two forms and uses the overall method to calculate the exchange of water between the crack networks and matrix regions as well as the exchange water between surface runoff and infiltration water. Results indicate that the SMC model has better performance compared with the traditional equivalent continuum model when those models are used to simulate the surface runoff movement and the soil water movement in the presence of cracks.


---

## Body

## 1. Introduction

Cracks are common in clay or expansive clay soils, which are recognized as one of the reasons to form preferential flow, involving water bypassing a large portion of the soil matrix and widely considered to be a common phenomenon in hydrology. Within those cracks, the infiltration water will reach the groundwater directly and quickly, and therefore reduce the irrigation efficiency. What is more, a variety of materials in soils, including heavy metals, radionuclides, pesticides, and contaminants, cause the environmental issues through preferential flow paths.When rainwater enters a cracked soil on a sloping surface, it flows in the network of cracks and seeps into the soil matrix and finally can influence hillside hydrology [1, 2]. As water moves in the network of cracks, there is an increase in the lateral flow rate through the soil as compared to flow through intact soils [3, 4]. Consequently, the occurrences of cracks in the soil influence the water movement in the soil. The dynamic parameters of overland flow are also influenced by the drying shrinkage cracks in the soil, leading to the reduction of the unit discharge and the number of Reynolds, and the overland flow decreased accordingly. In addition, the cracks can have a dramatic effect on the processes of surface water movement and flood dynamics by altering the partitioning of rainfall between infiltration and runoff, which is an important issue in modeling and forecasting flood events. Therefore, it is important to investigate the impact of soil crack networks on overland flow and water infiltration into soils, which will improve the understanding of physical basis of nonuniform water movement in soil.Some models have been used to simulate the combination of the surface hydrological and groundwater flow, such as the Soil Conservation Service (SCS) curve number method [5], the unit hydrograph approach [6], and transfer functions [7]. In the distributed hydrological models, overland flow is usually simulated with approaches based on the lowest neighbor routing algorithm [8] or the numerical solution of the Saint-Venant equations [9]. The Saint-Venant equations can be simplified using the kinematic wave [10] or diffuse wave approximations [11]. For example, MIKE-SHE [12, 13] is a model further developed based on the SHE model provided by Danish Hydraulic Institute in the early 1990s, which has been widely used in the resource and environmental issues related to surface water, groundwater, and their dynamic interactions [12, 14]. Singh and Bhallamudi [15] developed a mathematical model to simulate the slope runoff by combining surface and subsurface water flow. In their model, the Saint-Venant equation and simple explicit numerical method without swing (ENO) were used for surface runoff, while two-dimensional Richards’ equation and fully implicit finite difference numerical method were used for subsurface water flow. Therefore, their model can be used to investigate two-dimensional effects produced by nonuniform features of subsurface flow systems and to deal with complex subsurface conditions. Bai et al. [16] optimized the SWATMOD model by considering infiltration flow in the root layer calculated by SWAT, which was used as groundwater recharge of MODFLOW in their model. Gandolfi and Savi [17] built coupled model based on two-dimensional shallow water equation and one-dimensional Richards’ equation. Joint calculus and explicit finite difference were used to solve the coupled model, which was used to simulate dry soil initial surface conditions, with cases of changing initial moisture, as well as variation characteristics of spatial properties such as slope, roughness, and soil hydraulic properties. Panday et al. [18] established a hydrological evaluation which can be used in different scales. It coupled three-dimensional saturated-unsaturated groundwater flow and diffusion wave equation of surface flow, by taking into account a series of models of rivers, channels, and hydraulic structures. Trichakis et al. [19] developed an integrated model based on the LISFLOOD model, which was used to simulate the hydrological process at a pan-European scale.A complete description of field-scale water processes calls for an implementation of a dynamic simulation model including overland flow and subsurface flow components. However, in clay rich soils, inclusion of preferential flow as well as soil shrinkage and swelling processes is also crucial. Several earlier studies demonstrate that models lacking preferential flow processes are not well-suited for describing water flow in structured soils [20–23].Several models that simulate preferential flow have been developed in the past two decades [20, 23]. Those models are generally based on either the dual-permeability, kinematic wave, or two-stage conceptual approaches, which all assume mobile and immobile regions in soil [24]. The models differ in their description of flow and transport in the preferential flow domain.For example, dual-porosity and dual-permeability models both assume that the porous medium consists of two interacting regions, one associating with the interaggregate, macropore, or fracture system and the other one comprising micropores (or intra-aggregate pores) inside soil aggregates or the rock matrix. Dual-porosity models assume that water in the matrix is stagnant, but dual-permeability models assume that water can move in both the inter- and intra-aggregate pore regions [24, 25]. Available dual-permeability models differ mainly in how they implement water flow in and between the two pore regions, especially with respect to the degree of simplification and empiricism [20]. For the application of these models, it is necessary to know the properties of soil, and thus the uses of these models are obviously confined by irregular distribution of macropores and fractures in the soil.Dual-porosity models divide the total pore space into mobile and immobile (only water storage allowed) parts, while in dual- and multipermeability models, all the pore systems are mobile.According to Shao et al. [26], single pore system models do not emulate the true preferential flow well, while the multipermeability system approaches can be computationally demanding and difficult to be parameterized [27]. Gärdenäs et al. [28] showed dual-porosity approaches to be less realistic than dual-permeability approaches. Dual-permeability models can be further divided into conceptual [29], gravity-based, and Darcian approaches.Based on how preferential flow is simulated in the macropore system, the gravity-based approaches for preferential flow might be suitable for forest soils [30], but in agricultural fields, subsurface drains create pressure-driven flow fields which have to be handled with Darcian approaches [31].Several previous studies implicate that one-dimensional (1D) preferential flow and transport models are not comprehensive enough to simulate subsurface flow and solute transport in undulating fields and recommend the use of two-dimensional (2D) profile models or three-dimensional (3D) models instead [32–35]. However, there are a few 2D and 3D models that support simulation of preferential flow and transport. Five 3D soil and groundwater flow and substance transport models were found that simulate preferential flow and transport, including HILLFLOW [36], TOUGH2 [37], HYDRUS-3D [38], MIKE-SHE [12, 13], and a model by Laine-Kaulio et al. [39]. In addition, Blessent et al. [40] used the FRAC3Dvs numerical model to compare the dual porosity, equivalent porous medium, and the discrete fracture matrix diffusion conceptual model, and the study has an important role in decision-making for groundwater-related issues, while the research focus is the groundwater and the three models do not consider the surface runoff movement. Although soil swelling and shrinkage models have been included in several earlier models, most notably in FLOCR [41], MACRO [42], and SWAP [43], they are all limited to 1D or quasi 2D simulations.Water movement and solute transport in cracked soil have been widely investigated in the fields of geotechnical engineering, agriculture, and environmental sciences [1, 2]. The permeability of cracked soil is highly anisotropic due to the cracks. The permeability in the predominant direction of cracking is larger than that in other directions. Krisnanto et al. [44] used the model to predict the lateral flow rate through a network of cracks in the soils. Zhang et al. [45] simulated the water infiltration rate by considering the crack property. It is generally agreed that it will be difficult to quantitatively describe the anisotropic properties of cracked soil without an appropriate 3D crack network model. In the crack network model, the permeability tensor is often used to represent the anisotropic properties of a cracked medium [46, 47]. The permeability of cracked soil is often analyzed by means of continuum mechanics methods based on the permeability of the soil matrix and crack network [48], and thus the permeability tensor must be determined. Khoei et al. [49] have developed a numerical method to simulate the two-phase fluid flow through fractured porous media based on the linear momentum balance equation and the flow continuity equation. Li et al. [50], Krisnanto et al. [3], and Pierlot et al. [42] developed a method to determine whether a network of cracks can be considered as a continuum with anisotropic characteristics for the saturated coefficient of permeability.In this study, the dye-tracing experiments were conducted that simulated a series of natural rainfall amounts. In conjunction with laboratory experiments, numerical modelling was built to simulate coupled water movement in soil crack networks, surface water, and groundwater flow under the experimental conditions.


## 2. Materials and Methods

### 2.1. Field Experiment

#### 2.1.1. Experiment Sites

The experiments were carried out in Efficient Use of Water Resources in Arid Modern Agriculture, Ministry of Educational Engineering Research Center, western part of the Ningxia, Hui Autonomous Region, China, between July 2016 and August 2016. The climate in this area is temperate monsoon climate with a mean annual temperature of 9.5°C, an annual precipitation of 202.1 mm, and an average annual evaporation of 1947.1 mm. The soil physical properties of the experimental area are shown in Table1.Table 1
Soil physical properties in experiment area.
Soil depth (cm)
Particle size distribution (%)
Density (g·cm−3)
Porosity (%)
>50μm
2~50μm
<2μm
0~10
3.7
51.0
45.5
1.4253
45.2
>10~20
4.6
52.5
45.4
1.5445
45.0
>20~30
4.5
51.6
45.3
1.6876
44.8
>30~50
4.8
50.4
45.6
1.6424
44.3


#### 2.1.2. Experimental Design

According to different rainfall intensities and durations, the experiment is divided into five groups that are shown in Table2. Each experimental plot, with 100 cm long and 100 cm wide, is used to simulate a hillslope system with a relatively large slope angle of 10° (as shown in Figure 1), and the surface soil of 10 cm depth was removed to ensure that there is no large pore. Then, a piece of steel sheet with 0.5 cm thickness was squeezed into and then pull out the soil carefully in the experiment area to form the artificial rectangular crack networks which consist of 10 cracks, 30 cm in depth and 80 cm in length for each crack (as shown in Figure 1). At the beginning of the experiments, totally, 42 soil samples were collected at various depths to measure initial soil water content and the collector for rain was placed at the other side (lower boundary on the right hand side).Table 2
Experimental conditions of the plots.
Experiment label
Precipitation intensity (mm·h−1)
Duration of precipitation (min)
Existence of the preferential flow
S1
30
60
Yes
S2
50
30
Yes
S3
50
30
Yes
S4
30
60
No
S5
50
30
NoFigure 1
Sketch of experimental arrangement.The tracer of potassium iodide (KI) was used to investigate the flow patterns of infiltrated water, because the iodide ion (I−) can be easily visualized when mixed with starch. The iodide indication solution was prepared as described by Wang et al. [51], with a concentration of 20 g/L for KI.When the experiment began, the indication solution was sprayed above the plot area, resulting in a uniform irrigation system to mimic the outdoor rainfall shown in Figure1. When the amount of water in the gauge was increased to 100 mL, the time was recorded until the experiment was finished.The experiment area was then covered with plastic sheets to avoid evaporation in the following 12 hours when the indication solution penetrated the surface completely. Finally, the vertical soil profiles were excavated in each plot along thex-axis direction, with the interval spacing of 5–10 cm between two vertical profiles.For the tracer experiments, Fe3+ was used as a mild oxidizer for I−, which will not bleach the formed starch-iodide color. In addition, the effect of starch dissolved into the I− solution was better than that by applying starch powder directly to the soil profile, which could avoid the wetting process of the starch layer by I− solution in soil. The indication solution consists of 50 g/L of starch and 0.05 g/L iron(III) nitrate (Fe(NO3)3). The indication solution was added with a pressurized sprayer that can release a continuous and fine mist of indication solution. The indication solution was sprayed onto the surface as uniformly as possible by long, slow strokes. After spraying, it took about 10 min for a dark blue-violet color to develop as a result of the oxidizer iron(III) nitrate.Staining patterns of the soil profiles were taken using a CCD digital camera in daylight. Nine profiles were prepared and taken within 1 day, and procedure was taken following Klaus et al. [52]. A white sheet was placed over each plot to diffuse the light and to avoid direct radiation. The soil profiles were photographed with a Kodak color and gray scale. To adjust the color differences affected by the environment, photographs of a gray board were taken and used as a background.After taking photographs, soil samples were collected within the stained and unstained regions to measure soil water content distributions. More than fifty soil samples were collected in the stained and unstained regions, respectively, with each depth of 0, 0.10, 0.20, 0.30, 0.40, 0.50, 0.60, 0.70, and 0.80 m.


### 2.2. Image Processing

The most important sources of variability in soil images are usually with uneven illumination, geometric distortion, and the condition of the surface [53, 54]. These problems were minimized by fixed environmental conditions and carefully prepared soil face by hand prior to the acquisition of each image. Grey scales were positioned at the sides of the image for corrections of white balances. High-resolution (2048 × 1229 pixels) digital images were collected in the BMP file format. Corrections were made for white balance between images in order to standardize the images for comparison.In order to classify images into stained or unstained areas for the soil [55], a threshold algorithm was used to transform color images into black and white (binary) images shown in Figure 2. Specifically, red, green, and blue values in any pixels that fell within the specified range were classified into black or white accordingly, and the files are saved into binary format. The visual analysis of flow patterns was conducted qualitatively based on the generated binary images. Dye-stained areas were defined as the preferential flow region, whereas nonstained material represents the bypassed soil matrix.Figure 2
Heterogeneous dye coverage patterns at experiment plot before (a) and after (b) image processing.
(a)
(b)In order to quantify the degree of preferential flow, the binary data was used to calculate dye coverage (DC) which indicates the percentage ratio of the dye-stained area to the total profile area [51]. It should be noted that low values of the parameter indicate a low degree of preferential flow.


### 2.3. Model Description

The SMC (surface-matrix-crack) model can simulate water flow within discrete fracture as well as water flow in the soil matrix based on the dual-permeability concept using the equivalent freshwater formulation.The porous low-permeability matrix is represented by regular 3D blocks, while fractures of high permeability are represented by 2D rectangular planes. Fluid flow and transport equations are solved for both porous matrix and fractures. Vertical, horizontal, and inclined 2D fractures are incorporated into the 3D grid by superimposing 2D fracture elements into the 3D regular block elements. Faces and blocks share common nodes along the fracture walls. Thus, nodes at fracture locations receive fluid and mass contributions from both the block elements as well as from the fracture faces. Furthermore, for these mutual nodes, both hydraulic head and concentration at the fracture/matrix interface are assumed to be equal, ensuring continuity of these variables between fractures and matrix [56, 57].Mathematical modeling of overland flow involves the solution of the governing equations for both the surface flow and the groundwater flow by the interior physical mechanism. The interaction among hydrogeological condition, surface water, and groundwater should be considered. The surface flow equation is solved on a 2D finite element mesh stacked upon a subsurface grid when solving for both domains. For superposition, the grid generated for the subsurface domain is mirrored areally for the surface flow nodes, with surface flow node elevations corresponding to the top elevation of the topmost active layer of the subsurface grid.


#### 2.3.1. Soil Water Infiltration Movement

The water movement in the matrix area can be described using the Richards’ equation [20, 58], and the cubic law can be used to describe the preferential flow [59]:
(1)∂∂tθsSw=∂∂xiKskrw∂ϕ+Z∂xj+Q+QGf+QGs,and
(2)∂∂xi2bKfkrwf∂ϕf+Zf∂xj±QGf=2b∂Swf∂t,where θs is the saturated moisture content, L3L−3; Sw is the soil saturation; Swf is the water channel saturation, dimensionless; 2b is the flow channel width with b indicating the radius of the channel, L; ϕ and ϕf represent the soil water matrix potential and the flow channel water potential, L; Z is the elevation head and Zf is the elevation head within the fracture, L; Q is the root water uptake term or other sources and sinks, L3L−3T−1; xii=1,2,3 is spatial coordinates, x1=x,y denotes horizontal coordinate, and x3=z is the vertical coordinate and is positive upward; t is time, T; Ks and Kf represent the saturated hydraulic conductivity of soil and the water channel, LT−1; krw and krwf represent the relative hydraulic conductivity of soil and the water channel, dimensionless; QGf is the fracture water and soil water exchange capacity, L3L−3T−1; and QGs is the soil water and surface water exchange capacity, L3L−3T−1. The characteristic curve of soil water content and unsaturated hydraulic conductivity are represented based on the Mualem-van Genuchten model [60, 61].


#### 2.3.2. Surface Flow Model

Two-dimensional Saint-Venant equations with momentum term ignored (two-dimensional shallow water equations) were used to describe surface runoff movement [58]:
(3)∂hs∂t−∂∂xdsksx∂hs∂x−∂∂ydsksy∂hs∂y±dsQGs+Qs=0,where hs is the surface altitude, L; hs=Zs+ds; Zs is the ground altitude, L; ds is the water depth, L; Qs is the source-sink item, L3L−3T−1; ksx and ksy are land surface conduction coefficient in x and y directions. For surface runoff model parameters, ksx and ksy can be calculated, respectively, by the two formulas as follows:
(4)ksx=ds2/3nsx1∂hs/∂s1/2,ksy=ds2/3nsy1∂hs/∂s1/2,where nsx and nsy are Manning roughness coefficient in x direction and y direction, L−1/3T.


#### 2.3.3. Fracture Network Coupling with the Matrix Region

The discrete equation of soil water movement in the soil is described as follows:(5)θsSwL+1−SwLΔtI⋅VI=∑J∈ηIγIJλIJL+1⋅τIJ+QIL+1±QGfI±QGSI⋅VI+∫Beq∗⋅NI⋅dBe,where VI is a factor for impact volume node I, L; γIJ is the water potential difference between the nodes I and J matrix area, L; λIJ represents a relative unsaturated hydraulic conductivity between the node I and node J when using upstream weighting; ηI represents the order junction node I when the region is discrete; τIJ is the coefficient matrix; q∗ refers to the flow of flux on the border near the border node I; outflow is positive, for the internal node I, q∗=0; Be represents the study area boundary, L; and NI is the corresponding basis function in a unit.The discrete equation of preferential water movement within soil is described as follows:(6)2b⋅SwL+1−SwLΔt⋅AfI=2b⋅∑J∈ηfIγfIJλfIJL+1⋅τfIJ±QGfIL+1⋅AfI,where AfI is an impact factor for node fI, L2; ηfI represents the node order of junction connect I when region is discrete; γfIJ is head difference between the node I and node J, L; λfIJ indicates the relative conductivity between the node I to node J while using the weighted upwind method, dimensionless; and τfIJ is the coefficient matrix. QGfIL+1⋅AfI is the exchange water of the two regions.Replacing the exact value of the fluid exchangeQGfIL+1⋅AfI into (5) gives the following two-domain coupling equation:
(7)θsSwL+1−SwLΔtI⋅VI=∑J∈ηIγIJλIJL+1⋅τIJ+QIL+1±QGSI⋅VI+∫Beq∗⋅NI⋅dBe+2b⋅∑J∈ηfIγfIJλfIJL+1⋅τfIJ−2b⋅SwL+1−SwLΔt⋅AfI.


#### 2.3.4. Surface Coupling with the Subsurface

The discretized surface (7) is coupled with the 3D subsurface flow (5) via superposition. Fully implicit coupling of the surface and subsurface flow regimes provides an integral view of the movement of water, as opposed to the traditional division of surface and subsurface regimes. Flux across the land surface is, therefore, a natural internal process allowing water to move between the surface and subsurface flow systems as governed by local flow hydrodynamics, instead of using physically artificial boundary conditions at the interface. This provides significant robustness and accuracy over flux linkage techniques, often used when conjunctive modeling is required for the surface and subsurface regimes.From soil water flow movement discrete (5), the node I has double character of the two zones if the node I is shared by soil and surface.
(8)hSIL+1−hSILASIΔt=∑J∈ηSIλsSIJ+1/2L+1γSIJL+1⋅τSIJ+QSIL+1+QGfIL+1∓dsQGSIL+1⋅ASI.PresumehI=hSI, for the soil water flow movement discrete equation of the node I, QSIL+1+QGfIL+1∓dsQGSIL+1⋅ASI contains water exchange item of the two zones. But for the discrete equation of surface, dsQGSIL+1⋅ASI is the exchange item of the two zones; (8) can be shown as follows:
(9)dsQGSIL+1⋅ASI=∓hSIL+1−hSILASIΔt−∑J∈ηSIλsSIJ+1/2L+1γSIJL+1⋅τSIJ−QSIL+1+QGfIL+1⋅ASI.Substituting (9) into (8) leads to the integrated equation of two domains:
(10)θsSwL+1−SwLΔtI⋅VI=∑J∈ηIγIJλIJL+1⋅τIJ+QIL+1⋅VI+∫Beq∗NIdBe∓hSIL+1−hSILASIΔt−∑J∈ηSIλsSIJ+1/2L+1γSIJL+1⋅τSIJ−QSIL+1+QGfIL+1⋅ASI.


### 2.4. Initial and Boundary Conditions

The simulation area is a hexahedron consisting of parallel quadrilaterals with an angle of 80°, and the length of the parallelogram and the length perpendicular to the edge are both 1.0 m. The left and right boundaries were assumed to be the no-flow boundaries, and the top and bottom boundary was assumed to be the atmospheric boundary and critical depth boundary, respectively. Precipitation intensity and duration are shown in Table2, with distribution of the precipitation evenly. The discharge (Qs) per unit width at the critical depth boundary is given by the following:
(11)Qs=gds3.Surrounded boundaries were set as impermeable boundaries. Experimental data of soil water content was used as initial conditions for simulations.


### 2.5. Model Evaluation Method

In order to assess model performance, coefficient of determinationR2 and root mean square error (RMSE) are used to evaluate soil moisture prediction, while relative error (RE) is used to assess the surface runoff prediction performance. The mathematical expressions for R2, RMSE, and RE are
(12)SEobs=∑n=1NSwobs,i−Swobs¯2,SE=∑n=1NSwobs,i−Swsim,i2,R2=1−SESEobs,RMSE=SEn,RE=WSi−WmWSr×100%.In the (12), Swsim,i is the simulated value which is obtained through forward modeling by measured parameters or optimization parameters; Swobs,i is the observed value; Swobs¯ and Swsim¯ are the average values of Swobs,i and Swsim,i, respectively; N is the number of sample which is 9 in this paper; SEobs indicates the sum of squares of the total observations; SE is the sum of squares for error. WSi is the simulated value of the discharge of the surface runoff; Wm is the measured value of the discharge of the surface runoff; and WSr is the precipitation of the experimental plot.


### 2.6. Calibration and Parametrization of Models

Experiments S1 and S4 were chosen for the model calibration because of the completeness of the dataset. Specifically, experiment S1 was used to calibrate the saturation, while experiment S4 was to calibrate the surface runoff parameters. The calibration was conducted manually by adjusting the model parameters (one at a time) to reduce the residual between the predictive and the measured values. The model was forced by prescribing the observed rainfall. Model evaluation parameters which were calibrated against soil saturation and discharge indicate that the calibrated parameters could be reliable enough for model calibration, as shown in Figure3 and Table 3.Figure 3
Comparison of the measured and SMC-simulated value of experiment S1.
(a)
Comparison of soil saturation
(b)
Comparison of dischargeTable 3
Model evaluation parameters of the experiment S1.
Model calibration
R2
RMSE
RE
Value
0.87
0.039
0.2%The calibrated model parameter values (as shown in Table4) achieved the best which indicate satisfactory performances of the model as suggested by Zhu et al. [23].Table 4
Hydrodynamic parameters of the soil matrix and surface runoff parameters.
Bulk density
Porosity
Saturated hydraulic conductivity
Residual moisture
Saturated moisture
Shape parameter
Reciprocal of inhale
Manning coefficient
g·cm−3
%
Ks (m·h−1)
θr
θs
n
α (m−1)
min·cm1/3
1.58
46.0
0.005
0.034
0.45
1.37
0.01
0.0000001


## 2.1. Field Experiment

### 2.1.1. Experiment Sites

The experiments were carried out in Efficient Use of Water Resources in Arid Modern Agriculture, Ministry of Educational Engineering Research Center, western part of the Ningxia, Hui Autonomous Region, China, between July 2016 and August 2016. The climate in this area is temperate monsoon climate with a mean annual temperature of 9.5°C, an annual precipitation of 202.1 mm, and an average annual evaporation of 1947.1 mm. The soil physical properties of the experimental area are shown in Table1.Table 1
Soil physical properties in experiment area.
Soil depth (cm)
Particle size distribution (%)
Density (g·cm−3)
Porosity (%)
>50μm
2~50μm
<2μm
0~10
3.7
51.0
45.5
1.4253
45.2
>10~20
4.6
52.5
45.4
1.5445
45.0
>20~30
4.5
51.6
45.3
1.6876
44.8
>30~50
4.8
50.4
45.6
1.6424
44.3


### 2.1.2. Experimental Design

According to different rainfall intensities and durations, the experiment is divided into five groups that are shown in Table2. Each experimental plot, with 100 cm long and 100 cm wide, is used to simulate a hillslope system with a relatively large slope angle of 10° (as shown in Figure 1), and the surface soil of 10 cm depth was removed to ensure that there is no large pore. Then, a piece of steel sheet with 0.5 cm thickness was squeezed into and then pull out the soil carefully in the experiment area to form the artificial rectangular crack networks which consist of 10 cracks, 30 cm in depth and 80 cm in length for each crack (as shown in Figure 1). At the beginning of the experiments, totally, 42 soil samples were collected at various depths to measure initial soil water content and the collector for rain was placed at the other side (lower boundary on the right hand side).Table 2
Experimental conditions of the plots.
Experiment label
Precipitation intensity (mm·h−1)
Duration of precipitation (min)
Existence of the preferential flow
S1
30
60
Yes
S2
50
30
Yes
S3
50
30
Yes
S4
30
60
No
S5
50
30
NoFigure 1
Sketch of experimental arrangement.The tracer of potassium iodide (KI) was used to investigate the flow patterns of infiltrated water, because the iodide ion (I−) can be easily visualized when mixed with starch. The iodide indication solution was prepared as described by Wang et al. [51], with a concentration of 20 g/L for KI.When the experiment began, the indication solution was sprayed above the plot area, resulting in a uniform irrigation system to mimic the outdoor rainfall shown in Figure1. When the amount of water in the gauge was increased to 100 mL, the time was recorded until the experiment was finished.The experiment area was then covered with plastic sheets to avoid evaporation in the following 12 hours when the indication solution penetrated the surface completely. Finally, the vertical soil profiles were excavated in each plot along thex-axis direction, with the interval spacing of 5–10 cm between two vertical profiles.For the tracer experiments, Fe3+ was used as a mild oxidizer for I−, which will not bleach the formed starch-iodide color. In addition, the effect of starch dissolved into the I− solution was better than that by applying starch powder directly to the soil profile, which could avoid the wetting process of the starch layer by I− solution in soil. The indication solution consists of 50 g/L of starch and 0.05 g/L iron(III) nitrate (Fe(NO3)3). The indication solution was added with a pressurized sprayer that can release a continuous and fine mist of indication solution. The indication solution was sprayed onto the surface as uniformly as possible by long, slow strokes. After spraying, it took about 10 min for a dark blue-violet color to develop as a result of the oxidizer iron(III) nitrate.Staining patterns of the soil profiles were taken using a CCD digital camera in daylight. Nine profiles were prepared and taken within 1 day, and procedure was taken following Klaus et al. [52]. A white sheet was placed over each plot to diffuse the light and to avoid direct radiation. The soil profiles were photographed with a Kodak color and gray scale. To adjust the color differences affected by the environment, photographs of a gray board were taken and used as a background.After taking photographs, soil samples were collected within the stained and unstained regions to measure soil water content distributions. More than fifty soil samples were collected in the stained and unstained regions, respectively, with each depth of 0, 0.10, 0.20, 0.30, 0.40, 0.50, 0.60, 0.70, and 0.80 m.


## 2.1.1. Experiment Sites

The experiments were carried out in Efficient Use of Water Resources in Arid Modern Agriculture, Ministry of Educational Engineering Research Center, western part of the Ningxia, Hui Autonomous Region, China, between July 2016 and August 2016. The climate in this area is temperate monsoon climate with a mean annual temperature of 9.5°C, an annual precipitation of 202.1 mm, and an average annual evaporation of 1947.1 mm. The soil physical properties of the experimental area are shown in Table1.Table 1
Soil physical properties in experiment area.
Soil depth (cm)
Particle size distribution (%)
Density (g·cm−3)
Porosity (%)
>50μm
2~50μm
<2μm
0~10
3.7
51.0
45.5
1.4253
45.2
>10~20
4.6
52.5
45.4
1.5445
45.0
>20~30
4.5
51.6
45.3
1.6876
44.8
>30~50
4.8
50.4
45.6
1.6424
44.3


## 2.1.2. Experimental Design

According to different rainfall intensities and durations, the experiment is divided into five groups that are shown in Table2. Each experimental plot, with 100 cm long and 100 cm wide, is used to simulate a hillslope system with a relatively large slope angle of 10° (as shown in Figure 1), and the surface soil of 10 cm depth was removed to ensure that there is no large pore. Then, a piece of steel sheet with 0.5 cm thickness was squeezed into and then pull out the soil carefully in the experiment area to form the artificial rectangular crack networks which consist of 10 cracks, 30 cm in depth and 80 cm in length for each crack (as shown in Figure 1). At the beginning of the experiments, totally, 42 soil samples were collected at various depths to measure initial soil water content and the collector for rain was placed at the other side (lower boundary on the right hand side).Table 2
Experimental conditions of the plots.
Experiment label
Precipitation intensity (mm·h−1)
Duration of precipitation (min)
Existence of the preferential flow
S1
30
60
Yes
S2
50
30
Yes
S3
50
30
Yes
S4
30
60
No
S5
50
30
NoFigure 1
Sketch of experimental arrangement.The tracer of potassium iodide (KI) was used to investigate the flow patterns of infiltrated water, because the iodide ion (I−) can be easily visualized when mixed with starch. The iodide indication solution was prepared as described by Wang et al. [51], with a concentration of 20 g/L for KI.When the experiment began, the indication solution was sprayed above the plot area, resulting in a uniform irrigation system to mimic the outdoor rainfall shown in Figure1. When the amount of water in the gauge was increased to 100 mL, the time was recorded until the experiment was finished.The experiment area was then covered with plastic sheets to avoid evaporation in the following 12 hours when the indication solution penetrated the surface completely. Finally, the vertical soil profiles were excavated in each plot along thex-axis direction, with the interval spacing of 5–10 cm between two vertical profiles.For the tracer experiments, Fe3+ was used as a mild oxidizer for I−, which will not bleach the formed starch-iodide color. In addition, the effect of starch dissolved into the I− solution was better than that by applying starch powder directly to the soil profile, which could avoid the wetting process of the starch layer by I− solution in soil. The indication solution consists of 50 g/L of starch and 0.05 g/L iron(III) nitrate (Fe(NO3)3). The indication solution was added with a pressurized sprayer that can release a continuous and fine mist of indication solution. The indication solution was sprayed onto the surface as uniformly as possible by long, slow strokes. After spraying, it took about 10 min for a dark blue-violet color to develop as a result of the oxidizer iron(III) nitrate.Staining patterns of the soil profiles were taken using a CCD digital camera in daylight. Nine profiles were prepared and taken within 1 day, and procedure was taken following Klaus et al. [52]. A white sheet was placed over each plot to diffuse the light and to avoid direct radiation. The soil profiles were photographed with a Kodak color and gray scale. To adjust the color differences affected by the environment, photographs of a gray board were taken and used as a background.After taking photographs, soil samples were collected within the stained and unstained regions to measure soil water content distributions. More than fifty soil samples were collected in the stained and unstained regions, respectively, with each depth of 0, 0.10, 0.20, 0.30, 0.40, 0.50, 0.60, 0.70, and 0.80 m.


## 2.2. Image Processing

The most important sources of variability in soil images are usually with uneven illumination, geometric distortion, and the condition of the surface [53, 54]. These problems were minimized by fixed environmental conditions and carefully prepared soil face by hand prior to the acquisition of each image. Grey scales were positioned at the sides of the image for corrections of white balances. High-resolution (2048 × 1229 pixels) digital images were collected in the BMP file format. Corrections were made for white balance between images in order to standardize the images for comparison.In order to classify images into stained or unstained areas for the soil [55], a threshold algorithm was used to transform color images into black and white (binary) images shown in Figure 2. Specifically, red, green, and blue values in any pixels that fell within the specified range were classified into black or white accordingly, and the files are saved into binary format. The visual analysis of flow patterns was conducted qualitatively based on the generated binary images. Dye-stained areas were defined as the preferential flow region, whereas nonstained material represents the bypassed soil matrix.Figure 2
Heterogeneous dye coverage patterns at experiment plot before (a) and after (b) image processing.
(a)
(b)In order to quantify the degree of preferential flow, the binary data was used to calculate dye coverage (DC) which indicates the percentage ratio of the dye-stained area to the total profile area [51]. It should be noted that low values of the parameter indicate a low degree of preferential flow.


## 2.3. Model Description

The SMC (surface-matrix-crack) model can simulate water flow within discrete fracture as well as water flow in the soil matrix based on the dual-permeability concept using the equivalent freshwater formulation.The porous low-permeability matrix is represented by regular 3D blocks, while fractures of high permeability are represented by 2D rectangular planes. Fluid flow and transport equations are solved for both porous matrix and fractures. Vertical, horizontal, and inclined 2D fractures are incorporated into the 3D grid by superimposing 2D fracture elements into the 3D regular block elements. Faces and blocks share common nodes along the fracture walls. Thus, nodes at fracture locations receive fluid and mass contributions from both the block elements as well as from the fracture faces. Furthermore, for these mutual nodes, both hydraulic head and concentration at the fracture/matrix interface are assumed to be equal, ensuring continuity of these variables between fractures and matrix [56, 57].Mathematical modeling of overland flow involves the solution of the governing equations for both the surface flow and the groundwater flow by the interior physical mechanism. The interaction among hydrogeological condition, surface water, and groundwater should be considered. The surface flow equation is solved on a 2D finite element mesh stacked upon a subsurface grid when solving for both domains. For superposition, the grid generated for the subsurface domain is mirrored areally for the surface flow nodes, with surface flow node elevations corresponding to the top elevation of the topmost active layer of the subsurface grid.


### 2.3.1. Soil Water Infiltration Movement

The water movement in the matrix area can be described using the Richards’ equation [20, 58], and the cubic law can be used to describe the preferential flow [59]:
(1)∂∂tθsSw=∂∂xiKskrw∂ϕ+Z∂xj+Q+QGf+QGs,and
(2)∂∂xi2bKfkrwf∂ϕf+Zf∂xj±QGf=2b∂Swf∂t,where θs is the saturated moisture content, L3L−3; Sw is the soil saturation; Swf is the water channel saturation, dimensionless; 2b is the flow channel width with b indicating the radius of the channel, L; ϕ and ϕf represent the soil water matrix potential and the flow channel water potential, L; Z is the elevation head and Zf is the elevation head within the fracture, L; Q is the root water uptake term or other sources and sinks, L3L−3T−1; xii=1,2,3 is spatial coordinates, x1=x,y denotes horizontal coordinate, and x3=z is the vertical coordinate and is positive upward; t is time, T; Ks and Kf represent the saturated hydraulic conductivity of soil and the water channel, LT−1; krw and krwf represent the relative hydraulic conductivity of soil and the water channel, dimensionless; QGf is the fracture water and soil water exchange capacity, L3L−3T−1; and QGs is the soil water and surface water exchange capacity, L3L−3T−1. The characteristic curve of soil water content and unsaturated hydraulic conductivity are represented based on the Mualem-van Genuchten model [60, 61].


### 2.3.2. Surface Flow Model

Two-dimensional Saint-Venant equations with momentum term ignored (two-dimensional shallow water equations) were used to describe surface runoff movement [58]:
(3)∂hs∂t−∂∂xdsksx∂hs∂x−∂∂ydsksy∂hs∂y±dsQGs+Qs=0,where hs is the surface altitude, L; hs=Zs+ds; Zs is the ground altitude, L; ds is the water depth, L; Qs is the source-sink item, L3L−3T−1; ksx and ksy are land surface conduction coefficient in x and y directions. For surface runoff model parameters, ksx and ksy can be calculated, respectively, by the two formulas as follows:
(4)ksx=ds2/3nsx1∂hs/∂s1/2,ksy=ds2/3nsy1∂hs/∂s1/2,where nsx and nsy are Manning roughness coefficient in x direction and y direction, L−1/3T.


### 2.3.3. Fracture Network Coupling with the Matrix Region

The discrete equation of soil water movement in the soil is described as follows:(5)θsSwL+1−SwLΔtI⋅VI=∑J∈ηIγIJλIJL+1⋅τIJ+QIL+1±QGfI±QGSI⋅VI+∫Beq∗⋅NI⋅dBe,where VI is a factor for impact volume node I, L; γIJ is the water potential difference between the nodes I and J matrix area, L; λIJ represents a relative unsaturated hydraulic conductivity between the node I and node J when using upstream weighting; ηI represents the order junction node I when the region is discrete; τIJ is the coefficient matrix; q∗ refers to the flow of flux on the border near the border node I; outflow is positive, for the internal node I, q∗=0; Be represents the study area boundary, L; and NI is the corresponding basis function in a unit.The discrete equation of preferential water movement within soil is described as follows:(6)2b⋅SwL+1−SwLΔt⋅AfI=2b⋅∑J∈ηfIγfIJλfIJL+1⋅τfIJ±QGfIL+1⋅AfI,where AfI is an impact factor for node fI, L2; ηfI represents the node order of junction connect I when region is discrete; γfIJ is head difference between the node I and node J, L; λfIJ indicates the relative conductivity between the node I to node J while using the weighted upwind method, dimensionless; and τfIJ is the coefficient matrix. QGfIL+1⋅AfI is the exchange water of the two regions.Replacing the exact value of the fluid exchangeQGfIL+1⋅AfI into (5) gives the following two-domain coupling equation:
(7)θsSwL+1−SwLΔtI⋅VI=∑J∈ηIγIJλIJL+1⋅τIJ+QIL+1±QGSI⋅VI+∫Beq∗⋅NI⋅dBe+2b⋅∑J∈ηfIγfIJλfIJL+1⋅τfIJ−2b⋅SwL+1−SwLΔt⋅AfI.


### 2.3.4. Surface Coupling with the Subsurface

The discretized surface (7) is coupled with the 3D subsurface flow (5) via superposition. Fully implicit coupling of the surface and subsurface flow regimes provides an integral view of the movement of water, as opposed to the traditional division of surface and subsurface regimes. Flux across the land surface is, therefore, a natural internal process allowing water to move between the surface and subsurface flow systems as governed by local flow hydrodynamics, instead of using physically artificial boundary conditions at the interface. This provides significant robustness and accuracy over flux linkage techniques, often used when conjunctive modeling is required for the surface and subsurface regimes.From soil water flow movement discrete (5), the node I has double character of the two zones if the node I is shared by soil and surface.
(8)hSIL+1−hSILASIΔt=∑J∈ηSIλsSIJ+1/2L+1γSIJL+1⋅τSIJ+QSIL+1+QGfIL+1∓dsQGSIL+1⋅ASI.PresumehI=hSI, for the soil water flow movement discrete equation of the node I, QSIL+1+QGfIL+1∓dsQGSIL+1⋅ASI contains water exchange item of the two zones. But for the discrete equation of surface, dsQGSIL+1⋅ASI is the exchange item of the two zones; (8) can be shown as follows:
(9)dsQGSIL+1⋅ASI=∓hSIL+1−hSILASIΔt−∑J∈ηSIλsSIJ+1/2L+1γSIJL+1⋅τSIJ−QSIL+1+QGfIL+1⋅ASI.Substituting (9) into (8) leads to the integrated equation of two domains:
(10)θsSwL+1−SwLΔtI⋅VI=∑J∈ηIγIJλIJL+1⋅τIJ+QIL+1⋅VI+∫Beq∗NIdBe∓hSIL+1−hSILASIΔt−∑J∈ηSIλsSIJ+1/2L+1γSIJL+1⋅τSIJ−QSIL+1+QGfIL+1⋅ASI.


## 2.3.1. Soil Water Infiltration Movement

The water movement in the matrix area can be described using the Richards’ equation [20, 58], and the cubic law can be used to describe the preferential flow [59]:
(1)∂∂tθsSw=∂∂xiKskrw∂ϕ+Z∂xj+Q+QGf+QGs,and
(2)∂∂xi2bKfkrwf∂ϕf+Zf∂xj±QGf=2b∂Swf∂t,where θs is the saturated moisture content, L3L−3; Sw is the soil saturation; Swf is the water channel saturation, dimensionless; 2b is the flow channel width with b indicating the radius of the channel, L; ϕ and ϕf represent the soil water matrix potential and the flow channel water potential, L; Z is the elevation head and Zf is the elevation head within the fracture, L; Q is the root water uptake term or other sources and sinks, L3L−3T−1; xii=1,2,3 is spatial coordinates, x1=x,y denotes horizontal coordinate, and x3=z is the vertical coordinate and is positive upward; t is time, T; Ks and Kf represent the saturated hydraulic conductivity of soil and the water channel, LT−1; krw and krwf represent the relative hydraulic conductivity of soil and the water channel, dimensionless; QGf is the fracture water and soil water exchange capacity, L3L−3T−1; and QGs is the soil water and surface water exchange capacity, L3L−3T−1. The characteristic curve of soil water content and unsaturated hydraulic conductivity are represented based on the Mualem-van Genuchten model [60, 61].


## 2.3.2. Surface Flow Model

Two-dimensional Saint-Venant equations with momentum term ignored (two-dimensional shallow water equations) were used to describe surface runoff movement [58]:
(3)∂hs∂t−∂∂xdsksx∂hs∂x−∂∂ydsksy∂hs∂y±dsQGs+Qs=0,where hs is the surface altitude, L; hs=Zs+ds; Zs is the ground altitude, L; ds is the water depth, L; Qs is the source-sink item, L3L−3T−1; ksx and ksy are land surface conduction coefficient in x and y directions. For surface runoff model parameters, ksx and ksy can be calculated, respectively, by the two formulas as follows:
(4)ksx=ds2/3nsx1∂hs/∂s1/2,ksy=ds2/3nsy1∂hs/∂s1/2,where nsx and nsy are Manning roughness coefficient in x direction and y direction, L−1/3T.


## 2.3.3. Fracture Network Coupling with the Matrix Region

The discrete equation of soil water movement in the soil is described as follows:(5)θsSwL+1−SwLΔtI⋅VI=∑J∈ηIγIJλIJL+1⋅τIJ+QIL+1±QGfI±QGSI⋅VI+∫Beq∗⋅NI⋅dBe,where VI is a factor for impact volume node I, L; γIJ is the water potential difference between the nodes I and J matrix area, L; λIJ represents a relative unsaturated hydraulic conductivity between the node I and node J when using upstream weighting; ηI represents the order junction node I when the region is discrete; τIJ is the coefficient matrix; q∗ refers to the flow of flux on the border near the border node I; outflow is positive, for the internal node I, q∗=0; Be represents the study area boundary, L; and NI is the corresponding basis function in a unit.The discrete equation of preferential water movement within soil is described as follows:(6)2b⋅SwL+1−SwLΔt⋅AfI=2b⋅∑J∈ηfIγfIJλfIJL+1⋅τfIJ±QGfIL+1⋅AfI,where AfI is an impact factor for node fI, L2; ηfI represents the node order of junction connect I when region is discrete; γfIJ is head difference between the node I and node J, L; λfIJ indicates the relative conductivity between the node I to node J while using the weighted upwind method, dimensionless; and τfIJ is the coefficient matrix. QGfIL+1⋅AfI is the exchange water of the two regions.Replacing the exact value of the fluid exchangeQGfIL+1⋅AfI into (5) gives the following two-domain coupling equation:
(7)θsSwL+1−SwLΔtI⋅VI=∑J∈ηIγIJλIJL+1⋅τIJ+QIL+1±QGSI⋅VI+∫Beq∗⋅NI⋅dBe+2b⋅∑J∈ηfIγfIJλfIJL+1⋅τfIJ−2b⋅SwL+1−SwLΔt⋅AfI.


## 2.3.4. Surface Coupling with the Subsurface

The discretized surface (7) is coupled with the 3D subsurface flow (5) via superposition. Fully implicit coupling of the surface and subsurface flow regimes provides an integral view of the movement of water, as opposed to the traditional division of surface and subsurface regimes. Flux across the land surface is, therefore, a natural internal process allowing water to move between the surface and subsurface flow systems as governed by local flow hydrodynamics, instead of using physically artificial boundary conditions at the interface. This provides significant robustness and accuracy over flux linkage techniques, often used when conjunctive modeling is required for the surface and subsurface regimes.From soil water flow movement discrete (5), the node I has double character of the two zones if the node I is shared by soil and surface.
(8)hSIL+1−hSILASIΔt=∑J∈ηSIλsSIJ+1/2L+1γSIJL+1⋅τSIJ+QSIL+1+QGfIL+1∓dsQGSIL+1⋅ASI.PresumehI=hSI, for the soil water flow movement discrete equation of the node I, QSIL+1+QGfIL+1∓dsQGSIL+1⋅ASI contains water exchange item of the two zones. But for the discrete equation of surface, dsQGSIL+1⋅ASI is the exchange item of the two zones; (8) can be shown as follows:
(9)dsQGSIL+1⋅ASI=∓hSIL+1−hSILASIΔt−∑J∈ηSIλsSIJ+1/2L+1γSIJL+1⋅τSIJ−QSIL+1+QGfIL+1⋅ASI.Substituting (9) into (8) leads to the integrated equation of two domains:
(10)θsSwL+1−SwLΔtI⋅VI=∑J∈ηIγIJλIJL+1⋅τIJ+QIL+1⋅VI+∫Beq∗NIdBe∓hSIL+1−hSILASIΔt−∑J∈ηSIλsSIJ+1/2L+1γSIJL+1⋅τSIJ−QSIL+1+QGfIL+1⋅ASI.


## 2.4. Initial and Boundary Conditions

The simulation area is a hexahedron consisting of parallel quadrilaterals with an angle of 80°, and the length of the parallelogram and the length perpendicular to the edge are both 1.0 m. The left and right boundaries were assumed to be the no-flow boundaries, and the top and bottom boundary was assumed to be the atmospheric boundary and critical depth boundary, respectively. Precipitation intensity and duration are shown in Table2, with distribution of the precipitation evenly. The discharge (Qs) per unit width at the critical depth boundary is given by the following:
(11)Qs=gds3.Surrounded boundaries were set as impermeable boundaries. Experimental data of soil water content was used as initial conditions for simulations.


## 2.5. Model Evaluation Method

In order to assess model performance, coefficient of determinationR2 and root mean square error (RMSE) are used to evaluate soil moisture prediction, while relative error (RE) is used to assess the surface runoff prediction performance. The mathematical expressions for R2, RMSE, and RE are
(12)SEobs=∑n=1NSwobs,i−Swobs¯2,SE=∑n=1NSwobs,i−Swsim,i2,R2=1−SESEobs,RMSE=SEn,RE=WSi−WmWSr×100%.In the (12), Swsim,i is the simulated value which is obtained through forward modeling by measured parameters or optimization parameters; Swobs,i is the observed value; Swobs¯ and Swsim¯ are the average values of Swobs,i and Swsim,i, respectively; N is the number of sample which is 9 in this paper; SEobs indicates the sum of squares of the total observations; SE is the sum of squares for error. WSi is the simulated value of the discharge of the surface runoff; Wm is the measured value of the discharge of the surface runoff; and WSr is the precipitation of the experimental plot.


## 2.6. Calibration and Parametrization of Models

Experiments S1 and S4 were chosen for the model calibration because of the completeness of the dataset. Specifically, experiment S1 was used to calibrate the saturation, while experiment S4 was to calibrate the surface runoff parameters. The calibration was conducted manually by adjusting the model parameters (one at a time) to reduce the residual between the predictive and the measured values. The model was forced by prescribing the observed rainfall. Model evaluation parameters which were calibrated against soil saturation and discharge indicate that the calibrated parameters could be reliable enough for model calibration, as shown in Figure3 and Table 3.Figure 3
Comparison of the measured and SMC-simulated value of experiment S1.
(a)
Comparison of soil saturation
(b)
Comparison of dischargeTable 3
Model evaluation parameters of the experiment S1.
Model calibration
R2
RMSE
RE
Value
0.87
0.039
0.2%The calibrated model parameter values (as shown in Table4) achieved the best which indicate satisfactory performances of the model as suggested by Zhu et al. [23].Table 4
Hydrodynamic parameters of the soil matrix and surface runoff parameters.
Bulk density
Porosity
Saturated hydraulic conductivity
Residual moisture
Saturated moisture
Shape parameter
Reciprocal of inhale
Manning coefficient
g·cm−3
%
Ks (m·h−1)
θr
θs
n
α (m−1)
min·cm1/3
1.58
46.0
0.005
0.034
0.45
1.37
0.01
0.0000001


## 3. Results and Discussion

### 3.1. The Effect of the Fracture Network on Discharge

As shown in Figure4, the initial runoff-yielding times of experiments S2, S3, and S5 are 741 s, 819 s, and 350 s, respectively, and the discharges of the surface runoff are 11,790 mL, 6500 mL, and 19,650 mL, respectively, when the precipitation intensity is 50 mm/h. Because of the existence of fractures, infiltrated water will reach the groundwater faster along these cracks in experiments S1 and S2, thus reducing the discharge. When the precipitation intensity is 30 mm/h, the initial runoff-yielding times of experiments S1 and S4 are 3147 s and 1279 s and the discharges are 815 mL and 7960 mL, respectively. It shows that the existence of the fracture network delays the initial runoff yielding of slope flow and reduces the discharge.Figure 4
Surface runoff in comparison experiment.
(a)
Surface runoff for 30 mm/h
(b)
Surface runoff for 50 mm/h


### 3.2. The Effect of the Fracture Network on the Soil Water Movement

#### 3.2.1. The Effect of the Fracture Network on the Water Infiltration

In order to investigate the effect of the fracture network on soil water movement, two methods were used to calculate the accumulated infiltration, that is, water balance method and soil water content method. When the water content method was used, the difference between the experimental soil water content and the initial soil water content is obtained by layer accumulation to calculate the accumulated infiltration in the 9 sections. Then, they were added to the scatter diagram, with the fitting curve of the accumulated infiltration drawn. The integral of the curve is the experimental accumulated infiltration. Because the sampling points are relatively large, they cannot represent the water content of the whole section. In addition, due to the experimental error, the calculation results using the soil water content method are smaller than the results of water balance calculation.As shown in Table5, the accumulated infiltration of the experiments which had cracks (S2 and S3) was significantly higher than those in the noncrack experiments (S5) when the precipitation intensity is 50 mm/h, which indicates that the fracture network can increase the soil water infiltration. When the precipitation intensity is 30 mm/h, it leads to the same conclusion but the variable quantity of the accumulated infiltration was not as significant as that when the precipitation intensity is 50 mm/h. It can be concluded that the fracture network can increase water infiltration, with higher influence of the fracture network on infiltration due to larger precipitation intensity.Table 5
Accumulated infiltration of different methods at different experiments.
Experiment label
Water balance method (mL)
Water content method (mL)
R2
S1
29,185
24,665
0.7111
S2
13,200
11,681
0.6872
S3
18,500
12,852
0.9887
S4
22,040
15,146
0.9972
S5
5350
5247
0.9983


#### 3.2.2. The Influence of the Fracture Network on the Soil Water Infiltration Depth

The depth of staining of a section was obtained by calculating the number of pixel points on each vertical pixel of an infiltrating image, and the depth of infiltration on each section was averaged per 10 cm in the horizontal direction. The infiltration depth was the depth of penetration on each vertical pixel line, and then, the distribution of the infiltration depth of the experiment was obtained. Figure5 shows the spatial distribution of the infiltration depth on the Z-X plane for different experiments.Figure 5
The distribution of the infiltration depth.
(a)
Experiment S1
(b)
Experiment S2
(c)
Experiment S3
(d)
Experiment S4
(e)
Experiment S5The infiltration shape whose experiments had cracks (S1, S2, and S3) showed zonal distribution (as shown in Figure5). The average infiltration depth of experiments S2 and S3 is 33.16 cm and 32.49 cm, respectively, while the average infiltration depth of experiments S1, S4, and S5 is 18.45 cm, 13.69 cm, and 11.64 cm, respectively. It represents that the existence of the fracture network increases the depth of the average water infiltration. What is more, the depth of the average water infiltration is increasing as the precipitation intensity increases.The variance of infiltration depth was further analyzed to compare the infiltration uniformity of the sections. The variance is written as follows:(13)S=1n−1∑i=1nxi−x¯2,where xi is the infiltration depth of each vertical pixel line of the experiment and x¯ is the depth of the average infiltration of the experiment. The variance of experiments S1, S2, S3, S4, and S5, respectively, is 25.2, 74.9, 108.4, 5.6, and 8.6, which indicates that the uniformity of soil water infiltration is affected by cracks, and the values of variance for those experiments that have cracks are higher than those in the noncrack experiments. In addition, high precipitation intensity will decrease the uniformity of soil water infiltration.


#### 3.2.3. The Effect of the Fracture Network on the Amount of Stained Area per Depth (Dye Coverage)

As shown in Figure6, the soil water infiltration area is the area of a closed curve surrounded by curves and axes. The average dye coverage of experiments S1 to S5 was 30.7%, 55.8%, 54.2%, 22.8%, and 19.0%, respectively. These results show that the existence of the fracture network increases the soil water infiltration when the precipitation intensity is the same.Figure 6
The amount of stained area per depth (dye coverage) in different experiments; P1 to P9 represent sections 1 to 9.
(a)
Experiment S1
(b)
Experiment S2
(c)
Experiment S3
(d)
Experiment S4
(e)
Experiment S5Comparing different precipitation intensities (50 mm/h and 30 mm/h) for the dye coverage of the experiment which has cracks, it can be found that the fracture network obviously increases the dye coverage and the dye coverage is increasing with the increase of precipitation intensity, indicating that an increase of rainfall intensity may strengthen the flow channel, which is favorable for the formation of preferential flow.


### 3.3. Comparison between Coefficient of Infiltration and the Runoff Coefficient

Coefficient of infiltration refers to the ratio of the accumulation infiltration to the precipitation calculated by the moisture content, while the runoff coefficient is the ratio of the amount discharge to the amount precipitation. As shown in Table6, the coefficient of infiltration and the runoff coefficient of experiment S5 with no cracks are 0.786 and 0.2099, respectively, when the precipitation intensity is 50 mm/h, and the coefficient of infiltration increases to 0.4672 and 0.5141 for experiments S2 and S3, while the runoff coefficient decreases to 0.4716 and 0.26, respectively. The coefficient of infiltration and the runoff coefficient of experiment S4 (no cracks) are 0.2653 and 0.5049 when the precipitation intensity is 30 mm/h, and the runoff coefficient of the experiment S1 decreases to 0.0272, while the coefficient of infiltration increases to 0.8222. It indicates that the existence of the fracture network decreases the discharge and increases the soil water infiltration.Table 6
Runoff coefficient and the coefficient of rainfall infiltration.
Experiment label
Discharge (mL)
Infiltration capacity (mL)
Runoff coefficient
Coefficient of rainfall infiltration
S1
815
24,665
0.0272
0.8222
S2
11,790
11,681
0.4716
0.4672
S3
6500
12,852
0.26
0.5141
S4
7960
15,146
0.2653
0.5049
S5
19,650
5247
0.786
0.2099Furthermore, the average infiltration depth of the experiments with cracks was 21.34 cm higher than those with no cracks, and the average dyeing area increased by 36.0% when the precipitation intensity was 50 mm/h; in addition, they were 5.7 cm and 9.5%, respectively, when precipitation intensity was 30 mm/h, indicating that the coefficient of rainfall infiltration decreases as the rainfall intensity increases, but the average infiltration depth and the average dyeing area increase. This behavior shows that soil water flows mainly through the matrix if the intensity of rainfall or irrigation is low; otherwise, it flows preferentially through the cracks.


### 3.4. Comparison between Simulated and Measured Soil Moisture Content

In order to describe preferential flow motion clearly, the saturation variation between cross sectional saturation and the initial saturation indicates the amount of changes on saturation as the display indicator shows. Moreover, the simulation area of the five experiments was a hexahedron consisting of parallel quadrilaterals, and therefore, the soil water potential is different in different spatial dimensions. In order to further reduce the effect of the soil water potential, the same section of different simulated experiments was selected for comparison. Figure7 shows the result of each experiment using the SMC model when t = 12 h.Figure 7
Simulated flow pattern of each experiment by the SMC model whent = 12 h.
(a)
Experiment S1
(b)
Experiment S2
(c)
Experiment S3Of note, difference of the cross sectional saturation and the initial saturation is the saturation variation.In order to compare the SMC model and the traditional Richards’ model, results from the two models were compared in the same position through saturation variation. Figure8 shows that the SMC model is capable of simulating the rapid infiltration of water in the fracture network. Selecting the same section whose flow at x = 0.1 m of experiments S1, S2, and S3 for comparison (as shown in Figure 9 and Table 7), it shows that simulated soil water in the fracture network reached the groundwater quickly and directly. Compared with the traditional Richards’ model (1), the SMC model improves the model accuracy because it takes into consideration the existence of preferential flow. Figure 10 shows the simulated movement of water at different times during the experiment, which indicates that the model is capable of simulating preferential flow.Figure 8
Comparison of the flow pattern simulated by the Richards and SMC model whent = 12 h.
(a)
SMC model
(b)
Richards’ modelFigure 9
Comparison of saturation between measured values and simulated values of each experiment.
(a)
Experiment S1
(b)
Experiment S2
(c)
Experiment S3Table 7
Comparison of simulation results by different models.
Experiment label
S1
S2
S3
Statistic indexes
SMC model
Richards’ model
SMC model
Richards’ model
SMC model
Richards’ model
R2
0.80
0.47
0.33
0.14
0.72
0.59
RMSE
0.044
0.065
0.037
0.049
0.04
0.056Figure 10
SMC simulation of soil water movement of the experiment S1.
(a)
t = 1 h
(b)
t = 4 h
(c)
t = 8 h
(d)
t = 12 hIn general, the results indicate that the SMC model could simulate preferential flow more effectively and accurately than the model based on Richards’ equation.


### 3.5. The Comparison of Experimental Data and Simulated Curve of the Discharge

Figure11 shows that the SMC model can estimate the discharge at a high accuracy compared against the observed values, with relative error of 9.71% and 5.0% for experiments S2 and S3. This may result from the fact that our model considers the existence of the fracture networks. It also indicates that the model can simulate the overall flow of surface runoff, and the existence of the cracks may delay the production of surface runoff.Figure 11
Comparison of the measured and SMC-simulated value of the discharge.
(a)
Experiment S2
(b)
Experiment S3


## 3.1. The Effect of the Fracture Network on Discharge

As shown in Figure4, the initial runoff-yielding times of experiments S2, S3, and S5 are 741 s, 819 s, and 350 s, respectively, and the discharges of the surface runoff are 11,790 mL, 6500 mL, and 19,650 mL, respectively, when the precipitation intensity is 50 mm/h. Because of the existence of fractures, infiltrated water will reach the groundwater faster along these cracks in experiments S1 and S2, thus reducing the discharge. When the precipitation intensity is 30 mm/h, the initial runoff-yielding times of experiments S1 and S4 are 3147 s and 1279 s and the discharges are 815 mL and 7960 mL, respectively. It shows that the existence of the fracture network delays the initial runoff yielding of slope flow and reduces the discharge.Figure 4
Surface runoff in comparison experiment.
(a)
Surface runoff for 30 mm/h
(b)
Surface runoff for 50 mm/h


## 3.2. The Effect of the Fracture Network on the Soil Water Movement

### 3.2.1. The Effect of the Fracture Network on the Water Infiltration

In order to investigate the effect of the fracture network on soil water movement, two methods were used to calculate the accumulated infiltration, that is, water balance method and soil water content method. When the water content method was used, the difference between the experimental soil water content and the initial soil water content is obtained by layer accumulation to calculate the accumulated infiltration in the 9 sections. Then, they were added to the scatter diagram, with the fitting curve of the accumulated infiltration drawn. The integral of the curve is the experimental accumulated infiltration. Because the sampling points are relatively large, they cannot represent the water content of the whole section. In addition, due to the experimental error, the calculation results using the soil water content method are smaller than the results of water balance calculation.As shown in Table5, the accumulated infiltration of the experiments which had cracks (S2 and S3) was significantly higher than those in the noncrack experiments (S5) when the precipitation intensity is 50 mm/h, which indicates that the fracture network can increase the soil water infiltration. When the precipitation intensity is 30 mm/h, it leads to the same conclusion but the variable quantity of the accumulated infiltration was not as significant as that when the precipitation intensity is 50 mm/h. It can be concluded that the fracture network can increase water infiltration, with higher influence of the fracture network on infiltration due to larger precipitation intensity.Table 5
Accumulated infiltration of different methods at different experiments.
Experiment label
Water balance method (mL)
Water content method (mL)
R2
S1
29,185
24,665
0.7111
S2
13,200
11,681
0.6872
S3
18,500
12,852
0.9887
S4
22,040
15,146
0.9972
S5
5350
5247
0.9983


### 3.2.2. The Influence of the Fracture Network on the Soil Water Infiltration Depth

The depth of staining of a section was obtained by calculating the number of pixel points on each vertical pixel of an infiltrating image, and the depth of infiltration on each section was averaged per 10 cm in the horizontal direction. The infiltration depth was the depth of penetration on each vertical pixel line, and then, the distribution of the infiltration depth of the experiment was obtained. Figure5 shows the spatial distribution of the infiltration depth on the Z-X plane for different experiments.Figure 5
The distribution of the infiltration depth.
(a)
Experiment S1
(b)
Experiment S2
(c)
Experiment S3
(d)
Experiment S4
(e)
Experiment S5The infiltration shape whose experiments had cracks (S1, S2, and S3) showed zonal distribution (as shown in Figure5). The average infiltration depth of experiments S2 and S3 is 33.16 cm and 32.49 cm, respectively, while the average infiltration depth of experiments S1, S4, and S5 is 18.45 cm, 13.69 cm, and 11.64 cm, respectively. It represents that the existence of the fracture network increases the depth of the average water infiltration. What is more, the depth of the average water infiltration is increasing as the precipitation intensity increases.The variance of infiltration depth was further analyzed to compare the infiltration uniformity of the sections. The variance is written as follows:(13)S=1n−1∑i=1nxi−x¯2,where xi is the infiltration depth of each vertical pixel line of the experiment and x¯ is the depth of the average infiltration of the experiment. The variance of experiments S1, S2, S3, S4, and S5, respectively, is 25.2, 74.9, 108.4, 5.6, and 8.6, which indicates that the uniformity of soil water infiltration is affected by cracks, and the values of variance for those experiments that have cracks are higher than those in the noncrack experiments. In addition, high precipitation intensity will decrease the uniformity of soil water infiltration.


### 3.2.3. The Effect of the Fracture Network on the Amount of Stained Area per Depth (Dye Coverage)

As shown in Figure6, the soil water infiltration area is the area of a closed curve surrounded by curves and axes. The average dye coverage of experiments S1 to S5 was 30.7%, 55.8%, 54.2%, 22.8%, and 19.0%, respectively. These results show that the existence of the fracture network increases the soil water infiltration when the precipitation intensity is the same.Figure 6
The amount of stained area per depth (dye coverage) in different experiments; P1 to P9 represent sections 1 to 9.
(a)
Experiment S1
(b)
Experiment S2
(c)
Experiment S3
(d)
Experiment S4
(e)
Experiment S5Comparing different precipitation intensities (50 mm/h and 30 mm/h) for the dye coverage of the experiment which has cracks, it can be found that the fracture network obviously increases the dye coverage and the dye coverage is increasing with the increase of precipitation intensity, indicating that an increase of rainfall intensity may strengthen the flow channel, which is favorable for the formation of preferential flow.


## 3.2.1. The Effect of the Fracture Network on the Water Infiltration

In order to investigate the effect of the fracture network on soil water movement, two methods were used to calculate the accumulated infiltration, that is, water balance method and soil water content method. When the water content method was used, the difference between the experimental soil water content and the initial soil water content is obtained by layer accumulation to calculate the accumulated infiltration in the 9 sections. Then, they were added to the scatter diagram, with the fitting curve of the accumulated infiltration drawn. The integral of the curve is the experimental accumulated infiltration. Because the sampling points are relatively large, they cannot represent the water content of the whole section. In addition, due to the experimental error, the calculation results using the soil water content method are smaller than the results of water balance calculation.As shown in Table5, the accumulated infiltration of the experiments which had cracks (S2 and S3) was significantly higher than those in the noncrack experiments (S5) when the precipitation intensity is 50 mm/h, which indicates that the fracture network can increase the soil water infiltration. When the precipitation intensity is 30 mm/h, it leads to the same conclusion but the variable quantity of the accumulated infiltration was not as significant as that when the precipitation intensity is 50 mm/h. It can be concluded that the fracture network can increase water infiltration, with higher influence of the fracture network on infiltration due to larger precipitation intensity.Table 5
Accumulated infiltration of different methods at different experiments.
Experiment label
Water balance method (mL)
Water content method (mL)
R2
S1
29,185
24,665
0.7111
S2
13,200
11,681
0.6872
S3
18,500
12,852
0.9887
S4
22,040
15,146
0.9972
S5
5350
5247
0.9983


## 3.2.2. The Influence of the Fracture Network on the Soil Water Infiltration Depth

The depth of staining of a section was obtained by calculating the number of pixel points on each vertical pixel of an infiltrating image, and the depth of infiltration on each section was averaged per 10 cm in the horizontal direction. The infiltration depth was the depth of penetration on each vertical pixel line, and then, the distribution of the infiltration depth of the experiment was obtained. Figure5 shows the spatial distribution of the infiltration depth on the Z-X plane for different experiments.Figure 5
The distribution of the infiltration depth.
(a)
Experiment S1
(b)
Experiment S2
(c)
Experiment S3
(d)
Experiment S4
(e)
Experiment S5The infiltration shape whose experiments had cracks (S1, S2, and S3) showed zonal distribution (as shown in Figure5). The average infiltration depth of experiments S2 and S3 is 33.16 cm and 32.49 cm, respectively, while the average infiltration depth of experiments S1, S4, and S5 is 18.45 cm, 13.69 cm, and 11.64 cm, respectively. It represents that the existence of the fracture network increases the depth of the average water infiltration. What is more, the depth of the average water infiltration is increasing as the precipitation intensity increases.The variance of infiltration depth was further analyzed to compare the infiltration uniformity of the sections. The variance is written as follows:(13)S=1n−1∑i=1nxi−x¯2,where xi is the infiltration depth of each vertical pixel line of the experiment and x¯ is the depth of the average infiltration of the experiment. The variance of experiments S1, S2, S3, S4, and S5, respectively, is 25.2, 74.9, 108.4, 5.6, and 8.6, which indicates that the uniformity of soil water infiltration is affected by cracks, and the values of variance for those experiments that have cracks are higher than those in the noncrack experiments. In addition, high precipitation intensity will decrease the uniformity of soil water infiltration.


## 3.2.3. The Effect of the Fracture Network on the Amount of Stained Area per Depth (Dye Coverage)

As shown in Figure6, the soil water infiltration area is the area of a closed curve surrounded by curves and axes. The average dye coverage of experiments S1 to S5 was 30.7%, 55.8%, 54.2%, 22.8%, and 19.0%, respectively. These results show that the existence of the fracture network increases the soil water infiltration when the precipitation intensity is the same.Figure 6
The amount of stained area per depth (dye coverage) in different experiments; P1 to P9 represent sections 1 to 9.
(a)
Experiment S1
(b)
Experiment S2
(c)
Experiment S3
(d)
Experiment S4
(e)
Experiment S5Comparing different precipitation intensities (50 mm/h and 30 mm/h) for the dye coverage of the experiment which has cracks, it can be found that the fracture network obviously increases the dye coverage and the dye coverage is increasing with the increase of precipitation intensity, indicating that an increase of rainfall intensity may strengthen the flow channel, which is favorable for the formation of preferential flow.


## 3.3. Comparison between Coefficient of Infiltration and the Runoff Coefficient

Coefficient of infiltration refers to the ratio of the accumulation infiltration to the precipitation calculated by the moisture content, while the runoff coefficient is the ratio of the amount discharge to the amount precipitation. As shown in Table6, the coefficient of infiltration and the runoff coefficient of experiment S5 with no cracks are 0.786 and 0.2099, respectively, when the precipitation intensity is 50 mm/h, and the coefficient of infiltration increases to 0.4672 and 0.5141 for experiments S2 and S3, while the runoff coefficient decreases to 0.4716 and 0.26, respectively. The coefficient of infiltration and the runoff coefficient of experiment S4 (no cracks) are 0.2653 and 0.5049 when the precipitation intensity is 30 mm/h, and the runoff coefficient of the experiment S1 decreases to 0.0272, while the coefficient of infiltration increases to 0.8222. It indicates that the existence of the fracture network decreases the discharge and increases the soil water infiltration.Table 6
Runoff coefficient and the coefficient of rainfall infiltration.
Experiment label
Discharge (mL)
Infiltration capacity (mL)
Runoff coefficient
Coefficient of rainfall infiltration
S1
815
24,665
0.0272
0.8222
S2
11,790
11,681
0.4716
0.4672
S3
6500
12,852
0.26
0.5141
S4
7960
15,146
0.2653
0.5049
S5
19,650
5247
0.786
0.2099Furthermore, the average infiltration depth of the experiments with cracks was 21.34 cm higher than those with no cracks, and the average dyeing area increased by 36.0% when the precipitation intensity was 50 mm/h; in addition, they were 5.7 cm and 9.5%, respectively, when precipitation intensity was 30 mm/h, indicating that the coefficient of rainfall infiltration decreases as the rainfall intensity increases, but the average infiltration depth and the average dyeing area increase. This behavior shows that soil water flows mainly through the matrix if the intensity of rainfall or irrigation is low; otherwise, it flows preferentially through the cracks.


## 3.4. Comparison between Simulated and Measured Soil Moisture Content

In order to describe preferential flow motion clearly, the saturation variation between cross sectional saturation and the initial saturation indicates the amount of changes on saturation as the display indicator shows. Moreover, the simulation area of the five experiments was a hexahedron consisting of parallel quadrilaterals, and therefore, the soil water potential is different in different spatial dimensions. In order to further reduce the effect of the soil water potential, the same section of different simulated experiments was selected for comparison. Figure7 shows the result of each experiment using the SMC model when t = 12 h.Figure 7
Simulated flow pattern of each experiment by the SMC model whent = 12 h.
(a)
Experiment S1
(b)
Experiment S2
(c)
Experiment S3Of note, difference of the cross sectional saturation and the initial saturation is the saturation variation.In order to compare the SMC model and the traditional Richards’ model, results from the two models were compared in the same position through saturation variation. Figure8 shows that the SMC model is capable of simulating the rapid infiltration of water in the fracture network. Selecting the same section whose flow at x = 0.1 m of experiments S1, S2, and S3 for comparison (as shown in Figure 9 and Table 7), it shows that simulated soil water in the fracture network reached the groundwater quickly and directly. Compared with the traditional Richards’ model (1), the SMC model improves the model accuracy because it takes into consideration the existence of preferential flow. Figure 10 shows the simulated movement of water at different times during the experiment, which indicates that the model is capable of simulating preferential flow.Figure 8
Comparison of the flow pattern simulated by the Richards and SMC model whent = 12 h.
(a)
SMC model
(b)
Richards’ modelFigure 9
Comparison of saturation between measured values and simulated values of each experiment.
(a)
Experiment S1
(b)
Experiment S2
(c)
Experiment S3Table 7
Comparison of simulation results by different models.
Experiment label
S1
S2
S3
Statistic indexes
SMC model
Richards’ model
SMC model
Richards’ model
SMC model
Richards’ model
R2
0.80
0.47
0.33
0.14
0.72
0.59
RMSE
0.044
0.065
0.037
0.049
0.04
0.056Figure 10
SMC simulation of soil water movement of the experiment S1.
(a)
t = 1 h
(b)
t = 4 h
(c)
t = 8 h
(d)
t = 12 hIn general, the results indicate that the SMC model could simulate preferential flow more effectively and accurately than the model based on Richards’ equation.


## 3.5. The Comparison of Experimental Data and Simulated Curve of the Discharge

Figure11 shows that the SMC model can estimate the discharge at a high accuracy compared against the observed values, with relative error of 9.71% and 5.0% for experiments S2 and S3. This may result from the fact that our model considers the existence of the fracture networks. It also indicates that the model can simulate the overall flow of surface runoff, and the existence of the cracks may delay the production of surface runoff.Figure 11
Comparison of the measured and SMC-simulated value of the discharge.
(a)
Experiment S2
(b)
Experiment S3


## 4. Conclusions

In this paper, we investigated the surface runoff and the soil water movement in the presence of the fracture network. Our work leads to the following conclusions.(1)
Part of the water is infiltrated through cracks which put off the initial runoff-yielding time of the surface runoff, reducing the discharge of surface runoff and increasing the soil water infiltration. When the intensity of precipitation or irrigation is low, the soil water flow moves mainly in soil matrix; otherwise, it is mainly in terms of preferential flow through cracks. The nonuniformity of soil water infiltration and the depth of the average water infiltration increases as the precipitation intensity increases.
(2)
The complete coupling model (SMC model) is established, which includes the traditional Richards’ equation, the “cubic law” of fluid motion in cracks within smooth parallel plates’ equation, as well as the two-dimensional diffusion wave approximation to Saint-Venant equations. The model divides soil water infiltration into two forms and uses the overall method to calculate exchange water between the surface runoff and the soil water movement.
(3)
According to the coefficient of determinationR2 and the root mean square error (RMSE) (discussed in Section 3.4) and the simulated saturation which were, respectively, simulated by the SMC model and the Richards’ models as shown in Table 7 and Figure 9, the simulated saturation values of the SMC model were better because the SMC model considers the movement of preferential flow. It indicates that the SMC model has better performance compared with the Richards’ model. In addition, the SMC model can simulate the surface runoff movement, and therefore, the complete coupling model can be used to describe the water content distribution of all domains, as well as the dynamic process of infiltration and surface runoff.In this paper, the artificial soil fracture network and the hillslope system of a relatively large slope were utilized for our analysis. However, in field condition, the distribution of soil cracks is not clear and the soil surface is rough. Therefore, the extent of the fracture network, the spatial structure, and the soil surface roughness of the fracture will be considered in our further work.



---

*Source: 7056858-2018-07-15.xml*