How Wetting and Drainage Cycles and Wetting Angle Affect Capillary Air Trapping and Hydraulic Conductivity: A Pore Network Modeling of Experiments on Sand

Abstract

Entrapped air in porous media can significantly affect water flow but simulations of air entrapment are still challenging. We developed a pore-network model using quasi-static algorithms to simulate air entrapment during spontaneous wetting and subsequent drainage processes. The model, implemented in OpenPNM, was tailored to replicate an experiment conducted on a medium-sized unconsolidated sand sample. We started building the model with three types of relatively small networks formed by 54,000 pore bodies which we used to calibrate basic network topological parameters by fitting the model to the water retention curve and the saturated hydraulic conductivity of the sand sample. Using these parameters, along with X-ray image data (µCT), a larger network formed by over 250,000 pore bodies was introduced in the form of stacked sub-networks where topological parameters were scaled along the z-axis. We investigated the impact of two different contact angles on air entrapment. For a contact angle of 0, the model showed good agreement with the experimental data, accurately predicting the amount of entrapped air and the saturated hydraulic conductivity. On the contrary, for a contact angle of π/4, the model provided reasonable accuracy for saturated hydraulic conductivity but overestimated the amount of entrapped air. Overall, this approach demonstrated that a reasonable match between simulated and experimental data can be achieved with minimal computational costs.

1. Introduction

It has been experimentally demonstrated that entrapped air in a non-wetting phase strongly affects water flow through porous media [1,2]. In particular, the presence of entrapped air in preferential flow pathways of soil can create obstructions that significantly reduce the overall flow of water [3,4,5].

Recent advances in imaging techniques, such as X-ray microtomography, have provided deeper insights into the morphology of residual non-wetting phases [6]. For instance, Scanziani et al. [7] utilized this technique to visualize the dynamics of gas trapping in CO₂ storage scenarios, emphasizing the role of the pore structure in influencing the trapping efficiency. Wang et al. [8] demonstrated that both water-wet and CO₂-wet conditions enhance the trapping of CO₂ compared to weakly water-wet or neutral conditions.

Continuum models mainly implement entrapment and fate of air using the two-phase (or multiphase) flow approach [9,10]. It has frequently been assumed that the residual air content is constant, but this assumption is often not valid [11,12]. Alternatively, pore-scale models take the three-dimensional geometry of the pore network into account. Various pore-scale modeling approaches have been used in the past, such as the volume-of-fluid method [13,14], the lattice Boltzmann method [13,15], and pore-network modeling [16,17]. The volume-of-fluid and lattice Boltzmann methods are computationally demanding, and they have typically been used to describe only small domains of the order of several pore bodies in each direction [18].

Pore network models (PNMs) have been shown to successfully and efficiently predict numerous aspects of multiphase flow and transport processes in porous media [19,20,21,22,23]. Algorithms used in PNMs are often based on percolation theory [22,23], which simulate processes based on physical forces [24,25,26,27] that require little computational power and time [21,28,29]. A thorough overview of PNMs has been given in seminal works by Blunt et al. [19] and Sahimi [30].

The geometry and topology of the pore network are key parameters of PNMs [31]. The topology of the network determines the location of pore bodies and their connection through pore throats. Early studies used pore throats and pore bodies with circular cross sections [32,33]. Such approaches are being increasingly replaced by pore throats and pore bodies with angular cross sections (triangle, square, hyperbolic polygons, etc.), which allow simulations of corner flow [34]. The geometry and topology of the pore network can be constructed as random computer-generated structures [35,36,37] or can be derived directly from highly resolved image data of porous media [38,39]. The pores network generated based on high resolution 3D images, mostly obtained by X-ray imaging on synchrotrons provide the most realistic description [40,41,42], albeit only for sample volumes ranging from several cubic millimeters [41] to several cubic centimeters [43]. The small size of the domain is a disadvantage when investigating the hydraulic properties of complex porous media such as soils, which often require larger sample sizes to meet the prerequisites of a representative elementary volume.

The purpose of this study was to develop and parameterize a pore network model using the OpenPNM platform [44] to simulate air entrapment and water flow through coarse sand during repeated wetting and drainage cycles. Specifically, we sought to determine if a pore network constructed on (i) the water retention curve and (ii) the saturated hydraulic conductivity of the porous material, as well as (iii) the average diameter of entrapped air bubbles with depth obtained from X-ray image data could provide a representation of the experimental data on wetting and draining.

2. Materials and Methods

2.1. Model Background

2.1.1. Network Generation

The model presented is based on the OpenPNM framework [44]. This open-source package has previously been used to characterize diffusivity in the thin, porous layer [45] or to solve multiphase transport through the fuel cell [46,47].

The generated networks comprised interconnected pore bodies and pore throats, both a square cross-section. The network topology is defined by the coordinates of the pore bodies and a list of pairs of pore bodies connected by the corresponding pore throats. One of the parameters that describe the network is the coordination number z (-). The value z indicates the average number of pore throats connected to a single pore body. Another parameter, the aspect ratio, indicates the ratio between the pore body and the radius of the pore throat.

2.1.2. Physics Considered in the Simulation Model

We simulated both the drainage and wetting of water in the sand sample as a flow of an invading fluid phase towards a defending phase that is initially present in the network. This process is described by invasion percolation and simulated by considering two key parameters: the entry pressure and the conductivity of the individual network elements, namely the pore bodies and throats.

The entry pressure of the element corresponds to the capillary pressure (difference between the non-wetting and wetting phase pressures) at which the invading phase begins to enter the element filled by the defending phase. During drainage, when the capillary pressure drops, it eventually falls below the entry pressure of the element, then the element is filled by the non-wetting (invading) phase by a piston-like displacement. During the imbibition process, the capillary pressure increases, which leads to the filling of the element with the wetting phase when the entry pressure of the element is exceeded.

Our simulation model was based on a simplified representation of the porous medium using percolation theory. We assumed that only pore throats were relevant for the drainage process. The simplification is applicable because the entry pressure of the pore throat is generally lower than the entry pressure of the pore body because the size of the throat is smaller than the size of the pore body. Therefore, when a pore throat is drained, it means that the adjacent pore body also drains. However, capillary pressure increases during imbibition. In this way, elements with a smaller radius and lower entry pressure are saturated preferably, from which it follows that the pore bodies are more important for the imbibition process.

The entry pressure of a pore throat with a regular polygon cross section is described by Equation (1) [48]:

$$P_{e,t}={\displaystyle \frac{\sigma }{r_t}}\left\{\cos \theta +\sqrt{{\displaystyle \frac{\tan {\alpha }_c}{2}}\left(\sin 2\theta +\pi -2\left({\alpha }_c+\theta \right)\right)}\right\}$$

(1)

where r_t (m) is the inscribed radius of the pore throat, θ (rad) is the contact angle, and ${\alpha }_c$ (rad) is the half-angle of the corners (π/4 for square).

For the imbibition algorithm, the entry pressure of the pore body, P_e,p (Pa), depends on the geometrical characteristics which determine the radius of curvature of the interface R_n (m):

$$P_{e,p}={\displaystyle \frac{2\sigma }{R_n}}$$

(2)

where σ (N∙m⁻¹) is the interfacial tension. The radius of curvature was determined following [49], based on the calculation of the entry pressure for a spherical pore body. The relationship is identical for a cubic pore body since the interface geometry remains the same (see Figure 1a). The pore bodies filled with the non-wetting phase with a maximum angular pitch of their axes (α_p in Figure 1a) were used to calculate R_n. If there was more than one combination of throats with the same maximum angular pitch, then the combination with larger throats was used. The radius of curvature of the interface (R_n) for throats with equal radii was determined by Equation (3):

$$a{R}_n^2+b{R}_n+c=0$$

(3)

where:

$$\begin{gathered} a={\left(\mathit{cos}{\beta }_S+{\displaystyle \frac{\mathit{sin}{\beta }_S}{\mathit{tan}\frac{\alpha }{2}}}\right)}^2-1 \\ b=-2\left(\mathit{cos}{\beta }_S+{\displaystyle \frac{\mathit{sin}{\beta }_S}{\mathit{tan}\frac{\alpha }{2}}}\right)\left({\displaystyle \frac{r_t}{\mathit{sin}\frac{\alpha }{2}}}+R_p\mathit{cos}{\displaystyle \frac{\gamma }{2}}\right) \\ c={\left({\displaystyle \frac{r_t}{\mathit{sin}\frac{\alpha }{2}}}+R_p\mathit{cos}{\displaystyle \frac{\gamma }{2}}\right)}^2+{\left(R_p\mathit{sin}{\displaystyle \frac{\gamma }{2}}\right)}^2 \end{gathered}$$

(4)

where R_p (m) is pore body radius, γ (rad) is selected as α_p as discussed in [49], and ${\beta }_S$ = π/2 − θ − α_p (rad). When the pore throats were of different radii, the equivalent pore radius was determined according to [49] using the smaller of the two pore throat radii, and the value of angular space remained α_p. The entry pressure was then determined by Equations (2)–(4) with the equivalent parameters.

In addition, the capillary pressure within the pore network changes hydrostatically with elevation in the gravitational field. Hydrostatic pressure was accounted for by changing the entry pressure according to height.

$$P_{e\_fin,el}=P_{e,el}-{\rho }_{w}g{z}_{el}$$

(5)

where P_{e_fin,el} (Pa) is the final entry pressure of element el, ρ_w (kg∙m⁻³) is the density of the water, and z_el (m) is the height of the center of the element el from the bottom of the pores network.

Furthermore, the capillary pressure applied during drainage, $P_c^{*}$ (Pa), is suitable for calculating the radius of curvature of the water film in the corners of the pore throats and bodies (Figure 1b):

$$r_{nw}={\displaystyle \frac{\sigma }{P_c^{*}}}$$

(6)

The area of water in the cross-section, A_w (m²), is determined by

$$A_w=N_c{r}_{nw}^2\left[{\displaystyle \frac{\cos \theta }{\sin {\alpha }_c}}\cos \left({\alpha }_c+\theta \right)-{\displaystyle \frac{\pi }{2}}+{\alpha }_c+\theta \right]$$

(7)

where N_c (-) is the number of corners, and the air volume in the whole element, V_g (m³), is:

$$V_g=V_e-V_w=V_e-{\displaystyle \frac{A_w}{A_e}}{V}_e$$

(8)

where A_e (m²) is the area of the cross-section of the element, V_e (m³) is the volume of the element, and V_w (m³) is the volume of the water film. For the water film to form in the corners, the condition of θ <π/2 − α_c must be met.

We assumed the Poiseuille law in each pore throat to calculate its hydraulic conductivity [50]. The flow rate of the wetting phase from the pore body i to an adjacent pore body j, q_ij (m³∙s⁻¹), was determined as follows:

$$q_{ij}=g_{ij}∆P_{ij}$$

(9)

where ∆P_ij (Pa) is the difference between the pressure of the wetting phase in the pore bodies i and j, and g_ij (m³∙Pa⁻¹∙s⁻¹) is the throat hydraulic conductivity of the pore throat that connects i and j. The conductivity of the throat for single-phase flow, when the throat was filled by the wetting phase, was:

$$g_{ij}^{w,1-p}={\displaystyle \frac{\pi r_{eff}^4}{8{\mu }_{w}l}}$$

(10)

where r_eff (m) is the hydraulic radius, µ_w (Pa∙s) is the dynamic viscosity of the wetting phase, l (m) is the length of the throat, and the superscript w, 1 − p indicates the single-phase water flow. The hydraulic radius was determined following Bryant and Blunt [51]:

$$r_{eff}={\displaystyle \frac{\sqrt{A_t/\pi }+r_t}{2}}$$

(11)

where r_t (m) is the inscribed radius of the pore throat, and A_t is the area of the cross-section of the pore throat. The wetting phase conductivity of an unsaturated pore throat (water flow in corner water films) is given by [52]:

$$g_{ij}^{w,2-p}={\displaystyle \frac{A_w{r}_{nw}^2}{\beta {\mu }_{w}l}}$$

(12)

where the superscript w, 2 − p indicates two-phase water flow, and the resistance factor β (-) was determined by Equation (13) according [53]:

$$\beta ={\displaystyle \frac{12 {\sin }^2{\alpha }_c{\left(1-B\right)}^2{\left({\psi }_3+fB{\psi }_2\right)}^2}{{\left(1-\sin {\alpha }_c\right)}^2{B}^2{\left({\psi }_1-B{\psi }_2\right)}^2}}$$

(13)

where its coefficients are

$$\begin{gathered} {\psi }_1={\cos }^2\left({\alpha }_c+\theta \right)+\cos \left({\alpha }_c+\theta \right)\sin \left({\alpha }_c+\theta \right)\tan {\alpha }_c, { \psi }_2=1-{\displaystyle \frac{\theta }{\frac{\pi }{2}-{\alpha }_c}} \\ {{ \psi }_2=1-{\displaystyle \frac{\theta }{\frac{\pi }{2}-{\alpha }_c}}, \psi }_3={\displaystyle \frac{\cos \left({\alpha }_c+\theta \right)}{\cos {\alpha }_c}} \\ B=\left({\displaystyle \frac{\pi }{2}}-{\alpha }_c\right)\tan {\alpha }_c \end{gathered}$$

(14)

here the f = 0.

Appendix A describes both existing and new algorithms used in pore network modeling.

2.2. Experiments Carried Out on Sand Columns

We conducted two different experiments on columns packed with the same coarse sand. The first was the measurement of the primary drainage branch of the water retention curve on the sand column performed using the sandbox apparatus [54].

The second experiment [55] consisted of a series of drainage–imbibition cycles to study the relationship between the entrapped volumetric air content and the saturated and satiated hydraulic conductivity (Figure 2a).

The material used in this study was coarse sand (ST 03/08, Sklopisek Strelec, Czech Republic). Particle sizes ranged from 0.1 to 1.0 mm. The value of the uniformity coefficient (D₆₀/D₁₀,) was 2.2 and the effective particle size (D₁₀) was 0.29 mm. The particle density was 2.62 g∙cm³. The dry bulk density of the sample was calculated from the mass and volume of the dry sand sample (1.67 g∙cm³). Sample porosity was determined from particle density and the dry bulk density (0.354 cm³∙cm⁻³).

The water retention curve was used to fit the van Genuchten model [56], where we found the parameters α = 0.048, n = 3.976, m = 0.749.

The second experiment started with a sand sample that had been fully saturated in a vacuum chamber. The value of ω_r was assumed to be equal to 0 due to full saturation. The saturated hydraulic conductivity (K_s) was determined by repeating the fall-head method. (Figure 2b) with manual additions of 40 cm³ dose of water each time the water level dropped to the pin. The constant water level was kept around the sample by continuously removing excess water using a pump. The measurement was repeated ten times.

The mass of fully saturated sand was recorded and was later used as a reference for gravimetric measurements of the entrapped air content.

The drainage and spontaneous imbibition cycles were performed (Figure 2a). Drainage was carried out on a sand table with the given pressure head. The value of ω_init was obtained gravimetrically, by weighing the container with the sand column after drainage. Imbibition after each drainage step was performed by flooding the sample upwards through the bottom. By filling the outer container to the standard level, it was possible for each step to determine the mass of water replaced by the trapped air relative to the initial fully saturated state and to calculate ω_r. The satiated hydraulic conductivity (hydraulic conductivity at the actual residual gas saturation) K(ω) was measured after each imbibition using the same falling head method as used for the measurements of K_s. A total of ten drainage and imbibition cycles were completed. Furthermore, three-dimensional X-ray µCT images of the sample in the drained and flooded (saturated or satiated) states were recorded for two additional cycles. The experiments were detailed in [55].

2.3. Description of Numerical Simulations

Two contact angle values were considered for the simulations. The first contact angle θ = π/4 represents weakly wet conditions. This value was selected to represent a realistic value for the coarse sand. It is still small enough to prevent snap-off effects. Snap-off cannot occur for θ = π/4, due to the condition of θ <π/2 − α_c required for the creation of water film, because the water will not adhere in the corners [52]. The second value of θ was chosen to be 0, which well represents the state of very wet sand. Both θ-values were used for generating networks and simulations.

2.3.1. Small Networks Generation

A uniform network of 30 × 30 × 60 pore bodies was generated which has been shown to be a sufficient size for this purpose [31]. Because hydraulic characteristics mainly depend on the topology of the pore network and size distribution of the pore network elements, three different approaches for pore network generation were used.

The first network type was a cubic network with a connectivity equal to 6 and a spacing of 180 μm, with the truncated Weibull distribution of pore body radii:

$$R_p=\left(R_{max}-R_{\mathrm{m}\mathrm{i}\mathrm{n}}\right){\left(-\sigma \ln \left(x\left(1-e^{-{\displaystyle \frac{1}{\sigma }}}\right)+e^{-{\displaystyle \frac{1}{\sigma }}}\right)\right)}^{{\displaystyle \frac{1}{\delta }}}+R_{min}$$

(15)

where x (-) was a random number between 0 and 1. Parameters of this network were the minimal radius, R_min (m), the maximal radius, R_max (m), and the parameters σ (-) and δ (-).

The second network was a cubic network constructed with a connectivity equal to 26 and a spacing of 180 μm. The connectivity was then randomly reduced for the entire network to an integer value of the coordination number z between 4 and 8. Again, the Weibull distribution of the pore body radii was used (15). The parameters of this network type were R_min, R_max, σ, δ and z.

The third network type had a random spatial distribution of the pore body centers; the number of pores bodies was N_p (). The throats that connected the pore bodies were generated by the OpenPNM Gabriel algorithm. The maximum radius of the pore bodies was determined. Then the radius of each pore body was determined from the maximum radius by multiplying by the coefficient ε (-), which was constant for the entire network. The value of ε varied from {0.5, 0.7, 0.9}. The parameters of the network were ε and N_p.

The radius of the throat was determined from the radii of the pore bodies adjacent to the addressed throat (R_p1 and R_p2) as

$$r_t=C\min (R_{p1},R_{p2})$$

(16)

where the value of constriction factor C (-) was randomly selected from {0.80, 0.85, 0.90}.

2.3.2. Simulations in Small Networks

The drainage algorithm was used to simulate the sand retention curve to match it with the measured retention curve. The results of these backward simulations were two networks (parameters), one for the θ = 0 and the other for the θ = π/4. Hundreds of network realizations were generated for each network type and each contact angle.

The pressure head boundary was established at the bottom of the network, while air entered from the top. Saturated hydraulic conductivity was determined for each network (K_s,s). Six measured points of the retention curve (suction pressure at: 15, 25, 30, 35, 40, and 50 hPa) were used for the comparison of the measured and simulated curves. The objective function for parameter optimization φ (-) was:

$$\phi =\sum _{i=1}^6{({VWC}_s\left(h_i\right)-{VWC}_m\left(h_i\right))}^2+2{\left(p_s-p_m\right)}^2+2{\left(K_{s,s}-K_{s,m}\right)}^2$$

(17)

where VWC_s (-) and VWC_m (-) are simulated and measured volumetric water content, respectively, each corresponding to specific pressure head h_i. Parameters p_s (-) and p_m (-) are porosities, obtained from the simulation and from the measurement. The K_s,s and K_s,m are the (fully) saturated hydraulic conductivities, obtained from simulation and measurement.

The above-described network parameterization steps were carried out for two contact angles on all three considered network types. First, we parameterized three small network models using a contact angle θ = 0. Subsequently, we repeated the parameterization for a contact angle of θ = π/4, yielding another set of three small network models. After obtaining the network parameters, the retention curve based on the parameterized pore network models was calculated for the range of 5 to 50 hPa in small pressure head steps (≈2 hPa). A total of 20 unique network replicates were generated using the best-fitting parameters for each contact angle for the cubic regular approach. The retention curve was simulated for each network.

Furthermore, the relationships between the residual gas saturation of the sample, S_g,r (-), and initial gas saturation, S_g,init (-), were calculated and plotted for 5 unique networks with best-fitting parameter replicates. The values of S_g,r resp. S_g,init were obtained from the values of ω_r resp. ω_init by dividing porosity. The suction pressure values were {25, 28, 30, 32, 35, 40, 45, 50} hPa for networks generated for θ = 0 and {20, 21, 22, 25, 26, 27, 28, 30, 35, 40, 50} hPa for the networks generated for θ = π/4. First, a drainage simulation was performed with the lower pressure head boundary condition set at the bottom of the network and the air inlet located at the top of the network. The value of ω_init (and the derived value of S_g,init) was a result of this step. Subsequently, the imbibition algorithm with inflow from the bottom of the network was used. Then the trapping algorithm was used to obtain the value of ω_r (and the derived value of S_g,r). This was repeated for all pressure head values.

2.3.3. Large Network Generation

The generation of the network was aided by additional information provided by µCT of the sample. The µCT images showed that large globular air-filled pores were formed in the upper part of the sample during drainage–imbibition cycling probably by air bubbles, while the lower part remained mostly unaffected by bubble formation. Thus, the network parameters in the bottom most part (sub-network 1) were kept the same as in simulation in a small network (Section 2.3.2), while the pore network parameters were modified in the rest of the network. The resulting large pore network consisted of nine sub-networks, small pore networks stacked on top of each other, where pore-network elements size increased with increasing elevation (Figure 3a). These sub-networks represented sub-volumes of the sub-volumes of the µCT image of the sample. The size of the largest bubbles, and thus also the largest pore bodies, was determined from the µCT images for each sub-volume. This information was also used to build the pore network. An example of a cross section with region of interest in the upper half of the sand sample (height is from 25 to 50 mm) is shown in Figure 3a.

The network parameters R_min, R_max and spacing were kept constant within each individual sub-network but increased linearly with increasing distance of the sub-network above the base of the large network, in the same fashion as in the acquired sub-volumes of µCT images. The ratio between these parameters remained constant. The parameters σ and δ were kept constant throughout the pore network. Sub-networks 2–5 were more densely spaced to achieve a smoother transition between sub-network 1 and sub-network 6.

The histogram of µCT determined air-filled pores radii in the upper half of the sample is shown in Figure 3b. Air-filled pores radii were determined from the segmented µCT image as pore volumes with the assumption of the spherical shape of the pores. The maximal pore body radius was specified as 600 µm with the assumption that 1% represented outliers (Figure 3b). The radius of the pore throat (radius of the inscribed circle) was determined.

$$R_t=C^{*}\min (R_{p1},R_{p2})$$

(18)

where the constriction factor of the large network C* (-) was kept constant and R_p1, R_p2 were the radii of the pore bodies adjacent to the addressed pore throat.

2.3.4. Simulations in a Large Network

The simulation of the experiment involved two steps. The first step was the drainage simulation, performed on an initially fully saturated sample using a drainage algorithm with a boundary condition of −30 hPa pressure at the bottom.

The second step involved simulation of the imbibition using an algorithm with the input at the bottom. The order in which the pores were filled obtained from the imbibition algorithm was used as input for the entrapment algorithm, which eventually provided the set of pores occupied with trapped air. The air leak outlet was positioned at the top of the network and residual gas saturation was determined. Finally, the last algorithm was utilized to determine the satisfied hydraulic conductivity.

The result of the simulation was a histogram of the sizes of the air groups trapped. In addition, the frequency of the clusters that correspond to the ganglia was determined. Whether the given cluster represented a ganglion was decided using a topological parameter, the Euler number, which was calculated using the equation described by [31].

Figure 4 shows the schema of simulation processes for both small and large networks.

3. Results and Discussion

3.1. Generated Networks

The comparison of the results obtained using different networks for the numerical experiment in small networks is shown in Figure 5. The saturated hydraulic conductivity, porosity, and objective function were obtained for each of the networks. Each point in Figure 5 corresponds to the results obtained from one network. The points in Figure 5c,f form three separate groups that reflect the value of ε, shown in Figure 5c. Additionally, smaller subgroups, which correspond to different values of C, are visible in each group.

The regular cubic and Gabriel method-generated networks led to simulated hydraulic conductivities similar to the measured one. However, the use of a regular cubic approach led to the lowest value of the objective function. In the case of contact angle 0 the network generation parameters that led to the best fit of the pore network model to the measured retention curve and K_s were R_min = 6.4 µm, R_max = 86.6 µm, σ = 1.3931, and δ = 1.0308 which leads to simulated porosity 0.321 and K_s = 0.045 cm∙s⁻¹. In the case of the contact angle π/4 the optimal parameters were R_min = 14.6 µm, R_max = 86.1 µm, σ = 1.0071, and δ = 1.0073 leading to simulated porosity 0.353 and K_s = 0.058 cm∙s⁻¹. The optimized value of parameter C was 0.9 in both cases.

Network parameter sets were further used for generating sub-network 1 for the large pore network for the forward simulation. The values of the network parameters are listed in

Table 1. Parameters of 9 sub-networks of connected cubic networks. Sub-network 1 is the bottom layer and sub-network 9 is the top layer of the large network. The remaining parameters were C* = 0.93, σ = 1.3931, and δ = 1.0308 for θ = 0, and C* = 0.93, σ = 1.0071, and δ = 1.0073 for θ = π/4.

Sub-Network	Rmin		Rmax		Spacing	Height	Width
θ =	0	π/4	0	π/4	0/π/4	0/π/4	0/π/4
	(μm)		(μm)		(μm)	(pores)	(pores)
	Simulation in a small network
-	6.4	14.6	86.6	86.1	180	60	30
	Simulation in a large network
1	6.4	14.6	86.6	86.1	180	75	50
2	6.8	15.4	91.4	90.9	190	4	47
3	7.1	16.2	96.3	95.7	200	4	45
4	7.5	17.0	101.1	100.5	210	4	42
5	7.8	17.8	105.8	105.2	220	4	40
6	11.4	26.0	154.0	153.1	320	25	28
7	14.9	34.1	202.1	200.9	420	19	21
8	18.5	42.2	250.1	248.7	520	15	17
9	22.0	50.3	298.3	296.6	620	15	14

. The value of C* equal to 0.93 (Equation (18)) led to porosities of the whole sample equal to 0.35 for θ = π/4 and 0.32 for θ = 0. The pore network for this numerical experiment was constructed from 253,146 pore bodies and 744,360 pore throats.

Figure 6 shows examples of networks generated for backward and forward simulations with parameters used for numerical simulations. The small network, with sizes 5.4 × 5.4 × 10.8 mm³, contains 54,000 pores bodies and the large network, with sizes 9 × 9 × 50 mm³, contains over 250,000 pores.

3.2. Results of the Simulation in a Small Network

A comparison of the simulated and measured retention curves is shown in Figure 7. A better agreement was obtained for θ = π/4 with porosity close to the measured value.

Furthermore, the relationships between the residual gas saturation of the satiated sample, S_g,r (-), and the initial gas saturation, S_g,init (-), are shown in Figure 8a. The maximum of S_g,r was equal to 0.45 for θ = 0 and 0.52 for θ = π/4. The dotted lines show the fit of the quadratic model by [57]. There is an approximately linear growth of S_g,r(S_g,init) for lower values of initial gas saturation to the value over 0.5. Furthermore, the growth rate of the amount of entrapped air decreases significantly as the initial gas saturation increases. The general shape of the relationship is in good agreement with the findings of Pentland [58] and Kazemi [59]. Figure 8b shows the difference between the air distributions in the networks created for both contact angles for pressure head equal to 30 hPa.

3.3. Results of the Forward Simulation in a Large Network

The results are summarized in

Table 2. Summarized results for simulations in a large network. The results were obtained from five networks for each contact angle. SD indicates the standard deviation. The mean and SD values were calculated from these networks.

θ		K(ωg,r)	ωg,init	ωg,r	AR	SD(AR)
		cm∙s−1	cm3∙cm−3	cm3∙cm−3	-	-
0	mean	0.035	0.24	0.12	1.88	1.59
	SD	1.5 × 10−4	1.0 × 10−3	1.5 × 10−3	2.9 × 10−3	6.0 × 10−3
π/4	mean	0.035	0.34	0.18	1.55	0.82
	SD	3.8 × 10−4	1.1 × 10−3	1.3 × 10−3	7.3 × 10−4	1.1 × 10−3

, where SD indicates the standard deviation. The mean and SD in rows were calculated from five different networks. The mean value of the volumetric air content ω within the entire sample determined by the simulation with the network generated for the contact angle π/4 is 0.34 after drainage and 0.18 after spontaneous imbibition. The mean value of ω in the case with a network generated for the contact angle 0 is 0.24 after drainage and 0.14 after spontaneous imbibition. The ω determined gravimetrically was 0.244 after drainage and 0.040 after the infiltration experiment. A linear approximation was used instead of a specific measurement. A linear approximation with the equation S_g,r = 0.36∙S_g,init can be used as it was determined in [55]. Then, the residual saturation, S_g,r, the estimate would be 0.23, which corresponds to the value of ω_r equal to 0.083. This way, we can say that the values of ω in the case with the network generated for θ = 0 are relatively accurate.

The air distribution in the sample, represented as volumetric air content (ω), is shown in Figure 9. Figure 9a presents the simulation results, while Figure 9b shows the adjusted results where, to account for small pores missed by µCT, gas-filled pores and throats smaller than 170 µm were excluded. This air content was used to compare with µCT results. Figure 9a indicates a relatively uniform air distribution across the height, unlike µCT images, where higher ω_r values appear in the upper half. This difference is due to larger, easily detectable air bubbles in the upper part, while smaller bubbles in the lower part were not well captured by µCT. The pore network model, however, reflects the total air content.

The effect that can influence the value of ω_r is the mobility of entrapped air bubbles that was not considered in the model. It is reasonable to assume that the ganglia could be mobile and capable of escaping the network even after capture. Figure 10 shows the distribution of the air clusters by size; approximately 50% of the number of individual air clusters are single-pore bubbles. However, the total volume of single-pore bubbles represents only approximately one hundredth of the sample volume.

The overestimated simulated value of ω_r cannot be attributed to the neglected snap-off mechanism because the use of the mechanism generally leads to higher ω_r values [17] and the snap-off mechanism is suppressed due to the low aspect ratio value [60], where the aspect ratio means the relationship between the pore body radius and pore throat radius. The average aspect ratio for the network generated for θ = π/4 is 1.55 with a standard deviation 0.82 and for the network generated for θ = 0 the average aspect ratio is 1.88 with a standard deviation of 1.59.

Our simulation results align well with the microtomography experimental work on sandstone published by Herring et al. [61], who also observed an increase in residual non-wetting phase content with rising contact angles (from strongly water-wet to weakly water-wet conditions). The increase in residual non-wetting phase is attributed to the dominance of capillary over viscous and gravitational forces within the contact angle range of 0 to π/4. This dominance shapes the drainage patterns, resulting in a higher content of non-wetting fluid as compared to the initial condition. During imbibition, a significant amount of this non-wetting fluid becomes entrapped due to capillary forces acting on the non-uniform wetting front.

The value of $K_s$ was obtained from the simulation in a small network where K_s = 0.058 cm∙s⁻¹ for the network generated for θ = π/4 and K_s = 0.045 cm∙s⁻¹ for the network generated for θ = 0. The low K_s value in the latter case is probably related to a lower porosity. The satiated hydraulic conductivity for a pressure head of 30 hPa was determined from simulation in a large network. The value is 0.035 cm·s⁻¹ for the network generated for θ = π/4 and 0.035 cm∙s⁻¹ for the network generated for θ = 0. The measured value of K_s was equal to 0.060 cm∙s⁻¹. In our previous study [56], we experimentally determined K(ω_max) = 0.034 cm·s⁻¹ at a maximum obtained value of ω equal to 0.11. It can be seen in Figure 11 that the hydraulic conductivity for network generated for θ = 0 was close to the measured value. However, saturated hydraulic conductivity was underestimated in the case when the θ = 0 which is a consequence of lower porosity (see Figure 7).

4. Conclusions

The simple pore-network model based on the OpenPNM platform was presented to simulate physical experiments with air trapping conducted on a sand sample. In the first numerical experiment, with a small uniform network, the network parameters were optimized to fit the measured retention curve and saturated hydraulic conductivity. The regular cubic network led to the best agreement with the measurement.

In the second case, the infiltration experiment performed on the sand column that was drained and subsequently filled with water from the bottom was simulated in a large non-uniform network composed of nine sub-networks. This simulation was performed without any further parameter optimization other than considering the pore network element size variability obtained from µCT images.

The model overestimated the amount of entrapped air for network generated for the contact angle π/4. The saturated hydraulic conductivity was determined with reasonable accuracy.

The case of network generated for contact angle of 0 provided relatively good agreement between the model and the measured content of air within the sample, but the saturated hydraulic conductivity was underestimated, which is probably a consequence of the presence of smaller pores and overall lower porosity obtained for this contact angle than the measured porosity. However, the satiated hydraulic conductivity reached a value similar to the measured value.

Limitations of this study are the assumptions of absence of entrapped gas buoyancy and snap-off mechanism. However, the latter only becomes relevant for simulations with larger contact angles than the ones considered here. Our study has demonstrated that fast quasi-static algorithms and a network that was generated using independently measured hydraulic properties together with additional information from µCT successfully capture trends of air entrapment during imbibition of sand in relation to the initial air saturation.