On the importance of hysteresis and heterogeneity in the numerical simulation of unsaturated flow

Soil hydraulic properties in field scale porous media are distributed with large spatial heterogeneity. Thus, detailed understanding of micro-physical properties leading to macro-scale heterogeneity is important for generalizing purposes in modeling. Hysteretic water retention properties greatly affect solute transport under unsaturated condition. In this study, laboratory infiltration experiments were carried out to investigate the importance of hysteresis in numerical modeling of water flow under heterogeneous distribution of hydraulic conductivity. Soil water content was observed by time domain reflectometry using small printed circuit board probes (PCBP). The PCBP observations were compared to results from numerical simulations involving saturated-unsaturated seepage analysis and two-phase flow with/without hysteresis effects. Observations were compared to both uniform and heterogeneous parameter modeling. Differences between unsaturated-saturated and two-phase flow were quite small in the numerical results for soil water content. The uniform parameter model, however, could not properly reproduce experimental water content change. The numerical results of water content change including hysteresis matched better observations for the unsaturated-saturated and two-phase flow cases. Results displayed that it is important to take heterogeneity and hysteresis into consideration in the numerical modeling in both two-phase and saturated-unsaturated flow analysis. (Less)


INTRODUCTION
Field soils and aquifers are known to have a large heterogeneity. This heterogeneity affects water flow and solute transport through physical parameters such as hydraulic conductivity. In order to improve the understanding of heterogeneity effects on water flow, many researchers have performed theoretical and experimental studies on subsurface water flow with respect to physical properties. Wildenschild and Jensen (1999a, b) studied two-dimensional flow in a heterogeneous hydraulic conductivity field by experiments and numerical simulation. Nakagawa et al. (2003) conducted numerical simulation and laboratory experiment to investigate solute transport in a heterogeneous and unsaturated flow field. However, simulation results did not agree very well with experiments and more studies were necessary. Experimental data and numerical methods for differences between local and global water contents have not been fully examined. Possibly a way to improve simulation results may be to not only examine the effects of hydraulic conductivity distribution but also differences in local soil water content for the infiltration process. In many cases it is also important to visualize the physical processes. The property of water flow in the unsaturatedheterogeneous porous media is one such application. Saito and Kawatani (2005) showed that the effects of pore air were negligible during ponding and infiltration. However, during irrigation and subsequent drainage, it is important to include the hysteresis effects in the soil water retention curve for infiltration. Effects of hysteresis on unsaturated water content were investigated by, e.g., Kaluarachchi and Parker (1987). There are still, however, few studies on the hysteresis effects on the heterogeneous flow field.
Various techniques for investigating water and solute transport in heterogeneous soil have traditionally been used. Dye tracer is a common technique. Heilig et al. (2003) applied dye tracer near the soil surface. After 2-6 weeks, blue dye was examined by excavating a trench and photographing the exposed soil surface. Gish et al. (2004) used dye tracer to determine the relative importance of field scale matrix and preferential flow processes on transport of mobile chemicals under steady state infiltration. Akhtar et al. (2003a) used dye tracer to visualize flow pattern for unsaturated conditions. Kung et al. (2005) conducted field scale experiments by using an improved tile drain monitoring protocol to measure the mass flux breakthrough patterns of conservative tracers. Akhtar et al. (2003b) investigated the transport of Cl − and Li + in undisturbed soil columns. Boll et al. (1997) used non-adsorbing tracer application, wick, and gravity pan sampler data to obtain frequency distribution of water and solute transport properties. However, these studies are indirect methods and do not include direct observation of the infiltration pattern of water. Zevi et al. (2005) used bright field microscopy to determine the position and movement of water and colloids. DiCarlo et al. (1997) presented a technique which used synchrotron X-rays from the Cornell High Energy Synchrotron Source to measure three-phase fluid saturation on the time scale of seconds. Deinert et al. (2004) presented a method for using real-time neutron radiography to measure rapidly changing moisture profiles in porous media. In an interesting study, Wells et al. (2003) used needle penetrometer to measure the depth of the wetting front along a single transect of the soil, immediately after each rainstorm. Wet sand transmits more light than dry sand. Accordingly, Bauters et al. (1998) used light transmission to quantify moisture content. In the present study, infiltration behavior of the unsaturated heterogeneous zone was investigated using PCBP (Printed Circuit Board Probes; e.g., Nissen et al., 1999).
In view of the above, the objective of the present study was to investigate the importance of including hysteresis and heterogeneity in the numerical modeling for infiltrating water through heterogeneous hydraulic conductivity fields under step-wise ponding conditions. In the laboratory experiments, water content was measured by TDR (Time Domain Reflectometry) using PCBP. Results from the laboratory experiments were reproduced by numerical simulation with saturated-unsaturated seepage analysis. In order to improve the understanding of effects of air movement in the flow field, two-phase flow simulation analysis was also carried out. To evaluate the importance of hysteresis in the infiltration modeling, both simulations included this aspect. Additionally, to confirm the hydraulic conductivity heterogeneity effects on water flow, results from a uniform parameter model was compared with heterogeneous model results. For the above purpose, previous results (Nakagawa and Saito, 2008) were also used in this paper together with new numerical results.

LABORATORY EXPERIMENTS
Water infiltration through heterogeneous hydraulic conductivity fields Figure 1 shows a schematic illustration of the experimental flow tank where infiltration observations were made. Field soil was air-dried and passed through a 2 mm sieve and then separated into a diameter range of 5 classes (a: 0.1-0.2, b: 0.2-0.4, c: 0.4-0.6, d: 0.6-0.8, and e: 0.8-1.2 mm). The soil water content for each of these soil samples were adjusted to 10% water content. The saturated hydraulic conductivity for the samples for a, b, c, d, and e was 1.1 × 10 −3 , 2.3 × 10 −3 , 3.2 × 10 −3 , 5.6 × 10 −2 , and 8.3 × 10 −2 cm s −1 , respectively. Samples a, b, and c were soils taken from the Kagoshima University experimental field area and samples d and e were sand. They were re-packed into fixed positions in the flow tank with a density of 0.9 g cm −3 for samples a-c and 1.3 g cm −3 for samples d and e. The volumetric water content just after re-packing was 9% for samples a-c and 13% for samples d and e. The field generator in the groundwater modeling software PMWIN (Processing MODFLOW for Windows) (Chiang and Kinzelbach, 2001) was used to generate a field of randomly distributed 5 × 5 cm soil blocks having different saturated hydraulic conductivity in the flow tank ( Figure 1). The PCBPs were 6.0 cm long, 1.6 cm wide, and about 1 mm thick as shown in Figure 1. The length of the needles printed on the board was 15 cm in order to maintain a good measuring accuracy for the water content. The PCBPs were installed from the back surface of the flow tank at a right angle to the water flow so that the flow disturbance was minimized. The main water flow was vertically downward. The PCBPs were installed at 8 points in the flow tank (nos. 1 to 8 in Figure 1). TDR100 (Campbell Scientific, Inc.) was used as cable tester connected with the PCBPs by a multiplexer. The lowest point of the flow tank was set as a constant water table.
A total of 3 × 3000 cm 3 infiltration water was applied from the top of the flow tank with 80 min intervals. Infiltration patterns were photographed by a digital camera at prescribed time intervals. When the infiltrated water front reached the lowermost part of the flow tank, the experiment was halted.
Soil water content analysis using PCBP As discussed in the previous chapter, several techniques can be used for moisture and concentration observations in heterogeneous soils. The advantage of PCBPs is that water content patterns can be displayed at a small spatial scale. Since inexpensive probes can be mass-produced it is possible to take observations at many points. Consequently, the system only requires multiplexor and a TDR together with the probes. Although not introduced in the present study, since electrical conductivity can also be measured by the probe, it can potentially be used to monitor transport of both water and solutes. Even if not well suited in large-scale field studies, it is advantageously used in laboratory scale experiments and also in small-scale field studies. For example, probes could have been inserted along an excavated trench according to Boll et al. (1997). Although the probe is quite small, it is necessary to insert it so that the flow field is not disturbed. In total 256 discrete data points from one reflection wave were temporary saved on a data logger. These data were then transferred to a PC hard disk at 20 min intervals. The wave data were analyzed with the soil analysis software WinTDR  to obtain apparent dielectric constants. By using precedent derived calibration curves for all PCBPs, apparent dielectric constants were converted to volumetric water content.

Governing equations
The following equations for two-phase flow of air and water through porous medium were used for generating the finite element model (e.g., Meiri, 1981): (1) (2) where t = time; x i = Cartesian coordinates (x 3 being the vertical coordinate oriented positive upward); η = porosity; β a = formation volume factor of air (β a = , T = absolute temperature of air, T s = 273.2 K, p as = 1.0 atm = 101.325 kPa); ρ w = density of water; k rw , k ra = relative permeability of water and air, respectively; k ij = intrinsic permeability tensor; p a , p w = pressure of air and water, respectively; p c = capillary pressure (= p a − p w ); μ a , μ w = viscosities of air and water, respectively; S w = degree of water saturation; and g = acceleration of gravity. The boundary conditions for prescribed flux are: water: air: where n i denotes the components of the unit vector normal to the boundary; q wb and q ab are the prescribed discharges of out flowing water and air, respectively. The governing equation of saturated-unsaturated seepage analysis (onephase flow analysis) is given by ignoring Equation (2)

Water retention curve and relative permeability
Water retention curves (= drying phase) were estimated for the five kinds of soil samples which filled up the flow tank by the following method. At first, samples were filled in 20 cm long columns with the same density as experiments. The volumetric soil water content distribution was then measured by a drainage column test. For samples c, d, and e, the following van Genuchten (1980) equation was employed to estimate main drying curve and the parameters were identified by fitting to observed data: where S e is the effective saturation; S w is the degree of saturation; S r is the residual saturation; S f is the saturation at p c = 0; h c is the capillary pressure head (= p c /ρ w g); α, n, and m (= 1 − 1/n) are parameters. Here, S f is assumed to be 1.0, and judging from measured drying-phase curve, S r is conveniently assumed to be 0. Thus, α and n are targets of identification in this experiment. Due to difficulties to estimate main drying curves for samples a and b, these were estimated by a stochastic approach based on the grain size distribution curve of soil samples used by Kitamura et al. (1998), and VG (van Genuchten) parameters were identified by model fitting.
Main wetting curves were estimated by the method of Luckner et al. (1989) for samples c, d, and e. On the other hand, main wetting curves for samples a and b were determined based on the initial conditions of the experiment. The VG parameters and porosity for each sample type are listed in Table I and examples of water retention curves are shown in Figure 2.
The relative permeability and effective saturation are related as Mualem (1976); In general, the relationship between degree of saturation and capillary pressure depends on the initial water content or the history of wetting and drying process. The effects of hysteresis will develop as the wetting and drying cycle is repeated in a regular manner. In order to estimate the continuous wetting and drying curves (scanning curves), the method proposed by Scott et al. (1983) was employed.  Table II. In Run-0 and Run-0h, saturated hydraulic conductivity and water retention curves were unified, and the effects of pore air were ignored. The volumetric water content θ (= ηS w ) for each sample was adjusted from 11 to 13%. On the other hand, the measured initial volumetric water content of sample d was below 10%. Consequently, the initial S w was set to 20%, and the status 30 min later was considered as initial condition. In addition, the initial water pressure distribution was supposed to be on the main wetting curve and in response to initial saturation for each sample.    Run-1h 1 Hysteresis Run-2h 2 no. 5 improved the agreement substantially for the Run-1 and Run-2. The Run-2h is reproducing the observed value rather than the Run-1h slightly. A similar but smaller improvement occurred for point no. 7. According to Figure  2, the sample with lower permeability (sample a), showed a large hysteresis effects. Consequently, the hysteresis effect was quite large for no. 5 but relatively smaller for no. 7. In turn, including the hysteresis in the model improved the model behavior substantially for point no. 5 as seen from Figure 4. The scanning curves in water retention curves at point no. 5 (sample b) and no. 7 (sample d) is plotted on Figure 2. The sample with lower permeability shows large hysteresis effects and thus the discrepancy in Figure 4a becomes obvious. Figure 5 shows numerically obtained discharge with respect to time from the bottom of the flow tank Q b after the 2nd ponding (t = 80-160 min). It is clear that the discharge was affected by the pore-air flow. For the twophase flow simulations (Run-2 and Run-2h), there was an immediate sharp increase (0-25 min) caused by the increase of pore-air pressure just after ponding. This was due to that the water content in the flow tank decreased and thus the discharge gradually increased. Although the peak of Run-2 (ignoring hysteresis) was larger than Run-2h (including hysteresis), the difference between these was quite small. Thus, the influence of hysteresis for global behavior is smaller than for local behavior in two-phase flow simulation. In the one-phase flow simulations (Run-1 and Run-1h), the initial large discharge was not observed. The initial large discharge was caused by initially included pore-air. When the water was ponded at the soil surface, pore-air was pushed down lower water as a piston flow. Only two-phase flow models can simulate this phenomenon. One-phase flow simulations displayed a single peak after 100 min. When including hysteresis (Run-1h), the onset of increasing discharge was earlier than for the case of ignoring hysteresis (Run-1). The peak of including hysteresis (Run-1h) was higher than ignoring hysteresis (Run-1). In the case of ignoring hysteresis, the wetting curve was used as water retention curve. However, the situation before ponding represented more or less the drying phase. Thus, for small grain sizes, saturation may be underestimated as compared to the drying curve (see Figure 2a). If so, the global water content was smaller than for the situation including hysteresis. As a result, the discharge peak was delayed and lowered as compared to including hysteresis. In the uniform parameter model simulation (Run-0 and Run-0h), similar discharge behavior was seen in the one-phase flow simulation. According to the above discussion, the two-phase flow simulation will give most realistic results. Moreover, it is important to consider the pore-air flow in the case of the pore-air entrapment by ponding as in the present experiment because including hysteresis clearly improved volumetric water content change as seen in Figure 4. Supplement Movie S1 shows results from the two-phase flow simulation (Run-2h), which included hysteresis effects and permeability heterogeneity. Blue arrows show water velocity and red arrows pore-air velocity vectors, respectively. As seen, immediately after ponding started, infiltration through the soil surface was large, and much pore-air discharged upward through the soil surface. However, the simulations showed that the spatial variation of pore-air flow during infiltration was quite complex. Also, some downward flow was indicated during infiltration. Due to the complex hydraulic conditions the direction of pore air flow changed several times during the experiments. The main path of pore-air was the 3rd column from the righthand side in the experimental sand box. This column contained large grain size. This shows that pore-air and water flow display a complicated relationship during stepwise ponding. The induced preferential flow of air and water was determined by the hydraulic properties of the unsaturated media and not by the application rate. In the two-phase flow simulation, we can observe hysteresis loops for each soil sample. The sample with lower permeability displays a large hysteresis loop. Consequently, the unsaturated flow in this soil was influenced by large hysteresis effect (Nakagawa and Saito, 2008).

CONCLUSION
To investigate the importance of hysteresis and heterogeneity in numerical modeling, laboratory experiments and numerical simulations were performed. The experiments included step-wise ponding of water and subsequent infiltration into heterogeneous porous media. The soil water content at 8 spatial points in the experimental flow tank was investigated using PCBP with TDR. To analyze detailed properties of water infiltration and effects of pore-air, two-phase flow and saturated-unsaturated seepage analyses were performed and compared. The importance of hysteresis in numerical simulation was investigated by numerical experiments. Moreover, to investigate effects of heterogeneity, simulation results from heterogeneous and homogeneous cases were compared. The following results were obtained. (1) By including hysteresis in simulations, numerical results were in agreement with observations. The results showed that it is especially important to include hysteresis when infiltration and drainage are repeated as in the present experiment. (2) The inclusion of heterogeneity is important in the numerical simulation to obtain more realistic results reflecting real field-scale conditions. This was shown by comparing heterogeneous and homogeneous simulation results. (3) Preferential flow of pore-air was investigated in the numerical simulation. The pore-air exhibited a rapidly changing pattern and a complex behavior when repeated cycles of water infiltration and drainage were performed. Figure 5. Discharge from the bottom of the flow tank