Fingering and strain localization in porous media during imbibition processes

Fingered infiltration of a wetting fluid through a porous network is a widely studied subject in the field of fluid mechanics. However, the effect of this heterogeneous percolation on the response of granular materials, in particular fine-grained soils, is a poorly investigated and badly understood topic which deserves deep analysis, considering, among others, possible applications in soil remediation and underground energy storage. This paper presents a first application of a new formulation of unsaturated poromechanics based on a phase field approach that allows to characterize on the one hand the occurrence of fingering hydraulic instabilities and on the other one to capture their effects on the irreversible, and possible unstable, deformation of the solid skeleton

One of the most relevant characteristics of this process is the tendency of the interface between the injected wetting fluid and the receding non-wetting one to destabilize in the form of a fingered front, in other words to exhibit preferential paths, when a driving flux forces the fluid-fluid interface through the porous network. This phenomenon has been extensively studied in the past from the experimental point of view considering first of all the interaction between the two fluids in a thin gap separating two flat plates (Hele-Shaw cell), without any obstacle, see Lenormand 4 and Guo. 5 More recently important results have been obtained which clarify the role of wettability, summarized by the contact angle, in the transition from stable to fingered percolation, 6 and in enhancing/contrasting the efficient displacement of the defending fluid by the invading one. 7 The microscale mechanisms responsible of these two last counter-current effects have also been identified. The presence of gravity typically enhances the destabilization of the front when the injected fluid is more dense of the defending one, see for example, DiCarlo et al. 8 From the theoretical point of view the seminal paper of Saffman and Taylor 9 constitutes a milestone in the formulation of the problem followed by the study of Chuoke et al. 10 and Wang et al., 11 all based on the characterization of a linear instability criterion for the moving fluid-fluid interface. More recently an effort has been done to embed this study in the framework of the mechanics of porous media. On the one hand it was proven that the Richards equation, commonly adopted to describe bi-phasic flow through porous media, is unconditionally stable with respect to transversal perturbation and consequently unable to reproduce the fingering phenomenon. 12 On the other hand gradient and phase field models have been introduced to obtain a regularized representation of the fluid-fluid interface, studying its stability and modeling fingering occurrence. [13][14][15] The present study stems from these results in mechanical modeling and aims to face a quite different aspect of the imbibition process which has been poorly investigated in the past, say the action of heterogeneous percolation, in particular fingering, on the local rearrangement of grains in a fine-grained material. The idea is not just to consider the porous network as a sequence of rigid obstacles to the fluid flux but also to account, in the framework of continuum poromechanics and therefore adopting a purely macroscopic point of view, for the effect of this heterogeneous flux on the deformation of the solid skeleton and potential strain localization. As mentioned by several authors, see among others Li and Vanapalli 16 and recently Liu and Santamarina, 17 not just drainage but also imbibition processes can cause fabric changes in fine-grained sediments, which could be associated to contractant behavior (capillary collapse), as proven by Bruchon et al., 18 but also to swelling, because of the alteration of the yield surface and consequently of the reversibility domain. 19 To this purpose the present study develops a poromechanical model based on a phase field approach to partial saturation, 20 in order to overcome the above mentioned weakness of the formulation based on Richards' equation, endowed with an elasto-plastic model taking in due account the effect of saturation on the reduction of strength in a similar way as done by Tamagnini 21 and Rotisciani et al. 22 In particular the two above mentioned features of a fine-grained soil (a loamy sand), say its capability to contract or swell under hydraulic imbibition, are discussed as a consequence of the initial state of stress and the rate of change of the yield surface due to saturation.
The paper is organized as follows: in Section 2 a brief presentation of the phase-field approach to partial saturation is provided together with the implications of the extended Clausius-Duhem inequality in terms of the thermodynamical restrictions on both the behavior of the solid phase and the mixture of a wetting and a non-wetting fluid saturating the porous network. In Section 3, the governing equations of the problem are summarized and in Section 4 the corresponding weak formulation is derived together with the Algorithm 1 used to implement this new formulation within the Matlab Finite Element code written by Bonnet and Frangi. 23 Section 5 is devoted to the analysis of the pure hydraulic problem in order to get a preliminary estimate of the instability conditions that trigger fingering formation. Finally Section 6 presents the results of the coupled problem, in particular two cases corresponding to a dilatant and a contractant behavior of the same soil are discussed. Section 7 summarizes conclusions and perspectives of this study.

KINEMATICS AND THERMODYNAMICS
We consider a porous medium constituted by a deformable solid skeleton whose pores are saturated by a mixture of two immiscible fluids phases, a wetting liquid (water) and a non-wetting gas (air), this last directly in contact with the atmosphere. The liquid phase is assumed incompressible ( w = const) and the gaseous phase passive, say of infinite mobility. The gas density can therefore be neglected with respect to that of the liquid phase ( g ≃ 0). Rather than the classical approach to unsaturated poromechanics, 24 an alternative modeling scheme 20 is here adopted, consisting in considering the gas-liquid mixture filling the porous network as a non-uniform fluid in the sense of Cahn-Hilliard. 25 To this aim an order parameter, the phase field, is introduced which is uniform in the bulk phases and varies continuously through the interface between them. From now on the volume of the liquid phase with respect to the volume of the pores, say the degree of saturation S r will play the role of the phase field, being equal to 1 in the liquid and 0 in the gaseous phase. Being the gas passive, the non-uniform fluid density is f = w S r .
In a similar way as in the standard formulation of continuum poromechanics, 24 the pull-back of the mass balance of the non-uniform fluid in the reference configuration of the solid skeleton can be written, stating the small strains (small displacement gradients) assumption to hold true, as follows: where m f is the mass of the mixture per unit reference volume, M is the mass flow vector, the Lagrangian porosity, V f is the velocity of the fluid and V s that of the solid, this last being obviously the time derivative of a suitable displacement u.
It is worth to underline that in this case the fluid under consideration is the above mentioned mixture of the two phases. Following the seminal works of Cahn-Hilliard the distinction between the two phases is operated via an adapted form of the Helmoltz free energy of the fluid mixture, which is given by where Ψ f is a double well potential with two isopotential local minima which correspond to the equilibrium states of the mixture, say the above mentioned phases. The gradient, non-local, energy contribution Ψ nl incorporates a penalization associated to the diffuse interface formation and allows to convexify the problem. The coefficient C k , which controls the strength of the gradient-energy terms, is a higher order stiffness which has the same physical dimension as an energy per unit volume multiplied by a squared length . This last characterizes the thickness of the diffuse interface between the two phases. We underline that the gradient term does not involve just the saturation but the water content S r so as to account for the gradient of water distributed within a representative volume rather than just in its pores. The thermodynamic pressure P and the chemical potential of the fluid can be defined via the state equations in terms of the saturation degree S r as From now on we adopt the following simple functional form for the fluid bulk energy which assures the two isopotential minima of Ψ f associated to the gaseous and the liquid phase to be at S r = 0 and S r = 1, respectively. As depicted in Figure 1 the two minima of Ψ f correspond to the zeros of its non-monotonic first derivative . The physical meaning of will be clarified once the restrictions on the constitutive laws imposed by the first and the second principle of thermodynamics will be stated, see Section 2.4. In Equation (4), lg is the liquid-gas surface tension and R a characteristic length describing the effective size of the pore space through which the fluid flows; in the following the Leverett 26 estimate of R is considered, say: R = √ 0 , 0 being the initial porosity of the material and its intrinsic permeability. In a similar way as modeling of multi-phase flow in fluid mechanics, see for example, Kim 27 and references therein, this approach, whose extended formulation is provided by Sciarra, 20 allows to capture pattern formation in fluid flow through porous media, as for instance pinching-off, coalescence and fingering of a wetting fluid displacing a non-wetting one initially saturating the porous network.
As discussed by Sciarra 20 in order to account for the transition from partially to fully saturated conditions, the inequality S r ≤ 1 should be also taken into account. This can be done via a slack variable transforming the inequality constraint into an equality. The Lagrangian multiplier corresponding to the constraint is in this case the reactive chemical potential r , which will be zero in partially saturated conditions and different from zero when saturation is achieved. We refer to Sciarra 20 for more details. For the sake of simplicity all the deductions which are presented in the following are however only valid in the regime of partial saturation. A straightforward extension of these results can be obtained to incorporate in the model the above mentioned slack variable and the corresponding reactive chemical potential. Let us now move to thermodynamics and constitutive laws of the porous medium. Assuming the power of external forces to be a linear functional of the velocity field parametrized by bulk forces, tractions and double forces, 28 acting on the solid and the fluid mixture, and requiring the overall balance of momentum to hold true, allows to deduce an extended version of the Clausius-Duhem inequality from the first and the second principle of thermodynamics. We refer to Coussy 24 for the general framework of thermodynamics of porous media and to Sciarra 20 for the specific formulation valid in the case of second gradient continua, see Equation (27) therein. It is not the purpose of this paper to deduce again this result, the interested reader could refer to the above cited reference for further details. In order to further simplify the formulation, the hyper-stress of the solid is assumed to be the negative of that of the fluid so as to obtain a vanishing overall hyper-stress. Under isothermal conditions the Clausius-Duhem inequality can be finally formulated as: Here is the total stress of the overall porous medium, required to verify the overall momentum balance equation div + g = 0, being the total density of the porous medium including the solid skeleton and the fluid inside the pores; P c is the capillary pressure, is the fluid hyper-stress vector and b f the bulk force exerted on the fluid. Ψ s is the free energy of the porous solid defined as the difference of the free energy of the overall porous continuum and the bulk energy of the non-uniform fluid multiplied by the Lagrangian porosity, Ψ f . As the double well potential Ψ f accounts only for the energy of the fluid mixture, all the contributions relative to the solid-fluid and the fluid-fluid interfaces are stored into Ψ s . In the following an additive decomposition of the solid energy Ψ s will be proposed.
According to Coussy, 24 the dissipation inequality given by Equation (5) can be split into two separate non-negative contributions: Φ s relative to the solid and Φ f relative to the non-uniform fluid defined as: In the following the implications of the two contributions to the dissipation inequality are separately discussed.

Prescriptions on the constitutive law of the solid, with interfaces
We consider first of all the poro-elastic conditions, which correspond to the case when Equation (6) holds true as an equality. As usually done in rational thermodynamics, 29 assuming the dissipation to vanish allows to characterize the free energy as a state function, in this case of the state variables e , e , S r and ∇( S r ). e indicates the purely reversible part of the strain in the additive decomposition = e + p , valid within the assumed framework of small strains, while e corresponds to the reversible contribution to the current value of the Lagrangian porosity, also provided by an additive decomposition. In other words, the state of the porous skeleton is characterized at any place x and at any time t by the values of each of these variables. Hence, the skeleton Helmholtz free energy is assumed to have the following form: which depends not only on strain, porosity and saturation degree, as in standard unsaturated poromechanics, but also on the gradient of the wetting fluid content S r . Using the chain rule to develop the time derivative of Ψ s and replacing the obtained result into Equation (6) yields the following form for the skeleton dissipation: so that the skeleton state equations can be given by: According to the previous remarks on the free energy of the porous skeleton and adopting a similar approach as in classical unsaturated poromechanics to account for the retention properties of the solid, the following decomposition of Ψ s is proposed: where s represents the energy of the saturated porous solid and U is the so-called capillary energy, which only depends on the saturation degree and accounts for the energy stored within the solid-fluid interfaces. As already mentioned the energy of the fluid-fluid interfaces, intrinsic to the non-uniform fluid, is counted into Ψ nl . It is worth to notice that decomposition (11) implies that strain and porosity variation do not affect the capillary energy U and consequently the water retention curve. This choice assumed for the sake of simplicity must obviously be relaxed within the framework of a more general formulation, in particular when the soil volume significantly reduces. 30,31 Replacing Equations (11) into (6) allows to rephrase the third of the state Equation (10) in the usual form −P c = dU dS r . The classical van Genuchten 32 expression of the capillary pressure is assumed: where the dimensionless parameter m and the inverse of a characteristic length assign the retention properties of the porous medium, S res r being the residual saturation degree of the wetting fluid. The capillary energy is the integral of P c (S r ) between S r and 1.
Once moving from poro-elasticity to poro-elasto-plasticity we still consider the Helmoltz free energy Ψ s to be a function of the elastic strain, the elastic porosity, the saturation and, in this case, the gradient of the wetting fluid content; additional dependency on suitable internal (hardening) variables j is also considered: say Ψ s = Ψ s ( e , e , S r , ∇( S r ); j ). It is worth to underline that no hypothesis is formulated at this level concerning the tensorial order of these variables. Moreover considering that elastic deformations occur at a time scale much lower than that which characterizes dissipative processes, Equation (10) continue to hold true, which implies the solid dissipation Φ s to reduce to As usual in plasticity, this last is verified with a proper choice of the potential g providing the plastic strain rate (and the plastic porosity rate) as a function of the state of stress. 24 The plastic model adopted in the present study will be discussed in Section 2.3.

The effective stress
In order to further simplify the formulation, we refer here to the simplest deduction of the effective stress tensor in partially saturated conditions as reported for example, Coussy et al. 33 Assuming the solid grains, forming the matrix, to be incompressible and in the absence of any occluded porosity, the volumetric strain of the porous medium results from the porosity variation: Equation (14) can be verified, considering e = e − 0 anḋp =̇p.
In the poro-elastic case Equation (5) therefore reduces to with the effective stress ′ defined as the total stress plus the so-called equivalent pore pressure : and the energy per unit volume s just dependent on e . Apparently the state equations are therefore simplified into Using the same arguments as in the previous section, the reduced form (13) of the dissipation inequality in the poro-elasto-plastic case can be further simplified into

The elasto-plastic constitutive model
The constitutive relation for the effective stress is provided via an elasto-plastic constitutive model where, for the sake of simplicity, the elastic part is kept as the classical Hooke law, while the plastic one is assumed to be given by the modified Cam-Clay model (MCC). Let p ′ = − tr ′ 3 be the effective pressure and q ′ = √ 3 2 s ′ ∶ s ′ the deviator stress defined as the second invariant of the deviatoric (effective) stress s ′ = ′ + p ′ 1. The yield function f has the form where p c is the so-called preconsolidation pressure, say the highest mean effective stress ever experienced by the soil and M indicates the slope of the critical state line, which can be related to the friction angle via M = 6 sin ∕(3 − sin ). The preconsolidation pressure is, in the current formulation, the only non-vanishing among the previously cited hardening variables . A detailed presentation of the MCC model can be found in Nova. 34 The MCC being an associated plastic model, it is therefore characterized by the associated flow rulė prescribing the evolution of the plastic strain in terms of the derivative of the yield function with respect to the effective stress.Λ is the corresponding plastic multiplier. The presence of the hardening variable p c necessitates of a hardening law, which in the MCC is given by: wherėp v indicates the volumetric plastic strain rate, e is the void ratio, the compression index and the swelling index.
In the case of unsaturated soils it is of interest to include in the formulation of the plastic model, and more specifically in that of the hardening law, the effect of the progressive saturation of the porous network. Following a similar approach as that developed by Nova et al. 35 to describe plastic strains in bonded geomaterials, caused by bonding failure induced by chemical or mechanical degradation, a similar approach has been proposed 21 and exploited 22 to describe plastic deformations in geomaterials caused by failure of capillary bridges induced by saturation. In this framework the variation of the preconsolidation pressure depends not only on the volumetric plastic strain but also on suction or saturation degree. In particular a double hardening law 21 is introduced for the preconsolidation pressure which separately accounts for a parametrization ofṗ c by the plastic strain rate as well as by the variation of the saturation degree: whereṗ c(sat) is given by the original MCC hardening law, see Equation (21), andṗ c(unsat) describes the variation of the preconsolidation pressure with respect to the variation of the saturation degree: where is a new constitutive parameter to be calibrated. Equation (23) implies that the preconsolidation pressure decreases with increasing saturation, which causes the shrinkage of the yield surface and may induce plastification of the solid skeleton even at constant effective stress. As previously mentioned the elastic behaviour of the solid phase is assumed linear; and the isotropic material symmetry is also chosen. Accordingly the effective stress rate is given bẏ where K and G are the skeleton bulk and shear module, while J and K are the tensor product of two second order identity tensors and the fourth order identity tensor minus one third of the previous tensor J, respectively. The plastic strain rate, given by Equation (20) is determined calculating the positive plastic multiplier from the consistency condition, that is, requiring that a state of stress which at time t belongs to the yield surface f = 0 remains on it at time t + dt:

The generalized Darcy law
Positive definiteness of the fluid dissipation Φ f can be satisfied assuming this last to be a quadratic form of the Lagrangian mass flow vector M, so that: where is the already mentioned intrinsic permeability, the viscosity of the wetting fluid and k(S r ) a function that accounts for non-uniform mobility of the fluid mixture within the pore network. From now on we assume that it coincides with the relative permeability of the wetting phase. Accordingly, the van Genuchten or the simpler Leverett form of k(S r ) can be adopted. 36 Comparison of Equations (7) and (27), and taking in due account the state equations (17), yield the following generalized form of the isotropic Darcy law: which with respect to the standard one, that in partially saturated conditions just depends on the gradient of the capillary pressure, incorporates the contributions arising from the thermodynamic pressure P of the non-uniform fluid and the hyper-stress fluid vector . Using the definition of the thermodynamic pressure given by Equation (3) and the state equation (17) 3 implies Equation (28) to be rewritten as follows: It is worth to underline that in a similar way as discussed by Gurtin, 37 the argument of the ∇ operator in Equation (29) can be regarded as the variational derivative of Ψ mix , defined by Equation (2), plus U with respect to S r divided by , say . From now on the notion of pore-fluid energy is introduced as Ψ pf = Ψ f + U, so that, taking in due account the form of Ψ nl , a kind of effective chemical potential can be defined: which incorporates the so-called pore-fluid chemical potential, pf obtained as the derivative of the pore-fluid energy Ψ pf with respect to S r , and the divergence of the derivative of the non-local energy contribution, with respect to the gradient of the wetting fluid content. In the absence of bulk forces acting on the fluid, Equations (29) and (30) imply that the fluid mass flows from higher to lower effective chemical potentials. Apparently in the presence of conservative bulk forces an augmented chemical potential could be defined, as commonly done in the case of gravity. The pore-fluid energy Ψ pf describes the bulk energy of the wetting and the non-wetting phases stored into the porous network. Due to the form of the capillary energy, Ψ pf can possibly maintain a similar double well shape as Ψ f , however the two minima, if they exist, are no more isopotential. A double-tangent construction, equivalent to the well-known equal-area Maxwell rule can be adopted to identify the value of the pore-fluid chemical potential which re-equilibrates the two wells. From the physical point of view this is the same as saying that isopotential minima can exist when a suitable suction is applied at the boundary. Apparently the new isopotential minima will fall inside the (0, 1) interval, rather than being on its boundary, see Figure 2. When the double tangency condition is verified, and therefore its derivative with respect to S r is a non-monotonic function of the saturation, the current formulation is similar to the one based on a modified capillary pressure, introduced by DiCarlo et al. 38 which incorporates a hypo-diffusive term in addition to the capillary and gravity terms within the traditional Richards equation. The higher order stiffness C k is from now on prescribed in terms of the characteristic lengths −1 and , associated to capillary rise and fluid-fluid interface thickness, as C k = w ||g|| −1 2 , where ||g|| is the magnitude of the gravity acceleration g.

THE GOVERNING EQUATIONS
The equation governing the behavior of the fluid is obtained by substituting the generalized Darcy law Equation (29) into the Lagrangian mass balance of the fluid mixture Equation (1), which implies the following fourth order partial differential equation involving the current values of porosity and saturation degree: Even if Equation (31)  forms of instability of a gravity driven imbibition front propagating through a porous medium, even without any heterogeneity in the permeability or in the retention properties of the material. 14,15 A detailed discussion of instability conditions of a fluid-fluid interface with respect to transversal perturbation is reported in Section 5. According with the definition of the effective chemical potential provided in Equations (30), (31) can be rewritten as follows: Equation (32) has the same form as the well known Richards equation, 42 describing flow induced by gravity of a wetting fluid through an unsaturated porous medium, however in the present formulation the chemical potential eff , replacing the derivative of Richards capillary potential, is not directly prescribed via a state equation, function of the saturation degree and possibly of porosity, but it is the solution of the coupled boundary value problem defined by Equations (30) and (32), endowed with suitable boundary conditions. Apparently Equation (32) can be reduced to Richards equations assuming b f = w S r g, say the gravity force acting on the fluid, and neglecting the energy contribution of the non-uniform fluid Ψ f as well as the gradient, non-local, energy contribution Ψ nl . The proposed numerical implementation, that will be discussed in Section 4, is based on the above decomposition of Equations (31) into (30) and (32), considering therefore the degree of saturation as well as the chemical potential as the two primary unknowns of a system of two coupled second order partial differential equations. Equations (30) and (32) must be accompanied by suitable boundary and initial conditions. Adopting the usual distinction between essential (Dirichlet) and natural (Neumann) boundary conditions, the boundary of the integration domain could be differently subdivided when considering the boundary conditions relative to the above mentioned partial differential equations. In particular considering Equation (30), essential boundary conditions assign the value of the saturation degree S r while natural boundary conditions prescribe that of its normal derivative, which means providing the value of Young's contact angle. [43][44][45] Concerning Equation (32), essential boundary conditions assign the value of the effective chemical potential eff , natural boundary conditions its normal derivative. Taking in due account the form of Equation (32) this will imply providing the value of the fluid inflow. In Section 4 the weak formulation of the problem will allow to underline this fundamental distinction. To solve Equation (32) also the initial value of saturation must be prescribed.
In addition to Equations (30) and (32), with their own boundary and initial conditions, the momentum balance equation of the overall porous medium must be also considered, say where the total stress is given by Equation (16) coupled with the stress-strain constitutive relation of effective stress Equation (24), the flow rule Equation (20) and the consistency condition Equation (26). In this case prescribing essential boundary conditions will obviously mean assigning the value of displacement, while prescribing natural boundary conditions will obviously mean assigning the value of the traction applied on the overall porous medium. It is worth to underline that Equations (30), (32), and (33) constitute a fully coupled system of second order differential equations in the space variable and a first order differential equation in time.

INTEGRATION OF THE BOUNDARY VALUE PROBLEM
The coupled problem described by the mass balance Equation (31), endowed with the generalized Darcy law, Equation (29), the momentum balance Equation (33) and the constitutive law of the effective stress is here reformulated in a suitable weak form so that an adapted finite element method implementation strategy can be developed. Using the formulation of the fluid problem based on Equations (30) and (32) rather than the one directly based on Equation (31) implies, as already mentioned, the unknowns of the boundary value problem to be the two scalar fields S r and eff as well as the vector field u, representing the displacement of the porous medium; the plastic strain p which is also an unknown of the problem does not necessitate any boundary condition and is therefore treated as a local variable, as discussed in the following.
Let Ω be the partially saturated porous medium which constitutes the domain of the boundary value problem, and Γ its boundary. According to remarks concerning essential and natural boundary conditions for the balance of total moment and fluid mass stated in Section 3, Γ can be differently decomposed into two parts, Γ D and Γ N , where Dirichlet and Neumann boundary conditions are imposed respectively. Different decompositions of Γ can be considered when is chosen as the saturation degree S r , the effective chemical potential eff or the displacement u so that In Equations (34)-(36) T indicates the applied traction, q w the fluid inflow and q S r the so-called doubly normal double force. 28 The weak solution of the boundary value problem defined by Equations (30), (32), (33), (34)-(36) will be looked for in the functional space , R n being the vector space where displacement are observed. Test functions u, and S r , will be taken in the space H u0 × H 0 × H S r 0 , the 0 indicating that homogeneous Dirichlet boundary conditions are assumed. The weak formulation is obtained as usual multiplying Equation (32) by , Equation (30) by S r and Equation (33) by u and integrating over the domain Ω using the Green-Gauss Theorem to eliminate second-order derivatives. This leads to which must hold true for any u ∈ H u0 , ∈ H 0 and S r ∈ H S r 0 . Equation (37) represents the weak form of the momentum balance Equation (33) and can obviously be recognized as the principle of virtual working. In Equation (37) the test function u can therefore be interpreted as the virtual displacement, and consequently as the virtual strain. In Equation (38) the bulk force is assumed to be the weight of the wetting fluid as in Richards equation. It is worth to underline that Equations (37)- (38) must hold true at any time t. As usual when working with elasto-plastic problems the above mentioned weak formulation must be complemented by the equations that provide at any time step the current value of the effective stress as a function of the plastic strain. These are local equations which do not necessitate boundary conditions and are therefore just solved locally.
The spatial discretization is obtained by applying the usual Galerkin procedure with linear Lagrange finite elements. The unknowns u, eff and S r are therefore prescribed at any time t as a linear combination of piecewise linear functions weighted by the nodal values of the unknowns, say where N u , N , N s are the shape function matrices, say matrices assembling at the global level the contribution of the above mentioned piecewise linear functions that provide the approximate solution of the problem within each element, and u, , S r the vectors of nodal values. In a three dimensional space, the 3M-dimensional vector of nodal displacements, with M the number of nodes, is organized in the form In a similar way N and N s are line vectors of dimension M.
As a result of the spatial discretization, Equations (37)-(39) can be compactly written in the following form: where the l index indicates the nodal force working on the corresponding nodal virtual displacement in Equation (41), the nodal volume change of the liquid wetting phase working on the corresponding virtual chemical potential in Equation (42) and the nodal chemical potential working on the corresponding saturation in Equation (43).
The following definition of the quantities introduced in Equations (41)- (43) are assumed: Combining Equations (41), (42), and (43) the following coupled system of non-linear differential equations is obtained where the vector U and the matrices C(U) and K(U) and the vector N(U) are defined by Concerning time discretization, the implicit Euler scheme of the first order is used, which is unconditionally stable and commonly used in computational poromechanics literature. Following this method over a time step Δt = t n+1 − t n and using the approximation: Equation (44) can be rewritten at a time t n+1 as: The above system of non-linear coupled algebraic equations that allows to determine nodal displacements, effective chemical potential and saturation degrees will be solved adopting the Newton-Raphson method. As mentioned before the non-linear system (41)-(43) is coupled with the equations that allow to update, at any time step, and during the iteration process, the stress-strain constitutive relation. Being these equations local, they are required to hold true, for any finite element, just in the Gauss integration points used to obtain Equation (41) as the discretized form of Equation (37). A return mapping algorithm, 46,47 parametrized by the current approximation not only of the displacement solution but also of the saturation degree which controls the preconsolidation pressure is therefore implemented. A scheme of the overall integration algorithm is provided in Algorithm 1.

IMBIBITION INTO A RIGID POROUS SKELETON
In this section a brief analysis of the imbibition process of water into a relatively dry silt layer, driven by gravity is presented. The intention is to demonstrate the conditional stability of the imbibition front that results in fingering type instabilities. This is done, following Ommi et al. 14,15 by first resolving the traveling wave (TW) solutions of the imbibition boundary value problem (BVP), followed by a linear stability analysis against transverse perturbations. For the sake of simplicity the porous skeleton is assumed to be non-deformable ( = 0 = const). The sections that follow will be devoted to the analysis of the coupled hydro-mechanical problem where such an assumption is not done so that the influence of hydraulic instabilities on the mechanical response is revealed. Indeed, the current section serves as a precursor to the one following by providing an indication on the morphology of the fingering instability that can be expected for the given conditions.

One-dimensional TW-solutions
The imbibition fronts are assumed initially uniform along the transverse directions to direction of propagation and then are perturbed by a transverse harmonic field. These uniform solutions are of TW-type on an infinite domain of propagation with two uniform states of saturation degree connected by a diffused front that propagates with a constant speed, c, and thus are self-similar in nature. The rigorous proofs of existence and uniqueness of TW-solutions are out of scope of the current work. Instead, we resolve these solutions numerically by first posing the TW-problem corresponding to imbibition. This is done under the transformation, where = x − ct is the TW-coordinate moving with the speed of the front c. Introducing this transformation into the one-dimensional version of Equation (31) with b f = w S r g and g = ||g||e x gives, where (⋅) ′ is used to denote derivative of the function (⋅) with respect to s. The appropriate boundary conditions corresponding to imbibition fronts on an infinite domain are, where s + and s − represent the uniform saturation states on either side of the diffused front. s − is consistent with the injection flux at the boundary x = 0 m of the physical domain through the relation, Note that this is the same as the natural boundary condition derived in Equation (35) with ∇ eff vanishing, corresponding to the boundary condition of the TW-solution. Since we are posing the problem on an infinite domain, the speed of the front, c, is given by the Rankine-Hugoniot jump condition, To resolve the TW-solutions we use a second-order accurate central difference scheme and with uniform discretization on a sufficiently large finite domain ∈ [−L∕2, L∕2]. Further details regarding the above developments and numerical approximation can be found in Ommi et al. 14 Figure 3 shows these resolved TW-solutions for s + = 0.446, for different values of s − with soil hydraulic properties corresponding to loamy sand and L = 100 m. The values of the constitutive hydraulic model are reported in Table 1 of imbibition into dry sand by Ommi et al. 14 In essence, the two iso-potential minima, (S 1 c , S 2 c ), of the pore fluid energy, Ψ pf (S r ), determine the range within which the saturation degrees when used as boundary conditions cause overshoot (in the case of s − ) and undershoot (in the case of s + ) within the solution, respectively on the left and right-hand sides of the diffused front. For instance, in Figure 3, for s − < S 1 c an overshoot appears when the solution approaches the diffused front from the left. The small oscillation on the other hand when the solution leaves the diffused front towards s + is due to the complex nature of the stability properties of the associated equilibrium of the dynamical system governed by Equation (49). 14 In Figure 4 the function eff (S r ) is plotted for each TW-solution resolved earlier, along with the local part of the chemical potential pf (S r ). Since eff (S r ) has a non-local contribution, see Equation (30), it does not follow exactly the graph of pf (S r ). Instead, it follows a corrected path governed by the spatial gradient of the saturation degree, which is quite common in Cahn-Hilliard like fluid modeling.

Stability of solutions against transverse perturbations
In the work of Saffman and Taylor 9 the stability of a horizontal interface between two viscous fluids within a Hele-Shaw cell against transverse harmonic perturbations was studied. Stability of the said interface is understood as an exponential decay in time of the imposed perturbation for every possible wave number. Following this approach numerical study of the relationship between the maximal growth rate of the imposed perturbation and its corresponding wave number, the so-called dispersion relation, has been done in various works 12,13,15,48,49 for specific models intending to describe such fluid-fluid displacement in porous media. In this sense, the one-dimensional TW-solutions resolved earlier form the two-dimensional base solutions when extended uniformly along the transverse y-direction. These base solutions are represented s 0 ( ) which is by construction independent of y. Now, the solution of the Equation (31) in the moving coordinate system is assumed to be a perturbed one, that is composed of the base solution at the leading order and superposed perturbations of decreasing order, where represents the magnitude of the disturbance. Introducing such an expression into Equation (31) in the moving coordinate system results in a perturbed problem of order O( ) which governs the perturbation s 1 ( , y, t). Since the imposed perturbation is assumed harmonic and we intend to study its exponential growth/decay in time, the following form is assumed: wherek is characteristic wave number of the disturbance in y-direction,̂is the exponential growth factor in time andŝ( ) is the amplitude of the wave-like disturbance. Introducing Equation (54) into the perturbed problem for s 1 ( , y, t) results in a linear homogeneous ODE forŝ( ), where the spatially varying coefficients have the following expressions: and k 0 and pf 0 are respectively the functions k(s) and pf (s) evaluated at the base solution s 0 . For the order O( ) perturbed solution with the form assumed for s 1 ( , y, t) to be admissible,ŝ( ) needs to vanish uniformly at the boundaries in accordance with the boundary conditions of the base solution, Equation (50). This gives the appropriate boundary conditions that the solutions of the above ODE need to satisfy as, Upon discretization employing, for instance, a finite difference scheme one can construct an eigen value problem whose solution set for a givenk is the vector of growth-rates given by the spectrum, A , of the finite difference matrix with dimension determined by the discretization size. According to modal stability analysis for anyk ∈ R, if sup{ℜ( A )} > 0 then the corresponding perturbation grows exponentially in time according to Equation (54) and if sup{ℜ( A )} < 0 it decays exponentially.
The physical domain used for resolution of this ODE is the same as the one used for obtaining the TW-solutions, ∈ [−L∕2, L∕2], and the corresponding dimensionless domain is̃∈ [−0.5, 0.5] using a scaling equal to the length of the physical domain, L. In the current study, the scaling used for rendering the problem dimensionless and the numerical scheme for the resolution of the eigen value problem are adopted from Ommi et al. 15 Accordingly, we use a second order accurate finite difference scheme to discretize the derivatives in Equation (55) and a Krylov-Schur algorithm available in the MATLAB suite 50 to resolve the eigen value problem. The dispersion curves can then be plotted as the relation betweenk ∈ R and sup{ℜ(̃A)}, wherek,̃A represent the dimensionless counterparts ofk, A respectively. Figure 5 shows these curves for each base solution corresponding to the TW-solutions resolved in Section 5.1. As one can notice, base solutions for which the one-dimensional solution structure has an overshoot correspond to a dispersion relation with a range of wave numbers having a positive maximal growth rate. This indicates that for certain injection fluxes the fluid displacement is conditionally stable and the morphology of the ensuing fingering pattern is dependent on the range of unstable wave numbers. The critical wave number,k c , determines the unstable wave number range, (0,k c ), and the corresponding range of dimensionless wavelengths, (2 ∕k c , ∞). So within a physical domain of transverse dimension larger than the critical wavelength, 2 L∕k c , one would be able to excite unstable perturbations and observe their non-linear growth. These critical wavelengths for injection fluxes corresponding to s − = 0.85 and 0.80 are 10.46 and 5.23 m respectively.
In order to confirm the results of the linear stability analysis a two dimensional simulation over the same rectangular domain as the one used in Section 6, with the same hydraulic boundary conditions as those described in Figure 7, has been developed starting from an initial imbibition front perturbed by a transversal oscillation whose wavelength is consistent with the value ofk corresponding to the peak of sup{ℜ(̃A)} in Figure 5. The corresponding wavelength is consistent with four oscillations throughout the transversal section of the domain. In Figure 6 the initial imbibition front and the fully developed four fingers at the final stage of the simulation are depicted. F I G U R E 7 Schematic representation of the considered pre-stressed rectangular domain with an initial saturation S + . Boundary condition of the incremental hydraulic problem are also explicitly indicated.

WATER INFILTRATION INTO A DEFORMABLE SOIL
Starting from the weak formulation provided by Equations (37)- (39) and the corresponding numerical implementation given by Equation (47), the coupled problem is solved here, within a rectangular domain (width = 7.5 m, height = 20 m), for a pre-stressed loamy soil characterized by a non-vanishing initial water saturation, perturbing the reference state by a water injection through the top basis of the rectangle, see Figure 7.
The considered boundary conditions of the perturbed problem are therefore the following: (i) on the bottom basis of the rectangle the vertical displacement vanishes and the chemical potential of the fluid is assumed to constantly remain equal to its initial value + = pf (S + ), the normal derivative of the saturation degree is assumed zero; (ii) on the two lateral sides of the rectangle no traction increase with respect to the reference state Δ .n is applied, the fluid flux and the normal derivative of the saturation are zero as well; finally (iii) on the top basis on the rectangle zero traction is imposed and water is injected, again the the normal derivative of the saturation is assumed to vanish. Which means: Δ .e 1 (L, x 2 , t) = 0, k(S r (L, x 2 , t))( eff ,1 (L, x 2 , t) = 0, S r,1 (L, x 2 , t) = 0, (60)  Two different reference configurations which correspond to a heavily overconsolidated (ho) and to a lightly overconsolidated (lo) loamy soil are considered, depending on the value of the initial stress, assumed homogeneous all over the rectangular sample. The two initial states have been chosen so as to describe the different possible fabric changes of fine-grained sediments which can be associated to swelling and the contractant conditions, as discussed in the introduction.
In Table 2 the values of the constitutive parameters of the soil and those of the initial stress are reported, together with the initial saturation S + and the saturation which characterize the inflow at the top of the rectangular sample S − , according to Equation (51). The size of the rectangular domain has been chosen according with the results of the stability analysis conducted in Section 5.2. A perturbation of the initial condition characterized by a wavelength belonging to the unstable region has been introduced so that a single finger is expected to show up in the domain.
In the following the response of the considered loamy soil submitted to an imbibition process starting from the two above mentioned configurations is analyzed following the stress path of two characteristic points fixed in the rectangular domain.
Calculations are carried out using a finite element code derived by the research team from the one, implemented within a Matlab framework. 51 Mesh generation and post-treatment of the results are realized via Gmsh. 52

Imbibition of a heavily overconsolidated loamy soil
Consider the case when the homogeneous initial stress is on the left hand side of the critical state line q = Mp ′ in the (p ′ , q) plane, which means that the initial state of the material is heavily overconsolidated. From this initial configuration and assuming the above mentioned boundary conditions to hold true, the imbibition process is started with a transversal perturbation of the initial condition as the one reported in Figure 8A, as previously mentioned the wavelength of the profile belongs to the unstable range identified in Figure 5. We consider first of all the contour plots of the saturation degree, the volumetric plastic strain and the deviatoric plastic strain at three characteristic time steps, say t = 0, t = 40Δt and t = 60Δt. The unit time step used in the numerical simulations and also reported in Figures 8 and 10 is Δt = 10 4 s. Being the initial state of stress on the left hand side of the critical state line, a softening response is expected, which is responsible for the showing up, at the level of the sample, of a bifurcated solution of the mechanical problem corresponding to a couple of crossed shear bands which form beyond the front. As it is evident from Figure 8D-F the fluid fingering works as a trigger of the mechanical bifurcation. However, due to water accumulation at the top of the sample, a switch from a dilatant into a contractant behavior occurs which is evident in Figure 8E,H: the shear band initially accompanied by positive volumetric strains is progressively absorbed by the contractant zone appearing behind the front.
To better analyse this transition from dilatant to contractant behavior the stress path which characterizes the response of the material at the point P ≡ (3.75, 5) m has been drawn in Figure 9, the zero of the x 2 coordinate having been taken at the top basis of the rectangle. Apparently four different phases can be identified between the initial and the final state of the analysis.
-Phase 1 (running from ′ 0 to ′ 1 , Figure 9, from time step 0 to time step 30): this first phase of the stress path is purely elastic and not affected by the saturation which remains constant, see Figure 10A,B. It is characterized by a decrease of the effective pressure p ′ and a slight decrease of the deviator stress q, see Figure 10I,F. From the observation of the plastic strain charts, see Figure 8, one can observe that the decrease of the effective pressure and the corresponding elastic swelling (positive volumetric strain) in P are essentially consistent with the plastic dilatancy occurring beyond the front, this last being behind the monitoring point. At the same time the deviator stress decreases, this behavior being mainly due to that of the axial stresses in particular 11 and 22 which approach each other, see Figure 10K. During Phase 1 the elastic swelling, therefore, is not a straightforward consequence of the variation of the saturation degree, which indeed is still zero at point P, but it is due to the progressive plastification of the neighboring regions accompanied by their plastic dilatancy.  -Phase 2 (running from ′ 1 to ′ 2 , Figure 9, from time step 30 to time step 50): at the beginning of this phase the state of stress at P attains the yield surface and the volumetric plastic strain starts to grow in the positive domain. As the saturation profile is characterized by a slight undershoot, see Figure 10A, the unvaried shape of the yield surface is therefore due to the countercurrent effects of the desaturation and of the increase of the volumetric plastic strain, see Equation (22). Between time steps 30 and 40, p ′ changes from locally decreasing to locally increasing, the opposite occurs to q. The decrease of p ′ is associated to a positive volumetric elastic strain which sums up with the positive plastic strain, vice-versa the increase of p ′ is associated to a negative volumetric elastic strain as can be observed comparing Figure 10G,H, the growth of the volumetric plastic strain being more important than the one of the total volumetric strain. The behavior of q is still a consequence of that of the axial stresses. When the saturation starts to grow significantly, after the time step 40, two concurrent effects drive the shrinkage of the yield surface, softening accompanied by plastic dilatancy and negative hardening due to saturation, see Equation (22). This implies a strong increase of p ′ and a strong reduction of q. The last time step of the Phase 2 corresponds to the peak of the volumetric plastic strain and to the transition of the effective stress through the critical state line. It is worth to notice that, contrarily to what happens in purely mechanical problems, here the state of effective stress is allowed to cross the critical state line as the loading is driven by a different physical process, in this case the hydraulic imbibition. Finally one can also observe that the peak of the volumetric strain is anticipated with respect to that of the plastic strain because of the above mentioned increase of p ′ which induces an increase of the negative volumetric elastic strain (elastic shrinkage). -Phase 3 (running from ′ 2 to ′ 3 , Figure 9, from time step 50 to time step 58): this third part of the stress path is a straightforward continuation of the previous Phase 2, with p ′ progressively increasing and q decreasing; however here the rate of the volumetric plastic strain becomes negative, which means that the effect of the plastic strain rate and that of the progressive saturation are countercurrent: the first tending to induce swelling of the yield surface, the second its shrinkage. The rate of the total volumetric strain is already negative since the last time steps of Phase 2. As p ′ continues increasing from the previous Phase 2, a more and more negative volumetric elastic strain shows up, as a consequence the total volumetric strain passes from positive to negative before the plastic one. At the end of this Phase 3 the saturation degree has mainly attained its final value so that from now on any change of the preconsolidation pressure and consequently of the yield surface is just driven by the value of the plastic strain rate. The deviatoric plastic strain, as well as the total one, grow because of stress localization within the shear band. It is worth to notice that during imbibition, say during the transition from the initial saturation S + = 0.44 to the one corresponding to the imposed inflow S − = 0.8 the equivalent pore pressure decreases. This is a straightforward consequence of the non-monotonic profile of the pore fluid chemical potential, see Figures 2 and 4. On the other hand, the profile of the effective chemical potential is characterized by a kind of double well profile in the same time interval, which can be also followed along the curve of the pore fluid chemical potential of Figure 4. -Phase 4 (running from ′ 3 to ′ 4 , Figure 9, from time step 58 to time step 86): this fourth part of the stress path is just driven by the increase of the effective stress; the saturation, and consequently the equivalent pore pressure , have already attained their target values induced by the injected flow according to Equation (61). A two step growth of the deviator stress q reflects into a corresponding increment of the deviatoric plastic strain, see   Figure 10E, and of the deviatoric total strain, see Figure 10D; on the other hand the effective pressure p ′ first remains almost constant and then starts growing which implies the slope of the volumetric plastic strain to first remain constant and then to increase. The behavior of the total pressure clearly reflects that of the effective one, see Figure 10J.
In a similar way the stress path at point Q ≡ (0.5, 8) m is reported in Figure 11. It is worth to underline that within the interval of considered time steps this point has not yet been reached by the saturation front but conversely it is affected by deviatoric and volumetric plastic strains which arise within the shear band which nucleates beyond the saturation front. In this case the elastic phase (between ′ 0 and ′ 1 ) corresponds to an increase of the effective pressure p ′ and consequently to a negative (contractant) elastic volumetric strain, due to the contractant behavior of the layers above the monitoring point. Once plastic strains are attained p ′ increases and the stress path approaches the critical state (between ′ 1 and ′ 2 ), asymptotically tending to it with a decreasing deviator stress q (between ′ 2 and ′ 3 ). This behavior is definitely due to the shrinkage of the yield surface and the constraint for stress to remain on the left hand side of the critical state line.

Imbibition of a lightly overconsolidated loamy soil
Consider now the case when the homogeneous initial stress is on the right hand side of the critical state line q = Mp ′ in the (p ′ , q) plane, which means that the initial state of the material is lightly overconsolidated. In a similar way as previously done in Section 6.1, the imbibition process is started with a transversal perturbation of the initial condition, still characterized by a cosinusoidal profile whose wavelength belongs to the unstable range identified in Figure 5. We consider first of all the contour plots of the saturation degree, the volumetric plastic strain and the deviatoric plastic strain at three characteristic time steps, say t = 0, t = 40Δt and t = 60Δt, still with Δt = 10 4 s, see Figure 12. Being the initial state of stress on the right hand side of the critical state line, the two effects of negative hardening induced by saturation and mechanical hardening act in two opposite directions. In this case no bifurcation in the solution of the mechanical problem is expected and the fluid fingering instability just induces a localization of (plastic) strains within the finger itself. The mechanical response is always contractant. As in Section 6.1 the stress path which characterizes the response of the material at the point P ≡ (3.75, 5) m has been drawn in Figure 13. Apparently three different phases can be identified between the initial and the final state of our analysis.
-Phase 1 (running from ′ 0 to ′ 1 , Figure 13, from time step 0 to time step 32): as already seen in the case of heavily overconsolidated soils discussed in Section 6.1, Phase 1 is characterized by a purely elastic response, as the saturation degree remains almost unaltered with respect to initial conditions, see Figure 14A, exhibiting even a weak undershoot between time step 28 and time step 32. p and p ′ increase, a slight elastic compression is observed together with an elastic deviatoric strain, see Figure 14H,E. This behavior is driven by the plastic contraction that is occurring in the layers above the monitoring point which have been already saturated while the point P still maintains its initial saturation. The deviator stress q increases mainly because of the increase of  This means that the state of stress progressively moves towards the yield surface, this last starting to shrink mostly at the same time as plastic strains are triggered.
-Phase 2 (running from ′ 1 to ′ 2 , Figure 13, from time step 32 to time step 50): contractant plastic strains appear and progressively overwhelm the elastic ones which are an order of magnitude smaller than them and remain still contractant, see Figure 14I. A similar behavior can therefore observed for the total and the plastic volumetric strain, see Figure 14G,H. A different behavior is observed for the deviator strain: the plastic component increases while the total one decreases, see Figure 14E,D. This can be interpreted considering the corresponding decrease of the deviator stress which reflecting the reduction of the deviatoric stress implies the deviatoric strains to be negative. The decrease of q, see Figure 14F is mainly due to the evolution of the axial stresses in particular those along the x 1 and the x 2 directions, which approach one to the other, see Figure 14K. During this Phase 2 the saturation grows from its initial value to a value slightly larger than the one imposed at the boundary and a small overshoot is detected, see Figure 14A. -Phase 3 (running from ′ 2 to ′ 3 , Figure 13, from time step 50 to time step 86): now the saturation has attained its target value and only small strains variations are observed: p and p ′ slightly increase while q remains almost constant as the difference between the longitudinal stresses 11 and 33 , which mainly affect its evolution, is also constant.
In Figure 15 the stress path at point R ≡ (2, 5) m is reported which is definitely similar to the one at P. As a matter of fact in this case no perturbation of stresses and strains can be detected beyond the front, therefore the only difference resides in the time step which corresponds to the different transitions already discussed for point P.

CONCLUSIONS
In this paper, we present the first results in modeling the effects of hydraulic instabilities, in the form of fingering, on the mechanical response of a fine-grained soil. The analysis has been conducted in the particular case of the imbibition problem, thinking of possible applications to soil remediation from pollution, in particular from NAPLs. It has been proven that fingering instabilities, even in the absence of any heterogeneity of the permeability and the retention properties of the porous medium, are capable to trigger plastic strain localization, in particular dilatant shear bands, when the soil is heavily overconsolidated, and contractant regions reproducing similar patterns as those of saturation when the soil is lightly overconsolidated. To develop this analysis a generalized version of the Darcy law has been adopted deduced from a phase field description of partial saturation and incorporating a gradient regularizing contribution which introduced additional non-linearity but also a higher order diffusion term. The corresponding coupled problem has been implemented for the first time within a finite element code, adopting a mixed approach. Even if the numerical implementation is in some aspects quite elementary, we refer in particular to the choice of linear Lagrange elements for the whole set of nodal variables, which did not allow to account for the effects of porosity  gradients, interesting characteristics of coupled hydro-mechanical instabilities have been observed that deserve on the one hand further investigations from the numerical point of view and on the other one experimental validation. From the numerical point of view, in particular, the introduction of a strain gradient model will for sure be beneficial to model the showing up of shear band instabilities induced by saturation and to compare the regularizing effect of the porosity gradient with that of the strain gradient.
Moreover having in mind applications of this approach to the simulation of the drainage process that typically occurs at the interface between an aquifer rock, hosting stored hydrocarbons, and the tight caprock saturated by brine or the infiltration of hydrogen produced by corrosion of the canister hosting radioactive wastes against engineered barriers, an extended formulation of the model and a modified numerical implementation should be introduced in order to account for the compressibility of the gas, that up to now has just been considered as a passive phase.

ACKNOWLEDGMENTS
The authors would like to acknowledge the support of the French National Research Agency (ANR), project STOWENG (Project-ANR-18-CE05-0033).

DATA AVAILABILITY STATEMENT
The data that support the findings of this study will be made available to interested researchers by the corresponding author.