Formulation of a phase field model of multiphase flow in deformable porous media

In this contribution a new model for partially saturated porous media is presented, within the framework of gradient poromechanics, based on a phase field approach. While the standard retention curve is expected still to provide the intrinsic retention properties of the porous skeleton, depending on the porous texture, an enhanced description of the surface tension between the wetting and the non-wetting fluid, occupying the pore space, is stated considering a regularized phase field model based on an additional contribution to the overall free energy depending on the saturation gradient. This approach provides similar results as those of the model based on the concept of specific interfacial area. The dependence of the free energy on the gradient of the saturation also implies that Darcy law must be extended to become a fourth order partial differential equation.


Introduction
The constitutive characterization of partially saturated porous media became of interest at the middle of the last century when the scientific research started to face fundamental problems in geotechnics and petroleum engineering, concerning the response of partially imbibed soils, during drainage-imbibition cycles, or modeling the behavior of sedimentary reservoir rocks, when a multiphase fluid flows through the porous space.Starting from the analysis of basic static problems, it became clear that the balance between capillary and driving forces, in particular gravitational forces, would have been the central subject of the modeling efforts.This pushed the research in the direction of finding out a relation between the curvature of the wetting/non-wetting fluid interface and the average content of the wetting fluid, say the retention curve for the macro-scale capillary pressure.For long time this last has been the only relation used to describe the hydraulic flow through partially saturated porous media, being also the pivot of the hydromechanical coupling with the constitutive law of the porous skeleton, see e.g.Alonso et al. [1].
In the same period Cahn & Hilliard [2] established the basic framework of modeling multi-phase flow in terms of space and time evolution of a phase field, continuously varying over thin interfaces.Surface tension is recovered, in this context, considering the integral, through the thickness of the layer, of the phase field gradient.This approach progressively attracted more and more interest, in particular within fluid mechanics, because of its advantages for numerical calculations; however limited contributions attempted at incorporating these ideas into modeling of unsaturated porous media [3].
In this paper a novel general approach is developed which aims at merging phase field modeling of multiphase flow with unsaturated strain gradient poromechanics.While the standard retention curve is expected to still provide the intrinsic retention properties of the porous skeleton, depending on the porous texture, an enhanced description of surface tension between the wetting and the non-wetting fluid, occupying the pore space, is stated considering a regularized phase field model.

Thermodynamics
Thermodynamics of porous media has been summarized by Coussy [4] for two or more monophasic superimposed interacting continua, say the solid skeleton and the fluids saturating the porous space.In this case, the specific internal energies of the fluids are separately defined, whether they are in the liquid or in the gaseous phase; whilst the energy due to interfacial interactions between the fluids and among the solid and the fluids are incorporated into the macroscopic energy of the skeleton.Here a novel approach is adopted in order to account, within the framework of poromechanics, for the role of the interface between a liquid (water) and a gas (wet air), describing the mixture as a non-uniform diphasic fluid, which can reside in the liquid or in the gaseous phase.The fluid internal energy is given by: where fef is a double-well potential parametrized by the mass density, per unit volume of the mixture, f and the specific entropy sf, n being the Eulerian porosity.The a Corresponding author: giulio.sciarra@uniroma1.itnon-local energy f penalizes the formation of interfaces and provides a regularization of the non-convex energy fef.The state equation of the fluid defines the thermodynamic pressure P and the chemical potential , in terms of the the mass density f: however, because of equation ( 1), these quantities are not sufficient to completely characterize the constitutive law of the non-uniform fluid, which, on the other hand, necessitates of specifying the so-called fluid hyper-stress vector, related to the derivative of f with respect to (nf),  .In the following incompressibility of the liquid phase is assumed which implies the variation of f to be univocally determined by the variation of the saturation degree Sr, f = LSr, L being the intrinsic constant density of the liquid phase.The liquid (gaseous) phase therefore corresponds to Sr=1 (Sr=0).
From now on attention is focused on isothermal processes, a Duffing-like potential for the volumetric energy is assumed: where nw is the surface tension between the non-wetting and the wetting phase and R the average characteristic size of the channel through which the non-uniform fluid can pass.Let moreover the non-local term f be quadratic in the gradient of nSr.The first and the second principle of thermodynamics, together with the expression of the internal working relative to a strain gradient porous continuum [5], allow to state the constitutive prescriptions for the overall Piola-Kirchhoff stress and hyper-stress, Sij and Pijk, required to verify the overall balance of momentum (Sij -Pijk,k),j + bi = 0, as well as the generalized constitutive characterization of the capillary pressure Pc and the fluid hyper-stress vector k.After cumbersome calculations, see [7], one gets  (0) being the Lagrangian (initial) porosity, and J the determinant of the Jacobian matrix, relative to the deformation of the overall porous continuum.S'ij = s/Eij and P'ijk = s/Eij,k are the generalized effective stresses, prescribed in terms of the partial derivatives of the skeleton free energy s =  -f with respect to strain and strain gradient,  being the free energy of the overall porous continuum.Equations ( 4)-( 6) hold true within the assumption of small elastic strains.Thermodynamics also implies the fluid dissipation to be positive which yields the generalized Darcy law: -P,k/Sr + (Pc -(l/(nSr)),l),k + bk f /(Sr) = Akl  Ml/L (7) where Ml is the Lagrangian filtration vector and Akl the inverse of permeability.Introducing the capillary energy U, so that s/Sr= U/Sr, implies equation ( 7) to be rephrased as follows: The role of the capillary energy U is therefore that of modifying the double-well potential f which prescribes the free energy of the fluid, in order to account for the wetting properties of the solid skeleton.The new free energy f+U, which can be called effective energy of the pore-fluid, has not the same minima as f, as they are shifted inward the interval (0,1) from below or from above, whether the solid skeleton is gas or liquid wet.

Characterization of the pore-fluid
Providing a constitutive characterization of the pore-fluid, say of the fluid within the pore network, is generally achieved, in unsaturated poromechanics, assuming the capillary pressure Pc to be a prescribed function of Sr.
Here the liquid-gas mixture, which saturates the pore space, is regarded as a non-uniform fluid, the corresponding saturation ratio being used to characterize the state of the fluid at any current placement.No distinction is made explicit between the pressure of the liquid and the pressure of the gaseous phase, as well as between the corresponding chemical potentials.Moreover no explicit algebraic relation between the saturation degree and the capillary pressure or the effective chemical potential of the fluid  eff , defined in what follows, can a-priori be stated, without solving, at least at equilibrium, the generalized Richards equation:  eff = (f+U)/Sr -(s/Sr)l)l (10) Here the permeability of the porous medium has been assumed isotropic: (A -1 )kl = Ksatk(Sr)ij, Ksat and k(Sr) being the so-called saturated and relative permeabilities of the wetting phase, respectively.At the stationary state, equation ( 9) reduces to a fourth order partial differential equation in the space variable, which is definitely similar to the one prescribing the mass density distribution of a Cahn-Hilliard fluid at equilibrium.However a fundamental additional term is here accounted for, say the derivative of U with respect to Sr, which allows for describing the confining effect on the non-uniform fluid, due to the presence of the porous skeleton.The role of f and U in the characterization of the distribution of Sr,  and Pc deserves a deeper discussion.Consider the chemical potential of the pure fluid, prescribed by equations ( 2  This non-monotonic behavior of the pore-fluid chemical potential resembles that one postulated by DiCarlo [7] in order to explain the formation of gravity fingers in soils.However it is worth to underline that the appropriate form of the potential -U/Sr should be captured solving an inverse problem which stems from experimental evidence.
Notice that the effective generalized chemical potential  eff can be hydrostatic at equilibrium, and consistent with gravitational loading, even if the porefluid chemical potential  + U/Sr is not.

Micro-scale interpretation of the generalized constitutive law of capillary pressure
Equation ( 6) explicitly improves the basic constitutive prescription of the macroscale capillary pressure by a correction depending on the gradient of the Jacobian determinant J and the gradient of the water content nSr, once the quadratic form f for the nonlocal contribution to the energy of the fluid has been assumed.As a consequence when the hyperstress acting on the non uniform fluid vanishes, this correction vanishes as well; in other words the larger the gradient of the liquid content is, the wider the discrepancy between the standard and the enhanced constitutive prescriptions of Pc will be.The enhanced constitutive prescription of Pc is therefore expected to yield significant modifications of the capillary pressure, in the narrow subdomains of the current shape of the porous medium where significant gradients of the liquid content can be detected.This justifies an interpretation of this contribution in terms of the socalled specific interfacial area, see [8].
According with experimental data, reported among others in [9], and porenetwork numerical simulation, see e.g.[10], local variations of the macroscopic capillary pressure are accompanied not only by changes in the saturation degree but also by changes in the average interfacial area anw, which accounts for the local cumulative measure of the interfaces between the non wetting and the wetting phase (per unit volume of the Representative Volume Element, RVE).Even if changes in the interfacial area may not result into a (incoming/ outcoming) flow through the RVE, however, in this analysis, only such kind of variations of anw are considered.No account, on the other hand, is taken of those variations of the interfacial area which yield a pure remodeling of the internal structure of the RVE.This assumption will be clarified in the following when the homognization scheme adopted to upscale information from the micro to the macro scale is introduced and briefly discussed.Its rationale however resides in the identification of a relation between the interfacial area and the gradient of the saturation degree, to hold true only if the liquid particles are displaced to the boundary of the RVE.To verify the validity of this hypothesis a heuristic microscale analysis is developed deforming a prototype reference configuration B 0 of a RVE by a quasistatic loading path driven by the Lagrangian gradient of the macroscale Jacobian determinant J; the liquid is squeezed out of the porous chamber along the same direction of the gradient of J.The saturation degree Sr is assumed not to be affected by the considered microscale deformation, which means that the ratio between the volume of the liquid and the volume of the pores within the RVE does not vary during the deformation process.This assumption allows to capture the correction to the constitutive prescription of Pc provided by the gradient term of the fluid energy only.In Figure 2   Each RVE corresponds to the deformed current configuration B of a reference domain B 0 constituted by four identical beads among which a suitable amount of liquid water (the wetting phase) is trapped by capillary forces.The amount of liquid (in volume) is prescribed by the value of the degree of saturation; on the other hand its spatial distribution is a function of the local wetting properties of the beads.Here, for the sake of simplicity, the same equilibrium Young contact angle is postulated at each loading step, moreover the interface between the nonwetting and the wetting fluid is kept constantly circular, during the evolution process.Tuning the parametrizing value of Sr provides aposteriori the relation between the specific interfacial area and the saturation degree, and consequently allows to compare the results of the model with the benchmark experimental data, relative to drainageimbibition cycles.A suitable homogenization procedure, which stems from the averaging operations summarized in [11], can be developed in order to provide a motivation of the conjectured relation between the average interfacial area anw and the macroscopic capillary pressure.The average based homogenization scheme implies the corrective terms in equation ( 6) not to be zero only if the fluid trapped among the beads attains the boundary, during deformation.In this way the gradient of the liquid content does not fade away.Thus in order to prove a functional correlation between the gradient of Sr and the interfacial area, all the microscale deformations which do not drive the liquid to the boundary of the RVE are not taken into account.
Considering the geometrical data reported in the caption of Figure 2, the corrective term of capillary pressure, introduced in equation ( 6), is explicitly calculated using the above mentioned homogenization techniques.At the same time the average interfacial area consistent with the envisaged liquid distribution, can be estimated for each value of J ,1 .A contour plot of the constitutive law of Pc, parametrized by the value of anw, can therefore be drawn, see Figure 5, which qualitatively resembles that one traced in Figure 11 of [10], for the values of the saturation degree which have been considered in the present analysis.
The considered deformation and flow regimes are however totally different with respect to those simulated in [10].As a matter of facts they are driven by the gradient of the Jacobian J, at different values of Sr, whilst those ones of [10] are deduced simulating many scanning loops of drainage and imbibition.This is indeed an interesting point which corroborates the idea, introduced in [8], of resolving the hysteresis of the capillary pressure, between drainage and imbibition, introducing a suitable microstructural parameter, say anw, which tunes the value of Pc for each value of Sr.
It is worth to underline that the considered microscale analysis stems from the assumption of nonnegligible strain gradient, which in fact parametrizes the deformation process, even if only small strain of the solid skeleton are considered.If, on the other hand, this assumption were not valid, equation ( 6) whould abruptly simplify into the standard constitutive prescription of capillary pressure.

Conclusions
In this paper gradient theory of poromechanics endowed with phase field modeling have been used to describe partial saturation.The spatial distribution of saturation, and strain, is regularized at phase coexistence, within the nonuniform fluid, and in the neighbors of possible heterogeneities in the porous skeleton.To do this the free energy of the overall porous medium has been regarded as a function not only of strain and saturation but also on their (Lagrangian) gradients.
As usual in gradient theories, the governing partial differential equations relative, in this case, to the solid skeleton and the nonuniform fluid are deduced using integration by parts twice, which implies the equations to be, in general, of the fourth order.An enhanced version of classical Richards' equation has been deduced from generalized Darcy's law, which in the considered regime of partial saturation does not depend only on the capillary pressure, but also on the socalled generalized chemical potential.
A novel constitutive characterization of the capillary pressure is established and a microscale interpretation is provided, which is consistent with that originally formulated in [8].A comparison with experimental results also confirms that the contribution to capillary pressure due to the saturation gradient allows for recovering similar effects as those captured by the average interfacial area.

9 E 2016 -
)-(3), say = f/Sr = Lghf, hf = 2C (nw/(LgR)) Sr (1-3Sr+2Sr 2 )(11) and the derivative of the capillary energy, relative for DOI: 10.1051/ , instance to a sand, and given by the van Genuchten curve U/Sr = -LghU, hU = (1/){[(Sr-Sr res )/(1-Sr res )] -1/m -1}1/n  (12)Sr res being the residual value of saturation; the two curves are depicted in Figure1(upper panel, gray dotted and gray dashed lines, respectively).Both the chemical potential and the derivative of the capillary energy are expressed in head units adopting the classical Leverett scaling, which prescribes the characteristic length R in terms of the intrinsic permeability of the soil  and the initial porosity 0:R = (/0)1/2 The profile of the negative chemical potential of the pore-fluid, say -(+U/Sr), plotted against Sr, is also drawn (solid line).It exhibits a non-monotonic behavior and two additional zeros (the spots) with respect to the one at Sr=1, which can always be found, when considering only the capillary pressure Pc.This feature corresponds to the fact that the energy f + U maintains the double-well shape typical of f, see Figure1(lower panel).However the two minima of f + U are no more isopotential and the one associated to the smallest value of Sr is shifted inwards the interval (0,1).

Figure 1 .
Figure 1.In upper panel heads of the chemical potential  of the pure fluid (dotted gray line), of the derivative of the capillary energy (dashed gray line) and of the negative porefluid chemical potential -(+U/Sr) (solid black line) are depicted; in the lower panel the corresponding energies.The parameters which characterize the retention curve are those of a sand.
possible reference configurations of the RVE, parametrized by the degree of saturation are depicted.

Figure 2 .
Figure 2. Referential RVEs parametrized by the saturation degree: (a) Sr= 0.667298; (b) Sr= 0.834704; the dashed line indicates the boundary of the RVE.The distance between the beads, made dimensionless with respect to the characteristic size of the RVE, is L=0.25, whilst their radius is R=(1-L)/2.

Figure 3 .
Figure 3. Admissible distributions of the liquid phase within different RVEs, obtained as the current configurations of the reference shape of panel (a).The micro-scale displacement of the beads is parametrized by the value of J,1 via its dependence on the macro-scale second gradient of deformation.During the deformation the saturation ratio is kept constant: Sr= 0.667298.

Figure 4 .
Figure 4. Admissible distributions of the liquid phase within different RVEs, obtained as the current configurations of the reference shape on panels (a).The micro-scale displacement of the beads is parametrized by the value of J,1 via its dependence on the macro-scale second gradient of deformation.During the deformation the saturation ratio is kept constant: Sr= 0.834704.

Figure 5 .
Figure 5.The capillary pressure Pc is represented as a function of the saturation ratio Sr, within the interval delimited by the considered maximum and minimum value of Sr and parametrized by the interfacial area anw.