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

Zhu, Lei; Fan, Dongjun; Ma, Rong; Zhang, Yonggen; Zha, Yuanyuan

doi:https://doi.org/10.1155/2018/7056858

Geofluids

On this page

Abstract Introduction Materials and Methods Results and Discussion Conclusions Additional Points Conflicts of Interest Acknowledgments References Copyright Related Articles

Research Article | Open Access

Volume 2018 | Article ID 7056858 | https://doi.org/10.1155/2018/7056858

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

Lei Zhu,^1,2Dongjun Fan,^1,2Rong Ma,^1,2Yonggen Zhang,³and Yuanyuan Zha⁴

Academic Editor: Stefano Lo Russo

Received09 Oct 2017

Revised24 Apr 2018

Accepted20 May 2018

Published15 Jul 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.

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 Table 1.

2.1.2. Experimental Design

According to different rainfall intensities and durations, the experiment is divided into five groups that are shown in Table 2. 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).

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 Figure 1. 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 the -axis direction, with the interval spacing of 5–10 cm between two vertical profiles.

For the tracer experiments, Fe³⁺ 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(NO³)³). 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.

(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]: and where is the saturated moisture content, ; is the soil saturation; is the water channel saturation, dimensionless; is the flow channel width with indicating the radius of the channel, ; and represent the soil water matrix potential and the flow channel water potential, ; is the elevation head and is the elevation head within the fracture, ; is the root water uptake term or other sources and sinks, ; is spatial coordinates, denotes horizontal coordinate, and is the vertical coordinate and is positive upward; is time, ; and represent the saturated hydraulic conductivity of soil and the water channel, ; and represent the relative hydraulic conductivity of soil and the water channel, dimensionless; is the fracture water and soil water exchange capacity, ; and is the soil water and surface water exchange capacity, . 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]: where is the surface altitude, ; ; is the ground altitude, ; is the water depth, ; is the source-sink item, ; and are land surface conduction coefficient in and directions. For surface runoff model parameters, and can be calculated, respectively, by the two formulas as follows: where and are Manning roughness coefficient in direction and direction, .

2.3.3. Fracture Network Coupling with the Matrix Region

The discrete equation of soil water movement in the soil is described as follows: where is a factor for impact volume node , ; is the water potential difference between the nodes and matrix area, ; represents a relative unsaturated hydraulic conductivity between the node and node when using upstream weighting; represents the order junction node when the region is discrete; is the coefficient matrix; refers to the flow of flux on the border near the border node ; outflow is positive, for the internal node , ; represents the study area boundary, ; and is the corresponding basis function in a unit.

The discrete equation of preferential water movement within soil is described as follows: where is an impact factor for node , ; represents the node order of junction connect when region is discrete; is head difference between the node and node , ; indicates the relative conductivity between the node to node while using the weighted upwind method, dimensionless; and is the coefficient matrix. is the exchange water of the two regions.

Replacing the exact value of the fluid exchange into (5) gives the following two-domain coupling equation:

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 has double character of the two zones if the node is shared by soil and surface.

Presume , for the soil water flow movement discrete equation of the node , contains water exchange item of the two zones. But for the discrete equation of surface, is the exchange item of the two zones; (8) can be shown as follows:

Substituting (9) into (8) leads to the integrated equation of two domains:

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 Table 2, with distribution of the precipitation evenly. The discharge () per unit width at the critical depth boundary is given by the following:

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 determination and root mean square error () are used to evaluate soil moisture prediction, while relative error () is used to assess the surface runoff prediction performance. The mathematical expressions for , , and are

In the (12), is the simulated value which is obtained through forward modeling by measured parameters or optimization parameters; is the observed value; and are the average values of and , respectively; is the number of sample which is 9 in this paper; indicates the sum of squares of the total observations; is the sum of squares for error. is the simulated value of the discharge of the surface runoff; is the measured value of the discharge of the surface runoff; and 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 Figure 3 and Table 3.

(a) Comparison of soil saturation

(b) Comparison of discharge

The calibrated model parameter values (as shown in Table 4) achieved the best which indicate satisfactory performances of the model as suggested by Zhu et al. [23].

3. Results and Discussion

3.1. The Effect of the Fracture Network on Discharge

As shown in Figure 4, 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.

(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 Table 5, 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.

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. Figure 5 shows the spatial distribution of the infiltration depth on the Z-X plane for different experiments.

(a) Experiment S1

(b) Experiment S2

(c) Experiment S3

(d) Experiment S4

(e) Experiment S5

The infiltration shape whose experiments had cracks (S1, S2, and S3) showed zonal distribution (as shown in Figure 5). 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: where is the infiltration depth of each vertical pixel line of the experiment and 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 Figure 6, 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.

(a) Experiment S1

(b) Experiment S2

(c) Experiment S3

(d) Experiment S4

(e) Experiment S5

Comparing 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 Table 6, 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.

Furthermore, 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. Figure 7 shows the result of each experiment using the SMC model when = 12 h.

(a) Experiment S1

(b) Experiment S2

(c) Experiment S3

Of 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. Figure 8 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 = 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.

(a) SMC model

(b) Richards’ model

(a) Experiment S1

(b) Experiment S2

(c) Experiment S3

(a) = 1 h

(b) = 4 h

(c) = 8 h

(d) = 12 h

In 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

Figure 11 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.

(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 determination and the root mean square error () (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.

Additional Points

Highlights. (i) Dye-tracing experiments were conducted to investigate the surface runoff and the soil water movement in the presence of the fracture network. (ii) The effect of the fracture network on the soil water and overland flow was analyzed under different precipitation intensities and durations. (iii) Numerical model was built to simulate the coupled surface water flow and subsurface flow in soil crack networks and in soil matrix.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This research was financially supported by grants of the National Science Foundation of China (nos. 51569023 and 51369025) and the First Class Discipline Construction in Ningxia (no. NXYLXK2017A03). National Natural Science Foundation of China (51609173 and 51779179).

References

R. A. Stirling, S. Glendinning, and C. T. Davie, “Modelling the deterioration of the near surface caused by drying induced cracking,” Applied Clay Science, vol. 146, pp. 176–185, 2017.
View at: Publisher Site | Google Scholar
H.-z. Qiu, J.-m. Kong, R.-c. Wang, Y. Cui, and S.-w. Huang, “Response mechanism of post-earthquake slopes under heavy rainfall,” Journal of Seismology, vol. 21, no. 4, pp. 869–884, 2017.
View at: Publisher Site | Google Scholar
S. Krisnanto, H. Rahardjo, D. G. Fredlund, and E. C. Leong, “Water content of soil matrix during lateral water flow through cracked soil,” Engineering Geology, vol. 210, pp. 168–179, 2016.
View at: Publisher Site | Google Scholar
M. S. Kim, Y. Onda, T. Uchida, J. K. Kim, and Y. S. Song, “Effect of seepage on shallow landslides in consideration of changes in topography: case study including an experimental sandy slope with artificial rainfall,” Catena, vol. 161, pp. 50–62, 2018.
View at: Publisher Site | Google Scholar
S. Verma, S. K. Mishra, A. Singh, P. K. Singh, and R. K. Verma, “An enhanced SMA based SCS-CN inspired model for watershed runoff prediction,” Environmental Earth Sciences, vol. 76, no. 21, 2017.
View at: Publisher Site | Google Scholar
A. G. Awadallah, H. Saad, A. Elmoustafa, and A. Hassan, “Reliability assessment of water structures subject to data scarcity using the SCS-CN model,” Hydrological Sciences Journal, vol. 61, no. 4, pp. 696–710, 2016.
View at: Publisher Site | Google Scholar
C. Cruz and E. Miranda, “Evaluation of soil-structure interaction effects on the damping ratios of buildings subjected to earthquakes,” Soil Dynamics and Earthquake Engineering, vol. 100, pp. 183–195, 2017.
View at: Publisher Site | Google Scholar
J. F. O'Callaghan and D. M. Mark, “The extraction of drainage networks from digital elevation data,” Computer Vision, Graphics, and Image Processing, vol. 28, no. 3, pp. 323–344, 1984.
View at: Publisher Site | Google Scholar
S. Sayari, M. Rahimpour, and M. Zounemat-Kermani, “Numerical modeling based on a finite element method for simulation of flow in furrow irrigation,” Paddy and Water Environment, vol. 15, no. 4, pp. 879–887, 2017.
View at: Publisher Site | Google Scholar
H. Vergara, P.-E. Kirstetter, J. J. Gourley et al., “Estimating a-priori kinematic wave model parameters based on regionalization for flash flood forecasting in the conterminous United States,” Journal of Hydrology, vol. 541, pp. 421–433, 2016.
View at: Publisher Site | Google Scholar
L. Warsta, T. Karvonen, H. Koivusalo, M. Paasonen-Kivekäs, and A. Taskinen, “Simulation of water balance in a clayey, subsurface drained agricultural field with three-dimensional FLUSH model,” Journal of Hydrology, vol. 476, pp. 395–409, 2013.
View at: Publisher Site | Google Scholar
L. Ma, C. He, H. Bian, and L. Sheng, “MIKE SHE modeling of ecohydrological processes: merits, applications, and challenges,” Ecological Engineering, vol. 96, pp. 137–149, 2016.
View at: Publisher Site | Google Scholar
E. Dimitriou, E. Moussoulis, C. Díaz-Paniagua, and L. Serrano, “Hydrodynamic numerical modelling of the water level decline in four temporary ponds of the Doñana National Park (SW Spain),” Journal of Arid Environments, vol. 147, pp. 90–102, 2017.
View at: Publisher Site | Google Scholar
L. Locatelli, O. Mark, P. S. Mikkelsen et al., “Hydrologic impact of urbanization with extensive stormwater infiltration,” Journal of Hydrology, vol. 544, pp. 524–537, 2017.
View at: Publisher Site | Google Scholar
V. Singh and S. Murty Bhallamudi, “Conjunctive surface-subsurface modeling of overland flow,” Advances in Water Resources, vol. 21, no. 7, pp. 567–579, 1998.
View at: Publisher Site | Google Scholar
L. Bai, H. Li, and H. M. He, “Integrated simulation of surface water and groundwater in a high intensive coal mining area,” Journal of China Coal Society, vol. 40, no. 4, pp. 931–937, 2015.
View at: Google Scholar
C. Gandolfi and F. Savi, “A mathematical model for the coupled simulation of surface runoff and infiltration,” Journal of Agricultural Engineering Research, vol. 75, no. 1, pp. 49–55, 2000.
View at: Publisher Site | Google Scholar
P. K. Panday, C. A. Williams, K. E. Frey, and M. E. Brown, “Application and evaluation of a snowmelt runoff model in the Tamor River basin, Eastern Himalaya using a Markov Chain Monte Carlo (MCMC) data assimilation approach,” Hydrological Processes, vol. 28, no. 21, pp. 5337–5353, 2014.
View at: Publisher Site | Google Scholar
I. Trichakis, P. Burek, A. de Roo, and A. Pistocchi, “Towards a pan-European integrated groundwater and surface water model: development and applications,” Environmental Processes, vol. 4, Supplement 1, pp. 81–S93, 2017.
View at: Publisher Site | Google Scholar
J. Šimůnek, N. J. Jarvis, M. T. van Genuchten, and A. Gärdenäs, “Review and comparison of models for describing non-equilibrium and preferential flow and transport in the vadose zone,” Journal of Hydrology, vol. 272, no. 1-4, pp. 14–35, 2003.
View at: Publisher Site | Google Scholar
R. Jabakhanji and R. H. Mohtar, “A peridynamic model of flow in porous media,” Advances in Water Resources, vol. 78, pp. 22–35, 2015.
View at: Publisher Site | Google Scholar
R. D. Stewart, M. R. Abou Najm, D. E. Rupp, and J. S. Selker, “Modeling multidomain hydraulic properties of shrink-swell soils,” Water Resources Research, vol. 52, no. 10, pp. 7911–7930, 2016.
View at: Publisher Site | Google Scholar
L. Zhu, J. H. Chen, and D. D. Liu, “Preferential flow model coupling soil matrix with fracture network and its validation,” Transactions of the Chinese Society of Agricultural Engineering, vol. 32, no. 14, pp. 15–21, 2016.
View at: Google Scholar
K. Djabelkhir, C. Lauvernet, P. Kraft, and N. Carluer, “Development of a dual permeability model within a hydrological catchment modeling framework: 1D application,” Science of the Total Environment, vol. 575, pp. 1429–1437, 2017.
View at: Publisher Site | Google Scholar
E. T. De Vries, A. Raoof, and M. T. van Genuchten, “Multiscale modelling of dual-porosity porous media; a computational pore-scale study for flow and solute transport,” Advances in Water Resources, vol. 105, pp. 82–95, 2017.
View at: Publisher Site | Google Scholar
W. Shao, J. Ni, A. K. Leung, Y. Su, and C. W. W. Ng, “Analysis of plant root-induced preferential flow and pore-water pressure variation by a dual-permeability model,” Canadian Geotechnical Journal, vol. 54, no. 11, pp. 1537–1552, 2017.
View at: Publisher Site | Google Scholar
J. Bodin, P. Ackerer, A. Boisson et al., “Predictive modelling of hydraulic head responses to dipole flow experiments in a fractured/karstified limestone aquifer: insights from a comparison of five modelling approaches to real-field experiments,” Journal of Hydrology, vol. 454-455, pp. 82–100, 2012.
View at: Publisher Site | Google Scholar
A. I. Gärdenäs, J. Šimůnek, N. Jarvis, and M. T. van Genuchten, “Two-dimensional modelling of preferential water flow and pesticide transport from a tile-drained field,” Journal of Hydrology, vol. 329, no. 3-4, pp. 647–660, 2006.
View at: Publisher Site | Google Scholar
H. H. Gerke and J. Maximilian Köhne, “Dual-permeability modeling of preferential bromide leaching from a tile-drained glacial till agricultural field,” Journal of Hydrology, vol. 289, no. 1-4, pp. 239–257, 2004.
View at: Publisher Site | Google Scholar
A. Bronstert and E. J. Plate, “Modelling of runoff generation and soil moisture dynamics for hillslopes and micro-catchments,” Journal of Hydrology, vol. 198, no. 1-4, pp. 177–195, 1997.
View at: Publisher Site | Google Scholar
J. M. Kohne, S. Kohne, and J. Simunek, “Multi-process herbicide transport in structured soil columns: experiments and model analysis,” Journal of Contaminant Hydrology, vol. 85, no. 1-2, pp. 1–32, 2006.
View at: Publisher Site | Google Scholar
S. Hintikka, M. Paasonen-Kivekäs, H. Koivusalo, V. Nuutinen, and L. Alakukku, “Role of macroporosity in runoff generation on a sloping subsurface drained clay field — a case study with MACRO model,” Hydrology Research, vol. 39, no. 2, pp. 143–155, 2008.
View at: Publisher Site | Google Scholar
Y. Wang, S. A. Bradford, and J. Šimůnek, “Estimation and upscaling of dual-permeability model parameters for the transport of E. coli D21g in soils with preferential flow,” Journal of Contaminant Hydrology, vol. 159, pp. 57–66, 2014.
View at: Publisher Site | Google Scholar
R. P. Rechard, B. W. Arnold, B. A. Robinson, and J. E. Houseworth, “Transport modeling in performance assessments for the Yucca Mountain disposal system for spent nuclear fuel and high-level radioactive waste,” Reliability Engineering & System Safety, vol. 122, pp. 189–206, 2014.
View at: Publisher Site | Google Scholar
M. Turunen, L. Warsta, M. Paasonen-Kivekäs, and H. Koivusalo, “Computational assessment of sediment balance and suspended sediment transport pathways in subsurface drained clayey soils,” Soil and Tillage Research, vol. 174, pp. 58–69, 2017.
View at: Publisher Site | Google Scholar
R. Ruggenthaler, F. Schöberl, G. Markart, K. Klebinder, A. Hammerle, and G. Leitinger, “Quantification of soil moisture effects on runoff formation at the hillslope scale,” Journal of Irrigation and Drainage Engineering, vol. 141, no. 9, article 05015001, 2015.
View at: Publisher Site | Google Scholar
M. L. Sentis and C. W. Gable, “Coupling LaGrit unstructured mesh generation and model setup with TOUGH2 flow and transport: a case study,” Computers & Geosciences, vol. 108, pp. 42–49, 2017.
View at: Publisher Site | Google Scholar
M. Honari, A. Ashrafzadeh, M. Khaledian, M. Vazifedoust, and J. C. Mailhol, “Comparison of HYDRUS-3D soil moisture simulations of subsurface drip irrigation with experimental observations in the South of France,” Journal of Irrigation and Drainage Engineering, vol. 143, no. 7, article 04017014, 2017.
View at: Publisher Site | Google Scholar
H. Laine-Kaulio, S. Backnäs, T. Karvonen, H. Koivusalo, and J. J. McDonnell, “Lateral subsurface stormflow and solute transport in a forested hillslope: a combined measurement and modeling approach,” Water Resources Research, vol. 50, no. 10, pp. 8159–8178, 2014.
View at: Publisher Site | Google Scholar
D. Blessent, P. R. Jørgensen, and R. Therrien, “Comparing discrete fracture and continuum models to predict contaminant transport in fractured porous media,” Groundwater, vol. 52, no. 1, pp. 84–95, 2014.
View at: Publisher Site | Google Scholar
R. F. A. Hendriks, K. Oostindie, and P. Hamminga, “Simulation of bromide tracer and nitrogen transport in a cracked clay soil with the FLOCR/ANIMO model combination,” Journal of Hydrology, vol. 215, no. 1-4, pp. 94–115, 1999.
View at: Publisher Site | Google Scholar
F. Pierlot, J. Marks-Perreau, B. Réal et al., “Predictive quality of 26 pesticide risk indicators and one flow model: a multisite assessment for water contamination,” Science of the Total Environment, vol. 605-606, pp. 655–665, 2017.
View at: Publisher Site | Google Scholar
Q. de Jong van Lier, “Field capacity, a valid upper limit of crop available water,” Agricultural Water Management, vol. 193, pp. 214–220, 2017.
View at: Publisher Site | Google Scholar
S. Krisnanto, H. Rahardjo, D. G. Fredlund, and E. C. Leong, “Mapping of cracked soils and lateral water flow characteristics through a network of cracks,” Engineering Geology, vol. 172, pp. 12–25, 2014.
View at: Publisher Site | Google Scholar
Z. Y. Zhang, S. M. Fan, C. Wang, Y. Liu, and W. Qi, “Characteristics of unsaturated horizontal infiltration in desiccation cracks yellow cinnamon soil,” Transactions of the Chinese Society for Agricultural Machinery, vol. 48, no. 7, pp. 214–221, 2017.
View at: Google Scholar
P. S. Lang, A. Paluszny, and R. W. Zimmerman, “Permeability tensor of three-dimensional fractured porous rock and a comparison to trace map predictions,” Journal of Geophysical Research-Solid Earth, vol. 119, no. 8, pp. 6288–6307, 2014.
View at: Publisher Site | Google Scholar
J. Tang, Z. B. Peng, and Z. M. He, “Numerical analysis of fluid-solid coupling in fractured rock mass based on continuous medium,” Journal of Central South University of Science and Technology, vol. 47, no. 11, pp. 3800–3807, 2016.
View at: Google Scholar
F. Sheng, L. Y. Zhang, and D. Wu, “Review on research theories and observation techniques for preferential flow in unsaturated soil,” Transactions of the Chinese Society of Agricultural Engineering, vol. 32, no. 6, pp. 1–10, 2016.
View at: Google Scholar
A. R. Khoei, N. Hosseini, and T. Mohammadnejad, “Numerical modeling of two-phase fluid flow in deformable fractured porous media using the extended finite element method and an equivalent continuum model,” Advances in Water Resources, vol. 94, pp. 510–528, 2016.
View at: Publisher Site | Google Scholar
J. H. Li, L. M. Zhang, Y. Wang, and D. G. Fredlund, “Permeability tensor and representative elementary volume of saturated cracked soil,” Canadian Geotechnical Journal, vol. 46, no. 8, pp. 928–942, 2009.
View at: Publisher Site | Google Scholar
K. Wang, R. D. Zhang, and Z. H. Zhou, “Experiments on irrigation efficiency using iodine-starch staining,” Transactions of the Chinese Society of Agricultural Engineering, vol. 25, no. 12, pp. 38–44, 2009.
View at: Google Scholar
J. Klaus, E. Zehe, M. Elsner et al., “Controls of event-based pesticide leaching in natural soils: a systematic study based on replicated field scale irrigation experiments,” Journal of Hydrology, vol. 512, pp. 528–539, 2014.
View at: Publisher Site | Google Scholar
L. T. Zhang, G. Y. Jia, Y. M. Chen, and D. G. Fredlund, “Analytical solution for rainfall infiltration into infinite long slopes considering properties of unsaturated soil,” Chinese Journal of Geotechnical Engineering, vol. 32, no. 8, pp. 1214–1220, 2010.
View at: Google Scholar
S. Hua, S. Chen, and L. T. Liang, “Studying soil pore structure by using image filtering technology based on partial differential equation model,” Transactions of the Chinese Society of Agricultural Engineering, vol. 30, no. 3, pp. 78–85, 2014.
View at: Google Scholar
N. Otsu, “A threshold selection method from gray-level histograms,” IEEE Transactions on Systems, Man, and Cybernetics, vol. 9, no. 1, pp. 62–66, 1979.
View at: Publisher Site | Google Scholar
R. Therrien, R. G. McLaren, E. A. Sudicky, and S. M. Panday, HydroGeoSphere–A Three-Dimensional Numerical Model Describing Fully-Integrated Subsurface and Surface Flow and Solute Transport, Universite Laval and University of Waterloo, 2007.
T. Graf and R. Therrien, “Variable-density groundwater flow and solute transport in irregular 2D fracture networks,” Advances in Water Resources, vol. 30, no. 3, pp. 455–468, 2007.
View at: Publisher Site | Google Scholar
S. Panday and P. S. Huyakorn, “A fully coupled physically-based spatially distributed model for evaluating surface/subsurface flow,” Advances in Water Resources, vol. 27, no. 4, pp. 361–382, 2004.
View at: Publisher Site | Google Scholar
Z. Lei, R. Therrien, Y. Jinzhong, and Z. Qing, “Experimental and numerical simulation of water movement in soil,” Proceedings of the Institution of Civil Engineers - Water Management, vol. 165, no. 7, pp. 365–376, 2012.
View at: Publisher Site | Google Scholar
L. Zhu, J. C. Tian, and X. L. Sun, “A fully coupled numerical simulation for surface runoff and soil water movement,” Advances in Water Science, vol. 26, no. 3, pp. 322–330, 2015.
View at: Google Scholar
A. C. da Silva, R. A. Armindo, A. dos Santos Brito, and M. G. Schaap, “SPLINTEX: a physically-based pedotransfer function for modeling soil hydraulic functions,” Soil and Tillage Research, vol. 174, pp. 261–272, 2017.
View at: Publisher Site | Google Scholar

Copyright

Copyright © 2018 Lei Zhu et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

PDF Download Citation

Download other formats

Order printed copies

Views

1294

Downloads

1366

Citations