The unique experimental data of X-ray computed tomography and magnetic resonance were combined with three-dimensional numerical modeling to analyze the effect of air entrapment on steady state flow rate in two subsequent infiltration experiments. The change of the flow rate was explained by air entrapment in preferential flow regions.


Recurrent ponded infiltration experiments on undisturbed samples of coarse sandy loam have revealed a significant flow instability characterized by a decrease in the steady-state flow rate of the second infiltration run, conducted into wet soil, compared with the first infiltration run, conducted into drier soil. It has been hypothesized that this decrease was caused by air entrapment during the second run, with subsequent blocking of the preferential pathways. In this study, entrapped air distribution and its impact on water flow was studied through a novel combination of magnetic resonance (MR) imaging and numerical modeling. An undisturbed sample of coarse sandy loam was subject to recurrent ponded infiltration while being monitored by MR. Internal structure of the soil sample was visualized by x-ray computed tomography (CT). A parallel version of the three-dimensional water flow model based on Richards’ equation was used to simulate water flow through the heterogeneous soil sample. Information from the CT was used to describe internal heterogeneity of the soil sample via scaling factors. Magnetic resonance relaxometry imaging data were utilized to derive three-dimensional maps of entrapped air. These were implemented into a simulation of the second infiltration run as regions of no flow. The MR images were used to assess the water content distribution within the sample. The three-dimensional model was able to describe measured outflow rates and pressure heads and also to reproduce the heterogeneous distribution of water content within the sample. The results obtained support the assumption that the observed decrease in the outflow rate could be caused by entrapped air in large pores of the soil sample.

Current theory for mathematical predictions of water flow and solute transport in porous media assumes that the air phase plays a relatively minor role during unsaturated flow and transport (e.g., van Genuchten et al., 1999). There is a long-term evidence, however, of air-phase discontinuities occurring in porous media and leading either to water entrapment (Papafotiou et al., 2008) or to air entrapment effects on water flow (Christiansen, 1944; Faybishenko, 1995) and solute transport in soils (Geistlinger et al., 2005). For instance, entrapped air was found to be responsible for an increase in infiltration variability (Wang et al., 1998), temporal changes in hydraulic conductivity (Sakaguchi et al., 2005), and flow instabilities (Schultze et al., 1999; Wildenschild et al., 2001).

To study infiltration processes exhibiting air entrapment effects, Snehota et al. (2010) performed a recurrent ponded infiltration experiment on an undisturbed heterogeneous soil column. With the use of MR imaging, they detected significant air entrapment in large pores and associated this with the decrease in the quasi-steady-state outflow rate during infiltration into a wet sample. They presented only qualitative results because the MR pulse sequence used did not provide spatially resolved information about changes in water content.

Three-dimensional modeling of water flow through heterogeneous porous systems may serve as a valuable tool for detailed analysis of flow instabilities caused by entrapped air. The popularity of three-dimensional flow models has been increasing since parallel computing became available, which allows the introduction of a high degree of soil heterogeneity in the simulations (e.g., Hardelauf et al., 2007). In addition, the capability of one-dimensional models to reproduce flow in structured soils could be tested by means of three-dimensional modeling (Laloy et al., 2010). The high level of complexity of three-dimensional simulations is achieved by either introducing scaling factors to interpret the measured hydraulic properties (Javaux and Vanclooster, 2006) or using information obtained from MR imaging (Yoon et al., 2008).

A noninvasive visualization by means of MR imaging has proven to be useful in three-dimensional investigation of both water flow and solute transport in porous media (e.g., Deurer et al., 2004; Pohlmeier et al., 2009; Jelinkova et al., 2010). Votrubova et al. (2003) reported a decrease in the steady-state infiltration rate between two successive ponded infiltrations and used MR imaging to quantify changes in the water distribution within a sample. The internal structure of a soil sample is usually obtained using x-ray CT (Hopmans et al., 1994), which is another noninvasive visualization method. Cislerova and Votrubova (2002) concluded that CT provides information at the macroscopic level and therefore can be used for describing porous medium heterogeneity. Kasteel et al. (2000) used CT images to separate dense and less dense regions of the soil matrix; each of the two materials was then introduced into a three-dimensional numerical model with specific soil hydraulic characteristics.

The objectives of the present study were (i) to use the information from CT imaging to simulate transient water flow in a three-dimensional heterogeneous porous system, (ii) to compare the results of the simulations with both the infiltration experiment results and the MR imaging signal intensities, and (iii) to verify the assumption that the flow instabilities observed were caused by the air entrapped in the soil sample. More specifically, we analyzed the hypothesis that MR relaxometry imaging can be used to derive information about the volume and distribution of entrapped air and that this air entrapped in well-connected preferential pathways reduces the overall flow rate more effectively than the same volume of entrapped air localized in the soil matrix.

Background Information

Noninvasive Imaging Techniques

Noninvasive visualization techniques: (i) x-ray CT, (ii) MR imaging, and (iii) MR relaxometry imaging were used in this study. Brief background information about these techniques is presented below.

X-ray Computed Tomography

The internal pore structure of soil samples can be visualized by means of x-ray CT. The CT images are obtained using x-ray CT scanners. Then the three-dimensional scan image is mathematically reconstructed from acquired scans of consecutive horizontal slices (XY planes, perpendicular to vertical flow line). The scanner resolution in the XY plane and the slice thickness define the voxel (three-dimensional pixel) size. The intensity of the transmitted x-ray beam is usually expressed in Hounsfield units (HU), ranging between two extremes that represent either empty voxels or the most dense parts of the sample. Combined information about the pore size distribution of the soil studied, the voxel size, and x-ray imaging setup is used to determine the CT imaging outcome. In this study, the CT imaging provided macroscopic information predominantly reflecting the local bulk density distribution (Anderson et al., 1988). The water present in the soil during the CT imaging will cause a slight overestimation of the bulk density. It is possible, however, to make a reasonable assumption that primarily the low-conductivity regions (small pores) will be affected because the pores in highly conductive regions are empty due to their low air-entry value. Thus the CT imaging has the potential to provide information on the distribution of the low-density regions, i.e., large pores. These are the regions that form structural heterogeneity or potential preferential pathways.

Magnetic Resonance Method

The sample is placed in a large homogeneous static magnetic field and exposed to a series of radio frequency pulses. Subsequently, emitted radiation from H nuclei is detected by a receiver coil. The induced voltage defines the MR signal intensity I. Under idealized conditions, MR imaging directly refers to the spatial detection of water via measured signal intensity.

The decay of MR signals with time is a function of the longitudinal relaxation time T1 and the transverse relaxation time T2. Spatial resolution of relaxation times is subject of MR relaxometry imaging (e.g., a series of imaging using the spin-echo sequence). Determination of T1, T2, and I0 values is then achieved by fitting the acquired data series by the single exponential approximation of the decay function 
where t is time after the pulse, I0 is the spatially resolved amplitude value of signal intensity I, TR and TE are the repetition and echo times, respectively, of the applied sequence.

In general, MR relaxation in porous systems is enhanced by the presence of the air–solid–liquid interfaces. The factors controlling the relaxation process of water in a porous medium are the interface geometry (Brownstein and Tarr, 1979), the amount of the paramagnetic particles in contact with water (e.g., Keating and Knight, 2007), and differences in magnetic susceptibility at solid–liquid interfaces (e.g., Stingaciu et al., 2009). When dealing with MR imaging methods in natural structured soils these mentioned factors are combined in a complex way, which may result in nonvisibility of water in small pores (Cislerova et al., 2002). Furthermore, the relationships between the soil water content and MR-detected water content can be nonlinear and dependent on MR protocol and soil properties (Votrubova et al., 2000).

Modeling Approach

Governing Equation

Water flow through a three-dimensional variably saturated porous medium is described by Richards’ equation: 
where h is the pressure head (m), C is the specific water capacity (m−1), K is tensor of the hydraulic conductivity (m s−1), z is the vertical coordinate (m), and t is time (s). The hydraulic conductivity tensor was assumed to be isotropic and hence fully determined by a scalar function.

Hydraulic Characteristics

Hydraulic properties are described by a set of closed-form analytical expressions similar to those of van Genuchten (1980), who used the statistical pore-size distribution model of Mualem (1976) to determine the unsaturated hydraulic conductivity function. The parameterization was modified by Vogel and Cislerova (1988) and further extended by Vogel et al. (2000) to improve the description of the soil hydraulic characteristics near saturation. In the modified approach, the soil water retention function, θ(h), is expressed as 
where θr is the residual water content (m3 m−3), θm is the fictitious (extrapolated) parameter, and nVG (dimensionless) and α (m−1) are empirical shape parameters, respectively. The parameter m (dimensionless) is equal to 1 − 1/nVG and nVG is >1. This function is used to describe retention curves in the unsaturated pressure head range h ∈ (−∞, hs). For pressure heads larger than the “air-entry value” hs (m), the corresponding water contents become equal to the saturated water content θs. Combining Eq. [3] with Mualem’s model yields the unsaturated hydraulic conductivity function, K(θ): 
where Ks is the saturated hydraulic conductivity. The parameter θm is related to the air-entry value hs through the relationship 

Materials and Methods

Soil Sample

The soil under study is a coarse sandy loam from the experimental site Korkusova Hut (Sumava Mountain, Czech Republic); the soil is classified as a Dystric Cambisol. The internal structure of the soil is well developed with a broad range of pore sizes. Significant preferential flow effects were reported by Cislerova et al. (1990, 2002) for the same soil. The preferential flow can be attributed to conductive pathways along biopores and to the soil structure as well as to the spatial variability of local hydraulic properties (Cislerova et al., 1988). Three undisturbed soil columns were taken from the depth of 40 cm (B soil horizon). The sample diameter was 5.5 cm and the height was 11.2 cm.

Soil Sample Imaging

The CT images were obtained using a Siemens Somatom Definition medical scanner (voltage 140 kV and current 36 mA). The scanner resolution in the XY plane was 0.2 by 0.2 mm; the slice thickness was 0.6 mm. The CT imaging was performed on undisturbed soil columns. The sample showing the least damage to the internal structure was selected for further testing. Moreover, the information acquired prior to and after the water flow experiment was used to rule out pore structure changes. Local bulk densities derived from CT imaging before the experiment were implemented to the model.

The magnetic resonance measurements were performed on 4.7-T vertical bore (310-mm internal diameter) superconductive magnet system with an 80-mm i.d. radiofrequency coil tuned to 1H resonance frequency (Varian Inc.). Magnetic resonance imaging (MRI) was used to assess changes in the distribution of water within the sample during infiltration. Magnetic resonance relaxometry imaging (MRRI) was performed to obtain a three-dimensional map of the entrapped air in the soil sample. Details of the two MR methods used during the infiltration experiment is shown in Fig. 1. For both, the multi-echo–multi-slice sequence proposed by Edzes et al. (1998) was applied. Sequences mainly differed in the number of echoes, acquisition time, and orientation. Note that a single echo was applied in case of MRI. Detailed information on the MR setting is presented in Table 1. The MR sequence setting used (parameters TE and TR) was a compromise between time-scale resolution and soil pore visibility. While fitting Eq. [1] to the I(t,x,y,z) data in the case of MRRI, yielding spatially resolved values of I0, T1, and T2, only amplitude values I0 were used in this study. Derived I0 values associated with the respective performance measure of fitting provide more robust information than the signal intensity available from MRI. The MR measurement on the sample of interest was described in detail by Jelinkova et al. (2011).

Water Flow Experiments

The recurrent ponded infiltration–outflow experiment was performed while the sample was placed in the center of the MR scanner. The experiment included two consecutive runs of the ponded infiltration separated by a redistribution phase allowing gravity drainage (see Fig. 1). During each infiltration run, the surface of the sample was flooded by water, and a constant ponding depth of 2 cm was maintained at the top of the sample by a Marriotte bottle located outside the MR scanner. During the redistribution, the top of the sample was exposed to atmospheric pressure. At the bottom of the sample, the water flowed freely out of the sample through a perforated support plate.

The soil water pressure head at the depth of 8.8 cm below the sample surface was monitored using a modified UMS T5 tensiometer (UMS GmbH). Similarly to the experiment conducted by Snehota et al. (2010), the tensiometer body was substantially extended using nylon tubing so that the pressure transducer could be placed outside the MR scanner to prevent electrical interference between the pressure transducer circuit and the scanner. Unfortunately, inadvertent elasticity of the tensiometer assembly slowed down the tensiometer reaction. The outflow from the bottom of the sample was collected by a fraction collector located outside of the MR scanner. Similarly to the pressure transducer, the long tubing to the fraction collector led to a delay in the outflow record. Both delays were experimentally determined so the pressure head and outflow data could be corrected (the delay in outflow data measurement was constant during the experiment; the delay caused by elasticity of the tensiometer assembly was pressure-change dependent).

The first infiltration was performed into relatively dry soil with good connectivity of the air phase (a soil water pressure head of −25 cm before the start of the infiltration experiment). Thus when the infiltration was started, the air could escape freely from the bottom of the sample. Therefore, the first infiltration is supposed to be close to the conditions when the soil sample is without entrapped air. For the second infiltration run, the initial condition was noticeably different. The tensiometer showed a pressure head value of −2 cm, which means that the bottom of the sample was close to saturation and the air phase could easily become entrapped between the wetting front and the saturated bottom of the sample. Malfunction of the tensiometer resulted in an incomplete data series for the second infiltration run.

Water Flow Model

The governing Eq. [2] was solved numerically by the S3D model using the finite-element method on a general tetrahedral mesh. The S3D code is a follow-up version of the SWMS 3D code (Vogel, documentation of SWMS 3D code, CTU–Prague, internal report, 1993) developed jointly at the Czech Technical University in Prague and the Technical University of Liberec. The new version of the code uses the Portable, Extensible Toolkit for Scientific Computation (PETSc) library (Balay et al., 2008) for parallel solution of the nonlinear algebraic system that is solved at each time level. In particular, we used the inexact Newton method, the one-level additive Schwartz method for domain decomposition, and conjugate gradients preconditioned by incomplete LU factorization for the inner solver on each decomposed domain.

For time discretization, we used a backward Euler method with adaptive time stepping (Rathfelder and Abriola, 1994) and an error accounting method to achieve an accurate water balance at each time step. We treated the residual of the nonlinear solver as a water source for the next time step or across the time interval comparable to the actual time step.

Implementation of Computed Tomography Information

The original three-dimensional reconstruction of the CT horizontal slices with a resolution of 0.2 by 0.2 by 0.6 mm was resampled using simple averaging of neighboring values. This procedure was performed to obtain isotropic voxels of 0.6 by 0.6 by 0.6 mm and to reduce computational demands (hereafter voxels are considered as nodes for numerical modeling). In total, computational nodes amounted to 1,100,364.

Nodal porosities were calculated assuming a linear relationship between the measured HU values and the sample bulk densities (Hunt et al., 1988), similarly to Hopmans et al. (1994). Thus, the nodal porosity equaled 1 for voids, whereas the densest part of the sample corresponded to a nodal porosity of 0. Calculated porosities followed a nearly normal (Gaussian) distribution.

A scaling technique was used to introduce the spatial variability of the soil hydraulic properties. The scaling procedure was based on the concept of linearly variable hydraulic properties proposed by Vogel et al. (1991). According to this concept, the spatial distribution of soil hydraulic characteristics can be approximated using the set of dimensionless scaling factors of hydraulic conductivity αK, water content αθ, and pressure head αh. In this study, we assumed that the variability of hydraulic properties can be expressed based on the x-ray CT results using two scaling factors: 
where ni is porosity in the individual node (m3 m−3), forumla is the arithmetic mean of porosities in the whole sample, and Ksi and forumla denote the individual nodal saturated hydraulic conductivity and the sample mean value of saturated hydraulic conductivity (m s−1), respectively. The individual nodal saturated hydraulic conductivity Ksi could be defined by the modified Kozeny–Carman formula (Bear and Verruijt, 1987): 
where μ is the viscosity of water (kg s−1 m−1) and c is a constant depending on the pore geometry. Calculated conductivities were approximated by a lognormal distribution. The mean hydraulic conductivity could then be defined: 
where N is the total number of nodes. The value of forumla corresponds to the saturated hydraulic conductivity measured at the representative elementary volume. The mean values of both scaling factors αK and αθ are equal to 1.

Calculated values of αK used in the model were limited to the interval 0.001 to 43.75. Corresponding saturated hydraulic conductivity for stones with a low value of voxel porosity and for highly conductive pores were 0.0045 and 197 cm h−1, respectively. Inspection of the Kozeny–Carman formula reveals that the porosities >0.849 cm3 cm−3 were considered by the same value of αK (equal to 43.75). In analogy, the porosities <0.08 cm3 cm−3 were considered by the single scaling factor αK equal to 0.001. Similarly, the minimum value of αθ was limited to 0.001; the maximum value of αθ (=1.92) was kept unrestricted. The reason for the restrictions of the αK and αθ values was twofold: (i) to maintain laminar flow in tortuous soil pores; and (ii) to ensure convergence of the numerical model.

Note that the upper limit of αK values together with the reference value of Ks are two model parameters that control the maximum water flow velocity and flow field heterogeneity simulated within the sample. The resulting spatial distribution of scaling factors in the soil sample is shown in Fig. 2. Locations with low values of αθ (i.e., stones) are shown in Fig. 2a, while highly conductive pathways characterized by high values of αK are shown in Fig. 2b.

The CT information was also used for setting computational nodes inside the tensiometer body and ceramics as a no-flow region. Several computational nodes in the vicinity of the tensiometer ceramics were set as observation points to allow comparison with pressure head values monitored by the tensiometer.

Implementation of Magnetic Resonance Information

The MRRI was used to obtain information about the volume and distribution of the entrapped air phase. Due to extremely fast MR relaxation, the soil water in small pores produces unreliable MR signals (Votrubova et al., 2003); however, “free” water in large pores can be detected (Hall et al., 1997). At the end of redistribution (before the second infiltration run), pressure heads ranged from 0 cm at the bottom to −11.2 cm at the top of the sample. Maximum curvature of menisci calculated at the air–water interface from these pressure values implies that water, which drained between two different MR measurements, drained from pores with a radius >300 μm. Therefore, water content changes took place in large pores containing “MR visible” water.

The change in water volume between the first and second infiltration run, δVMR (m3), was calculated from I0 maps of the individual stages of the experiment (Jelinkova et al., 2011): 
where forumla and forumla are the sums of I0 values for entire sample volume during the steady-state flow of the first and second infiltration runs, respectively; forumla is the sum of I0 values for the entire volume of the sample measured during the final equilibrium state of the redistribution phase; and Vgdw is the measured volume of gravitationally drained water during the redistribution phase (Vgdw = 7.75 mL). Jelinkova et al. (2011) reported that the change in “MR visible” soil water content between the steady state of the first and the second infiltration runs was equal to 2.2% (δVMR divided by the volume of the soil sample). This volumetric water content decrease was assumed to correspond to the volume of entrapped air in the soil sample during the second infiltration run.
To map the spatial locations where the water content decrease took place, I0diff values for each voxel were calculated as 

The detection limit of I0 mapping was estimated to be 0.19 (based on coupled noise filtering and the performance measure assessment of Eq. [1]). Note that for calculating I0 differences, only the relevant values of I0 (based on the detection limit) were considered. The histogram of I0diff in the whole soil sample, shown in Fig. 3, illustrates that the number of highly positive I0diff values predominates over negative ones. Maps of entrapped air filling 1.1, 2.2, and 4.4% of voxels with the greatest positive values of I0diff were considered (note that the limit used for the entrapped air mapping of 2.2% was almost 50 times higher than the estimated detection limit).

Three different volumes of air entrapment were applied in the analysis to account for a relatively high uncertainty level associated with the MR measurement (see below). For instance, MR-visible changes in water content could correspond to voxels not fully drained for scenarios based on 1.1% and 2.2% of voxels. On the contrary, the scenario with 4.4% of voxels represent the possibility that the porosity of the voxels with entrapped air can be <1. The corresponding entrapped air volumes for the considered voxel percentages are presented in Table 2. To apply the detected entrapped air structures in the numerical model, their distribution map with the original MR image resolution of 1 by 1 by 5 mm was resampled using a bilinear interpolation method to match the resolution of the computational mesh (i.e., 0.6 by 0.6 by 0.6 mm). The entrapped air structures were then introduced to the computational mesh as no-flow regions. The map of no-flow regions representing the entrapped air as detected by MRRI is shown in the Fig. 4.

Three-Dimensional Numerical Realizations of Air Entrapment Structures

To examine the effect of varying the entrapped air distribution on the water flow, 30 different numerically generated realizations of the three-dimensional spatial arrangement of air entrapped within the sample were considered. The sequential Gaussian conditioned kriging approach implemented in the SGSIM module of the GSLIB program (Deutsch and Journel, 1992) was used to generate these realizations of the spatial distribution of the I0 differences. The spatial correlation structure model of the I0 differences was based on the MR relaxometry imaging; the experimental variograms suggested a correlation length of 5 to 10 mm for the horizontal plane and 5 mm for the vertical plane. The generated three-dimensional I0diff fields were then sampled to obtain the MR-detected amount of entrapped air voxels (2.2%). Note that the obtained three-dimensional I0diff fields were fully independent from the CT images.

Model Application

The numerical simulation of soil water flow was performed for the whole experiment—the first infiltration run together with the redistribution phase and the second infiltration run (the simulation was interrupted between runs to include the air entrapment structures). The initial pressure head distribution was derived from the tensiometric observation at the start of the first infiltration run. For the second infiltration run, the pressure head distribution at the end of the redistribution phase was used. The upper boundary of the soil column was exposed to: (i) ponded infiltration using a Dirichlet-type boundary condition with a 2-cm depth of ponding, and (ii) a zero flux Neumann-type condition during the redistribution phase. The lower boundary condition was the seepage face, which assumes zero flux when the bottom part of the sample is unsaturated and allows water to drain from the sample when the bottom layer becomes saturated.

The main drying and wetting retention curve branches determined by standard laboratory equilibrium techniques on 100-cm3 undisturbed soil samples were very similar except near saturation. In fact, the natural soil under study proved that the uncertainty associated with the process of retention curve measurement is usually greater than the difference determined between the main retention curve branches (e.g., Kodesova and Gribb, 2006). Therefore, the hysteretic behavior in the simulated case was associated exclusively with the inclusion of entrapped air structures into the second infiltration run.

The saturated hydraulic conductivity values determined on 1000-cm3 soil cores and by two independent infiltration tests performed on the sample under study (before the infiltration experiment was simulated) were similar (9.0 cm h−1 for 1000-cm3 soil cores; 8.2 and 5.5 cm h−1 for the soil sample under study). The primary estimate of the saturated hydraulic conductivity was then adjusted so that the three-dimensional simulation with the heterogeneity introduced via a scaling factor of hydraulic conductivity (see above) matched the measured “steady-state” outflow rate during the first infiltration run. The soil hydraulic parameters used for the simulations were θr = 0.03 cm3 cm−3, θs = 0.52 cm3 cm−3, α = 0.11 cm−1, nVG = 1.21, Ks = 4.5 cm h−1, and hs = −0.5 cm.

The CT-based scaling factors of hydraulic conductivity (αK) and water content (αθ) were used in the simulations of both runs. The air entrapment at the beginning of the first infiltration run could not be determined using Eq. [11] because entrapped air is associated with the change in “MR-visible” water content between two successive infiltration runs. The presence of entrapped air was thus neglected during the first infiltration run. The simulation of the second infiltration run made use of the inclusion of air entrapment structures (either measured or geostatistically generated) as no-flow regions. Note that model parameters were not changed between the first and second infiltration runs.


Simulation Results Comparison (First Infiltration Run)

Figure 5 depicts the measured and simulated outflow rate as a function of time for the first infiltration run and a part of the redistribution phase. The outflow from the sample occurred approximately 190 s after the beginning of ponding. The onset of measured outflow was captured reasonably well by the model (see inset in Fig. 5). The measured outflow rate initially reached 12 cm h−1, then it gradually decreased to 9 cm h−1. While the gradual decrease in the outflow rate could not be reproduced by the model based on Richards’ equation, the later phase of quasi-steady state was predicted appropriately. Note that the simulated steady-state outflow rate remained higher than Ks. This was caused by the sample heterogeneity implemented via the scaling factors of hydraulic conductivity (more specifically by the spatial arrangement of the local Ks values). When the ponding ceased at 20,445 s, the redistribution phase of the experiment started. Water continued draining from the bottom of the sample until 22,165 s. Figure 5 shows a good agreement between the simulated and measured outflow rates during this phase of the experiment. The redistribution of soil water content within the sample continued for the next 23 h.

The pressure head data from the tensiometer indicated that the wetting front reached the tensiometer shortly after the ponding started (Fig. 6). The first contact of the tensiometer with the wetting front was recorded at t = 120 s. Figure 6 further shows the simulated pressure heads for different observation points in the vicinity of the tensiometer ceramics. The figure illustrates a relatively high variability of simulated pressure heads near the tensiometer. The first contact of the tensiometer with the infiltrating water was predicted in the range between 140 and 175 s after the ponding started.

Overall, the outflow rates in Fig. 5 as well as the pressure heads in Fig. 6 indicate that the simulation results of the first infiltration run exhibited a high level of coherency with the experimental data.

The MR images taken during the transient phase of the infiltration experiment show distribution of “free” (MR-detectable) water within the sample. These images were used to assess the general agreement between the measured and predicted wetting front advancements in the sample. First, the MR signal intensity differences (Idiff) were calculated between each actual MR image and the one taken at the initial stage of the infiltration (t = 0 s). The same procedure was repeated for the simulated water content to obtain the difference Δθ. The differences (Idiff and Δθ) for two selected time levels are shown in Fig. 7. The wetting front detected by MR imaging showed more heterogeneous patterns than simulated one. The simulation did reproduce reasonably well the wetting of high-porosity regions. This can be illustrated by MR-detected wetting of these regions, which are indicated by arrows in Fig. 7.

Simulation Results Comparison (Second Infiltration Run)

During the second infiltration run, the measured steady-state outflow rate dropped to 3 cm h−1 (Fig. 8). The outflow from the sample occurred 169 s after the beginning of the second run. Figure 8 clearly shows that entrapped air introduced in the simulation of the second infiltration run significantly decreased the steady-state outflow rate. Compared with simulation without entrapped air, an outflow rate decrease of 0.7 cm h−1 was obtained for the 1.1% voxels of introduced entrapped air, while the inclusion of 4.4% led to an outflow rate decrease of 3.0 cm h−1. Nevertheless, none of the four scenarios shown in Fig. 8 was able to decrease the outflow rate to the measured values. Possible reasons are discussed below. The simulated onset of outflow occurred 110 s earlier than measured. The difference in timing of the beginning of outflow simulated based on different air entrapment percentages was negligible (<5 s).

Numerical Analyses

Sensitivity of Volume of Entrapped Air

Analysis of the gradual increase in the entrapped air volume introduced into the simulations revealed that the same small increment of entrapped air voxels had a different impact on the outflow rate decrease (Fig. 9). The figure was constructed by gradually varying the air entrapment (by 0.55%) in the sample while registering the decrease in the steady-state outflow rate. For instance, the increase of air entrapment from 1.1 to 1.65% would lead to a decrease of 0.32 cm h−1 in the steady-state outflow rate (Fig. 9). Note that different amounts of entrapped air voxels were implemented in the simulations based on MR measurements, i.e., entrapped air structures spatially followed the measured MR information. The impact was relatively small (<0.36 cm h−1) for the amount smaller than about 2% because for these scenarios preferential pathways were probably not completely blocked by entrapped air. By increasing the air entrapment, the effect became significant as preferential pathways were substantially or fully blocked. The largest decrease in the steady-state outflow rate was obtained for about 3.8% of entrapped air (see Fig. 9). A larger volume of entrapped air had a less noticeable impact on the outflow rate decrease. This could be caused by the fact that air entrapment grew outside of the preferential pathways, thus starting to block also the surrounding soil matrix, which had a small impact on the rate of outflow from the sample.

Projection of Air Entrapment Structures into the Soil Hydraulic Properties

The results of this analysis demonstrate that the MR-derived sample regions potentially blocked by entrapped air had considerably higher local porosity and conductivity than the rest of the sample (Table 2). Considering the entrapped air structures with total voxels of 2.2%, the porosity of the related locations was 0.64 cm3 cm−3, while the average porosity of the whole sample was 0.52 cm3 cm−3. Similarly, the saturated hydraulic conductivity for these structures was three times higher than the reference value of the sample. The reason for the decrease in porosity and conductivity values for the scenario based on 4.4% of entrapped air (Table 2) was probably similar to that stated above: the growth of air structures outside the preferential pathways (to the soil matrix). Note that the interconnection of high-porosity regions, which seemed to have a crucial impact on the water flow regime, could be assessed by three-dimensional model only.

Sensitivity to Different Locations of Entrapped Air: Geostatistical Simulations

The average steady-state outflow rate given by 30 simulations based on different spatial arrangements of entrapped air was 8.49 cm h−1 (minimum 8.01 cm h−1, maximum 8.77 cm h−1). The average value corresponded to 111% of the steady-state outflow rate predicted based on MR-detected air entrapment of 2.2%. Moreover, the majority of these realizations had steady-state outflow rates higher even than the simulation with MR-detected air entrapment voxels of 1.1% (8.24 cm h−1). The difference in the outflow timing of these 30 simulations was negligible (57–65 s).

Degree of Heterogeneity

Initially, we used αK = 43.75 as the upper limit in the simulations. The decrease in the outflow rate between the first and second infiltration runs for the simulation based on 2.2% of entrapped air was about 26%. As the MR-derived air entrapment suggested blocking within high-porosity regions with high αK (controlled by the αK limit), the simulated decrease in the outflow rate was affected by the αK limit used. To study this effect, several simulations using less restricted values of αK were additionally tested in which the upper limit of αK was relaxed up to a value of 700. For a limited number of nodes having such high αK values, it was expected that laminar flow in soil pores would still be maintained. The decrease in the outflow rate between the first and second infiltration runs for the simulation with an upper limit of αK = 700 was about 50%.


In the present study, air entrapment was included in the single-phase model as no-flow regions. This simplified approach is intended neither to simulate the actual mechanism of air trapping nor to capture complex interactions between the entrapped air phase, the soil matrix, and soil water (e.g., air dissolution in the water phase, dissolved gas transport, etc.). Nevertheless, we believe that the approach principally allows the verification of the initial hypothesis about the effect of air entrapment on outflow rates. As far as can be ascertained, the approach combining CT and MR measurements with three-dimensional numerical modeling seems to be unique and promising with respect to analyzing transport processes in heterogeneous porous systems. A two-phase flow model with complex air-phase parameterization may serve as an alternative to the approach used in this study. The objective of this analysis, however, was to evaluate the decrease in the outflow intensity caused by the presence of a discontinuous air phase at specific MR-detected regions.

When interpreting differences between the MR-detected and simulated wetting fronts, the inherent difficulties of MR imaging in natural structured soils should be kept in mind: water in small pores is not visible (Cislerova et al., 2002), and even though the MR acquisition time of 19 s seems to be relatively short, it is still long considering the velocity of the wetting front advancement (see Fig. 7, t = 10 s). The correction of the possible artifacts caused by rapid water content change was not feasible; however, it is reasonable to assume that the upper part of the sample was captured more reliably then the voxels near the wetting-front edge.

Flow simulations for different entrapped air arrangements based on geostatistical simulations suggested that the generated entrapped air structures were distributed in the soil matrix, hence not effectively acting as blockage of larger pores with high values of αK. This further implies that the spatial distribution of air entrapment structures played a crucial role in predicting the steady-state outflow rate.

Although the volume of entrapped air affected the simulated steady-state flow rate significantly, the time of outflow occurrence from the bottom of the sample remained almost the same for all considered scenarios (Fig. 8). This was caused by both the seepage face boundary condition on the bottom and the volume of the soil sample. The simulation results suggested that outflow starts when the whole sample is close to its full saturation. The decrease in the sample storage capacity (caused by an increase in the entrapped air volume) led to earlier outflow from the sample. More pronounced differences in outflow timing for the scenarios were suppressed by the relatively low volume of entrapped air in the soil sample.

Because only water in the large pores could be detected by MRRI, possible air entrapment in the smaller pores was not included in the model as no-flow regions. Therefore, the simulated decrease in the steady-state flow rate of the second infiltration run could resemble the measured decrease better if the entrapped air in these pores was considered appropriately. The entrapped air detection might possibly be improved by the use of advanced MR techniques such as single-point ramped imaging with T1 enhancement (SPRITE) (Pohlmeier et al., 2008), which allows monitoring of water in small soil pores, or by the use of other high-resolution noninvasive techniques such as neutron imaging or high-energy x-ray CT (e.g., Oswald et al., 2008; Wildenschild et al., 2005).


The presented experimental data include a unique data set collected during concurrent MR measurement and an infiltration–outflow experiment on an undisturbed heterogeneous soil sample of submacroscopic scale.

A parallel version of three-dimensional numerical code was used to simulate water flow in the heterogeneous soil sample. The soil heterogeneity detected by CT imaging was implemented into the model by means of scaling factors of hydraulic conductivity and water content for each individual computational node in the three-dimensional domain. The simulation results were compared with the data from the infiltration experiment and signal intensities from MR imaging.

Implementation of MR information into the simulations supported the initial hypothesis that the measured flow instability (the decrease in the outflow rate between two successive infiltration runs) could be caused by a discontinuous air phase entrapped in large pores of the soil sample. Through this numerical analysis, air entrapment was found to have a negligible effect on the beginning of the outflow from the soil sample.

Analysis of simulations based on geostatistically generated air entrapment structures proved that MR-detected entrapped air was located at positions with a significant effect on water flow (e.g., preferential pathways). Moreover, our results indicated that the spatial distribution of entrapped air within the soil sample had a greater impact than its overall volume.

The research was supported by the Czech Science Foundation, Projects no. 103/08/1552 and 205/09/P567. The comments of Milena Cislerova, Jana Votrubova, and Tomas Vogel (CTU in Prague) are gratefully acknowledged. We are also grateful to Andreas Pohlmeier and Dagmar van Dusschoten, who made possible the magnetic resonance measurements in Forschungszentrum Jülich. The Supercomputing Center of the CTU in Prague deserves our special thanks for the possibility of running the simulations on a supercomputer network.

All rights reserved. No part of this periodical may be reproduced or transmitted in any form or by any means, electronic or mechanical, including photocopying, recording, or any information storage and retrieval system, without permission in writing from the publisher.
Open access article.