# Vadose Zone Journal

- Copyright © by the Soil Science Society of America, Inc.

## 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.

In general, mass exchange may occur between the solid phase and the pore fluids. Thus, variables *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):where is the total mass exchange rate between the solid phase and pore fluid *f*.

Hence, it is proposed herein that

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

According to the second law of thermodynamics (see Appendix), it is necessary thatwhere η^{α} 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

According to the Bowen (1982) definition, the total chemical potential tensor of a bulk phase is given by

It is straightforward to prove that

Here, the diffusion terms including and have been neglected, since these terms are generally small compared to the magnitude of In addition, one has

Let be the chemical potential tensor of a species. can be related to through

According to the classical definition of electrochemical potential (Atkins and dePaula, 2002), it is immediately apparent that multiplier λ is equal to , where is Faraday’s constant, and ξ is the *local* electric potential. Dropping all the diffusion terms, one can cast Eq. [7] into

The residual dissipation inequality is given by Eq. [A46] in the Appendix. To address chemical reactions in the porous media, consider several stoicheiometric chemical reactions (including dissolution–precipitation between the two bulk phases) occur in the pores, through which the solid matrix and the pore liquid (water) exchange mass. These chemical reactions can be generally represented by , where [α_{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

Substituting these two equation in Eq. [A46], one obtains the residual inequality aswhere δ^{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

Inequality [19] includes a set of independent variables {*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.

Several common results can be deduced from Eq. [25]. First, for a general chemical reaction occurring in pore fluids, δ^{sf} = 0, and Eq. [2.25] yields

Now consider a species *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] thatand for (*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.

With Eq. [24], eliminating the exchange terms, , from the linear momentum balance equations of pore fluids (Eq. [A6] in the Appendix), one can derivewhere **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

From Eq. [26], it follows that (*k* = 1, 2, ..., *Z* − 1). Noting that **a*** *^{f} = 0 and ∇*T* = 0 at *static* equilibrium, one obtains

Integrating Eq. [32] yieldswhere *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

In analyzing the behavior of a porous medium, it is crucial to characterize its evolution process. It is assumed herein that during its evolution process, the porous medium deviates only slightly from thermodynamic equilibrium. With this assumption, and according to Inequality [19], it is assumed thatwhere **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).

Substituting Eq. [34–38] in Eq. [19] yields

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.

Now consider a system consisting of an infinitesimal amount of fluid in . Let the system be transported, in a thermodynamically reversible manner, from to the same sate as the REV. To maintain constant ρ* *^{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 obtainswhere *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

To derive a generic formulation of chemical potential, consider a species in the pore fluid, denoted as *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)Substituting Eq. [40] in Eq. [41], one obtainswhere is the chemical potential of a species in the bulk fluid free of surface forces, andAccording to its classical definition (Atkins and dePaula, 2002), chemical potential is explicitly given bywhere is the activity of the species that is a function of *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

Hence, the chemical potential can be expressed asand the electrochemical potential is given by

To gain more insights into the composition and character of Ω* *^{f}, recall the state equation for , i.e., Eq. [20]. Using Eq. [40], one obtains

Integrating Eq. [47] yieldswhere is the fluid volume fraction at full saturation, and for the wetting fluid, is approximately equal to porosity *n*; is considered as a function of *T* and *n *^{f}. Equation [48] can be recast intowherewhere 0^{+} 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.

Consider Point A in the equilibrium solution in the well, which is at the same elevation as Point B in the soil. The condition of electrical neutrality requires that in the equilibrium solution,where *c*_{0} [mol/m^{3}] is the total number of negative charges in a unit volume of equilibrium solution. In the pore water,

At equilibrium, according to Eq. [33],

According to Eq. [46], the chemical potential of the solvents in the well and the soil arewhere is the mass density of pure water, which is assumed to be constant; is the pressure of pure water at some reference state; is the chemical potential of pure water.

For ions γ (= *X *^{m−}, *Y *^{n+}), one haswhere and are the mass density and chemical potential, respectively, of species γ at the pure state. Here, is assumed to be constant. It is also noted that the local electric potential is zero in the well, where ions are uniformly distributed in the equilibrium solution.

Using Eq. [54], [57], and [58], it is straightforward to prove thatwhere is the generalized osmotic pressure, and *O*(Π, Ω^{l}) is given by

By substituting Eq. [55] and [56] in Eq. [54], one can derive

Comparing Eq. [60] and [61], one can see that *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, .

It follows from Eq. [61] thatwhere Π_{D} is the Donnan osmotic pressure, given by

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}.

Consider a point (denoted as C) in the vadose zone, whose vertical coordinate is *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 aswhere . In derivation, use is made of and .

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

For unsaturated soils, the effect of surface forces on *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

By substituting Eq. [65] in Eq. [66] and introducing Eq. [47] and [63], one deriveswhere is the matric suction or the soil water characteristic function. For a dilute solution, . By integrating Eq. [67] from *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.

It is remarkable that Eq. [67] also applies to the fully saturated condition, at which

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}.

According to Hilf (1956), if the saturation remains unchanged, the difference between the air pressure and pore water pressure is constant, i.e., . Hence, one has

The applicability of the axis-translation technique depends mainly on the validity of Hilf’s assumption as stated above. Application of this technique has been seriously criticized by Baker and Frydman (2009). These authors argued that the pore water in a natural soil should cavitate when the pore water pressure drops down to the cavitation pressure (about −100 kPa), and artificially elevating the pore air pressure through the axis-translation technique could intervene in this natural cavitation process, rendering a pseudo behavior of the unsaturated soil. It is noted, however, that according to Eq. [65], the true pore water pressure (*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.

At equilibrium, the chemical potential of the solvent is continuous across the interface between the pore water and the pore gas, i.e., . is given by Eq. [56], dropping subscript *B*, whilewhere is the partial vapor pressure in the pore gas. If takes the value of the saturated partial pressure of the vapor at *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

To determine the equilibrium solution pressure (i.e., ), one can substitute Eq. [65] into Eq. [73], and, after some rearrangements, obtainwhere is equal to the relative humidity, , of the mixed gas in equilibrium with the equilibrium solution under the same *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

To model the constitutive behavior of unsaturated porous media, it is crucial to define the effective stress tensor properly. To explore this issue, recall the definition of the total stress tensor for a porous medium as a whole (e.g., Bowen, 1982); that is,where all the diffusion velocity terms have been omitted.

In general, the effective stress can be defined only for the porous media at equilibrium. Combining state equations, Eq. [5] and [6], and equilibrium condition, Eq. [23], as well as Eq. [75], it is straightforward to prove thatwhere *p*_{eq} is the equivalent pore pressure, given by

As discussed in the previous sections, . As a consequence, *p*^{g} = *p*^{s} and . Hence, it follows from Eq. [76] thatwhere *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).

Equation [79] can be cast into the form as *p*_{eq} = *p*^{g} + σ_{s}, where σ_{s} is the so-called suction stress (Lu and Likos, 2006), and defined by

It follows from Eq. [76] thatwhere is the intergranular stress tensor (e.g., Mitchell and Soga, 2005, see Chapter 7), which is the part of the total stress transferred through grain–grain contacts.

At full saturation, Π = Π_{D} − ρ^{l}Ω^{l} = *s*_{M} = AEV and , so thatwhere Π, Π_{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.

In the classical soil mechanics textbook, Terzaghi’s effective stress tensor **σ′** 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

Using energy balance equations (Eq. [A10–A18] in the Appendix) and Eq. [34], one can derive the heat conduction equation aswhere , and

*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

Using Eq. [36] to eliminating from the linear momentum balance equation of a fluid (Eq. [A6] in the Appendix), one obtains

For the pore gas, and thus

For the pore water, using Eq. [47] and acknowledging that Ω^{l} is a function of *T* and *n*^{l}, one can cast Eq. [86] into

Noting that for a dilute solution, . Using Eq. [63] and [65], one obtains

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).

Interestingly, Eq. [89] predicts the existence of a capillary fringe immediately above the groundwater table in the vadose zone (Fig. 1), where the soil is fully saturated and is negative. To make this point transparent, consider the unsaturated soil at equilibrium under the atmospheric pressure and the isothermal condition, Eq. [89] becomes . Assuming and integrating this equation from *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 aswhere is the equilibrium solution pressure at plane *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

For the gaseous phase, . Using Eq. [35] and [47], one obtains

Assuming that Eq. [63] and [65] are valid in transient conditions, one can derive

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

Species diffusion occurs only in pore fluids. Multiplying Eq. [A2] by , and then subtracting the resultant from Eq. [A1], it follows after some algebraic manipulations that

Substituting Eq. [17], [18], and [38] in Eq. [92], one obtainswhere can be determined once if the chemical reaction is specified.

### Mass Exchange and Chemical Reactions

As discussed in earlier sections, the mass exchange processes among various bulk phases can be addressed through specifying several typical stoicheiometric chemical reactions. To derive a general equation governing the mass exchange and chemical reactions, we cast Eq. [37] into the Arrhenius’ type aswhere is the activation energy of reaction *J*; Γ^{J} is a coefficient accounting for the rate of reaction; is the chemical affinity of the reaction, given bywhere if ; δ^{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

For a vaporizing–condensing processes, it follows from Eq. [95] that

In a vaporizing process, , whereas in a condensing process. Using Eq. [56] and [72], it is straightforward to prove thatwhere is the relative humidity of the pore gas, and is the relative humidity of the equilibrium solution at *T* and *p*^{g} (as in the pore gas). At equilibrium,

Here, subscript 0 denotes the equilibrium state. Linearizing Eq. [97] and using Eq. [98], one obtains

#### Pressure Solution–Precipitation

Pressure dissolution and precipitation occur only between the solid matrix and the pore water in unsaturated soils. Without loss of generality, consider a group of chemical reaction processes, which can be represented by

As an instance, the quartz dissolution is represented by .

Accordingly, the chemical affinity is given bywhere is a function of (or *n*_{r}), defined by Eq. [21]; and represents the deviations of the chemical potentials from their equilibrium counterparts.

The chemical potential of a species in the solid phase is generally dominated by the volumetric strain energy of the solid skeleton. Thus, one can writewhere *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 aswhere all the terms including are dropped out, since they are generally much smaller than the bracketed term.

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

The chemical affinity of a chemical reaction in a pore fluid is

Particularly, in the pore water, the chemical potential is given by

During a chemical reaction, both the mass and the number of electric charges are conserved. Hence, substituting Eq. [105] in Eq. [104], one obtainswhere is the molar volume of species *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.

## Acknowledgments

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

The mass balance equation for a species is given by

Summing Eq. [A1] over *i* (= 1, 2, ..., *Z*) yields the mass balance equation for a bulk phase; that is,where , and **v**^{α} is the velocity the mass center of the α phase, defined by

Mass conservation requires that

#### Linear Momentum Balance

The linear momentum balance equation for a species is given by

Summing up Eq. [A5] over *i* (= 1, 2, ..., *Z*) yields the linear momentum balance equation for a bulk phase; that is,where **t**^{α}, the Cauchy stress tensor of α phase, is given by is the diffusion velocity. Equation [A5] and [A6] are restricted by

#### Energy Balance

The energy balance equation for a species is given by

Summing up Eq. [A10] yields the energy balance equation for a bulk phasewhere

Equation [A10] and [A11] are restricted by

#### The Second Law of Thermodynamics

The second law of thermodynamics requires that the rate of total entropy production Λ in the unsaturated soils due to dissipation be nonnegative; that is,

The entropy balance requires that

It has been well recognized that, in deriving a general thermomechanical theory of porous media, it is more convenient to use the Helmholtz free energy instead of the internal energy function. By Legendre transformation, the Helmholtz free energy density for a species and a bulk phase are defined, respectively, as

Using Eq. [A12] and noting that , one can derive

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

Replacing in Eq [A19] by and *A*^{α}, one can derive the following expression for the second law,

### Additive Decomposition of Porosity

Introducing the mass balance equation of the solid phase (Eq. [A2]), one can derive

It is clear that the variation of porosity includes two contributions, which are induced by mechanical actions (including solid material compression and skeletal deformation) and mass exchanges, respectively. The variation of the porosity due to chemical reactions is described bywhere the subscript “r” denotes the chemical reactions.

### Constraints of the Second Law

Substituting Eq. [2] and [3] into Inequality [A25] can yield an expression of *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,

Taking the spatial gradient on both sides, one obtains

In addition, the species are electrically charged. The electrical neutrality requires that the total number of charges in the porous media be zero. According to this requirement, one haswhere ν_{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

The constraints represented by Eq. [A28], [A29], and [A31] can be removed by using the Lagrange-multiplier method. To this end, the right-hand sides of Eq. [A28], [A29], and [A31] are first multiplied by multipliers , , and λ, respectively, and then added into the right-hand side of Eq. [A25]. Finally, introducing Eq. [2] and [3], one can derivewhere λ is a scalar, a vector, a symmetric second-order tensor, and

In a general thermodynamic process, variables D^{α}*T/*D*t*, ∇**v**^{α}, and can vary arbitrarily. Hence, for Inequality [A32] to be constantly satisfied, it requires that

It follows from Eq. [A43] thatand from Eq. [A42] that

Using Eq. [A27] and Eq. [A39–A43], one obtains the residual dissipation inequality as

## Footnotes

All rights reserved. No part of this periodical may be reproduced or transmitted in any form or by any means, electronic or mechanical, including photocopying, recording, or any information storage and retrieval system, without permission in writing from the publisher.

- Received July 17, 2013.

Open access