Effects of the Grain Size on Dynamic Capillary Pressure and the Modified Green-Ampt Model for Infiltration

(1)

Research Article

Effects of the Grain Size on Dynamic Capillary Pressure and the

Modified Green–Ampt Model for Infiltration

Yi-Zhih Tsai ,

1

Yu-Tung Liu,

2

Yung-Li Wang,

1

Liang-Cheng Chang,

3

and Shao-Yiu Hsu

1

1_{Department of Bioenvironmental Systems Engineering, National Taiwan University, Taipei 10617, Taiwan} 2_{Department of Water Resources Engineering and Conservation, Feng Chia University, Taichung 40724, Taiwan} 3_{Department of Civil Engineering, National Chiao Tung University, Hsinchu 30050, Taiwan}

Correspondence should be addressed to Shao-Yiu Hsu; syhsu@ntu.edu.tw

Received 20 April 2018; Accepted 12 July 2018; Published 26 August 2018 Academic Editor: Ching Hung

Copyright © 2018 Yi-Zhih Tsai 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. 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 alsofind 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

Inﬁltration involves gas and liquid ﬂows in porous media and occurs during precipitation or when liquid contaminants leak underground or onto the soil surface. This contributes

to runoﬀ generation, crop irrigation, and transport of

nutri-ents and contaminants. Inﬁltration is a main process in

sub-surface hydrology and plays a crucial role in geohazards, such

as landslides, ﬂooding, and groundwater contamination.

Richards’ equation is the most general model to address such ﬂows with spatially and temporally variable saturation [1, 2]. Nevertheless, computational time is an issue of numerically solving Richards’ equation in a large-scale simulation. There-fore, 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 hydrolog-ical 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

inﬁltration based on the GAM remains a challenging task.

The classical GAM for downward inﬁltration can be

expressed as an ordinary diﬀerential equation (ODE) [3, 11]:

θs− θi

K_s dl_f

dt lf = hw+ sf+ lf, 1

whereθsis the saturated water content,θiis the initial water

content,K_sis the saturated hydraulic conductivity,l_f is the

length of the porous medium with inﬁltrated liquid (the

wet-ting front depth), νf= dlf/dt is the wetting front velocity,

sf = Pc/ρg is the equilibrium suction head, Pc is capillary

pressure,ρ is the water density, g is the gravitational

acceler-ation, andh_wis the ponding height. In the model, a wetting

front moves downward as a sharp moving boundary during

downward inﬁltration. The porous media above and below

the wetting front are assumed to have water saturated with

θsand unsaturated withθi, 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

(2)

movement is neglected. The capillary pressurePc 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 con-ducted under nonequilibrium conditions have indicated that the soil water retention curve may depend on the dynamics

of waterﬂow [14, 15]. Water content measured under

tran-sient ﬂow conditions has been shown to be signiﬁcantly

diﬀerent 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 ﬂow. This dynamic component

depends on both ﬂow dynamics and the process of either

drainage or imbibition. In addition, studies conducting

col-umn experiments have shown that during downward in

ﬁltra-tion, the water pressure head at the gas-water or liquid-liquid

interfaces changes as theﬂow velocity changes [11, 23–28].

The aforementioned studies suggest that the capillary pres-sure under dynamic conditions can be less than that under static conditions during inﬁltration. Therefore, inﬁltration can be better described by a GAM with a nonequilibrium

suction headsf.

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 inﬁltration [25]. Hsu and Hilpert [29] showed

that the MGAM is better than the classical GAM at

describ-ing both transient capillary rise and downward inﬁltration in

dry sands. For downward inﬁltration, the MGAM can also be

expressed as an ordinary diﬀerential equation [11]:

θs− θi K dl_f dtlf=hw+sf− γ dρgα η γ dl_f dt β +lf, 2

whereγ is the interfacial tension, d is the grain size, α = αε

ϵ, θi is the lumped parameter related to a nondimensional

function (ε) of porosity ϵ and θi,α and β are model

parame-ters, andη is the dynamic viscosity of water. The additional

term γ/dρg α η/γ dlf/dt βis added into GAM to

con-sider the dynamic eﬀect 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 eﬀect 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 aﬀect the dynamic eﬀect

on the soil water retention curve [20, 35–37]. DiCarlo [36]

showed that when a saturation overshoot occurs during the

downward inﬁltration, the length of the gravitational ﬁnger-ing 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 satu-ration 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

poros-ity and grain size. Based on (2), the dynamic eﬀect on the suc-tion head should be inversely proporsuc-tional to the grain size. However, the eﬀect of the grain size and pore size distribu-tion of the sand column on the nonequilibrium sucdistribu-tion 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 inﬁltration experiment to

systemati-cally investigate the eﬀect 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 inﬁltration depths were

measured during a series of downward inﬁltration

experi-ments in dry sand columns with diﬀerent grain sizes, subject

to diﬀerent 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 mate-rials. In Section 3, we describe the results of the experi-ments, and in Section 4, we compare the predictions from

classical GAM and MGAM as well as the values of ﬁtting

parameters. We also discuss the eﬀects of the grain size on

the dynamic capillary pressure and theﬁtting of the

equilib-rium 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 inﬁltration 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 inﬂuence of the aggregates on the observed dynamic cap-illary pressure. Moreover, the grain size distribution (GSD) of the sands was determined through sieve analysis (ASTM

C136) of eight diﬀerent mesh screens.

2.2. Downward Inﬁltration Experiment. The inﬁltration

experiments were conducted using glassﬁlter 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

fin-gering flow for downward infiltration in sands are mostly

ranged from 3 to 15 cm. Theﬁngering ﬂow 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ﬁngering ﬂows. 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

(3)

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

tub-ing 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 mea-sured 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 (K_s) is an important

param-eter of the GAM, as it dparam-etermines theﬂow rate of water in the

saturated soil. K_s in this study was determined by the

constant-head method. The measurement procedure was the same as that for the inﬁltration experiment previously described. Following the inﬁltration experiments, the sand

in the column was at saturation, and theﬂow rate gradually

reached equilibrium. We measured the equilibrium

inﬁltra-tion rate, and Kswas estimated based on Darcy’s law:

K_s=QL

Ah, 3

whereQ is theﬂow rate, L is the length of the sand column, A

is the total cross-section area, andh is the constant hydraulic

head. The ﬂow 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−2cm and was followed in order by B60 (distributed

between 1 and 6 × 10−2cm), B150 (ranging from 5 × 10−3to

2.3 × 10−2cm), and B320 (which had the ﬁnest grain size

and was distributed between 1 × 10−4and 2 × 10−2cm).

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 (d₅₀). Table 1 shows the geometric mean grain sizes of

sands B35 (4.51 × 10−2cm), B60 (3.48 × 10−2cm), B150

(1.38 × 10−2cm), and B320 (1.3 × 10−3cm).

The sands were uniform. Speciﬁcally, the uniformity coeﬃcients of sands B35, B60, B150, and B320 were 1.33, 1.46, 1.55, and 3.85, respectively. In the discussion that

fol-lows, we used d50to quantify the eﬀect of grain size on the

inﬁltration model and dynamic capillary pressure.

3.2. Results of Inﬁltration Experiments. Table 1 shows the

measured porosity, bulk density, and saturated hydraulic

conductivity. In each of the 36 inﬁltration replicate

experi-ments, we controlled the bulk density and porosity of the soil column. The average bulk density of sand B35 column was

1.36 g/cm3and 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

poros-ity 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 θ = θsbehind the wetting front and an initial water

content θi ahead of the front. Mobile water content (Δθ)

indicates the moisture change θs− θi at the wetting front

and can be estimated from the total cumulative inﬁltration

and total volume of the sand column, that is, Δθ = Ftotal/

V_total [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−3cm/s, respectively. Mariotte_′s bottle Analytical balance Computer Valve Column Wetting front

Figure 1: Schematic of the experimental setup for the downward inﬁltration experiments. 100 90 80 70 60 50 40 30 20 10 0 C u m u la ti ve p er cen ta ge b y w eig h t passin g (%) 10−1 ₁₀−2 ₁₀−3 ₁₀−4 Particle diameter (cm) B150 B320 B35 B60

Figure 2: Grain size distribution curves of sands B35, B60, B150, and B320.

(4)

We used the analytical balance to measure the mass

(mw) of water inﬁltrating the sand column. In addition, the

cumulative inﬁltration (F) was calculated from the mass

(mw), water density (ρ), and soil column cross-sectional area

(A), such that F = mw/ρA. We plotted the relationship

betweenF and time t according to diﬀerent water ponding

heights and soil sizes. Figure 3 shows F versus t in our

experiments. Each plotted point represents the average result of three replicate experiments. However, one of the inﬁltration 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 F increased as the water

ponding height increased.

Table 1: Geometric mean grain size, porosity, bulk density, and saturated hydraulic conductivity for column experiments with diﬀerent ponding heights.

Materials Geometric mean grain size, d₅₀× 10−4(cm) Water ponding height, h_w(cm) Bulk density, ρb(g·cm−3) Porosity,

ϵ (cm3_·cm−3₎ Mobile water content,_{Δθ (cm}3_·cm−3₎

Saturated hydraulic conductivity, K_s× 10−2(cm·s−1) B35 10 1.36 ± 0.006 0.46 ± 0.003 0.44 ± 0.011 19.68 ± 0.683 451 20 1.36 ± 0.002 0.46 ± 0.001 0.42 ± 0.006 18.35 ± 0.132 40 1.36 ± 0.017 0.46 ± 0.007 0.43 ± 0.006 18.37 ± 1.591 B60 10 1.43 ± 0.003 0.43 ± 0.001 0.41 ± 0.014 6.53 ± 0.202 348 20 1.42 ± 0.005 0.43 ± 0.002 0.39 ± 0.003 6.42 ± 0.079 40 1.42 ± 0.008 0.43 ± 0.003 0.39 ± 0.003 6.50 ± 0.393 B150 10 1.38 ± 0.004 0.45 ± 0.002 0.40 ± 0.001 1.21 ± 0.024 138 20 1.38 ± 0.007 0.45 ± 0.003 0.40 ± 0.005 1.16 ± 0.055 40 1.39 ± 0.000 0.45 ± 0.000 0.42 ± 0.021 1.42 ± 0.391 B320 10 1.36 ± 0.004 0.46 ± 0.002 0.39 ± 0.007 0.08 ± 0.003 13 20 1.34 ± 0.013 0.47 ± 0.005 0.44 ± 0.040 0.11 ± 0.026 40 1.33 ± 0.006 0.47 ± 0.002 0.42 ± 0.036 0.10 ± 0.028 25 20 15 10 5 0 0 20 40 60 80 Time (s) C u m u la ti ve infil tra tio n (cm) 25 20 15 10 5 0 0 50 100 150 200 250 Time (s) (a) (b) (c) (d) C u m u la ti ve infil tra tio n (cm) Time (s) 0 200 400 600 800 1000 C u m u la ti ve infil tra tio n (cm) 30 25 20 15 10 5 0 hw = 10 cm hw = 20 cm hw = 40 cm Time (s) 0 500 1000 1500 2000 2500 C u m u la ti ve infil tra tio n (cm) 14 12 10 8 6 4 2 0

(5)

3.3. Wetting Front Depth and Velocity versus Water Ponding

Height and Grain Size. Wetting front depth l_f can be

cal-culated fromF andΔθ, such that lf=F/Δθ. Figure 4 shows

lf versus t in our experiments. In general, the lf

monoto-nously increased with t. In addition, water ponding height

hw= 40 cm had the largest lf value at the same time,

followed byhw= 20 cm andhw= 10 cm. For diﬀerent grain

sizes, the wetting front of B35 with the largest grain size and

fastest inﬁltration rate spent 55–80.8 s to arrive lf= 60 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

slow-est inﬁltration rate spent 2040–2430 s to arrive lf = 30 cm.

We divideddl_f bydt to obtain the wetting front velocity

νf and plotted the relationship ofνf andt in Figure 5. To

compare the wetting front velocities in diﬀerent experiments

witht, we normalized the X axis through t divided by the

inﬁltration end time tie . Figure 5 showsν_fversus

standard-ized time t/tie . We can see that νf was the fastest at the

beginning of inﬁltration, gradually decreased with t, and ﬁnally approached an equilibrium value. Regardless of the

water ponding height,νfincreased as the grain size increased.

The main purpose of this study is to examine the

relationship between ﬂow velocity and capillary pressure

during transient inﬁltration. The combination of ponding

depths and grain sizes can lead to a range of ﬂow rates

and front velocities for examining our hypothesis. In

addition, our experimental results shown in Figures 4 and

5 conﬁrm that our experimental design in accord with the

basic physical assumptions of inﬁltration theory and the

GAM. It implies that theﬁngering ﬂow did not occur

dur-ing our experiments.

3.4. Inﬁltration 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 inﬁltration conditions. GAM

required only a parameter sf for simulation, whereas other

parametersKs andΔθ were inputted as ﬁxed values (listed

in Table 1). Parameters used by MGAM included sf,α, β,

d, Ks, andΔθ, where KsandΔθ were inputted as ﬁxed values

(listed in Table 1). Here, grain size d was represented by

the geometric mean grain size (d50).

The ﬁtting parameters were inversely determined from

the measured data using a MATLAB nonlinear ﬁtting

function. The nonlinear ﬁtting function uses a

Levenberg-Marquardt nonlinear least squares algorithm to solve bound-constrained optimization problems [42]. MATLAB

nonlinear ﬁtting function must specify the initial value of

s_f,α, β, and d. Based on the studies of Hsu et al. [11] and

Pellichero et al. [26], we setα = 100 and β = 0 3 as the initial

values for modelﬁtting. The initial value of d of sands B35,

60 50 40 30 20 10 0 W et tin g f ro n t dep th (cm) 0 20 40 60 80 Time (s) 60 50 40 30 20 10 0 W et tin g f ro n t dep th (cm) 0 50 100 150 200 Time (s) (a) (b) (c) (d) 60 50 40 30 20 10 0 W et tin g f ro n t dep th (cm) 0 200 400 600 800 Time (s) hw = 10 cm hw = 20 cm hw = 40 cm 30 25 20 15 10 5 0 W et tin g f ro n t dep th (cm) 0 500 1000 1500 2000 Time (s)

Figure 4: Wetting front depth of sand versus time: (a) B35, (b) B60, (c) B150, and (d) B320. The column length of sands B35, B60, and B150 was 60 cm and of sand B320 was 30 cm.

(6)

B60, B150, and B320 was d50= 4 51 × 10−2, 3.48 × 10−2,

1.38 × 10−2, and 1.3 × 10−3cm, respectively. The initial value

ofsf was obtained by

sf =

2γ cos φ

rρg , 4

whereγ is interfacial tension (in which the surface tension of

water and air isγ = 7 2 × 10−2N/m),g is gravity acceleration

(g = 9 81 m/s2_),_{ρ is water density (ρ = 1000 kg/m}3_),_{φ is}

con-tact angle, and r is eﬀective pore radius. In this study, the

contact angle was φ = 30°; 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 d into eﬀective pore size r :

r = d

2 8/3 × 1 52/8ϵ2×1 5, 5

whereϵ is porosity. The aforementioned initial values were

substituted in (1) and (2) to simulate GAM and MGAM for

downward inﬁltration, respectively.

The results of inﬁltration simulation of GAM and

MGAM are shown in Figure 6, and the modelﬁtting

param-eters are listed in Table 2. In Figure 6, the black circles are the wetting front depths, which were calculated from the

measured values of mw. The remaining four curves show

the simulation results of GAM and MGAM. The orange-dotted curve (GAM) was plotted by GAM with the

calcu-lated sf based on (4). The yellow-dashed curve (MGAM)

was simulated by MGAM with α = 100, β = 0 3, sf (which

was calculated by (4)), andd = d50. The purple-dashed curve

(optimized GAM) used the nonlinearﬁtting method to

sim-ulate downward inﬁltration by GAM with sf as a ﬁtting

parameter. Finally, the green-solid curve (optimized MGAM)

was represented by MGAM with α, β, sf, andd as ﬁtting

parameters. Both the purple-dashed and green-solid curves

with ﬁtting parameters showed better simulation results;

these curves could approach the black circles of the wetting

front depth. The orange-dotted curve with aﬁxed sf was

far-thest from the black circles, followed by the yellow-dashed

curve withﬁxed α, β, sf , and d in the middle. This indicates

that the performance of MGAM with multiparameters for

simulating downward inﬁltration was superior to that of

the classical GAM.

Regardless of the grain size of the sand, the optimized

GAM and optimized MGAM with the ﬁtting parameters

had good simulation results. For GAM and MGAM with

no ﬁtting parameters, the smaller the grain size, the more

the simulation results deviated from the black circles, partic-ularly the sand B320.

Table 2 shows that two of theﬁtting suction heads sf of

sand B35 based on GAM were negative. Because the grain 4 vf (cm/s) 3 2 1 0 0 0.2 0.4 0.6 0.8 1 t/tie 4 vf (cm/s) 3 2 1 0 0 0.2 0.4 0.6 0.8 1 t/tie (a) (b) 6 vf (cm/s) 4 2 0 0 0.2 0.4 0.6 0.8 1 t/tie (c) B35 B60 B150 B320

Figure 5: Wetting front velocity of sands B35, B60, B150, and B320 under ponding height of: (a) h_w= 10 cm, (b)hw= 20 cm, and (c) h_w= 40 cm.

(7)

size of the B35 was larger than that of other sands, the

equi-libriumsf 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 sf was underestimated by using classical

GAM during a downward inﬁltration with high wetting front 60 50 40 30 20 10 0 W et tin g f ro n t dep th (cm) 0 10 20 30 40 50 60 70 80 Time (s) hw = 10 cm B35 60 50 40 30 20 10 0 W et tin g f ro n t dep th (cm) 0 10 20 30 40 50 60 70 Time (s) hw = 20 cm 60 50 40 30 20 10 0 W et tin g f ro n t dep th (cm) 0 10 20 30 40 50 Time (s) (a) (b) (c) (d) (e) (f) (g) (h) (i) (j) (k) (l) hw = 40 cm 60 50 40 30 20 10 0 W et tin g f ro n t dep th (cm) B60 0 50 100 150 200 Time (s) 60 50 40 30 20 10 0 W et tin g f ro n t dep th (cm) 0 20 40 60 80 100 120 140 160 Time (s) 60 50 40 30 20 10 0 W et tin g f ro n t dep th (cm) 0 20 40 60 80 100 120 Time (s) 60 50 40 30 20 10 0 W et tin g f ro n t dep th (cm) 0 200 400 600 800 Time (s) B150 60 50 40 30 20 10 0 W et tin g f ro n t dep th (cm) 0 100 200 300 400 500 600 700 800 Time (s) 60 50 40 30 20 10 0 W et tin g f ro n t dep th (cm) 0 100 200 300 400 500 600 Time (s) 0 500 1000 1500 2000 Time (s) Optimized GAM Optimized MGAM 30 25 20 15 10 5 0 W et tin g f ro n t dep th (cm) Measured data GAM MGAM B320 0 500 1000 1500 2000 Time (s) 30 25 20 15 10 5 0 W et tin g f ro n t dep th (cm) 0 500 1000 1500 2000 Time (s) 30 25 20 15 10 5 0 W et tin g f ro n t dep th (cm)

(8)

velocity. According to Chow et al. [12], thesf 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

tex-ture class of B320 can be classiﬁed as silt loam. Its sf 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 sf 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 diﬀerence between the measured wetting front depth

(lf(m)) and the calculated wetting front depth (lf(p)) was

expressed in RMSE as:

RMSE = ∑ N i=1 lf m − lf p 2 N , 6

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 inﬁltration 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 sf

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 the-oretical relationship based on (4) and (5) between the grain Table 2: Model parameters and RMSE of GAM and MGAM for downward inﬁltration simulation.

Materials hw(cm) GAM MGAM s_f(m) RMSE s_f (m) α β d (×10−4cm) RMSE B35 10 0.011 0.007 0.105 86.138 0.305 425 0.004 20 −0.009 0.011 0.094 93.316 0.297 456 0.008 40 −0.052 0.009 0.09 100.741 0.283 442 0.006 B60 10 0.053 0.002 0.112 85.391 0.339 361 0.003 20 0.025 0.005 0.133 98.828 0.3 351 0.002 40 0.015 0.007 0.135 93.037 0.29 361 0.004 B150 10 0.134 0.008 0.326 97.797 0.288 138 0.003 20 0.147 0.004 0.323 98.817 0.3 138 0.001 40 0.099 0.006 0.31 98.585 0.291 137 0.002 B320 10 0.528 0.003 1.063 95.888 0.335 16 0.003 20 0.467 0.002 1.007 97.124 0.335 17 0.004 40 0.279 0.009 1.496 109.98 0.305 13 0.002 2 1.5 1 0.5 0 0 200 400 600 Grain size (×10−4_cm) Sf (m)

Estimated from (4) with contact angle = 30° Estimated from (4) with contact angle = 63° Estimated from (4) with contact angle = 84° B35 from MGAM B60 from MGAM B150 from MGAM B320 from MGAM B35 from GAM B60 from GAM B150 from GAM B320 from GAM

Figure 7: Wetting front suction head s_f versus grain size. In all sands,s_f from MGAM was larger than s_f from GAM. The solid curve estimatess_f from a static contact angle of 30°, the dotted line represents the best contact angle (63°) withs_f from MGAM, and the dashed-dotted line represents the contact angle = 84°with s_f from GAM.

(9)

sizes of the sand columns and the equilibrium suction heads with the assumed eﬀective 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 ﬁtted relationship between grain sizes

and the equilibrium suction heads estimated by GAM and

MGAM with the ﬁtted eﬀective 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 eﬀective contact angle estimated by the GAM was also systematically greater than that by

the MGAM. The overestimated eﬀective contact angle from

GAM was because we neglected the dynamic eﬀect when

estimating the suction head. The eﬀective contact angle

esti-mated by GAM was the average dynamic contact angle dur-ing the inﬁltration process. Nevertheless, the contact angle

(63°) estimated by MGAM was still greater than that (30°)

measured in previous studies, even when the dynamic eﬀect was considered. One explanation for this is related to the contact angle hysteresis: the contact angle no longer corre-sponds to the static equilibrium angle but is larger when the liquid is advancing and smaller when the liquid is

reced-ing [45]. In general, durreced-ing the inﬁltration, the water

advances the air in the sand. Therefore, the contact angle during the inﬁltration should be in the form of an advancing contact angle, which is commonly greater than the equilib-rium contact angle.

4.2. MGAM Fitting Parameters for Dynamic Capillary Pressure. Hsu and Hilpert [29] pointed out that the parame-tersα and β in the MGAM are related to the dynamic eﬀect on the capillary pressure, and their values might not be the same for diﬀerent initial water contents, porosities, and grain

sizes. Figure 8 shows two trends in which theﬁtted α is

pro-portional to the porosity and inversely propro-portional to the grain size. From the empirical relationship proposed by

Stauﬀer [20] relating the dynamic coeﬃcient to measurable

porous media andﬂuid 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 sta-tistically increases as the grain size decreases. Therefore, our fitted values of the parameters were consistent with those of the aforementioned studies.

Limited numbers of Darcy-scale studies focus on

deter-mining the β values in sand. Weitz et al. [28] yielded β =

0 5 ± 0 1 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 ﬂat surfaces and in circular capillary tubes, which can

110 100 90 80 aˆ 0 100 200 300 400 500 Grain size (×10−4_cm) B60 B35 B150 B320 Y = −0.02X + 101.06 (a) 110 100 90 80 aˆ _B60 B35 B150 B320 0.42 0.44 0.46 0.48 Porosity (cm3_/cm3₎ Y = 189.09X + 11.32 (b)

Figure 8: MGAM parameter α versus (a) grain size and (b) porosity.

0.34 0.32 0.3 0.28 0.26 0.24 B320 B150 B60 B35 0 100 200 300 400 500 Grain size (×10−4_cm) 훽 Y = 0.3

Figure 9: MGAM parameter β versus grain size. The ﬁtted β was approximately 0.3.

(10)

be described by cosφ_eq− cos φ_dyn=αCaβ [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β = 0 3 resulted in the best ﬁt in downward

inﬁltration experiments with the same sand but varying

ponding heights. The sameβ value also worked well in

sim-ilar experiments but with various initial water contents [11].

In this study, β was set as one of the ﬁtting parameters.

Figure 9 shows that theﬁtted β 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ﬁtted β in

this study agreed with the values derived or determined from the aforementioned Darcy and pore-scale studies. We con-ﬁrmed again that 0.3 is a universal value for β for Darcy’s scale sand experiments.

5. Conclusions

We performed a series of downward inﬁltration experiments and compared the predictions derived from classical GAM

and MGAM as well as the values of ﬁtting parameters. In

this study, cumulative inﬁltration monotonously increased

with time, and a greater water ponding height showed more

cumulative inﬁltration. This kind of result was also reﬂected

in the wetting front depth. Furthermore, the sand column

with smaller grain size had a lower inﬁltration rate and

wetting front velocity, and a greater amount of time was spent on inﬁltration. 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 inﬁltration was superior to that of the classical GAM.

Because the classical GAM ignores the dynamic eﬀect on

estimating the suction head, the eﬀective contact angle

esti-mated 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 inﬁltration. Therefore, the equilibrium contact angle estimated by MGAM was greater than the intrinsic

con-tact angle. In addition, theﬁtted α was proportional to the

porosity and inversely proportional to the grain size in our

study, indicating that the dynamic eﬀect of capillary pressure

increased with decreasing grain size and increasing porosity.

Furthermore, theﬁtted β was approximately 0.3, which is

con-sistent with the values obtained from most of Darcy’s and pore-scale studies. This study conﬁrmed that 0.3 can be the

universal value ofβ for Darcy-scale sand experiments.

Data Availability

The data used to support theﬁndings of this study are

avail-able from the corresponding author upon request.

Conflicts of Interest

The authors declare that there is no conﬂict 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 ﬁnancially supporting this research

under Grant MOST 106-2628-M-002-009-MY3 and MOST 106-2923-E-009-001-MY3.

References

[1] D. A. Barry, J.-Y. Parlange, G. C. Sander, and M. Sivaplan, “A class of exact solutions for Richards' equation,” Journal of Hydrology, vol. 142, no. 1–4, pp. 29–46, 1993.

[2] G. Sposito,“Recent advances associated with soil water in the unsaturated zone,” Reviews of Geophysics, vol. 33, no. S2, pp. 1059–1065, 1995.

[3] W. H. Green and G. A. Ampt,“Studies on soil physics part I– the ﬂow of air and water through soils,” Journal of Agricul-tural Science, vol. 4, pp. 1–24, 1911.

[4] S. M. Hsu, C.-F. Ni, and P.-F. Hung,“Assessment of three inﬁl-tration formulas based on modelﬁtting on Richards equation,” Journal of Hydrologic Engineering, vol. 7, no. 5, pp. 373–379, 2002.

[5] G. Liu, J. R. Craig, and E. D. Soulis, “Applicability of the Green-Ampt inﬁltration model with shallow boundary condi-tions,” Journal of Hydrologic Engineering, vol. 16, no. 3, pp. 266–273, 2011.

[6] A. Van den Putte, G. Govers, A. Leys, C. Langhans, W. Clymans, and J. Diels,“Estimating the parameters of the Green–Ampt inﬁltration equation from rainfall simulation data: why simpler is better,” Journal of Hydrology, vol. 476, pp. 332–344, 2013.

[7] A. D. Feldman, “Hydrologic modeling system HEC-HMS: technical reference manual,” in US Army Corps of Engineers, Hydrologic Engineering Center, 2000.

[8] J. M. Laﬂen, W. Elliot, D. Flanagan, C. Meyer, and M. Nearing, “WEPP-predicting water erosion using a process-based model,” Journal of Soil and Water Conservation, vol. 52, no. 2, pp. 96–102, 1997.

[9] P. Tuppad, K. R. Douglas-Mankin, T. Lee, R. Srinivasan, and J. G. Arnold, “Soil and water assessment tool (SWAT) hydrologic/ water quality model: extended capability and wider adoption,” Transactions of the ASABE, vol. 54, no. 5, pp. 1677–1684, 2011. [10] F. Bouraoui and T. A. Dillaha,“ANSWERS-2000: runoﬀ and sediment transport model,” Journal of Environmental Engi-neering, vol. 122, no. 6, pp. 493–502, 1996.

[11] S.-Y. Hsu, V. Huang, S. Woo Park, and M. Hilpert,“Water inﬁltration into prewetted porous media: dynamic capillary pressure and Green-Ampt modeling,” Advances in Water Resources, vol. 106, pp. 60–67, 2017.

[12] V. T. Chow, D. R. Maidment, and L. W. Mays, Applied Hydrol-ogy, McGraw-Hill, 1988.

[13] S. P. Neuman,“Wetting front pressure head in the inﬁltration model of Green and Ampt,” Water Resources Research, vol. 12, no. 3, pp. 564–566, 1976.

[14] R. S. Mokady and P. F. Low,“The tension-moisture content relationship under static and dynamic conditions1,” Soil Sci-ence Society of America Journal, vol. 28, no. 4, 1964. [15] G. C. Topp, A. Klute, and D. B. Peters, “Comparison of

(11)

steady-state, and unsteady-state methods1,” Soil Science Society of America Journal, vol. 31, no. 3, p. 312, 1967. [16] E. Diamantopoulos and W. Durner, “Dynamic

nonequilib-rium of waterﬂow in porous media: a review,” Vadose Zone Journal, vol. 11, no. 3, 2012.

[17] W. C. Lo, C. C. Yang, S. Y. Hsu, C. H. Chen, C. L. Yeh, and M. Hilpert,“The dynamic response of the water retention curve in unsaturated soils during drainage to acoustic excita-tions,” Water Resources Research, vol. 53, no. 1, pp. 712–725, 2017.

[18] D. M. O'Carroll, T. J. Phelan, and L. M. Abriola,“Exploring dynamic effects in capillary pressure in multistep outflow experiments,” Water Resources Research, vol. 41, no. 11, 2005. [19] T. Sakaki, D. M. O'Carroll, and T. H. Illangasekare, “Direct quantification of dynamic effects in capillary pressure for drainage-wetting cycles,” Vadose Zone Journal, vol. 9, no. 2, 2010.

[20] F. Stauﬀer, “Time dependence of the relations between capil-lary pressure, water content and conductivity during drainage of porous media,” in IAHR Symp. on Scale Eﬀects in Porous Media. Thessaloniki Greece 29 Aug.-1. IAHR, Madrid, Spain, 1978.

[21] L. Zhuang, S. M. Hassanizadeh, C. Z. Qin, and A. de Waal, “Experimental investigation of hysteretic dynamic capillarity eﬀect in unsaturated ﬂow,” Water Resources Research, vol. 53, no. 11, pp. 9078–9088, 2017.

[22] S. M. Hassanizadeh and W. G. Gray,“Mechanics and thermo-dynamics of multiphaseﬂow in porous media including inter-phase boundaries,” Advances in Water Resources, vol. 13, no. 4, pp. 169–186, 1990.

[23] T. Annaka and S. Hanayama,“Dynamic water-entry pressure for initially dry glass beads and sea sand,” Vadose Zone Jour-nal, vol. 4, no. 1, pp. 127–133, 2005.

[24] D. A. DiCarlo,“Stability of gravity-driven multiphase ﬂow in porous media: 40 years of advancements,” Water Resources Research, vol. 49, no. 8, pp. 4531–4544, 2013.

[25] S. L. Geiger and D. S. Durnford,“Inﬁltration in homogeneous sands and a mechanistic model of unstableﬂow,” Soil Science Society of America Journal, vol. 64, no. 2, 2000.

[26] E. Pellichero, R. Glantz, M. Burns, D. Mallick, S. Y. Hsu, and M. Hilpert, “Dynamic capillary pressure during water inﬁl-tration: experiments and Green-Ampt modeling,” Water Resources Research, vol. 48, no. 6, 2012.

[27] N. Weisbrod, T. McGinnis, M. L. Rockhold, M. R. Niemet, and J. S. Selker,“Eﬀective Darcy-scale contact angles in porous media imbibing solutions of various surface tensions,” Water Resources Research, vol. 45, no. 4, 2009.

[28] D. A. Weitz, J. P. Stokes, R. C. Ball, and A. P. Kushnick, “Dynamic capillary pressure in porous media: origin of the viscous-ﬁngering length scale,” Physical Review Letters, vol. 59, no. 26, pp. 2967–2970, 1987.

[29] S. Y. Hsu and M. Hilpert,“Incorporation of dynamic capillary pressure into the Green-Ampt model for inﬁltration,” Vadose Zone Journal, vol. 10, no. 2, 2011.

[30] T. Tabuchi,“Inﬁltration and capillarity in the particle pack-ing,” Records Land Reclam Res., vol. 19, pp. 1–121, 1971. [31] P. G. D. Gennes, X. Hua, and P. Levinson,“Dynamics of

wet-ting: local contact angles,” Journal of Fluid Mechanics, vol. 212, no. 1, p. 55, 1990.

[32] P. Joos, P. van Remoortere, and M. Bracke,“The kinetics of wetting in a capillary,” Journal of Colloid and Interface Science, vol. 136, no. 1, pp. 189–197, 1990.

[33] T. E. Mumley, C. J. Radke, and M. C. Williams,“Kinetics of liquid/liquid capillary rise,” Journal of Colloid and Interface Science, vol. 109, no. 2, pp. 398–412, 1986.

[34] E. Schäﬀer and P.-z. Wong, “Contact line dynamics near the pinning threshold: a capillary rise and fall experiment,” Physi-cal Review E, vol. 61, no. 5, pp. 5257–5277, 2000.

[35] G. Camps-Roach, D. M. O'Carroll, T. A. Newson, T. Sakaki, and T. H. Illangasekare, “Experimental investigation of dynamic eﬀects in capillary pressure: grain size dependency and upscaling,” Water Resources Research, vol. 46, no. 8, 2010. [36] D. A. DiCarlo, “Experimental measurements of saturation overshoot on inﬁltration,” Water Resources Research, vol. 40, no. 4, 2004.

[37] D. Wildenschild, J. W. Hopmans, and J. Simunek,“Flow rate dependence of soil hydraulic characteristics,” Soil Science Soci-ety of America Journal, vol. 65, no. 1, 2001.

[38] M. Hilpert, “Velocity-dependent capillary pressure in theory for variably-saturated liquid inﬁltration into porous media,” Geophysical Research Letters, vol. 39, no. 6, 2012.

[39] Z. Wang, J. Feyen, and D. E. Elrick,“Prediction of ﬁngering in porous media,” Water Resources Research, vol. 34, no. 9, pp. 2183–2190, 1998.

[40] R. J. Glass, J. Y. Parlange, and T. S. Steenhuis,“Immiscible dis-placement in porous media: stability analysis of three-dimen-sional, axisymmetric disturbances with application to gravity-driven wetting front instability,” Water Resources Research, vol. 27, no. 8, pp. 1947–1956, 1991.

[41] M. D. Fredlund, D. G. Fredlund, and G. W. Wilson,“An equa-tion to represent grain-size distribuequa-tion,” Canadian Geotechni-cal Journal, vol. 37, no. 4, pp. 817–827, 2000.

[42] G. A. F. Seber and C. J. Wild, Nonlinear Regression, Wiley-Interscience, Hoboken, N.J, 2003.

[43] S. P. Friedman,“Dynamic contact angle explanation of ﬂow rate-dependent saturation-pressure relationships during tran-sient liquid ﬂow in unsaturated porous media,” Journal of Adhesion Science and Technology, vol. 13, no. 12, pp. 1495– 1518, 1999.

[44] P. W. Glover and E. Walker,“Grain-size to eﬀective pore-size transformation derived from electrokinetic theory,” Geophys-ics, vol. 74, no. 1, pp. E17–E29, 2009.

[45] L. Makkonen, “A thermodynamic model of contact angle hysteresis,” The Journal of Chemical Physics, vol. 147, no. 6, article 64703, 2017.

[46] L. H. Tanner,“The spreading of silicone oil drops on horizon-tal surfaces,” Journal of Physics D, vol. 12, no. 9, pp. 1473– 1484, 1979.

(12)

Hindawi www.hindawi.com Volume 2018 Journal of

Chemistry

Archaea

Hindawi www.hindawi.com Volume 2018

Marine Biology

Journal of Hindawi

www.hindawi.com Volume 2018

Biodiversity

International Journal of

Hindawi

Ecology

International Journal of

Hindawi www.hindawi.com Volume 2018 Hindawi www.hindawi.com

Applied &

Environmental

Soil Science

Volume 2018

Forestry ResearchInternational Journal of Hindawi www.hindawi.com Volume 2018 Hindawi www.hindawi.com Volume 2018

Geophysics

Environmental and Public Health Journal of Hindawi www.hindawi.com Volume 2018 Hindawi www.hindawi.com Volume 2018 International Journal of

Microbiology

Hindawi www.hindawi.com Volume 2018

Agriculture

Advances in Hindawi www.hindawi.com Volume 2018

Agronomy

Hindawi www.hindawi.com Volume 2018 International Journal of Hindawi www.hindawi.com Volume 2018

Meteorology

Advances in

Hindawi Publishing Corporation

http://www.hindawi.com Volume 2013 Hindawi www.hindawi.com

World Journal

Volume 2018 Hindawi

www.hindawi.com Volume 2018 Hindawiwww.hindawi.com Volume 2018

Hindawi

Geological ResearchJournal of

Analytical Chemistry International Journal of

Hindawi

Effects of the Grain Size on Dynamic Capillary Pressure and the Modified Green-Ampt Model for Infiltration

Research Article