Derivation of homogenized equations ( a ) Deriving the two fluid Cahn – Hilliard equation

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.


Introduction
The macroscopic flow of multiple fluid phases in porous media, for example soil, is often described by Richards' equation [1][2][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][5][6] and the numerical solution [7,8].
Mathematically it has been shown, using the method of homogenization [9,10]  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][16][17] and vuggy porous structures [18][19][20][21][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][25][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][28][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][33][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][39][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.

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· 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 S and a fluid part B with boundary ∂B . We assume that B is connected and that ∂B is smooth. The geometry of the porous structure is of typical sizeL x and comprises a series of regularly repeating units of size (0,L y ) 3 , whereL y /L x = 1. To model the interaction between two fluids we use the Cahn-Hilliard model [27][28][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λ which we refer to as the interface thickness. We define the viscosity and density of fluid j asη (j) andρ (j) , respectively, for j = {0, 1}. We consider the caseρ (1) ρ (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][28][29][30][31]. Specifically we write the bulk and surface free energies of the fluid as 2 ,γ is the surface tension and α = 6 √ 2 scales the free energy such thatγ is the total excess free energy. It can be shown [40] 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 RayleighianR whereF =F b +F s is the total fluid free energy,ũ (0) andũ (1) are the fluid velocities of fluid 0 and 1, respectively,ũ = φũ (1) is the combined velocity,ζ is the drag coefficient between the two fluids,σ =η[(Ṽũ) + (Ṽũ) T ] is the combined stress tensor,g is the acceleration due to gravity andη(φ) is the phase-dependent viscosity which takes the valuesη(0) =η (0) andη(1) =η (1) in fluid 0 and fluid 1, respectively. Similarly,ρ(φ) is the phase-dependent density which takes the values ρ(0) =ρ (0) andρ(1) =ρ (1) in fluid 0 and fluid 1, respectively. Finally,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 ∂φ ∂t = −Ṽ · (φũ (1) ) and 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 wheren is a unit vector normal to the surface of the porous medium. We minimize the Rayleighian with respect top,ũ (1) andũ (0) and, after some algebra and application of the divergence theorem, obtain the Cahn-Hilliard-Stokes system of equations and a zero flux condition to ensure the total volume of each phase is conserved Hereμ is the capillary pressure which is defined in equation (2.5d) andê 3 is a unit vector in the directionx 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 φ(x, 0) =φ(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 scaleL y such thatx =L y y and V =L yṼ . Using this scaling we define the unit cell Y = (0, 1) 3 composed of a fluid part B with boundary ∂B. We introduce the non-dimensional velocity We also introduce the dimensionless capillary number (Ca), Peclet number (Pe), Cahn number (λ) and scaled gravitational force (g): where we have usedρ (0) /ρ (1) ∼ 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 [36,40]. Using these equations the scaled and dimensionless Cahn-Hilliard-Stokes equations are ∂φ ∂t where σ = η[(∇u) + (∇u) T ] and M = φ 2 (1 − φ) 2 . We solve these equations subject to the boundary conditions and the initial condition φ(y, 0) =φ(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 −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 §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 scaleL x and describe the movement of these two fluids averaged over the length scaleL y , whereL y /L x = 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 lengthscaleL 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 V = V y + V x , where V x denotes the gradient operator on the macroscopic length scale and V 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 scaleL 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 λ∇ 2 φ balances λ −1 f (φ). Hence, we expand the velocity, pressure and phase only in powers of , (2.10) We also expand the stress tensor, σ = σ 0y + O( ), and the mobility and and We now substitute (2.10) and (2.11) into (2.9) and solve for ascending powers of .
Substituting equations (2.11) into (2.9) and collecting terms of order −1 we obtain with the boundary conditionsn the initial condition φ 0 (x, y, 0) =φ(y) and 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, B μ 0 = B μ 0 dy, where B = B dy. We note in passing that, at steady state, φ 0 is a function of both x and y which we write as 0 is the modulated part of φ 0 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 μ 0 ∼ μ 0 [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 0 (x) + Ca −1 φ 0 (y)μ 0 (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: Heren φ 0 is a unit vector normal to the fluid-fluid interface. Hence we can write s is the specific pressure in the jth fluid and, using equation (2.15), we find p 0 = p (0) s , i.e. 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,φ. For now we shall assume that we knowφ and, hence, μ 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): (2.17b) Before we proceed we note that μ 0 ∼ μ 0 (x) such that V 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 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 τ −1 dependence of φ 1 if we denote the volume average over the unit cell as · = B · dy and integrate in τ −1 between 0 and T −1 we obtain where we have taken 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 : the correction to the phase can be found using the following equation for φ 1 , andn · λV y φ 1 +n · λV x φ 0 + h (φ 0 )φ 1 = 0, y ∈ ∂B, (2.20g) and are 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, y and S. Specifically μ 0 is a function of saturation and, hence, x, p 0 is a function of x only and φ 0 is a function of y and 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 φ 0 = S(x) + φ (m) 0 (y) the terms of the form V x · V y φ 0 = 0. We write the solution in the form where χ The remaining functions χ g , κ g , ω g and ψ g satisfy the cell problem originating from g: (2.24b) κ g = 0,n · M 0 V y χ g = 0, y ∈ ∂B, (2.24c)

(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 with the relevant boundary conditions where φ 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 V 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, that we could also have derived the mixed form of Richards' equation simply by leaving equations (2.26) and (2.27) in terms of μ 0 . The functions φ 0 , κ μ k , κ p k and κ 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.

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][33][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 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 φ (i) . Neglecting terms O(λ) we find f (φ (o−) ) = f (φ (o+) ) = 0. For the inner region we rescale ρ = λ −1 (r − r b ) to obtain Expanding to order λ we obtain for ρ → ±∞. Equation (3.5) is a linear equation of the form Lφ In order that a solution exists we require that ψ, S = 0 for ψ ∈ ker(L † ), where the superscript † denotes the adjoint. In this case L is self adjoint and its kernel is spanned by the function ∂φ (i) 0 /∂ρ. 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 B for low saturation and V = (1 − S) B 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 = 0 and equations (2.22) simplify to the Stokes problem where B 1 is the region in which φ 0 = 1 and B 0 is the region in which φ 0 = 0. We note that at this point χ μ k and ψ μ k are undefined. Equations (3.8) require a pair of boundary conditions at the interface between B 0 and B 1 . To determine these conditions we rescale space with the interface thickness, V y = λ −1 V z , to obtain 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 twodimensional 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 μ 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  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η (0) = 20 × 10 −6 Pa s andη (1) = 10 −3 Pa 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),K(S) andb(S), are shown in figure 5.
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   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) andK(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. The value ofb(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) andK(S). Hence, in the intermediate saturation region, the final equations (2.26) and (2.27) could be simplified to neglect K(S) andb g (S) such that equation (2.30) becomes valid under the constant pressure assumption.

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