Effects of Pore-Scale Geometry and Wettability on Two-Phase Relative Permeabilities within Elementary Cells

Abstract

We study the relative role of the complex pore space geometry and wettability of the solid matrix on the quantification of relative permeabilities of elementary cells of porous media. These constitute a key element upon which upscaling frameworks are typically grounded. In our study we focus on state immiscible two-phase flow taking place at the scale of elementary cells. Pressure-driven two-phase flow following simultaneous co-current injection of water and oil is numerically solved for a suite of regular and stochastically generated two-dimensional explicit elementary cells with fixed porosity and sharing main topological/morphological features. We show that the relative permeabilities of the randomly generated elementary cells are significantly influenced by the formation of preferential percolation paths, called principal pathways, giving rise to a strongly nonuniform distribution of fluid fluxes. These pathways are a result of the spatially variable resistance that the random pore structures exert on the fluid. The overall effect on relative permeabilities of the diverse organization of principal pathways, as driven by a given random realization at the scale of the elementary cell, is significantly larger than that of the wettability of the host rock. In contrast to what can be observed for the random cells analyzed, the relative permeabilities of regular cells display a clear trend with contact angle at the investigated scale.

1. Introduction

Multi-phase flows within natural and engineered porous structures are ubiquitous [1,2]. In the context of groundwater quality and subsurface energy management and exploitation, the migration of NAPLs (Non-Aqueous Phase Liquids) through the vadose zone and multi-phase dynamics in oil bearing formations and their feedback with overlying fresh groundwater bodies constitute an example of the interest in the issues associated with our understanding of such flow settings. Despite recent progress in studying these processes, knowledge of them is still far from being complete. Multi-phase flows are characterized by a multiscale nature, being distributed over a multiplicity of scales ranging from pore- to continuum-scales. This has suggested the development and application of upscaling frameworks to derive macroscopic models able to capture key features of processes occurring at the pore-scale. Upscaled models are typically formulated in terms of effective medium properties that are evaluated on the basis of micro-scale characteristics. The latter are commonly evaluated within elementary (or unit) cells (e.g., [3,4,5,6,7]).

Upscaling procedures based on homogenization via multiple scale expansion [4,8] or volume averaging techniques [9,10] have been proposed to analyze immiscible continuous or dispersed (e.g., [6]) two-phase flow in porous media as well as non reactive (e.g., [3,5]) and reactive (e.g., [11] and references therein) transport processes. The upscaling of immiscible two-phase flow has also been proposed through the thermodynamically constrained averaging theory by Gray and Miller [12]. Auriault [4] and Whitaker [9,10] derived effective models for two-phase immiscible fluid flow that include terms representing the effect of viscous coupling between flowing phases. When the latter phenomenon is negligible or the same pressure gradient is applied to both fluids, these models reduce to the traditional Darcy-Buckingham equation, i.e., $q_{\alpha }={K{K}_{r,\alpha }/{\mu }_{\alpha }\nabla p}_{\alpha }$, with subscript α indicating fluid phase, and $q_{\alpha }$, ${\mu }_{\alpha }$, $\nabla p_{\alpha }$, and K respectively being Darcy velocity, dynamic viscosity, pressure gradient across fluid phase α, and intrinsic permeability of the medium. Here, $K_{r,\alpha }$ indicates the relative permeability of fluid phase α, lumps all complexities of the interfacial fluid dynamics behavior and strongly depends on the pore-scale features. Experimental (e.g., [13,14]) and numerical (e.g., [15,16,17]) pore-scale studies focus on the way $K_{r,\alpha }$ is influenced by the viscosity ratio between the two fluids, the capillary number, medium wettability and initial configuration of the fluids in the system. Contact angles, which are associated with medium wettability, are recognized as key parameters in controlling macroscopic features of multiphase flow, such as fluid arrangement in the pore space and relative permeabilities. Nowadays there is a strong potential for their direct evaluation through advanced imaging techniques ([18,19] and reference therein). This element is also evidenced in Trojer et al. [20], wherein the displacement pattern of two migrating fluid phases is analyzed for diverse contact angles in a radial Hele-Shaw cell packed with glass beads. It has also been suggested that the micro-scale structure of the void space might play a critical controlling role on relative permeabilities and on the spatial arrangement of the flowing phases [21] at diverse spatial scales of interest. In this context, Liu et al. [22] presented a numerical study of an unsteady displacement problem for a two-dimensional setting, in which the solid phase is modeled through a set of non-overlapping circles associated with a random distribution of their diameter and location of their center, representing disaggregated porous media. Random realizations of pore-scale topological features have been employed in the context of pore network modeling of immiscible two-phase flows by Jiang et al. [23]. These authors investigated the effects on flow properties of random microscopic rock features, highlighting the key role of the connectivity function.

The choice of the geometrical features of the pore space of elementary cells significantly impacts effective coefficients evaluated with upscaling procedures (e.g., [24,25,26]). Icardi et al. [27] pointed out the need to include uncertainty in pore-scale physics to be used in upscaling models.

The assessment of the joint role that the uncertain geometry of an elementary cell and medium wettability have on relative permeabilities is still an open and challenging problem. Here we tackle this issue and study the role of wettability and pore-scale structure of realistic randomly generated and regular elementary cells in fluid phase distribution and relative permeabilities characterizing steady state immiscible two-phase flows. We do so through a set of numerical simulations of two-phase fluid flow performed in a diffusive-interface framework within two-dimensional explicit pore spaces representing elementary cells which can be employed for upscaling studies.

2. Materials and Methods

2.1. Generation of Regular and Random Elementary Cells

We consider water and oil as immiscible displacing fluids, i.e., the fluid phase α = w (water) or o (oil). We explore a suite of twelve two-dimensional elementary cells: (i) three media associated with an explicit regular structure of the pore space (hereafter termed as settings HDη, with η = a, b, c); and (ii) nine stochastically generated explicit pore spaces (hereafter termed Ei, with i = 1, …, 9). Each case is characterized by the same macro-scale characteristics, i.e., size of the medium along x-, L_x = 3.98 × 10⁻³ m, and y-, L_y = 3.45 × 10⁻³ m directions (see Figure 1), and connected porosity, n ≈ 0.5.

Systems in which solid elements are regularly distributed in space, as in HDη, are commonly employed in upscaling studies of flow and transport processes [6,7]. The elementary cell HDa (see Figure 1a) is composed of circular grains of diameter D_g = 0.48 mm with throat diameter D_t = 0.11 mm. D_g is chosen to guarantee that the Specific Surface Area (SSA) of the regular geometries is equal to the average SSA of the nine random samples, while D_t controls the sample porosity. Systems HDb and HDc are characterized by elliptical and square solid elements with round corners, respectively (see Figure 1b,c).

The random elementary cells we consider are generated on a rectangular lattice of size L_x/2 and L_y/2 using the procedure implemented in Hyman et al. [28] and in Siena et al. [29], to which we refer for additional details about the generation techniques and code as well as the associated generation parameters. The synthetic generation approach employed has been shown to produce realistic pore spaces and pore space geometric observables such as porosity and SSA to be controlled through the appropriate tuning of generation parameters [30]. Virtual pore spaces generated using this method were used to investigate the influence of (a) porosity on transport properties [28]; (b) porosity and mean hydraulic radius on single-phase fluid flow permeability [31]; (c) pore wall geometry and network topology on local mixing [32]; and (d) pore size distributions on fluid velocity statistics [29], as well as to assess the way in which statistics of porosity and SSA are affected by scale [33,34].

The generated samples are then mirrored along the x- and y-directions to obtain the periodic heterogeneous pore structures required by upscaling approaches. Figure 1d–l depicts the nine heterogeneous elementary cells we consider. As listed in

Table 1. Connected porosity (n), Specific Surface Area (SSA), intrinsic permeability (K) and maximum likelihood (ML) estimates of range and sill (${\widehat{r}}_i$ and ${\widehat{s}}_i$; i = x, y) of directional variograms of the indicator function of the pore space.

Sample Identifier	n	SSA (10−3 m−1)	K (10−10 m2)	${\widehat{\mathit{r}}}_{\mathit{x}}$	${\widehat{\mathit{s}}}_{\mathit{x}}$	${\widehat{\mathit{r}}}_{\mathit{y}}$	${\widehat{\mathit{s}}}_{\mathit{y}}$
HDa	0.48	8.63	4.97	-	-	-	-
HDb	0.50	8.63	11.67	-	-	-	-
HDc	0.52	8.63	10.47	-	-	-	-
Average HDa–HDc	0.50	8.63	9.06	-	-	-	-
E1	0.48	8.39	1.25	0.12	0.23	0.12	0.21
E2	0.49	7.84	2.42	0.16	0.25	0.17	0.23
E3	0.48	8.61	3.43	0.17	0.25	0.12	0.22
E4	0.51	9.11	3.84	0.11	0.21	0.13	0.23
E5	0.49	9.24	4.43	0.18	0.26	0.14	0.24
E6	0.50	8.20	2.29	0.16	0.25	0.21	0.24
E7	0.49	7.69	1.01	0.12	0.22	0.25	0.27
E8	0.49	9.44	3.11	0.13	0.24	0.16	0.27
E9	0.49	9.13	1.04	0.14	0.27	0.13	0.24
Average E1–E9	0.49	8.63	2.54	0.14	0.24	0.16	0.24

, the connected porosity of all analyzed samples is virtually the same and equal to the value set for the regular media (≈0.5). Table 1 also lists values of sample SSA. This quantity contributes to control transport dynamics and does vary considerably in all of the cases.

2.2. Pore Space Spatial Correlation

All random cells are characterized by the same structure of the pore space spatial correlation. The latter has been described by analyzing the directional variograms ${\gamma }_i(h_i)$ ($h_i$ being spatial separation distance rescaled by the domain size L_i with i = x, y) of the indicator function, I (I = 1 in the solid matrix and I = 0 otherwise), of the pore space. These variograms are interpreted through a Maximum Likelihood (ML) approach (details not shown) via a spherical model, i.e., ${\gamma }_i(h_i)=s_i\left(3h_i/2r_i-h_i^3/2r_i^3\right)$ for 0 ≤ $h_i$ ≤ $r_i$ and ${\gamma }_i(h_i)=s_i$ for $h_i$ > $r_i$, $s_i$, and $r_i$ respectively representing variogram sill and range (rescaled by L_i). ML estimates of $s_i$, ${\widehat{s}}_i$, and $r_i$, ${\widehat{r}}_i$, for all random samples are listed in Table 1. The degree of correlation of the void spaces, as quantified by ${\widehat{r}}_i$, does not vary significantly with direction (indicating isotropy of the correlation structure of the pore spaces) and across samples (i.e., the spatial correlation structure of void spaces and solid matrix is statistically equivalent for all geometries). The variogram sill is fairly constant for all geometries and close to the theoretical value $n\left(1-n\right)=0.25$.

2.3. Topological/Morphological Analysis of the Elementary Cells

Both regular and random elementary cells are associated with comparable topological and morphological features. These characteristics are not controlled by the generation algorithm we employ and can be assessed a posteriori. We investigate the structure of the pore spaces within a pore-network framework using the maximal ball algorithm described in Dong and Blunt [35]. Each pore is associated with a coordination number, Cn, representing the number of pores connected to the pore under consideration. Figure 2 depicts the relative frequency of Cn evaluated in the nine random samples. The patterns of the shapes of the distributions associated with different geometries indicate that the systems share similar topological characteristics.

Table 2. Topological and morphological information of the elementary cells: total number of pores, N_pores; average coordination number, $\langle Cn\rangle$ ; ratio $\langle Cn\rangle /N_{pores}$; area of the solid phase, M₀; boundary length between solid and fluid, M₁; and Euler number, M₂.

Sample Identifier	Npores	$\mathbf{\langle }\mathit{C}\mathit{n}\mathbf{\rangle }$	${\mathbf{\langle }\mathit{C}\mathit{n}\mathbf{\rangle }\mathbf{/}\mathit{N}}_{\mathit{p}\mathit{o}\mathit{r}\mathit{e}\mathit{s}}$ (%)	M0 (10−6 m2)	M1 (10−3 m)	M2
HDa	112	3.00	2.68	8.10	69.91	−173.25
HDb	112	3.00	2.68	8.53	73.19	−181.00
HDc	112	3.00	2.68	10.73	80.77	−64.25
Average HDa–HDc	112	3.00	2.68	9.12	74.62	−139.50
E1	64	2.63	4.10	7.20	60.45	−145.25
E2	72	2.72	3.78	7.03	55.10	−133.25
E3	80	2.70	3.38	7.16	61.61	−147.50
E4	76	2.37	3.12	6.74	61.39	−160.75
E5	100	2.68	2.68	6.95	64.24	−151.00
E6	80	2.55	3.19	6.87	56.28	−134.25
E7	84	2.48	2.95	7.06	54.30	−135.75
E8	100	2.32	2.32	9.44	65.62	−162.50
E9	104	2.35	2.26	9.13	63.93	−157.75
Average E1–E9	84.4	2.53	3.09	6.99	60.32	−147.56

lists the total number of pores, N_pores, the average coordination number, $\langle Cn\rangle$, and the ratio $\langle Cn\rangle /N_{pores}$ computed for the regular structures (HDa, HDb, HDc) and for each of the random samples considered. The ratio $\langle Cn\rangle /N_{pores}$, which is a measure of the global connectivity of the pore network, is constant for the three regular domains HDη and varies weakly among the nine stochastic pore spaces, the coefficient of variation of $\langle Cn\rangle /N_{pores}$ being about 20%. Note that, on average and in seven of nine cases, the heterogeneous structures appear to be more connected than the regular ones.

To further characterize the morphology and the topology of the pore structures, we evaluate the Minkowski functionals by relying on the algorithm proposed by Mantz et al. [36]. These functionals are geometric measures defined for binary structures (see for details [37,38,39]). In a two-dimensional setting they correspond to the area of the solid phase (M₀), the boundary length between solid and fluid (M₁) and the Euler number (M₂). While M₀ and M₁ are morphological measures, M₂ is a topological metric describing connectivity. The latter is calculated by counting the number of isolated objects from which the number of redundant connections is subtracted and is negative for well connected structures and positive otherwise. The results listed in Table 2 show that all of the considered structures are well connected. Note also that, according to M₂, the regular structure HDc appears to be characterized by the lowest connectivity.

2.4. Single Phase Flow Simulations

All flow simulations are performed using the finite element software COMSOL Multiphysics [40] by relying on a triangular element mesh. A body force of magnitude ΔP/L_x, with ΔP = 39.8 Pa, acts along the x-direction (from left to right in Figure 1). This magnitude has been selected to be comparable to a body force produced by the gravity (ρ|g|, where ρ = 1000 kg m⁻³ is the fluid density and |g| is the absolute value of gravitational acceleration).

We start by solving the steady state single phase fluid problem to evaluate the intrinsic permeability $K=Q{\mu }_{\alpha }{L}_x/\left(L_y\Delta P\right)$, where the flow rate, Q, is obtained by integration of the steady state fluid velocity component along the y-axis at the system outlet (i.e., at x = L_x). The values of K associated with each scenario are listed in Table 1. These vary significantly across all elementary cells (either regular or random), the coefficient of variation of K for the random media being about 50%. The degree of heterogeneity of each cell is further quantified in terms of tortuosity values. We do so by considering the trajectories of N tracer particles which are released uniformly across the domain inlet (at x = 0) and advective displacement of which is tracked through the fluid domain. The tortuosity of a particle trajectory i is computed as ${\tau }_i=l_i/L_x$, l_i being the trajectory length. Particles reaching stagnation regions are removed from the analysis, which then considers N_P particles percolating through the medium. Relying on N = 5000 particles enabled us to obtain stable basic statistics of $\tau$ (details not shown). Depending on test case, 0.46 ≤ N_P/N ≤ 0.98, as listed in

Table 3. Number of percolating particles, N_P, and basic statistics of tortuosity, $\tau$.

Sample Identifier	Np	Mean	Variance	Minimum	Maximum	Skewness	Kurtosis
HDa	4906	1.18	1.5 × 10−4	1.06	1.21	0.44	1.00
HDb	4947	1.14	6.5 × 10−5	0.96	1.16	−9.33	197.76
HDc	4935	1.21	3.3 × 10−4	1.13	1.27	1.27	3.87
Average HDa–HDc	4929	1.18	1.8 × 10−4	1.05	1.21	−2.54	67.54
E1	4726	1.30	5.3 × 10−3	1.21	1.51	0.56	2.64
E2	2305	1.32	2.0 × 10−2	1.15	1.92	0.71	3.46
E3	2422	1.34	5.4 × 10−2	1.12	2.14	1.23	3.43
E4	3729	1.17	3.7 × 10−3	1.09	1.42	0.63	2.81
E5	3471	1.42	2.7 × 10−2	1.18	1.74	0.54	1.97
E6	3730	1.80	2.6 × 10−2	1.15	2.11	−1.32	7.85
E7	3151	1.26	8.7 × 10−3	1.15	1.66	1.21	3.67
E8	4480	1.21	1.2 × 10−2	1.06	1.42	0.06	1.59
E9	4697	1.17	3.7 × 10−3	1.09	1.42	0.64	2.76
Average E1–E9	3635	1.33	8.1 × 10−3	1.10	1.50	0.64	2.67

. Figure 3 depicts the histogram of tortuosity, $\tau$, for HDa and three selected random elementary cells. The mean, variance, skewness, and kurtosis of $\tau$, are listed in Table 3. The tortuosity of HDη (η = a, b, c) varies over a narrow range of values and its histogram is characterized by a single peak. Otherwise, histograms of $\tau$ obtained in the randomly heterogeneous geometries exhibit a multimodal behavior and are characterized by an appreciable range of variability within and across samples (see Table 3 and Figure 3). This behavior is related to the complex spatial distribution of the void spaces characterizing the heterogeneous samples.

Colorized images of particle trajectories are depicted in Figure 3. Particles migrating in the randomly heterogeneous media tend to focus within areas of high velocity, forming trajectory bundles termed principal pathways [28]. The multimodal distribution of tortuosity for the heterogeneous cases can be related to the grouping of particle trajectories along principal pathways during percolation through the media. The emergence of these features is an intrinsic characteristic of the flow field taking place in such intricate pore structures, as opposed to what can be observed in regular pore spaces, where flow is spatially distributed according to a mostly uniform pattern. Note the occurrence of pore space regions that are not visited by particles during the time frame of the simulation (white regions in Figure 3). These observed behaviors will affect the results of our two-phase flow simulations, the ensuing fluid distribution in the pore space and the magnitude of relative permeabilities as a function of fluid saturation as discussed in Section 3.

2.5. Two-Phase Flow Simulations

Two-phase flow within the elementary cells is solved through a diffusive-interface formulation which simulates the evolution of the fluid/fluid interface within an energy based framework [41,42]. Within this approach the interface between the two phases is approximated through a small but finite transition region inside which the two components are mixed. The thickness of this region is controlled by a capillary length scale parameter, ε [41,42]. This methodology has been applied to describe steady state two-phase flows and displacement processes in porous media by Ahmadlouydarab et al. [17] and Amiri and Hamouda [43,44]. It requires the joint solution of the Cahn-Hilliard (1) and Navier-Stokes (2) equations

$$\frac{\partial \varphi }{\partial t}+v\cdot \nabla \varphi =m{\nabla }^{2}G$$

(1)

$$\rho \frac{\partial v}{\partial t}+\rho v\cdot \nabla v=-\nabla p+\nabla \cdot \left[\mu \left(\nabla v+\nabla v^{\mathrm{T}}\right)\right]+G\nabla \varphi +{\Delta P/L}_x;\nabla \cdot v=0$$

(2)

Here φ is the phase variable (φ = 1 in water and φ = −1 in oil), and v, p, m, and G respectively indicate fluid velocity vector, pressure, the Chan-Hilliard mobility constant and the chemical potential of the system. The latter arises from the variation of the system free energy with respect to φ and is defined as $G=\lambda \left[-{\nabla }^2\varphi +\varphi {({\varphi }^2-1)/\epsilon }^2\right]$, λ being the mixing free energy density that arises from the molecular interaction between the two phases inside the transition region [41].

At equilibrium for a planar fluid/fluid interface the ratio λ/ε is proportional to the interfacial tension, $\sigma$, as $\sigma =2\sqrt{2}\lambda /\left(3\epsilon \right)$ [42,43], the quantity $G\nabla \varphi$ in (2) being the diffusive-interface representation of $\sigma$. Fluid density, ρ, and viscosity, μ, are defined as the average between the corresponding values associated with each fluid phase weighted by their volume fraction as $\xi ={\xi }_w\left(1+\varphi \right)/2+{\xi }_o\left(1-\varphi \right)/2$ with $\xi$ = ρ, μ. The fourth-order Chan-Hilliard Equation (1) is decomposed in two second-order equations through the auxiliary function ψ, i.e.,

$$\frac{\partial \varphi }{\partial t}+v\cdot \nabla \varphi =\frac{m\lambda }{{\epsilon }^2}{\nabla }^2\psi ;\psi =-{\epsilon }^2{\nabla }^2\varphi +\varphi \left({\varphi }^2-1\right)$$

(3)

Computations are therefore performed using the set of four independent variables {v, p, φ, ψ}. Bi-periodic boundary conditions are applied to the external sides of the domain for all variables of the system. The following boundary conditions are employed at the solid-fluid boundaries

$$v=0;n\cdot \nabla \psi =0;n\cdot \frac{\nabla \varphi }{\left|\nabla \varphi \right|}=\cos \theta$$

(4)

with n being the unit normal vector pointing from the fluid to the solid. For the purpose of the present study we assume no slip condition on the solid matrix, ensured by the first condition of (4). The contact line motion is then due solely to Chan-Hilliard diffusion. Zero flux through the solid matrix is imposed via the second of (4), while the third of (4) corresponds to the diffuse-interface formulation of the contact angle, θ.

In our settings the oil phase is considered to reside initially in the central part of the system according to a given saturation, S_o, as exemplified in Figure 4 for HDa and E4 with S_o = 0.5. The initial configuration of the two phases has been selected considering that one of the assumptions underlying upscaling of two-phase flow is that the volume occupied by each fluid is connected [4]. One or both phases start flowing in the system under the action of the constant driving force ${\Delta P/L}_x$ until steady state conditions are reached and the resulting equilibrium configuration of the water-oil interface is observed. During the simulation the saturation of the two fluids is controlled by the bi-periodic boundary conditions applied to the external sides of the domain.

In typical experimental conditions the flow rate of the two fluid phases is imposed at the domain inlet and the pressure gradient, $\Delta P$, between inlet and outlet is measured (see, e.g., Tallakstad et al. [14]). One can note that under these conditions the steady-state regime is commonly not identified in terms of a static equilibrium configuration of the water-oil interface and is considered to be attained in an average sense when the values of $\Delta P$ appear to fluctuate around a constant value.

Benchmarking of the numerical approach has been performed through comparison against the two-phase Poiseuille flow (details not shown). A key issue is the proper selection of the Chan-Hilliard model parameters, i.e., the capillary width, ε, and the mobility constant m [45,46] as well as the grid size h. The capillary width is associated with the Cahn number, $Cn=\epsilon /D_t$, for the definition of which we select the pore throat of HDa as the characteristic length scale. The transition region (and therefore the associated Cn) needs to be very small to capture the physics of the two-phase flow problem [41,42]. The convergence of the numerical results with variations of h, ε and m has been verified for all settings (details not shown). On these bases and according to Amiri and Hamouda [42], we set $Cn=0.1$ a maximum value of h equal to 2 ε and m = ε².

The simulations are performed by considering ${\mu }_w$ = 0.001 Pa s and ${\mu }_o$ = 0.002 Pa s (i.e., a viscosity ratio ${\mu }_o/{\mu }_w$ = 2), water and oil density equal to 1000 kg m⁻³, and interfacial tension $\sigma$ = 0.04 N m⁻¹, these values being representative of a typical oil/water system at reservoir condition [43,47]. Medium wettability is defined through the contact angle, θ, measured with respect to the water fluid phase. We set θ = 90°, 60° and 120° to mimic neutral wet (NW), water wet (WW) and oil wet (OW) systems, respectively. The distribution of pore-scale velocities is then computed for diverse values of S_o.

The magnitude of the imposed body force (the same used in single-phase simulations, Section 2.4) controls the flow rate of the two phases at steady-state and the capillary number $C{a}_{\alpha }$ = ${\mu }_{\alpha }{v}_{\alpha }^M/\sigma$ ($v_{\alpha }^M$ being the largest velocity of phase α along the x-direction).

We note that the largest value of the capillary number, is observed for the oil phase, and its magnitude is of the order of 10⁻⁴. This number represents the ratio of viscous to capillary forces, and the value we employ can be considered representative of standard oil/water systems at reservoir condition [44]. Albeit of interest, the study of the influence of diverse values of the capillary number (e.g., [14]) on relative permeability curves is outside the scope of the present contribution.

3. Results and Discussion

3.1. Equilibrium Two-Phase Distribution Patterns

Figure 5 depicts patterns of equilibrium (steady state) spatial distributions of the two fluid phases and streamlines in HDa and in three selected heterogeneous elementary cells for θ = 90° and S_o ≈ 0.5 (the actual values of S_o at steady state are listed in the caption of Figure 5). Pore space geometry significantly affects the steady state distribution of the interface between the two fluids and the flow field characteristics. The pattern of the flow fields in the random elementary cells is characterized by several regions not contributing to the total flow rate (clearly identifiable in Figure 5). For this degree of oil saturation both phases percolate in HDa (see Figure 5a) and in the heterogeneous samples E4 (Figure 5c) and E8 (Figure 5d), the water phase being blocked in sample E1 (Figure 5b). This behavior is related to the complex topology of the investigated geometries within which the two fluids are displaced and experience spatially variable resistances that are differently distributed depending on the realization analyzed. The percolation of the two fluids is strongly linked to (a) the number of inlet/outlet openings from which a fluid can be injected into (or exit from) each sample and (b) the distribution of principal pathways within the samples. The two fluids may flow simultaneously in cells characterized by more than one inlet and multiple principal pathways favoring the occurrence of several distinct percolation paths in the system (see Figure 3c,d). One or both fluids are immobile (at steady state) when the medium has only one inlet and/or one principal pathway (see Figure 3b). The emergence of the general patterns observed in Figure 5 has been found for all wettability values analyzed (see Figure 6).

3.2. Relative Permeabilities

The flow rate of the α phase, Q^α, for a given oil saturation is computed by integrating the steady state fluid velocities along the portion of the outlet boundary through which the α phase flows. Normalization of Q^α by the total flow rate computed when the medium is fully saturated by the α phase leads to an estimate of the relative permeability, K_rα. We further note that the solution of closure equations associated with effective models [4,8,9,10] and leading to the evaluation of a full relative permeability tensor which includes the coupling phenomenon, is equivalent to the solution of (1)–(4), where flow is sustained by an imposed body force acting on both fluids at steady state conditions.

Relative permeabilities, $K_{r,\alpha }^X$, versus oil saturation are depicted in Figure 7 in regular media and in three selected heterogeneous elementary cells under diverse wettability conditions (here the superscript X = NW, WW and OW). It can be observed that in the regular media $K_{r,\alpha }^X$ varies almost linearly with S_o for all wettability conditions analyzed (see Figure 7a–c). The simultaneous flow of both phases is observed for a wide range of saturations. Otherwise, the relative permeability values for a given saturation level consistently display a clear trend with contact angle. In this case the oil permeability under water wet conditions is higher than that corresponding to neutral or oil wet settings, the opposite behavior being documented for relative permeabilities of water. The relative permeabilities in the randomly heterogeneous cells are strongly influenced by the pore structure configuration (see Figure 7d–f) and vary in a highly non-linear fashion with fluid saturation. For instance, water relative permeability, $K_{r,w}^X$, remains constant for all wettability values when 0.4 ≤ S_o ≤ 0.6, as in the case depicted in Figure 7e. This behavior is markedly controlled by the number of principal pathways and by the distribution of the two phases in the heterogeneous domains. The resulting fluid configuration could have a significant feedback with the occurrence of regions (Figure 5c) which do not contribute to water mass transport and could be linked to the number of inlet/outlet openings from which a fluid can be injected into (or exit from) a given sample. To further explore this point, we analyze the sample probability distributions (PDFs) of the horizontal component of the velocity field, u_α (α = o, w, respectively for oil and water phase) normalized by the mean velocity of the corresponding phase. Results of these investigations are summarized in the Appendix A. It can be seen from these results that sample PDFs are not affected by saturation values for the homogeneous domains analyzed. For heterogeneous geometries the peak and the positive tail of the sample PDF of oil (or water) velocities tend to be respectively higher and longer with increasing saturation of the corresponding phase (i.e., S_o or S_w). This behavior is related to the observation that heterogeneous domains are characterized by the occurrence of regions with very limited contribution to oil/water mass transport.

Contrary to what is observed in the regular systems, the wettability conditions of a given random medium do not affect systematically the relative permeabilities of soil and water, the variability of which is mostly driven by the underlying random pore space configuration. This behavior is associated with the topological complexity of the heterogeneous cells, which shadows the influence of the contact angle on the final distribution pattern of the two fluid phases.

The role of the random pore-size distribution patterns on $K_{r,\alpha }^{NW}$ is further emphasized in Figure 8a where the $K_{r,\alpha }^{NW}$ curves evaluated for the regular settings and for the random samples are juxtaposed. The corresponding results for water and oil wet conditions are depicted in Figure 8b,c respectively. Comparison of Figure 7 and Figure 8 allows the dominant effect of random pore space geometry of the elementary cells to be recognized, as opposed to the impact of wettability on relative permeabilities. This is clearly evidenced by considering that the relative permeability values in Figure 8 are characterized by a significantly higher spread than that is evidenced in Figure 7 for a given saturation level. On the other hand, it is noted that the relative permeabilities associated with the three diverse regular settings analyzed can be barely distinguishable for a given wettability value.

Water relative permeability is generally smaller in the heterogeneous cells we analyze than in the HDη setting. This behavior is a consequence of the occurrence of regions where virtually no flow is observed for the water phase in Ei. The location of principal pathways within the medium in relation to the spatial pattern of fluid injection along the inlet section has a remarkable controlling role in the relative magnitude of oil and water relative permeabilities observed for the heterogeneous settings Ei as opposed to the regular pore distribution, where no distinct principal pathways are observed. For example, the relative permeability of oil for heterogeneous systems is always larger than that associated with the regular pore space patterns in cases where the segment along which oil injection takes place corresponds to the location of the principal pathways (e.g., samples E1, and E4; Figure 3b,c), the opposite behavior emerging when oil is injected through portions of the inlet that are not associated with principal pathways (e.g., sample E8; Figure 3d). We note that our work does not address the question of the existence or of the relevance of a representative scale for the assessment of the relative permeability curves of natural porous materials. Our findings are framed in the context of upscaling studies based on, e.g., volume averaging [48] or homogenization [49], wherein one relies on the solution of pore-scale physics at a scale that is much smaller than that of an investigated porous system.

In such a framework, results obtained at the micro-scale are used to evaluate macroscopic parameters (in our case relative permeabilities) which are intended to characterize a continuum-scale description of the system. When conditions of the applicability of these approaches are warranted, our results can then potentially be transferred to characterize larger scale systems.

The geometry of the pore space at the scale of the elementary cell is typically uncertain and this aspect impacts the quantification of macroscopic quantities such as relative permeabilities. Our results therefore support the need for considering such uncertainties when dealing with characterization of relative permeabilities based on upscaling frameworks.

3.3. Qualitative Comparison of the Numerical Results against Experimental Data

Finally, we provide a preliminary assessment of the ability of our generated pore structures to provide relative permeability curves which can be considered compatible/consistent with those arising from experimental observations. We do so by performing a qualitative comparison of our simulation results against experimental data collected within two-dimensional micromodels of porous media [13,50] under steady state conditions.

These authors analyze steady state two-phase flows at diverse capillary numbers and using diverse fluid pairs by way of two-dimensional micromodels forming a connected system of pores and throats. We collected a set of data from this study, which appears to be amenable to a qualitative comparison against our numerical results obtained under water wet conditions. The physicochemical properties of the fluids used in the experiments and in our computational analyses are listed in

Table 4. Physicochemical properties of fluids employed by Avraam and Payatakes [13], and Chang et al. [50] and used in this study.

Property	Experiments from Avraam and Payatakes [13]	Experiments from Chang et al. [50]	Our Computations
Non wetting fluid	n-Dodecane	Diesel Fuel/Tetrachloroethane	Oil
Wetting fluid	Water	Water	Water
${\mu }_o$ (Pa s)	1.36 × 10−3	not specified	2.0 × 10−3
${\mu }_w$ (Pa s)	0.94 × 10−3	not specified	1.0 × 10−3
$M={\mu }_o/{\mu }_w$	1.45	not specified	2.0
${\rho }_o$ (kg/m3)	730	820/1600	1000
${\rho }_w$ (kg/m3)	995	997	1000
$\sigma$ (mN/m)	25	48.7/38.75	40.0
$\theta$ (°)	40	0	60
Capillary Number	1.19 × 10−5	not specified	10−4

µ_α: dynamic viscosity of phase α (α = o or w for oil and water phase, respectively), ρ_α: density of phase α (α = o or w), σ: interfacial tension, θ: contact angle.

.

Figure 9a,b juxtapose the selected numerical results depicted in Figure 8b (curves) and the relative permeability values estimated in the considered experimental studies (symbols). Note that the experimental and numerical analyses have been performed under different conditions, e.g., in terms of the physicochemical properties of the fluids employed (see Table 4). It can be seen that the results associated with our random elementary cells are qualitatively consistent with the behavior displayed by the available relative permeability information. In this sense, it can be recognized that our relative permeability results for the random geometries display features that can be observed also in the experimental studies, namely: (a) the lack of connectivity of the oil phase for values of S_o smaller than a given threshold; and (b) relative permeability curves that vary in a highly non-linear fashion with fluid saturation.

When viewed also in this context, the results of our analysis are promising and suggest the worth of further systematic studies on the assessment of the way the uncertain geometry of elementary cells propagates to the macroscale description of porous media constructed through such systems.

4. Conclusions

We illustrate the role of complex micro-scale pore geometry and wettability on oil and water relative permeabilities within elementary cells under immiscible two-phase flow conditions. Two types of micro-scale elementary cells, characterized by comparable macro-scale properties and topological and morphological features, have been considered: (a) regular samples, constructed starting from a solid element of fixed geometry which is then distributed in the domain according to a regular pattern, and (b) random geometries, which share similarities with the complex pore structure of natural rocks. In the regular elementary cells, the two phases are well connected for several oil/water saturation values. This leads to relative permeabilities which vary (approximately) linearly with phase saturation and exhibit a clear trend with contact angle. In the random elementary cells, the relative permeability curves vary in a highly non-linear fashion with fluid saturation due to the occurrence of stagnant flow regions.

The key driver for the variability of relative permeability curves is the random geometry of the pore space of the elementary cell. This is clearly evidenced by our results (see Figure 7 and Figure 8), indicating that the relative permeability values associated with random pore spaces (for a given wettability) are characterized by a markedly higher spread than their counterparts related to diverse wettabilities (for a given random pore space). As our knowledge of the details of pore space structure and wettability is always incomplete, our findings: (a) have an impact on the assessment of the uncertainty associated with the quantification of relative permeabilities of porous materials stemming from the use of upscaling techniques which are based on the solution of closure equations formulated on elementary cells; and (b) suggest the need to perform systematic upscaling studies in a stochastic context to propagate the effects of uncertain pore space geometries to a probabilistic description of relative permeability at the continuum scale.

All data used in the paper will be retained by the authors for at least five years after publication and will be available to the readers upon request.