The Impact of Capillary Trapping of Air on Satiated Hydraulic Conductivity of Sands Interpreted by X-ray Microtomography

: The relationship between entrapped air content and the corresponding hydraulic conductivity was investigated experimentally for two coarse sands. Two packed samples of 5 cm height were prepared for each sand. Air entrapment was created by repeated inﬁltration and drainage cycles. The value of K was determined using repetitive falling-head inﬁltration experiments, which were evaluated using Darcy’s law. The entrapped air content was determined gravimetrically after each inﬁltration run. The amount and distribution of air bubbles were quantiﬁed by micro-computed X-ray tomography (CT) for selected runs. The obtained relationship between entrapped air content and satiated hydraulic conductivity agreed well with Faybishenko’s (1995) formula. CT imaging revealed that entrapped air contents and bubbles sizes were increasing with the height of the sample. It was found that the size of the air bubbles and clusters increased with each experimental cycle. The relationship between initial and residual gas saturation was successfully ﬁtted with a linear model. The combination of X-ray computed tomography and inﬁltration experiments has a large potential to explore the e ﬀ ects of entrapped air on water ﬂow.


Introduction
Interaction between water and gas needs to be considered for water fluxes in near-saturated and saturated soil because the entrapped air may block pores and thereby affect the hydraulic conductivity [1][2][3]. Entrapped air is defined as a discontinuous non-wetting phase that is not connected to the atmosphere [1]. Entrapped air may appear in the form of single pore bubbles, clusters or ganglia [4]. With growing size, such ganglia become more mobile [5] and may escape from the soil to the atmosphere. The pressure in the discontinuous air phase usually differs from the atmospheric pressure [6] and can be measured from the curvature of the water-air interface.
We distinguish between three states of saturation. If the pressure in the soil water is lower than the air entry value, the soil is denoted as unsaturated soil. In this state, it contains a continuous air phase with connection to and in equilibrium with the atmosphere. Quasi-saturated [1] or satiated [7] soil contains air in the form of entrapped air bubbles or air ganglia with water pressures higher than the air entry value [6]. Fully saturated soil, which is rarely found under natural conditions, does not contain any air bubbles or ganglia. The amount of entrapped air in porous media is commonly expressed as a residual gas saturation, S g,r (-), or an entrapped volumetric air content [1], ω (-). The residual gas Snehota et al. [3] illustrated and quantified air redistributions during infiltration experiments in a repacked sand sample using neutron imaging. The sand sample contained blocks of fine sand surrounded by coarse sand. Snehota et al. [3] found that entrapped air was transferred from the fine into the coarse sand while the flow rate decreased from 0.400 cm·min −1 to 0.296 cm·min −1 . Similar results were obtained in an experiment with a similarly designed sand sample [2] on which repeated infiltration runs were performed. Also here, the redistribution of entrapped air from the finer to the coarser sand caused the infiltration rate to slow down. The flow rate decreased from 3.7 cm·h −1 to 1.0 cm·h −1 during the first infiltration run and during the second infiltration run from 3.6 to 1.6 cm·h −1 . Thereafter, an infiltration run with degassed water was performed during which the flow rate increased to a maximum value of 10.5 cm·h −1 .
Hydraulic conductivity is a crucial parameter in all components of water cycle modeling, rainfall excess determination, water and solute transport in vadose zone as well as in recharge of groundwater [26,[45][46][47]. The impact of entrapped air on measured values of satiated hydraulic conductivity means that instead of a well-defined constant parameter, its value may vary in dependence on many factors, some of them described in this paper.
This study aimed to unravel the relationship between the entrapped air content and the dependence of the air trapping on the initial air saturation in two coarse sands. Recent experimental studies [3,48] on air entrapment involved the same sands as components of more complex, heterogeneous samples. To enable numerical modeling of these experiments in future, the K(ω) relationship of coarse sands is needed. An additional goal of the current study was to investigate variations of the entrapped air bubbles in space and time during repeated infiltration drainage cycles. Experiments were conducted on packed samples of coarse sand representing homogeneous media with narrow pore size distribution to rule out the effects of heterogeneity. The aim was to use X-ray computed micro-tomography to characterize the sizes, shapes and abundance of entrapped air bubbles.

Material Characterization
We prepared four packed samples, each two with different sands that we denote with abbreviations ST and FH, respectively. The ST sand is a commercially available technical quartz sand-with a SiO 2 content 99.40% and Fe 2 O 3 content 0.04% (technical designation ST 03/08, Sklopisek Strelec, Czech Republic). The FH sand is a commercially produced quartz sand from with a SiO 2 content of larger than 99% (technical designation FH31, Quarzwerke Frechen, Germany). We evaluated the grain-size distribution of the two different sands by dry sieving. Here we employed mesh-sizes of 0.20, 0.25, 0.40, 0.50, 0.63, 0.80, and 1.00 mm. The particle densities were measured using a pycnometer and turned out to be identical 2.62 g·cm −3 for all samples. The dry bulk densities of the samples were calculated from the mass and volume of dry sand samples. The samples' porosities were determined from the particle densities and the dry bulk densities. The grain size distributions of all samples are shown in Figure 1. Table 1 contains the respective values of dry bulk density, porosity as well as packing height and sample volume. The FH sand is finer and has a lower porosity than the ST sand. The value of the uniformity coefficient (D 60 /D 10 , where D x represents the size for which x% particles are smaller) for ST sand was 2.2 and for FH sand the value was 1.7. Both sands are considered to be poorly graded, but the FH sand is more uniform. The effective size (D 10 ) of the ST sand is 0.29 mm and of the FH sand is 0.22 mm.

Experimental Setup
The saturated and satiated hydraulic conductivity of sand were measured with repetitive fallinghead infiltration experiments. The respective setup is presented in Figure 2. The sample container consisted of an acrylic glass (PMMA) cylinder (7.2 cm inner diameter, 7.0 cm height) and a PMMA ring with a filter-mesh (size 99 µm) attached to its bottom. An o-ring made of nitrile rubber was used as a hydraulic seal, ensuring no water could escape between the ring and cylinder. The sample container was filled with the sand up to the height of approximately 5.0 cm. The container was placed into a larger, outer PMMA vessel (inner diameter 12.3 cm), where it was mounted onto five PMMA retainers. An overflow was inserted into the wall of the outer vessel (WADI-PZ, M12×1.5, 3.0-6.5 mm, Hugro Plastic Cable Gland, Waldkirch, Germany) through which a constant water level could be kept. To achieve a sufficiently large flow rate through the overflow, we attached a tube connected to a peristaltic pump (see Figure 2). The vessel was placed on a balance with 0.1 g precision to monitor the amount of water and gas in the outer vessel and the sample.

Experimental Setup
The saturated and satiated hydraulic conductivity of sand were measured with repetitive falling-head infiltration experiments. The respective setup is presented in Figure 2. The sample container consisted of an acrylic glass (PMMA) cylinder (7.2 cm inner diameter, 7.0 cm height) and a PMMA ring with a filter-mesh (size 99 µm) attached to its bottom. An o-ring made of nitrile rubber was used as a hydraulic seal, ensuring no water could escape between the ring and cylinder. The sample container was filled with the sand up to the height of approximately 5.0 cm. The container was placed into a larger, outer PMMA vessel (inner diameter 12.3 cm), where it was mounted onto five PMMA retainers. An overflow was inserted into the wall of the outer vessel (WADI-PZ, M12×1.5, 3.0-6.5 mm, Hugro Plastic Cable Gland, Waldkirch, Germany) through which a constant water level could be kept. To achieve a sufficiently large flow rate through the overflow, we attached a tube connected to a peristaltic pump (see Figure 2). The vessel was placed on a balance with 0.1 g precision to monitor the amount of water and gas in the outer vessel and the sample.

Sample Preparation
Two samples (ST1 and ST2) were packed with ST sand, while another two samples (FH1 and FH2) were prepared with FH sand. The samples were packed to a height of approximately 5.0 cm (the exact heights are given in Table 1) into cylindrical containers with an inner diameter of 7.2 cm. The sand packing was performed in degassed, deionized water, to minimize the introduction of entrapped air bubbles. The sand was added in 20 steps. After each addition step, the sand in the column was compacted by a small manual compactor, consisting of a rubber stopper and rod. The saturation by vacuum was used to ensure full saturation of the sand in the column [49]. The vessel with a submerged sample container was inserted into a vacuum chamber. The vacuum chamber was then evacuated for approximately 20 min at −95 kPa (relative to the atmospheric pressure). To minimize changes in pore geometry due to the resulting air bubble ebullition, a perforated plate loaded with 1.24 kg weight was placed on the sample's surface during this process.

Determination of Saturated Hydraulic Conductivity
To determine the value of saturated hydraulic conductivity, the first infiltration experiment with the repetitive falling-head was performed on the fully saturated sample, directly after it had been removed from the vacuum chamber. We assumed that the entrapped air content ω was zero at this point. The water used in infiltration experiments was kept in an open bottle one day before the experiment in order to equilibrate the temperature and dissolved air content to the laboratory conditions.
The workflow of the infiltration experiments is shown in Figure 3. At the first stage (Figure 3a), the sample container was submerged in degassed water. Then a simple repetitive falling-head infiltration experiment was carried out to determine the volumetric flux through the sample at known changes of hydraulic gradient. Here, a pin inserted in the sand (Figure 3b) was used to indicate a reference water level. A batch of water (40 cm 3 ) was manually added on top of the column each time the water level had dropped to the tip of the pin. The water was allowed to flow freely through the filter-mesh at the bottom of the sample. The saturated hydraulic conductivity, KS (cm•s −1 ), was determined in each interval by Darcy's law adapted for falling-head method [50]:

Sample Preparation
Two samples (ST1 and ST2) were packed with ST sand, while another two samples (FH1 and FH2) were prepared with FH sand. The samples were packed to a height of approximately 5.0 cm (the exact heights are given in Table 1) into cylindrical containers with an inner diameter of 7.2 cm. The sand packing was performed in degassed, deionized water, to minimize the introduction of entrapped air bubbles. The sand was added in 20 steps. After each addition step, the sand in the column was compacted by a small manual compactor, consisting of a rubber stopper and rod. The saturation by vacuum was used to ensure full saturation of the sand in the column [49]. The vessel with a submerged sample container was inserted into a vacuum chamber. The vacuum chamber was then evacuated for approximately 20 min at −95 kPa (relative to the atmospheric pressure). To minimize changes in pore geometry due to the resulting air bubble ebullition, a perforated plate loaded with 1.24 kg weight was placed on the sample's surface during this process.

Determination of Saturated Hydraulic Conductivity
To determine the value of saturated hydraulic conductivity, the first infiltration experiment with the repetitive falling-head was performed on the fully saturated sample, directly after it had been removed from the vacuum chamber. We assumed that the entrapped air content ω was zero at this point. The water used in infiltration experiments was kept in an open bottle one day before the experiment in order to equilibrate the temperature and dissolved air content to the laboratory conditions.
The workflow of the infiltration experiments is shown in Figure 3. At the first stage (Figure 3a), the sample container was submerged in degassed water. Then a simple repetitive falling-head infiltration experiment was carried out to determine the volumetric flux through the sample at known changes of hydraulic gradient. Here, a pin inserted in the sand (Figure 3b) was used to indicate a reference water level. A batch of water (40 cm 3 ) was manually added on top of the column each time the water level had dropped to the tip of the pin. The water was allowed to flow freely through the filter-mesh at the bottom of the sample. The saturated hydraulic conductivity, K S (cm·s −1 ), was determined in each interval by Darcy's law adapted for falling-head method [50]: where ∆H (cm) is a hydraulic head difference at the moment when the water level touches the pin, V is batch volume of added water and A is cross-section area of sample (V/A ≈ 1 cm), t is time of infiltration of the batch volume and L (cm) is the sample height. A range of ∆H was from 1.8 cm to 2.4 cm. After the infiltration experiment, the outer vessel with the sample was flooded precisely to the reference level indicated by the tip of the bent wire ( Figure 3c). It was then weighed to obtain the reference mass of the fully saturated sample. where ∆H (cm) is a hydraulic head difference at the moment when the water level touches the pin, V is batch volume of added water and A is cross-section area of sample (V/A ≈ 1 cm), t is time of infiltration of the batch volume and L (cm) is the sample height. A range of ∆H was from 1.8 cm to 2.4 cm. After the infiltration experiment, the outer vessel with the sample was flooded precisely to the reference level indicated by the tip of the bent wire ( Figure 3c). It was then weighed to obtain the reference mass of the fully saturated sample.

Determination of Satiated Hydraulic Conductivity and Entrapped Air
To determine the satiated hydraulic conductivity for varying contents of entrapped air, several cycles of ponded infiltration runs followed by drainage runs were performed. After each infiltration run, the sample was placed onto a standard laboratory sand box of the design similar to [49]. The box contained saturated fine sand with a high air entry value. The water within sand box was connected to a burette by a flexible hose. Tension was produced by the height difference between the sand surface in the sand box and the water level in the burette. The sample was gradually drained to a specific pressure head. Ten drainage-satiation runs were conducted, each providing one satiated hydraulic conductivity. The chosen pressure heads, as well as the corresponding and drainage durations, are given in Table 2. A perforated plate loaded with a weight of 0.62 kg was placed onto the surface of the sand sample during the drainage process to minimize possible pore geometry changes within the sample. The sample mass, Msamp,init, was determined after each drainage run to calculate the gas saturation of the drained sample, Sg,init. The calculation is described in Appendix A. Then, a ponding infiltration experiment was performed in the same way as described in Section 2.4, this time to obtain the satiated hydraulic conductivity, K(ω). A range of ∆H was from 2.0 cm to 2.5 cm.
The entrapped air content ω was calculated from the difference between the mass of the outer vessel with a satiated sample and the mass of the outer vessel with a degassed, fully saturated sample.
The obtained satiated hydraulic conductivity values were fitted with (i) the relationship published in Faybishenko [1]:

Determination of Satiated Hydraulic Conductivity and Entrapped Air
To determine the satiated hydraulic conductivity for varying contents of entrapped air, several cycles of ponded infiltration runs followed by drainage runs were performed. After each infiltration run, the sample was placed onto a standard laboratory sand box of the design similar to [49]. The box contained saturated fine sand with a high air entry value. The water within sand box was connected to a burette by a flexible hose. Tension was produced by the height difference between the sand surface in the sand box and the water level in the burette. The sample was gradually drained to a specific pressure head. Ten drainage-satiation runs were conducted, each providing one satiated hydraulic conductivity. The chosen pressure heads, as well as the corresponding and drainage durations, are given in Table 2. A perforated plate loaded with a weight of 0.62 kg was placed onto the surface of the sand sample during the drainage process to minimize possible pore geometry changes within the sample. The sample mass, M samp,init , was determined after each drainage run to calculate the gas saturation of the drained sample, S g,init . The calculation is described in Appendix A. Then, a ponding infiltration experiment was performed in the same way as described in Section 2.4, this time to obtain the satiated hydraulic conductivity, K(ω). A range of ∆H was from 2.0 cm to 2.5 cm.
The entrapped air content ω was calculated from the difference between the mass of the outer vessel with a satiated sample and the mass of the outer vessel with a degassed, fully saturated sample.
The obtained satiated hydraulic conductivity values were fitted with (i) the relationship published in Faybishenko [1]: and (ii) with the van Genuchten-Mualem model (VGM) for unsaturated hydraulic conductivity [18,21]: where K S (cm·s −1 ) is the saturated hydraulic conductivity, n (-), m(-) and l (-) are exponents, the S (−) is the water saturation (S = 1-S g,r ) and K 0 (cm·s −1 ) is the hydraulic conductivity obtained for a particular ω max (−). The parameter ω max is the maximum possible entrapped air content, which depends on the method chosen to saturate the sample. According to Faybishenko [1], it is not a soil property but rather a state variable. Faybishenko's equation is based on a power-law relationship from Averyanov [51], used for two-phase flow, where the residual saturation was replaced by satiated saturation. The equations were verified on high columns (0.5-5 m in height) of soil during long-term infiltration.
Mualem [18] derived a general equation of relative hydraulic conductivity by the model of pore size distribution. Mualem then used Brooks and Corey's equation of retention curve [20] and the measurements of 45 soils from literature to estimate exponent l, which can be positive or negative. He obtained the value of l = 0.5, which is commonly used. van Genuchten [21] combined Mualem's model of pore distributions with an equation for the retention curve, and verified the equation on five samples with good agreement. The equation is known as the Mualem-van Genuchten equation.
We moreover investigated how residual gas saturations during each infiltration run were related to the initial gas saturation. Here we evaluated relationships proposed by Land [14] and by Aissaoui [52]. Land [14] uses the simple term to describe the relationship between initial gas saturation S g,init (−) and residual gas saturation S g,r (−): 1 where S g,r,max (−) is a maximum value of the S g,r obtained for S g,init = 1. Since the S g,init value was in most cases less than 1, we propose a modified version of Land's term: where S g,r,max is maximum value of the S g,r and the S g,init,max is appropriate value of the S g,init . The term proposed by Aissaoui reads: S g,r = S g,r,max S g,C S g,init when S g,init < S g,C (6a) S g,r = S g,r,max when S g,init ≥ S g,C where S g,C (−) is the critical gas saturation corresponding to the point where the maximum trapped saturation, S g,r,max, is reached. Land derived the relationship on the published experimental data on six sandstone samples [14]. He presented figures of six plots with good agreement between fitted functions and data. Aissaoui [52] proposed the linear relationship between S g,r and S g,init where after the reaching of the critical gas saturation S g,C the S g,r stays constant. According Pentland [53], Aissaoui's equation gives a better prediction of entrapped air saturation for unconsolidated soil than the Land's equation.

X-Ray Computed Tomography of Samples
In addition to the ten drainage and infiltration cycles specified in Table 2, two drainage and infiltration experiments were conducted in combination with X-ray computed micro-tomography (CT) imaging. In this fashion, we obtained information on the spatial distribution of entrapped air bubbles. The imaging was conducted using a GE Phoenix v|tome|x m industrial X-ray scanner with a 240 kV X-ray tube and a tungsten target, located at the Swedish University of Agricultural Sciences in Uppsala (Sweden). The result is a stack of 16-bit horizontal slices with pixel size 47 µm and 47 µm slice thickness. The samples were first fully saturated in a vacuum chamber, after which they were X-rayed while submerged inside the water-filled outer vessel. The samples were then drained at a tension of −30 hPa in the case of material ST and −35 hPa in the case of material FH. The drained ST samples were scanned in air-filled outer vessels. Then an infiltration experiment was applied to all samples followed by immediate imaging. The samples were still submerged inside the water-filled outer vessel during imaging. The second drainage cycle was carried out under a pressure of −50 hPa before a second infiltration run and more X-ray scans were conducted. Details of the experiment schedule are summarized in Table 3 and Figure 4. a 240 kV X-ray tube and a tungsten target, located at the Swedish University of Agricultural Sciences in Uppsala (Sweden). The result is a stack of 16-bit horizontal slices with pixel size 47 µm and 47 µm slice thickness. The samples were first fully saturated in a vacuum chamber, after which they were X-rayed while submerged inside the water-filled outer vessel. The samples were then drained at a tension of −30 hPa in the case of material ST and −35 hPa in the case of material FH. The drained ST samples were scanned in air-filled outer vessels. Then an infiltration experiment was applied to all samples followed by immediate imaging. The samples were still submerged inside the water-filled outer vessel during imaging. The second drainage cycle was carried out under a pressure of −50 hPa before a second infiltration run and more X-ray scans were conducted. Details of the experiment schedule are summarized in Table 3 and Figure 4.

Image Preprocessing
The CT images were preprocessed and analyzed using the software Fiji (ImageJ) [54]. Preprocessing includes the elimination of the noise, edge enhancement and artifact removal [55].
First, all horizontal image layers in the region of interest were selected. The location of the topmost layer depended on the shape of the surface because after the infiltration experiment the surface was not exactly horizontal. The bottommost selected layer was defined by the location of the sample retainer. Then a 3D-median filter with radius 2 voxels and an unsharp mask was used.
The tomograms were further adjusted to minimize the effect of beam hardening. Beam hardening refers to the non-linear filtering of soft (or low-frequency) X-rays during the beam's passage through a sample, leaving only harder (or-high-frequency) X-rays in the beam [43]. For soil samples with a cylindrical shape, beam hardening manifests itself as a bright ring artifact along the column's perimeter [44].
More specifically, the beam-hardening artifacts in this study were detected in two forms: First, along with the entire sample, a general darkening of the image towards the horizontal center of the sand sample. Second, below the PMMA ring, the gray values were attenuated due to the extra X-rayfiltering effect of the ring. Schedule of the experiment with X-ray micro-tomography imaging. CT imaging (shown in red color) was conducted after the first infiltration (I1) and after the second infiltration (I2). CT imaging was also performed after drainage (D1 and D2) experiment and character D indicates CT scanning after drainage. The duration of scan was 1 hour and duration of infiltration experiment was from 10 to 20 minutes.

Image Preprocessing
The CT images were preprocessed and analyzed using the software Fiji (ImageJ) [54]. Preprocessing includes the elimination of the noise, edge enhancement and artifact removal [55].
First, all horizontal image layers in the region of interest were selected. The location of the topmost layer depended on the shape of the surface because after the infiltration experiment the surface was not exactly horizontal. The bottommost selected layer was defined by the location of the sample retainer. Then a 3D-median filter with radius 2 voxels and an unsharp mask was used.
The tomograms were further adjusted to minimize the effect of beam hardening. Beam hardening refers to the non-linear filtering of soft (or low-frequency) X-rays during the beam's passage through a sample, leaving only harder (or-high-frequency) X-rays in the beam [43]. For soil samples with a cylindrical shape, beam hardening manifests itself as a bright ring artifact along the column's perimeter [44].
More specifically, the beam-hardening artifacts in this study were detected in two forms: First, along with the entire sample, a general darkening of the image towards the horizontal center of the sand sample. Second, below the PMMA ring, the gray values were attenuated due to the extra X-ray-filtering effect of the ring.
To remove the artifact caused by the PMMA ring, the intensity of each slide was multiplied by a correction coefficient α i . The average image intensity of a 60 pixels wide region of interest located within the wall of the sample cylinder (Figure 5a) was measured for each slice (I i ). The value of correction coefficient α i for the slice i was then: where I avg is the average of all the values of I i and β is a correction factor. The effect of the parameter β on profiles of average image intensity in the central part of the sample (shown in red in Figure 5c) is depicted in Figure 5d-f. A factor of β = 0.5 led to the best correction results. To remove the artifact caused by the PMMA ring, the intensity of each slide was multiplied by a correction coefficient αi. The average image intensity of a 60 pixels wide region of interest located within the wall of the sample cylinder (Figure 5a) was measured for each slice (Ii). The value of correction coefficient αi for the slice i was then: where Iavg is the average of all the values of Ii and β is a correction factor. The effect of the parameter β on profiles of average image intensity in the central part of the sample (shown in red in Figure 5c) is depicted in Figure 5d-f. A factor of β = 0.5 led to the best correction results. The next adjustment was related to the second artifact, i.e., the darkening towards the horizontal center of the images. This was done by creating vertical cross-sections through the 3D images. Then, a 9.4 mm (200 pixels) wide strip along the diameter was selected and a profile of intensity was evaluated (Figure 6a,b). This profile was fitted with a center-axis parabola. The ratio between the value on the parabola and the edge values were given by factor ρ. Then each of the images was modified by dividing the intensity by ρ, given the distance from the center of the sample. The next adjustment was related to the second artifact, i.e., the darkening towards the horizontal center of the images. This was done by creating vertical cross-sections through the 3D images. Then, a 9.4 mm (200 pixels) wide strip along the diameter was selected and a profile of intensity was evaluated (Figure 6a,b). This profile was fitted with a center-axis parabola. The ratio between the value on the parabola and the edge values were given by factor ρ. Then each of the images was modified by dividing the intensity by ρ, given the distance from the center of the sample.

Image Analysis
The images were binarized into two classes-one class depicting the gas phase and the other all other phases, namely the water and solid phases. An automatic binarization was not successful in this case due to poorly recognizable peaks on the intensity histogram. The threshold values were obtained by visual inspection for each sample. Then a morphological opening with a kernel radius of 1 voxel was carried out using the ImageJ plugins 3D Erode and 3D Dilate. We then removed 'holes' in the segmented air bubbles, which represented segmentation artifacts, using the 3D Fill Holes plugin [56].
We determined the imaged gas volume content as a function of the height of the sand samples as the ratio between the area of gas and the column cross-section for each image slice. These results were also used to determine the air volume content of the whole sample. The plugin BoneJ [57] was used to evaluate the thickness of the air bubbles, defined as the diameter of the largest sphere that fits within individual bubbles.

Initial Gas Saturation
The relationship between initial and residual gas saturation as fitted using Aissaoui's and Land's equations is shown in Figure 7 together with the results of Pentland et al. We observed a positive correlation between Sg,init and the residual saturation Sg,r. The maximum of the initial gas saturation ranged between 80% to 90%. The evaluated parameters for Land's and Aissaoui's models are shown in Table 4. In our case the growth of Sg,r value did not converge to an upper limit as it was supposed by Aissaoui. On the basis of experiments with octane used as a non-wetting phase, Pentland et al.

Image Analysis
The images were binarized into two classes-one class depicting the gas phase and the other all other phases, namely the water and solid phases. An automatic binarization was not successful in this case due to poorly recognizable peaks on the intensity histogram. The threshold values were obtained by visual inspection for each sample. Then a morphological opening with a kernel radius of 1 voxel was carried out using the ImageJ plugins 3D Erode and 3D Dilate. We then removed 'holes' in the segmented air bubbles, which represented segmentation artifacts, using the 3D Fill Holes plugin [56].
We determined the imaged gas volume content as a function of the height of the sand samples as the ratio between the area of gas and the column cross-section for each image slice. These results were also used to determine the air volume content of the whole sample. The plugin BoneJ [57] was used to evaluate the thickness of the air bubbles, defined as the diameter of the largest sphere that fits within individual bubbles.

Initial Gas Saturation
The relationship between initial and residual gas saturation as fitted using Aissaoui's and Land's equations is shown in Figure 7 together with the results of Pentland et al. We observed a positive correlation between S g,init and the residual saturation S g,r . The maximum of the initial gas saturation ranged between 80% to 90%. The evaluated parameters for Land's and Aissaoui's models are shown in Table 4. In our case the growth of S g,r value did not converge to an upper limit as it was supposed by Aissaoui. On the basis of experiments with octane used as a non-wetting phase, Pentland et al. [53] assumed that the upper limit of S g,r corresponds to the occupancy of all large pores by the non-wetting phase. They also assumed that small pores continue to fill causing further increase of S g,init , while S g,r does not increase, because gas is displaced from the small pores. In our study the upper limit of S g,r was not detected, probably due to the apparently unstable pore geometry in the upper portion of the sample, where pores enlarged with increasing S g,init . This effect will be discussed in detail in the next paragraph. Yet, Aissaoui's equation results in a better fit of the experimental values than the Land's equation, which is in agreement with Pentland's study.
Water 2019, 11, x FOR PEER REVIEW 11 of 20 [53] assumed that the upper limit of Sg,r corresponds to the occupancy of all large pores by the nonwetting phase. They also assumed that small pores continue to fill causing further increase of Sg,init, while Sg,r does not increase, because gas is displaced from the small pores. In our study the upper limit of Sg,r was not detected, probably due to the apparently unstable pore geometry in the upper portion of the sample, where pores enlarged with increasing Sg,init. This effect will be discussed in detail in the next paragraph. Yet, Aissaoui's equation results in a better fit of the experimental values than the Land's equation, which is in agreement with Pentland's study.   Table 5 lists the volumetric air contents ω of the whole samples, obtained gravimetrically and by X-ray image analysis. The amount of entrapped air, after infiltrations I1 and I2, was higher for FH sand. The values of ω obtained gravimetrically were higher than the values obtained by X-ray imaging. This was probably caused by the fact that bubbles under certain sizes were not detected being below the image resolution. The field of view which determines the image resolution is given by the size of the sample. The size of the sample was designed to be larger than the expected representative elementary volume of sand, which led to an image resolution that neglected the smaller bubbles. Table 5 shows that the gravimetrically measured ω values were significantly larger after drainage than the ω values measured after each infiltration experiment. However, this effect is not reproduced in the X-ray imaging results. On the contrary, we observed the opposite. This could be due to the presence of predominantly small bubbles after drainage that remained too small to be detected in the image data, as can be seen in Figure 8, where vertical cross-sections of the binarized images of the sample ST2 are shown. After drainage (D1) there was a large number of small bubbles, mostly located in the upper half of the sample. We can assume that in this case an additional large   Table 5 lists the volumetric air contents ω of the whole samples, obtained gravimetrically and by X-ray image analysis. The amount of entrapped air, after infiltrations I1 and I2, was higher for FH sand. The values of ω obtained gravimetrically were higher than the values obtained by X-ray imaging. This was probably caused by the fact that bubbles under certain sizes were not detected being below the image resolution. The field of view which determines the image resolution is given by the size of the sample. The size of the sample was designed to be larger than the expected representative elementary volume of sand, which led to an image resolution that neglected the smaller bubbles. Table 5 shows that the gravimetrically measured ω values were significantly larger after drainage than the ω values measured after each infiltration experiment. However, this effect is not reproduced in the X-ray imaging results. On the contrary, we observed the opposite. This could be due to the presence of predominantly small bubbles after drainage that remained too small to be detected in the image data, as can be seen in Figure 8, where vertical cross-sections of the binarized images of the sample ST2 are shown. After drainage (D1) there was a large number of small bubbles, mostly located in the upper half of the sample. We can assume that in this case an additional large number of small bubbles existed, which remained undetected because they were below the resolution of the X-ray images. After the first infiltration experiment (I1), larger pores were created by displacing nearby grains of sand (Figure 8b) and the air had aggregated into larger bubbles, probably due to the effect of Ostwald ripening [58]. After the second infiltration experiment (I2), the bubbles became even larger (Figure 8c). It seems that the higher amount of the visible air bubbles in the sample occurred due to the enlargement of pores. number of small bubbles existed, which remained undetected because they were below the resolution of the X-ray images. After the first infiltration experiment (I1), larger pores were created by displacing nearby grains of sand ( Figure 8b) and the air had aggregated into larger bubbles, probably due to the effect of Ostwald ripening [58]. After the second infiltration experiment (I2), the bubbles became even larger (Figure 8c). It seems that the higher amount of the visible air bubbles in the sample occurred due to the enlargement of pores.  The air saturation profile along the height of the sample is shown in Figure 9. It should be noted that only relatively large bubbles were detected by x-ray CT. It is evident that the air is predominantly distributed in the upper part of the sample, from 25 mm upwards, especially after infiltration experiments (I1, I2). The accumulation of air in large bubbles at the upper half of the sample was probably caused by the buoyancy force in combination with Ostwald ripening.

Air Content and Air Distribution
The amount of air in the upper part is approximately doubled between I1 and I2. The air in the lower 5 mm is predominantly present in newly created cracks. The air saturation profile along the height of the sample is shown in Figure 9. It should be noted that only relatively large bubbles were detected by x-ray CT. It is evident that the air is predominantly distributed in the upper part of the sample, from 25 mm upwards, especially after infiltration experiments (I1, I2). The accumulation of air in large bubbles at the upper half of the sample was probably caused by the buoyancy force in combination with Ostwald ripening.
The amount of air in the upper part is approximately doubled between I1 and I2. The air in the lower 5 mm is predominantly present in newly created cracks.
The frequency of the thickness of the bubbles can be seen in Figure 10. To eliminate the effect of the air-filled fractures, only parts of the image 2 cm above the bottom were analyzed. In all the cases, the air bubble thickness increased between stages I1 and I2. More specifically, the thickness of the smallest 50% of the bubbles changed from 300 to 405 µm for ST1, from 370 to 470 µm for ST2, from 310 to 405 µm for FH1, and from 370 to 470 µm for FH2. The maximum thickness of all detected air bubbles was approximately 2000 µm. The frequency of the thickness of the bubbles can be seen in Figure 10. To eliminate the effect of the air-filled fractures, only parts of the image 2 cm above the bottom were analyzed. In all the cases, the air bubble thickness increased between stages I1 and I2. More specifically, the thickness of the smallest 50% of the bubbles changed from 300 to 405 µm for ST1, from 370 to 470 µm for ST2, from 310 to 405 µm for FH1, and from 370 to 470 µm for FH2. The maximum thickness of all detected air bubbles was approximately 2000 µm.  Figure 11 shows the relationship between the satiated hydraulic conductivity and the air content in the sand. Also depicted are fits with Faybishenko's (2) and VGM functions (3). The respective model parameters n and m were obtained by the least squares method. The corresponding values are shown in Table 6. It can be seen that the results of the infiltration run on the ST samples conducted next to the X-ray scanner followed the same trend as observed for the other 10 experiments. This was not the case for the FH samples, especially not for FH2. One possible explanation could be the smaller  Figure 11 shows the relationship between the satiated hydraulic conductivity and the air content in the sand. Also depicted are fits with Faybishenko's (2) and VGM functions (3). The respective model parameters n and m were obtained by the least squares method. The corresponding values are shown in Table 6. It can be seen that the results of the infiltration run on the ST samples conducted next to the X-ray scanner followed the same trend as observed for the other 10 experiments. This was not the case for the FH samples, especially not for FH2. One possible explanation could be the smaller number of measured points and the experiment FH2 can be considered as an outlier. The coefficients of determination for both models are given in Table 6. The coefficients of determination for Faybishenko's model for batch with CT scanning are high due to fitting three points only. K S of the coarser sand ST is almost twice as high as the K S of the finer sand FH. The mean K S of FH sand (0.035 cm·s −1 ) was similar to the value reported for the same sand    Faybishenko's model led to similar fits of the K(ω) as the VGM equation. An exception was sample FH1, where the overall shape of the K(ω) was abnormal. In comparison, Marinas et al. [24]

Conclusions
For the purpose of simulation of two-phase flow, the relationship between entrapped air content and satiated hydraulic conductivity K(ω) of two coarse sands was studied. Simultaneously with falling head infiltration and drainage experiments, the X-ray computed tomography imaging was performed to reveal temporal changes of the entrapped air distribution and the bubble shapes within the sample.
The positive correlation between initial air saturation S g,init and residual air saturation S g,r was found; this result is in agreement with Land [14] or Mansoori et al. [59]. The maximum value of S g,r was 0.31 for ST sand and 0.38 for FH sand. A good fit of the S g,r (S g,init ) relationship was obtained using Aissaoui's equation.
As expected, the saturated hydraulic conductivity of coarser sand ST was consistently higher than the hydraulic conductivity of finer sand FH. The hydraulic conductivity of samples at full saturation was approximately twice as high as the hydraulic conductivity of samples with maximal entrapped air content. The relationship between air volume content and satiated hydraulic conductivity was successfully fitted with Faybishenko's equation. Faybishenko's model seems to be appropriate to express the K(ω) relationship.
X-ray computed tomography revealed that the size of air bubbles and clusters increased significantly with each infiltration-drainage cycle. The final shapes of observed bubbles appeared to be elongated in the horizontal direction. The increase in bubble size was not spatially uniform within an individual sample. Instead, air bubbles grew larger close to the surface of the sample. Bubble ebullition probably also led to the formation of the larger voids in the sand packing.
The obtained K(ω) relationship can be used as a characteristic in newly developed numerical models of two-phase flow, such as the model by Fucik and Mikyska [60]. It is evident that the entrapped air significantly influences measured values of hydraulic conductivity even in samples of simple material and geometry. The presented methodology can be used to determine the K(ω) in other natural and artificial porous materials.