Homogenization of two fluid flow in porous media

1. Introduction

The macroscopic flow of multiple fluid phases in porous media, for example soil, is often described by Richards' equation [1–3]. This equation describes the local saturation under the influence of saturation and pressure gradients and is parametrized by the water release curve and the saturation-dependent hydraulic conductivity, both are measured experimentally [4]. Richards' equation offers challenges in terms of both parametrization [4–6] and the numerical solution [7,8].

Mathematically it has been shown, using the method of homogenization [9,10], that single phase flow in porous media can be approximated by Darcy's Law [11]. This equation can be derived from Stokes equations for single phase flow in the pore space and is parametrized by the hydraulic conductivity [3,11,12]. Such techniques have been applied in single porosity materials [3,11,13,14], dual porosity materials [15–17] and vuggy porous structures [18–22]. However, the homogenization process in partially saturated porous media is less well defined and relies on assumed knowledge of the interface location (ch. 5 in [3]). Knowledge of the air–water interface is often hard to obtain. Studies using X-ray computed tomography have been carried out [23]. However, these are computationally intensive and require scans to be carried out across the whole range of saturation. To overcome the need for repeated scanning, various researchers have suggested different empirical or approximate formula to describe the water release curve [2,3,24–26]. However, the water release curve exhibits multi-branched hysteresis loops [6] and needs to be parametrized with experimental measurements that can take months to gather even for a single branch of the hysteresis curve [5].

In order to derive the water release curve and saturation-dependent hydraulic conductivity, the dynamic interaction of the two fluid phases must be considered. One way to capture the physics of the two fluids is to use the Cahn–Hilliard phase field model coupled with Stokes equations [27–29]. The Cahn–Hilliard model can be derived through a free energy minimization [30,31] and has widely been applied to study moving contact lines and bubbles in a two fluid system [32]. This model overcomes the difficulties of a sharp interface by using the assumption that the interface between the two phases is of finite thickness which is assumed to be small in comparison with the geometry considered. By considering the limit when the interface thickness approaches zero the phase field model can be shown to reduce to standard free–boundary problems [32–34]. Homogenization of the Cahn–Hilliard model in porous media has previously only been studied for the case in which the interface thickness is comparable with the characteristic pore size [35,36]. This results in an effective Cahn–Hilliard equation where the interface mobility is derived through a set of cell problems.

In this paper, we consider the case of two immiscible fluids separated by an interface of finite width. This width is assumed small relative to the pore size and, hence, may be considered the smallest length scale in the problem. We use the method of homogenization to derive a coupled set of equations which describe the macroscopic flow properties of these fluids in partially saturated porous media. These equations, which are based on fundamental physical assumptions, are shown to reduce to Richards' equation in an appropriate parameter regime, i.e. the gas pressure is assumed constant. We assume that, to leading order, the interface positions and, hence, the water release curve are determined by capillary forces. The hydraulic conductivity is then determined by studying the perturbation due to a weak pressure gradient. The result is that the water release curve and the saturation-dependent hydraulic conductivity are determined through a series of cell problems which capture the porous geometry and the effects of the pressure gradients.

To capture the physics associated with two phase flow in porous media, the correct behaviour of the contact angle between the solid and the two fluid phases must be included in the Cahn–Hilliard formulation. The contact angle is one of the many factors associated with the hysteresis observed in the water release curve [37]. There have been numerous studies on the effect of the contact angle in the Cahn–Hilliard formulation [38–40]. Here, we make use of the more recent boundary conditions which are derived geometrically [40].

This paper is arranged as follows: in §2, we derive the Cahn–Hilliard fluid model using a free energy minimization and apply the method of homogenization to derive the two fluid model. The cell problems, which result from the calculations, allow the hydraulic conductivity and water release curves to be calculated based entirely on the underlying geometry. In §3, we study the solution properties of the cell problems and consider the limit when the interface thickness approaches zero. In §4, we solve the cell problems for the example case of air and water flowing in soil. Finally, in §5, we discuss our results and draw conclusions.

2. Derivation of homogenized equations

(a) Deriving the two fluid Cahn–Hilliard equation

We consider the interaction of two single component fluids, for example air and water, in a porous geometry as illustrated in figure 1. Throughout this derivation we use $\tilde{\cdot }$ to denote that ⋅ is a dimensional quantity. A list of the dimensional quantities used in this derivation is given in table 1. We consider a porous domain Ω comprising a solid matrix $\mathfrak{S}$ and a fluid part ${\mathfrak{B}}^{ϵ}$ with boundary $\mathrm{\partial }{\mathfrak{B}}^{ϵ}$. We assume that ${\mathfrak{B}}^{ϵ}$ is connected and that $\mathrm{\partial }{\mathfrak{B}}^{ϵ}$ is smooth. The geometry of the porous structure is of typical size ${\tilde{L}}_x$ and comprises a series of regularly repeating units of size ${(0,{\tilde{L}}_y)}^3$, where ${\tilde{L}}_y/{\tilde{L}}_x=ϵ\ll 1$.

• Download figure
• Open in new tab
• Download PowerPoint

To model the interaction between two fluids we use the Cahn–Hilliard model [27–29]. We define the phase field ϕ, a dimensionless variable that takes the value ϕ=1 in fluid 1 and ϕ=0 in fluid 0. At the interface between the two fluids ϕ changes smoothly from ϕ=0 to ϕ=1 over a distance $\tilde{\mathrm{\lambda }}$ which we refer to as the interface thickness. We define the viscosity and density of fluid j as ${\tilde{\eta }}^{(j)}$ and ${\tilde{\rho }}^{(j)}$, respectively, for j={0,1}. We consider the case ${\tilde{\rho }}^{(1)}\gg {\tilde{\rho }}^{(0)}$ corresponding to, for example, air and water. However, we note that this is not a limitation of the model, it is instead used to simplify the notation and algebra used in the remainder of the paper.

To derive the equations that describe the interaction between these two fluids, we write an appropriate fluid free energy which we will then minimize [27–31]. Specifically we write the bulk and surface free energies of the fluid as

where f(ϕ)=ϕ²(1−ϕ)², $\tilde{\gamma }$ is the surface tension and $\alpha =6\sqrt{2}$ scales the free energy such that $\tilde{\gamma }$ is the total excess free energy. It can be shown [40] that h(ϕ) satisfies $h^{\prime }(\varphi )=\sqrt{2}\cos (\theta )\sqrt{f(\varphi )}$, where θ is the contact angle between the fluid–fluid interface and the surface of the porous matrix which is assumed constant and h′(ϕ) is the variational derivative δh/δϕ. The fluid free energy consists of two terms: the bulk free energy which has minima at ϕ=0 and ϕ=1 and the interface energy which acts to minimize the total fluid volume over which ϕ is changing. We now proceed as in [30,31] and derive the momentum equations describing the fluid motion by minimizing the Rayleighian where $\tilde{\mathcal{F}}={\tilde{\mathcal{F}}}_b+{\tilde{\mathcal{F}}}_s$ is the total fluid free energy, ${\widetilde{\mathit{\text{u}}}}^{(0)}$ and ${\widetilde{\mathit{\text{u}}}}^{(1)}$ are the fluid velocities of fluid 0 and 1, respectively, $\widetilde{\mathit{\text{u}}}=\varphi {\widetilde{\mathit{\text{u}}}}^{(1)}+(1-\varphi ){\widetilde{\mathit{\text{u}}}}^{(0)}$ is the combined velocity, $\tilde{\zeta }$ is the drag coefficient between the two fluids, $\tilde{\sigma }=\tilde{\eta }[(\tilde{\mathbf{\nabla }}\widetilde{\mathit{\text{u}}})+{(\tilde{\mathbf{\nabla }}\widetilde{\mathit{\text{u}}})}^{\mathrm{T}}]$ is the combined stress tensor, $\tilde{g}$ is the acceleration due to gravity and $\tilde{\eta }(\varphi )$ is the phase-dependent viscosity which takes the values $\tilde{\eta }(0)={\tilde{\eta }}^{(0)}$ and $\tilde{\eta }(1)={\tilde{\eta }}^{(1)}$ in fluid 0 and fluid 1, respectively. Similarly, $\tilde{\rho }(\varphi )$ is the phase-dependent density which takes the values $\tilde{\rho }(0)={\tilde{\rho }}^{(0)}$ and $\tilde{\rho }(1)={\tilde{\rho }}^{(1)}$ in fluid 0 and fluid 1, respectively. Finally, $\tilde{p}$ is a Lagrange multiplier, effectively a combined or reduced fluid pressure, used to enforce incompressibility of the overall mixture. Assuming conservation of mass for each fluid, we also have Differentiating equation (2.1) with respect to time, using equation (2.3), and assuming the fluid velocity vanishes on the porous structure boundary, we find where $\widehat{\mathit{\text{n}}}$ is a unit vector normal to the surface of the porous medium. We minimize the Rayleighian with respect to $\tilde{p}$, ${\widetilde{\mathit{\text{u}}}}^{(1)}$ and ${\widetilde{\mathit{\text{u}}}}^{(0)}$ and, after some algebra and application of the divergence theorem, obtain the Cahn–Hilliard–Stokes system of equations combined with the no slip condition on the surface of the porous matrix the contact angle boundary condition derived from (2.4) and a zero flux condition to ensure the total volume of each phase is conserved Here $\tilde{\mu }$ is the capillary pressure which is defined in equation (2.5d) and ${\widehat{\mathit{\text{e}}}}_3$ is a unit vector in the direction ${\tilde{x}}_3$. We shall discuss the relationship between the capillary pressure and the specific fluid pressures in §2c. The additional condition, equation (2.5h), ensures the conservation of mass of each fluid on the boundary of the porous structure. Equations (2.5) combined with the initial condition $\varphi (\widetilde{\mathit{\text{x}}},0)=\overline{\varphi }(\widetilde{\mathit{\text{x}}})$ describe the two fluid behaviour in the porous structure.

(b) Non-dimensional equations

We non-dimensionalize equations (2.5) by first scaling space with the microscopic length scale ${\tilde{L}}_y$ such that $\widetilde{\mathit{\text{x}}}={\tilde{L}}_y\mathit{\text{y}}$ and $\mathbf{\nabla }={\tilde{L}}_y\tilde{\mathbf{\nabla }}$. Using this scaling we define the unit cell Y =(0,1)³ composed of a fluid part $\mathfrak{B}$ with boundary $\mathrm{\partial }\mathfrak{B}$. We introduce the non-dimensional velocity $\mathit{\text{u}}={[\mathit{\text{u}}]}^{-1}\widetilde{\mathit{\text{u}}}$, pressure $p={[p]}^{-1}\tilde{p}$, capillary pressure $\mu ={[\mu ]}^{-1}\tilde{\mu }$ and time $t={[t]}^{-1}\tilde{t}$, where

We also introduce the dimensionless capillary number (Ca), Peclet number (Pe), Cahn number (λ) and scaled gravitational force (g): where we have used ${\tilde{\rho }}^{(0)}/{\tilde{\rho }}^{(1)}\sim O(ϵ)$, physically this is equivalent to neglecting the influence of gravity on the phase ϕ=0. Finally, we define the phase-dependent viscosity We note that the definition of the Peclet number, which relates the diffusive motion of the interface to the advection by the fluid, is not, strictly speaking, a conventional Peclet number. However, as this is widely used in the literature we have chosen to keep this terminology [40,36]. Using these equations the scaled and dimensionless Cahn–Hilliard–Stokes equations are where σ=η[(∇u)+(∇u)^T] and M=ϕ²(1−ϕ)². We solve these equations subject to the boundary conditions and the initial condition $\varphi (\mathit{\text{y}},0)=\overline{\varphi }(\mathit{\text{y}})$. Equations (2.9) form a complete non-dimensional description of the two fluid motion in the porous material. These will form the starting point of the homogenization procedure. We have scaled equations (2.9) such that Ca∼O(1) and Pe∼O(1) and the only small parameters in the final equations are ϵ and λ, i.e. the interface is narrow and we are considering a porous geometry with well-defined micro and macro scales. The result of the scaling given in equations (2.6) is that a unit change in μ drives a fluid velocity of order ϵ⁻¹. This velocity corresponds to the movement of the fluid–fluid interface which decays rapidly to zero. Hence, the first non-zero contribution to the scaled velocity is order 1.

We are considering a problem with two different small parameters ϵ, the ratio of the microscopic and macroscopic length scales, and λ, the ratio of the interface thickness to the microscopic length scale, see figure 1. Before we proceed it is useful to discus the role of these two parameters. The first of these, ϵ, is standard in the homogenization literature [9] and will form the basis of the asymptotic expansion we will use in §2c to derive the averaged equations. The second small parameter λ is the non-dimensional interface thickness which must be small with respect to the minimum radius of curvature of the porous structure. We shall show in §3 that, in the limit λ→0, this model reduces to existing models where the fluid–fluid interface location is known (ch. 5 in [3]).

(c) Homogenizing the Cahn–Hilliard fluid equation

We consider the geometry illustrated in figure 1. This geometry consists of a solid structure surrounded by pore space which contains two fluids. Our aim is to derive a set of macroscopic equations which are applicable on the length scale ${\tilde{L}}_x$ and describe the movement of these two fluids averaged over the length scale ${\tilde{L}}_y$, where ${\tilde{L}}_y/{\tilde{L}}_x=ϵ\ll 1$. Due to the separation in length scales, and the periodicity of the geometry, the behaviour of the two fluids is, to first approximation, assumed periodic on the pore scale. Using this assumption and considering the solution to equations (2.9) in successive powers of ϵ we will derive a set of equations which describe the slow variation of these periodic solutions on the lengthscale ${\tilde{L}}_x$.

We define two different spatial coordinate systems; x denotes position on the macroscopic length scale and y denotes position on the microscopic length scale. The key assumption used in the following homogenization procedure is that these two length scales may be treated as independent [9]. Hence, we expand ∇=∇_y+ϵ∇_x, where ∇_x denotes the gradient operator on the macroscopic length scale and ∇_y denotes the gradient operator on the microscopic length scale. We also consider a set of different time scales τ₋₁=ϵ⁻¹t, τ₀=t and τ₁=ϵt. These time scales correspond to the fast equilibration of the fluid–fluid interface, the medium time scale movement of the fluid–fluid interface on the scale ${\tilde{L}}_y$ and the slow variation in saturation due to applied pressure gradients, respectively.

Intuitively it may seem natural, as a first approximation, to neglect terms of order in λ and write λ in terms of ϵ before expanding as in, for example, [41]. However, if we do this, then the leading order solution is ϕ=const and multiple phases cannot co-exist. This is because the terms of order λ multiply the highest derivatives in equation (2.9) resulting in a singular perturbation scheme [42]. In order to accommodate the two phases there must be a region of thickness λ in which the function ϕ changes rapidly. As the interface position can change, it is not straightforward to construct an analytic solution in the interface region and match it to the solution in the regions of constant phase. Therefore, to leading order we must consider terms of order λ such that λ∇²ϕ balances λ⁻¹f′(ϕ). Hence, we expand the velocity, pressure and phase only in powers of ϵ,

We also expand the stress tensor, σ=σ_0y+O(ϵ), and the mobility M=M₀+ϵM₁+O(ϵ²), where and μ=μ₀+ϵμ₁+O(ϵ²), where and We now substitute (2.10) and (2.11) into (2.9) and solve for ascending powers of ϵ.

(i) O(ϵ⁻¹) problem

Substituting equations (2.11) into (2.9) and collecting terms of order ϵ⁻¹ we obtain

and with the boundary conditions and the initial condition ${\varphi }_0(\mathit{\text{x}},\mathit{\text{y}},0)=\overline{\varphi }(\mathit{\text{y}})$ and p₀, μ₀ and ϕ₀ are periodic with period 1. Physically equations (2.12) describe the location of the fluid–fluid interface and are satisfied at steady state for any μ₀ and p₀ which are constant in y, hence, $\parallel \mathfrak{B}\parallel {\mu }_0={\int }_{\mathfrak{B}}{\mu }_0\hspace{0.17em}\mathrm{d}y$, where $\parallel \mathfrak{B}\parallel ={\int }_{\mathfrak{B}}\hspace{0.17em}\mathrm{d}y$. We note in passing that, at steady state, ϕ₀ is a function of both x and y which we write as ${\varphi }_0=S(\mathit{\text{x}})+{\varphi }_0^{(m)}(\mathit{\text{y}})$ where ${\varphi }_0^{(m)}$ is the modulated part of ϕ₀ with zero average. Hence, the saturation is defined as where S takes value S=1 for a fully saturated region and S=0 for a fully unsaturated region. Therefore, we write μ₀∼μ₀[S(x)], where S(x) varies only on the macroscopic scale.

Finally, we observe that, by defining the fluid pressure p_s(x,y)=p₀(x)+Ca⁻¹ϕ₀(y)μ₀(x), we can rewrite equation (2.12b) as

Following the method outlined in [32] we integrate equation (2.14) over a cylinder of height 2hλ, where h≫1 centred about the interface and, after some algebra, obtain the Young–Laplace equation relating the capillary pressure to the pressure drop across the interface: Here ${\widehat{\mathit{\text{n}}}}_{{\varphi }_0}$ is a unit vector normal to the fluid–fluid interface. Hence we can write where $p_s^{(j)}$ is the specific pressure in the jth fluid and, using equation (2.15), we find $p_0=p_s^{(0)}$, i.e. p₀ represents the specific pressure of phase 0.

In order to obtain a macroscopic theory which is valid for all saturation levels we will have to solve equations (2.12) for all possible initial saturation values. In reality this can be achieved using a discrete set of different saturation values and interpolating. It is also clear that the resulting value μ₀(S) is dependent not only on the initial saturation, but also on the initial conditions, $\overline{\varphi }$. For now we shall assume that we know $\overline{\varphi }$ and, hence, μ₀(S) can be determined and will revisit this point in §4.

(ii) O(ϵ⁰) problem

To proceed we collect terms of O(ϵ⁰) from the expansion of equation (2.9). First we consider the expansion of equation (2.9a) and the corresponding boundary condition (2.9g):

and Before we proceed we note that μ₀∼μ₀(x) such that ∇_yμ₀=0 and the terms involving M₁ in equations (2.17) vanish. We now check for solvability by integrating equations (2.17) over $\mathfrak{B}$ and applying the divergence theorem. Hence, by the Fredholm alternative [9] for a solution to equations (2.17a) to exist, we require This is an equation for the τ₋₁ dependence of ϕ₁ if we denote the volume average over the unit cell as $⟨\cdot ⟩={\int }_{\mathfrak{B}}\cdot \hspace{0.17em}\mathrm{d}y$ and integrate in τ₋₁ between 0 and T₋₁ we obtain where we have taken T₋₁≫1 such that ϕ₀ has been in steady state for a long time. In order that ϕ₁ does not grow linearly in time we require that ∂_τ0ϕ₀=0 and, hence, ϕ₀ is independent of τ₀ and ϕ₁ is independent of τ₋₁. The result is a set of equations for u₀, p₁, ϕ₁ and μ₁: with boundary conditions and the correction to the phase can be found using the following equation for ϕ₁, and and are u₀, p₁, μ₁ and ϕ₁ periodic with period 1. We note, however, that we do not need to explicitly calculate ϕ₁ in order to obtain the averaged equations. Equations (2.20) are linear in u₀, p₁, ϕ₁ and μ₁ and depend on x, y and S. Specifically μ₀ is a function of saturation and, hence, x, p₀ is a function of x only and ϕ₀ is a function of y and S. In order to find u₀, p₁, ϕ₁ and μ₁ we consider the effects of the pressure and saturation gradients in equations (2.20) separately. We note that, as ${\varphi }_0=S(\mathit{\text{x}})+{\varphi }_0^{(m)}(\mathit{\text{y}})$ the terms of the form ∇_x⋅∇_yϕ₀=0. We write the solution in the form where ${\chi }_k^j$, ${\mathit{\kappa }}_k^j$, ${\omega }_k^j$ and ${\psi }_k^j$ for j={μ,p,g} are assumed periodic with period 1. We note that the velocity written in the form of equation (2.21) is effectively a generalized Darcy's law. The functions ${\chi }_k^{\mu }$, ${\mathit{\kappa }}_k^{\mu }$, ${\omega }_k^{\mu }$ and ${\psi }_k^{\mu }$ satisfy the cell problem originating from ∇_xμ₀ where ${\sigma }_k^{\mu }=\eta [({\mathbf{\nabla }}_y{\mathit{\kappa }}_k^{\mu })+{({\mathbf{\nabla }}_y{\mathit{\kappa }}_k^{\mu })}^{\mathrm{T}}]$. The functions ${\chi }_k^p$, ${\mathit{\kappa }}_k^p$, ${\omega }_k^p$ and ${\psi }_k^p$ satisfy the cell problem originating from ∇_xp₀ where ${\sigma }_k^p=\eta [({\mathbf{\nabla }}_y{\mathit{\kappa }}_k^p)+{({\mathbf{\nabla }}_y{\mathit{\kappa }}_k^p)}^{\mathrm{T}}]$. The remaining functions χ^g, κ^g, ω^g and ψ^g satisfy the cell problem originating from g: where σ^g=η[(∇_yκ^g)+(∇_yκ^g)^T].

(iii) O(ϵ¹) problem

We now expand the phase equation and the incompressibility condition to O(ϵ) and use the results of the previous expansions to write

and with the relevant boundary conditions and where ϕ₂ and μ₂ are periodic with period 1. As in the O(ϵ⁰) case the final terms in equation (2.25a) which contain multiples of ∇_yμ₀ vanish as μ₀∼μ₀(x). Considering the solvability condition for equations (2.25), as in §2c(ii), we require that ϕ₁ and ϕ₂ do not grow linearly with time on any scale. Hence, the total fluid volume is conserved. Integrating equation (2.25a), applying the divergence theorem, and using equations (2.25b), (2.25c), (2.25d) with equations (2.20c) and (2.20d) we find, after some algebra, Similarly, we integrate equation (2.25b) over $\mathfrak{B}$ and, using the divergence theorem, obtain where the saturation-dependent parameters in equation (2.26) are We have simplified equations (2.28) by taking the limit λ→0. Hence, as M₀ is only non-zero in a region of width λ, ${\int }_{\mathfrak{B}}{M}_0\hspace{0.17em}\mathrm{d}y\sim \mathrm{\lambda }\to 0$. The parameters in equation (2.27) are Equations (2.26) and (2.27), with (2.28) and (2.29), form a pair of coupled equations for the macroscopic saturation and pressure. These are parametrized by the cell problems, equations (2.22), (2.23) and (2.24), which must be solved for a range of saturation.

We note that if ϕ₀=1 everywhere, corresponding to full saturation, then $b(S)=\overline{b}(S)$ is the hydraulic conductivity of phase 1 and equations (2.26) and (2.27) become degenerate. Hence, we obtain Darcy's law for single phase flow. Similarly if ϕ₀=0 everywhere then b(S)=0 and $\overline{b}(S)$ is the hydraulic conductivity of phase 0. Again, we obtain Darcy's law for single phase flow. By applying the assumptions used in deriving Richards' equation, i.e. that the pressure of phase 0 is constant in Ω we can write

Equation (2.30) is the saturation form of Richards' equation which is valid assuming $\overline{K}(S)\ll \overline{b}(S)$ and ${\overline{b}}_g(S)\ll \overline{b}(S)$ such that equation (2.27) is approximately satisfied for constant pressure. We note that we could also have derived the mixed form of Richards' equation simply by leaving equations (2.26) and (2.27) in terms of μ₀. The functions ϕ₀, ${\mathit{\kappa }}_k^{\mu }$, ${\mathit{\kappa }}_k^p$ and κ^gδμ₀/δS must be found for all saturation values. In reality it is enough to compute them for a subset of values and interpolate between them. The function μ₀ is the scaled capillary pressure which is a function of geometry, contact angle and is history dependent. We investigate these effects using two and three dimensional examples in §4.

3. Analysis of homogenized equations

Equations (2.26) and (2.27) are valid across the whole range of saturation values. However, before we study the numerical solution of these equations we use matched asymptotics [42] to consider the limit of the Cahn–Hilliard–Stokes model when λ→0. It has been previously shown that, in this limit, the model reduces to standard free boundary problems for two fluid flow [32–34]. Here, in §3a, we use matched asymptotics to consider the limit of high and low saturation. Then, in §3b, we use these techniques to reduce the cell problems to those derived for fixed interfaces [3,22].

(a) High and low saturation limit

Physically the high and low saturation limits correspond to the case where either the water or the air collects in the sampled geometry to form bubbles which are attached to the soil aggregate surface through capillary pressure. These bubble solutions are found by considering a single droplet of fluid attached to the porous structure. These droplets are assumed to be sufficiently small such that the surface to which they are attached may be considered planar, i.e. the radius of curvature of the droplet is much smaller than the radius of curvature of the porous structure. We start by following the method presented in [34] and considering a steady state radially symmetric solution to equations (2.11d), for μ₀=C and C is a constant, in the absence of the porous structure. The solution is then patched onto the porous structure such that the contact angle which the droplet makes with the surface is correct and, hence, equations (2.12) are solved. Rewriting equation (2.11d) in radial coordinates we obtain

where N is the number of dimensions considered, i.e. N=2 for 2D or N=3 for 3D. We consider the case of a bubble of dimensionless radius r_b, where λ≪r_b. We solve equation (3.1) using the method of matched asymptotics [42]. We define two outer regions, r<r_b and r>r_b and an inner region r∼r_b. In order to obtain the full solution we are required to solve equation (3.1) in each region and match the solutions. However, as we are only interested in obtaining the value of the constant C in terms of radius, it is enough to consider only the leading order solution and the first order correction in the inner region.

We denote the solutions in the outer regions ϕ^(o−) for r<r_b and ϕ^(o+) for r>r_b. The solution in the inner region is denoted ϕ⁽ⁱ⁾. Neglecting terms O(λ) we find f′(ϕ^(o−))=f′(ϕ^(o+))=0. For the inner region we rescale ρ=λ⁻¹(r−r_b) to obtain

with the boundary conditions ϕ⁽ⁱ⁾→ϕ^(o−) as $\rho \to -\mathrm{\infty }$ and ϕ⁽ⁱ⁾→ϕ^(o+) as $\rho \to \mathrm{\infty }$. We expand ${\varphi }^{(i)}={\varphi }_0^{(i)}+\mathrm{\lambda }{\varphi }_1^{(i)}+O({\mathrm{\lambda }}^2)$ and substitute into equation (3.2) to obtain the planar interface solution at leading order, i.e. ${\varphi }_0^{(i)}$ satisfies with ${\varphi }_0^{(i)}\to {\varphi }^{(o-)}$ as $\rho \to -\mathrm{\infty }$ and ${\varphi }_0^{(i)}\to {\varphi }^{(o+)}$ as $\rho \to \mathrm{\infty }$. Solving for ${\varphi }_0^{(i)}$ we find Expanding to order λ we obtain with ${\varphi }_1^{(i)}\to {\varphi }_1^{(o\pm )}$ for $\rho \to \pm \mathrm{\infty }$. Equation (3.5) is a linear equation of the form $\mathcal{L}{\varphi }_1^{(i)}=S$. In order that a solution exists we require that 〈ψ,S〉=0 for $\psi \in \text{ker}({\mathcal{L}}^{†})$, where the superscript † denotes the adjoint. In this case $\mathcal{L}$ is self adjoint and its kernel is spanned by the function $\mathrm{\partial }{\varphi }_0^{(i)}/\mathrm{\partial }\rho$. Therefore, after some algebra, we find The partial bubble solution, where the fluid droplet is attached to a planar surface with contact angle θ, can be derived using simple geometric arguments. The droplet is assumed to be a partial sphere of volume $V=S\parallel \mathfrak{B}\parallel$ for low saturation and $V=(1-S)\parallel \mathfrak{B}\parallel$ for saturation close to 1 hence We will make use of the bubble solution as the limit of high and low saturation in §4.

(b) Comparison with existing models

We now compare our model with the existing two fluid homogenized models which are, for example, presented in ch. 5 in [3]. In these models the interface location is assumed to be known and the interface width is zero. The fluid velocity at the interface is assumed to be continuous in the directions tangential to the interface and zero in the direction normal to it. Finally, these equations are closed by assuming an unknown capillary pressure across the interface. These models can be homogenized using standard methods presented in, for example, [9]. The resulting cell problems can then be used to determine the macroscopic flow properties. In order to compare the equations we have derived to these models we consider the limit of (2.22), (2.23) and (2.24), as λ→0.

We consider cell problem (2.22) in detail as the same procedure can be applied easily to the cell problems (2.23) and (2.24). Far from the interface M₀=0 and equations (2.22) simplify to the Stokes problem

where ${\mathfrak{B}}_1$ is the region in which ϕ₀=1 and ${\mathfrak{B}}_0$ is the region in which ϕ₀=0. We note that at this point ${\chi }_k^{\mu }$ and ${\psi }_k^{\mu }$ are undefined. Equations (3.8) require a pair of boundary conditions at the interface between ${\mathfrak{B}}_0$ and ${\mathfrak{B}}_1$. To determine these conditions we rescale space with the interface thickness, ∇_y=λ⁻¹∇_z, to obtain We find a balance by scaling ${\omega }_k^{\mu }=\mathrm{\lambda }{\overline{\omega }}_k^{\mu }$ and ${\chi }_k^{\mu }={\overline{\chi }}_k^{\mu }\mathrm{\lambda }$. In order to obtain the appropriate interface conditions, we integrate equation (3.9b) over a cylinder of volume V and height 2hλ, where h≫1, centred about the interface. Applying the divergence theorem and taking the limit λ→0 we obtain the condition where we recall ${\widehat{\mathit{\text{n}}}}_{{\varphi }_0}$ is a unit vector normal to the fluid–fluid interface and I is the identity matrix. Similarly, we integrate equation (3.9a) over V to obtain We denote the end facets of the cylinder S0 in the region ϕ₀=0 and S1 in the region ϕ₀=1 and the cylinder side surface ΔS. Applying the divergence theorem, using equation (3.9b) and M₀(1)=M₀(0)=0 gives us the condition where we have used continuity of velocity to simplify the integrals. Finally, observing that variation in ${\chi }_k^{\mu }$ tangential to the boundary is small and may be neglected at this order we obtain the boundary condition Hence, if we define the fluid–fluid interface as $\mathrm{\partial }{\mathfrak{B}}_I$ we obtain the cell problem in the limit λ→0 where ${\hat{\mathit{\tau }}}_q$ for q={1,2} is a unit vector tangent to the interface $\mathrm{\partial }{\mathfrak{B}}_I$. Similarly, we find that ${\mathit{\kappa }}^g=Ca{\mathit{\kappa }}_3^{\mu }$, ${\omega }^g=Ca{\omega }_3^{\mu }$ and ${\chi }^g=Ca{\chi }_3^{\mu }$ and we obtain the cell problem for the pressure driven velocity field

Equations (3.14) and (3.15) correspond to the cell problems derived for a fixed interface in [22]. They can also be derived directly from the microscale equations for a fixed interface presented, for example, in ch. 5 in [3]. Hence, we see that taking the limits with respect to λ and ϵ commute for a fixed interface location.

4. Example

In this section we solve the cell problems (2.12), (2.22) and (2.23) to obtain the macroscopic parameters used in equations (2.26) and (2.27). As shown in §3, in the limit λ→0, the solution to the cell problem (2.24) is proportional to the solution to (2.22), hence we do not consider the solution to this cell problem explicitly. We consider, as an example, the flow of air and water in soil.

We recall μ₀ is the leading order capillary pressure which we derive for a simplified two-dimensional example. The advantage to studying a simplified two-dimensional geometry is that it allows us to easily relate the calculated properties of the capillary pressure to the geometrical features. However, in this case we are interested only in the capillary pressure as the air and water phases can each easily form a plug for a large range of saturation resulting in zero hydraulic conductivity. In order to understand how contact angle affects the capillary pressure, we do this for a range of contact angles. The air–water contact angle with soil has been measured to be approximately 90° [43]. Hence, we consider contact angles of 70°, 90° and 110° as this range of angles provides a clear picture of how the capillary pressure varies near 90°. We then consider a simplified three-dimensional geometry with 90° contact angle and obtain the capillary pressure and hydraulic conductivities.

We note that the solution to the equations (2.12) is dependent on the initial condition. Therefore, the interface location and the capillary pressure are history dependent. In order to accommodate this, we calculate the capillary pressure curves for increasing saturation. We start from the water bubble solution and integrate equations (2.12) to steady state. The saturation is then increased by weakly perturbing the solution and equations (2.12) are solved again using the perturbed solution as the initial condition. By slowly increasing the saturation till the air bubble solution is reached the whole capillary pressure curve can be calculated. The equations are implemented in comsol multiphysics 4.3 using a combination of coefficient form PDEs and computational fluid dynamics modules. The equations are solved using a direct PARDISO method on a single 16 core node of the Iridis 4 supercomputing cluster at the University of Southampton. For both the two-dimensional and the three-dimensional case, the total simulation time is less than 20 h.

(a) Two-dimensional soil

We derive water release curves, which show capillary pressure as a function of saturation, for the geometry shown in figure 2a for a range of different contact angles using the following method. We start by considering an initial partial air bubble attached to the soil. The size of the initial bubble is chosen to give 5% saturation. Equations (2.12) are solved to find the steady-state interface profile and, hence, value of μ₀ corresponding to 5% saturation. The saturation is then increased by 1% and the process is repeated until 95% saturation is reached. The drying curve is captured by reversing the process and decreasing the saturation to 5%.

• Download figure
• Open in new tab
• Download PowerPoint

The water release curve which results from this process is shown in figure 3 for contact angles of 70°, 90° and 110°. These curves are used to parametrize equations (2.26) and (2.27) through the parameter (2.28) and are valid on a timescale much slower than τ₋₁. The bubble solution, equation (3.6), is used for <5% or >95% saturation. The corresponding interface profiles as calculated using equations (2.12) are shown for the 90° contact angle in figure 2. It can be seen that at high and low saturation the water release curve follows the 1/r_b dependence that we would expect for a partial bubble solution.

• Download figure
• Open in new tab
• Download PowerPoint

In the geometry-dependent part of the water release curve, there are several discontinuities, shown with a dashed line. These jumps correspond to the saturation values at which the interface changes topology. As an example we follow the drying curve for the 90° contact angle in figure 3. At high saturation S>0.75, the interface forms a half bubble on the surface of the porous structure, see figure 2a. As the saturation is decreased below S=0.75 the volume of water contained in the bubble becomes to large to fit in the pore and the topology of the solution changes to the one shown in figure 2b. This solution remains valid for 0.75>S>0.43. For S<0.43 the solution topology changes again giving rise to the one shown in figure 2c. This solution remains valid until the air film becomes too thin and the solution collapses to the air bubble solution shown in figure 2d. In our simulations, we have taken this point to be S=0.95.

Increasing the contact angle away from normal incidence acts to increase the capillary pressure whilst maintaining the same set of topological solutions. The overall shape of the water release curve is unaffected by these changes.

(b) Three-dimensional soil

For the three-dimensional case, the same algorithm is used as in the two-dimensional case to determine the water release curve. Once the steady-state interface location has been derived, the cell problems given in equations (2.22) and (2.23) are solved to obtain the hydraulic conductivities. We use ${\tilde{\eta }}^{(0)}=20\times {10}^{-6}\hspace{0.17em}\text{Pa}\hspace{0.17em}\text{s}$ and ${\tilde{\eta }}^{(1)}={10}^{-3}\hspace{0.17em}\text{Pa}\hspace{0.17em}\text{s}$ corresponding to the viscosity of air and water, respectively. For simplicity, we take Ca=Pe=1 and θ=90°.

The water release curve is calculated as in the two-dimensional case. Starting from the bubble solution at 95% saturation we find five topologically different solutions; these are shown in figure 4. The corresponding water release curve is also shown in figure 4 and the diffusivity curves, K(S), b(S), $\overline{K}(S)$ and $\overline{b}(S)$, are shown in figure 5.

• Download figure
• Open in new tab
• Download PowerPoint

• Download figure
• Open in new tab
• Download PowerPoint

As in the two-dimensional case, we see discontinuities in the water release curve where the interface location jumps between topologically different solutions. Following the drying curve we see that the bubble solution, figure 4a, is valid for a small range of saturation values S>0.91 before the droplet spans the gap between adjacent soil particles. At this point the solution topology changes to the one shown in figure 4b, which is valid for 0.91>S>0.68. In contrast to the two-dimensional case, an additional topological state is observed for 0.68>S>0.53, shown in figure 4c. This is observed when the air trapped between a pair of particles expands so much that it interacts with the air trapped between an adjacent pair of particles. As the saturation continues to decrease, the solution changes again to the one shown in figure 4d for 0.53>S>0.05 before collapsing to the bubble solution at low saturation S<0.05.

The corresponding diffusivity curves, neglecting terms of order λ, are plotted in figure 5 for the wetting and drying cycle. The functions b(S), K(S) and $\overline{K}(S)$ are zero for S=0 and increase to give the saturated hydraulic conductivity of the ϕ=1 phase at S=1. The function $\overline{b}(S)$ moves from the hydraulic conductivity of the ϕ=0 phase at S=0 to the hydraulic conductivity of the ϕ=1 phase at S=1.

As with the water release curve, the diffusivities are discontinuous corresponding to the solution switching between topologically different solutions. We note that the diffusivities are non-zero, providing that there is a connected flow pathway for either phase. To illustrate this, we follow the drying curves for all four diffusivities. At full saturation S=1 we obtain simply the diffusivity of phase 1 in all directions and all four diffusivity values are identical. Decreasing the saturation, we observe the half bubble solution figure 4a and the diffusivity in the x direction begins to differ from the solution in the y and z directions. Decreasing further we move through the solution shown in figure 4b,c with discontinuities visible in the diffusivity curves. These are most visible in the K(S) and $\overline{K}(S)$ curves; however, they are present in all four curves. In all four cases, the diffusivity increases as the lower viscosity phase, phase 0, acts to lubricate the flow of the higher viscosity phase, phase 1. At S=0.53 the solution changes to the one shown in figure 4c. At this point, neither phase 0 nor phase 1 is connected in the x direction and the diffusivity values all drop to zero in this direction. The corresponding values of b(S), $\overline{K}(S)$ and K(S) decrease monotonically from this point to 0 at S=0. The function $\overline{b}(S)$ increases monotonically from this point to reach the phase 0 diffusivity at S=0.

The value of $\overline{b}(S)$ is much larger than the other diffusivity curves and in the region 0.05<S<0.95 will dominate the behaviour of the flow. In the high and low saturation regimes, this will not be the case as a(S) grows rapidly will act to amplify the role of K(S) and $\overline{K}(S)$. Hence, in the intermediate saturation region, the final equations (2.26) and (2.27) could be simplified to neglect $\overline{K}(S)$ and ${\overline{b}}_g(S)$ such that equation (2.30) becomes valid under the constant pressure assumption.

5. Discussion

In this paper we have used the method of homogenization to derive, for the first time, a set of macroscopic equations for coupled saturation and pressure-driven fluid flow based entirely on the underlying geometry. The presence of multiple phases is described by the Cahn–Hilliard equation and fluid flow is governed by Stokes equations. The final equations are parametrized by the water release curve and four different diffusivity curves. These curves are captured by solving a series of different cell problems over a range of different saturation values. We have shown that these cell problems, and the resulting macroscopic equations, reduce to standard, simplified homogenizaiton models where the fluid–fluid interface is assumed to be known.

We have used several key assumptions in developing our model. The first is that the porous medium may be considered as a periodic structure. This approximation may be valid for man-made porous structures. However, for natural ones, such as soil, this is clearly only an approximation with the structure being quasi-periodic at best. This assumption has been tested for the case of two fluid flow in imaged soil geometries [23], for the case in which the air–water interface is obtained directly via imaging. In this case, the cell problems were solved on geometries of increasing size and the hydraulic properties were seen to converge. It is expected that the same will apply in this case and, hence, quasi-periodicity is enough for the model to remain valid. Secondly, based on the scaling used in equations (2.9) we require that the capillary forces dominate flow such that there is a separation in time scales between the movement of the fluid–fluid interface and the fluid velocity. This assumption is valid for sufficiently small pores as the fluid velocity shrinks with pore radius whilst the capillary pressure grows. Thirdly, we have modelled the interface between the two fluids using the Cahn–Hilliard equation in which the interface width is assumed non-zero. We have shown that, as the interface width tends to zero, the cell problems derived in §2a converge to those traditionally used for two fluid flow (ch. 5 in [3], [22]). For this approximation to be valid we have implicitly assumed that the interface width used in the numerical calculations is significantly less than the smallest geometrical feature observed. This assumption neglects the effects of small-scale surface roughness which could induce contact angle hysteresis. These effects could, in principle, be included through an effective contact angle-dependent on the small-scale surface geometry.

Using these assumptions we have been able to capture the main features of two fluid flow and, for a given periodic geometry, predict the water release and diffusivity curves. There are three distinct regions observed in these curves, the low- and high-concentration regimes, in which we find an approximate analytic solution for the water release curve, and an intermediate region, in which the water release curve becomes geometry dependent. In this intermediate region the water release curve is discontinuous due to the topologically different solutions obtained at different saturation. Even in simple cases such as those considered in §4 the simulations we have done show that the macroscopic features are strongly related to the geometry studied.

The resulting complexity and discontinuities in both the water release curve and the diffusivity requires some attention. The solutions obtained in §4 are based around a chosen initial condition. However, it is clear that, even for the two-dimensional case, the solutions pictured are not the only possible ones. A different initial condition would have resulted in a different solution structure. For example, the solution shown in figure 2 would be equally valid if the images were simply rotated by 90°. Similarly, the solutions pictured for the three-dimensional case, see figure 4, are equally valid if the solution is rotated by 90° about any of the three coordinate axis. Hence, in general, we expect that the calculated solution will be a combination of topologically different states and that the observed properties of the two fluids will be an average of all possible states.

The water release and diffusivity curves obtained for the sample geometries exhibit hysteresis. This hysteresis is entirely due to the non-linear behaviour of the fluid–fluid interface, i.e. for a given saturation there are multiple stationary solutions. For increasing saturation a different solution is obtained to decreasing saturation. In principle, other sources of hysteresis such as contact angle hysteresis could be included in the model. However, excluding these does not prevent the model from capturing the main observable effects of the two fluid–fluid interaction in a porous geometry.

Finally, we note that the ability to directly predict the water release curves directly from a porous geometry, for example soil, enables a much more precise set of macroscopic equations to be derived without the need for time consuming measurements. While we have demonstrated this method in using parameters appropriate to the flow of air and water in soil, it is applicable to a variety of two fluid systems. In the context of soil, this method, combined with image-based modelling, can be used as a tool to study the effects of different microscopic soil properties on the macroscopic behaviour. This in turn will directly feed back into how soil structure may be optimized in order to maximize flow and transport properties.

Data accessibility

The corresponding author may be contacted for the Comsol implementation of the computational method described in this work.

Author contributions

K.R.D. co-authored the mathematical derivation, conducted all the numerical simulations and wrote the manuscript. T.R. designed the study, co-authored the mathematical derivation of the equations/theory and co-drafted the manuscript. All authors gave final approval for publication.

Funding statement

K.R.D. is funded by BBSRC BB/J000868/1 and T.R. was funded by the Royal Society University Research Fellowship.

Conflict of interests

The authors have no conflict of interests to declare.

Footnotes