Water and Vapor Movement with Condensation and Evaporation in a Sandy Column

soil hydraulic property model correctly describes retention and conductivity properties in the dry water content range.

One of the most widely used functions for describing un-saturated hydraulic properties is van Genuchten’s (1980) set of closed-form equations for the soil water retention curve, coupled with Mualem’s (1976) pore-size distribution model for the unsaturated hydraulic conductivity function. The van Genuchten (1980) model sometimes fails to describe soil water retention properties, however, and underestimates the unsatu-rated hydraulic conductivity at low water contents. Tuller and Or (2001) showed a relation between the unsaturated hydrau-lic conductivity and the pressure head based not only on the pore-size distribution, but also ? lm ? ow on the surface of soil particles, a process that dominated at low saturations. F ayer and Simmons (1995) modi? ed the van Genuchten soil water retention model to better describe soil water retention at low water contents and described the unsaturated hydraulic con-ductivity function by coupling their model with the Mualem (1976) model. The Mualem model contains the pore-connec-tivity coef? cient that accounts for the effects of pore connec-tivity and tortuosity, which can have a signi? cant impact on the gradient of the unsaturated hydraulic conductivity with respect to the pressure head. Although Mualem (1976), based on an analysis of 45 soil samples, suggested the average value of the pore-connectivity coef? cient to be 0.5, some recent stud-ies had utilized other values (e.g., Romano and Santini, 1999; Hopmans et al., 2002; Schaap et al., 2001).

Although numerical simulations of the simultaneous movement of liquid water, water vapor, and heat in soils with medium or high water contents have been performed in the past (e.g., Nassar et al., 1992b), few simulations have been performed for relatively dry soils. This is largely because of uncertainty about unsaturated hydraulic conductivity. The ob-jective of this study was to evaluate the vapor condensation experiments performed by Miyazaki (1976) using the Philip and de Vries (1957) theory and the Fayer and Simmons (1995) hydraulic property model implemented in the HYDRUS-1D code (?im?nek et al., 2008a,b; Saito et al., 2006). These laboratory experiments involved coupled movement of liquid water, water vapor, and heat in sandy columns under an im-posed temperature gradient. During the experiments, water vapor diffused from the hot and humid end of the columns and subsequently condensed at the cold dry end. In our study, the water retention curve parameters were ? rst ? tted to inde-pendently measured water retention curve data using the Fayer and Simmons (1995) model. The numerical model was then calibrated against experimental data of volumetric water con-tent pro? les by inversely estimating the pore-connectivity co-ef? cient to determine the unsaturated hydraulic conductivity function for low water contents while using independently ? t-ted water retention curve parameters. Further analyses of vari-ous processes were then performed using the calibrated model. Second, four components of the total water ? ux, including liquid and vapor ? uxes due to pressure head and temperature gradients, were evaluated using calculated pressure head and temperature pro?les. Third, condensation and evaporation rates were calculated from the mass balance and liquid water and water vapor ? uxes. Fourth, the impact of the enhancement factor on the coupled transport involving condensation and evaporation processes was investigated. Finally, the uncertainty of the inverse parameter estimation procedure using the avail-able data set was determined and the need for additional data to guarantee a unique solution was recognized. MATERIALS AND METHODS

Condensation Experiments

Miyazaki (1976) conducted one-dimensional laboratory column experiments to observe the diffusion of water vapor from hot and hu-mid air into cold and dry soil and to evaluate the various mechanisms involved in water vapor ? ow. The experiments were performed using Hamaoka dune sand, which had the typical particle size distribution of dune sand. It contained 3% of soil particles in the range of 0.002 to 0.02 mm, 7% of soil particles in the range of 0.02 to 0.1 mm, and 90% of soil particles >0.1 mm. The saturated hydraulic conductivity of 34.6 m d?1 was measured using the falling-head method. Hamaoka dune sand with an initial volumetric water content of 0.0045 m3 m?3 was packed uniformly at a dry bulk density of 1.6 g cm?3 in six acrylic columns with 10-cm height and 10-cm diameter. The sand columns were then placed vertically in a chamber with a constant air tempera-ture of 37°C and relative humidity of 85 to 90%. The tops of the columns were exposed to the moist air, while the closed bottoms of the columns were maintained at 20°C by circulating constant-tem-perature water along the base of the columns using a water pump. The walls of the columns were insulated using foam polystyrene to achieve one-dimensional heat ? ow. Since differences between the ob-served temperatures at the center and close to the wall of the column were <0.6°C (or 12°C m?1), the experiment was assumed to have one-dimensional heat ? ow.

Because of the difference in the vapor density between the two ends of the columns, water vapor diffused from the hot and humid end of the column toward the cold and dry end, where it condensed. Subsequently, the condensed liquid water moved upward due to the pressure head gradients. Temperature pro? les in the columns were monitored with thermocouples located at depths of 0, 2, 5, 8, and 10 cm at the center of the columns. The amount of cumulative water vapor diffusion into the column was measured based on the column weights. Identical experiments were repeated on all six columns to ob-tain water content pro? les at different times (2, 5, 10, 16, 21, and 30 d). Each sand column at those times was sectioned into seven 1- and 2-cm-thick layers to gravimetrically determine the water contents. Numerical Model

Unsaturated Water Flow

Darcy’s law for liquid water ? ux, q L, under time-variable tem-perature conditions has to consider both the ? ow driven by the tem-perature gradient, accounting for the temperature dependence of the surface tension, as well as the term driven by the pressure head gradi-ent (Philip and de Vries, 1957). Since the vapor density is also a func-tion of the pressure head and the temperature, Fick’s law for the vapor ? ux, q v, also consists of two (thermal and isothermal) components. Thus, the total water ? ux, q Total, in a variably saturated soil is the sum of the liquid water ? ux and the water vapor ? ux:

Total L v L L v v

L L v v

=1

h T h T

h T h T

q q q q q q q

h T h T

K K K K

z z z z =+=+++

????

??

?+???

??

????

??

[1]

708 SSSAJ: Volume 73: Number 3 ? May–June 2009

710

SSSAJ: Volume 73: Number 3 ? May–June 2009

()0.5L 12L 3L

b b b λθ=+θ+θ [14]where b 1, b 2, and b 3 are constants (for sands, b 1 = 0.228 W m ?1 K ?1,

b 2 = ?2.406 W m ?1 K ?1, and b 3 = 4.909 W m ?1 K ?1, as implemented in the HYDRUS-1D code). The two terms on the left-hand side of Eq. [13] represent changes in the energy content and the latent heat of the vapor phase. The terms on the right-hand side of Eq. [13] repre-sent soil heat ? ow by conduction, the convection of sensible heat with ? owing water, the transfer of sensible heat by diffusion of vapor, and the transfer of latent heat by diffusion of vapor. The volumetri

c heat capacity of the soil, C p [M L ?1 T ?2 K ?1] is de? ne

d as th

e sum o

f the volumetric heat capacities of the solid C n (= 1.92 MJ m ?3 K ?1), liquid C w (= 4.18 MJ m ?3 K ?1), and air C a (= 6.3 kJ m ?3 K ?1) phases multi-plied by their respective volumetric fractions θ (de Vries, 1958):

p n n w L a a C C C C =θ+θ+θ [15]

where θn is the volumetric fraction of the solid phase [L 3 L ?3]. The third term is usually neglected because it is signi? cantly smaller than the other two terms.

Initial and Boundary Conditions

Initial and boundary conditions were based on the condensation experiment. The initial volumetric water content was assumed to be uniform throughout the soil column (i.e., θL = 0.0045). The corre-sponding initial pressure head, calculated from the water retention curve (Eq. [5]), was approximately ?7 × 105 cm. The upper boundary condition for water ? ow was given by the total incoming water ? ux:

()()()Total L L v v 00, 0h T h T z q t q q q q t ==+++> [16]

The total surface ? ux, q Total (0,t ), was assumed to be equal to the change in weight of the entire column, thus re? ecting the condensa-tion of water vapor in the column. We ? tted a quadratic curve to the observed increase in weight of the column (cumulative condensation), and determined q Total (0,t ) from the time derivative of this curve.

It was assumed that there was no water ? ux at the lower bound-ary describing the closed end of the column:

()()()Total L L v v 1010,0 0h T h T z q t q q q q t =??=+++=> [17]

The sum of the liquid water ? ux (the sum of q L h and q L T ) and the vapor ? ux (the sum of q v h and q v T ) must be zero to satisfy the no-? ux boundary condition.

Since only measured temperature pro? les for the 6th and 16th days were available from Miyazaki (1976), the initial temperature was assumed to be a constant room temperature of 23.5°C throughout the soil pro? le. Dirichlet boundary conditions with speci? ed tempera-tures were considered for heat transport calculations. As it is unlikely that the boundary temperatures of the soil are equal to imposed tem-peratures (Farlow, 1993, p. 19–26), the upper and lower boundary temperatures were set equal to the average temperatures observed at those locations nearest to both boundaries, i.e., 36.3°C at the soil surface and 23.5°C at the soil bottom.

Governing Eq. [2] and [13] describing coupled movement of liquid water, heat, and water vapor were solved numerically using the HYDRUS-1D code that included vapor transport (?im ?nek et al., 2008a,b). The soil column was discretized uniformly into ? nite ele-ments of 0.1 cm. The time step was allowed to vary between the initial and maximum allowed time steps of 1 × 10?10 and 0.01 d, respectively.

Hydraulic Parameters

F igure 1a shows water retention

curves for the Hamaoka dune sand measured using both drainage and imbibition processes. Since Miyazaki (1976) did not provide the water re-tention curve, both processes were ob-served using a hanging water column. A time domain re? ectometry (TDR) probe (three 6.5-cm-long metallic rods of 0.15-cm diameter with 0.5-cm separations) and a tensiometer (having a 1.5-cm-long and 0.6-cm-diameter porous cup) with a pressure trans-ducer were installed horizontally in a soil sample with a height of 2 cm and

an internal diameter of 10 cm. After draining from saturation to a pressure head of ?60 cm, the imbibition pro-Table 1. Parameters used in numerical simulations.

Parameter

Equation?

Reference

Surface tension of soil water (γ), g s ?2γ = 75.6 ? 0.1425T ? 2.38 × 10?4T 2

Hillel (1971)Density of liquid water (ρw ), kg m ?3ρw = 1 ? 7.37 × 10?6(T ? 4)2+3.79 × 10?8(T ? 4)3

Hillel (1971)Saturated vapor density (ρvs ), kg m ?3ρvs = exp(31.37 ? 6014.79T abs ?1 ? 7.92 × 10?3T abs )T abs ?1 × 10?3Campbell (1985)Diffusivity of vapor in air (D a ), m 2 s ?1D a = 2.12 × 10?5(T abs /273.15)2Campbell (1985)

Relative humidity (H r )

H r = exp(hMg /RT abs )

Philip and de Vries (1957)Latent heat of vaporization of water (L w ), J kg ?1L w = 2.501 × 106 ? 2369.2T

Monteith and Unworth (1990)

? T , temperature; T abs , absolute temperature; h , pressure head; M , molecular weight of water; g , gravitational acceleration; R , universal gas constant.

Fig. 1. (a) The soil water retention curve of the Hamaoka dune sand, with the Fayer and Simmons (1995) model ? tted to the imbibition data; and (b) unsaturated hydraulic conductivity functions evaluated using pore-connectivity coef? cient (l ) values of 0.5, 4.91, 6.00, and 6.17.

SSSAJ: Volume 73: Number 3 ? May –June 2009 711

cess was subsequently performed by connecting the soil sample with a Mariotte tank and returning it to saturation. Finally, the drainage curve was measured and the ? nal water content at ?60 cm was deter-mined gravimetrically. The TDR measurements were then calibrated based on the ? nal water content and the cumulative amount of drain-age. Using this calibration, dielectric constants measured during the imbibition process were then converted to volumetric water contents. We used the water content values estimated by TDR in Fig. 1a. The water retention curve for soils with lower water contents was then measured using the pressure-plate apparatus for pressure heads lower than ?100 cm (see F ig. 1a) and assuming that hysteresis could be ignored for these low pressure heads.

Since the Hamaoka dune sand displayed a relatively important hysteretic behavior and because the condensation experiment rep-resented a wetting process during the entire period, the imbibition water retention curve was used in the calculations. The Fayer model (Eq. [5]) was ? tted to the observed data, resulting in values of θa = 0.027, θs = 0.325, α = 0.0656 cm ?1, and n = 4.71. Note that θs was about 0.06 smaller than the porosity value of 0.385 calculated from bulk and particle densities. This was probably partly due to the inevitably entrapped air during the drainage and imbibition processes. Furthermore, since we estimated θs based on the ? nal water content at h = ?60 cm and the cumulative drainage amount, θs might be slightly underestimated. Since the condensation experiment was performed in the pressure head range below ?20 cm, we believe that these possible experimental errors near saturation did not affect our calculations.

Figure 1b shows the relation between the pressure head, h , and the unsaturated hydraulic conductivity, K L h , derived by substituting Fayer’s water retention curve model (Eq. [5]) into Mualem’s model (Eq. [7]). The observed value of the saturated hydraulic conductivity (K s = 34.6 m d ?1) and the value of the pore-connectivity coef? cient suggested by Mualem (l = 0.5) were used to obtain the unsaturated hydraulic conductivity function. The gradient, dlog K Lh /dlog h , is larg-er for the higher pressure head range (?10 ≥ h ≥ ?100 cm) than for the lower pressure head range (h < ?100 cm).

When the a parameter in the enhancement factor is assumed to be known (i.e., equal to 8 (Eq. [11]), the pore-connectivity coef? cient l is the only remaining unknown parameter (when not assumed to be 0.5). Therefore, the value of l was optimized using an inverse analysis (e.g., ?im ?nek et al., 1998), during which the numerical solution was ? tted to volumetric water contents observed in the vapor condensa-tion experiment. Parameters θa , θs , α, n , and K s were kept constant at the values discussed above during the inverse analyses. The objective function, Φ, was de? ned as

()()2

L L 11

*,,j i

n n j i j i j i z t z t ==??Φ=θ?θ??

∑∑ [18]where θL *(z j , t i ) and θL (z j , t i ) were the observed and calculated water

contents at time t i for the spatial coordinate z j , respectively, n i was the number of measurement times, and n j was the number of measure-ment locations. Minimization of the objective function, Φ, was ac-complished using the Levenberg–Marquardt nonlinear minimization method (Marquardt, 1963).

The impact of the enhancement factor, η, on the optimized l pa-rameter in the unsaturated hydraulic conductivity, K L h , was evaluated by performing the model calibration with different values of the a pa-rameter (= 5, 8, or 15) in Eq. [11] (Fig. 2). We also demonstrated how

water content pro? les and water ? uxes change when different values

of the enhancement factor (with a = 5, 8, or 15 or η = 1) were used.

RESULTS AND DISCUSSION

Unsaturated Hydraulic Conductivity

The liquid water content pro? les were ? rst calculated with the original values of the pore-connectivity coef? cient (l = 0.5) suggested by Mualem (1976) and the enhancement factor (a = 8) suggested by Cass et al. (1984). Figure 3 compares the calculated and observed water content pro? les. Similar to the observed water contents, the calculated water contents also in-creased from the bottom due to vapor condensation; however, the shape of the simulated water content pro? les was quite dif-ferent from those observed. It was found that the simulated wa-

ter contents were much smaller than the observed values near

712

SSSAJ: Volume 73: Number 3 ? May–June 2009

the bottom of the soil column. For example, while the calculat-ed θL was 0.03, the observed value was 0.096 at a depth of 9.5 cm after 10 d. On the other hand, calculated water contents were overestimated in the middle of the column (e.g., calculat-ed θL = 0.025 vs. observed θL = 0.0075 at a depth of 5 cm after 10 d). This disagreement between observations and calculations might be the result of either overestimation of the upward liquid water ? ow or underestimation of the downward vapor ? ow.

Since our focus was mainly on the uncertainty in the un-saturated hydraulic conductivity, K L h , for lower water contents, we ? xed the water vapor ? ow parameters (i.e., a = 5, 8, and 15 in the enhancement factor) and attempted to decrease the upward liquid water ? ow by adjusting K L h through changing the pore-connectivity coef? cient, l . Optimized values of l us-ing the inverse analysis of the observed water contents were 6.00, 6.17, and 4.91 for a = 5, 8, and 15, respectively. The

relation between the pressure head, h , and the unsaturated hydrau-lic conductivity, K L h , is shown in Fig. 1b. Higher values of l lead to a faster decrease in K L h with de-creasing h . F igure 4a shows water content pro? les calculated with an optimized value of l = 6.17 (for a = 8). In this case, the calculated water contents agreed very well with the observed water contents. Figure 4b shows the corresponding pressure head pro? les. While pressure heads increased signi? cantly from the initial value (h = ?7 × 105 cm) in the lower part of the soil column (?30 ≤ h ≤ ?20 cm), the condensation process caused them to remain very low (?1 × 105 ≤ h ≤ ?3 × 104 cm) in the upper part.

Since similar results were ob-tained with the other two combina-tions of parameters a and l , these two parameters were mu-tually correlated and could not be estimated simultaneously from available data, i.e., water content pro? les only. For a given experimental setup, only one of these two parameters can be optimized, while the other parameter has to be estimated inde-pendently and ? xed during calibration. Notice, however, that differences between the optimized hydraulic conductivity func-tions for different enhancement factors were relatively small and probably below the precision of the available measurement techniques for low water contents, for which hydraulic con-ductivity measurements are extremely dif? cult to do.

The agreement between the measured and calculated water contents indicates that the relation of K L h and h in the Hamaoka dune sand at lower water contents can be described well using Mualem’s (1976) pore-size distribution model with the Fayer water retention curve model. The K L h (h ) function shows different gradients, dlog K L h /dlog h , for higher and lower pressure heads. This shape for the K L h (h ) function is similar to that proposed by Tuller and Or (2001), which considers ? lm ? ow for lower water contents as well as pore water ? ow for higher water contents. While the sand retains water in soil pores by capillary forces at higher water contents, at lower wa-ter contents water is retained on soil particles as water ? lms. These different retention mechanisms may alter the depen-dency of the hydraulic conductivity on the pressure head. A similar argument was used for solute mixing in terms of the hy-drodynamic dispersion in the dune sand (Toride et al., 2003). Solutes mix well in larger pores with higher pore-water veloci-ties and higher water contents. On the contrary, slow solute mixing, primarily by transverse diffusion, becomes dominant for lower water contents.

F igure 5 shows the unsaturated hydraulic conductiv-ity functions, K L h (h ), calculated using the van Genuchten–Mualem (VG) model (van Genuchten, 1980) with various values of l (?1 ≤ l ≤ 5). Although the V

G model shows differ-ent values of dlog K L h /dlog h for different values of l , it never resembles the shape of the Fayer model, which has different values of dlog K L h /dlog h

for different ranges of pressure heads.

SSSAJ: Volume 73: Number 3 ? May –June 2009

713

This indicates that the Fayer model not only modi? es the VG model for the soil water retention curve, but it also has a pro-found effect on the unsaturated hydraulic conductivity func-tion for lower water contents.

Since the numerical model that uses the unsaturated hy-draulic conductivity function K L h (h ) with the pore connectiv-ity parameter l = 6.17 that was optimized for the enhancement factor with a = 8 (Fig. 1b and 5) describes the experimental data (Fig. 4) well, the following discussion will be based on calcula-tions with this unsaturated hydraulic conductivity function. It needs to be emphasized here, however, that the optimized un-saturated hydraulic conductivity function may not be unique and could not be independently validated against directly mea-sured hydraulic conductivities at low water contents.

Temperature Pro? les

Figure 6 presents the simulated and observed temperature pro? les obtained during the experiment. A linear decrease in temperature between the upper and lower boundary values was reached after only 0.1 d. Observed temperatures after 6 and 16 d closely followed this linear distribution, except for small convex deviations at a depth of 9 cm after 6 d and 7 cm after 16 d. This small deviation from the linear pro? le can be ex-plained by the difference in thermal conductivities between the lower and upper part of the column as a result of their different water contents. For example, the convex deviation occurred at a depth of 7 cm after 16 d, which agreed well with the depth of the moisture front (Fig. 4a). The thermal conductivity cal-culated using Eq. [14] for the lower part of the column (θL > 0.075) was >1.4 W m ?1 K ?1, while it was <0.7 W m ?1 K ?1 for the upper part of the column (θL < 0.01).

Liquid Water and Water Vapor Fluxes

Once the pressure head and temperature pro? les were ob-tained, it was possible to evaluate the liquid and vapor ? uxes due to pressure head (isothermal) and temperature (thermal) gradients using Eq. [1]. Figure 7 shows the calculated total water ? ux pro? les, as well as their four components, after 6 and 16 d. Positive values represent upward movements, while negative val-ues represent downward movements.

The total water ? ux, q Total , was always directed downward (negative) throughout the col-umn from the warm end to the cold end, similarly as in previ-ous studies of nonisothermal ? ow in column experiments (Jones and Kohnke, 1952; Nassar et al., 1992a). While the total water ? ux, q Total , was large in the upper part of the soil column above the mois-ture front (?0.017 cm d ?1 after 6 d and ?0.02 cm d ?1 after 16 d), it was much smaller in the lower part of the column with higher water contents (Fig. 4a), and it was zero at the bottom of the column. The downward water movement oc-curred mainly as thermal water va-por ? ux, q vT , due to the large downward temperature gradient

(Fig. 6). The isothermal vapor ? ux, q vh , was small except near the moisture front (in the 9-cm depth after 6 d and the 7-cm depth after 16 d), where upward vapor ? ux was calculated due to the large positive pressure head gradient (Fig. 4b). Since q vh was much smaller than q vT , the actual vapor ? ux, q v , was down-ward after both 6 and 16 d.

Upward water movement occurred mainly as isothermal liquid water ? ux, q Lh , due to the upward pressure head gradi-ent, which developed at the bottom of the soil column as a result of vapor condensation. The isothermal liquid water ? ux had a peak at both times near the moisture front due to the extremely large pressure head gradients. The isothermal liquid water ? ux was also large at the bottom of the column after 16 d (0.067 cm d ?1) because of the large unsaturated hydraulic con-ductivity, K L h , for high water contents. As the water content gradually increased at the column bottom and the moisture front moved upward, so did the maxima and absolute values of q L h and q v T . The thermal liquid water ? ux, q L T

, showed

714

SSSAJ: Volume 73: Number 3 ? May–June 2009

downward water movement near the bottom of the column. Although q L T was small after 6 d, it increased after 16 d (?0.06 cm d ?1 at the bottom) due to a larger thermal hydraulic conductivity, K L T (Eq. [8] with the large K L h ). Because q L h was larger than q L T at every depth, the actual liquid water ? ux, q L , was always upward. The value of q L was negligibly small at depths above the moisture front after both 6 and 16 d.

Miyazaki (1976) calculated the increase in water contents at the bottom of the columns by evaluating the water vapor ? ux due to the temperature gradient using the Philip and de Vries theory and the Penman theory (Penman, 1940), which was the original theory of water vapor ? ow without the enhancement factor. The water content calculations signi? cantly overesti-mated the observed values, however, because of the dif? culty of calculating the simultaneous upward liquid ? ux resulting from the water content increase. Our numerical simulations of the coupled movement of water vapor, liquid water, and heat could consider simultaneously all four ? ux components and thus lead to better agreement between simulated and observed water contents.

Surface Boundary Condition

Both Milly (1984) and Saito et al. (2006) evaluated surface

water ? ux. The former used meteorological variables, while the latter utilized the surface energy balance. Since evaluating wa-

ter movement inside the soil column was the primary purpose of this study, the observed amount of water vapor diffused into the column, q Total (0,t ), was determined from weight changes in the entire soil column and used as the surface boundary con-dition for water ? ow. We assumed that q Total (0,t ) was equal to the sum of liquid and vapor ? uxes passing through the surface boundary, as described by Eq. [16]. The four components of the total ? ux at the soil surface were then determined using the pressure head and temperature gradients, as described by Eq. [1]. Each ? ux may have a different magnitude and direction depending on the pressure head and temperature gradients at the soil surface. Since the temperature in the soil column quickly established steady-state conditions (during about 0.1 d as shown in Fig. 6), the pressure head at the soil surface was the only other variable that could adjust q Total to be equal to the

observed diffusion ? ux of water vapor into the column during calculations.

F igure 8a presents the calcu-lated surface pressure heads dur-ing the experiment. Although not clearly visible in Fig. 8a, at the be-ginning of the experiment the sur-face pressure head ? rst decreased quickly from the initial pressure head of ?7 × 105 cm, and then in-creased because of the rapid change in the surface temperature from the initial value of 23.5°C to the im-posed boundary value of 36.3°C. After that, the surface pressure head started increasing steadily un-til it reached a value of around ?1.5 × 105 cm. The increasing pressure head re? ects vapor condensation at

the soil surface. As the amount of condensed water was rela-tively small, water content increases were barely noticeable.

Figure 8b shows the four components of the surface ? ux q Total (0,t ) as a function of time. Liquid water ? uxes (q L h and q L T ) were signi? cantly smaller (and in fact negligible) than vapor ? uxes (q v h and q v T ) because the hydraulic conductivi-ties for liquid water were several orders of magnitude smaller (K L h ≈ 10?19 cm d ?1 and K L T ≈ 10?15 cm 2 K ?1 d ?1) than vapor hydraulic conductivities (K v h ≈ 10?6 cm d ?1 and K v T ≈ 10?1 cm 2 K ?1 d ?1). Although the direction of q v h was changing at the beginning of the experiment, due to oscilla-tions of the surface pressure head in time (Fig. 8a), downward q v T and upward q v h dominated the overall water ? ux q Total (0, t ), and de? ned the water vapor diffusion into the column.

Condensation and Evaporation Rates

The evaporation rate (E ) at each internal discretization node can be calculated with a discretized form of Eq. [3]:

11/21/2L L L 1/2L 1/2j j j j j

i i

i i i q q E t z ????+θ?θ?=?+ΔΔ [19]

where the subscript i represents the position of a lattice point and the superscripts j and j ? 1 denote the current and previous time points, respectively.

The evaporation rate at locations i = 1 + 1/4 and i = N ? 1/4 can be evaluated using a method similar to Eq. [19], where i = 1 and N represent the lower and upper boundary nodal points, respectively:

11/21/2L11/4L11/4L1L11/211/4

2j j j j j q q E

t z ???++++θ?θ?=?+

ΔΔ [20]11/21/2

L 1/4L 1/4L L 1/214

2

j j j j j

N N N N N q q E

t z ???????θ?θ?=?+

ΔΔ [21]Values at i = 1 + 1/4 and i = N ? 1/4 (e.g., θL1+1/4) can be obtained using a linear interpolation of adjacent node values. Note that q L1 does not have to be zero because, as described by Eq. [17], only the total boundary ? ux is assumed to be zero.

SSSAJ: Volume 73: Number 3 ? May –June 2009

715

The evaporation rate at the surface boundary, E 1, is evaluated by lin-early extrapolating E 1+1/4 and E 2 toward the surface. A similar lin-ear extrapolation from E N ?1/4 and E N ?1 is made to estimate E N .

Figure 9a displays evaporation rate pro? les calculated using Eq. [19–21]. Note that negative evapo-ration represents vapor condensa-tion. The maximum evaporation rate was located at a depth of 9 cm after 6 d, and moved upward to a depth of 7 cm after 16 d. This lo-calized evaporation resulted in the maximum downward thermal va-por ? ux, q v T , behind the moisture front as shown in Fig. 7. Figure 9b shows relative humidity pro? les

(Table 1) at 6 and 16 d. The relative humidity decreased from almost unity near the bottom of the column to 0.98 at the moisture front. The location of the maximum evaporation rate corresponded exactly with the position of the abrupt decrease in the relative humidity. As shown in Fig. 4b, a large pressure head drop occurred at the moisture front (h = ?2 × 104 cm), resulting in a decrease in the relative humidity.

Vapor condensation started just below the maximum evaporation peak, where the relative humidity increased to al-most unity. Vapor had to condensate at this position because the air could not hold all the water vapor coming from the up-per location. As shown in Fig. 6, vapor moved farther down and eventually condensed at the cold end. The maximum condensa-tion (minimum evaporation) rate occurred at the bottom of the column: ?0.123 d ?1 after 6 d and ?0.085 d ?1 after 16 d (Fig. 9a).

The total amount of vapor condensation in an entire soil column can be obtained by integrating the condensation rate pro? le (negative evaporation rate in Fig. 9a) with depth and time. Figure 10 shows the cumulative amount of condensation as a function of time. The total quantity of condensation in the soil pro? le includes not only the condensation at the bottom, but also that in the rest of the soil pro? le. The condensation at the bottom dominated initially and accounted for 40% of the total condensation during the ? rst 5 d. Condensation inside the column subsequently increased and reached about 75% of the total condensation (25% at the bottom) after 20 d. The observed cumulative amount of water vapor diffused into the column was also plotted in Fig. 10. Note that the total amount of condensation in the column was greater than the cumula-tive amount of water vapor diffused into the column due to internal evaporation (Fig. 9a). In other words, the difference between the total amount of condensation and the cumulative amount of diffusion into the column was equal to the amount of evaporation inside the column. The amount of internal evaporation accounted for around 25% of the total condensa-tion, regardless of the time.

Enhancement Factor

The impact of the enhancement factor, η, can also be demonstrated using water content pro? les. F igure 11 shows

the observed water content pro? les as well as those calculated with different enhancement factors after 30 d with an opti-mized value of l = 6.17. In calculations without the enhance-ment factor (η = 1), calculated water contents were underesti-mated in the lower part of the column and overestimated in the upper part. More vapor condensed in the upper part of the column due to underestimation of the vapor ? ow. This simula-tion clearly demonstrates the importance of considering the enhancement factor. Calculated water content pro? les were similar when different values of the a coef? cient (a = 5, 8, and 15) were used to calculate the enhancement factor (Fig. 11). Figures 12a and 12b show pro? les of the liquid water and water vapor ? uxes and the evaporation rate, calculated with a = 8 and 15. Both the downward vapor ? ux and the upward liquid water ? ux increased as the enhancement factor increased, and so did both the maximum evaporation and condensation rates. The amount of water circulating in the soil column increased as the enhancement factor increased. These two opposing effects

(evaporation and condensation) compensated for each other,

716

SSSAJ: Volume 73: Number 3 ? May–June 2009

making the impact of the enhancement factor less important for the water content pro? le in the condensation process.

We have also evaluated the impact of η on the estimation of the pore-connectivity coef? cient, l . Optimized l values were 6.00, 6.19, and 4.91 for the enhancement factor a parameters of 5, 8, and 15, respectively. As shown in Fig. 1b, there is only a relatively small range of uncertainty in the unsaturated hydraulic conductiv-ity function due to differences in η for lower pressure heads. As discussed above, however, these two parameters are mutually high-ly correlated for given experimental data and only one of them can be optimized independently. Additional data need to be measured so that both the pore-connectivity parameter l and the enhance-ment factor η can be estimated simultaneously.

CONCLUSIONS

We have numerically analyzed the vapor condensation experiments reported by Miyazaki (1976). These experiments involved vapor diffusion and condensation in laboratory sand columns due to an imposed temperature gradient. The move-

ment of water vapor and liquid water in the sand column were evaluated using the Philip and de Vries (1957) model. The cou-pled transport of water vapor, liquid water, and heat were calcu-lated using the HYDRUS-1D code that included vapor ? ow.

The numerical model was calibrated against experimen-tal data by optimizing the pore-connectivity coef? cient, l , of the unsaturated hydraulic conductivity function for different values of the enhancement factor η. The Fayer and Simmons (1995) model, describing soil hydraulic properties with an estimated l parameter of 6.17 (for the enhancement factor with a = 8), provided a very good agreement between simu-lated and observed water contents. The rate of decrease in the unsaturated hydraulic conductivity with decreasing pressure head, described by the Fayer model, differed for high and low water contents, re? ecting capillary pore water ? ow and ? lm ? ow, respectively. It was concluded that the Fayer model with the calibrated l parameter of Mualem’s pore-size distribution model worked well for simulating water ? ow in sandy soils at low water contents. Due to the lack of directly measured unsaturated hydraulic conductivities, however, the calibrated hydraulic conductivity function could not be independently validated. Mutual correlation between the pore-connectivity coef? cient and the enhancement factor also raises a question of uniqueness of the optimized unsaturated hydraulic conductiv-ity function.

Additionally, we quantitatively evaluated the four compo-nents of the total water ? ux: liquid and vapor ? uxes driven by pressure head and temperature gradients. The evaporation and condensation rates inside the soil column were described further by analyzing these water ? uxes. Water vapor that en-tered the soil column from the hot surface moved downward due to the temperature gradient and condensed at the cold bottom end of the column. Liquid water subsequently moved upward due to the pressure head gradient and evaporated at the moisture front where the relative humidity decreased from 1 to 0.98. Vapor generated by evaporation then moved down-ward together with vapor coming from the soil surface. A cer-tain amount of vapor condensed just below the location with maximum evaporation at the moisture fronts, and the maxi-mum condensation occurred at the cold bottom of the column

(25–40% of the total condensation in the column). Liquid water and water vapor circulated between the bottom of the column and the moisture front, accompanied by condensation and evaporation pro-cesses in the sand column.

The impact of the enhance-ment factor on vapor ? uxes inside the soil column was also demon-strated. Downward vapor ? ux was signi? cantly underestimated when the enhancement factor was ne-glected. Calculated water content pro? les were very similar when the a parameter from the function of the enhancement factor was >5. Higher enhancement factors pro-

duced larger upward liquid water

? uxes and evaporation rates, as well as downward vapor ? uxes and condensation rates. These two opposing effects compen-sated for each other, making the water content pro? les in the condensation experiment less sensitive to the enhancement factor. This indicates that water content pro? les do not pro-vide enough information to simultaneously estimate both the unsaturated hydraulic conductivity function and the enhance-ment factor. Additional data are needed so that both the pore-connectivity parameter l and the enhancement factor η can be estimated simultaneously. ACKNOWLEDGMENTS

We would like to thank Dr. Tsuyoshi Miyazaki from the University of Tokyo for his kindness in offering us his experimental data as well as the Hamaoka dune sand that was used in his column experiments for additional measurements of the water retention curve. This study was supported in part by the Terrestrial Sciences Program of the Army Research Of? ce (Terrestrial Processes and Landscape Dynamics and Terrestrial System Modeling and Model Integration), by the National Science F oundation (NSF) Biocomplexity programs no. 04-10055 and NSF DEB 04-21530.

REFERENCES

Cahill, A.T., and M.B. Parlange. 1998. On water vapor transport in ? eld soils.

Water Resour. Res. 34:731–739.

Campbell, G.S. 1985. Soil physics with BASIC. Elsevier, New York.

Cass, A., G.S. Campbell, and T.L. Jones. 1984. Enhancement of thermal water vapor diffusion in soil. Soil Sci. Soc. Am. J. 48:25–32.

Chung, S.O., and R. Horton. 1987. Soil heat and water ? ow with a partial surface mulch. Water Resour. Res. 23:2175–2186.

de Vries, D.A. 1958. Simultaneous transfer of heat and moisture in porous media. Trans. Am. Geophys. Union 39:909–916.

Eching, S.O., and J.W. Hopmans. 1993. Optimization of hydraulic functions from transient out? ow and soil water pressure data. Soil Sci. Soc. Am. J.

57:1167–1175.

F arlow, S.J. 1993. Partial differential equations for scientists and engineers.

Dover, New York.

Fayer, M.H., and C.S. Simmons. 1995. Modi? ed soil water retention functions for all matric suctions. Water Resour. Res. 31:1233–1238. Fujimaki, H., and M. Inoue. 2003. A transient evaporation method for determining soil hydraulic properties at low pressure. Vadose Zone J. 2:400–408. Hillel, D. 1971. Soil and water: Physical principles and processes. Academic Press, New York.

Hopmans, J.W., J. ?im?nek, and K.L. Bristow. 2002. Indirect estimation of soil thermal properties and water ? ux using heat pulse probe measurements: Geometry and dispersion effects. Water Resour. Res. 38(1):1006, doi:10.1029/2000WR000071. Inoue, M., J. ?im?nek, J.W. Hopmans, and V. Clausnitzer. 1998. In situ estimation of soil hydraulic functions using a multistep soil-water extraction technique. Water Resour. Res. 34:1035–1050.

Jones, H.E., and H. Kohnke. 1952. The in? uence of soil moisture tension on vapor movement of soil water. Soil Sci. Soc. Am. Proc. 16:245–248. Klute, A., and C. Dirksen. 1986. Hydraulic conductivity and diffusivity: Laboratory methods. p. 687–734. In A. Klute (ed.) Methods of soil analysis. Part 1. Physical and mineralogical methods. 2nd ed. Agron.

Monogr. 9. ASA and SSSA, Madison, WI.

Marquardt, D.W. 1963. An algorithm for least-squares estimation of nonlinear parameters. SIAM J. Appl. Math. 11:431–441.

Millington, R.J., and J.M. Quirk. 1961. Permeability of porous solids. Trans.

Faraday Soc. 57:1200–1207.

Milly, P.C.D. 1984. A simulation analysis of thermal effects on evaporation.

Water Resour. Res. 20:1087–1098.

Miyazaki, T. 1976. Condensation and movement of water vapor in sand under temperature gradient. (In Japanese.) Trans. Jpn. Soc. Irrig. Drain. Rural Eng. 61:1–8.

Monteith, J.L., and M.H. Unworth. 1990. Principles of environmental physics.

2nd ed. Edward Arnold, London.

Mualem, Y. 1976. A new model for predicting the hydraulic conductivity of unsaturated porous media. Water Resour. Res. 12:513–522. Nassar, I.N., A.M. Globus, and R. Horton. 1992a. Simultaneous soil heat and water transfer. Soil Sci. 154:465–472.

Nassar, I.N., and R. Horton. 1989. Water transport in unsaturated nonisothermal salty soil: II. Theoretical development. Soil Sci. Soc. Am. J. 53:1330–1337. Nassar, I.N., and R. Horton. 1992. Simultaneous transfer of heat, water, and solute in porous media: I. Theoretical development. Soil Sci. Soc. Am. J. 56:1350–1356. Nassar, I.N., and R. Horton. 1997. Heat, water, and solute transfer in unsaturated porous media: I. Theory development and transport coef? cient evaluation. Transp. Porous Media 27:17–38.

Nassar, I.N., R. Horton, and A.M. Globus. 1992b. Simultaneous transfer of heat, water, and solute in porous media: II. Experiment and analysis. Soil Sci. Soc. Am. J. 56:1357–1365.

Nimmo, R.J. 1990. Experimental testing of transient unsaturated ? ow theory at low water content in a centrifugal ? eld. Water Resour. Res. 26:1951–1960. Nimmo, J.R., and E.E. Miller. 1986. The temperature dependence of isothermal moisture vs. potential characteristics of soils. Soil Sci. Soc.

Am. J. 50:1105–1113.

Noborio, K., K.J. McInnes, and J.L. Heilman. 1996. Two-dimensional model for water, heat, and solute transport in furrow-irrigated soil: I. Theory.

Soil Sci. Soc. Am. J. 60:1001–1009.

Penman, H.L. 1940. Gas and vapor movement in soil: I. The diffusion of vapors in porous solids. J. Agric. Sci. 30:437–462.

Philip, J.R., and D.A. de Vries. 1957. Moisture movement in porous materials under temperature gradients. Trans. Am. Geophys. Union 38:222–232. Romano, N., and A. Santini. 1999. Determining soil hydraulic functions from evaporation experiments by a parameter estimation approach: Experimental veri? cations and numerical studies. Water Resour. Res. 35:3343–3359. Rossi, C., and J.R. Nimmo. 1994. Modeling of soil water retention from saturation to oven dryness. Water Resour. Res. 30:701–708.

Saito, H., J. ?im?nek, and B.P. Mohanty. 2006. Numerical analysis of coupled water, vapor, and heat transport in the vadose zone. Vadose Zone J. 5:784–800. Scanlon, B., K. Keese, R.C. Reedy, J. ?im?nek, and B.J. Andraski. 2003. Variations in ?ow and transport in thick desert vadose zones in response to paleoclimatic forcing (0–90 kyr): F ield measurements, modeling, and uncertainties.

Water Resour. Res. 39(7):1179, doi:10.1029/2002WR001604. Schaap, M.G., F.J. Leij, and M.Th. van Genuchten. 2001. ROSETTA:

A computer program for estimating soil hydraulic parameters with

hierarchical pedotransfer functions. J. Hydrol. 251:163–176.

?im?nek, J., M. ?ejna, H. Saito, M. Sakai, and M.Th. van Genuchten. 2008a. The HYDRUS-1D software package for simulating the movement of water, heat, and multiple solutes in variably saturated media. Version 4.0. HYDRUS Softw. Ser. 3. Dep. of Environ. Sci., Univ. of California, Riverside.

?im?nek, J., M.Th. van Genuchten, and M. ?ejna. 2008b. Development and applications of the HYDRUS and STANMOD software packages, and related codes. Vadose Zone J. 7:587–600.

?im?nek, J., O. Wendroth, and M.Th. van Genuchten. 1998. Parameter estimation analysis of the evaporation method for determining soil hydraulic properties. Soil Sci. Soc. Am. J. 62:894–905.

Toride, N., M. Inoue, and F.J. Leij. 2003. Hydrodynamic dispersion in an unsaturated dune sand. Soil Sci. Soc. Am. J. 67:703–712.

Tuller, M., and D. Or. 2001. Hydraulic conductivity of variably saturated porous media. Film and corner ? ow in angular pore space. Water Resour.

Res. 37:1257–1276.

van Genuchten, M.Th. 1980. A closed-form equation for predicting hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44:892–898.

SSSAJ: Volume 73: Number 3 ? May –June 2009 717