Dynamic Pore-Scale Model of Drainage in Granular Porous Media: The Pore-Unit Assembly Method.

Dynamics of drainage is analyzed for packings of spheres, using numerical experiments. For this purpose, a dynamic pore-scale model was developed to simulate water flow during drainage. The pore space inside a packing of spheres was extracted using regular triangulation, resulting in an assembly of grain-based tetrahedra. Then, pore units were constructed by identifying and merging tetrahedra that belong to the same pore, resulting in an assembly of pore units. Each pore unit was approximated by a volume-equivalent regular shape (e.g., cube and octahedron), for which a local capillary pressure-saturation relationship was obtained. To simulate unsaturated flow, a pore-scale version of IMPES (implicit pressure solver and explicit saturation update) was employed in order to calculate pressure and saturation distributions as a function of time for the assembly of pore units. To test the dynamic model, it was used on a packing of spheres to reproduce the corresponding measured quasi-static capillary pressure-saturation curve for a sand packing. Calculations were done for a packing of spheres with the same grain size distribution and porosity as the sand. We obtained good agreement, which confirmed the ability of the dynamic code to accurately describe drainage under low flow rates. Simulations of dynamic drainage revealed that drainage occurred in the form of finger-like infiltration of air into the pore space, caused by heterogeneities in the pore structure. During the finger-like infiltration, the pressure difference between air and water was found to be significantly higher than the capillary pressure. Furthermore, we tested the effects of the averaging, boundary conditions, domain size, and viscosity on the dynamic flow behavior. Finally, the dynamic coefficient was determined and compared to experimental data.

Introduction

The capillarity phenomenon is fundamental to the flow of fluids in porous materials. For example, the rise of water into a hydrophilic porous material is due to capillary action. Capillarity is a key process in many natural and industrial processes, for example, not only in unsaturated soils (Hassanizadeh et al., 2002; Likos & Jaafar, 2013; Nikooee et al., 2013; Nuth & Laloui, 2008; Oh & Lu, 2014), carbon sequestration, oil and reservoir engineering but also in water transport in thin porous media such as tissues (e.g., Sun et al., 2015), fuel cells (e.g., Qin, 2015), and hygienic products (Diersch et al., 2010).

Capillarity is often described as the relation between capillary pressure and saturation under equilibrium conditions. The capillary pressure‐saturation relationship of air and water is often experimentally characterized where experiments may take hours to weeks. However, in many applications, flow of water may occur in much shorter time scales. For example, during dynamic drainage experiments by Camps‐Roach et al. (2010) on a 20 cm long column, it took just 30 min for the majority of drainage to occur. Moreover, imbibition of ink in paper occurs within a second (Kettle et al., 2010). Relatively fast drainage and imbibition processes give rise to the concept of dynamic capillarity effect. Under fast transient drainage or imbibition, the difference in fluid pressures is not equal to capillary pressure anymore. Many examples exist in the soil physics and two‐phase flow literature that describe the dynamic capillarity effect (see e.g., Bottero et al., 2011; DiCarlo, 2005; Topp et al., 1967). For comprehensive reviews of dynamic capillarity effects, see Hassanizadeh et al. (2002) and Diamantopoulos and Durner (2012).

Here we would like to employ an equation based thermodynamic considerations by Hassanizadeh and Gray (1993), who have conjectured the following dynamic capillarity relationship: where and denote the macroscopic pressure of air and water, respectively, denotes the capillary pressure, is a dynamic coefficient, and is the rate of change of saturation with time. Experiments have been conducted to explore equation (1) and to determine the value of . It has been shown that for transient drainage, at any given saturation, is larger than ; the larger the imposed pressure difference the more negative and thus the larger the dynamic effect (e.g., Bottero et al., 2011; Camps‐Roach et al., 2010). During dynamic imbibition, is smaller than (see e.g., Bottero et al., 2006). O'Carroll et al. (2005) conducted Multi‐Step‐Outflow (MSO) experiments on perchloroethylene (PCE) and water, to investigate the effect of dynamic capillarity. In MSO, parameters are estimated by inverse modeling of experimental data. They found that the goodness of the fit of inverse modeling of the experiments improved significantly when the dynamic effect was included.

Dynamic capillarity effects have also been studied using pore‐scale models. Pore‐scale models can be divided into two types, namely, (i) models that discretize the pore space into pore bodies and pore throats, such as pore‐network models (Joekar‐Niasar et al., 2010; Qin, 2015; Thompson, 2002) and pore‐unit assembly models (Mason & Mellor, 1995) and (ii) models that employ direct, or fully resolved, simulations of flow within the pore geometry, such as Lattice‐Boltzman (Vogel et al., 2005), smooth particle hydrodynamics (e.g., Kunz et al., 2016), and computational fluid dynamic (CFD) simulations (e.g., Ferrari & Lunati, 2012; Rabbani et al., 2017). Direct simulations of two‐phase flow are typically applied to a small number of pores or limited to two‐dimensional domains, while pore‐network and pore‐unit assembly models can handle a very large number of pores in two and three dimensions. On the other hand, pore‐network and pore‐unit assembly models employ very simplified pore geometries and flow equations. Pore‐network models have been mainly employed to study the effect of various parameters on two‐phase flow. For example, Mogensen and Stenby (1998) developed a single‐pressure algorithm to simulate dynamic imbibition. They studied the effects of contact angle and capillary number on the residual oil saturation after a water flooding event. Thompson (2002) developed a comprehensive pore‐network model for dynamic imbibition to study water movement in fibrous porous materials, where both the nonwetting and wetting pressures were solved implicitly and saturation was updated explicitly. They included corner flow, capillary pressure in a pore body and flow dynamics of two viscous fluids. Joekar‐Niasar et al. (2010) developed a pore‐network model of two‐phase flow based on a two‐pressure algorithm similar to that in Thompson (2002), but saturation was updated semi‐implicit in order to enhance numerical stability at low flow rates. To evaluate their model, they conducted simulations on a three‐dimensional lattice of cubic pore bodies, with a coordination number of six (i.e., each pore body was connected to six other pore bodies). For a cubic pore body, the local capillary pressure‐saturation curve was analytically derived. The model has been used to study the effect of various parameters, such as viscosity ratio and capillary number, on the dynamic capillarity effect (Joekar‐Niasar & Hassanizadeh, 2011).

However, pore‐network models assume the pore geometry to be rigid, despite of the common occurrence of coupled flow and deformation. For example, a deformation can be induced by applying a confining stress to a porous sample, causing the pore space to decrease in size, which can invoke a change in water saturation. In addition, soils can deform during a change of water saturation (Kharaghani et al., 2011; Lins & Schanz, 2005); ink penetration in paper can induce swelling of the hydrophilic binders (Rosenholm, 2015); and fluid distribution in foods occurs alongside changes in the polymeric chains (Takhar, 2014). Furthermore, in swelling materials, fluid absorption by the solid phase causes deformation of the pore space and thus causes fluid redistribution, which is the case in super absorbent polymer particle beds in hygienic products (Diersch et al., 2010).

An appropriate pore‐scale model that can include deformation and flow is the development based on coupling of pore‐unit assembly method (PUA) with the Discrete Element Method (DEM). DEM is a particle‐based model that is capable of simulating particle‐particle interactions during deformation of a granular material (Cundall & Strack, 1979). In the PUA method, the pore geometry of spherical particle packings can be extracted using a regular triangulation; four neighboring particles form the vertices of a tetrahedron, which encloses a so‐called pore unit. PUA has been applied on Finney packings to construct imbibition curves for varying contact angles (Gladkikh & Bryant, 2005). More recently, PUA has been coupled to the Discrete Element Method (DEM) to enable simulation of hydromechanical coupling in saturated granular materials (Catalano et al., 2014; Chareyre et al., 2012). PUA in combination with DEM has also been used to construct equilibrium capillary pressure‐saturation curves for packing of sand and glass‐beads. First, the size distribution of spheres and the porosity of the packing were matched to those of a granular material. Next, the equilibrium capillary pressure‐saturation curve was constructed, see for example, Yuan et al. (2016), Sweijen et al. (2017a, 2016), and Mahmoodlu et al. (2016). In those studies, the capillary pressure‐saturation computations were based on binary saturation values (i.e., a pore unit is either empty or fully saturated). They also calculated entry pressures for pore throats, which were based on pore throat radii following the Mayer‐Stowe and Princen method or the hemisphere method. See Sweijen et al. (2017a) for a comparison between the two methods.

Based on the coupling of PUA with DEM, we introduce in this work a pore‐scale model for unsaturated flow to simulate dynamic drainage in a packing of spheres. We have developed an algorithm for simplifying the complex pore geometry of pore units by fitting regular shapes into pore units. This algorithm accounts for varying coordination numbers of pore units. By simplifying the pore geometry, a local capillary pressure‐saturation relationship for each pore unit could be determined, enabling dynamic flow computations with pore units being partially saturated, rather than using a binary saturation. Dynamic drainage is simulated by employing a single‐pressure algorithm to solve for water pressure, based on a finite difference scheme. We refer to this algorithm as a pore‐scale version of IMPES (implicit pressure solver and explicit saturation update).

The aims of this research are (i) to develop a dynamic pore‐scale model that captures the main geometrical information from packings of spheres and to simulate unsaturated flow during dynamic drainage, (ii) to test and verify the model by reproducing the quasi‐static capillary pressure‐saturation curve, and (iii) to investigate dynamic capillary effects in packings of spheres and to study the effects of averaging, boundary conditions, pore‐scale processes, domain size, and viscosity and by comparing the dynamic coefficient in simulations to that of experiments.

The outline of this study is as follows. First, we explain the algorithm to fit regular shapes into pore units. Then, we evaluate the single‐pressure solver and the local rules for two‐phase flow. Finally, the model is tested by simulating dynamic drainage experiments.

Numerical Model

Packings of spheres are generated using the Discrete Element Method (DEM), which was first introduced by Cundall and Strack (1979). In this study, we employ open‐source software Yade‐DEM, which is an extendable and three‐dimensional model (Šmilauer et al., 2015). Packings of spheres are generated for a predefined porosity value and particle size distribution, following the methodology of Chareyre et al. (2002). From the particle packing, the pore structure is extracted using a regular Delaunay triangulation, resulting in an assembly of tetrahedra. Afterwards, the pore space is simplified by replacing the pore geometry inside each tetrahedron by that of a regular shape, which is often referred to as a platonic solid (e.g., tetrahedron, cube, and octahedron). Thus, each regular shape represents a pore body, which we refer to as a pore unit. Finally, we obtain a pore‐unit assembly. In the following sections, we explain the methodology of generating an assembly of regular shapes.

Representing Pore Units by Regular Shapes

The pore space from a spherical particle packing is extracted using a regular triangulation, such that the pore space is subdivided into tetrahedra (see Figure 1a). Each tetrahedron is enclosed by four spheres, where the vertices of the tetrahedron are located at the centers of spheres. One tetrahedron is connected to four other tetrahedra via narrow openings at the facet of the tetrahedra. These facets thus represent zero‐volume pore throats that offer resistance to flow (see Figure 1a). The tetrahedra are referred to as grain‐based tetrahedra, such that the whole pore space is represented by an assembly of grain‐based tetrahedra. In this research, we employ the regular triangulation which accounts for differences in particle sizes unlike the Delaunay triangulation. Regular triangulation was implemented into Yade‐DEM by Chareyre et al. (2012).

Figure 1

Two‐dimensional illustration of the algorithm to find regular shapes to represent pore units: (a) grain‐based tetrahedra from triangulation, (b) pore‐based tetrahedra having the same volume as the void space inside the grain‐based tetrahedra, note the indication of and , and (c) regular shapes that were found after merging of pore‐based tetrahedra.

The pore space inside an individual grain‐based tetrahedron has a complex shape that is different for each grain‐based tetrahedron. To enable unsaturated flow simulations, the pore space inside each grain‐based tetrahedron is simplified by representing it by a pore‐based tetrahedron, which is placed with its vertices into the pore throats and which has a volume equal to the pore volume of a grain‐based tetrahedron (see Figure 1b). The size of a pore throat is characterized by radius of the inscribed circle between three particles that make up a facet of a grain‐based tetrahedron, while the size of the pore space is described by radius of the inscribed sphere inside the pore‐based tetrahedron (see Figure 1b).

In some cases, the pore throat radius of pore throat ij can be larger than the radius and/or of the inscribed sphere of pore‐based tetrahedra i and j, following Al‐Raoush et al. (2003) and Sweijen et al. (2016). Per definition, the pore throat size should be smaller than the pore size. However, this is not the case for larger pores that contain two or more grain‐based tetrahedra. This happens in the case of pore units that are surrounded by more than four spheres. To alleviate this issue, we apply a merging algorithm with the aim of identifying grain‐based tetrahedra that belong to the same pore unit. This algorithm is explained further in section 2.2.

When two grain‐based tetrahedra are merged, another regular shape is chosen to represent the pore space inside the merged tetrahedra such that the vertices of a regular shape are fitted into the pore throats, see for example, Figure 1c. Each regular shape represents a pore unit and has a void volume and adjacent pore throats. For example, a set of merged grain‐based tetrahedra that has six adjacent pore throats is represented by an octahedron, which has six vertices. Here we do not consider the exact relative location of those pore throats within the pore unit, but we merely assume them to match the location of pore throats (i.e., vertices) of a regular shape.

Merging Algorithm

The aim of the merging algorithm is to identify grain‐based tetrahedra that belong to the same pore unit, while considering the regular shape that is chosen to represent the merged grain‐based tetrahedra. Initially, we replace the pore space of each grain‐based tetrahedron by a pore‐based tetrahedron, having four vertices (Figure 1b). To check whether the pore throats are smaller than the radius of pore‐based tetrahedra, we compute the radius of inscribed sphere of each pore‐based tetrahedron. The radius of the inscribed sphere ( ) of a pore‐based tetrahedron (in fact any regular shape) is given by in which is a shape factor, whose value is different for different regular shapes, and is the volume of the regular shape (and thus of the pore unit). Values of are given in Table 1. In the case that one (or both) of the inscribed spheres is smaller than the pore throats connecting them, we merge the two spheres into one big pore unit. In order to avoid very small differences between and or playing a role, we employ a factor to control the strictness of merging. Thus, when one of the following criteria are not satisfied, we merge pore units i and j: where is a user‐defined parameter less than unity. In this research was set to 0.90. The method of merging grain‐based tetrahedra is similar to the work by Al‐Raoush et al. (2003), who merged grain‐based tetrahedra together, based on their inscribed spheres. They showed that merging of tetrahedra is an efficient way of generating a pore network but yields slightly different pore statistics than image analysing techniques provide. Typically, the pore body radii and pore throat radii are larger when using merging of tetrahedra, compared to radii from image analyzing techniques. Bakke and Øren (1997) also employed a technique similar to merging; they used an imaging method (i.e., thinning) to draw Voronoi polyhedral around grains, from which the vertices were combined until the inscribed sphere with its center on a vertex was always bigger than a pore throat.

Table 1

Geometrical Constants for Various Regular Shapes

Name	N_i	χ	λ	θ	N_edges	κ
Tetrahedron	4	0.42	2.04	1.23	6	3.87
Octahedron	6	0.53	1.62	1.91	12	8.71
Cube	8	0.50	1.00	12π	8	6.83
Octahedron + cube	10	0.44	1.19	12π and 1.91	12	5.15
Icosahedron	12	0.58	0.77	2.41	30	24.11
Dodecahedron	20	0.57	0.51	2.03	30	22.87

The merging procedure starts with a pore‐based tetrahedron that has the largest value of ; it will be merged to its neighboring pore‐based tetrahedron. After merging, a new regular shape is fitted into the merged tetrahedra, and are subsequently updated. Then, a pore unit with the next largest value of is merged to its neighbor; this is repeated until conditions in equation (3) are satisfied for all pore units and throats. To avoid long chains of merged pore units, the maximum coordination number is set to 20, as otherwise pore units with very large coordination numbers may arise. Finally, a small number of unresolved pore throats remain (typically <1% of the pore throats). For those, we change their to . Finally, we ensure that only one pore throat exists between two pore units, as it may occur that more throats exist. If more throats exist, the throats are combined to one throat by taking an arithmetic average of pore throat properties (e.g., pore throat radii). The final set of merged tetrahedra is referred to as an assembly of pore units.

Benchmark testing revealed that our merging algorithm works well. For example, we generated different packings for which we knew in advance how the pore structure should be. First, a packing of spheres was generated having tetrahedral pore units (i.e., pore units have a coordination number of four). Second, a packing of particles was generated having octahedral pore units (i.e., pore units have a coordination number of eight), see Figure 2a. Finally, spheres were put in a triangular prism (i.e., pore units have a coordination number of 5), see Figure 2b.

Figure 2

Three dimensional illustrations of packings of spheres from which the pore space is represented by a regular shape that is known prior to applying our merging algorithm: (a) octahedral assembly of spheres, whose pore unit has a coordination number of 8 and (b) a triangular prism assembly of spheres, whose pore unit has a coordination number of 5.

Pore‐Scale Version of IMPES

To simulate flow of water in a pore‐unit assembly, we employ a pore‐scale algorithm of a finite difference scheme to solve implicitly for pressure and explicitly for saturation, which is often referred to as IMPES (see e.g., Chen et al., 2004; Huber & Helmig, 1999). Pressure is only solved for water, because air is assumed to be infinitely mobile compared to water. In other words, the two‐phase flow of air and water is characterized by a viscosity ratio of 0.02 (ratio of air over water viscosity), hence its pressure is assumed to be constant at . This assumption was also employed by Kibbey et al. (2016) to study the effect of drainage rate on residual saturation in a pore‐network model. Let us consider a pore unit having a water saturation , water pressure and air saturation . Per definition, we have

In pore unit , a local capillary pressure is defined that under no‐flow conditions is related as follows:

When solving for water pressure, is not strictly imposed as the capillary pressure may be different than . The volumetric flux of water though pore throat ( ) is assumed to be linearly proportional to the pressure gradient (e.g., Chareyre et al., 2012; Joekar‐Niasar et al., 2010): in which is a pore throat conductivity. For saturated pore units, a volume balance is written as follows: which yields the following equation, using equation (6):

For partially saturated pore units, a different set of equations is employed for computing the water pressure. In what follows, we consider the air pressure to be equal to zero such that equation (5) reduces to . For partially saturated pore units, the volume balance is given by

Equations (8) and (9) are can be written into a discretized form (see Appendix A) as follows: where coefficients , , and are defined in Table 2. Equation (10) provides a linear set of equations that is solved for water pressure, using implicit solver Sparse supernodal LU factorization in Eigen libraries (Eigen, 2017). Once the water pressure is solved, the saturation is updated explicitly following:

Table 2

Parameters to Solve the Linear System in Equation 7.10

	Name	a	b	c
siw=1	Saturated	∑j=1N_ikijw	kijw	0
0<siw<1	Partially saturated	∑j=0N_ikijw+V_iΔt∂siw∂piw	kijw	V_iΔt∂siw∂piwpit

Local Capillary Pressure Saturation Curve

A set of equations is derived to describe the local capillary pressure‐saturation relationship for partially saturated pore units, following the work by Joekar‐Niasar et al. (2010). See Appendix B for detailed and illustrated information. We define a threshold saturation, , for pore unit i that corresponds to a saturation at which air fills the inscribed sphere of the pore unit (Joekar‐Niasar et al., 2010; Thompson, 2002):

When the saturation in a pore unit is smaller than , water resides in corners and edges of a pore unit. Thus, following the work by Joekar‐Niasar et al. (2010), we can write where and is the volume of water residing in corners and edges, respectively. The volume of water in the corners can be computed, using equation (2): where is the radius of curvature of the air‐water interface, which is given by Young's Laplace equation , with being the interfacial tension of air‐water (0.072 N m−1). The volume of water along the edges is determined as follows. Consider an edge that has a dihedral angle . In the plane that is perpendicular to the length of an edge, the area occupied by water is given by, following Joekar‐Niasar et al. (2010):

Note that the radius of curvature in an edge is given by , because the second principal curvature is assumed to be flat. The length of the edges inside a regular shape is given by where is a geometrical parameter that determines the length of one edge as a function of the volume of a regular shape (see Table 1 for its values). We substitute equations (15) and (16) into , where is the number of edges in the regular shape. The total water volume in a pore unit, , is given by

Thus, water saturation inside pore unit can be written as where and are geometrical parameters defined by

Equations (18) and (19) are fitted by the following equation, following Joekar‐Niasar et al. (2010): where is a geometrical constant; for more information about the fitting of equation (20), see Appendix C.

Geometrical Parameters for a Range of Coordination Numbers

The pore space does contain pore units that can be fitted not only with regular shapes but also with irregular shapes that have other coordination numbers than regular shapes. To overcome this problem, we fit an equation to relevant geometrical properties, namely: , and , as a function of ; see Figure 3. When we can replace a pore unit by a regular shape, we will use the values of , and that are derived directly for that regular shape, otherwise we will use fitting equations to find values for those parameters. The fitted equations should yield values that correspond to a sphere for converging to infinity. This was checked and all equations given in Table 3 satisfy this test.

Figure 3

Pore‐scale parameters needed for the local capillary pressure‐saturation relationship for partially saturated pore units. Symbols indicate the values for various regular shapes, given in Table 1, and the dashed lines indicate the fitting by equations given in Table 3.

Table 3

Fitted Equations for Geometrical Parameter of Regular Shapes, Which Are Used to Estimate Geometrical Parameters of Irregular Shapes

Parameter	Fitting equation
κ	1.32N_i
λ	7.12N_i^−0.89
χ	34π31−e^−0.170N_i

Local Drainage

The entry pressure of pore throat , , defines the air‐water pressure difference at which air can invade a water‐saturated pore unit. To estimate the entry pressure, we assume that the air‐water interface has the shape of a hemisphere. The entry pressure is given by the Young's Laplace equation, assuming a perfectly water wet solid phase:

Equation (21) is an approximation for the entry pressure for sake of computational time, but of course a more rigorous method such as the level‐set method would yield more accurate entry pressure values (Prodanović & Bryant, 2006). During drainage, air can only invade from pore unit to its neighboring pore unit if

Air‐Water Interfacial Area

The air‐water interfacial area within a pore unit ( ) is determined as follows. For pore units with , the shape of the air‐water interface is unknown. Therefore, we assume it to be equal to that of a sphere having the same volume as air, following Joekar‐Niasar et al. (2010). For stable partially saturated pore units ( ), the water is residing in the edges and in the corners, for which a derivation of air‐water interfacial area is given in Appendix D. We approximate the air‐water interface as follows: where is the number of pore throats that are not invaded by air, yet.

Hydraulic Conductivity

The conductivity is assumed to depend on the local capillary pressure in pore units and and thus on the water saturation. Following Chareyre et al. (2012), we assume that can be approximated by in which is the surface area of water that resides inside the pore throat, is a scaling parameter to adjust the overall permeability of the pore‐unit assembly (following Chareyre et al., 2012), is the viscosity of water, is the distance between centers of pore units and , and is a representative radius that is assumed to depend on local capillary pressure. We define as follows: in which is the hydraulic radius of pore throat that is computed for a fully saturated pore throat, see Chareyre et al. (2012) for its definition. Moreover, depends on the local capillary pressure in pore units and such that is given by in which is the area of pore throat . Note that in cases where , we set equal to 0. Finally, we define is an effective curvature in pore throat , which is defined as

Time Step and Stability

Controlling the time step is essential for the numerical model to be accurate and stable. For each unsaturated pore unit, a time step is computed, which is an estimation of the time needed for the next event to occur in that pore unit. Following Joekar‐Niasar et al. (2010), we compute that is required to fully empty or fill a pore unit during local drainage or imbibition, respectively. In addition, we have added a time step that is needed to reach a saturation of to ensure that the local capillary pressure‐saturation behavior was reproduced. We define as follows:

Note, that in determining the absolute value of saturation differences, e.g., , should be larger than a truncation value of 10−6 in order to avoid converging to zero. The minimum saturation that a pore unit can have is set to 10−6, at that saturation, the pore unit was considered to be empty. Then, the smallest value of of all pore units is used and multiplied by a safety factor of 0.2.

The numerical scheme presented above is a pore‐scale version of IMPES, which is a highly nonlinear system of equations that is known to have stability problems for capillary dominated flow (Chen et al., 2004; Thompson, 2002). To stabilize the solution, we evaluate and backward in time. Thus, we linearly extrapolate to the next time step. A similar approach was employed by Joekar‐Niasar et al. (2010), who linearly extrapolated the gradient of capillary pressure over a pore throat, namely , over time.

Simulation Setup

Simulating the Experimental Setup

To test our model, we simulated the fast drainage experiments reported by Zhuang et al. (2017). In those experiments, fast drainage was studied for a filter sand having a particle diameter of 0.0–0.50 mm (with an average diameter of 0.20 mm) and porosity value of 0.39. The sample size in that experiment was 3 × 3 × 2 cm³. To reproduce the filter sand in our model, a packing of 4,000 spheres was generated having the same particle size distribution and porosity value as in experiments, following the approach of Chareyre et al. (2002). The final size of the packing was 3.2 × 3.2 × 3.2 mm3. In this work, the viscosity of water was set to 10−3 Pa s, the surface tension of air‐water to 0.072 N m−1 and the contact angle to zero.

As our model is based on smooth spherical particles, it cannot represent irregular sand grains. So even though we use the grain size distribution of the sand in experiments, the resulting pore space does not have the same pore size distribution as the sand. To be able to match the entry capillary pressure measured in the experiments (which has a link to pore sizes), we have modified the grain size distribution. We have multiplied the grain size distribution of the spheres by a factor of 1.28 as a first order modification. Consequently, results compared very well with measured capillary pressure‐saturation curve (see Figure 4). Note that the quasi‐static capillary pressure‐saturation curve was based on the pore‐unit assembly as described in this work, with partially saturated pore units, in contrast to previous work where a binary saturation was assumed (e.g., Sweijen et al., 2016; Yuan et al., 2016).

Figure 4

Quasi‐static capillary pressure‐saturation curve from experiments (symbols) and PUA simulations (solid line). For simulations, the particle radii had to be increased by a factor of 1.28.

The permeability of the packing of spheres was calculated to be 2.11 × 10−10 m2 and was found to be 12.4 times larger than the measured value (1.7 × 10−11 m2). This is a consequence of using spherical particles rather than the actual particle shape (Sweijen et al., 2017b; Torskaya et al., 2014).

Initial and Boundary Conditions

The top boundary of the packing was assumed to be an air reservoir at pressure and the bottom boundary was assumed to be a water reservoir at a constant pressure . The difference between and is referred to as . All other boundaries were considered impermeable. In this research, we used as the reference pressure and was set to . As initial conditions, we assumed all pores to be fully water saturated, except for the first row of pore units along the air reservoir. All pore units were assigned an initial water pressure of – Pa, which was just below the entry pressure of the sample. Thus, along the air boundary, pore units were partially saturated at a saturation value of where was obtained from equation (20). Then, the water reservoir pressure was decreased in order to initiate drainage.

Averaging Procedure

The following three variables have been considered for each pore unit: water pressure, saturation, and air‐water interfacial area. Macroscale equivalents for those three variables were defined as their averages over all pore units. First, the macroscale water saturation is given by where is the total number of pore units in the packing. The volume‐averaged water pressure, , was computed by

The volume‐averaged air pressure is equal to its boundary value , since air pressure is assumed to be constant. The macroscale capillary pressure, , was obtained by averaging the capillary pressure in pore units over their respective air‐water interfacial area (Joekar‐Niasar et al., 2010):

Hereafter, we define as the capillary pressure, while is referred to as the pressure difference between air and water, which under equilibrium conditions is the same as , but not under dynamic conditions. Under dynamic conditions, the relationship between capillary pressure and the air‐water pressure difference is given by equation (1).

Results and Discussion

In the supporting information, sections S1 and S2, we evaluate the minimum number of particles that are required in order to have a representative elementary volume and we verify our simulations. In what follows we discuss the results of our dynamic drainage simulations.

Dynamic Drainage: Formation of Air Fingers

Simulations of dynamic drainage showed the formation of air fingers in an initially water‐saturated domain. The volume‐averaged pressure difference between air and water, , and capillary pressure, , as function of saturation are presented in Figure 5. Furthermore, to illustrate the migration of air, two‐dimensional cross sections are shown in Figure 6 that depict water pressure and water saturation distributions. For these simulations, the value of was set to 7 kPa.

Figure 5

Averaged pressure difference as a function of saturation for dynamic drainage simulations with all sides open to water (denoted by blue lines) and with closed sides (denoted by red lines).

Figure 6

Distribution of water pressure and water saturation obtained using dynamic drainage simulations, where was set to 7 kPa. The sides are closed and the bottom represent a water reservoir while the top represents an air reservoir. Results for three different times, at which the average saturation values are (a) 1.0, (b) 0.84, and (c) 0.11.

Initially, as most pore units were saturated, the water pressure was almost linearly distributed, varying from just less than −3 kPa at the air reservoir to −7 kPa at the water reservoir; see Figure 6a. As drainage progressed, air infiltrated the modeling domain via channels (so‐called air fingers). Figure 6b shows an air finger that has extended from the air boundary to the water reservoir. Consequently, the water pressure of pore units closer to the water reservoir increased significantly, as they were invaded by the air and the local capillary pressure corresponding to a local residual saturation was reached. The water pressure in many pore units was close to −4 kPa, except for pore units near the water boundary, which had lower values of water pressure (–4 to −7 kPa). Toward the end of drainage as the water flow rate decreased and no‐flow conditions were approached, the water pressure in all pore units connected to the water reservoir gradually reached −7 kPa.

The unsaturated flow illustrated in Figure 6 resulted in a peculiar behavior of when plotted against saturation (see Figure 5). Initially, the value of was higher than the value of . Then, the value of decreased until it reached the value of at saturation 0.82, which corresponded to the moment that air fingers reached the water reservoir. Drainage for saturations lower than 0.82 occurred with values of close to the quasi‐static capillary pressure‐saturation curve. The value of deviated from the quasi‐static capillary pressure, which can be related to the approximation of surface area in pore units that is presented in equation (23).

To investigate the effect of boundary conditions, dynamic drainage was also simulated for the case that all sides were open to water, having the same pressure as the water reservoir. Thus, Dirichlet boundary conditions were imposed for the side walls, using a water pressure of – . Figure 7 shows two‐dimensional distributions of the saturation and water pressure at three different saturation values. Similar to results in Figure 6, air fingers were formed in the granular material. However, the air fingers did not infiltrate much farther than the center of the packing, whereas for closed side walls, the air fingers migrated down to the water reservoir. Consequently, the water pressure at saturation 0.80 is much lower (i.e., higher values of pressure differences) for simulations with open side walls compared to simulations with closed side walls (compare Figures 6b and 7b). The difference in between the two simulations is shown in Figure 5, where simulations with open side boundaries show higher values of and a smaller effect of air‐finger formation. Note that the higher values of for open side boundaries compared to that of impermeable side boundaries is merely an artifact of an increased number of pore units along an water reservoir (even though the boundary pore units were not included in obtaining averaged quantities).

Figure 7

Cross sections of dynamic drainage simulations for water pressure and water saturation, where is 7 kPa. The sides and bottom represent a water reservoir while the top boundary represents an air reservoir. Results for three different times, at which the average saturation values are (a) 1.0, (b) 0.82, and (c) 0.32.

Air fingers were a consequence of preferential pathways in the pore structure due to small‐scale heterogeneities. Reconstruction of quasi‐static capillary pressure‐saturation curves showed the same preferential pathways that were present in dynamic flow. The formation of air fingers during primary drainage is a common phenomenon, as for example, shown in the modeling efforts by Aker et al. (1998). They studied the infiltration patterns during primary drainage for viscosity ratios ranging from 10−3 to 102. Their drainage simulation, using a low viscosity ratio, clearly showed finger formations. Their results are compatible with our results, as an air‐water system has a low viscosity ratio of 0.02. Similarly, using their dynamic pore‐network model, Joekar‐Niasar et al. (2010) showed with their dynamic pore‐network model that the pattern of the nonwetting fluid is much more irregular under low values of viscosity ratio (as shown for a viscosity ratio of 0.1).

The Dynamic Coefficient

Effect of Averaging

In Figure 5, the averaged pressure differences, , are much higher than the corresponding capillary pressure values for saturations higher than 0.82. As explained above, this is because the domain initially comprises mostly of saturated pore units, whose number and thus their contribution to decreases as drainage progresses. So, the large values of are a result of the contribution of saturated pore units to and it is not a capillarity effect. However, in column experiments, where only macroscale quantities are observed and water occupancy of pores is unknown, the difference between and is attributed to dynamic capillarity effect.

Thus, obtaining by averaging over a domain that comprises mostly saturated pores is questionable and various studies show that we get pseudo‐dynamic effects due to averaging. This issue has been discussed by Nordbotten et al. (2007, 2008) and it is attributed to the fact that the domains occupied by air and water phases have different centroids, which are also different from the centroid of the averaging domain. A proper way of calculating average values of phase pressures for such domains has been suggested by Nordbotten et al. (2008). Therefore, we have also computed the pressure difference while accounting for location of the centroid of air in the averaging domain, using Korteland et al. (2010, equation (2)). Figure 8 shows the effect of averaging on , where the pressure difference for saturation values higher than 0.82 are slightly smaller (maximal 500 Pa smaller) when using a centroid‐corrected averaging. In addition, we plotted in Figure 8 the values of based on volume averaging of only partially saturated pore units. This curve lacks the enhanced values of for saturation values higher than 0.82, because the water pressure inside a partially saturated pore unit is close to the local capillary pressure (following equation (5)). By only averaging over partially saturated pore units, the pressure difference is close to the capillary pressure, but this pressure is not a true average of the sample; it only represents specific pore units and it is virtually impossible to measure experimentally.

Figure 8

Investigation of the effect of averaging on the pressure difference between air and water. The dashed line is the quasi‐static capillary pressure‐saturation curve, while the black solid line is the centroid‐corrected pressure difference and the red solid line is the volume‐averaged pressure difference.

Thus, by averaging over all pore units, is not equal to , which results in a pseudo‐dynamic effect and a value for , following equation (1). Therefore, the dynamic effect that we determine by using volume‐averaged pressures can be partially contributed to the existence of saturated pore units inside the pore‐unit assembly; in the beginning of primary drainage, the contribution of saturated pore units to the dynamic effect is obvious, because saturated pore units form a distinct saturated domain having a layer of unsaturated pore units at the air reservoir. But, as drainage progresses, saturated pore units become smaller in number and also randomly dispersed inside the unsaturated pore‐unit assembly, yet they still contribute to a dynamic effect. In the following sections, we still use the volume‐averaged pressures.

Effect of Boundary Pressure

To quantify the effect of boundary pressure values , we simulated dynamic drainage for three different values of : 7, 10, and 14 kPa. In general, a higher value of resulted in a higher value of average pressure difference, for saturation values between 0.82 and 1.0 (see Figure 9a). However, for saturation values lower than 0.82, the average pressure difference curves converged to the quasi‐static capillary pressure‐saturation curve. During the saturation interval of 0.82–1.0, the rate of saturation change is relatively large (see Figure 9b). Figure 9c shows the value of as function of saturation, which quantifies the dynamic effect. For saturations of 0.82 or higher, was relatively high but it decreased during drainage to almost zero. The decrease in was a result of the average pressure difference that converged to the capillary pressure (see Figure 9a). The value of increased again after a saturation of 0.5. This was related to an increase in the deviation of the pressure difference of air and water and the capillary pressure (see Figure 9a) as well as the low value of , which converged to zero when the saturation reached the residual saturation. Note that the value of for low saturations is considered to be not reliable because converges to zero.

Figure 9

Investigation of dynamic effect for three dynamic drainage simulations using of 14 kPa (black line), 10 kPa (red line), and 7 kPa (blue line). Results are reported as (a) average pressure difference versus saturation, (b) rate of saturation change versus saturation, and (c) as a function of saturation.

Effect of Domain Size

We conducted dynamic simulations for 7kPa for five different domain sizes, with number of particles varying from 4,000 to 12,000. Results are shown in Figure 10. In general, we found that the larger the domain size, the larger the value of . The increase of is thus related to the increase in flow length ( ), which was experimentally confirmed by Abidoye and Das (2014). According to Stauffer (1978), Manthey et al. (2005, 2008), Dahle et al. (2005), and Bottero et al. (2011), is linearly proportional to . To test this proportionality in our simulations, we plotted values of for various domain sizes as a function of saturation in Figure 10b. We can see that the curves approximately coincide with each other. This indicates that in our simulations, .

Figure 10

(a) Dynamic coefficient as a function of saturation for different number of particles and (b) values of normalized to .

Illustration of water in a cube for a water saturation lower than , (a) illustration of the whole cube, (b) illustration of water in a corner, (c) smaller cube that contains the eight corners in Figure B1b, (d) illustration of an edge of the cube, and (e) cross‐sectional section of the edge shown in Figure B1d.

Effect of Permeability

To study the effect of permeability on , we varied the value of parameter in equation (24) from 1 to 10−4. Thus, we decreased the permeability of the packing from 2.11 × 10−10 to 2.11 × 10−14 m2, respectively, without changing the pore structure, thus the entry pressures are not affected. Detailed results are reported in the supporting information, section S3. Results showed that averaged pressure differences were not affected by changing the permeability. However, the permeability affects all fluxes through the pore throats ( , see equation (6)). This leads to a proportionally smaller rate of change of saturation with decreasing permeability, which in turn results in larger values of . Thus, we found that is inversely proportional to the permeability. This relation was also suggested by Stauffer (1978) and by Manthey et al. (2005).

In case that drainage is fast, the assumption of a constant air pressure may not hold; air pressure gradients may become significant in experiments while in simulations they are assumed to be constant. Of course, a constant air pressure is not always a valid assumption; air pressure gradient may arise if flow rate is high, if a rate limiting membrane is added to the air reservoir boundary (see Hou et al., 2012, 2014), and in the case of heterogeneous layered samples (Sakaki et al., 2011).

Simulations Versus Experiments

Dynamic effects in our simulations are small (less than 5.5 Pa s) for saturation values higher than 0.82 whereas for saturation values lower than 0.82, dynamic effects barely exist (i.e., is close to zero). However, experimental studies typically obtain larger values of , for example, Zhuang et al. (2017) reported a value of in the range of 5 × 104 to 3 × 106 Pa s. But in experiments the sample size is larger than that of our simulations. The domain size in our model is 3.2 × 3.2 × 3.2 mm3 whereas the sample size in experiments by Zhuang et al. (2017) was 30 × 30 × 20 mm3.

The small value of in our simulations agrees with other studies on pore‐scale simulations, which typically are conducted on small domains. For example, Joekar‐Niasar et al. (2010) found values between 100 and 300 Pa s, based on their dynamic pore‐network simulations for a viscosity ratio of 0.1 and a domain size of 6 × 6 × 6 mm³. Dahle et al. (2005) reported a value of 274 Pa s based on simulations of a bundle of tubes that had a characteristic length scale of 10−3 m. Similarly, Mumford and O'Caroll (2011) found values of between 30 and 230 Pa s using a bundle of tubes model for varying contact angels.

To compare the magnitude of values from our simulations with those of experiments, we scaled values ( ) by domain size and permeability, based on the findings of previous subsections: where and are the permeability values of simulations and experiments, respectively, and denotes the length of the domain in flow direction. The value of in simulations (less than 5.5 Pa s) was scaled to (less than 6.000 Pa s), which is closer to, but still lower than, the values in experimental data.

Of course, care must be taken when drawing conclusions on the comparison between experiments and simulations. As given by equation (1), for determining the pressure difference between the nonwetting and wetting phase pressure are needed. In this study, the function of averaged pressure difference and saturation is nonmonotonic for dynamic drainage simulations, which is a direct result of flow patterns and the averaging theorem as discussed in section 4.2.1; while air fingers progress toward the water reservoir, the pressure differences decreases toward the capillary pressure. The effect of breakthrough was also discussed in Joekar‐Niasar et al. (2010), where a decrease in pressure difference was obtained after breakthrough. In our simulations, dynamic effects disappear after air fingers reached the water reservoir. The saturation at which air reaches the water reservoir can be decreased by using a larger modeling domain, such that more fingers can exist and the effect of one finger is substantially smaller. However, modeling a larger domain is not possible with our code because of the small time steps that are required for dynamic drainage in an irregular pore‐unit assembly.

In literature, various observations have been made on the shape of the pressure difference‐saturation curves. Using pore‐scale models, Joekar‐Niasar et al. (2010, 2011) reported nonmonotonic curves, whereas Mumford and O'Carroll (2011) obtained monotonic curves. Similar observations are made in experimental studies, where monotonic curves are reported such as in Sakaki et al. (2011), Camps‐Roach et al. (2010) and in Zhuang et al. (2017), while Bottero et al. (2006) reported a nonmonotonic curve. To conclude, the magnitude of in pore‐scale simulations can be scaled from small modeling domains to larger experimental samples (using equation (33)). But in our case, upscaling would also entail increasing the number of fingers that can exist within a domain, thus upscaling will also affect the shape of the pressure difference‐saturation curve, and thus the functionality of with saturation.

In addition, when comparing pore‐scale simulations of an idealized pore structure with experiments on sand samples, the difference in pore structure may play a significant role in discrepancies values obtained from simulations and experiments. Sands are typically irregular in shape and will have surface roughness, which is of course simplified when using an assembly of regular shapes. The complex solid surface in real porous materials may lead to (i) interfaces getting trapped more easily, thus reducing the rate of saturation change and consequently larger values, (ii) reduction or enhancement of air‐finger formation, which can affect as a function of saturation, or (iii) a significant change in the local capillary pressure‐saturation curve that is idealized in simulations (such as by equation (20)).

Summary and Conclusions

This work presents a pore‐scale simulation technique for dynamic drainage in a packing of spheres, using the pore‐unit assembly approach. The pore space inside the packing of spheres was extracted by a sequence of meshing techniques: first, a triangulation was applied to extract grain‐based tetrahedra, which were then merged to find pore units. Subsequently, the pore space of each pore unit was replaced by a regular shape, for which the capillary pressure‐saturation curves were determined. Using the assembly of pore units, a pore‐scale version of IMPES (implicit pressure and explicit saturation scheme) was employed to solve for water pressure during dynamic drainage. The dynamic model was then verified by simulating slow dynamic drainage on a pore‐unit assembly. The resulting capillary pressure‐saturation curve was close to the simulated quasi‐static capillary pressure‐saturation curve of that same pore‐unit assembly.

Dynamic drainage simulations revealed that air infiltrated via preferential pathways in the pore structure, resulting in air fingers. During finger infiltration, dynamic effects occurred; the pressure difference between air and water was much larger than the capillary pressure. However, the air reached the water reservoir at relatively high values of saturation (0.82), after which further drainage occurred with the pressure difference between air and water being close to the global capillary pressure value. The enhanced values of pressure difference were mainly due to the employed averaging method, which gives much weight to saturated pores. When one averages over partially saturated pore units only, the air‐water pressure difference gets close to the capillary pressure. The value of obtained from simulations was found to be low (<5.5 Pa s), which is inherent to pore‐scale modeling of dynamic drainage in a small domain.

Nomenclature

parameter to adjust the overall permeability of the porous media.

area of water in the cross section of an edge in a regular shape, L 2.

air‐water interfacial area in a pore unit, L 2.

surface area of water in the cross section of pore throat , L 2.

cross‐section area pore throat , L 2.

water conductivity of pore throat , T L 4 M −1.

macroscale permeability, L 2.

length, L.

length of an edge of a regular shape, L.

distance between the centers of pore units and , L.

number of edges in a regular shape.

number of pore units connecting to pore unit .

number of saturated pore units adjacent to pore unit .

total number of pore units.

volumetric water flux through pore throat , L 3 T −1.

radius of curvature of the air‐water interface in a pore unit, L.

hydraulic radius of pore throat , L.

radius of the inscribed sphere of pore unit , L.

pore throat radius, L.

water saturation in pore unit .

air saturation in pore unit .

saturation corresponding to the inscribe sphere of a pore unit.

macroscale water saturation.

air pressure, M T −2 L −1.

water pressure in pore unit , M T −2 L −1.

capillary pressure in pore unit , M T −2 L −1.

entry pressure of pore throat , M T −2 L −1.

pressure difference between air and water reservoir, M T −2 L −1.

macroscale water pressure, volumetric averaged, M T −2 L −1.

macroscale capillary pressure, interfacial area averaged, M T −2 L −1.

volume of pore unit , L 3.

water volume in pore unit , L 3.

water volume in corners of a pore unit, L 3.

water volume in edges of a pore unit, L 3.

centroid of averaging domain, L.

surface tension of water, M T −2.

dihedral angle of an edge, rad.

geometrical parameter of a regular shape.

geometrical parameter of a regular shape.

water viscosity, M T −1 L −1.

maximum ratio of pore throat radius to pore‐unit radius.

geometrical parameter that relates with .

dynamic coefficient, M T −1 L −1.

time step, T.

void volume associated with pore throat , L 3.

solid surface area in pore throat , L 2.

Supporting Information

Supporting Information S1

Click here for additional data file.