## Abstract

The macroscopic behaviour of air and water in porous media is often approximated using Richards' equation for the fluid saturation and pressure. This equation is parametrized by the hydraulic conductivity and water release curve. In this paper, we use homogenization to derive a general model for saturation and pressure in porous media based on an underlying periodic porous structure. Under an appropriate set of assumptions, i.e. constant gas pressure, this model is shown to reduce to the simpler form of Richards' equation. The starting point for this derivation is the Cahn–Hilliard phase field equation coupled with Stokes equations for fluid flow. This approach allows us, for the first time, to rigorously derive the water release curve and hydraulic conductivities through a series of cell problems. The method captures the hysteresis in the water release curve and ties the macroscopic properties of the porous media with the underlying geometrical and material properties.

## 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 *Ω* comprising a solid matrix

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 *j* as *j*={0,1}. We consider the case

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
*f*(*ϕ*)=*ϕ*^{2}(1−*ϕ*)^{2}, *h*(*ϕ*) satisfies *θ* 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
*c*. 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

### (b) Non-dimensional equations

We non-dimensionalize equations (2.5) by first scaling space with the microscopic length scale *Y* =(0,1)^{3} composed of a fluid part *Ca*), Peclet number (*Pe*), Cahn number (λ) and scaled gravitational force (*g*):
*ϕ*=0. Finally, we define the phase-dependent viscosity
*σ*=*η*[(**∇ u**)+(

**∇**)

*u*^{T}] and

*M*=

*ϕ*

^{2}(1−

*ϕ*)

^{2}. We solve these equations subject to the boundary conditions

*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

*ϵ*

^{−1}. 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 §2*c* 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 *ϵ* we will derive a set of equations which describe the slow variation of these periodic solutions on the lengthscale

We define two different spatial coordinate systems; ** x** denotes position on the macroscopic length scale and

**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***∇**=

**∇**

_{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

*τ*

_{−1}=

*ϵ*

^{−1}

*t*,

*τ*

_{0}=

*t*and

*τ*

_{1}=

*ϵ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

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 λ∇^{2}*ϕ* balances λ^{−1}*f*′(*ϕ*). Hence, we expand the velocity, pressure and phase only in powers of *ϵ*,
*σ*=*σ*_{0y}+*O*(*ϵ*), and the mobility *M*=*M*_{0}+*ϵM*_{1}+*O*(*ϵ*^{2}), where
*μ*=*μ*_{0}+*ϵμ*_{1}+*O*(*ϵ*^{2}), where
*ϵ*.

#### (i) *O*(*ϵ*^{−1}) problem

Substituting equations (2.11) into (2.9) and collecting terms of order *ϵ*^{−1} we obtain
*p*_{0}, *μ*_{0} and *ϕ*_{0} 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 *μ*_{0} and *p*_{0} which are constant in ** y**, hence,

*ϕ*

_{0}is a function of both

**and**

*x***which we write as**

*y**ϕ*

_{0}with zero average. Hence, the saturation is defined as

*S*takes value

*S*=1 for a fully saturated region and

*S*=0 for a fully unsaturated region. Therefore, we write

*μ*

_{0}∼

*μ*

_{0}[

*S*(

**)], where**

*x**S*(

**) varies only on the macroscopic scale.**

*x*Finally, we observe that, by defining the fluid pressure *p*_{s}(** x**,

**)=**

*y**p*

_{0}(

**)+**

*x**Ca*

^{−1}

*ϕ*

_{0}(

**)**

*y**μ*

_{0}(

**), we can rewrite equation (2.12b) as**

*x**h*λ, 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:

*j*th fluid and, using equation (2.15), we find

*p*

_{0}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 *μ*_{0}(*S*) is dependent not only on the initial saturation, but also on the initial conditions, *μ*_{0}(*S*) can be determined and will revisit this point in §4.

#### (ii) *O*(*ϵ*^{0}) problem

To proceed we collect terms of *O*(*ϵ*^{0}) from the expansion of equation (2.9). First we consider the expansion of equation (2.9a) and the corresponding boundary condition (2.9g):
*μ*_{0}∼*μ*_{0}(** x**) such that

**∇**

_{y}

*μ*

_{0}=0 and the terms involving

*M*

_{1}in equations (2.17) vanish. We now check for solvability by integrating equations (2.17) over

*τ*

_{−1}dependence of

*ϕ*

_{1}if we denote the volume average over the unit cell as

*τ*

_{−1}between 0 and

*T*

_{−1}we obtain

*T*

_{−1}≫1 such that

*ϕ*

_{0}has been in steady state for a long time. In order that

*ϕ*

_{1}does not grow linearly in time we require that ∂

_{τ0}

*ϕ*

_{0}=0 and, hence,

*ϕ*

_{0}is independent of

*τ*

_{0}and

*ϕ*

_{1}is independent of

*τ*

_{−1}. The result is a set of equations for

*u*_{0},

*p*

_{1},

*ϕ*

_{1}and

*μ*

_{1}:

*ϕ*

_{1},

*u*_{0},

*p*

_{1},

*μ*

_{1}and

*ϕ*

_{1}periodic with period 1. We note, however, that we do not need to explicitly calculate

*ϕ*

_{1}in order to obtain the averaged equations. Equations (2.20) are linear in

*u*_{0},

*p*

_{1},

*ϕ*

_{1}and

*μ*

_{1}and depend on

**,**

*x***and**

*y**S*. Specifically

*μ*

_{0}is a function of saturation and, hence,

**,**

*x**p*

_{0}is a function of

**only and**

*x**ϕ*

_{0}is a function of

**and**

*y**S*. In order to find

*u*_{0},

*p*

_{1},

*ϕ*

_{1}and

*μ*

_{1}we consider the effects of the pressure and saturation gradients in equations (2.20) separately. We note that, as

**∇**

_{x}⋅

**∇**

_{y}

*ϕ*

_{0}=0. We write the solution in the form

*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

**∇**

_{x}

*μ*

_{0}

**∇**

_{x}

*p*

_{0}

*χ*

^{g},

*κ*^{g},

*ω*

^{g}and

*ψ*

^{g}satisfy the cell problem originating from

*g*:

*σ*

^{g}=

*η*[(

**∇**

_{y}

*κ*^{g})+(

**∇**

_{y}

*κ*^{g})

^{T}].

#### (iii) *O*(*ϵ*^{1}) problem

We now expand the phase equation and the incompressibility condition to *O*(*ϵ*) and use the results of the previous expansions to write
*ϕ*_{2} and *μ*_{2} are periodic with period 1. As in the *O*(*ϵ*^{0}) case the final terms in equation (2.25a) which contain multiples of **∇**_{y}*μ*_{0} vanish as *μ*_{0}∼*μ*_{0}(** x**). Considering the solvability condition for equations (2.25), as in §2c(ii), we require that

*ϕ*

_{1}and

*ϕ*

_{2}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,

*M*

_{0}is only non-zero in a region of width λ,

We note that if *ϕ*_{0}=1 everywhere, corresponding to full saturation, then *ϕ*_{0}=0 everywhere then *b*(*S*)=0 and *Ω* we can write
*μ*_{0}. The functions *ϕ*_{0}, *κ*^{g}*δμ*_{0}/*δ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 *μ*_{0} 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 *μ*_{0}=*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
*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 *ϕ*^{(i)}. Neglecting terms *O*(λ) we find *f*′(*ϕ*^{(o−)})=*f*′(*ϕ*^{(o+)})=0. For the inner region we rescale *ρ*=λ^{−1}(*r*−*r*_{b}) to obtain
*ϕ*^{(i)}→*ϕ*^{(o−)} as *ϕ*^{(i)}→*ϕ*^{(o+)} as *ψ*,*S*〉=0 for *θ*, can be derived using simple geometric arguments. The droplet is assumed to be a partial sphere of volume

### (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}=0 and equations (2.22) simplify to the Stokes problem
*ϕ*_{0}=1 and *ϕ*_{0}=0. We note that at this point **∇**_{y}=λ^{−1}**∇**_{z}, to obtain
*V* and height 2*h*λ, where *h*≫1, centred about the interface. Applying the divergence theorem and taking the limit λ→0 we obtain the condition
*I* is the identity matrix. Similarly, we integrate equation (3.9a) over *V* to obtain
*S*0 in the region *ϕ*_{0}=0 and *S*1 in the region *ϕ*_{0}=1 and the cylinder side surface Δ*S*. Applying the divergence theorem, using equation (3.9b) and *M*_{0}(1)=*M*_{0}(0)=0 gives us the condition
*q*={1,2} is a unit vector tangent to the interface

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 *μ*_{0} 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 2*a* 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 *μ*_{0} 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%.

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 *τ*_{−1}. 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.

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 2*a*. 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 2*b*. 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 2*c*. This solution remains valid until the air film becomes too thin and the solution collapses to the air bubble solution shown in figure 2*d*. 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 *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*),

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 4*a*, 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 4*b*, 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 4*c*. 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 4*d* 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 *S*=0 and increase to give the saturated hydraulic conductivity of the *ϕ*=1 phase at *S*=1. The function *ϕ*=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 4*a* 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 4*b*,*c* with discontinuities visible in the diffusivity curves. These are most visible in the *K*(*S*) and *S*=0.53 the solution changes to the one shown in figure 4*c*. 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*), *K*(*S*) decrease monotonically from this point to 0 at *S*=0. The function *S*=0.

The value of *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

## 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.

## Acknowledgements

The authors acknowledge the use of the IRIDIS High Performance Computing Facility, and associated support services at the University of Southampton, in the completion of this work. The authors would also like to thank Prof. J. Yeomans FRS for helpful discussions relating to this work.

- Received July 23, 2014.
- Accepted January 29, 2015.

© 2015 The Authors. Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited.