Bridging Effective Stress and Soil Water Retention Equations in Deforming Unsaturated Porous Media: A Thermodynamic Approach

The finite deformation of an unsaturated porous medium is analysed from first principles of mixture theory. An expression for Bishop’s effective stress is derived from (1) the deformation-dependent Brooks and Corey’s water retention curve and (2) the restrictions on the constitutive relationships of an unsaturated medium subject to finite deformation. The resulting expression for the effective stress parameter \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\chi $$\end{document}χ is reasonably consistent with experimental data from the literature. Hence, it is shown that Bishop’s equation is constitutively linked to water retention curves in deforming media.


Introduction
Unsaturated soil mechanics is of utmost importance in many applications. Examples are rainfall-induced landslides and slope failures, the settlement of foundations on (and bearing capacity of) the unsaturated soils, and wetting-and drying-induced volume changes in the expansive and collapsible soils. This explains why the mechanical behaviour of unsaturated soils has been a subject of interest in recent decades.
In order to construct a fundamental framework for modelling the mechanical behaviour of unsaturated soils, a clear understanding of coupling between hydraulics and stress states of unsaturated soils is necessary. On the one hand, the volume fraction of different phases (air, water and solid) appears in the definition of stress measures (stress state variables). On the other hand, the variation of volume fractions themselves depends on the current stress level in a deformable unsaturated porous medium. This coupling is commonly referred to as hydromechanical coupling. Recent experimental studies in unsaturated soil mechanics have pointed to such coupling. Also, a number of modelling studies have tried to take it into account (Gallipoli et al. 2003a, b;Sun et al. 2007Sun et al. , 2008Sheng and Zhou 2011;Zhou et al. 2012;Sun and Sun 2012;Casini et al. 2013). The hydromechanical coupling frameworks introduced to date for unsaturated soils are mainly of an empirical nature, and the need for a rigorous thermodynamic framework to address the hydromechanical coupling is still being felt. As the term hydromechanical coupling explicitly points out, it is a two-way coupling. Therefore, an essential question is that as one can measure and obtain the volume fraction of different phases in unsaturated soils at certain stress levels, how to formulate the stress state variables which account for the stress-level dependency of the amount of different phases. A comprehensive formulation for this purpose can only be obtained based on the principles of thermodynamics, by means of balance laws, and through a rigorous mathematical framework rather than a heuristic approach. The starting point is to understand how volume fractions of different phases change in an unsaturated soil and how they can be defined based on other physical state variables. Moreover, it is necessary to know how such relationship is affected by the stress level. One major characteristic of unsaturated soils is the soil water retention curve (SWRC), which specifies capillary pressure value associated with the water saturation degree. In recent years, experimental evidences have clearly shown that soil water retention curves are dependent on the stress level (Romero and Vaunat 2000;Karube and Kawai 2001;Gallipoli et al. 2003b;Tarantino and Tombolato 2005;Nuth and Laloui 2008;Tarantino 2009;Uchaipichat 2010). Tarantino (2009) has considered the dependency of air entry value on stress level and has offered a SWRC equation for deformable media. He has modified the van Genuchten equation for SWRC to account for soil deformation and SWRC dependency on the stress level. Van Genuchten equation, in its original form, reads (Genuchten 1980): where α, n and m are fitting parameters and Θ = S r −S ro 1−S ro is the effective degree of saturation which is defined based on S r (water saturation) and residual degree of saturation S ro (the water saturation at very high suction values which is mainly in the form of water films surrounding the particles). Θ can be replaced by S r for coarse-grained soils as a good approximation. In fact, α can be related to the inverse of air entry capillary pressure, while n is related to pore-size distribution. Parameter m is linked to n and is usually set to n−1 n . Tarantino (2009) considered the dependency of α (or air entry pressure) on stress level by relating it to the void ratio. He considered zero residual saturation for simplicity. Consequently, his modified formulation of van Genuchten equation reads (Tarantino 2009): where e is void ratio, a, n and b are fitting parameters. Recently, Salager et al. (2010Salager et al. ( , 2013 experimentally studied the void ratio dependency of retention curves and introduced the water retention surface, which clearly illustrates the change of saturation degree-capillary pressure relationship caused by a change in void ratio. Nuth and Laloui (2008) proposed a constitutive relationship between air entry value and void ratio as well as the elasto-plastic analogy in the degree of saturation versus capillary pressure relationship (an analogy with mechanical hysteresis) to model hysteretic water retention curves in deformable soils. By means of an empirical equation for the effective stress, Masin (2010) formulated and investigated the void ratio dependency of degree of saturation for various capillary pressure values. All these studies point to the fact that soil water retention curve changes as a result of a change in void ratio. Thus, stress-level dependency of the soil water retention curve has to be taken into account for a proper hydromechanical modelling. As indicated before, we expect not only to observe the stress-level dependency of the amount of different phases in unsaturated soils but also the stress measures have to be dependent on the amount of different phases. Such dependency has commonly been introduced to the current hydromechanical modelling frameworks using soil water retention (capillary pressure-saturation) curves.
There are different formulations for relating the state of stress in unsaturated soils to capillary pressure. Usually, this is done by defining an effective stress (Gallipoli et al. 2003a;Khalili et al. 2004;Gens et al. 2006;Sun et al. 2007Sun et al. , 2008. Among different stress state variables that have been proposed for modelling unsaturated soils, the effective stress has proved to work satisfactorily for modelling unsaturated soils (Khalili and Khabbaz 1998;Khalili et al. 2004;Lu 2008). It provides a consistent framework for both saturated and unsaturated conditions and is able to incorporate hydromechanical coupling. In most effective stress formulations, Bishop's proposal for the effective stress equation (3) is employed along with different empirical equations for the effective stress parameter (Khalili et al. 2008;Uchaipichat 2010).
where σ is the total stress tensor, p air is the pore air pressure, χ is the so-called effective stress parameter, and p c is the matric suction (capillary pressure) and I denotes the unit tensor. In all such empirical equations, the effective stress parameter is explicitly expressed either in terms of the soil saturation or in terms of capillary pressure and parameters of soil water retention curve. For instance, Khalili and Khabbaz (1998) have introduced the following form for the effective stress parameter χ: in which the p c,ae is the air entry value and Ω is a fitting parameter, which is suggested to be taken equal to 0.55 by Khalili and Khabbaz (1998) based on a large number of experimental data points they assessed. A similar formula exists for the soil water retention curve which was proposed by Brooks and Corey (1964): Comparison of Eqs. (4) and (5) indicates that the effective stress parameter and the effective degree of saturation are related. In fact, in the literature we find various equations for the effective stress parameter as a function of (or simply being equal to) effective degree of saturation (Lu et al. 2010). Vanapalli et al. (1996), for instance, expressed the contribution of the suction to shear strength of unsaturated soils τ us = χ p c tan φ , as follows: where κ is a material parameter and φ is the effective friction angle, the same as the one for saturated soils. In some other cases, the effective stress parameter is expressed alternatively as a function of capillary pressure and air entry pressure as previously discussed (Khalili and Khabbaz 1998). From this brief review, we can conclude that the stress measures in unsaturated soils are closely connected to SWRC equation. But, a number of questions remain unanswered. One may ask how the stress-level dependency of soil water retention curve would be introduced to the current effective stress measures, which are empirically formulated. Moreover, it is not either known or proved how the coefficients appearing in the soil water retention curve are related to the empirical coefficients appearing in the proposed equations for the effective stress. The effect of net stress (stress level) on the soil water retention properties and consequently on the effective stress parameter has been mostly studied experimentally (Oh and Lu 2014) rather than through rigorous theoretical frameworks such as mixture theory. While some researchers such as Lu et al. (2010) and Nikooee et al. (2013) have resorted to an energetic and thermodynamic approach to formulate the effective stress in unsaturated soils, they have not looked into the stress-level dependency of hydraulic parameters directly and possible higher-order couplings it may bring about. A theoretical study is therefore essential to identify the scope of applicability of current formulations and to arrive at a most comprehensive formulation for the stress measures accounting for the all aspects of hydromechanical coupling.
In the current study, hence, we employ the mixture theory in order to study the mutual relationship between SWRC equation and a stress measure for an unsaturated porous medium. In fact, in this study we are limited to current hydraulic formulations that account for the hydraulic hysteresis (different equations for drying, wetting and scanning branches need to be introduced). Also, the contribution of interfaces to the mechanical behaviour of the unsaturated porous medium has been neglected. This can be a reasonable assumption when coarse-grained soils are considered. In fine-grained soils especially in scanning curves and/or in low saturation degrees, the contribution of interfaces can be considerable (Nikooee et al. 2013;Likos 2014;Gray and Schrefler 2001;Dangla and Coussy 2002). In what follows, first the theoretical framework is introduced. Then, the expected coupling is derived based on the definition of the strain energy function of the mixture. Next, coupling/connection between the stress measure and soil water retention curve is discussed and compared to current formulations. Finally, concluding remarks and some suggestions for future studies are provided.

Basic Assumptions
We consider three phases: a solid (superscript s), a wetting fluid (superscript w) and a nonwetting fluid (superscript nw). The solid deforms from an initial configuration X to a current configuration x. The deformation gradient tensor of the solid In Eq. (7), the usual superscript s is left out because in the subsequent equations we do not consider deformation gradient tensors of other constituents. The determinant of the deformation gradient tensor represents the volume change of the mixture

Conservation Laws and Constitutive Laws
Excluding mass transfer between phases, the mass balance of each phase is then written as: in which ρ α is the apparent density and v α the velocity of component α. In analogy to Vankan et al. (1996), Huyghe and Janssen (1997) and Huyghe and Bovendeerd (2004), we refer current descriptors of the mixture with respect to an initial state of the porous solid. If we introduce the Lagrangian apparent density [compare to eq. (4.1), in Bowen (1980)] per unit initial volume of the mixture, we can rewrite the mass balance equation (9) as follows: In order to simplify our task, we assume the solid and the wetting phase (usually water) to be incompressible, so that the intrinsic density 1 2 ) is constant. φ α is the volume fraction of phase α. For the non-wetting phase (often air), we assume the phase to be drained, i.e. the pressure of the air is atmospheric everywhere. Because of the constant pressure and the assumption of isotherm, the intrinsic density of the air is constant in time and space. Because of Eqs. (12) and (13), the mass conservation equations reduce to volume conservation equations: and D α Dt is the time derivative for an observer fixed to the constituent α. The total mass balance is obtained by adding the three mass balances together: Neglecting body forces and inertia, the momentum balance takes the form: in which σ α is the partial stress tensor of constituent α and π α is the momentum imparted upon phase α by the other phases, which after summation over the three components yields: if use is made of the balance condition: Although the non-wetting phase is drained and non-wetting pressure is atmospheric pressure, the momentum interaction π nw is nonzero if gradients of non-wetting volume fraction exist. However, the gradient in non-wetting volume fraction is a self-equilibrated term that does not induce flow. Balance of moment of momentum requires that the stress tensor σ be symmetric.
If no moment of momentum interaction between components occurs, the partial stresses σ α also are symmetric. In this paper, we assume all partial stresses to be symmetric. Under isothermal conditions, the entropy inequality for a unit volume of mixture reads Bowen (1980): in which F α is the Helmholtz free energy of constituent α per unit mass constituent. We introduce the strain energy function [compare to Biot (1972, eq. (3.22)) and Bowen (1980, eq. (4.2))] as the Helmholtz free energy of a mixture volume which in the initial state of the solid equals unity. ψ α is the Helmholtz free energy of constituent α per unit mixture volume. Rewriting the inequality (20) for the entropy production per initial mixture volume-i.e. we multiply inequality (20) by the relative volume change J -we find: The entropy inequality should hold for an arbitrary state of the mixture, complying with the balance laws. We introduce into the entropy inequality (22) the total mass balance (16) using a Lagrange multiplier p: with σ eff the effective stress: The Lagrange multiplier p has obviously the dimension of a pressure. We choose as independent variables the Green strain E, the Lagrangian form of the volume fractions of the fluid and air Φ w and Φ nw , and of the relative velocity We apply the principle of equipresence, i.e. all dependent variables depend on all independent variables, unless the entropy inequality requires otherwise: We apply the chain rule for time differentiation of W : in which μ w is the chemical potential of the wetting fluid: and μ nw is the chemical potential of the non-wetting phase: Equation (31) should be true for any value of the state variables. Close inspection of the choice of independent variables and the inequality (31), reveals that the first term of (31) is linear in the solid velocity gradient ∇v s , the second term linear in D s Dt v ws , the third term linear in D s Dt v nws , the fourth term linear in the relative velocity gradients ∇(v w − v s ) and the fifth term linear in the relative velocity gradients ∇(v nw − v s ). Therefore, by a standard argument, we find: leaving as inequality: Equation (34) indicates that the effective stress of the mixture can be derived from a strain energy function W which represents the free energy of the mixture. From Eqs. 24 and 34, one can infer that the pressure p can be interpreted as the pressure present at the interface between solid and wetting fluid. Indeed, the stress in the solid has a hydrostatic component that does not induce deformation of the incompressible solid. This is the pressure exerted on the solid by the wetting fluid. Equation (34) is the constitutive law of the solid skeleton, which unlike the solid itself is compressible as its porosity can change during deformation.
Equations (35) and (36) show that the strain energy function cannot depend on the relative velocities. Thus, the stress of an unsaturated medium can be derived from a regular strain energy function, which physically has the same meaning as in solid materials, but which can depend on both strain and composition of the medium. According to Eqs. (37) and (38), the partial stress of the fluid and the gas is a scalar. Transforming the relative velocity to its Lagrangian equivalent, we find instead of (39): in which ∇ 0 = F T · ∇ is the gradient operator with respect to the initial configuration. Because Reynolds numbers are usually very low, we assume that the system is not too far from equilibrium. Hence, we can express the dissipation (40) associated with relative flow of fluid as a quadratic function of the relative velocities: B βγ is a positive definite matrix of frictional coefficients. Substituting Eq. (37) into Eq. (41) yields Lagrangian forms of the classical equations of irreversible thermodynamics:

Corey-Brooks Relationship for Capillary Pressure and Bishop's Equation
The pores of the incompressible solid are filled with an incompressible wetting fluid and air at atmospheric pressure. The sample as a whole can change its volume by expelling wetting or non-wetting fluid. The sample is in contact with a water reservoir at zero pressure. At the interface between the reservoir and the sample, the chemical potential μ w of the wetting phase is the same inside and outside: the chemical potential of the non-wetting phase is: The capillary pressure difference between the wetting fluid and non-wetting fluid is [Bowen 1980[Bowen , p. 1141]: in which p c is the capillary pressure of the wetting fluid and Φ w the Lagrangian volume fraction of the wetting fluid: The saturation is defined as: According to Gallipoli (2012), Brooks-Corey can be written for deformable media as, for p c ≥ p c,ae (49) in which we assume that the effective saturation is equal to the saturation and p c,ae is the air entry pressure dependent on the void ratio yielding a water retention curve for deformable media: Inverting this relationship yields or, after substitution into Eq. (45): Differentiating the above equation with respect to strain yields: In Eq. (57), we use the identity for any tensor A. This results in the following expression for the stress: in which σ sat is an integration constant dependent on the strain E, representing the effective stress in the saturated case. Equation (60) is rewritten in the form: or, after integration: Considering that σ sat = σ − p c,ae I we can rewrite Eq. (62) Making use of Eq. (50) Fig. 1 Soil water retention curve at different net stress levels (weathered granite), experimental data are from Lee et al. (2005) soils samples associated with each net stress level (the same as what we consider for our proposed equation), and for the exponent of their equation, we keep it constant and equal to their suggested best fit value of 0.55. This will make Ω a material independent parameter (this is in line with the assumption of previous researchers, e.g. that of Masin (2010), who investigated the void ratio dependency of soil water retention curves using an effective stress approach). Lee et al. (2005) made a novel device to determine the soil water retention curve at different net stress levels. However, the current experimental practice in unsaturated soil mechanics is solely to determine the soil water retention curve using plate pressure apparatus at zero net stress level. Of course, the water retention properties at different net stress levels can also be obtained from the triaxial testing apparatus, but the measurement of the water volume change in this device is not as accurate as plate pressure apparatus, and therefore, there are limited data in the literature on the soil water retention properties at different net stress levels. Hereafter, in order to explore different experimental datasets in the literature, we present a systematic approach to use the proposed formula for the determination of the effective stress parameter at different net stress and suction values. For weathered granite, the experimental data on the water retention curve at four net stress levels, namely 0, 100, 200, and 300 kPa, are available (Fig. 1). Therefore, first using the effective parameter values at different suction levels, and at zero net stress the value of coupling parameter is obtained, as the value of air entry value at higher net stress levels can be directly obtained from soil water retention curves experimentally measured by Lee et al. (2005) at higher net stress levels. For this purpose, the Brooks-Corey model (Eq. 48) has been used to simulate the experimental data and to obtain the air entry value. The resulting effective stress parameter for weathered granite obtained from the proposed equation (Eq. 67) is depicted in Fig. 2. Table 1 presents the values of different parameters required for predicting the effective parameter for weathered granite. The air entry values in this approach were all determined directly from the available  experimental data, and the slope of SWRC was considered to be almost the same under different net stress levels (consistent with the assumptions in the previous studies [see, e.g. Tarantino (2009)]. Therefore, the coupling parameter was the only parameter which was found from calibrating the proposed equation with the experimental data at zero net stress level. Then for higher net stress levels, the proposed equation was evaluated, which resulted in a fairly well agreement with the experimental data. For Kidston tailing and Talybont clay, as the parameters of soil water retention curve at higher net stress levels were not available in the experimental results, a systematic approach was employed. For Kidston tailing, first the slope of soil water retention curve was determined directly from soil water retention curve measured at zero net stress level. Then, the slope of the SWRC and the proposed equation for the effective stress parameter were used to determine the air entry value at 30 kPa and coupling term. Having the coupling term, and the slope of soil water retention curve, then for higher net stress levels, namely, 125 and 250 kPa, the sole parameter which needed to be found was air entry value which could directly be obtained by the use of the proposed formula for the effective stress parameter, while the performance of the proposed equation could also be assessed (see Fig. 3; Table 2, for a comparison of the proposed equation and the experimental data for Kidston tailing as well as the values of different parameters). For Talybont clay, the values of slope of SWRC at zero net stress and the effective stress parameter at zero net stress level were used to find the coupling parameter, and then, the air entry value at 204 kPa (30 psi) was obtained by fitting the proposed equation to the effective stress values at higher net stress level. The comparison between the proposed equation and  the experimental data of Talybont clay (Fig. 4, Table 3) shows some discrepancy which can root in the existence of micro-and macro-pores (a double porous structure). In order to be able to address the presence of double porous structure, the use of more advanced retention equations may be necessary which account for the presence of steps (plateaus) in the water retention curves, and accordingly in the effective stress parameter versus suction. However, the overall agreement of the proposed equation with the experimental data is good enough to consider it as a proper approach to address the variation of the effective stress parameter with the net stress and capillary pressure. It is noteworthy that the equation proposed by Khalili and Khabbaz (1998) works perfectly for Kidston tailing at different net stress levels (Fig. 3), but it shows less favourable results for weathered granite and Talybont clay. As we just discussed for Talybont clay, a double porous structure can be the possible reason for a poor prediction by both approaches. The impreciseness in prediction for weathered granite as compared to our approach roots in the empirical nature of the equation proposed by Khalili and Khabbaz (1998) and the assumptions involved.

Conclusion
As indicated in the current study, a direct link between soil water retention curve and the effective stress formulation for unsaturated soils can be established via principles of thermodynamics and mixture theory. The resulting formula contains coefficients that have physical meaning and can be obtained from soil water retention curve. Because hydromechanical cou-  pling is introduced through a theoretical framework, the stress-level dependency of soil water retention curve is explicitly taken into account in the resulting effective stress formulation and one does not need to heuristically examine its coefficients for new conditions. It is, therefore, superior to current formulations in the literature, which are empirically obtained. The proposed equation holds as long as its assumptions hold, and thus, its scope of application can be easily explored, whereas for the empirical equations of the effective stress parameter, it is difficult to notice the range of their coefficients, physical meaning and connection to other parameters and variables. It is noteworthy that the proposed thermodynamic framework allows us to connect the effective stress equation to that of soil water retention curve in each hydraulic path (i.e. wetting, drying and scanning) which means theoretically it can capture the hydraulic hysteresis. In this study, we, however, focused on the variation of effective stress parameter in one hydraulic path. Future study is recommended to examine the efficacy and applicability of the proposed equation in different hydraulic paths. Last but not the least, the presence of disconnected phases and/or interfaces needs to be taken into account in most general case which we neglected in this study to simplify the framework. Accounting for hydromechanical coupling while considering the presence of all interfaces and phases will be a next step complementing this study. Of course, in order to be able to look into different phases and interfaces and their variation under different loading conditions, further progress in the current experimental techniques is essential.
of Economic Affairs (Project 12538). EN and SMH authors would like to thank the European Research Council for ERC advanced grant which has supported the fundamental porous media research at Utrecht University and their contribution to this study (Agreement No. 341225).
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.