Physicochemical effects remain elusive in the theory of porous media. A continuum theory is presented for modeling the chemomechanical behavior of unsaturated soils with chemomechanical coupling. The new theory effectively explains many salient phenomena occurring in soils and describes the multiple coupled processes in the vadose zone.

## Abstract

A theoretical framework is presented for modeling the chemomechanical behavior of multiphase porous media, in general, and unsaturated soils, in particular, which can address skeletal deformation, fluid flow, heat conduction, solute diffusion, chemical reaction, and phase transition in a consistent and systematic way. A general expression is derived for the electrochemical potential of a fluid species with explicitly accounting for the effects of osmosis, capillarity, and adsorption. The equilibrium behavior of porous media is investigated, and the composition of pore water pressure is identified. Explicit formulations are developed for the effective stress and intergranular stress, with consideration of physicochemical effects. It is shown that the negative water pressure measured by a conventional transducer can be significantly different than the true pore water pressure. It is also theoretically revealed that, other than the soil water characteristic function, a new pressure (or potential) function accounting for the physicochemical effects is generally required in analyzing the coupled chemomechanical processes in unsaturated soils. The new theory is capable of effectively explaining many salient phenomena occurring in water-saturated porous media with a degree of saturation varying from an extremely low value to 100%, including Donnan osmosis, capillary fringe, air entry value, initial hydraulic head during seepage, and pressure solution. The new theory can be used to analyze the multiple coupled physical and chemical processes in the vadose zone.

**Chemomechanical behavior** of porous media with multiphases and multispecies is of great interest in many diverse fields of science and engineering. To mention a few, this includes the industries of nuclear, hazardous, and municipal waste isolation; petroleum and gas extraction; technologies of methane gas hydrates exploitation; CO_{2} sequestration; and weak soil reinforcement, landform stability assessment, structural material durability, and weathering of rock masses. Understanding, controlling, and predicting the long-term effect of physicochemical processes on the mechanical performance of geomaterials is becoming an indispensable part of environmental impact assessment and performance assessment analysis. Thus, there is a clear need to develop comprehensive chemohydromechanical mathematical modeling capacities that are able to address realistically the reactive multispecies transport, multiphase flow, chemical reaction, phase transition, chemically induced deformation, and other related physicochemical processes in deformable soils.

The coupling of multiple processes in porous media, including skeletal deformation, seepage, diffusion, and heat conduction, has been extensively investigated, and the poromechanic theory that is capable of describing these coupled physical processes has been very well developed (e.g., Coussy, 2004). Thus far, however, the behavior of porous media with multiphase and multispecies where physicochemical effects come into play remains elusive. Earlier research and practices were mainly concerned with reactive flow and transport, and numerous flow and transport models have been developed, without addressing mechanical issues (e.g., Spycher and Sonnenthal, 2003; Steefel et al., 2005; Xu and Pruess, 2001). With regard to the chemomechanical behavior of porous media, a few research studies have been performed with a focus on the chemoplasticity of fully saturated soils (e.g., Hueckel, 2002; Loret et al., 2004; Witteveen et al., 2013). In addition, the cross-scale effect and multiprocesses coupling of multiphase systems has been an active field of research on experimental, theoretical, and numerical bases (e.g., Di Maio et al., 2002; Coussy, 2010).

Over the past two decades, numerous efforts have been devoted to developing general theoretical frameworks for modeling multiphase and multispecies porous media. Among these, we mention the biological tissues-oriented theory of porous media developed by Huyghe and his coworkers (Huyghe and Janssen, 1997; Huyghe et al., 2004, 2007) and the hybrid theory of swelling porous media (Bennethum and Cushman, 1996a,b; 2002a,b; Bennethum et al., 1996, 2000; Bennethum, 2007, 2012). Through these efforts, general macroscopic governing equations have been developed that can be used to address the physicochemical and electrochemical coupling in multiphasic and multispecies systems. In addition, the significances and implications of pressures and potentials have been well addressed.

Despite of this progress, a comprehensive theoretical framework is still lacking for modeling the chemomechanical behavior of unsaturated soils with variable saturation. Here our hypothesis is that when the saturation varies from 100% to a low extreme, the composition and concentration of a pore fluid are variable, resulting in intensive physicochemical effects in the soils. Indeed, for the unsaturated soils with low saturation, the pore water pressure (PWP) has not been very well defined, and in current practice, the methods for PWP measurement are used quite ambiguously, even without distinguishing the difference between the potential and the pressure. To overcome such difficulties, several fundamental issues have to be resolved, including:

How do we reconcile the concepts of potentials and stresses (or pressures) in a unique framework that is used to address the coupled thermo-hydro-chemo-mechanical processes in unsaturated soils?

How do we characterize the chemical potential of a species in the unsaturated soils with multispecies, where osmotic, capillary, and adsorptive effects are important?

What is the explicit expression of the effective stresses when compositions and concentrations of pore water are variable?

Solving the first problem depends largely on the answers to the last two issues. Based on fundamental thermodynamic principles and an averaging procedure, Nitao and Bear (1996) had rigorously derived mathematical formulations for the potentials of unsaturated soils and provided a complete set of governing equations for flow and transport processes in the soils. Although this framework describes the flow and transport processes in a multiphase system, it does not address the mechanical issues.

As one of the fundamental concepts in the classical soil mechanics, the effective stress concept plays a crucial role in describing the mechanical behavior of saturated soils. After Terzaghi’s proposal for it, intensive efforts have been made to extending the effective stress concept to unsaturated soil problems (e.g., Bishop, 1959; Kohgo et al., 1993; Khalili et al., 2004; Laloui and Nuth, 2009). Perhaps the most commonly used formulation of the effective stress for unsaturated soils is the one proposed by Bishop (1959). This equation states that the effective stress of unsaturated soils is equal to the total stress minus the averaged pore pressure that is the average of pore air and pore water pressure weighted by factor χ. The validity of Bishop effective equation have been recently examined by Gray and Miller (2007), Gray and Schrefler (2007), Gray et al. (2009), and Nikooee et al. (2012), based on the principles of thermodynamics and a local averaging procedure with explicit consideration of interfacial effects. These authors have shown that Bishop’s effective stress can be theoretically recovered only if the interfacial tension terms are neglected.

In a somewhat intuitive way, Lu and Likos (2006) and Lu et al (2010) have proposed a new formulation for the effective stress of unsaturated soils, which includes a stress term called the “suction stress.” Remarkably, the suction stress accounts for the effect of surface tension, negative pore pressure, electrochemical interactions, and other factors related to unsaturated soils. The applicability of the suction stress concept has been validated based on the shear strength properties of unsaturated soils extensively collected from the literature. Thus far, however, a solid theoretical basis has yet to be built for this concept.

On the basis of the approach of the mixture hybrid theory, a detailed derivation of a comprehensive theoretical framework is presented here that can be used to describe the chemomechanical behavior of unsaturated soils with saturation varying from 100% to a low extreme value. Within this context, the effective stress concept is examined, and the general expression for the chemical potential of a species is developed. To shed new insights into the unsaturated soil behavior, equilibrium behavior of the soils is investigated in detail.

## Thermodynamic Constraints

### Preliminary

The unsaturated soils under consideration are viewed as the porous continua composed of a solid matrix (*s*) with interconnected pores saturated with two immiscible fluids, namely, a liquid (*l*) and a mixed gas (*g*). Unless otherwise specified, symbol α or β represents a bulk phase, i.e., α or β = *s*, *l*, *g*, and *f* denotes a pore fluid (i.e., *l* and *g*). For the mathematical convenience, it is assumed that every bulk phase is composed of the same set of electrically charged species, denoted by *i* = 1, 2, ..., *Z*, respectively. Let be the mass concentration of species *i* in phase α, i.e., , where is the partial mass density of the *i* species in the α phase and ρ^{α} () is the intrinsic mass density of the bulk phase. With this notation, if a species is absent in a bulk phase, its concentration is zero, i.e., . For convenience, all the other variables are summarized in the Appendix.

Formally, the following theoretical derivation follows the procedure adopted by Bennethum and Cushman (1996b) and Bennethum et al. (1996, 2000), with the assumption that the interfaces among bulk phases possess no thermodynamic properties. To avoid redundancy, only key assumptions and general results are given, and useful details of derivation can be found in the above-cited references. Nonetheless, for convenience of reference, the balance equations are summarized in the Appendix.

### Constitutive Assumptions

The derivation starts with the assumptions regard to the specific Helmholtz free energy densities of three bulk phases, namely, *A ^{s}*,

*A*,

^{l}*A*, respectively. To this end, we introduce the principle of phase separation (Passman et al., 1984), stating that the Helmholtz free energy of a phase depends solely on its independent state variables. Based on this principle,

^{g}*A*

^{s}can be assumed as a function of temperature

*T*, volume fraction

*n*

^{s}, intrinsic mass density ρ

^{s}, deformation gradient

**F**

^{s}and species mass concentration (

*k*= 1, 2, …, Z − 1), while

*A*as a function of

^{f}*T*,

*n*, ρ

^{f}*, and (*

^{f}*k*= 1, 2, …, Z − 1). Here it is noted that among all , only Z − 1 components are independent. As argued by Passman et al. (1984), use of the phase separation principle is justified by the fact that individual bulk phases are physically separated in porous media.

*n*

^{s}, ρ

^{s}, and

**F**

^{s}are independent (recalling the solid mass balance, Eq. [A2] in the Appendix). Kinematically, the variation of soil porosity can be additively decomposed into two components: one is due to the skeletal deformation and/or the solid grain compression and the other solely due to the mass exchange (via, say, dissolution and precipitation). As a consequence, variation of

*n*

^{s}also has two contributions, one of which can be solely attributed to the mass exchange (denoted as ). The evolution of is governed by (see Appendix):

*f*.

Inclusion of Lagrangian strain tensor **E**^{s} instead of deformation gradient **F**^{s} in Eq. [2] is due to the requirement of objectivity. Clearly, all the arguments in *A ^{s}* and

*A*are independent. From Eq. [2] and [3], it is noted that the solid matrix is assumed to be elastic and the pore fluids are of Newtonian-type. Viscous effect of pore fluids can be addressed by assuming a dependence of

^{f}*A*on ∇

^{f}**v**

*. For simplicity, however, ∇*

^{f}**v**

*is excluded in Eq. [3].*

^{f}### Constraints by the Second Law

^{α}is the entropy density;

**τ**

^{s}the intrinsic intergranular stress tensor, transited through grain–grain contacts;

*p*

^{α}the thermodynamic pressure of α phase;

**1**is the unit isotropic tensor with components of δ

_{ij}(i.e., Kronecker’s delta); the electrochemical potential tensor of α

_{i}species; the relative chemical potential of species the diffusion velocity; ν

_{i}the valence number of species

*i*; the molar mass of species

*i*; and λ is a Lagrangian multiplier, associated with the electrical neutrality. Variables

**τ**

^{s},

*p*

^{α}, and are defined by

*local*electric potential. Dropping all the diffusion terms, one can cast Eq. [7] into

_{i}] represents the molar concentration [mol/m

^{3}] of species α

_{i}, and represents the reaction coefficients of the

*J*th reaction (

*J*= 1, 2, ...,

*M*). If a species (say α

_{i}) does not participate in the reaction, ; otherwise, for a reactant, and for a product. Let be the rate of the

*J*th reaction, then

^{sα}= 1, if α =

*s*; otherwise, δ

^{sα}= 0;

is the configurational pressure, which is sometimes called the swelling pressure for expansive soils (Bennethum and Weinstein, 2004; Bennethum, 2012); represents the effort (with the same unit as pressure) required to break the chemical bonds of the solid material in the non-deformed condition.

### Thermodynamic Equilibrium

*z*

_{m}}with components ∇

*T*, D

^{s}

*n*

^{f}/D

*t*,

**v**

^{f,s}, and . It is straightforward to prove that the sufficient and necessary conditions for the porous media to attain equilibrium are

Equation [22] is the condition for *thermal* equilibrium, Eq. [23, 24, and 26] are the conditions for *hydromechanical* equilibrium, and Eq. [25] is the conditions for *chemical* equilibrium. Remarkably, for a system in equilibrium, all the thermal, hydromechanical, and chemical equilibrium conditions must be satisfied.

^{sf}= 0, and Eq. [2.25] yields

*j*that can exchange its mass among the solid, the pore liquid, and the pore gas. Symbolically, this exchanging process can be represented as . For exchange , one can assume , and it follows from Eq. [27] that

*f*=

*l*or

*g*), setting , one obtains

An equation similar to Eq. [29] was derived by Bennethum et al. (1996, 2000), although the breaking energy term (i.e., ) appears here for the first time.

Let be the electrochemical potential, and then . Equation [28] states that the chemical potential of any species in a multiphasic system in equilibrium is continuous across the interfaces between bulk phases, which is consistent with the classical result of thermodynamics (Atkins and dePaula, 2002). Equation [29] is a new result, implying that the chemical equilibrium between the solid matrix and a pore fluid can be intervened by the intrinsic intergranular stress. In the left-hand side of Eq. [29], accounts mainly for the elastic energy of the solid matrix, the second term is the elastic energy of the solid material, and the third term represents the energy required to break the chemical bonds of the solid material.

**a**

^{f}is the acceleration of the fluid. Using the state equations for

*G*

^{f},

*p*

^{f}, , and

**μ**

^{f}, Eq. [30] can be cast into

*k*= 1, 2, ...,

*Z*− 1). Noting that

**a**

^{f}= 0 and ∇

*T*= 0 at

*static*equilibrium, one obtains

*i*= 1, 2, ...,

*Z*; sym( ) represents the symmetrical part of a second-order tensor;

**x**and

**x**

_{0}are the spatial coordinates of the point of interest and the reference point, respectively;

*c*is independent of spatial coordinate

**x**. Clearly, at equilibrium, the sum of electrochemical and gravitational potentials is spatially uniform for any species.

### Linear Dissipative Processes

**k**

_{q}, θ

_{f},

**ξ**

_{f}, ζ

_{J}, and (

*k*= 1, 2, ...,

*Z*− 1) are material coefficients. Here, for simplicity, all the cross effects have been neglected. Equation [34] is the generalized Fourier law for porous media with multiphase and multispecies; Eq. [35] describes the rate of the fluid volume fraction; Eq. [36] represents the constraint on the linear momentum exchange between bulk phases; Eq. [37] accounts for the mass exchange processes; Eq. [38] is the generalized Fick’s law for porous media with multiphase and multispecies.

It is remarkable that the second law does not exert any control on the linear momentum exchange of species, i.e., . This result is different from the previous results (e.g., Bennethum and Cushman, 1996b; Bennethum et al., 1996, 2000). It is clear that within this context porous media are treated as highly reactive ones in the sense of Bataille and Kestin (1977).

Coefficient tensor **k**_{q}, **ξ**_{f}, and are symmetric and positively definite, and the scalar coefficients ς_{J} and θ_{f} are non-negative.

## Chemical and Electrochemical Potentials

As discussed above, the chemical potential, , (or electrochemical potential, ) of the species in the pore fluids plays a crucial role in characterizing the behavior of a porous medium with multiphase and multispecies. In the following, general expressions of and are developed for the pore fluids.

### Decomposition of the Specific Helmholtz free energy of pore fluids

In Eq. [3], *A *^{f} is assumed to be a function of *T*, ρ* *^{f}, and (*k* = 1, 2, ..., *Z* − 1), as well as *n *^{f}. The dependence of *A *^{f} on *n *^{f} accounts for the effect of the surface forces associated with interfaces, implying that a solution in the pores differs from that in the free bulk state in that the former is under the influence of the surface forces (Nitao and Bear, 1996). These surface forces include the surface tension on the interfaces and the adsorptive forces stemming from the physicochemical interactions among different bulk phases, including electrostatic forces, Van der Waals attraction, double-layer repulsion, and so on. In general, porous media (e.g., geomaterials) are electrically charged. Hence, it is expected that significant physicochemical interactions occur among mineral surfaces, water dipoles, and electrically charged species. These interactions can modify the potential energy of pore fluids, so that a fluid in the pores has a potential energy different from the fluid free of the surface forces at the same thermodynamic condition (temperature, mass density, and mass fractions).

To gain insight into the character of the specific Helmholtz free energy defined by Eq. [3], we adopted a procedure by Nitao and Bear (1996) to analyze the effect of surface forces on the energy potential. Consider a representative volume (REV) of a pore fluid at point **x**, at temperature *T*, average mass density ρ* *^{f} and mass fraction . Consider also a reservoir outside the porous medium domain containing the same fluid as the REV at the same thermodynamic conditions (*T*, ρ* *^{f}, and ). The reservoir is at the same reference elevation as the REV. It is important to note that, unlike the fluid in the reservoir, the pore fluid is subjected to the surface forces.

^{f}and

*T*during the movement, the same amount of surface forces as in the REV has to be applied to the system so that the system (the transported fluid) has the same potential field as in the REV. As a consequence, the system changes by a pressure Δ

*p*

^{f}and amount of specific heat Δ

*q*(=

*T*Δη

^{f}). In applying the surface force, some amount of work must be externally input to the system. According to the principle of energy conservation, under the condition of constant

*T*and ρ

^{f}, the variation of specific internal energy,

*E*

^{f}, of the system is given by Δ

*E*

^{f}= −Δ

*w*

^{f}+

*T*Δη

^{f}, or equivalently,

*E*

^{f}− = −Δ

*w*

^{f}+

*T*(η

^{f}− ), where Δ

*w*

^{f}is the specific work done against the surface forces, and and are the specific internal energy and specific entropy of the system before applying the surface forces, that is, in . Let Ω

^{f}be the surface energy potential, and Ω

^{f}= − Δ

*w*

^{f}< 0. Then, introducing the Legendre transformation (Eq. [A23] in the Appendix), one obtains

*A*

^{f}and are the specific Helmholtz free energies of the pore fluid in the REV and the fluid in , respectively.

Remarkably, function as defined is equal to the specific Helmholtz free energy of the system, i.e., the fluid in . Hence, once is defined, as a function of the same *T*, ρ^{f} , and as in , Ω^{f} can be fully specified as a function of on *T* and *n*^{f} only. Decomposition of energy potential has usually been exercised in addressing the effect of pore water films on the behavior of unsaturated porous media (e.g., Coussy, 2004, 2010). It is noted, however, that unlike the previous approaches (Coussy, 2004, 2010), here the surface energy potential Ω^{f} is considered part of the free energy potential of the pore fluid, which is subjected to the surface forces.

### A Generic Expression for Chemical Potential

*f*

_{i}. Consider a small representative volume, d

*V*, of a deforming porous medium, whose initial volume before deformation is d

*V*

_{0}. Apparently, d

*V*can be related to d

*V*

_{0}through d

*V*=

*J*d

*V*

_{0}, where

*J*is the Jacobian of the deformation gradient,

**F**, of the solid matrix, i.e.,

*J*= det(

**F**) > 0. Hence, the total specific energy of the fluid with regard to d

*V*

_{0}is

*Jn*

^{f}ρ

^{f}

*A*

^{f}. Let be the specific mass of a species

*f*

_{i}with regard to d

*V*

_{0}, i.e., . According to its very definition, the chemical potential of species

*f*

_{i}can be defined as (e.g., Bennethum and Cushman, 2002b)

*T*,

*p*

^{f}, and molar fraction ; ;

*R*is the universal gas constant; the chemical potential of the species

*f*

_{i}at the pure state. Activity is defined in such a way that it approaches as for a solute species, and for the solvent. is related to through

^{f}, recall the state equation for , i.e., Eq. [20]. Using Eq. [40], one obtains

*n*; is considered as a function of

*T*and

*n*

^{f}. Equation [48] can be recast into

^{+}is a small positive quantity, approaching zero. For the pore liquid, 0

^{+}represents the smallest fluid volume fraction at which the porous medium is considered “dry.”

Assuming the dependence of A ^{f} on is to highlight the fact that the surface energy potential depends on the porosity only through the fluid volume fraction, *n *^{f}, which is equal to multiplied by the degree of saturation, *S*_{f}* *. Noticeably, in Eq. [48] can be viewed as the total surface energy potential of the pore fluid fully occupying the pore space in the porous medium, and is equal to the work required to overcome both the adsorptive and capillary forces in transporting, in a thermodynamically reversible way, the fluid from a reservoir to fully saturate the porous medium at the same *T*, ρ* *^{f}, and . Because the integral in Eq. [50] represents the surface energy potential associated with the capillary forces, it is clear from Eq. [50] that A* *^{f} is indeed the surface energy potential accounting only for the adsorptive forces in a fully saturated porous medium.

As described by Eq. [49], the total surface energy potential has two contributions, i.e., A* *^{f} and C* *^{f}. In the literature (e.g., Tuller et al., 1999), the surface energy (or matric) potential is usually decomposed into two components, which accounts for the effect of adsorptive forces and the effect of capillary forces, respectively. As implied by Eq. [49], it is generally difficult, if not impossible, to clearly distinguish A* *^{f} and C* *^{f}, since both account for the effect of the adsorptive forces. Nevertheless, it is clear that only C* *^{f} depends on the degree of saturation (through *n *^{f}) and accounts for the effects of both the adsorptive and capillary force, whereas A* *^{f} is associated only with the effect of adsorptive forces in the porous media at full saturation.

## Equilibrium Properties

Because is a function of *T*, *p *^{f}, *n *^{f}, and , enforcing equilibrium equations, Eq. [27–29] and Eq. [33] may yield rich colligative properties such as osmotic pressure, undercooling, superheating, and Kelvin effect in porous media. These properties describe the correlations among *T*, *p *^{f}, *n *^{f}, and . Because of its central importance, the Donnan osmotic phenomenon is discussed in detail in the following, In addition, a detailed analysis of different procedures for measuring the negative pore water pressure in unsaturated soils is also given.

### Donnan Osmotic Phenomenon

Consider a soil layer in the vadose zone, in which a groundwater observation well (denoted as W) is drilled down below the water table, as schematically shown in Fig. 1. The soil solution in the observation well keeps in contact, and in equilibrium, with the pore water in the soil. In the following, to distinguish the water solutions in the well and in the soil pores, they are called “equilibrium solution” and “pore water,” respectively. Suppose that the soil has a fixed charge density of *c*_{fix} [mol/m^{3}], which represents the total number of electric charges fixed in a unit volume of the soil. For the purpose of illustration, it is assumed that the fixed charges are *negative* and that both the soil solution and the pore water are composed of a solvent (H_{2}O) and two charged species (*X *^{m−} and *Y *^{n+}), whose valence numbers are −*m* and *n*, respectively.

*c*

_{0}[mol/m

^{3}] is the total number of negative charges in a unit volume of equilibrium solution. In the pore water,

*X*

^{m−},

*Y*

^{n+}), one has

*O*(Π, Ω

^{l}) is given by

*O*(Π, Ω

^{l}) is a small quantity, which is negligible in Eq. [59]. For a dilute solution, , and . Hence, one obtains

This equation can be solved for and , provided that *c*_{0}, *c*_{fix}, and *n *^{l} are specified.

Remarkably, unless the soil has no fixed charges (i.e., *c*_{fix} is zero), is generally different than . This phenomenon is usually called Donnan’s effect (Mitchell and Soga, 2005), which can be attributed to the existence of the fixed charges in the porous media. The example discussed here represents a special case, in which only two solutes (i.e., ions *X *^{m−} and *Y *^{n+}) exist in the water solution and the soil (at Point B) is fully saturated. In a more general case, where the water solution consists of multiple species, it can be expected that the mass concentrations of species are different in the pores than in the well; that is, .

Equation [63] is a general expression for the osmotic pressure in porous media, implying that Π has two contributions, namely, the Donnan osmotic pressure (Π_{D}) and the pressure induced by the surface forces . Figure 2 depicts the dependence of Π_{D} on *c*_{0}, *c*_{fix}, and *n*^{l} in a clayey soil saturated with a NaCl solution (i.e., *i* = H_{2}O, Na^{+}, Cl^{−}). It can be seen that Π_{D} could be significant, depending on *c*_{0}, *c*_{fix}, and *n*^{l}.

*x*

_{C}(Fig. 1). At equilibrium, , where assumes a form similar to ,

*x*

_{0}is the vertical coordinate of the groundwater table, and

*b*the gravitational acceleration. Now it is straightforward to prove that the true pore water pressure at point C is given by . Because point C is arbitrarily chosen, one obtains a general expression for the true pore water pressure as

Clearly, *p*^{l} differs from by amount of Π. If a conventional transducer (say a tensiometer) is used to measure the “pore water pressure” in a soil, the measured value is indeed equal to . In the following, is called the “equilibrium solution pressure” to highlight the fact that can be measured using a conventional transducer. Remarkably, both *p*^{l} and are mechanical pressures, and thus they can be used to describe the mechanical behavior of unsaturated soils. In constitutive modeling, however, it is critical to distinguish which pressure is used. Unfortunately, such an important issue traditionally has been overlooked.

### Determination of Surface Energy Potential

*A*

^{g}is generally small, and the gaseous phase behaves as an ideal mixed gas. Hence,

*A*

^{g}is assumed to be a function of

*T*,

*p*

^{g}, and

*C*

^{g}only. According to Eq. [23],

*p*

^{g}

*− p*

^{s}= 0, and = 0. It follows that

*n*

^{l}to , one obtains

Noticeably, while *s*_{M} is a function of *T* and *n*^{l}, Π_{D} can be considered as a function of *T* and *n*^{l} only to the extent that both *c*_{0} and *c*_{fix} are specified. In general, Π_{D} is a function of *T*, *n*^{l}, and , which has yet to be specified. Such a development goes beyond the scope of this paper and may become a topic for future research.

Here it is recalled that porosity *n* is equal to . Equation [69] is the evolution equation for Ω^{l} in fully saturated *deforming* soils. At infinitesimal deformation, the derivative term is negligible, so that . Equation [63] yields Π(*T*, *n*) = *s*_{M}(*T*, *n*). According to the very definition of the air entry value (AEV) of a soil, *s*_{M}(*T*, *n*) simply equals to the AEV. Hence, at full saturation, the generalized osmotic pressure is equal to the air entry value of the soil.

### Measurement of Negative Pore Water Pressure

In the following, three typical types of the techniques for measuring the negative equilibrium solution pressure (i.e., ) are introduced, including the tensiometer, the vapor equilibrium technique, and the axis-translation technique. Briefly discussed here are the working principles of these techniques, and for a detailed analysis of their advantages and limitations, interested readers can refer to related references (e.g., Ridley and Wray, 1996; Lu and Likos, 2004).

#### Tensiometer

A tensiometer consists of a ceramic filter with a high air entry value, a strain gauged diaphragm, and a small reservoir. The reservoir is located behind the ceramics and at the front of the strain gauged diaphragm. Before measurement, the ceramics and the reservoir have to be fully saturated. Then the tensiometer is inserted into the point of interest in the soil, at which the pore water contacts with the water in the reservoir of the tensiometer (i.e., the equilibrium solution) through the ceramics. Measurement begins after equilibrium is achieved.

The working principle of a tensiometer is quite similar to that of a traditional PWP. Unlike the traditional PWP transducer, however, a tensiometer includes a high air-entry value ceramics, which is used to prevent the air enters into the water reservoir in the tensiometer during the measurement. The pressure measured by a tensiometer is indeed , as discussed above. The true pore water pressure can be recovered using Eq. [65]. Due to the water cavitation, a tensiometer can only be applied to measure a negative pore water pressure down to −100 kPa.

#### Axis-Translation Technique

To prohibit the water cavitation, the so-called axis-translation technique can be applied. Consider a representative volume of the unsaturated soil exposed to an ambient air pressure *p *^{g}, at which the equilibrium solution pressure and the true pore water pressure are and *p*^{l}, respectively. To avoid the cavitation of pore water, *p*^{g} is increased to during the measurement. As a consequence, the true pore water pressure rises to . Now consider a reservoir containing the soil solution in equilibrium with the pore water, and the solution pressure in the reservoir (denoted as ) remains constant during the measurement. For convenience, can be chosen as the atmospheric pressure, *p*_{atm}.

*p*

^{l}) is larger than by amount of Π. If the cavitation occurs at

*p*

^{l}= −100 kPa, then = −(Π + 100) kPa. Hence, the range of that ensures the validity of the axis-translation technique is given by

Based on the experimental data compiled in Baker and Frydman (2009), the cavitation tensions of many unsaturated soils, which is represented by the value of at cavitation, are all significantly higher than 100 kPa. For example, the cavitation tension of Barcelona silt (Gens et al., 1995) is up to 1 to 2 MPa. Hence, it is suggested here that the axis-translation technique should be applicable when is significantly lower than −100 kPa, and perhaps even in the whole range of matric suction in which this technique is commonly applied.

#### Vapor Equilibrium Technique

At low saturation, the aqueous phase in the pores is disconnected, and the water is adsorbed onto the grain surfaces, forming water films. To measure , the vapor equilibrium technique can be applied, in which the relative humidity of the pore gas is controlled and monitored.

*B*, while

*T*and

*p*

^{g}, quantity is indeed the relative humidity of the pore gas in the soil, which can be measured using either a psychrometer or filter papers (e.g., Lu and Likos, 2004). At the reference state, . Hence, it follows that

*T*and

*p*

^{g}as in the pores. Noticeably, determination of requires the extraction of equilibrium solution from the soil.

## Formulation of the Effective Stress Tensor

### Intergranular and Effective Stresses

*p*

_{eq}is the equivalent pore pressure, given by

*p*

^{g}=

*p*

^{s}and . Hence, it follows from Eq. [76] that

*p*

^{l}is given by Eq. [65]; that is, . Substituting in Eq. [78], one can derive

Note that, for a dilute aqueous solution, . Because Π_{D} is a function of *T* and *n*^{l} to the extent that the mass fractions of species (or, both *c*_{0} and *c*_{fix}) are specified, *p*_{eq} is generally a function of *p*^{g}, *T*, and *n*^{l}, as well as (*k* = 1, 2, ..., *Z* − 1).

*p*

_{eq}=

*p*

^{g}+ σ

_{s}, where σ

_{s}is the so-called suction stress (Lu and Likos, 2006), and defined by

_{D}− ρ

^{l}Ω

^{l}=

*s*

_{M}= AEV and , so that

_{D}, and Ω

^{l}are functions of

*T*,

*n*, and , as discussed above. From the above derivations, it is clear that Eq. [81] can smoothly convert to Eq. [82] as the soil switches from a partially saturated state to full saturation.

**σ′**for the fully saturated soils is defined as . Clearly,

**σ**′ can be related to

**σ**

_{I}

**′**through

In the current practice of soil mechanics, Π is usually neglected, so that the effective stress and intergranular stress are used interchangeably. However, such a simplification is not always permissible, since the bracketed term or the generalized osmotic pressure in Eq. [83] may vary significantly when the physicochemical effect comes into play.

### Evaluation of Pressure or Potential Functions

To apply the proposed theory, several pressure or potential variables have yet to be determined, including matric suction (or soil water characteristic function) *s*_{M}, surface energy potential function Ω^{l}, Donnan osmotic pressure Π_{D} and generalized osmotic pressure Π. In general, these variables are functions of temperature, pore water content (or saturation), and concentrations. Because of Eq. [63] and [68], only two of these variables are independent. Suction can be routinely measured (see above sections), while Function can be determined indirectly. One of the indirect methods for determining Π is to measure the shear strength of soils, incorporating the Mohr–Column criterion with Eq. [81] (e.g., Lu and Likos, 2006), from which Π_{D} and Ω^{l} are determined with *s*_{M} being given.

For some cases, one needs to measure *s*_{M} only, and other variables can be calculated using Eq. [63], [64], and [68]. As an example, consider a drying process of a water-saturated clayey soil through vaporization, whose *s*_{M} has been determined. Initially, the soil is fully saturated by a dilute NaCl solution, with a concentration corresponding to *c*_{0} = 100 mol/m^{3}. The soil has a porosity of 0.6, and a fixed charged density of 30 mol/m^{3}. Now the soil is slowly dried through vaporization (by reducing the humidity of the pore air) at room temperature (∼296 K) and atmospheric pressure, without any applied loads. As a first approximation, it is assumed that only the solvent (H_{2}O) in the pore solution can vaporize, and other species (Na^{+} and Cl^{−}) remain in the pores. As such, one can calculate the variation of Π_{D} with saturation by following the procedure stated in the “Donnan Osmotic Phenomenon” section above. Variables Π and Ω^{l} are then calculated using Eq. [63] and [68]. Because the initial state is at full saturation, is conveniently chosen as zero.

The variations of Π, *s*_{M}, Π_{D}, and Ω^{l} with the liquid saturation are depicted in Fig. 3. It can be seen that both Π and Π_{D} increase with the decrease of saturation. In particular, as the saturation decreases, Ω^{l} increases, implying that it becomes more and more difficult to remove the pore water, as is commonly observed in the laboratory. Figure 4 illustrates the variation of the mean intergranular stress, , during the drying process. Clearly, as the degree of saturation decreases, the intergranular *pressure* (i.e., −*p*_{intGR}) increases, and the soil becomes stronger and stronger. Noticeably, the variation of *p*_{intGR} depends also on the fixed charge density and the initial concentration of pore water. As an illustration, the results of *c*_{0} = 10 and 500 mol/m^{3} are also given in Fig. 4 (other conditions remain unchanged). Unlike those at lower concentrations, for *c*_{0} = 500 mol/m^{3}, the intergranular *pressure* first increases and then decreases. These results are consistent with those of the suction stress concept (Lu and Likos, 2006; Lu et al., 2010). Indeed, it has been well recognized that, during a drying process, some clayey soils may shrink first and become more and more aggregated at low saturation.

## Theoretical Framework

In general, several coupled physical and chemical processes may occur in an unsaturated soil, including heat conduction, fluid flow, capillary relaxation, diffusion, phase transition and chemical reactions as well as irreversible deformation.

### Heat Conduction Equation

*T*Λ_{M} is the rate of heat release associated with the physical and chemical processes occurring in porous media. Equation [84] takes the same form as derived by Coussy (2004). Noticeably, the dissipation associated with irreversible skeletal deformation is excluded in *T*Λ_{M}, since only the elastic deformation is considered here for clarity. Within this framework, however, the irreversible deformation can be incorporated in a straightforward way.

## Pore Fluid Flow Equations

^{l}is a function of

*T*and

*n*

^{l}, one can cast Eq. [86] into

Equation [88] and [89] are the general equations governing the pore fluid transport in the soils. Traditionally, the term ∇Π_{D} and ∇Ω^{l} have been neglected, while and *p*^{l} are used indiscriminately. As implied by Eq. [88] and [89], however, such a practice must be exercised with caution.

In fully saturated porous media, Π_{D} depends on *T* and (through *c*_{0} and *c*_{fix}). Hence Eq. [89] implies that if are heterogeneously distributed, there should exist a threshold of , under which no filtration can occur. Indeed, in carefully prepared samples in which the fluid concentrations are uniformly distributed, such a threshold phenomenon does not occur (Mitchell and Soga, 2005).

*x*

_{0}(the water table) to

*x*

_{LL’}, where x

_{LL’}is the vertical coordinate of the plane separating the fully saturated and unsaturated zones (Fig. 1), one obtains the height of the capillary fringe as

*LL*′ and equals the AEV. As a consequence, , that is, the true pore water pressure at

*x*

_{LL’}is zero.

The steady-state pressure profiles in the vadose zone deduced from Eq. [88] and [89] are illustrated in Fig. 1, which highlights the difference between *p*^{l} and . Remarkably, the profile is different from the hydrostatic pressure profile, which is the case, since the gradients of *T* or generally exist in the vadose zone.

## Capillary Relaxation

Equation [90] and [91] is the evolution equation of *n*^{l} that can be used to address the dynamic effect of capillarity (Hassanizadeh et al., 2002; O’Carroll et al., 2005; Manthey et al., 2008). It has been shown that the dynamic effect of capillarity can be attributed to local fluid flow in locally heterogeneous porous media (Wei and Muraleetharan, 2006, 2007).

### Species Diffusion

### Mass Exchange and Chemical Reactions

*J*; Γ

^{J}is a coefficient accounting for the rate of reaction; is the chemical affinity of the reaction, given by

^{sα}= 1 if α =

*s,*and otherwise δ

^{sα}= 0.

In soils, the solid grains are generally coated by the aqueous phase (i.e., the pore water). Therefore, the reactions of main concern here may include vaporization–condensation between the pore water and pore gas, pressure solution–precipitation between the solid grains and the pore water, and the chemical reactions in pore fluids.

#### Vaporization–Condensation Process

*T*and

*p*

^{g}(as in the pore gas). At equilibrium,

#### Pressure Solution–Precipitation

As an instance, the quartz dissolution is represented by .

*n*

_{r}), defined by Eq. [21]; and represents the deviations of the chemical potentials from their equilibrium counterparts.

*p*

_{I}(= −(

**τ**

^{s}:

**1**)/3) is the intergranular pressure. Noting that Δε

_{v}= Δ

*p*

_{I}/

*K*, where

*K*is the secant bulk modulus of the solid matrix, Eq. [101] can be approximated as

As implied by Eq. [103], a dissolving–precipitating process in the pores can be intervened by the intergranular pressure applied to the porous media. In general, however, the applied pressure and temperature do not exert simple controls on the pressure solution, since is also involved in. For dissolution, , whereas for precipitation. In the early beginning of the dissolution, few chemical bonds break, and thus . From Eq. [103], one can show that for the dissolution to occur, it requires that *p*_{I} > *K*. Hence, there exists a threshold of the applied pressure, below which no pressure dissolution can occur. Clearly, as the pressure increases, more and more chemical bonds will break, i.e., becomes more and more negative, and the dissolution will cease at . These theoretical results are consistent with experimental observations (e.g., Taron and Elsworth, 2010).

Variable is defined by Eq. [21]. Due to the current shortage of experimental data, it is impossible to develop here a realistic constitutive equation for . As a first approximation, one can assume based on the linearization of Eq. [21], where Θ is a positive material parameter. Provided that the mass exchange rates are given by Eq. [18], can be determined by integrating Eq. [1].

## Chemical Reactions in a Pore Fluid

*l*

_{i}.

Noticeably, for a specified , the first term in the right-hand side of Eq. [106] could be positive, zero, or negative, depending on the specific type of chemical reaction. Because , one can conclude that pressurization and the surface forces in the pores could have different controls on the chemical reactions in porous media.

## Summary and Conclusions

Developed here is a continuum theory of unsaturated soils, in which the soils are viewed as the porous media composed of multiphase and multispecies. Within the proposed framework, various coupling processes in unsaturated soils, including skeletal deformation, heat conduction, pore fluid flow, diffusion, phase transition, chemical reactions, and capillary relaxation, are addressed in a systematic and consistent way. The conditions for the thermal, mechanical, and chemical equilibrium of unsaturated soils are established, and the driving forces for various dissipative processes are identified and characterized. A general expression is developed for the chemical potential of a species in a fluid phase, which takes into account the thermal, pressure, osmotic, capillary, adsorptive, and electrostatic effects of unsaturated soils.

The behavior of unsaturated soils at equilibrium is investigated by enforcing the equilibrium conditions with incorporation of the proposed chemical potential formula. In particular, the Donnan osmotic phenomenon is discussed in detail, providing deep insights into the conception, composition, and measurement of pore water pressure, and a general formulation is developed for calculating the osmotic pressure. It is found that the true pore water pressure has three contributing components, including the equilibrium solution pressure, Donnan osmotic pressure, and the pressure induced by the surface forces. While the first component can be measured based on the conventional methods, the last two components cannot be directly measured. In addition, the experimental techniques used to measure the equilibrium solution pressure are discussed, and it is theoretically shown that the axis-translation technique can be applied, without water cavitation, down to a negative pressure much lower than expected.

Formulations of effective and intergranular stress tensors are developed for both fully saturated and unsaturated soils. It is suggested that it is the intergranular stress, instead of the effective stress, that controls the chemomechanical behaviors of both unsaturated and saturated soils. The new theory provides solid theoretical evidence for the recent development of the suction stress concept. In particular, a closed-form equation is developed for the suction stress, which can be used to describe the chemomechanical behavior of unsaturated soils.

A complete set of general governing equations are developed for addressing the chemomechanical behavior of unsaturated porous media in general, and unsaturated soils in particular, with a saturation ranging from extremely low to 100%. Remarkably, it is shown that other than the soil water characteristic function, another pressure (or potential) function Π(*T* ,*n*^{l} , (or Π_{D} or Ω^{l}) is generally required in analyzing the coupled chemomechanical processes in unsaturated soils. This conclusion is consistent with the recent development of the suction stress concept. The new theory is capable of explaining many salient phenomena occurring in water-saturated porous media, including Donnan osmosis, capillary fringe, air entry value, initial hydraulic head during seepage, pressure solution, and more.

This research was funded by one of National Basic Research Programs of China under Grant 2012CB026102, the National Science Foundation of China (NSFC) under Grants 51239010 and 11372078, and the Natural Science Foundation of Guangxi under Grant 2011 GXNSFE018004.

### Appendix

#### Notation

#### Balance Equations and the Second Law

The balance equations and the second law of thermodynamics, originally developed by Hassanizadeh and Gray (1979a,b) and Gray and Hassanizadeh (1989) for multiphase porous media based on a local averaging procedure, and later generalized by Bennethum and Cushman (1996a; 2002a) to multispecies porous media, are introduced here. In the following derivations, no thermodynamic properties are assigned to the interfaces between bulk phases, and thus only the balance equations of species and bulk phases are introduced.

##### Mass Balance

*i*(= 1, 2, ...,

*Z*) yields the mass balance equation for a bulk phase; that is,

**v**

^{α}is the velocity the mass center of the α phase, defined by

##### Linear Momentum Balance

*i*(= 1, 2, ...,

*Z*) yields the linear momentum balance equation for a bulk phase; that is,

**t**

^{α}, the Cauchy stress tensor of α phase, is given by

##### Energy Balance

##### The Second Law of Thermodynamics

Remarkably, *A*^{α} includes a contribution due to the diffusion of the species.

#### Additive Decomposition of Porosity

#### Constraints of the Second Law

*T*Λ, which is a linear function of D

^{α}

*T/*D

*T*, ∇

**v**

^{α}, ∇

*T*,

**v**

^{f,s}, , , and (

*i*= 1, 2, ...,

*Z*). These variables are not independent, because ’s and ’s are correlated. Namely,

_{i}and are the valence number and molar mass of charged species α

_{i}, respectively. ν

_{i}< 0 (or > 0) for an anion (or a cation), and ν

_{i}= 0 for a species with zero charge. Applying operator D

^{s}

*/*D

*t*to Eq. [A30], and introducing Eq. [A2], one can prove that

^{α}

*T/*D

*t*, ∇

**v**

^{α}, and can vary arbitrarily. Hence, for Inequality [A32] to be constantly satisfied, it requires that