Upscaling and Effective Behavior for Two-Phase Porous-Medium Flow using a Diffuse Interface Model

We investigate two-phase flow in porous media and derive a two-scale model, which incorporates pore-scale phase distribution and surface tension into the effective behavior at the larger Darcy scale. The free-boundary problem at the pore scale is modeled using a diffuse interface approach in the form of a coupled Allen-Cahn Navier-Stokes system with an additional momentum flux due to surface tension forces. Using periodic homogenization and formal asymptotic expansions, a two-scale model with cell problems for phase evolution and velocity contributions is derived. We investigate the computed effective parameters and their relation to the saturation for different fluid distributions, in comparison to commonly used relative permeability saturation curves. The two-scale model yields non-monotone relations for relative permeability and saturation. The strong dependence on local fluid distribution and effects captured by the cell problems highlights the importance of incorporating pore-scale information into the macro-scale equations.


Introduction
Flow through porous media, especially in multi-phase systems, is of interest in a variety of applications from oil recovery and CO 2 sequestration to fuel cells and biological systems.Traditionally these are modeled using extended Darcy's law for two phases, which lacks the mathematical derivation from pore-scale information available for the single-phase Darcy's law.Single-phase Darcy's law can be derived from a pore-scale model of (Navier-)Stokes equations through upscaling procedures.For two phases Darcy's law has been extended by introducing empirically derived relative permeability saturation curves to capture fluid interactions.The effective behavior depends only on the averaged saturation and the model fails to capture different pore-scale effects.Further empiric modifications such as play-type hysteresis have been added to account for missing behavior.In this work we start from a two-phase-flow model at the pore scale and through appropriate assumptions and a periodic homgenization approach, we derive a two-scale model of macroscopic Darcy's law-type equations with effective parameters computed from solutions of resulting pore-scale cell problems.
The flow of the two fluids is modeled on the pore scale using quasi-compressible Navier-Stokes equations with phase-dependent density and viscosity as well as an additional term capturing surface tension forces at the fluid-fluid interface.For the pore-scale model with resolved phase distribution we use a diffuse interface approach in the form of an Allen-Cahn phase-field model [5].Where a sharp-interface model separates the pore space into subdomains for each fluid phase, our approach captures the interface implicitly through a smooth phase field function.All unknowns are defined and modified equations hold on the entire domain of the pore space.This removes the need to track or capture the interface and decompose the domain in the numeric simulation.To resolve the transition zone between the two bulk phases a fine discretization is required near the interface, requiring high computational effort in the absence of adaptive local mesh refinement.Furthermore the Allen-Cahn model used in this work is not conservative, including a curvature driven motion of the interface.While the conservative Cahn-Hilliard model requires solving a fourth order partial differential equation, we consider the second order Allen-Cahn equation, which offers a much simpler numerical implementation, and we investigate the two-scale model derived from it.While there are different approaches to ensuring conservation of the phase field [7] or counteracting this curvature-driven motion entirely [26], we apply the scaled mobility approach presented in [1] to eliminate the curvature-driven motion only in the sharp-interface limit.A core advantage of the diffuse interface approach is the ease of upscaling the equations, as they are defined on a stationary domain encompassing both fluid phases.In contrast, upscaling a sharp-interface model requires special attention to the evolving domains [25].
The model is presented for two phases but extends well to more evolving phases [14,20], only a non-neutral contact angle at internal triple points between phases requires more attention [21].
To bridge the scale gap between pore scale and the averaged Darcy scale we employ the method of formal asymptotic homogenization [13].Representing the porous medium by a periodic arrangement of scaled reference cells, one introduces an additional spatial coordinate for this smaller scale.Here the macroscopic scale corresponds to the Darcy scale and the microscopic scale is the pore scale.These two scales are sufficiently separated, with representative length scales L and ℓ respectively defining a length scale ratio ϵ = ℓ/L ≪ 1.Assuming asymptotic expansions of the unknowns in terms of ϵ, considering the limit ϵ → 0 and grouping by orders of ϵ, one obtains equations that are defined on the macro scale and form micro-scale cell problems.The former contains effective parameters which are computed from the cell problem unknowns, linking the two scales.Through the cell problems, these parameters and the effective behavior depend additionally on the distribution of fluids and are able to capture more detailed effects on the fluid velocities.
The homogenization leads to macroscopic equations similar to the extended Darcy's law but with effective parameters computed from solutions of cell problems instead of being given by empiric relations.These effective parameters can be understood as total phase mobilities and for isotropic pore geometries a relative permeability can be computed.Through the cell problems these relative permeabilities in the two-scale model depend on the local phase distribution at the pore scale.This motivates a comparison of computed relative permeabilities to commonly used functional relations, enabling us to investigate under which conditions the models agree and which additional pore-scale effects are captured by the two-scale model.We investigate numerically the effective parameters computed for different geometries and fluid properties.Constructing fluid distributions of varying saturations, we solve the cell problems for the velocity contributions to obtain a relative permeability curve.The models are implemented using a staggered finite volume discretization in DuMu x [15], with Dune-SPGrid [19] and Dune-Subgrid [12] used for the grid.
In [10], [18], and [22] a Cahn-Hilliard model was used to derive a similar two-scale model.We instead investigate the applicability of the simpler Allen-Cahn model and the effects of micro-scale information on the effective flow parameters.The derivation of a two-scale model in [18] additionally assumes a separation of time scales and a different asymptotic expansion to arrive at an instationary problem for the phase distribution.We use the full expansion in the spatial separation ϵ instead and afterwards introduce an artifical evolution to obtain a local phase distribution.In [10] three time scales are used to separate interface equilibration, macroscopic flow and saturation changes.In [22] a solute-dependent surface tension is considered and simulations of the coupled two-scale model are presented.We instead focus on a systematic investigation of the effective parameters computed from the cell problems.A similar comparison of effective parameters is included in [16], where a sharp-interface Stokes model for two-phase flow is upscaled using volume averaging.However, [16] does not account for a three-phase contact line or slip velocities at the solid, which are accounted for in the current work.
In this contribution we aim to derive a two-scale model for two-phase flow in porous media, which captures effects of pore-scale fluid distribution and surface tension on the effective behavior at the macro-scale.Our goal is to determine conditions under which the computed effective parameters show significantly different behavior to commonly used saturation-dependent curves and highlight the importance of resolving pore-scale behavior through a two-scale model for two-phase porous-medium flow.
The structure of this paper is as follows.We present our pore-scale model for two-phase flow using an Allen-Cahn formulation coupled to a modified Navier-Stokes system in Section 2. In Section 3 we perform the upscaling and obtain a two-scale model using periodic homogenization.We investigate the dependence of computed effective parameters on pore-scale conditions using DuMu x in Section 4. Section 5 summarizes the derived model and relates it to the extended Darcy's law.

Pore-scale modeling
We begin by presenting the pore-scale model describing two-phase flow.The model is first presented in a sharp-interface formulation before we introduce the diffuse-interface description for the phase-field model.
The general domain is a stationary void space inside a porous medium.We consider a domain Ω = P ∪ G decomposed into a domain P capturing the pore space and the solid matrix G.Here the distribution of two fluid phases needs to be captured and evolved according to the flow of the fluids, including effects of surface tension at the fluid-fluid interface.

Sharp-interface model
The model is first formulated with a sharp interface separating the fluid phases and corresponding evolving subdomains Ω α (t) ⊂ P, α = 1, 2 of the stationary pore space P (Figure 1).The interfaces between the fluid phases, denoted by Γ f (t) = Ω 1 (t) ∩ Ω 2 (t), as well as the fluid-solid interfaces Γ α = Ω α ∩ G evolve accordingly.Within each subdomain the flow is modelled using incompressible Navier-Stokes equations, with constant viscosity µ α and density ρ α .
At the interfaces Γ α with the solid matrix we prescribe a Navier-Slip [23] boundary condition with slip length λ ≥ 0 for the velocity v α where t is tangential to the fluid-solid interface ∂G.At the fluid-fluid interface Γ f we prescribe continuity of velocities with a no-slip condition enforcing zero tangential velocity, with interface velocity V Γ f in the direction of the interface normal n Γ f pointing from fluid 1 into fluid 2, The interface itself moves due to the local velocity of the fluids and with surface tension γ and curvature H Γ f we have a jump condition where [ψ] = ψ 1 − ψ 2 denotes the jump over the interface.At the three-phase contact point the fluid-fluid-interface meets the solid with a contact angle θ eq (Figure 1).In a sharp-interface formulation this could be incorporated using a boundary condition for a level set equation, when a level set is used to track the location of the fluid-fluid-interface.

Phase-field model
To model two-phase flow at the pore scale we use a diffuse-interface approach, capturing the distribution of phases with a phase-field function u defined on the total void space u : T × P → [0, 1] with u = 0 and u = 1 corresponding to the two distinct phases α = 2 and α = 1 respectively.The phase-field variable evolves according to an advective Allen-Cahn equation, a second order partial differential equation (PDE) with a non-linear source term.The multi-phase system is then modeled as one fluid with varying properties, depending on u(t, x) to account for phase distribution.
We use the compressible Navier-Stokes equations for non-constant viscosity, derived from conservation laws using Stokes hypothesis ζ = 0.A surface tension flux is added to the momentum equation, coupling it to the Allen-Cahn equation as presented for the stationary Stokes equation in [3].In [2] the surface tension term is introduced to a Navier-Stokes equation with a phase field evolved according to the Cahn-Hilliard equation.
The resulting Navier-Stokes equations in conservative form are We focus on the case of each fluid being incompressible, treating the multiphase system as a quasi-incompressible fluid with density and viscosity dependent on the present phase in a linear manner.Denoting mass and viscosity ratios R = ρ 1 /ρ 2 and M = µ 1 /µ 2 , we write This is combined with the advective Allen-Cahn equation, additionally coupled via the advecting velocity and the surface tension momentum flux.The phase-field equation is given as where we consider the scaled diffusivity σξ suggested in [1] in order to obtain a more favorable sharp-interface limit.The diffusivity σ and diffuse-interface width ξ prescribe how strongly a certain profile is enforced for the transition zone between bulk phases and the width of this zone respectively.The double-well potential, P (u), encourages a separation of phases through stable minima at u = 0 and u = 1.The (no-)slip boundary conditions on solid boundaries apply also here, with interpolation based on the phase-field u if phase-dependent slip lengths are desired.Additionally for the phase-field a Neumann boundary condition encoding the contact angle θ eq measured through fluid 1, corresponding to u = 1 (Fig 2 ), is prescribed [11,17].
Note that for vanishing slip length λ = 0 and for neutral contact angle θ eq = π 2 , the boundary conditions reduce to homogeneous Dirichlet and Neumann conditions for velocity and phase field respectively.The Appendix A includes the sharp-interface limit of the coupled Navier-Stokes Allen-Cahn system presented in this section, yielding again the sharp-interface model presented in Section 2.1.by the cell problems, arriving at a coupled two-scale problem.As will be shown, these effective parameters can be interpreted as phase mobilities.
To achieve this we assume a sufficient separation of scales, with characteristic length scales ℓ for the pore scale and L for the Darcy scale at which we are interested in the effective behavior.Assuming a scale separation ϵ = ℓ/L ≪ 1, we first non-dimensionalize the model.We then follow the approach of periodic homogenization, modeling the porous medium as a periodic arrangement of scaled reference cells ϵY , Y = [0, 1] d with dimension d, where d is 2 or 3, (see Figure 3).The domain of the porous medium is thus decomposed into Ω ϵ = ∪ w∈WΩ ϵ(w + Y ), with W Ω ∈ Z d a set of indices.With local pore space and solid matrix denoted as Y = P ∪ G and the outer boundary of the reference cell given as ∂Y we denote the interior boundary between solid and fluid as The global pore space and its boundary with the solid matrix in the non-dimensional case is then given as While their fluid content can generally vary between the cells, we here assume that the porous matrix is constant both in time and space in order to facilitate the derivation of the two-scale model.As the reference cell is the unit cube, we can define the constant porosity φ = |P|.The conditions inside the pore spaces still vary, with fluid saturations and distributions changing in time and between cells.We use the scale separation to introduce a micro-scale coordinate, rewrite spatial derivatives accordingly, assign scalings in terms of ϵ for the dimensionless numbers and assume all unknowns to have an asymptotic expansion in ϵ.Inserting the expansions into the model equations and gathering terms of equal order in terms of ϵ we obtain a new set of equations containing derivatives with respect to coordinates of only one of the spatial scales.This procedure yields macroscopic equations capturing the effective behavior and cell problems defined on the reference cell.Effective parameters are defined through cell problems and integrals of their solutions, linking the two sets of equations and capturing detailed effects of local phase distribution on effective behavior.In addition to pressure driven flow, a separate cell problem captures flow due to surface tension forces, reflected in an added velocity contribution on the effective scale.

Non-dimensionalization
In preparation of upscaling by periodic homogenization we non-dimensionalize the micro-scale model.The homogenization requires a separation of scales between representative length ℓ at the pore scale and the length scale of interest L at the macro scale, quantified by the small number ϵ = ℓ/L ≪ 1. Defining reference values with dimensions and letting one can rewrite the phase-field model equations using non-dimensional variables and parameters In the following all variables and parameters are non-dimensional and the overline notation is dropped for convenience.Using dimensionless numbers one obtains (see Appendix B for details) For the phase-field equation the non-dimensionalization yields (see Appendix B for details) with boundary condition Phase-constant density and viscosity Using the above equations ( 7) for density and viscosity and the reference values ( 12), the non-dimensional viscosity and density in (13) are

Periodic Homogenization
Given the scale separation ϵ we introduce the micro-scale coordinate y = ϵ −1 x and assume all unknowns can be written as an asymptotic expansion in ϵ, depending on both x and y with periodicity in the unit cell Y .For an unknown ψ ∈ {p, v, u} we introduce Y -periodic ψ k (t, x, y) such that The spatial derivatives are rewritten as For the upscaling we consider a flow regime where Darcy's law is considered valid, with laminar flow driven by the pressure drop and capillary forces, and where advective and diffusive time scales are of the same order.For the phase distribution the phase-field diffusivity σ should be comparable to the micro-scale advection, captured by S ∈ O(ϵ 0 ).As will be shown in Remark 1, with S ∈ O(ϵ 1 ) the advective term separates in the upscaling process and yields a restrictive pore-scale equation as well as introducing mixed-scale derivatives into the phase-field equation.
and Fr ∈ O(ϵ 0 ) the leading terms for ϵ → 0 are given as follows (see Appendix C for details).
For the phase dependent fluid properties we observe and for the double-well potential of the phase-field equation, due to its polynomial structure, Denoting Eu := ϵ 2 Eu, we have for (13a), (13b) and (15) in and for boundary conditions ( 14) and ( 16) on Γ

Phase field
From the leading order term of (20c) we therefore obtain a stationary equation involving only local derivatives.
together with boundary condition (20e) Remark 1.If the phase-field equation is dominated by advection (S ∈ O(ϵ 1 )) the leading order advective term of (20e) would be isolated at O(ϵ −2 ) as Applied to the leading order term in (20a) this yields a divergence free velocity field and reduces to resulting in a strong limitation on modeled problems.At the order O(ϵ −1 ) the equation (20c) would instead contain the advective terms of the next order, yielding with an undesireable mix of scales that leads to neither a pore-scale cell problem nor an effective equation for the Darcy scale.A balance of the two terms is expressed by S ∈ O(1), which yields (21) instead of ( 23) as the leading order equation.
This local phase-field equation ( 21) is underconstrained, admitting among others the trivial solutions of u 0 = 0 and u 1 = 1 for divergence free velocity fields.The local cell problem does not offer a way to compute the saturation of fluid 1 as the mean integral of u 0 but rather requires it as a constraint as in [22].
From the next order O(ϵ 0 ) terms of the phase-field equation (20c) we obtain We would like to use this instationary equation to update the saturation to in turn constrain the stationary phase-field equation obtained from the leading order term of the asymptotic expansion.After integrating over the constant local periodicity cell P and using the periodicity with respect to y, (26) reduces to While the first two terms yield a macroscopic equation in saturation and phase-specific velocity, the last term includes the additional unknown u 1 .Assuming a homogeneous porous medium, with P not depending on x, one can use a solvability constraint to show that this term is zero and obtain the desired saturation equation.
Integrating the leading order terms (21) and applying the divergence theorem and periodicity conditions, one is left with only the potential derivative.
One now views the third term in (27) as a Fredholm operator of index zero L(u 0 ) applied to u 1 .
Applying it instead to ∇ x u 0 , using the chain rule and the assumption of P not depending on x, we obtain Together with the previously derived information about P ′ from (28), one sees that ∇ x u 0 is an element of the kernel of L(u 0 ).Rewriting (27 this yields the solvability constraint To avoid trivial behavior (∇ x u 0 ≡ 0, ∇ x u 0 independent of x), the integral over the local porespace P must disappear.Using again the assumption of a stationary and homogeneous pore space P with constant porosity φ, this yields the saturation equation Introducing averaged quantities for the saturation of fluid 1, S (1) = φ −1 P u 0 dy, as well as velocities we obtain The second order term of the mass conservation equation (20a) yields, using (19b) and after integration over P, a second conservation equation.Inserting (34) into (35) yields ∇ x • v = 0 and with S (2) = 1 − S (1) a saturation equation for the second phase, corresponding to u 0 = 0, These macroscopic equations (34), (36), through the saturation, yield an integral constraint for the local phase-field.Together with the stationary equation ( 21) and boundary condition (22) we obtain the cell problem on Γ , u 0 is Y -periodic and P u 0 = φS (1) .
(37) While this ensures mass-conservation it does not fully prescribe a distribution of the phases.The stationary phase-field equation (37) still defines the profile of the diffuse transition zone but the position is not clear.This phase distribution however is central to the equations for fluid flow.To obtain a meaningful solution to this stationary problem one could introduce an artifical time evolution, taking care to introduce additional source terms to the non-conservative Allen-Cahn equation in order to enforce the saturation constraint.That is, one would instead solve (1)  in P , The third term on the right-hand side is introduced in order to make the Allen-Cahn equation conservative [7].The fourth term on the right-hand side moves the system towards the desired saturation and can be localized to the interface using an interface indicator function

Fluid flow
The derivation of the local cell problems for the fluid velocity follows the common approach of separating derivatives of the two scales and using the linear structures of the equations to obtain different problems for the effect of macroscopic pressure gradients, and additionally one to capture the contribution of the surface tension at the interface.
From the leading order term of the momentum equation (20b) we obtain ∇ y p 0 = 0 and from the next order term we get or equivalently, using µ(u 0 (t, x, y)) = 1 + u 0 (t, x, y)(M − 1), Due to the linear structure in the unknowns p 1 and v 0 as well as the contributions of p 0 (t, x) and u 0 (t, x, y), we can find Y -periodic functions w j (t, x, y), Π j (t, x, y), j = 0, . . ., d such that with a function p1 (t, x) independent of y and thus not relevant for (39).These velocity contributions and local pressures can be obtained from the following auxiliary cell problems, where we denote the symmetrized gradient 2ε y (w on Γ , Π j , w j are Y -periodic and P Π j dy = 0, (42) for j ∈ {1, . . ., d} and on Γ , Π 0 , w 0 are Y -periodic and P Π 0 dy = 0, (43) From these solutions we define tensors capturing the effective behavior of the fluid, denoting the phase indicator of phase k as and the components of w j as w j,i .
Multiplying (40) with u (k) and integrating over P, we obtain the macroscopic velocity equations containing the effective parameters. (46)

Two-scale model
In summary we arrive at a micro-macro model, coupling a Darcy-scale flow problem reminiscent of the extended two-phase Darcy's law to d + 2 pore-scale cell-problems on Y = [0, 1] d at every point of the domain.We solve for five macroscopic unknowns p 0 (t, x), S (1) (t, x), S (2) (t, x), v(1) (t, x), v(2) (t, x) with 1 = S (1) + S (2) , and for every set of cell problems the local unknowns u 0 and w i , i = 0, . . .d.
The macro-scale model reminds of Darcy's law with one shared pressure unknown as well as additional effective parameters M (k) which, just as K (k) , are computed from cell-problem variables.For convenience we drop the subscript x on the macroscopic gradient.
We see that the effective parameters K (k) represent the phase-specific effective mobilities, containing information about both the absolute permeability of the porous medium and the interactions between the two fluids.Both contributions are in general anisotropic and can't easily be separated.For isotropic geometries the intrinsic permeability κ abs is a scalar and the effective mobility can be written as κ abs K (k) rel (u 0 )/µ k with the relative permeability rel ∈ R d×d depending on the local phase distribution u 0 through the cell problems (49).The second effective parameter M (k) captures effective flow due to surface tension forces between the two fluids.
In return, the cell problems depend on the local value of the macroscopic saturation through a constraint on the phase-field function.
on Γ , u 0 is Y -periodic and P u 0 dy = S (1) , (48) with v 0 depending on the global pressure gradient and local flow functions w j as defined in (40).These in turn depend on the phase field and are computed through the cell problems defined in (42) and ( 43) on Γ , Π j , w j are Y -periodic and P Π j dy = 0, (49) Figure 4: Staggered discretization for j ∈ {1, . . ., d} and on Γ , Π 0 , w 0 are Y -periodic and P Π 0 dy = 0. (50)

Numeric Investigation
In the following we present a simulation using the tightly coupled pore-scale model from Section 2.2, as well as the investigation of the cell problems for velocity components w j and the behavior of the computed effective parameters.For the former we simulate the movement of two fluids through a channel geometry modeling a single pore.For the latter we consider two exemplary geometries, denoted "cross" and "obstacle", and investigate the relative permeability saturation relation for different fluid properties and fluid distributions.

Implementation
For the numerical investigation all model equations are implemented in DuMu x [15] , using Newton's method as a nonlinear solver and backward Euler for temporal discretization in transient problems.The equations are discretized in space with the finite volume method on a regular rectangular grid, using Dune-SPGrid [19] and an adapted Dune-Subgrid [12] to model the considered geometries with periodic boundaries.For the phase-field and pressure unknowns a cell-centered approach is used, with control volumes equal to grid cells and degrees of freedom placed at their centers.Fluxes over control-volume faces are approximated using the two adjacent values.For the velocities a staggered discretization is used (see Figure 4), with separate degrees of freedom for each velocity component placed at the edges between grid cells and control volumes centered around them.
While the cell problem for the phase field uses only the cell-centered discretization, the pore-scale problem and the cell problems for velocity contributions require a coupled system of cell-centered pressures and staggered velocities.In DuMu x this is achieved using a multidomain formulation with a coupling manager handling the volume coupling of the different discretizations of a shared domain [15].For the pore-scale model the existing coupling manager for Navier-Stokes models was extended to additionally exchange phase-field data.
For the phase-field equation a boundary condition prescribing the flux implements the contact angle condition.In the velocity problems a no-flux condition for the mass conservation and homogeneous Dirichlet conditions for the normal velocity are used at the fluid-solid interfaces.To implement Navier slip boundary conditions a solution-dependent flux can be described.In the case of the cell problem capturing surface tension effects the additional flux can be given based on the prescribed contact angle.On the outer boundary of the reference cell Y periodicity conditions are prescribed.

Pore-scale simulation
We present an examplary simulation of the fully coupled and non-dimensionalized pore-scale model from Section 3.1.In a channel geometry containing both phases the system conforms to a prescribed contact angle and the interface is advected by a pressure gradient.This could be applied to simulating two-phase flow through a narrowing in a porous medium.
We consider a two-dimensional domain Ω = (0, 0.2) × (0, 0.1).At the inlet (x = 0) and outlet (x = 0.2) we prescribe fixed pressures and fluxes for the phase-field equation.The top and bottom boundaries are walls with slip boundary conditions for the velocity and no-flux boundary conditions for the density.For the phase field we apply the contact angle boundary condition (16) with contact angle θ eq = π/3.
We initialize the fluid distribution with fluid 1 on the left and fluid 2 on the right with a curved and diffuse interface between them, see Figure 5.We prescribe a non-dimensional slip length Figure 6 shows the fluid distribution at t = 40s.Due to very weak surface tension forces in this setup, the velocity field does not deviate significantly from a parabolic flow profile.In the center of the channel the interface is advected by these higher velocities, leading to an inversion of the curvature.At the three-phase contact point the prescribed contact angle is successfully maintained despite the near-parabolic profile.

Cell problems
We consider the cell problems (49) and (50) for the velocity components w j and investigate the behavior of the resulting effective parameters.The derived two-scale model is similar to twophase Darcy's law, with effective parameters computed from cell problems rather than prescribed relative permeability saturation curves.We investigate how the computed parameters compare to commonly used curves and under which conditions core assumptions such as a monotone relation is violated.Anisotropic permeabilities are not the focus of this investigation and we choose simple geometries such that the intrinsic permeability of the pore geometry can be separated from fluidfluid interactions.This allows us to isolate the influence of the fluid distribution on the relative permeability in a simple manner.
For two different geometries we vary the local phase field to determine the effects of saturation and phase distribution and the impact of the dimensionless ratios for viscosity and density as well as the surface tension.Both the slip length and the contact angle influence the results, see e.g. the investigation of dynamic contact angles in [17], but will not be considered here.Instead these effects will be controlled by a no-slip condition for the velocity components and a homogeneous neumann condition for the phase field, encoding a neutral contact angle.Since our goal is to address behavior of the effective parameters, we do not couple the flow to the cell problem for the phase distribution (48).Instead we manually prescribe a phase-field function corresponding to a certain saturation and run a few timesteps of the instationary version of the phase-field equation (48) without advection in order to enforce the neutral contact angle.This fluid distribution is then used to solve (49) and (50).
The two domains under consideration vary in geometry and porosity, see Figure 7.The first setup, "obstacle", contains a square of side length 0.45 placed in the center of the domain, leaving a void fraction of 79.75%.The second investigated geometry, denoted "cross", is a cross of channels with diameter 0.3, resulting in a porosity of 51%.
We consider four core scenarios of which the first three investigate the effective mobility K (k) for different ratios of viscosity M and density R. In the last setup the surface tension tensor M (k) is analyzed, with identical fluid properties for both fluids.The dimensionless numbers in the cell problems are chosen as 1.
For all four cases we fix the center of a circle with varying radius r and assign cells inside this circle to phase 1, with the rest of the void space being filled with fluid 2. This initial data already contains a diffuse interface according to the following radial function.
This initial data is evolved under the transient version of the Allen-Cahn equation ( 48) without advection to enforce a neutral contact angle (see Figure 8).The resulting phase field is fed into the different cell problems for flow velocity components (49), (50), which are then solved.
In the "cross" cell geometry we observe recirculation patterns (Figure 14), and e.g. for a horizontal pressure gradient forcing flow from left to right the recirculation near the top and bottom of the domain leads to flow from right to left at these edges of the periodic cell.Depending on the fluid distribution one of the phases may contain primarily these flows in the opposite direction to the main flow, leading to negative integrated velocities.As the effective mobility should capture net flow through the cell, we exclude these recirculations based on the computed flow field, avoiding negative mobilities.The computed velocity field is loaded into ParaView [4], computing streamlines from the boundaries and marking affected cells.This process is not performed for the flow problem driven by surface tension (50).
After postprocessing we compute the weighted integrals over marked cells to obtain the effective parameters.The saturation and an approximation to the specific interfacial area are determined by integrals over the entire domain.
For the effective mobility K (k) , the computed entries are divided by the absolute permeability and multiplied with the viscosity, yielding a relative permeability.The results are plotted over the saturation of phase 1.

Equal fluid properties
In the first case we consider two fluids with equal viscosity and density (M = R = 1).The cell problems for velocity contributions driven by the global pressure gradient treat the two-phase system as single-phase flow.The effective mobilities are obtained through integration weighted by the phase-indicator u (k) and the resulting relative permeabilities always sum to 1, see Figure 9.
Due to the symmetries in the fluid distribution the relative permeability is approximately a diagonal matrix and we show only the dominant diagonal entries.For the isotropic setup in the "obstacle" geometry the entries are equal and the relative permability reduces to a scalar value.
Due to the choice of fluid distribution for every radius r the domain corresponding to fluid 1 contains that for lower radii r ′ < r and thus lower saturation.Equivalently, for phase-field  functions u 1 < u 2 the corresponding saturations fulfil S 1 < S 2 .Under this condition the relative permeability is monotone with respect to the saturation.For differing fluid distributions this monotonicity is not guaranteed and a higher saturation concentrated at areas of lower velocities can result in a lower permeability.Figure 10 shows changing relative permeabilities for fixed saturation and varying fluid distributions obtained by changing the center of the initial circular subdomain for fluid 1.This information about local phase distribution and its impact on macroscopic flow can only be captured by solving the cell problems.

Influence of viscosity differences
Next we investigate the effect of the viscosity ratio by considering M = 2, the investigated cell problems now correspond to an incompressible fluid with varying viscosity.
In stark contrast to typically used relative permeability curves we observe non-monotone behavior for the case "obstacle" with relative permeabilities of the more viscous fluid reaching values above 1 for saturations above 0.8, up to K (1) xx = 1.039 at a saturation of 0.95 (see Figure 11).Due to the chosen fluid distribution the less viscous fluid 2 is located on the surface of the solid matrix at low saturations (high saturation for phase 1).This reduces the resistance exerted  on the first fluid and results in higher velocities, especially for higher viscosity ratios M .Such lubricating effects have been observed experimentally in different settings [6].

Influence of density differences
In the case of non-trivial density ratio R, the cell problem (49) corresponds to Stokes equation for a quasi-compressible fluid.For ratio R = 2 (Figure 12) no non-monotone behavior is observed but the curves are notably different to the one obtained for the case of equal fluid properties.For density ratio R = 10 (Figure 13) non-monotonicity is observed for a horizontal pressure gradient in the "cross" setup, with the velocity in the lighter fluid increasing as the main flow path is filled with the more dense fluid.The fluid distribution and computed flow field for highest permeability at a saturation of 0.45 are shown in Figure 14.In Figure 15 the results for the "obstacle" geometry are compared to those for equal fluid properties, which highlights the fact that the sum of the relative permeabilities is less than one.

Surface tension tensor
For the isotropic geometry and fluid distribution considered above the velocities driven by surface tension cancel out in the integral and the computed effective parameter is equal to 0 (Figure 16, left).For asymmetric fluid distributions (Figure 16, right) a non-zero contribution is obtained, but no visible trends are observed as the saturation changes.The redistribution of fluids driven by surface tension forces leads to a net flow for the two phases, which corresponds to M (k) being non-zero.The size and direction of this net flow and hence the size and sign of the effective parameter M (k) depends highly on the fluid distribution.Such effects can only be accounted for by solving the cell problem (50).

Discussion
For equal fluid properties and fixed fluid distributions we obtain monotone relative permeability saturation curves comparable to commonly used relations such as Brooks-Corey [8] or van Genuchten [24].However, even in this simple case the relative permeability depends not only on saturation but also on the distribution of the fluids.As observed in Section 4.3, even in this simple case no relation can be given without accounting for the local distribution of the fluids.
For differing viscosities the model is able to reproduce non-monotone relative permeabilities with values above 1 as discussed in [6].When the less viscous fluid coats the solid surface it can reduce the resistance exerted on the more viscous fluid, resulting in higher velocities.Due to the different fluid distribution such a coating is not observed to the same extent for the "cross" setup, leading to monotone relations for the relative permeability.These effects depend strongly on the local phase distribution and cannot be captured solely by the saturation without additional assumptions.
For moderate density differences the combined permeability is reduced and for higher density ratios non-monotone behavior can be observed.We remark that for very high differences in fluid properties, this should be considered in the derivation of the cell problems and the twoscale problem, which would lead to a different appearance of the effective equations.To obtain equations similar to the often used two-phase Darcy's law, we use the assumptions presented in Section 3.2.
The derived two-scale model contains an additional effective parameter which accounts for the influence of surface tension between the fluid phases.This is not a common part of the extended Figure 16: Symmetric and asymmetric case with center at (0.35, 0.3).
Darcy's law and there is no counterpart to compare it to.A similar cell problem and resulting effective parameter has been investigated by [22], incorporating the influence of solute-dependent surface tension.Here only constant surface tension is considered.For symmetric distributions the effective parameter capturing surface tension effects disappears as expected.For anisotropic phase distributions the flow driven by surface tension does not cancel out and the associated cell problem is able to capture significant contributions to effective flow behavior.The size and direction of flow resulting from surface tension and thus the magnitude and sign of the effective parameter depend strongly on the local distribution of fluids.
The numeric investigation of the effective parameters highlights how different fluid distributions with equal saturation can result in very different net flow and effective behavior.Increased relative permeability due to viscosity differences and significant surface tension effects further emphasize the importance of resolving the local fluid morphology.

Conclusion
We derived a two-scale model for two-phase flow in porous media using homogenization and investigated the dependence of effective parameters on the fluid distribution in the pore-scale cell problems.
The model consists of macro-scale equations similar to the extended Darcy's law with effective parameters computed from solutions to cell problems defined at every macroscopic point.One of the effective parameters can be understood as an effective mobility, while the second effective parameter represents effects of interfacial tension on the effective flow.On the pore-scale the phase distribution is captured by a phase field, using an advective Allen-Cahn formulation.For the velocity field we use a stationary Stokes equation with phase-dependent fluid properties and an additional source term incorporating surface tension forces.
The two-scale structure allows to investigate numerically the effects of pore-scale information on effective behavior.We construct local fluid distributions for the pore scale and solve the cell problems associated with pressure driven velocity contributions.
The similarity to two-phase Darcy's law motivates a comparison of the dependence of relative permeabilities on saturation.By selecting isotropic pore geometries, relative permeability tensors can be extracted from the effective parameters K (k) in our model and for isotropic fluid distributions these simplify to diagonal matrices or scalar values.As one of the core qualities of relative-permability-saturation-curves we investigate the conditions under which a monotone relationship is obtained.If a common fluid morphology is maintained for different saturations, then for equal fluid properties monotone curves, similar to commonly used functions like van Genuchten and Brooks-Corey, are obtained.The calculated curves however depend strongly on the fluid distribution and change with the pore geometry, highlighting the need to account for pore-scale information.For different viscosities the model can produce non-monotone relationships and relative permeabilities greater than 1.This can be related to a lubricating effect of the less viscous fluid coating the surface of the solid matrix, an effect observed experimentally.For moderate density ratios the effective total permeability decreases and no non-monotone behavior was observed.
We also performed numeric simulations for the cell problem corresponding to the effective parameter M (k) , capturing effects of surface tension forces.For isotropic geometry and fluid morphologies these forces result in no net flow and a vanishing effective parameter.For asymmetric fluid distributions we obtain significant net flows captured by M (k) .
Our two-scale model for two-phase flow captures pore-scale fluid-fluid interactions related to surface tension forces as well as differences in viscosity and density.Through numeric investigation of the effective parameters we highlight the importance of local fluid distribution for accurately capturing effective flow behavior.such that x = yξ(t, s) + rnξ(t, s) . (54) Note than ∂ t r = −v n , with v n the velocity of the interface in direction of nξ corresponding to the sharp-interface velocity V Γ f .Furthermore it can be shown (see [9]) that for mean and Gaussian interface curvatures κ and Π |∇r| = 1, ∇r • ∇s i = 0, ∇ 2 r = κ + 2Πr 1 + κr + Πr 2 .

A.1 Outer expansions
Due to the polynomial structure of P ′ it can be expanded as Substituting the outer expansions (53) into the phase-field equation ( 8) yields The leading order equation P ′ (u out 0 ) = 0 has solutions u out 0 ∈ {0, 1/2, 1}, of which 0 and 1 are stable minimizers of P .This allows a decomposition of the unified P into bulk domains Ω out 1 and Ω out 2 for fluid 1 and 2, respectively, through We consider the remaining outer expansions for Ω out i .For the mass balance equation (6a) we obtain by using the linear dependence of density on the phase-field variable Denoting phase velocities we recover the individual mass conservation equations and, due to constant phase densities ρ i , the divergence-free condition (1).The momentum balance (6b) yields, using additionally the linear relation for viscosity, With phase pressures p i given by restricting p out 0 to the bulk domain Ω out i we recover the momentum balance (2).
Inserting the outer expansions into the boundary conditions yields the conditions for the fluid-solid interfaces of the sharp-interface model.From (9a) we obtain corresponding to the slip condition (3).

A.2 Inner expansions
To obtain the boundary conditions at the fluid-fluid interface we insert the inner expansions (55), rewriting derivatives according to (56), and use the matching conditions (57).
Phase-field equation From the phase-field equation ( 8) we obtain due to the polynomial form of P Multiplying the leading order terms with ∂ z u in 0 and integrating with respect to z yields a condition for the interface velocity where the terms in the first integral are equal to 0 due to the matching conditions.Taking the limit in z, we obtain With v in 0 • n 0 independent of z, the leading order terms reduce to Multiplying the remaining terms with ∂ z u in 0 and integrating with respect to z yields, applying the chain rule and matching conditions with u in 0 (t, y 1/2− ) = 0 and P (0) = 0, By definition of the curvilinear coordinates we have ∂ z u in 0 ≥ 0 and Using u out (t, 0, s) = 1/2 one can obtain a profile for the phase-field, namely Mass balance Inserting the inner expansions into the mass balance equation (6a) yields 0 Integrating with respect to z and applying the matching conditions we obtain 0 which is fulfilled for v 1 = v 2 = v n,0 at the interface.
Momentum balance Inserting the inner expansions into the momentum balance equation (6b) yields with leading order and second order Using the definition of the outer product, the leading order terms at O( ξ−2 ) can be written as Let t 0 the tangential vector, with t 0 • n 0 = 0. Taking the dot product of the above (74) with Integrating with respect to z and using matching conditions, we obtain Since µ ̸ = 0 we have ∂ z v in 0 • t 0 = 0, corresponding to the sharp-interface boundary condition and together with the assumption ∂ z v in 0 • n 0 = 0 the velocity v in 0 is independent of z.Using also the definition of the outer product as well as ∇ Γ ψ • n 0 = 0 , ∇ Γ • n 0 = κ 0 and ∇ Γ u in 0 = 0, the second order terms simplify to Integrating with respect to z and applying the matching conditions yields where [ψ] = ψ 1 − ψ 2 denotes the jump of ψ across the interface and we used the boundedness of ∇ Γ v in 0 as well as ∂ z u in 1 .As the integral evaluates to 2/(3L) and the phase velocities are equal to v n,0 at the interface, we recover the boundary condition (5).

B Non-dimensionalization
As described in Section 3.1, in addition to the two length scales ℓ and L with ϵ = ℓ/L, we define reference values with dimensions Using the non-dimensional numbers and dropping the overline notation this is written as Dividing by ξ2 / t, using the dimensionless number S and dropping the overline notation the equation becomes For the contact-angle boundary condition we obtain

C Homogenization
We introduce the micro-scale coordinate y = ϵ −1 x with ϵ = ℓ/L the scale separation and assume asymptotic expansions for ψ ∈ {p, v, u} with each ψ k y-periodic on the reference cell Y .Rewriting spatial derivatives according to We note that due to the dependence of ρ on u we have and the equation can be written as Analogously for the viscosity µ we have µ 0 = µ(u 0 ) and µ 1 = u 1 (M − 1).Denoting Eu := ϵ 2 Eu ∈ O(ϵ 0 ), the momentum equation (13b) yields the following terms.

Figure 2 :
Figure 2: Contact angle boundary condition for the diffuse-interface phase-field model.

Figure 3 :
Figure 3: Separation of scales and domain definitions

Figure 5 :
Figure 5: Initial condition of pore-scale simulation of two-phase flow through a channel.u = 1 (red) corresponds to fluid 1, u = 0 (blue) to fluid 2.

Figure 8 :
Figure 8: Preprocessing to enforce a neutral contact angle.Initial condition from (51) (left), and the resulting phase field with desired contact angle (right).

Figure 10 :
Figure 10: Impact of phase distribution on relative permeability for equal fluid properties with radius r = 0.15 (left) and r = 0.3 (right).

Figure 14 :
Figure 14: Phase distribution and velocity component for density ratio R = 10.