Abstract
Estimating runoff and infiltration in slopes is essential to predict flooding and landslide triggering. An approach jointly using a dual-permeability model and a kinematic wave model was proposed for integrative analyses of infiltration and surface runoff of rainwater for slopes with macropores in soils subjected to rainfalls. Laboratory tests were conducted to evaluate the performance of the proposed approach. Analyses showed that the presence of macropores greatly facilitates rainwater movement in the slope and the maximum depth of wetting front could be one order of magnitude larger than in cases without macropores. The influence of rainfall intensity, initial pore water pressure, coefficient of permeability of soil matrix, and slope angle was also examined.
1. Introduction
Slope stability failures due to short-duration high-intensity rainfalls such as seasonal typhoons and rainstorms are not uncommon in coastal areas of southeastern China. Rainfall-induced slope stability failures (e.g., landslides) are mainly attributed to the infiltration of rainwater which causes changes in pore water pressures and seepage forces [1]. Soil suctions in an unsaturated slope dissipate, totally or partially, with increasing degree of saturation of that slope. The loss of soil suctions results in reduction of shear strength of soils and could eventually lead towards failures of that slope.
Rainwater infiltration and flow within a slope have been extensively studied, for example, Akan and Yen [2], Abbott et al. [3], Garg and Sen [4], Bautista et al. [5], and Banti et al. [6]. When a slope is subjected to a rainfall event, at early stage the rainwater is absorbed by soils at and near the slope surface because of their high infiltration capacity at dry conditions. At this stage, the infiltration rate is equal to the rainfall intensity. As the infiltration continues, the infiltration capacity of the surface soils decreases rapidly and then is exceeded by the rainfall intensity. As such, slope surface runoff takes place and the infiltration rate is now equal to the infiltration capacity of the surface soils under saturated conditions. Models have been developed to describe this rainfall-infiltration-runoff process for slopes where the infiltration capacity of soils across a slope is more or less similar; for example, differences in coefficient of permeability are usually within one order of magnitude [7–10].
However, this type of models was later found inadequate in explaining the occurrence of shallow landslides for a slope subjected to a short-duration high-intensity rainstorm [1]. This is because there are always macropores in a slope commonly due to surface cracking, channels formed by roots of trees, pores caused by animals, and structural cavities in soils due to drying-wetting or freezing-thawing cycles [11]. These macropores form channels where the magnitude of soil permeability is several orders higher than that of the soil matrix (i.e., micropores). The presence of macropores allows rainwater to reach a much greater depth into the slope within a short time. This largely increases the chance for shallow landslides to take place even under short-duration rainfall conditions.
Various models have been developed to simulate rainwater flow in soils with macropores, including dual-porosity model [12], dual-permeability model [13], multidomain model [14], and two-phase model [15]. The dual-permeability model is adopted in this paper. In a dual-permeability model, soils are assumed to consist of a less permeable micropore domain (soil matrix) and a highly permeable macropore domain. The flows of rainwater in both domains are governed by two Richards’ equations. The two equations are connected through sharing a common term describing the water transfer between the two pore domains. The ability of dual-permeability models in modeling water infiltration and water flow in soils with micro- and macropores has been demonstrated (e.g., [16, 17]).
When a slope is subjected to a short-duration high-intensity rainstorm, the rainfall intensity can be reasonably assumed to be much larger than the saturated permeability of the slope surface soils. Under this situation, both flow of rainwater in soils and surface runoff take place simultaneously. The focus of this paper is on integrative analyses of surface runoff and rainwater flow in soil for slopes under rainfall conditions. The surface runoff and the rainwater flow in soil are described using one-dimensional kinematic wave equation [18] and dual-permeability model equation [16] reported in the literature, respectively. The equations are numerically solved using the finite difference method. The ability of the two equations to simulate the rainfall-infiltration-runoff process is demonstrated by the good agreement between the numerical analysis outcomes and the experimental outcomes in laboratory. A parametric study is also carried out to show the influences of rainfall intensity, initial pore water pressure, saturated soil permeability, and slope angle on the rainfall-infiltration-runoff process.
2. Analytical Models for Slope Surface Runoff and Macropore Flow
2.1. Kinematic Wave Equation for Slope Surface Runoff
Slope surface runoff takes place when the rainfall intensity exceeds the infiltration capacity of the slope. The surface runoff can be modeled using the one-dimensional kinematic wave equation [18], which is expressed aswhere is the surface ponding depth (m); is time (s); is the flow rate per unit width (m2/s) computed as () with being Manning’s roughness coefficient of slope surface; is the distance between a certain point at slope surface and the toe of the slope (m); is the rainfall intensity (m/s); is the angle of the slope from the horizontal; and is the infiltration rate of slope soil (m/s). In the brackets are the units for the variables in (1). Readers interested in details of derivations of (1) and its application are directed to, for example, Walker and Skogerboe [19] and Huang and Yeh [20].
2.2. Dual-Permeability Model for Water Flow in Micropore and Macropore Domains
Dual-permeability models are widely used to simulate water flow in soils with both micro- and macropores (e.g., [21–23]). There are in general two types of dual-permeability models based on model assumptions: (1) water is mobile in macropores while it is immobile in micropores and (2) water is mobile in both micro- and macropores. In this study, dual-permeability model with the second assumption is adopted, which is expressed as [13]
Here, , , and are volumetric water content (dimensionless), unsaturated permeability coefficient (m/s), and pressure head (m) in macropores, respectively; , , and are volumetric water content (dimensionless), unsaturated permeability coefficient (m/s), and pressure head (m) in soil matrix (i.e., micropores), respectively; is depth taken to be positive downward (m); is ratio of volume of macropores to volume of the total soil (dimensionless); is water transfer rate between soil matrix and macropores (1/s).
Gerke and van Genuchten [16] proposed a first-order expression for the water transfer rate () assuming is proportional to the difference in pressure head between the macropore and soil matrix domains. The expression of is [16, 17]where is geometric shape factor of soil aggregates, that is, for cubic soil aggregates and for spherical soil aggregates; is adopted in the analyses that are presented later in this study. The parameter is half width of the soil matrix or thickness of the matrix mantle in the case of hollow cylindrical geometry (m); is a nondimensional scaling factor which can be typically taken as 0.4 [17]; and is effective permeability coefficient at the macropore-soil matrix interface, which is commonly defined as an arithmetic average relating to both and as follows [24, 25]:
3. Numerical Method to Solve the Slope Surface Runoff and Macropore Flow Models
It is generally difficult to derive analytical solutions for (1), (2a), and (2b) even under simple initial and boundary conditions. In this study, the finite difference method was adopted to numerically solve the equations. The following presents the discretization of the equations, the convergence criteria for numerical calculations, and the steps to carry out an integrative analysis. The numerical calculations were performed using Matlab™.
3.1. Discretization of Equations
The kinematic wave equation of slope surface runoff (1) is discretized based on time and space using the Lax-Friedrichs method [26]. The discretized equation can be written as where and are computation time and space steps, respectively. Note that in the discretization is equal to () as noted earlier.
Similarly, the macropore flow equations ((2a) and (2b)) are discretized using the Alternating Direction Implicit Method [27], which is expressed as [28]
3.2. Convergence Criteria for Numerical Computation
The philosophy of selecting convergence criteria is that the computation outcomes should be sufficiently accurate while the computation should be numerically efficient. Previous studies [7] have shown that a satisfactory convergence criterion can be set as the difference in the computation outcomes of two successive computation steps is no larger than 0.1%. This convergence criterion is adopted here and mathematically, it is expressed aswhere is the th computation step.
3.3. Steps to Perform Integrative Analysis of Slope Surface Runoff and Macropore Flow
Analyses considering both slope surface runoff and macropore flow are conducted based on infiltration capability and water supply rate on the slope surface. The water supply rate depends on rainfall intensity and existing ponding depth. When the infiltration rate of the slope surface soil is larger than the water supply rate, all the water is absorbed by the soil. In this case water only flows in both macropores and soil matrix (no surface runoff) and (2a) and (2b) are employed for analysis. As the water supply rate exceeds the infiltration rate, slope surface runoff also takes place and (1), (2a), and (2b) are needed for analysis. Steps to perform integrative analyses are outlined as follows.
Step 1 (determine soil infiltration rate, , at slope surface). Since the slope soil consists of both domains of macropore and micropore, the total infiltration rate of the soil is the sum of infiltration rates of both domains. Based on the discretization of the macropore flow equation ((2a) and (2b)), the soil infiltration rate can be determined as [25]
Here, , , and are infiltration rate for soil matrix, macropore, and the total soil at the slope surface, respectively. Parameters and are permeability coefficients of soil matrix at Node 0 and Node 1, respectively. Node 0 is on the slope surface and Node 1 is the first node along depth (i.e., direction) as shown in Figure 1. Parameter is pressure head of soil matrix at Node 1. For the macropore domain, parameters of the same meaning are denoted as , , and . Parameter is available water from rain and existing surface ponding on the slope surface which is calculated as
At initial conditions, for example, at the beginning of a rainfall, there is no surface ponding and is equal to 0.
Step 2 (calculate water supply rate, ). The water supply rate () is the available water over a small time interval, which is mathematically expressed as
Step 3 (set boundary conditions). The boundary conditions are set by comparing the soil infiltration rate, , to the water supply rate, . If , then the actual infiltration is controlled by . In this case, water only flows in macropores and soil matrix and the boundary condition for analysis using (2a) and (2b) is where is the pressure head at the boundary node, that is, Node 0; is the pressure head at Node 1; is the permeability coefficient of soil matrix under saturated conditions; is the permeability coefficient of soil matrix at Node 1.
If , then the infiltration is . The boundary condition for (2a) and (2b) can be defined as
Step 4 (compute the volumetric water contents, and , given initial values of and ). and are soil water characteristic curves relating and to and , respectively.
If , compute at each node using the boundary condition of (11) and (2b) (discretized form (6b)).
If , then (1) is used to determine the ponding depth with boundary conditions described as (12). Equations (2a) and (2b) or equivalently the discretized forms (6a) and (6b) are adopted to calculate and .
Step 5 (check the satisfaction of convergence criteria of (7a) and (7b)). If both (7a) and (7b) are satisfied simultaneously, the analysis continues to next time step until the end of rainfall; otherwise, rerun the analysis from Step 4 by assuming different initial values of and .
The steps to carry out integrative analyses of slope surface runoff and macropore flow are summarized in the flow chart shown in Figure 2.
4. Evaluation of Performance of the Surface Runoff and Macropore Flow Model
4.1. Laboratory Tests
Laboratory tests were conducted to simulate the slope surface runoff and macropore flow process. Monitoring data from laboratory tests were then used to verify the proposed method of jointly using (1), (2a), and (2b) for analysis. Figure 3 shows the setting of laboratory tests. A typical silty clay sampled from a slope site in Wuyishan, China, was used as the soil matrix. A coarse sand was used to simulate the macropore domain. Both soils were placed in a plexiglass box which was 100 cm (length) × 50 cm (width) × 70 cm (depth) (Figure 3(a)). The slope angle was set as 15°. The variation of soil water contents was monitored and recorded using TDR (time-domain reflectometer) probes.
(a)
(b)
(c)
A rainfall simulator system similar to that used by Wu et al. [29] was built. The rainfall simulator system consisted of a pump, valves, PVC (Polyvinyl Chloride) pipes, pressure gauges, and raindrop generator as shown in Figure 3(b). To assess the sprinkling homogeneity of the rainfall simulator system, three graduated cylinders were placed at points A, B, and C at the bottom of the plexiglass box, as shown in Figure 3(c). The system was then initiated to simulate a rainfall with “nominal” intensity of 6.25 mm/h for one hour. The water amounts collected by the three graduated cylinders were recorded. This test was repeated for five times and the average values for the cylinders at points A, B, and C were 6.4, 6.7, and 6.0 mm/h, respectively, which were basically close to the nominal value of 6.25 mm/h. Hence, the performance of the rainfall simulator system was considered satisfactory from the perspective of sprinkling homogeneity. Test procedures to simulate the slope surface runoff and macropore flow process are described in detail in the following.
Step 1. The plexiglass box was sealed to be impervious for all sides. While at the bottom of the box, small holes with a diameter of 5 mm were drilled at a spacing of 10 cm in both horizontal directions. A layer of geotextile was then placed at the bottom. This was done to simulate a free-drainage boundary condition.
Step 2. The box was filled with silty clay layer by layer with each of 5 cm thick. The initial volumetric water content was about 15%. Each clay layer was compacted to about 85% degree of compaction. Here, the degree of compaction is the ratio of dry soil unit weight to the maximum dry soil unit weight. The maximum dry soil unit weight is typically determined through the standard proctor test [30]. A PVC pipe with a diameter of 2 cm was placed at the middle of the box (Figure 3(a)) when the filling height reached half of the box height. After the placement of the pipe, filling and compaction of the silty clay were continued. A total of six TDR probes were embedded in the silty clay to measure the water content.
Step 3. The PVC pipe was carefully pulled out after completion of the silty clay filling and compaction. Then, the hole was filled with coarse sand which was used to simulate the macropore domain.
Step 4. A rainfall intensity of 6.25 mm/h (equivalently 150 mm/d) was chosen since this is a typical rainfall intensity in typhoon-rainstorm prone area. The rainfall duration was taken as one hour (60 minutes) in the experiment.
4.2. Determination of Input Parameters
The saturated and unsaturated coefficients of permeability must be known before integrative analyses of surface runoff and macropore flow can be conducted. These parameters in this study were determined based on the soil water characteristic curve (SWCC) obtained through laboratory experiments.
The van Genuchten model [31] for description of SWCC curves (i.e., water content versus pressure head) was adopted in this paper, which is expressed aswhere and are residual and saturated volumetric water contents, respectively. Parameters and are empirical constants dependent on soil types. Parameter is pressure head as defined earlier in this paper. The filter paper method [31] was used to determine the and values and (13) was fitted to experiment results to determine the values for and . The outcomes are summarized in Table 1.
The saturated coefficients of permeability for the silty clay and coarse sand were measured using falling head and constant head permeability tests, respectively. The results are also shown in Table 1. Unsaturated coefficients of permeability for both soils can be described using the van Genuchten and Mualem model as [32, 33]where for soil matrix and for macropore; is the effective degree of saturation defined as
4.3. Comparisons between Experimental Results and Model Solutions
The water content variation with time (rainfall duration) was recorded by the six TDR probes (i.e., A1, A2, A3, B1, B2, and B3) at depths of 5, 15, and 25 cm from the slope surface. The data were collected at a sampling rate of one data point per minute. Figure 4 shows the measured water contents with time. The volumetric water contents at all depths increase with increasing rainfall duration. The initial in the micropore domain were about 15% at all depths (time = 0 minutes) but increased to about 40%, 35%, 30%, 28%, 27%, and 24% at points A1, A2, A3, B1, B2, and B3 at time of 60 minutes, respectively. The water content was higher at shallower depth than at greater depth.
(a)
(b)
For validation purpose, the established model ((1), (2a), and (2b)) was used to numerically compute the water contents at the same points and same depths with time using the input parameters given in Table 1. The steps to solve the model have been summarized in Figure 2. The predicted water contents using the established model were also plotted in Figure 4 for comparison. The trend in measured values is well captured (at least visually) by the predicted values at all depths. The good agreements were confirmed quantitatively using a statistical approach as shown in Figures 5 and 6.
(a)
(b)
(a)
(b)
Figure 5(a) shows the plots of measured values from laboratory tests versus paired predicted values at points A1, A2, and A3 using the analytical model established in the present paper. Similar plots are shown in Figure 5(b) for points B1, B2, and B3. All datasets in the figures follow closely the 1 : 1 correspondence line. The mean values of bias were close to 1.0 at all points and depths. Here, bias is defined as the ratio of measured to predicted value. The close-to-one mean values for the three bias datasets indicate that the predictions using the proposed model are satisfactorily accurate on average, regardless of depths. The respective coefficients of variation (COV) of bias values were basically between 3% and 7%. The COV of bias values were all small, implying small spreads in prediction accuracy.
The bias values were then plotted against the predicted values as shown in Figures 6(a) and 6(b). Visual inspections suggest that there are no statistical correlations between bias and predicted values. The observations were confirmed quantitatively using Spearman’s rank correlation test. Spearman’s rank correlation test is widely used to measure the monotonic relationship (both linear and nonlinear) between two datasets. The values of Spearman’s rank correlation tests were larger than 0.05. This means the rejection of the null hypotheses that the bias values and the predicted values are correlated at a level of significance of 5%. Hence, based on the mean of bias, COV of bias, and Spearman’s rank correlation test outcomes, the performance of the proposed model is judged to be satisfactory. However, it is also important to note that the high accuracy of the analytical model in this study is largely attributed to the fact that the measured values are obtained under well-controlled laboratory conditions. The accuracy of the model can be expected to be less when applied to simulate in situ slope surface runoff and macropore flow.
5. Parametric Study
The parametric study presented in this section was based on a slope with geometry shown in Figure 1. Input parameters and boundary conditions are the same as given in Table 1. Influences of presence of macropores, rainfall intensity, initial pore water pressure, and coefficient of permeability of soil matrix on the maximum depth of wetting front in soil and ponding depth on slope surface are examined. The baseline case is set as rainfall intensity = 12.5 mm/h (equivalently 300 mm/d), initial pore water pressure of −10 m (initial pressure head), and permeability coefficient of = 5.70 × 10−5 cm/s for soil matrix. In the following analyses, the water content is always for the soil matrix.
5.1. Effect of Presence of Macropores
The presence of macropores in soils greatly increases the infiltration depth of wetting front within the same time period compared to the case with soil matrix only, as shown in Figure 7. For example, the wetting front reaches a maximum depth of about 6 cm and 8 cm after 1.2 hours and 2.4 hours of rainwater infiltration, respectively, if there are no macropores in the soils. However, with the existence of macropores the maximum depth of wetting front would reach about 58 cm for 1.2 hours and 76 cm for 2.4 hours. The difference between these two cases is about one order of magnitude based on maximum depth of wetting front. This can be expected since the coefficient of permeability for macropores is three orders larger than that of the soil matrix (Table 1).
The ponding depth on slope surface is smaller if macropores are present in soil, as shown in Figure 8(a). However, the difference is not significant. This can be attributed to the occurrence of surface runoff on an inclined slope, which also explains the observation of larger ponding depth at lower part of the slope. The effect of presence of macropores can be expected to be larger on smoother slopes, that is, smaller slope angle.
(a)
(b)
Since the ponding depth varies along the slope surface, the boundary pressure in (11) also varies at each point “”. This results in different infiltration rates and thus depths of wetting front corresponding to different locations, as shown in Figure 8(b). Within the same time period, the closer to the toe of the slope (i.e., smaller ), the higher the ponding depth, and hence the deeper the wetting front. However, the difference in depths of wetting front is small due to the fact that the difference in boundary pressures caused by different ponding depths is small.
Figure 9 shows the accumulation of slope surface ponding at m and m with increasing duration of rainfall. At the early stage of the rainfall event, there is basically no surface ponding since all rainwater is absorbed into the soil. As the rainfall lasts, the soil near the slope surface becomes saturated and surface runoff takes place. It takes longer time for surface runoff to take place on a slope with macropores; however, the difference of time for initiation of surface runoff on slopes with and without macropores is not significant in the example used for parametric analyses in this study.
5.2. Sensitivity Analysis for Pressure Head (Depth of Wetting Front) in the Slope
The influences of rainfall intensity, initial pressure head, saturated permeability coefficient of soil matrix, and slope angle on the pressure head profile in the slope are examined in this section. Their influences on the surface ponding depth are presented in the next section.
Keeping other input parameters constant, the rainfall intensity was varied from 3.33 mm/h (80 mm/d) to 12.5 mm/h (300 mm/d) to examine its effect on computed maximum depth of wetting front and surface ponding. Figure 10(a) shows the maximum depths of wetting front corresponding to different rainfall intensities for slopes with soil matrix only and with both soil matrix and macropores. The effect of rainfall intensity on wetting front depth is practically negligible for slopes without macropores, that is, with soil matrix only. This is due to the fact that the three rainfall intensities used for analyses were all much higher than the coefficient of permeability of the soil matrix under saturated conditions. For slopes with macropores, larger maximum depth of wetting front was found for higher intensity of rainfall, for instance, about 50 cm for rainfall intensity of 3.33 mm/h but about 76 cm for rainfall intensity of 12.5 mm/h.
(a)
(b)
(c)
(d)
The initial conditions of pore water pressure expressed as pressure head was set at −0.5 m, −5 m, and −10 m. The initial pore water pressure of −0.5 m was intended to simulate slopes at a state close to fully saturated conditions, while the value of −10 m was used for slopes under dry conditions. The plots of pressure head in the depth direction are shown in Figure 10(b) for analyses using different initial pore water pressures. For slopes with macropores, soils within a depth of about 60 cm become saturated after 2.4 hours of rainfall, regardless of the initial pore water pressures.
Figure 10(c) shows the depths of wetting front with respect to three different coefficients of permeability assumed for the soil matrix, that is, = 1.2 × 10−5 cm/s, 5.7 × 10−5 cm/s, and 1.2 × 10−4 cm/s. The wetting front depth was greatly influenced by . As increases, the soil matrix is more permeable and the rainwater infiltrates into the soils more rapidly. The effect of is insignificant within the range under consideration for the case of slopes with soil matrix only.
Figure 10(d) shows the effects of slope angle on the infiltration of rainwater. The depth of wetting front decreases from about 80 cm to about 60 cm as the slope angle correspondingly increases from 5° to 25° based on a = 2.4 h rainfall period. This is because the surface runoff flows at higher speeds on steeper slopes, which results in less surface ponding depths and thus shallower depths of the wetting front.
5.3. Sensitivity Analysis for Surface Ponding Depth along the Slope
The surface ponding depth increases with increasing rainfall intensity as shown in Figure 11(a). The surface ponding depths appear to be roughly the same for slopes with and without macropores, although larger ponding depth is obtained for the latter case.
(a)
(b)
(c)
(d)
Figure 11(b) showed that smaller initial pore water pressures lead to larger surface ponding depths. At = 10 m, the difference in surface ponding depths is negligible, while at = 1 m, the difference reaches about 3 mm. In general, the surface ponding depths are similar regardless whether there are macropores in slopes or not.
The surface ponding depth was much less influenced by in this analysis, although smaller value resulted in larger ponding depth. No noticeable difference was found between slopes with and without macropores in this case, as shown in Figure 11(c).
Figure 11(d) shows the effects of slope angle on the surface runoff patterns. The steeper the slope (e.g., increasing from 5° to 25°), the less the surface ponding depth (accordingly decreasing from about 11 mm to about 4 mm), given all other conditions unchanged. This observation is consistent with that from Figure 10(d).
6. Summary and Conclusions
The paper presents integrative analyses of rainwater infiltration and surface runoff for slopes subjected to short-duration high-intensity rainfall events. Soils in slope were considered to be a two-domain system, including micropore domain (soil matrix) and macropore domain. A dual-permeability model proposed in the literature was adopted to describe the rainwater movement and exchange between the two domains. A kinematic wave model was employed to describe the slope surface runoff of rainwater. The two models were jointly used to carry out an integrative analysis of rainwater infiltration inside slope soils and runoff on slope surface based on boundary conditions on slope surface. Laboratory tests were then conducted to evaluate the ability of using the two models for integrative analysis. The main conclusions drawn from this study are as follows.
(1) The approach proposed in this study of using dual-permeability model coupled with kinematic wave model is applicable for description of rainwater movement inside slope soils and runoff on slope surface. The approach is satisfactorily accurate on average and the spread in prediction is very small from a statistical point of view based on comparisons between predicted outcomes and laboratory testing outcomes.
(2) The presence of macropores in slope would change the infiltration of rainwater dramatically based on the maximum depth of wetting front. With the existence of macropores, the maximum depth of wetting front can be one order of magnitude larger than that for slopes without macropores. This highlights the need to take macropores into account in analyses of rainwater infiltration and movement in slopes and in slope stability analyses under rainfall conditions. Also, for slopes without the presence of macropores, the great majority of rainwater would become surface runoff which increases the risk of flooding hazards, whereas for slopes with macropores, a significant amount of rainwater would infiltrate into the soils, which decreases the risk of flooding hazards but increases the risk slope stability failures in the meanwhile.
(3) The influences of rainfall intensity, initial pore water pressure, coefficient of permeability of soil matrix, and slope angle are significant on the movement of rainwater inside slope soils but are practically insignificant on depth of surface ponding on slope surface.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
The authors are grateful for financial support from the National Natural Science Foundation of China (nos. 41302215 and 41772297) and the Natural Sciences and Engineering Research Council of Canada (NSERC).