Phase field modeling of partially saturated deformable porous media

Phase field modeling of partially saturated deformable porous media

Giulio Sciarra Dipartimento Ingegneria Chimica Materiali Ambiente,
Università di Roma La Sapienza, Via Eudossiana 18, 00184 Rome, Italy

A poromechanical model of partially saturated deformable porous media is proposed based on a phase field approach at modeling the behavior of the mixture of liquid water and wet air, which saturates the pore space, the phase field being the saturation (ratio). 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 surface tension between the wetting (liquid water) and the non-wetting (wet air) fluid, occupying the pore space, is stated considering a regularization of the phase field model based on an additional contribution to the overall free energy depending on the saturation gradient. The aim is to provide a more refined description of surface tension interactions.

An enhanced constitutive relation for the capillary pressure is established together with a suitable generalization of Darcy’s law, in which the gradient of the capillary pressure is replaced by the gradient of the so-called generalized chemical potential, which also accounts for the “force”  associated to the local free energy of the phase field model. A micro-scale heuristic interpretation of the novel constitutive law of capillary pressure is proposed, in order to compare the envisaged model with that one endowed with the concept of average interfacial area.

The considered poromechanical model is formulated within the framework of strain gradient theory in order to account for possible effects, at laboratory scale, of the micro-scale hydro-mechanical couplings between highly-localized flows (fingering) and localized deformations of the skeleton (fracturing).

Phase field, Poromechanics, Strain gradient, Capillarity


The constitutive characterization of partially saturated porous media became of interest at the beginning of the last century when scientific research started to face fundamental problems in geotechnics and petroleum engineering, concerning modeling the response of partially imbibed soils, during imbibition/drainage cycles (see Buckingham (1907); Richards (1931)), or modeling the behavior of sedimentary reservoir rocks, when a multi-phase fluid flows through the porous network. 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 modeling efforts. This pushed the research in the direction of finding out simple relations between the curvature of the wetting/non-wetting fluid interface and the average content of the wetting fluid, over a suitably defined Representative Volume Element (RVE). The main ideas of this identification have been clearly sketched by (Coussy, 2010, Chapter 6) for a partially saturated truncated conical pore, specifying the infinitesimal variation of the interfacial energy in terms of the infinitesimal variation of the volume occupied by the wetting fluid (saturation). This naturally provided the well-known definition of macro-scale capillary pressure, as the derivative of a macro-scale capillary energy. In Leverett (1940), Brooks and Corey (1964), in the summarizing contribution by Bear (1972), in van Genuchten (1980) etc., different semi-empirical relations, between macro-scale capillary pressure and saturation have been stated, in order to specify the retention properties of the porous skeleton. These retention curves typically exhibit hysteresis during imbibition/drainage cycles, so as to fit with experimental evidence.

At the same time the pioneering papers by Cahn and Hilliard Cahn and Hilliard (1958); Cahn (1959); Cahn and Hilliard (1959) established the basic framework within which modeling of multi-phase fluid flow is formulated in terms of space and time evolution of the mass concentration or, in the general case, of a phase field which can vary continuously over thin interfacial layers. Beyond the original formulation due to Cahn and Hilliard (1958) one can refer to Jacqmin (2000) and Kim (2012) for general reviews, as well as to Lowengrub and Truskinovsky (1998), Boyer and Lapuerta (2006), Boyer et al. (2010) or Lamorgese et al. (2011). Surface tension is recovered, in this context, considering the integral, through the thickness of the layer, of the concentration gradient (or the gradient of the phase field). This approach progressively attracted more and more interest, in particular within the framework of fluid mechanics, because of its advantages for numerical calculations which do not necessitate adaptive interface fitting grids, see e.g. Jacqmin (2000); Kim (2012); Gomez and Hughes (2011). Surprisingly however limited contributions attempted at incorporating these ideas into modeling of unsaturated porous media, one can refer for instance to Papatzacos (2002) and Papatzacos and Skjoeveland (2004) and more recently to Cueto-Felgueroso and Juanes (2008, 2009b, 2009a) and Gomez et al. (2013).

This paper follows the research path traced by these authors to model partially saturated deformable porous media, considering the pore network infused with a two-phase fluid. Also the constitutive characterization of the solid will be generalized with respect to the classical poromechanical model, in order to describe the coupling between highly localized flows and (possibly localized) strains.

For long time the above mentioned constitutive characterization of the macro-scale capillary pressure in terms of the corresponding hysteretic retention curve has been the only relation used for describing the hydraulic flow through partially saturated porous media, being also the pivot of the hydro-mechanical coupling with the constitutive law of the porous skeleton. Due to the coarse simplification provided by this model, however, it was finally recognized the existence of several problems in which the relation between capillary pressure and saturation is not sufficient to describe at the same time surface tension, between the wetting and the non-wetting fluid, and the retention properties of the porous skeleton, due to its texture. As a matter of fact what is understood by this constitutive law is, as already mentioned, that no variation of the capillary pressure can be observed if no change of the saturation is induced by external forces. This is due to the rough upscaling rule, from the micro to the macro-scale, which assumes the variation of the area of interfaces proportional to the volume occupied by the fluids, see Coussy (2010). On the other hand, as originally pointed out by Morrow (1970), redistribution of the fluids within a porous medium can occur also keeping constant the saturation, and instabilities can even arise from fluid interfaces that are unable to change curvature smoothly with variations in pressure (“Haines jumps”). In other words, different micro-scale configurations of interfaces between a wetting and a non-wetting phase are possible for a given saturation.

In more refined macro-scale formulations the specific cumulative measure of interfaces, between the wetting and the non-wetting phase, say the specific interfacial area, is also introduced to account for the micro-scale features of the fluids within the RVE, see Hassanizadeh and Gray (1993) and the related literature, e.g. Reeves and Celia (1996); Papatzacos (2002); Niessner and Hassanizadeh (2008). In order to clarify the meaning of this additional macro-scale state parameter, one can think of the effect of the saturation degree, on the macro-scale capillary pressure, as that which coarsely accounts for the retention characteristics of the porous skeleton, depending on its texture, and consider on the other hand the effect of the cumulative measure of interfaces as a suitable corrective term, which allows for describing different admissible coexistence configurations of the saturating fluids. Within this framework, several experimental campaigns have been carried out, based on laboratory tests, see e.g. Dalla et al. (2002); Culligan et al. (2004); Costanza-Robinson et al. (2011), and also pore network micro-scale numerical simulations have been implemented, see e.g. Joekar-Niasar et al. (2010), in order to characterize the constitutive relation among capillary pressure, saturation and specific interfacial area.

Consider now the behavior of the porous skeleton, it is almost standard in continuum poromechanics to introduce a constitutive prescription of stress in terms of strain and saturation, both in the case of elastic and elasto-plastic deformations; one can refer among others to the seminal paper by Alonso et al. (1990) and the related literature, see e.g. Olivella et al. (1996); Alonso et al. (1999), to the systematic formulation by Lewis and Schrefler (1998), see also Sanavia et al. (2002, 2006), to the approach based on the generalized Bishop effective stress, see e.g. Houlsby (1997); Gens et al. (2006); Buscarnera and Einav (2012); Jommi (2000); Tamagnini (2004) etc., or to some recent contributions within the range of finite deformations of the porous skeleton, see Borja et al. (2013) and Song and Borja (2014). Possible modifications of these constitutive models, in order to account for the effect of the specific interfacial area, concern the improvement of the retention curve provided within the approach to multi-phase mechanics proposed by Hassanizadeh and Gray (1990, 1993), in which velocities of phases and interfaces are employed to describe fluid flow and skeleton deformation, see also Gray and Schrefler (2001); Gray et al. (2009); Nikooee et al. (2013). Within this framework, however, several phenomena, related to micro-structure remodeling, can not be captured, exactly as several phenomena related to the spatial distribution of fluids within the pores could not be described only in terms of saturation evolution. The typical case is that of localized deformations which have been observed in laboratory scale experimental tests, since the beginning of the 70s, initially using medical scanners for imaging sand bodies, see e.g. Vardoulakis et al. (1978) and Desrues et al. (1985), and later with the aid of RX computerized tomography and Digital Image Correlation, see e.g. Desrues and Viggiani (2004); Lenoir et al. (2007); Hall et al. (2010) etc. Thanks to the improvements in the resolving power of tomographs as well as to the progressively increasing skills of scientists in the treatment of experimental data, dry and (partially) saturated specimens have been imaged under triaxial loading conditions to observe localized deformations and porosity change, see e.g. Andò et al. (2012a, b) and Desrues and Andò (2015).

In order to model the effects at the scale of the specimen of these micro-scale grain displacements, and in particular to obtain the regularization of localized deformation patterns, the gradient theory of elasto-plasticity, as well as the gradient approach to damage and fracture mechanics, have been adopted in the literature. Typical references are Sulem and Vardoulakis (2004), Chambon et al. (2001) and Matsushima et al. (2002), as well as Francfort and Marigo (1998) and Bourdin et al. (2000, 2012). In all these cases however the hydro-mechanical coupling, if any, still remains restricted to the constitutive prescription of the effective or net stress, see e.g. Collin et al. (2006) and Sieffert et al. (2009).

As already mentioned, in this paper a novel general approach is developed, within the framework of continuum poromechanics, which aims at merging phase field modeling of multi-phase fluid flow with unsaturated strain gradient poromechanics. 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 surface tension between the wetting and the non-wetting fluid, occupying the pore space, is stated considering a regularized phase field model (in the sense of Cahn and Hilliard (1958)). Similarly to the model which introduces the specific interfacial area as a corrective term of the standard retention curve, this alternative formulation allows for detecting variations of capillary forces even when no change in saturation is caused by the external loading. Indeed a similar approach has already been addressed in the literature, attempting at modeling gravitational fingering and saturation overshooting, through an undeformable porous skeleton, or viscous fingering through a Hele-Shaw cell, see Cueto-Felgueroso and Juanes (2009a, 2012, 2014). However no account has been taken of the possible coupling with the deformations of the porous skeleton. In the model which is going to be presented in this paper, on the other hand, the behavior of the solid skeleton is described by means of a strain gradient model, to predict localized strains and fracture, when considering irreversible processes. To the best of author’s knowledge, no general model accounting for both the above mentioned items, say localization in fluid flow and strain, can be retrieved from the literature, except a preliminary study by the author Sciarra (2013), stemming from previous results on modeling of porous media saturated by quasi-incompressible fluids, see Dell’Isola et al. (2003); Sciarra et al. (2007, 2008) and Cirillo et al. (); Cirillo et al. (2010, 2011, 2012, 2013). Within this enhanced framework, the hydro-mechanical coupling is therefore responsible not only for the effects of average variations of the hydraulic regime within the RVE on the skeleton deformation, and vice-versa, see among others Olivella et al. (1996); Alonso et al. (2010); Tamagnini (2004); Casini (2012); Rotisciani et al. (2014, 2015), but it can also account for the effects of capillary fingering on damaging and fracturing of the skeleton and vice-versa for the effects of strain localization on permeability variations and heterogeneous/anisotropic fluid flow.

The paper is organized as follows: in section 1 the basics of kinematics of (partially) saturated porous media are summarized. In section 2 the expression of the virtual working of external forces, relative to a second gradient multi-phase continuum, is manipulated to obtain the virtual working of internal forces, which is consistent with the balance of the overall momentum. In section 3 the first and the second principle of thermodynamics are introduced in order to get the generalized form of the Clausius-Duhem inequality, prescribing the solid, fluid and thermal dissipations. In section 4, the constitutive laws for the effective stress and hyper-stress are established, considering the contribution of strain and saturation gradients to the overall free energy of the porous continuum; the generalized Darcy law is also deduced. Section 5 is devoted to discuss the properties of the so-called pore-fluid, constitutively characterized in terms of a suitable double well potential, specifying the two phases of the fluid within the pores. In section 6 the generalized prescription of the macro-scale capillary pressure is placed, accounting for its dependence on the saturation gradient. In section 7 the governing equations, say the overall balance law and the generalized Darcy law, are deduced by means of a variational approach which also provides the proper boundary conditions. Finally in section 8 some conclusions are stated.

1 Kinematics

A macroscopic description of kinematics is adopted, treating the partially saturated porous medium as the superposition of a skeleton and a binary mixture of liquid water and wet air. The current placement of the skeleton particles is provided by the deformation of a reference configuration . The partial saturation of the pore space is accounted for thinking of the fluid mixture as a non-uniform fluid in the sense of Cahn and Hilliard (1958), i.e. a fluid possibly having a spatial variation in one of its intensive scalar properties.

Let be the -th component of the deformation (placement) of the skeleton particles with respect to a fixed orthonormal frame in the Euclidean space of positions; is the reference configuration of the solid constituent and a time interval. The image of under the deformation map is the current configuration of the porous medium, say . Moreover let indicate the displacement of the solid particles, the deformation gradient, along the -th direction of an orthonormal frame in the reference configuration of the solid, and the associated strain tensor, being the Kronecker delta. From now on Greek indices label the components of any -th order tensor with respect to a frame in the Euclidean space, whilst Latin indices label the components of any -th order tensor with respect to a frame into the reference configuration of the skeleton . Finally let be the -th component of the velocity of the skeleton particles. As usual within the framework of poromechanics, the notions of Eulerian porosity and Lagrangian porosity are introduced, as the current pore volume density per unit volume of the porous continuum, and the current pore volume density per unit reference volume, respectively. Clearly ; being the determinant of .

Concerning the non-uniform fluid the intensive scalar property, used from now on to characterize the biphasic characteristic of the fluid, is not the mass concentration but the volume density of the liquid phase with respect to the volume of the pores, say the degree of saturation. The liquid phase of the mixture (water) is assumed incompressible and the gaseous phase (wet air) is assumed to be passive (that is, have infinite mobility) so that its density can be neglected with respect to that of the liquid, ; the apparent density of the fluid mixture, say the mass density of the non-uniform fluid per unit volume of the porous medium, is therefore definitely prescribed in terms of the volumetric liquid content . If indicates the saturation degree (saturation ratio), measuring the current volume occupied by the liquid per unit volume of the pores, the liquid content is and the fluid apparent density is . The saturation degree ranges within the interval , where and indicate the gaseous and the liquid phase, respectively. The kinematics of the non-uniform fluid is specified as in Sciarra et al. (2008); the reference placement of the non-uniform fluid particles is therefore prescribed by means of a regular map defined on the reference configuration of the skeleton. It identifies the fluid material particle within a reference domain , which occupies, at time , the same current place as the solid particle, chosen in . Accordingly the fluid velocity coincides with the time derivative of , which means for every reference particle and every current place :


where ; the Einstein summation for repeated indices is understood in equation (1). Similarly to the deformation gradient of the solid skeleton also the gradient of can be defined in terms of the gradient of and as .

Both the saturation degree and the liquid content are naturally defined on the current configuration of the porous medium, however the corresponding pull-back in the reference configuration of the skeleton can easily be defined by means of the inverse of the deformation map . The mass conservation of the fluid mixture can be stated, following Coussy (2004), in the form


where is the Lagrangian fluid mass content and is the Lagrangian filtration vector; the time derivative is computed keeping fixed the reference placement of the solid particle. According with previous remarks, the incompressibility of the liquid, say , implies equation (2) to reduce to:


In equations (2)-(3), for the sake of simplicity, no explicit dependence of the considered fields on the current position or the corresponding placement in the skeleton configuration has been specified. It is worth to underline that even if the liquid is incompressible the non-uniform fluid is not, consequently no restriction on the fluid velocity is placed.

2 External & strain working: from a mixture model towards a Biot-like model

Here an approach much similar to that developed by Sciarra et al. (2007) is adopted, in order to deduce the strain working relative to a porous medium within the framework of gradient continuum mechanics. Starting from the standard formulation of the mechanics of superimposed continua, and using the almost classical formulation of gradient theories introduced by Toupin (1962), Mindlin (1964) and Germain (1973), the external working of the solid-fluid mixture is defined as a continuous linear functional of the velocity fields and :


where is the boundary of the current configuration of the porous continuum, assumed differentiable almost everywhere, and is the union of the edges of , on which the normal to the boundary suffers a jump. Equation (4) implies that not only bulk forces and surface tractions , but also double forces and tractions per unit line on both constituents , are expected to be be balanced by proper internal tractions. This balance can be achieved considering the extended Cauchy theorem, see Germain (1973); Dell’Isola et al. (2009), which for the -constituent reads


and are tensorial quantities defined over the current configuration, which represent the stress and the so-called hyper-stress, acting on the solid and the fluid. Let , indicate the local parametrization of the boundary , the tensor field is the projection tensor onto the tangent space of , while the partial derivatives with respect to -coordinates indicate the surface-gradient; are the components in the surface coordinate system of the Darboux tangent-normal vector to each edge of . Finally indicates the jump through the edge. Equations (5) imply that surface tractions depends on the curvature of , say the surface gradient of the normal unit vector , and, in the limit when this curvature tends to infinity, say over the edges of , that a line density of forces must arise. Thinking of the pure solid, double forces generalize the concept of skew-symmetric couples, which is typical of Cosserat or couple-stress theories, bearing into account not only surface density of moment but also surface density of symmetric couples. The former, working on a combination of differential shearing, bending and torsion, the latter, on differential elongation. On the other hand, thinking of non-viscous but non-uniform fluids, characterized by internal capillarity, say Cahn-Hilliard fluids, double forces are associated only to symmetric couples working on differential liquid content or, within our framework, saturation ratio, which allows to describe surface tension effects, see e.g. Seppecher (1989, 1996). In this case, which is the one considered in this paper, the hyper-stress tensor relative to the fluid reduces to .

In order to deduce a Biot-like poromechanical model the overall balance of momentum must hold true and the corresponding strain working must be characterized in terms of strains and stresses defined over the reference configuration of the skeleton. Using equations (5) and (5) together with the orthogonal decomposition of the velocity gradients, given by , the external working (4) reduces to the following form:


Let and be the -th and the -th components of the overall stress and hyper-stress tensor, respectively, let moreover be the -th component of the overall bulk force per unit volume; in order for the overall porous medium to be balanced the following equation must hold true


for every place . A consistent definition of the strain working of the porous skeleton can deduced from equation (6):


where, from now on, the following positions are assumed: and .

Following Dell’Isola et al. (2009) the generalized Piola-Kirchhoff overall stress tensors and are defined as the pull-back in the reference configuration of the porous skeleton of the current stress and hyper-stress tensors:


Accordingly the strain working is re-written as


where indicates the time derivative, keeping the placement of the solid particles fixed in the reference configuration. As already mentioned the fluid is a non-uniform Cahn-Hilliard fluid, consequently the fluid hyper-stress and its pull-back in the reference configuration of the skeleton reduce to vectors: and ; which implies , with an obvious definition of the third order tensor . On the other hand the fluid stress is split into a spherical and a deviatoric part, say , where this last does not a-priori vanish.

Consider the first three terms of equation (10), the following chains of equalities hold true:


once considered that for every vector field one has . Identities (11) imply the strain working (10) to be rewritten in terms of the pull-back of stresses in the reference configuration as:


where is the -th component of the pull-back of the bulk forces acting on the fluid, in the reference configuration of the skeleton.

3 Thermodynamics

Thermodynamics of porous media has been summarized by Coussy (2004, 2010) for two or more monophasic superimposed interacting continua, say the solid skeleton and the fluids saturating the porous space. In that case, the specific internal energies of the fluids, which fill the pores, 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. This contribution to the energy accounts for the micro-scale adhesion properties among the constituents, coherently with the Young-Dupré equation, which at the micro-scale states the equilibrium of the triple line keeping in contact the solid wall (of the pore), with the wetting and the non-wetting phase, see e.g. de Gennes (1985). Considering a geometrically simple configuration of the porous space and a rough upscaling rule, as that already mentioned in the Introduction, which assumes the variation of the area of interfaces proportional to the volume occupied by the phases, a macroscopic interfacial energy, per unit volume, can be explicitly separated from the energy of the skeleton and regarded as a function of the saturation ratio only, see Coussy et al. (2004). Afterwards the micro-scale surface tension can be related to the derivative of this energy, with respect to , which is interpreted a-posteriori as the so-called macro-scale capillary pressure, see e.g. Bear (1972); van Genuchten (1980) and Coussy (2004). The macro-scale capillary pressure, defined in this way, accounts at the same time for two distinct features, say the effect of surface tension between the liquid and the gas, stored within any possible reservoir, and the effect of the retention characteristics of the porous material.

Several criticisms have been moved to this definition of the macroscopic capillary pressure; in particular Gray and Hassanizadeh Gray and Hassanizadeh (1989); Hassanizadeh and Gray (1990, 1993) demonstrated that the spatial distribution of interfaces within a multi-phase system is fundamental in order to characterize the intrinsic state of the system.

Here a novel approach is adopted in order to incorporate, into the macroscopic constitutive prescription, the role of the interface between the two fluids stored in the pores, describing the fluid mixture as a non-uniform diphasic fluid in the spirit of Cahn and Hilliard (1958). A phase-field model with possible diffuse interface is considered, which is analogous to that one of (more than) two immiscible components incompressible flows, see Lowengrub and Truskinovsky (1998), Boyer and Lapuerta (2006), or Kim (2012) for a general review. As already noticed in § 1 and following Cueto-Felgueroso and Juanes (2012, 2014) the phase field here is the degree of saturation .

The internal energy of the fluid is given by


where the term is a double-well potential depending on the specific density and parametrized by the specific entropy . The non-local energy penalizes the formation of interfaces and provides a regularization of the non-convex energy contribution. As usual the state equation of the fluid defines a relation between conjugate variables; in particular the so-called fluid thermodynamic pressure , and the fluid chemical potential can be defined in terms of the the density :


in equation (14) is the absolute temperature, which is conjugate to the specific entropy by . The double-well volumetric energy resembles that of Van der Waals’ model and allows for describing coexistence of the immiscible phases, even if no mass exchange is considered. No information on the domains occupied by the phases in the current configuration is supplied by the model if this contribution to the fluid energy is the only non-vanishing one. According to Maxwell’s rule, an affine term, with respect to , can be added to the Van der Waals-like energy to make arise two isopotential phases, characterized by the same (vanishing) value of the chemical potential , see Figure 0(a) and Coussy (2010).

(a) Van der Waals’ like potential
(b) Chemical potential
Figure 1: The liquid and the gaseous phase coexist at equilibrium as they are isopotential minima of the fluid free energy

The non-local contribution allows for governing the coarsening of the domains occupied by the fluid phases and pattern formation. Indeed it provides just a correction of the bulk energy, when considering the so-called flat interface limit, see Cahn and Hilliard (1958), whilst it constitutes the main part of the energy of the non-uniform fluid, when describing phenomena in which a clear scale separation can not be assumed, as in topological transitions in multi-phase fluids Lowengrub and Truskinovsky (1998) or gravity driven fingering through porous media Cueto-Felgueroso and Juanes (2009a).

Because of the constitutive law (13), the following prescriptions on the fluid stress and hyper-stress hold true


being the identity tensor and the dimension of the Euclidean space . The thermodynamic pressure and the chemical potential given in (14) are just a part of these constitutive prescriptions. represents the deviatoric stress acting on the fluid, which is non-vanishing because of the gradient contribution to energy.

As already noticed incompressibility of the liquid phase means that the variations of the density are univocally determined by the variations of the degree of saturation, so that equation (15) can be rephrased in terms of the saturation ratio , see e.g. Cueto-Felgueroso and Juanes (2012). In particular the double well shape of the energy potential can be prescribed assuming the fluid energy per unit volume of the porous medium to be a kind of Duffing potential


whilst the non-local term is typically assumed quadratic in the gradient of :


Here is the surface tension between the non-wetting and the wetting phase, whilst is the characteristic size of the channel through which the fluid can pass, see e.g. Boyer and Lapuerta (2006). In Figure 1 the energy per unit volume and the corresponding chemical potential are plotted. The energy of the non-uniform fluid could also be prescribed in terms of a different phase field as for instance the mass concentration of the liquid without its double well feature being modified. With an abuse of notation from now on will indicate the derivative of with respect to rather than to , as indicated in equation (14).

As saturation can not overwhelm the limit , a suitable constraint should be stated when formulating the minimization of the fluid energy potential with respect to , say which corresponds to , with a slack variable used to transform the inequality constraint into equality. This yields the introduction of a Lagrangian multiplier, in the functional to be made stationary, representing the reactive chemical potential, , which allows to account for the transition from partially-saturated to fully-saturated conditions. In partially-saturated states implies and consequently ; on the other hand when saturation is attained vanishes and the reactive chemical potential plays the role of the liquid pressure in standard saturated poromechanics. Looking at equation (17) it is worth to notice that in the case of saturation the non-local contribution to energy still does not vanish, but reduces to a function of the porosity gradient.

The effect of confinement of the air-water mixture into the porous space, say the retention characteristics of the porous material, will be accounted for by an additional energy . Within the considered framework, no a-priori constitutive characterization is assumed on this functional, but suitable restrictions on it are deduced from the first and the second principle of thermodynamics, summarized in the generalized form of the Clausius-Duhem (dissipation) inequality.

In the following paragraphs the explicit form of the dissipation relative to the partially saturated porous medium is stated, taking in due account the expression of the strain working given by equation (12). The corresponding restrictions on the constitutive law of the porous skeleton and the local diffusion fluid mass flux are deduced.

3.1 The first principle of thermodynamics

Modifying the form of the energy equation, stated by (Coussy, 2004, equation (3.14)) for a standard porous continuum, so as to separate the interaction energy between the solid and the non-uniform fluid implies the first principle of thermodynamics to read as


where indicates the time derivative following the motion of the particle . In a similar way to the fluid, the energy of the solid is defined in terms of the corresponding intrinsic energy as . Moreover the rate of change of the coupling energy is associated partly to the motion of the solid grains, partly to the motion of the molecules of the fluid, say the liquid and the gas. accounts for the rate of heat externally supplied and can be prescribed in terms of a surface rate of heat due to conduction,


finally is the strain working, given by equation (12). The pull-back of the energy balance (18) into the reference configuration of the skeleton yields


where is the overall Lagrangian density of internal energy; indicates the time derivative following the motion of the solid particle fixed in the reference configuration of the skeleton. Equations (12) and (20) imply the following Lagrangian local energy equation to hold true:


3.2 The second principle of thermodynamics & the Causius-Duhem inequality

In the same way as the first principle, the second principle of thermodynamics (entropy balance) for the porous continuum is stated in the form:


where stands for the specific entropy of the -th constituent and again indicates the surface rate of heat due to conduction. The pull-back of the entropy balance in the reference configuration of the skeleton reads therefore:


being the overall Lagrangian entropy. Using the Legendre transform ; the local form of equation (23) can be written as:


which using equation (21) delivers the Clausius-Duhem inequality for a porous material within the framework of gradient poromechanics, see Sciarra et al. (2007) for similar results:


3.3 Dissipation

A characterization of the dissipation relative to the solid and the non-uniform fluid is obtained, starting from the Clausius-Duhem inequality (25), separating, within the overall free energy and the overall entropy , the contributions of the bulk fluid, say and , from the residual terms relative to the energy and the entropy of the bulk skeleton and the interfaces, say and :


Remind that is the Legendre transformation of , such that . This approach stems from Biot (1972), who defined the wetted solid as a system composed of the solid skeleton and a thin layer of fluid attached to the internal walls of the pore. The aim is proving that is a state function which prescribes the constitutive behavior of the skeleton in thermo-poroelasticity, say in the case when just reversible processes are involved by the deformation of the porous continuum. In this case the residual part of the dissipation will depend just on the Lagrangian filtration vector and the gradient of the absolute temperature.

The fluid mass conservation (2) together with the constitutive law relative to the fluid stress and hyper-stress (15), allow to rephrase the Clausius-Duhem inequality (25) in the following form:


The coefficients of , , and reported in equation (27) are explicitly deduced from equation (25) in Appendix A. Following Coussy (2004) the dissipation (27) can be split into three terms, one related to the solid skeleton, one to the non-uniform fluid and the last one referred to thermal effects. These three contributions are separately assumed non-negative:


As previously remarked, equation (28) states that within thermo-poroelasticity the Helmholtz free energy of the skeleton is a state function:


which depends not only on the strain of the skeleton, the Lagrangian porosity and the degree of saturation, as in standard unsaturated poromechanics, but also on the gradient of strain, , and the gradient of Lagrangian water content, . It is worth to notice that equation (31) could be also formulated considering a free energy function depending on the gradient of strain, , and the gradient of water content, , in addition to the dependence on , so modifying the definition of conjugate variables.

As the fluid dissipation (29) does not vanish except at equilibrium, the filtration force must be related to the filtration vector itself, in a way to fulfill the dissipation inequality during evolution. Finally equation (30) recalls the well-known result that heat flows through materials along the direction of the negative gradient of temperature (from higher to lower).

Now, assuming isothermal conditions, the poro-elastic constitutive equations of the solid skeleton and the generalized Darcy law relative to the fluid are deduced and compared with those ones typically adopted in modeling the behavior of partially saturated porous media.

4 Poroelastic constitutive relations

According with equation (28) and its corollary remarks, the poroelastic constitutive relations of the solid skeleton must fulfill the following restriction:


if no frozen contribution to the free energy has been taken into account, see e.g. Collins (2005) and Coussy et al. (2010). Within this framework one gets the constitutive relations for the overall stresses and ,


the constitutive prescription for the coupling energy term ,


and the constitutive constraints to be satisfied by the skeleton free energy once the thermodynamic pressure and the fluid hyper-stress have been assigned according with equation (16), say


When the gradient contribution to the fluid energy vanishes, equation (34) reduces to the prescription of the macro-scale capillary pressure, postulated in classical unsaturated poromechanics, see Coussy (2004), say , if one assumes


As a consequence equations (34) and (35) can be rephrased as the generalized constitutive characterizations of the macro-scale capillary pressure and the fluid thermodynamic pressure , which account for the spatial distribution of interfaces, within the porous network, considering the additional contribution of the fluid hyper-stress:


Equation (33) is standard, within the framework of gradient theories of continuum mechanics; on the other hand equations (37) deserve a deeper discussion in order to understand in which sense they can be considered an enhanced constitutive prescription of the macro-scale capillary pressure and the thermodynamic pressure of the fluid. Developing a micro-scale analysis this model will be compared, in § 6, with the micro-structured one which accounts for the dependence of the capillary pressure on the so-called interfacial area between the non-wetting and the wetting phase, see e.g. Hassanizadeh and Gray (1990, 1993); Joekar-Niasar et al. (2010); Reeves and Celia (1996).

4.1 Generalized effective stresses

Some additional remarks can be deduced from the thermodynamical restrictions stated by equation (28) if the solid grains, which constitute the matrix of the porous medium, undergo negligible volume change. In this case the Lagrangian porosity is definitely prescribed in terms of volume changes of the porous elementary volume, say: , being the reference value of porosity, see e.g. (Coussy, 2004, Chapter 3). The solid dissipation can therefore be rewritten and the generalization of Bishop’s effective stress naturally arises, within the framework of the considered gradient model. Assuming small strains, but allowing the porous skeleton to suffer strain gradients of order , with respect to the considered small perturbation parameter, the simplified expression of reads as follows:


Here and are the generalized Bishop stress and hyper-stress defined by


As a consequence the corresponding free energy of the solid skeleton can be regarded as a state function of strain, strain gradient, saturation degree and gradient of liquid content, only. The general form of equations (39), when no restriction to small strains is assumed, is reported in Appendix B.

If one aims at investigating non-reversible processes affecting the behavior of the solid skeleton, different dissipative mechanisms should be analyzed, which could concern interactions among the grains of the solid skeleton, or between the grains and the fluid. In the first case inelastic constitutive relations for the solid stress and (eventually) hyper-stress should be introduced, in the second one hysteresis during wetting-drying cycles should modify the prescriptions of the retention properties, see e.g. Hassanizadeh et al. (2002). It is not the goal of this paper to address modeling inelastic processes, however this can be done specifying the way in which dissipation depends on a suitable set of internal variables.

4.2 Generalized Darcy law

As usual in poromechanics the dissipation inequality (29) relative to the fluid can be satisfied requiring to be a quadratic function of , which means: This implies the coefficient of in equation (29) to be constitutively constrained to


which reads as a generalization of classical Darcy’s law. is the inverse of permeability, its dimension is therefore . According to equations (14) and (37) equation (40) can be reformulated as follows


where is the specific free energy of the fluid per unit volume of the pores and therefore coincides with the double-well energy potential , possibly describing coexistence of the liquid and the gaseous phases, see equation (16). For the sake of simplicity, the capillary energy has been introduced, following Coussy et al. (2004), so that .

According with classical arguments proposed within the framework of diffuse interface models in fluid mechanics, see e.g. Jacqmin (2000); Jamet et al. (2001); Boyer et al. (2010), the first two terms in equation (41) correspond to the so-called generalized chemical potential of the non-uniform fluid defined as


the whole quantity can therefore be interpreted as the generalized chemical potential of the pore-fluid, say of the non-uniform fluid within the pore network. This last accounts, on one hand, for surface tension effects which are typical of a diphasic fluid on the other one for the retention properties of the skeleton, due to its texture, by means of .

Apparently the role of the capillary energy is that of modifying the double-well potential which prescribes the free energy of the fluid, in order to account for the presence of a confining surface densely distributed within the porous continuum. The new free energy , which can be called effective energy of the pore-fluid, has not the same minima as , as they are shifted inward the interval from below or from above whether the solid skeleton is gas or liquid wet, see e.g. Papatzacos (2002); Papatzacos and Skjoeveland (2004). In § 5 the effects of combining these two energy contributions are discussed, considering the retention properties of different soils. A similar behavior has been discussed by Cueto-Felgueroso and Juanes (2014), modeling two-phase flow in a Hele-Shaw cell, introducing a suitable symmetry-breaking function.

5 Constitutive 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 a proper prescription of the retention curve, which means prescribing the capillary pressure , or the capillary energy , as a function of . As already noticed in § 3 however, the air-water mixture, which saturates the pore space, is regarded here 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 therefore 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 algebraic relation between the saturation degree and the capillary pressure or the chemical potential of the fluid can a-priori be stated, without solving, at least at equilibrium, the following generalized Richards equation:


Here the permeability of the porous medium has been assumed isotropic: , being the so-called saturated permeability and the relative permeability of the wetting phase, see e.g. Coussy (2004, 2010). As usual in unsaturated poromechanics equation (43) which provides the distribution in space and time of the saturation degree is deduced merging the mass conservation law (2) and the Darcy law, which in this case takes the form given by equation (41).

At stationary conditions, say when , equation (43) 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 with respect to , which allows for describing the confining effect on the non-uniform fluid, due to the presence of the porous skeleton. Equation (43) looks also similar to that one stated in (Cueto-Felgueroso and Juanes, 2008, equation (2)) and (Cueto-Felgueroso and Juanes, 2009a, equation (13)), who adopted a phase field model for describing gravity fingers and saturation overshoot, during water infiltration through a non-deformable soil. In this case the difference resides in the account for the chemical potential in (43), allowing for possible coexistence between the phases of the non-uniform fluid. Moreover equation (43) is here definitely coupled with the balance equation of the overall continuum.

Following the same procedure which induced the definition of the generalized chemical potential (42) when simulating multi-phase flows, see e.g. Boyer et al. (2010), equation (43) can be rephrased in terms of the functional derivative of the free energy of the pore-fluid


with respect to , as follows:


where the generalized effective chemical potential of the pore-fluid has been introduced. The energy of the skeleton depends on the gradient of the liquid content, so that the gradient contribution to the constitutive law of the non-uniform fluid is not modified, with respect to equation (15). For the sake of completeness the effective chemical potential of the pore-fluid, accounting only for the derivative of the effective energy is defined as . As already noticed is different from µ, because of the term involving the derivative of the capillary energy , but it also differs from the flow potential introduced by (Cueto-Felgueroso and Juanes, 2009a, equation (12)), because of the term involving the derivative of the fluid energy .

The role of and in the characterization of the spatial distribution of , , and deserves a deeper analysis in order to underline the main differences between the proposed modified Richards equation (45) and that one stated by Cueto-Felgueroso and Juanes (2009a) for describing gravity fingering in soils. In Figure 2 the chemical potential of the pure fluid prescribed, according with equation (16), as


and the derivative of the capillary energy, given by a van Genuchten-like curve


are depicted (gray dotted and gray dashed lines, respectively), considering different soil textures and therefore different choices for the relative partition of sand, clay and silt in the soil. The profile of the negative effective chemical potential of the pore-fluid with respect to is also drawn in each panel of Figure 2 (solid lines). Both the chemical potential and the derivative of th capillary energy are expressed in head units, via the two fields and , adopting the classical Leverett (1940) scaling, which prescribes the characteristic length in terms of the intrinsic permeability (of the soil): . The different values of the parameters , and (), which characterize the retention curve (47) for the considered soil textures, are listed in the panels.

(a) Sand
(b) Loamy sand
(c) Sandy loam
(d) Loam
(e) Silt
(f) Silt loam
(g) Sandy clay Loam
(h) Clay loam
(i) Silty clay loam
(j) Silty clay
(k) Clay
Figure 2: Heads relative to the chemical potential of the pure fluid (dotted gray line), the derivative of the capillary energy (dashed gray line) and the negative effective chemical potential (solid black line). The values of the parameters which characterize the retention curve (47) are listed in each panel, together with the intrinsic permeability and the referential porosity , relative to the portion of sand, clay and silt which constitute the soil. The red spots in panels (a), (b), (c), (g) and (i) highlight the real zeros of , different from , which corresponds to the highest relative minimum and the saddle point of the effective free energy