Abstract
Darcy-scale capillary pressure is traditionally assumed to be constant. By contrast, a considerable gap exists between the measured and equilibrium capillary pressures when the same moisture saturation is considered with a high flow rate, and this gap is called the dynamic effect on the capillary pressure. In this study, downward infiltration experiments of sand columns are performed to measure cumulative infiltration and to calculate the wetting front depth and wetting front velocity in sands with different grain sizes. We estimate the equilibrium capillary pressure head or suction head at the wetting front using both the classical Green–Ampt (GAM) and modified Green–Ampt (MGAM) models. The results show that the performance of MGAM in simulating downward infiltration is superior to that of GAM. Moreover, because GAM neglects the dynamic effect, it systematically underestimates the equilibrium suction head in our experiments. We also find that the model parameters and of MGAM are affected by the grain size of sands and porosity, and the dynamic effect of the capillary pressure increases with decreasing grain size and increasing porosity.
1. Introduction
Infiltration involves gas and liquid flows in porous media and occurs during precipitation or when liquid contaminants leak underground or onto the soil surface. This contributes to runoff generation, crop irrigation, and transport of nutrients and contaminants. Infiltration is a main process in subsurface hydrology and plays a crucial role in geohazards, such as landslides, flooding, and groundwater contamination. Richards’ equation is the most general model to address such flows with spatially and temporally variable saturation [1, 2]. Nevertheless, computational time is an issue of numerically solving Richards’ equation in a large-scale simulation. Therefore, the simple one-dimensional Green–Ampt model (GAM) [3] is an attractive alternative [4, 5]. Actually, the GAM has been widely incorporated in large-scale hydrological processes and erosion models [6] such as HEC-HMS [7], WEPP [8], SWAT [9], and ANSWERS [10] models. However, describing the transient behavior during early infiltration based on the GAM remains a challenging task.
The classical GAM for downward infiltration can be expressed as an ordinary differential equation (ODE) [3, 11]: where is the saturated water content, is the initial water content, is the saturated hydraulic conductivity, is the length of the porous medium with infiltrated liquid (the wetting front depth), is the wetting front velocity, is the equilibrium suction head, is capillary pressure, is the water density, is the gravitational acceleration, and is the ponding height. In the model, a wetting front moves downward as a sharp moving boundary during downward infiltration. The porous media above and below the wetting front are assumed to have water saturated with and unsaturated with , respectively. In addition, the air pressure in the porous medium is assumed to be constant and, therefore, the viscous pressure drop due to the air movement is neglected. The capillary pressure across the moving wetting front is commonly assumed to be constant and determined from the soil water retention curve [12, 13].
Nevertheless, during the past few decades, studies conducted under nonequilibrium conditions have indicated that the soil water retention curve may depend on the dynamics of water flow [14, 15]. Water content measured under transient flow conditions has been shown to be significantly different from that measured under static and steady-state conditions [16–21]. Considering the soil water retention curve based on thermodynamics, Hassanizadeh and Gray [22], postulated the existence of a dynamic component in the unsaturated water flow. This dynamic component depends on both flow dynamics and the process of either drainage or imbibition. In addition, studies conducting column experiments have shown that during downward infiltration, the water pressure head at the gas-water or liquid-liquid interfaces changes as the flow velocity changes [11, 23–28]. The aforementioned studies suggest that the capillary pressure under dynamic conditions can be less than that under static conditions during infiltration. Therefore, infiltration can be better described by a GAM with a nonequilibrium suction head .
The modified Green–Ampt model (MGAM) was first developed in Hsu and Hilpert [29]. They verified the model by the experimental data from upward infiltration [30], and downward infiltration [25]. Hsu and Hilpert [29] showed that the MGAM is better than the classical GAM at describing both transient capillary rise and downward infiltration in dry sands. For downward infiltration, the MGAM can also be expressed as an ordinary differential equation [11]: where is the interfacial tension, is the grain size, is the lumped parameter related to a nondimensional function (ε) of porosity ϵ and , and are model parameters, and is the dynamic viscosity of water. The additional term is added into GAM to consider the dynamic effect of capillary pressure.
The physical concept of the MGAM and the additional parameters and are based on the theory of the pore-scale dynamic contact angle. The upscaling of dynamic contact angle to the dynamic effect of capillary pressure has been derived and discussed in Hsu and Hilpert [29]. Based on the studies of dynamic contact angle, is in the range of 1/2 to 1/3, but the value of highly relies on the solid surface structure [31–34].
Hsu et al. and Pellichero et al. [11, 26] showed that the MGAM approach can avoid the initial unphysical infinite rate of infiltration and better describe the experimental results than can the classical GAM for downward infiltration at an early stage in both dry and prewetted sand columns with varied initial water contents. Indications are that the grain size and pore size distributions affect the dynamic effect on the soil water retention curve [20, 35–37]. DiCarlo [36] showed that when a saturation overshoot occurs during the downward infiltration, the length of the gravitational fingering tips varies with the grain size and distribution. Hilpert [38] proposed a theory based on velocity-dependent capillary pressure that correctly predicts the degree to which the saturation overshoot depends on both the liquid contents and grain size. Hsu and Hilpert [29] pointed out that the model parameters and β in MGAM might be related to the porosity and grain size. Based on (2), the dynamic effect on the suction head should be inversely proportional to the grain size. However, the effect of the grain size and pore size distribution of the sand column on the nonequilibrium suction head as well as the model parameters and in (2) of the MGAM have not yet been systematically investigated. Therefore, we performed a series of infiltration experiment to systematically investigate the effect of the grain size of the sand column on the nonequilibrium suction head as well as the model parameters and .
In this study, the cumulative infiltration depths were measured during a series of downward infiltration experiments in dry sand columns with different grain sizes, subject to different constant ponding heights. The remainder of this paper is organized as follows. In Section 2, we describe our experimental method and preparation of the porous materials. In Section 3, we describe the results of the experiments, and in Section 4, we compare the predictions from classical GAM and MGAM as well as the values of fitting parameters. We also discuss the effects of the grain size on the dynamic capillary pressure and the fitting of the equilibrium suction head.
2. Experimental Setup
2.1. Experimental Material Properties and Grain Size Measurement
In this study, sands (glass beads) with four grain sizes (labeled B35, B60, B150, and B320) were used for our infiltration experiments. The sands had the same shape, with a particle density of 2.5 g/cm3. Because glass beads after conditioning hardly aggregate, we could exclude the influence of the aggregates on the observed dynamic capillary pressure. Moreover, the grain size distribution (GSD) of the sands was determined through sieve analysis (ASTM C136) of eight different mesh screens.
2.2. Downward Infiltration Experiment
The infiltration experiments were conducted using glass filter columns (inner diameter = 2.6 cm, depth = 60 and 30 cm, cross-section area = 5.3 cm2). Based on Wang et al. [39], the sizes of the fingering flow for downward infiltration in sands are mostly ranged from 3 to 15 cm. The fingering flow with a width of less than 3 cm was only observed by Glass et al. [40]. In this study, our sand column with the inner diameter of 2.6 cm is small enough to avoid 2D fingering flows. All the sands were conditioned by being rinsed with deionized water to remove impurities and dried overnight in an oven at 100°C. We packed sands B35, B60, and B150 into the 60 cm column and sand B320 into the 30 cm column. To achieve uniform packing, the sands poured into the column were maintained at a constant 3 cm distance between the supply funnel and the top of the sand packing. In addition, a rubber hammer was used to tap the top of the sand in the column to obtain an even more homogeneous packing.
Figure 1 is a schematic representation of the experiment setup. A Mariotte’s bottle was connected through Tygon tubing and a valve to the sand column to maintain the hydraulic head. The bottle was placed on the top of an analytical balance (Sartorius TE1502S) to record the weight at 0.2 s intervals and automatically send the data to a computer. When the valve was opened, the water in the Mariotte’s bottle flowed into the column packed with dry sand, infiltrated the sand and reached the bottom of the column. The weight change measured by the analytical balance was used to calculate the cumulative infiltration and infiltration rate. All infiltration experiments in this study were conducted in triplicate.
2.3. Saturated Hydraulic Conductivity (Ks) Measurement
Saturated hydraulic conductivity (Ks) is an important parameter of the GAM, as it determines the flow rate of water in the saturated soil. Ks in this study was determined by the constant-head method. The measurement procedure was the same as that for the infiltration experiment previously described. Following the infiltration experiments, the sand in the column was at saturation, and the flow rate gradually reached equilibrium. We measured the equilibrium infiltration rate, and Ks was estimated based on Darcy’s law: where is the flow rate, is the length of the sand column, is the total cross-section area, and is the constant hydraulic head. The flow rate value was obtained by calculating the weight change of water per unit time.
3. Results
3.1. Grain Size Distribution of Experimental Materials
Figure 2 shows the grain size distributions of the sands. Sand B35 was the coarsest in grain size, ranging from 2 to 7 × 10−2 cm and was followed in order by B60 (distributed between 1 and 6 × 10−2 cm), B150 (ranging from 5 × 10−3 to 2.3 × 10−2 cm), and B320 (which had the finest grain size and was distributed between 1 × 10−4 and 2 × 10−2 cm).
The curves in Figure 2 are plotted based on a grain size distribution equation provided by Fredlund et al. [41]. We used the same equation to calculate the geometric mean grain size (d50). Table 1 shows the geometric mean grain sizes of sands B35 (4.51 × 10−2 cm), B60 (3.48 × 10−2 cm), B150 (1.38 × 10−2 cm), and B320 (1.3 × 10−3 cm).
The sands were uniform. Specifically, the uniformity coefficients of sands B35, B60, B150, and B320 were 1.33, 1.46, 1.55, and 3.85, respectively. In the discussion that follows, we used d50 to quantify the effect of grain size on the infiltration model and dynamic capillary pressure.
3.2. Results of Infiltration Experiments
Table 1 shows the measured porosity, bulk density, and saturated hydraulic conductivity. In each of the 36 infiltration replicate experiments, we controlled the bulk density and porosity of the soil column. The average bulk density of sand B35 column was 1.36 g/cm3 and the average porosity was 0.46 cm3/cm3. The bulk density of sand B60 was 1.42–1.43 g/cm3, B150 was 1.38–1.39 g/cm3, and B320 was 1.33–1.36 g/cm3. The porosity of sands B60, B150, and B320 was 0.43, 0.45, and 0.47 cm3/cm3, respectively.
Under the assumption of GAM, the sand has a water content behind the wetting front and an initial water content ahead of the front. Mobile water content (Δθ) indicates the moisture change at the wetting front and can be estimated from the total cumulative infiltration and total volume of the sand column, that is, [11]. In our study, the Δθ of sand B35 column was obtained from 0.42 to 0.44 cm3/cm3, B60 from 0.39 to 0.41 cm3/cm3, B150 from 0.40 to 0.42 cm3/cm3, and B320 from 0.39 to 0.44 cm3/cm3.
The results of the saturated hydraulic conductivity from additional experiments show that B35 had the largest value, followed by B60, B150, and B320, at 1.84–1.97 × 10−1, 6.42–6.53 × 10−2, 1.16–1.42 × 10−2, and 0.8–1.1 × 10−3 cm/s, respectively.
We used the analytical balance to measure the mass (mw) of water infiltrating the sand column. In addition, the cumulative infiltration (F) was calculated from the mass (mw), water density (ρ), and soil column cross-sectional area (A), such that . We plotted the relationship between and time according to different water ponding heights and soil sizes. Figure 3 shows versus in our experiments. Each plotted point represents the average result of three replicate experiments. However, one of the infiltration experiments of soil B320 column at a water ponding height of 20 cm failed because of water leakage, a fact that was not mentioned in our subsequent discussion. Figures 3(a)–3(d) clearly show that increased as the water ponding height increased.
3.3. Wetting Front Depth and Velocity versus Water Ponding Height and Grain Size
Wetting front depth can be calculated from and , such that . Figure 4 shows versus in our experiments. In general, the monotonously increased with . In addition, water ponding height cm had the largest value at the same time, followed by cm and cm. For different grain sizes, the wetting front of B35 with the largest grain size and fastest infiltration rate spent 55–80.8 s to arrive cm, B60 spent 125–217.6 s, and B150 spent 597.2–980 s. Finally, the wetting front of B320 with the smallest grain size and slowest infiltration rate spent 2040–2430 s to arrive cm.
We divided by to obtain the wetting front velocity and plotted the relationship of and in Figure 5. To compare the wetting front velocities in different experiments with , we normalized the X axis through divided by the infiltration end time . Figure 5 shows versus standardized time . We can see that was the fastest at the beginning of infiltration, gradually decreased with , and finally approached an equilibrium value. Regardless of the water ponding height, increased as the grain size increased.
The main purpose of this study is to examine the relationship between flow velocity and capillary pressure during transient infiltration. The combination of ponding depths and grain sizes can lead to a range of flow rates and front velocities for examining our hypothesis. In addition, our experimental results shown in Figures 4 and 5 confirm that our experimental design in accord with the basic physical assumptions of infiltration theory and the GAM. It implies that the fingering flow did not occur during our experiments.
3.4. Infiltration Simulation of GAM and MGAM
In this study, experimental data were modeled by GAM and MGAM. We used a MATLAB implementation of ODE solvers to provide a numerical solution of GAM (1) and MGAM (2) in downward infiltration conditions. GAM required only a parameter for simulation, whereas other parameters and Δθ were inputted as fixed values (listed in Table 1). Parameters used by MGAM included , where and were inputted as fixed values (listed in Table 1). Here, grain size was represented by the geometric mean grain size (d50).
The fitting parameters were inversely determined from the measured data using a MATLAB nonlinear fitting function. The nonlinear fitting function uses a Levenberg-Marquardt nonlinear least squares algorithm to solve bound-constrained optimization problems [42]. MATLAB nonlinear fitting function must specify the initial value of . Based on the studies of Hsu et al. [11] and Pellichero et al. [26], we set and as the initial values for model fitting. The initial value of of sands B35, B60, B150, and B320 was , 3.48 × 10−2, 1.38 × 10−2, and 1.3 × 10−3 cm, respectively. The initial value of was obtained by where is interfacial tension (in which the surface tension of water and air is ), is gravity acceleration (), is water density (), is contact angle, and is effective pore radius. In this study, the contact angle was ; this was the static contact angle of quartz and water suggested by Friedman [43]. We used the empirical formula of Glover and Walker [44] to convert grain size into effective pore size : where is porosity. The aforementioned initial values were substituted in (1) and (2) to simulate GAM and MGAM for downward infiltration, respectively.
The results of infiltration simulation of GAM and MGAM are shown in Figure 6, and the model fitting parameters are listed in Table 2. In Figure 6, the black circles are the wetting front depths, which were calculated from the measured values of . The remaining four curves show the simulation results of GAM and MGAM. The orange-dotted curve (GAM) was plotted by GAM with the calculated based on (4). The yellow-dashed curve (MGAM) was simulated by MGAM with , , (which was calculated by (4)), and . The purple-dashed curve (optimized GAM) used the nonlinear fitting method to simulate downward infiltration by GAM with as a fitting parameter. Finally, the green-solid curve (optimized MGAM) was represented by MGAM with , , , and as fitting parameters. Both the purple-dashed and green-solid curves with fitting parameters showed better simulation results; these curves could approach the black circles of the wetting front depth. The orange-dotted curve with a fixed was farthest from the black circles, followed by the yellow-dashed curve with fixed and in the middle. This indicates that the performance of MGAM with multiparameters for simulating downward infiltration was superior to that of the classical GAM.
Regardless of the grain size of the sand, the optimized GAM and optimized MGAM with the fitting parameters had good simulation results. For GAM and MGAM with no fitting parameters, the smaller the grain size, the more the simulation results deviated from the black circles, particularly the sand B320.
Table 2 shows that two of the fitting suction heads of sand B35 based on GAM were negative. Because the grain size of the B35 was larger than that of other sands, the equilibrium was relatively small or even close to zero (Figure 7 and Table 2). Nevertheless, to have a negative suction head for a sand column packed with glass beads is still abnormal. It implies that the was underestimated by using classical GAM during a downward infiltration with high wetting front velocity. According to Chow et al. [12], the is distributed from 0.97 cm (0.0097 m) to 156.5 cm (1.565 m) from sand to clay. From the grain size distribution in Figure 2, the texture class of B320 can be classified as silt loam. Its value ranges from 2.92 cm (0.0292 m) to 95.39 cm (0.9539 m). However, theoretical values of sf based on (4) are over 2 m. The value for sand B320 is larger than that from Chowet al. [12], but it should still be in the range of silt loam or soils slightly smaller than silt loam for practical problems.
The performances of GAM and MGAM were evaluated using root mean square error (RMSE). For each simulation, the difference between the measured wetting front depth (lf(m)) and the calculated wetting front depth (lf(p)) was expressed in RMSE as: where N is the number of observed samples. By statistically verifying RMSE, we could evaluate the reliability of the model simulation; when the RMSE was small, the model was more reliable. The results of RMSE for infiltration simulation of GAM and MGAM are listed in Table 2. In general, the RMSEs of most MGAMs (0.001–0.008) in this study were less than those of GAMs (0.002–0.011), indicating that the simulation results of MGAM were superior to those of GAM. Except for two sets of simulation results for B60 and B320, the RMSE of GAM was smaller than that of MGAM. Fortunately, the RMSE of these two sets of results were also very small, which meant that the performance of MGAM was also good.
4. Discussion
4.1. Relationship between Wetting Front Suction Head and Grain Size
The relationship between the suction head and grain size can be linked by using (4) and (5). Figure 7 shows the relationship between the grain sizes of the sand columns and the equilibrium suction heads estimated by GAM and MGAM. The solid line in Figure 7 shows the theoretical relationship based on (4) and (5) between the grain sizes of the sand columns and the equilibrium suction heads with the assumed effective contact angle between air, water, and the glass bead of 30°. Using the same equations, we determined the effective pore radius and contact angles by fitting the effective contact angle to the grain size and the determined equilibrium suction heads. The dashed lines in Figure 7 show the fitted relationship between grain sizes and the equilibrium suction heads estimated by GAM and MGAM with the fitted effective contact angle of 84° and 63°, respectively.
The equilibrium suction heads estimated by the GAM were systematically lower than those estimated by the MGAM. Consequently, the effective contact angle estimated by the GAM was also systematically greater than that by the MGAM. The overestimated effective contact angle from GAM was because we neglected the dynamic effect when estimating the suction head. The effective contact angle estimated by GAM was the average dynamic contact angle during the infiltration process. Nevertheless, the contact angle (63°) estimated by MGAM was still greater than that (30°) measured in previous studies, even when the dynamic effect was considered. One explanation for this is related to the contact angle hysteresis: the contact angle no longer corresponds to the static equilibrium angle but is larger when the liquid is advancing and smaller when the liquid is receding [45]. In general, during the infiltration, the water advances the air in the sand. Therefore, the contact angle during the infiltration should be in the form of an advancing contact angle, which is commonly greater than the equilibrium contact angle.
4.2. MGAM Fitting Parameters for Dynamic Capillary Pressure
Hsu and Hilpert [29] pointed out that the parameters and in the MGAM are related to the dynamic effect on the capillary pressure, and their values might not be the same for different initial water contents, porosities, and grain sizes. Figure 8 shows two trends in which the fitted is proportional to the porosity and inversely proportional to the grain size. From the empirical relationship proposed by Stauffer [20] relating the dynamic coefficient to measurable porous media and fluid properties, the equation suggests that the dynamic effects would be greater for fine-grained soil with high air entry pressure and high porosity. Camps-Roach eat al. [35] observed that the dynamic coefficient statistically increases as the grain size decreases. Therefore, our fitted values of the parameters were consistent with those of the aforementioned studies.
(a)
(b)
Limited numbers of Darcy-scale studies focus on determining the values in sand. Weitz et al. [28] yielded for water displacing decane in a glass bead column. Hsu and Hilpert [29] tested a range of values from 1/5 to 1 based on the studies of pore-scale dynamic contact angles on flat surfaces and in circular capillary tubes, which can be described by [33, 34, 46]. The best value was approximately 0.31 for the experiments performed by Geiger and Durnford [25]. Pellichero et al. [26] showed that resulted in the best fit in downward infiltration experiments with the same sand but varying ponding heights. The same value also worked well in similar experiments but with various initial water contents [11]. In this study, was set as one of the fitting parameters. Figure 9 shows that the fitted was approximately 0.3, and the error bars overlapped for most grain sizes except for those of B320, which had the highest value (up to 0.34) with the smallest grain size. Nevertheless, the value of the fitted in this study agreed with the values derived or determined from the aforementioned Darcy and pore-scale studies. We confirmed again that 0.3 is a universal value for for Darcy’s scale sand experiments.
5. Conclusions
We performed a series of downward infiltration experiments and compared the predictions derived from classical GAM and MGAM as well as the values of fitting parameters. In this study, cumulative infiltration monotonously increased with time, and a greater water ponding height showed more cumulative infiltration. This kind of result was also reflected in the wetting front depth. Furthermore, the sand column with smaller grain size had a lower infiltration rate and wetting front velocity, and a greater amount of time was spent on infiltration. Regarding the model simulation, the RMSEs of MGAM were less than those of GAM, indicating that the performance of MGAM with multiparameters for simulating downward infiltration was superior to that of the classical GAM.
Because the classical GAM ignores the dynamic effect on estimating the suction head, the effective contact angle estimated by GAM was systematically greater than that by MGAM, meaning that the equilibrium suction from GAM was systematically less than that from MGAM. Moreover, the contact angle is typically in the form of an advancing one during infiltration. Therefore, the equilibrium contact angle estimated by MGAM was greater than the intrinsic contact angle. In addition, the fitted was proportional to the porosity and inversely proportional to the grain size in our study, indicating that the dynamic effect of capillary pressure increased with decreasing grain size and increasing porosity. Furthermore, the fitted was approximately 0.3, which is consistent with the values obtained from most of Darcy’s and pore-scale studies. This study confirmed that 0.3 can be the universal value of for Darcy-scale sand experiments.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this paper.
Acknowledgments
The authors thank the editors and anonymous referees for their thoughtful comments and suggestions. The authors also want to thank the Ministry of Science and Technology (MOST), Taiwan, for financially supporting this research under Grant MOST 106-2628-M-002-009-MY3 and MOST 106-2923-E-009-001-MY3.