For unsaturated sands and silts, spatial distributions of capillary water bridges and intergranular fabrics roughly coincide. Microscale capillary-suction-induced interparticle forces lead to a macroscale capillary-suction-induced intergranular stress tensor with similar spatial characteristics. Its physical parameters follow from a macroscale experimental database by using inverse modeling.

## Abstract

For the saturated case with only one pore fluid, either water or air, the roles of both the intergranular stress tensor and the pore fluid stress can be distinguished easily. In the unsaturated case, the capillary water is recognized to induce capillary suction in the pores and capillary-suction-induced interparticle forces. At the macroscale, volume averaging of these forces would lead to the capillary-suction-induced intergranular stress tensor. In its approximate formulation, the concept of the fabric stress tensor is applied, enabling the effect of the spatial distribution of the intergranular fabric on the capillary water bridges as occurring in the drier pendular saturation phase to be accounted for. Subsequently, the combined intergranular stress tensor and the combined pore fluid stress tensor can be derived directly. The constitutive relation of a granular skeleton, composed of elastic particles with mainly frictional interaction, like quartz sands and silts, is considered to remain independent of the degree of saturation. Under such restrictive conditions, only the additional physical parameters of the capillary-suction-induced intergranular stress tensor need to be determined, which can be achieved by means of inverse modeling, taking advantage of all macroscale experimental data and physical modeling for the whole unsaturated range. For clays and peats, with potential physicochemical and biochemical actions and double porosity and/or fibrous microstructures, the constitutive models can be expected to be physically more complicated, thus involving more physically relevant parameters. Hence, clays and peats must be considered to fall outside the scope of the proposed model framework.

**The physical relevance** of the continuum mechanical measures of stress, deformation, and flow of the pore fluids and stress and deformation of the solid skeleton forms the basis of any constitutive modeling and subsequent application for predictions in geomechanics.

The physical relevance of the applied continuum measures of stress, deformation, and flow is a reflection of the physical concepts as applied in their descriptions. For instance, for the fluid-saturated case, the calculation of the deformation of the solid skeleton due to a changing pore fluid stress can be achieved by applying the macroscale isotropic pore fluid stress tensor *p***I**, at least if the substance composing the particles remains elastic. For the solid skeleton of granular materials, the intergranular stress tensor **σ*** and a potential microstructure tensor are derived using micromechanics in combination with volume averaging, irrespective of the degree of saturation. For the unsaturated case, the effects of the two simultaneous pore fluids on the solid skeleton are limited by the conditions for the granular skeleton that its deformation remains identically dependent on the tensors of the intergranular stress and a combined measure of both pore fluid stresses, as for both saturated cases. In fact, these conditions will be used for the quantification of the combined unsaturated measure of both pore fluid stresses. In addition, the potential of experimental methods for quantifying these measures may also help to appreciate their physical significance, irrespective of whether they are obtained directly by physical observation or by discrete element modeling.

Furthermore, the measures of stress, deformation, and flow must be consistent with the basic conservation laws of mass, momentum, and energy and the non-negative entropy production law of thermodynamics. The model formulation by Gray and Schrefler (2007) satisfies these conditions and implicitly derives the types of required experimental data needed to finally arrive at applications. However, while some experimental capabilities for solid geomaterials do reach the three-dimensional nanometer scale (e.g., Desbois et al., 2011), for loose granular materials and multiple pore fluids even the three-dimensional micrometer scale is still rather demanding. Most observations concern macroscale quantities, increasingly supplemented with some microscale data, which are widely used in phenomenological constitutive models in engineering attempts to capture the most essential phenomena as reproducible in element tests, such as oedometer and triaxial tests.

The main motivation of this study was to understand why and how the saturated stress measures can be extended to keep the unsaturated stress measures physically justified and the amendments to the constitutive relation minimal or none, thus enabling maximum advantage to be taken of the earlier findings on the study of the mechanical characteristics for the saturated case.

First, the saturated case of granular materials with only one pore fluid, albeit a liquid, e.g., water, or a gas, e.g., air, is reviewed. These two cases form the limits for the unsaturated case, in which the pore space is simultaneously filled by two fluids, namely pore water with dissolved air and pore air with water vapor. The microstructure and the measures of stress and deformation are reviewed, starting from a description at the microscale of a packing of spheres. The resulting macroscale quantities involve the intergranular fabric tensor, the intergranular stress tensor, the scalar pore fluid stress, and the constitutive model of the granular material, relating the material rates of stress and strain. Consequently, field measures of stress can be decomposed into the intergranular stress tensor field and the isotropic pore fluid stress field. Both stress measures occur in the constitutive relation for describing the deformation of the solid skeleton.

For the unsaturated case of granular materials then, first the concept of a capillary-suction-induced interparticle force is introduced, which is shown to depend on matric suction. Furthermore, considering that the microscale anisotropic spatial distribution of the water bridges at the interparticle contacts in the drier pendular saturation phase corresponds approximately to the intergranular fabric tensor, the expression of the capillary-suction-induced intergranular stress tensor can be derived.

The constitutive relation of a granular skeleton, composed of elastic (or rigid) particles with mainly frictional interaction, like quartz sand, is considered to remain independent of the degree of saturation. Under such restrictions for the unsaturated case, we demonstrate that the effects on the solid skeleton by both pore fluids can be accounted for appropriately through amendments to mainly both saturated stress measures, leading to a consistent framework for this major class of unsaturated granular materials. Nevertheless, this proposed approach does not imply that microscale observations could not further validate and/or strengthen its physical basis, in particular concerning the spatial distributions of the microstructures of both the granular skeleton and the pore water and their evolution.

However, for soils with more complicated microstructure characteristics and/or interactions at interparticle contacts, like clays, and/or the biochemistry of the organic fibers and bulk substance, like peats, the proposed approach is expected to not be suitable because the constitutive relations of the solid skeleton will also depend on the degree of saturation.

The main objective of this study was to show that the formulation with two stress tensors, the intergranular stress tensor **σ*** and the second-order pore fluid stress tensor **p**, is physically justified and applicable for irreversible saturated and unsaturated granular materials such as sands and silts.

We start with a review of the existing stress measures for saturated and unsaturated granular materials. Then the micromechanical characteristics of the solid skeleton of simplified granular materials as applied in the discrete element method are described and the corresponding macroscale measures as fabric and intergranular stress tensors are derived. Next, for unsaturated granular materials, a microscale description of capillary suction and capillary-suction-induced interparticle forces is provided, followed by the derivation of the corresponding macroscale intergranular and combined pore fluid stress tensors.

Finally, the intended way of application of the derived capillary-suction-induced intergranular and combined pore fluid stresses is discussed. For the short term, the application of the current understanding may allow further progress by applying inverse modeling for the determination of the described model characteristics of sands and silts as affected exclusively by capillary action. However, for clays and peats, additional non-capillary physical properties will have to be accounted for as well.

## Review of Stress Measures of Saturated Granular Materials

The sign convention of continuum mechanics is applied here, with tension and stretching being positive.

*p*and normal total stress σ, Fillunger (1915) reported that in unjacketed tensile tests the water pressure inside and outside of the pores had no influence on the strength of the solid skeleton. Terzaghi (1923, 1925 [p. 50–52], 1943) proposed, as the effective normal stress σ′ for the constitutive properties of a solid skeleton composed of incompressible particles, the normal stress difference σ′ = σ −

*p*, which expressed in tensor format leads to Terzaghi’s effective stress tensor:

being the difference between the total stress tensor **σ** and the isotropic pore fluid stress tensor *p***I**.

*A*, average force per intergranular contact

*P*, and average number of contacts per unit of cross-sectional area

*N*. Bishop indicated the effective contact area of the soil particles per unit cross-sectional area by contact fraction

*a*=

*AN*. Then the resulting average intergranular force per unit of cross-sectional area of the contact plane indicated by normal stress σ

_{i}′ is expressed by (Bishop, 1959, Eq. [4])

*p*over the whole of its surface, it undergoes no distortion but a small decrease in volume. Hence it is only that part of the local contact stress (

*P*/

*A*= σ

_{i}′/

*a*[authors’ addition]) which is in excess of

*p*that causes deformation of the soil structure.” This reasoning is supported by the decomposition of the tractions on a particle, as illustrated on the right side of Fig. 1. Then continuing: “By summing the corresponding components of the intergranular force per unit of cross-sectional area, [thus (σ

_{i}′/

*a*–

*p*)

*a*, authors’ addition], an expression is obtained for normal stress σ

_{c}′, defined as that part of the normal stress which controls volume change due to the deformation of the soil structure,” leading to (Bishop, 1959, Eq. [6])

_{c}′, the controlling volume change due to the deformation of the soil structure, in accordance with Bishop (1959, Eq. [7]):

On this basis, Bishop concluded that “…although the average intergranular force per unit area depends on the magnitude of (contact fraction) *a*, volume changes due to deformation of the soil structure depend simply on the stress difference σ − *p*, whatever the value of *a*.” Consequently, Bishop’s theoretical normal stress σ_{c}′ in Eq. [4] must be the normal component of the intergranular stress tensor because it controls the volume change due to the deformation of the soil structure. We conclude that Bishop’s theoretical intergranular stress tensor in the soil structure is identical to Terzaghi’s proposal of the effective stress tensor **σ′** in Eq. [1].

*C*and

*C*

_{s}are the compressibilities of the skeleton of the porous material and the solid substance comprising the particles, respectively. Note that in Eq. [5] we have replaced the small increments Δ by super dots as material time derivatives for allowing coefficients

*C*and

*C*

_{s}to represent tangential quantities, enabling the representation of severely nonlinear material characteristics. Subsequently, aiming for the definition of the corresponding effective stress (Skempton, 1960, sixth and seventh expressions on p. 12), Eq. [5] was rearranged as follows:

where η = 1 – *C*_{s}/*C* is Biot’s coefficient. Consequently, the effective stress according to Eq. [7] will be indicated here for convenience as Biot’s effective stress.

Lade and de Boer (1997) showed experimentally that in one-dimensional compression of quartz and gypsum sands up to extremely high pressures, the Biot’s coefficient η equals unity at zero pressure and decreases increasingly from unity with increasing pressure, as illustrated in Fig. 2 (based on Lade and de Boer, 1997, Fig. 14). For the calculation of Biot’s coefficient η, published data on the compressibility coefficients of quartz and gypsum were used. Lade and de Boer (1997) concluded on that basis that at higher stresses significant deviations from Terzaghi’s effective stress principle (Eq. [1]) occur.

Application of the Biot’s effective stress (Eq. [7]) is mainly appropriate if both the solid skeleton and the solid substance comprising the particles are linear elastic, which was also the basic assumption behind Eq. [5].

*C*in Eq. [5], and the often linear elastic properties of the substance composing the particles, like

*C*

_{s}in Eq. [5] to be simultaneously accounted for. Consequently for saturated granular materials, we propose to apply as complementary stress measures the isotropic pore fluid stress tensor

*p*

**I**and the intergranular stress tensor:

with normal components equal to Bishop’s stress σ_{c}′ in Eq. [4], thus controlling the volume change due to the deformation of the soil structure.

in which the mechanical characteristic of the solid skeleton with compressible particles, *C*{**σ***, ,*e*_{v},…}, depends at least on the intergranular stress tensor **σ***, the strain rate tensor , the instantaneous void ratio *e*_{v}, and possibly several other scalar and tensorial state parameters, while the mechanical characteristics of the substance comprising the particles, *C*_{s}, may often (but not necessarily) remain a simple constant scalar material parameter, except for large intergranular stresses inducing particle visco-plasticity and/or crushing.

**σ***(Eq. [9]) equals Terzaghi’s phenomenological proposal for an effective stress tensor

**σ′**(Eq. [1]) for incompressible particles. The differences between the intergranular and Terzaghi’s and Biot’s effective stress tensors occur as a consequence of the compressibility of the individual particles (Skempton, 1960), leading to Eq. [5] and [6]. A comprehensive demonstration of their difference is obtained by equating the first term of Eq. [6] and [10], giving

*C*

_{s}/

*C*of particle substance and solid skeleton:

From Eq. [12] it follows directly that both the rates of Biot’s effective stress and the intergranular stress become equal if the compressibility ratio *C*_{s}/*C* approaches zero, *C*_{s}/*C* → 0, thus when the compressibility of the particle substance *C*_{s} becomes negligible with respect to the compressibility *C* of the solid skeleton. Biot’s effective stress according to Skempton’s Eq. [7] and its alternative form Eq. [12], obtained after substituting Eq. [9], implicitly combines the effects of the intergranular and pore fluid stresses, while both right-side terms of Eq. [11] account for the volume change of the intergranular structure and the volume change of the substance composing the particles due to a change in the pore fluid stress, respectively.

## Review of Effective Stress for Unsaturated Granular Materials

*p*

^{w}in the liquid phase, which is lower than that in the gaseous phase

*p*

^{g}, will act only on a reduced area. This notion enabled Bishop to propose the following expression for the normal effective stress:

which is the sum of the net stress σ − *p*^{g} and the product of fraction χ and the matric suction *p*^{w} − *p*^{g}. This fraction χ equals unity for saturated soils and decreases with decreasing degree of saturation *S*^{w} to zero. Following this proposal of normal effective stress, its role in the constitutive modeling of the solid skeleton was disputed while distinguishing the effects of each stress measure occurring in Eq. [13], net stress and matric suction.

Bishop and Blight (1963) demonstrated the determination of the magnitudes of the χ parameter by combining measured shear strengths of saturated and unsaturated soils in terms of effective stresses. However, a severe shortcoming of the proposed effective stress definition σ′ was also recognized (e.g., Jennings and Burland, 1962; Burland, 1965), namely the impossibility of explaining pore collapse, which is understood to involve the sudden instability phenomenon of volumetric contraction during wetting by reduction in the matric suction.

*S*

^{w}, thus χ =

*S*

^{w}, has been found to be a first rough approximation for sands and to some extent also for silts but decreasingly applicable for clays with increasing plasticity index. The following proposal for fraction χ was investigated (Brooks and Corey, 1964; Khalili and Khabbaz, 1998):

which was applied for matric suctions *p*^{w} − *p*^{g} larger than the matric suction at air entry (*p*^{w} − *p*^{g})_{e}, for which two values were recognized, namely one for air entry during drying and the other for air expulsion during wetting. Khalili et al. (2004) demonstrated experimentally that with this approximation, good predictions of not only the shear strength and volume change of sands and silts are obtained but also of the critical state line, independent of the degree of saturation. Consequently, these results are in support of Bishop’s effective stress proposal (Eq. [13]), which significantly simplifies constitutive modeling.

*S*

_{r}

^{e}(Alonso et al., 2010, Eq. [12]) of these macropores, the microstructural threshold parameter

*S*

_{r}

^{m}is introduced, quantifying the lower limit of the degree of saturation

*S*

_{r}

^{e}= 0 of the macropores for the case in which the macropores are dry, containing only gas, while the intra-aggregate micropores are still pore water saturated. Following the definition of this degree of saturation

*S*

_{r}

^{e}of the macropores, subsequently fraction χ is assumed equal to the degree of saturation

*S*

_{r}

^{e}of the macropores:

Experimental evidence, involving the strength and elastic stiffness of the solid skeleton of high-plasticity clays, has demonstrated that the single macroscale Bishop type of effective stress measure (Eq. [13]) is still very useful for engineering practice if the additional microscale characteristic as expressed by the microstructural threshold parameter *S*_{r}^{m} is accounted for, too.

## Microstructure, Stress, and Deformation of Granular Packing of Spheres

For the development of macroscale stress measures for unsaturated granular materials, first we review the microscale descriptions of the structures of the solid skeleton and the pores, the interparticle forces, and the capillary microscale stresses in the pore water and pore gas. The major part of these descriptions have been developed and applied for the numerical simulation of the behavior of granular materials by means of the discrete element method (DEM) (e.g., Cundall and Strack, 1979). The analytical and numerical investigations have already increased insight into the roles of various microscale characteristics and the corresponding macroscale quantities, such as the evolution of the fabric tensor, intergranular stress, and deformation (e.g., Thornton, 2000; Rothenburg and Kruyt, 2004; Sun and Sundaresan, 2010).

*r*on each side of the interparticle contact point indicated by the sequence number

*k*for

*k*= (2

*r*– 1,2

*r*), respectively. The outward normal unit contact vectors of the material on both sides are indicated by

**n**

^{(2r−1)}and

**n**

^{(2r)}and the corresponding radius vectors by

**R**

^{(2r−1)}and

**R**

^{(2r)}, while the opposing interparticle contact force vectors are

**f**

^{(2r−1)}and

**f**

^{(2r)}, respectively. The microstructure of the packing of spheres can be quantified in terms of the fabric tensor

**A**(Satake, 1978, 1982, 2004; Kanatani, 1984; Oda et al., 1985; Subhash et al., 1991; Li et al., 2009):

*C*= max(

*r*) represents the number of direct contacts of the packing. Then the corresponding average intergranular stress tensor

**σ***is expressed by (Goddard, 1977; Rothenburg and Selvadurai, 1981; Christoffersen et al., 1981; Nemat-Nasser, 1982; Thornton and Barnes, 1982, 1986; Kruyt and Rothenburg, 1996; Thornton, 2000)

where *V* is the sphere packing volume. For the saturated case, the pore space is filled by a pore fluid with isotropic stress tensor *p***I**.

The effects of the changes in both stress measures are illustrated in Fig. 4. The left-hand side illustrates the densification of the particles due to the increase in both the pore fluid pressure and the isotropic part of the intergranular stress tensor **σ***, while on the right-hand side the interparticle shearing due to an increase of the deviatoric part of the intergranular stress tensor **σ*** is depicted. It may be noted that Fig. 4 has some similarity with Fig. 1, namely that the pore fluid pressure acting across the whole surface of each particle, indicated in the middle sketch of Fig. 1, also contributes to the densification of each particle, as indicated on the left side of Fig. 4. Nevertheless, each particle is also densified due to an increase in the isotropic pressure of the intergranular stress tensor **σ***, represented above by Bishop’s intergranular normal stress σ_{c}′ (Eq. [4]).

where the contribution by the objective material rate of the intergranular stress tensor **σ*** (Zaremba, 1903; Jaumann, 1906), indicated by superscript *J*, is represented by the double dot product with the possibly anisotropic nonlinear elastic visco-plastic fourth-order tangent stiffness tensor function **D** of the solid skeleton.

We consider this constitutive fourth-order tensor function **D** to be limited to a dependence on scalar and tensorial state parameters of all components of the solid skeleton exclusively, thus including the bulk and shear stiffness of the elastic particles, and to be formulated in an observer-independent form. This macroscale fourth-order intergranular tensor function could be partly phenomenological, based on experimental macroscale observations, when embedded within an appropriate thermodynamic framework. However, the constitutive tensor function could, in principle, also be derived by volume averaging of a microscale model, e.g., in the form of a DEM, with appropriate evolution terms of the state parameters, if the resulting macroscale behavior would correspond sufficiently to micro-macroscale observations.

The contribution in Eq. [18] by the material rate of the isotropic pore fluid stress tensor *p***I** depends on the tangent bulk stiffness *K*^{sf} of the material comprising the solid particles. In the case of elastic particles, the right-side terms of Eq. [18] may be rearranged to arrive at Biot’s effective stress tensor , occurring in the rate relation of stress and strain, namely = **D**^{−1}: (Biot, 1941; Verruijt, 1982, 1984, 2010), which for the simple volumetric strain rate model Eq. [5] as described by Skempton (1960) would reduce to _{vol} = C in Eq. [6]. Biot’s effective stress tensor would approach the intergranular stress tensor **σ*** in Eq. [9] when the tangent elastic bulk stiffness *K*^{sf} of the substance comprising the particles would approach infinity with respect to the equivalent tangent bulk modulus of the solid skeleton *K* = (**I**:**D**:**I**)/9, thus if *K*/*K*^{sf} → 0, which was concluded already above when considering the compressibility ratio *C*_{s}/*C* → 0 following Eq. [12]. It should be noted that for mineral granular materials like quartz sands and silts and for mean intergranular stress levels occurring in geotechnical engineering applications, the assumption of elasticity may be rather realistic. However, the effects of grain shape and roughness of the particle surface may lead to additional micromechanical parameters like the average particle orientation for best-fitted ellipsoidal particles. For the description of the microstructure of clays, composed of plate-like particles, more complicated fabric tensors will be needed, while the interaction at the macroscale may become elasto-visco-frictional, based on the combined mechanical and chemical interactions at the microscale particle surface.

For geomaterials like fibrous peat, a macroscale multicomponent organic material model may be needed, distinguishing between organic fibers with a spatial distribution and a homogeneous organic bulk, each of which may require an elasto-visco-plastic macromodel.

**N**

^{(2r)}=

**n**

^{(2r)}

*N*

^{(2r)}and tangential component vector

**T**

^{(2r)}=

**m**

^{(2r)}

*T*

^{(2r)}of the interparticle contact force vector

**f**

^{(2r)}, in which the unit vector

**m**

^{(2r)}indicates the direction of

**T**

^{(2r)}. Substitution of the normal and tangential component vectors enables the intergranular stress tensor

**σ***in Eq. [17] to be decomposed into the intergranular stress tensors

**σ**

^{N}and

**σ**

^{T}due to the normal and tangential interparticle contact force vectors, respectively, namely

**σ**

^{N}due to the normal interparticle contact force vectors

**N**

^{(2r)}(Thornton, 2000), adopting the DEM, introduced the concept of the fabric stress tensor

**σ**

^{A}, which is defined in the second equation, namely

and **A** is the fabric tensor (Eq. [16]) of the packing of spheres. From Eq. [20] follows that the fabric stress tensor **σ**^{A} is the product of the mean normal interparticle forces (Eq. [21]) and the fabric tensor **A** (Eq. [16]), thus involving the same quantities as defining the intergranular stress tensor **σ**^{N} due to the normal interparticle contact force vectors but without the physical correlations embedded in **σ**^{N}.

Thornton (2000) demonstrated the physical relevance of each quantity by means of Fig. 6, illustrating various types of deviatoric intergranular stress states in the π plane in the principle stress space at peak stress for purely deviatoric loading of an initially dense dry packing of spheres. In Fig. 6, the larger locus illustrates the deviatoric peak strength. The locus by the deviatoric part of the intergranular stress tensor **σ**^{N} (defined according to the first expression of Eq. [20]), due to the normal interparticle contact force vectors, is shown to be only marginally smaller than the outer locus. It can be concluded that the contribution by **σ**^{T} due to the tangential interparticle forces **T**, represented by the difference between the two larger loci, remains relatively small. Furthermore, Fig. 6 demonstrates the significance of the intergranular fabric in terms of the deviatoric part of the fabric stress tensor **σ**^{A} (Eq. [20]), represented by the size of the inner locus, involving both the fabric tensor **A** (Eq. [16]) and the mean scalar measure (Eq. [21]) due to the average normal interparticle force, which is already about 33% of the outer locus.

Figure 7 depicts the corresponding evolution of the deviatoric part of the fabric tensor **A** as a function of the deviatoric strain measure ε_{a} − ε_{r}. The deviatoric fabric–strain curve shows a peak of about *A*_{a} – *A*_{r} ≈ 0.12 at a similar strain ε_{a} − ε_{r} ≈ 0.07 as that at which the deviatoric peak stress is occurring. At the critical state, the deviatoric fabric value reduces to *A*_{a} – *A*_{r} ≈ 0.085. This pattern is similar to the corresponding deviatoric stress–strain behavior shown in Fig. 8. The illustrated decomposition of the deviatoric stress into the contributions by the normal and tangential interparticle force components demonstrates the significance of the normal components as depicted in Fig. 6.

## Suction-Induced Interparticle Forces

In unsaturated granular materials for low degrees of saturation, the air phase is continuous while the liquid pore water occurs both as water bridges between neighboring particles, illustrated in Fig. 9a, and as isolated pockets of adhesive water at the particle surfaces. This saturation phase is known as the *pendular phase* (Bear, 1972). Figure 9b illustrates the detailed actions by a water bridge between two rough approximately spherical particles with equal radii *R* in direct contact. These actions are all induced by the surface tension σ along the interface between the water and air. One major consequence of this surface tension σ is the pore water suction *p*^{w} occurring in the pore water composing the water bridge, while in the surrounding air the gas stress *p*^{g} occurs. A second consequence of the surface tension σ is the line force acting along the intersection curve of the water surface and the surface of a particle. The orientation of the water surface and the coinciding surface tension σ along the intersection curve, indicated by contact angle θ, depends on the direction of the motion of the intersection curve along the particle surface, either wetting or drying. Both the induced pore water suction *p*^{w} and surface tension σ acting along the intersection curve pull the particles toward each other and as such, due to equilibrium at the interparticle contact point, cause the suction-induced interparticle force vector **f** as a reaction. We note that, due to the surface roughness *s* of each particle surface, the approximately representative spheres in contact are kept at a relative distance of 2*s*.

The resulting interparticle force vector **f** can be expressed in dimensionless form as illustrated in Fig. 10, namely in terms of **F** = **f**/(σ*R*) (Molenkamp and Nazemi, 2003a) as a function of dimensionless capillary (matric) suction Ψ = ψ*R*/σ = (*p*^{w} − *p*^{g})*R*/σ and dimensionless surface roughness *S* = *s*/*R*, while still depending also on the contact angle θ.

The curve for smooth particles with zero dimensionless surface roughness, *S* = 0, seems to correspond to the observed behavior of, e.g., soft clays during severe drought, where during continued drying, solid lumps of clay particles are formed, separated by a developing pattern of cracks due to the shrinking of the lumps. On the other hand, drying of initially saturated sand leads finally to dry sand grains, such as occur higher up sandy beaches, without any suction-induced interparticle forces. Consequently, the curves for the rough particles *S* ≥ 0.001 in Fig. 10 with the dimensionless interparticle force **F** approaching zero seem to represent the characteristics of sand.

Despite the apparent suitability of the role of surface roughness, we believe that the lack of actual physics of surface phenomena at the microscale, occurring at interparticle contacts during drying of sands and clays, warrants further microscale investigation.

## Capillary-Suction-Induced Intergranular Stress

^{x}is provisionally defined by

in which the generalized pore air stress *p*^{x} is introduced for all phases of unsaturation. Subsequently for the phases of unsaturation, Cases a and b are distinguished:

*p*

^{g}occurs in the continuous air phase throughout the porous network and is usually in open contact with the atmosphere:

*insular saturation phase*(Bear, 1972). The pore air stress

*p*

^{b}is the background pore air stress, with a virtual magnitude that would occur if the local pore air would be in open hydrostatic equilibrium with the atmosphere:

In unsaturated granular materials, the definition of an unsaturated state involves at least the generalized capillary suction ψ^{x} (Eq. [22]), the degree of saturation *S*^{w} = *V*^{water}/*V*^{pores}, the contact angle θ, and the magnitudes and spatial distribution of the capillary-suction-induced intergranular force vectors **f**^{k}, being normal contact forces similar to **N**^{k} in Eq. [19]. These force vectors **f**^{k} concern a large series of microscale vectors, while all other quantities are scalars defined at the macroscale. The particle surface roughness *s* is a material parameter.

The aim of this section is to replace the effects of the capillary-suction-induced microscale intergranular force vectors **f**^{k} by a representative macroscale capillary-suction-induced intergranular stress tensor **ψ***. Such a capillary-suction-induced intergranular stress tensor has already been derived analytically for homogeneous unsaturated pyramidal packing of rough spheres (Molenkamp and Nazemi, 2003b) for the pendular saturation phase (e.g., Bear, 1972; Nitao and Bear, 1996).

These unsaturated state parameters are mutually dependent, thus distinguishing between functionally independent and dependent quantities is merely a matter of convenience.

*S*

^{w}and its material rate , enabling us to distinguish between wetting and drying through , namely for wetting, thus

^{x}and the macroscale capillary-suction-induced intergranular stress tensor

**ψ*** can be symbolically expressed by

**f**

^{k}with approximately equal magnitude at all water bridges between particles in direct contact. Here, for the low degrees of the pendular saturation phase, the relative humidity plays an equilibrating role (e.g., Molenkamp and Nazemi, 2003a). In such an idealized simple case, the scalar normal magnitudes

*f*

^{k}of the volume-averaged suction-induced interparticle (normal) force vectors

**f**

^{k}can be approximated by the application of the volume-averaged scalar measure of the normal interparticle forces (Eq. [21]), while replacing the magnitudes

*N*

^{k}of the normal intergranular forces by the magnitudes

*f*

^{k}of the approximately normal capillary-suction-induced intergranular force vector

**f**

^{k}, leading to the following scalar volume-averaged capillary-suction-induced contact force measure:

*f*

^{k}, illustrated in Fig. 10, and combining this with the functional form of the contact angle θ according to Eq. [25] leads, for the functional form of the capillary-suction-induced intergranular normal force

*f*

^{k}, to

*s*of the particles surfaces has been assumed to be constant and subsequently for simplicity left out of the following considerations. The corresponding functional form of the volume-averaged suction-induced contact force measure then becomes

**A**, acting here as an independent state parameter of the granular microstructure. In that particular case, the suction-induced intergranular stress tensor

**ψ*** will have the same spatial structure as the fabric stress tensor

**σ**

^{A}according to Eq. [20], leading, for the functional form of the capillary-suction-induced intergranular stress tensor

**ψ***, to

for which Eq. [30] has been substituted in the second expression of Eq. [20], while also accounting for the fabric tensor **A** as an additional tensorial state parameter.

**ψ*** with the functional form according to Eq. [31], first the suction-induced part of Eq. [13] according to Bishop (1959) and Bishop and Blight (1963) is considered:

where the generalized capillary suction ψ^{x} according to Eq. [22] is combined with Bishop’s dimensionless isotropic suction-induced intergranular stress function , while **I** is the second-order isotropic tensor. Note that if the generalized capillary suction ψ^{x} = *p*^{w} – *p*^{x} (Eq. [22]) is positive, and thus the pore water suction *p*^{w} is larger (thus tensile) than the generalized pore air stress *p*^{x}, then the capillary-suction-induced intergranular stress **ψ*** is negative (thus compressive), justifying the minus sign in Eq. [32].

**A**

^{dev}of the fabric tensor

**A**. To this end the fabric tensor

**A**according to Eq. [16] and occurring in Eq. [31] is first decomposed into the isotropic part

**A**

^{iso}=

*A*

^{iso}

**I**and the deviatoric part

**A**

^{dev}by

**ψ*** to

in which the additional dimensional scalar function is introduced to account for the relative magnitude, depending on the degree of saturation, of the deviatoric capillary-suction-induced intergranular stress, defined in more detail below.

**A**

^{dev}= 0, thus ignoring the deviatoric part of the fabric tensor, gives

**ψ***, to

**X**has been introduced:

After recognizing the physical reality of the existence of a deviatoric part of the capillary-suction-induced intergranular stress tensor **ψ***, its practical significance should also be considered, which for unsaturated loose granular materials should at least include its effect on the occurrence of wetting-induced pore collapse.

On this topic we expect as a provisional speculative guess, before considering wetting, that in general terms the real deviatoric part of the combined intergranular stress tensor may be significantly larger than when completely ignoring the deviatoric part of the capillary-suction-induced intergranular stress tensor **ψ***. Then, during subsequent wetting, not only the isotropic part of the combined intergranular stress tensor may decrease but also its deviatoric part, which for loose granular materials may increase the occurrence of volumetric contraction. This induced-deformation mode, together with the ongoing rate of resaturation, may increasing the possibility of the occurrence of a spontaneous instability, recently investigated by Buscarnera (2010) and Buscarnera and Nova (2011), in a somewhat similar way as the spontaneous instability leading to the static liquefaction of saturated loose sand.

## Relations between Capillary-Suction-Induced Quantities

*d*

_{50}/2 represents the mean grain radius. The synthetic relations have been composed while taking account of the earlier microscale characteristics according to Fig. 4, 5, and 7 together with some experimental data on unsaturated granular materials, particularly concerning the dimensionless relations . In addition, for the corresponding dimensionless capillary-suction-induced intergranular stress tensor

**Ψ***, the following isotropic component Ψ

_{p}* and deviatoric component Ψ

_{q}* are distinguished:

in which subscripts *a* and *r* indicate the axial and radial directions, respectively.

Figure 11 illustrates both the dimensionless isotropic suction-induced intergranular stress function (Bishop, 1959; Bishop and Blight, 1963) and the corresponding dimensionless fabric-related deviatoric suction-induced function for drying and wetting. We note that the isotropic Bishop’s functions increase from zero at the dry state with increasing saturation to reach unity at full saturation, while the deviatoric fabric functions start at unity for the dry state and decrease to zero with increasing saturation.

In Fig. 12, the corresponding dimensionless relations of the capillary pore water suction (horizontal to the right side) (Eq. [39]), the capillary-suction-induced isotropic intergranular stress component (vertical downward) (Eq. [40]), and the capillary-suction-induced deviatoric intergranular stress component (horizontal to the left side) (Eq. [40]) are shown. These figures are mainly self-explanatory, although the left-bottom capillary-suction-induced intergranular stress path should be noted to pass through a significant deviatoric stress level during repeated cycles of drying and wetting. From Eq. [37], the latter can be noted to be due to the magnitude of ||**A**^{dev}||/*A*^{iso} for which, based on Fig. 7, ||**A**^{dev}||/*A*^{iso} ≈ 0.33 has been applied.

## Stress Measures for Unsaturated Granular Materials

**σ**equals the sum of the combined intergranular stress tensor

**σ*** and the combined pore fluid stress tensor

**p**, thus (see also Eq. [9] for saturated case)

**σ*** is the sum of the intergranular stress tensor due to the boundary tractions and the capillary-suction-induced intergranular stress tensor

**ψ*** (Eq. [37]), thus

**σ**(Eq. [43]) and the local normal outward unit vector

**n**of the boundary:

*p*

^{x}(Eq. [23]) are the externally exerted intergranular boundary traction and pore gas stress

*p*

^{g}, respectively. This combined traction corresponds to the total stress tensor

**σ**, being the sum of the pore gas stress tensor

*p*

^{x}

**I**(Eq. [23]) and the intergranular stress tensor due to the boundary tractions, thus

while the corresponding total stress tensor **σ** on the boundary follows from substituting Eq. [42] and [54] into Eq. [41], reproducing the total stress tensor **σ** (Eq. [53]) on the boundary.

**ψ*** (Eq. [37]) into the expressions of both the combined intergranular stress tensor

**σ*** in Eq. [42] and the combined pore fluid stress tensor

**p**in Eq. [54] leads, after also substituting Eq. [53] and [22], respectively, to

*p*

^{x}=

*p*

^{g}∩

**X**= 0 (see Fig. 11 and Eq. [23]), it follows from Eq. [55] and [56] after substituting Eq. [53], that

*p*

^{x}=

*p*

^{b}∩

**X**=

**I**(see Fig. 11 and Eq. [24]) and substitution of Eq. [22] and [53] into Eq. [55] and [56], that

We conclude that for both limits of unsaturation, the unsaturated stress measures of the combined intergranular stress tensor **σ*** in Eq. [55] and the combined pore fluid stress tensor **p** in Eq. [56] coincide with the respective expressions for saturated geomechanics.

**σ*** in Eq. [55] and the combined pore fluid stress tensor

**p**in Eq. [56] are accounted for. The corresponding constitutive characteristics are represented by the tangent elasto-visco-plastic flexibility tensor function

**D**

^{−1}of the solid skeleton and the tangent elastic flexibility matrix (

**D**

^{e_sf})

^{−1}of the substance comprising the particles. We note that the tangent isotropic linear elastic flexibility matrix (

**D**

^{e_sf})

^{−1}replaces the compressibility 1/

*K*

^{sf}of the substance comprising the particles, occurring in Eq. [18], while the isotropic pore fluid stress tensor

*p*

**I**in Eq. [18] is replaced by the combined second-order pore fluid stress tensor

**p**(Eq. [56]). The isotropic linear elastic flexibility matrix of the substance comprising the particles is, in index format,

*K*

^{sf}and

*G*

^{sf}are the average bulk modulus and shear modulus of the particles, respectively. Substituting Eq. [61] into the last term of Eq. [60] leads, for the corresponding part of the strain rate due to the elastic flexibility of the particles, to

Note that the first term on the right-hand side of Eq. [62] depends on the bulk stiffness *K*^{sf} of the particles and is caused by the material rates of both the pore air stress (Eq. [23–24]) and the isotropic part of the capillary-suction-induced intergranular stress **ψ***. The last term on the right-hand side of Eq. [62] depends on the shear modulus *G*^{sf} of the particles and is due to the material rate of the deviatoric part of the capillary-suction-induced intergranular stress **ψ***. This term is expected to usually remain negligible compared with the deviatoric strain rate due to the interparticle deformation mode of the granular skeleton as calculated from the first right-side term of Eq. [60].

From Eq. [60] we conclude that for a solid skeleton of granular materials, composed of elastic particles with mainly frictional particle interaction, both developed capillary-based unsaturated stress measures should be sufficient to account for the mechanical effects on the solid skeleton for the full range of unsaturation, thus including both the dry and fully saturated limits. The concept is in general agreement with recently published DEM simulations (e.g., Chareyre et al., 2009; Kim et al., 2012), although the evolution of the spatial distributions of interparticle contact points and water bridges during the loading history may not coincide precisely. For instance, during interparticle slip, the local water bridge will distort at least due to the dependence of the local contact angle θ on the wetting and drying motions of the capillary interface in contact with the solid surface around the interparticle contact point. This will affect the magnitude and direction of the local capillary-suction-induced interparticle force.

## Potential Approach for Progress toward Applications

It is understood that, in principle, the intergranular capillary suction structure tensor function **X** (Eq. [38]) of the described model framework for unsaturated granular materials, like sands and silts, could be determined directly on the basis of microscale observations of the evolution of the intergranular structure and the geometries of the interfaces between the pore water and pore air and the particle surfaces. These microscale observations should be sufficient to calculate both the microscale capillary pore fluid stresses and the corresponding suction-induced interparticle forces. Then subsequent volume averaging of these microscale data should enable the calculation of the macroscale quantities occurring in function **X** (Eq. [38]). Unfortunately, at present the development of such microscale experimental capabilities seems still far away from geotechnical engineering practice. For instance, Fig. 13 shows a (state of the art) two-dimensional picture of microsized spheres of about 10-μm diameter with concave and convex water menisci formed simultaneously at contacts for a pore air relative humidity of 88% (Lourenço et al., 2012). Although this type of data is very important to observe essential features of wetting and drying, for the quantification of the complete microstructure three-dimensional geometrical data are needed for both the granular skeleton and the interfaces of water, air, and grains.

Therefore another more approximate and less demanding and therefore probably much quicker approach for advancing the development of predictive tools for unsaturated sands and silts is considered here, involving the application of inverse modeling, for which the described model framework for unsaturated granular materials is considered to be potentially useful as a basis. The method of determination of the χ parameter of Bishop’s effective stress (Eq. [13]), combining measured shear strengths of saturated and unsaturated soils as demonstrated by Bishop and Blight (1963), can be seen as an early form of such inverse modeling. Inverse modeling by means of the ensemble Kalman filter (e.g., Evensen, 2003, 2006) is expected to allow the simultaneous determination of the macroscale combined capillary characteristics of both pore fluids and the mechanical characteristics of the solid skeleton, provided the existence of an experimental database containing both saturated and unsaturated data on the mechanical behavior of the granular skeleton and both pore fluids. In appropriate inverse modeling, all differences between the experimental observations and the model predictions are minimized. For the required modeling predictions occurring in this minimization process, all stress measures in Eq. [55] and [56] can be used as volume averages. The corresponding deformation modes need to be compatible and thermodynamically justified, for equilibrium minimizing their combined internal energy, in which in principle also the capillary energy as stored in the liquid–gas interface must be accounted for at least approximately.

It should be recognized that the major quantity of this framework for unsaturated conditions concerns the capillary-suction-induced intergranular stress tensor **ψ*** (Eq. [37]). Its macroscale quantities are the degree of saturation *S*^{w}, the generalized capillary matric suction ψ^{x} (Eq. [22]), and the normalized deviatoric fabric tensor **A**^{dev}/*A*^{iso} (Eq. [33]), while the (Bishop) isotropic and deviatoric fitting relations and depend also on the wetting and drying modes. However, without microscale three-dimensional geometrical data on the microstructure, the fabric tensorial quantity **A**^{dev}/*A*^{iso} remains unknown. Nevertheless application of the proposed model framework in inverse modeling can be expected to enable the approximate determination of **A**^{dev}/*A*^{iso} for the following three reasons:

The constitutive relations of the solid skeleton of sands and silts are considered to be independent of the degree of saturation. Consequently, the mechanical characteristics of the granular skeleton can be determined by means of inverse modeling based exclusively on the experimental data on saturated behavior.

The anisotropic effect of the microstructure, described by the normalized deviatoric fabric tensor **A**^{dev}/*A*^{iso}, on the mechanical behavior of the granular skeleton is maximum for the pendular saturation phase and zero for the insular saturation phase. Consequently, this parameter can be derived by means of inverse modeling from unsaturated experimental data on exclusively the pendular saturation phase, while applying the mechanical characteristics of the granular skeleton according to Point 1.

The (Bishop) isotropic and deviatoric fitting curves and can be derived by means of inverse modeling using experimental data on the whole range of unsaturation.

Therefore a database of controlled observations is needed for both fully saturated and unsaturated laboratory tests of the same granular materials and for the full range of initial void ratios occurring in engineering practice. In these experiments, the controlled boundary conditions will concern both pore fluid stresses and total stress, although for the stability of the control (e.g., Buscarnera, 2010; Buscarnera and Nova, 2011) other types of control may also be needed to complete the behavioral overview.

Having arrived at the best-fitting mechanical characteristics of the solid skeleton and the capillary-suction-induced intergranular stress tensor for unsaturated granular materials, the corresponding combined intergranular stress tensor **σ*** (Eq. [55]) and the combined pore fluid stress tensor **p** (Eq. [56]) can finally be calculated.

## Evaluation of Capabilities and Limitations of Proposed Formulation

**σ*** (Eq. [55]). Both are recalled here for convenience:

*p*

^{g}and as the tensor

**σ**−

*p*

^{x}

**I**, while noting

*p*

^{x}=

*p*

^{g}(Eq. [23]) for the pendular and funicular phases and

*p*

^{x}=

*p*

^{b}(Eq. [24]) for the insular phase. The last parts of both Eq. [13] and [55] concern the effects of the matric suction

*p*

^{w}−

*p*

^{g}on the intergranular stress as a normal stress χ(

*p*

^{w}−

*p*

^{g}) and as a tensor (

*p*

^{w}–

*p*

^{x})

**X**. Clearly their main difference concerns the scalar fraction χ and the second-order intergranular capillary suction structure tensor

**X**, which is

Note that in Eq. [38] Bishop’s scalar fraction also occurs, multiplied by the second-order isotropic tensor **I**, thus **I**, and the additional second-order deviatoric tensor **A**^{dev}/*A ^{iso}*. Consequently, the difference between Eq. [13] and [55] concerns only this last deviatoric contribution, representing the anisotropic effect of the intergranular fabric on the capillary-suction-induced intergranular stress tensor

**ψ***.

Therefore it can be concluded that both Eq. [13] and [55] involve the same isotropic capillary-suction-induced intergranular stress, and additionally Eq. [55] also represents a potential deviatoric contribution.

Thus, for evaluating Eq. [55], the significance of this deviatoric contribution should be considered. Current experience, involving Bishop’s effective stress (Eq. [13]), has shown good predictions of shear strength, volume change, and critical state lines for sands and silts (Khalili et al., 2004). At first sight, this would strongly suggest that the significance of the deviatoric part of Eq. [55] remains relatively small and therefore can be neglected.

However, Khalili et al (2004) also expected it to be possible to model pore collapse on wetting by means of a plasticity model by accounting for the hardening effect of matric suction.

On this point, we would like to remark that the very need for matric suction as a constitutive quantity of the solid skeleton brings Bishop’s effective stress formulation outside the framework as considered by us, which is summarized as follows: “The constitutive models of the granular skeleton are considered to depend exclusively on the state parameters of the solid skeleton, thus being independent of, among others, the degree of saturation or the types of pore fluids (or the matric suction).” We expect that the additional deviatoric part of Eq. [55] will enable the modeling of the contraction of the solid skeleton during wetting because, through this deviatoric term, wetting will also induce the reversal of deviatoric loading, increasing the tendency of loose sand to contract.

Consequently, while recognizing the major modeling capabilities of Bishop’s effective stress (Eq. [13]) for sands and silts, we expect that through its deviatoric part, Eq. [55] may enable the modeling of some additional physical features, like pore collapse, without requiring any additional stress measures. In general, for saturated granular materials, the loss of controllability, inducing collapse, occurs for specific combinations of stress and deformation (Nova, 1994; Buscarnera and Nova, 2011).

More complicated granular materials with constitutive models of their solid skeletons, also depending on the degree of saturation, fall outside the scope of the proposed model framework. Nevertheless, for clays with physical-chemical particle interactions, Bishop’s effective stress (Eq. [13]) has been shown to be very useful after accounting for its double-porosity microstructures (e.g., Alonso et al., 2010, 2013). However, the corresponding micro-physics, explaining this microstructure, is still lacking.

For peats, composed of organic, partly decomposed, fibrous and bulk substances (Landva and Pheeney, 1980), illustrated in Fig. 14 (Boylan, 2008), with potential biochemical actions both within and between the solid components and with the pore water and pore air, the proposed model framework must be considered to be not physically justifiable.

In the longer term, appropriate microscale observations may become possible. Expecting for sands and silts that the capillary phenomenon will remain governing, detailed microscale observations of the granular fabric, the water–gas interfaces, and their contacts with the particle surfaces may be sufficient for arriving at an effective and physically fully justified model. However, for clays and peats, the required microscale observations and the physical modeling can be expected to involve additional non-capillary physical characteristics.

## Summary and Conclusions

First, the derivations of the expressions for the microstructure and the intergranular stress for saturated granular materials were reviewed. For the constitutive modeling of solid skeletons composed of elastic (or rigid) particles with mainly frictional particle interaction, the intergranular stress tensor follows from the volume averaging of the microscale interparticle contact forces, while the microstructure is described by a fabric tensor formulated for the packings of spheres. The pores are saturated by air or water, each requiring the definition of a pore fluid stress. The strain rate of the saturated granular material was shown to depend on the rates of both stress measures simultaneously, the intergranular stress and the pore fluid stress.

For unsaturated granular materials, the additional action of capillary water was recognized as inducing capillary suction in the pores and capillary-suction-induced interparticle forces at the microscale. At the macroscale, volume averaging of these forces leads to the capillary-suction-induced intergranular stress tensor **ψ***. In its formulation, the concept of the fabric stress tensor was applied, enabling the effect of the spatial distribution of the capillary water bridges at the particle contacts to be approximately accounted for in the drier pendular saturation phase. To this end, the intergranular fabric tensor was decomposed into its isotropic and deviatoric parts, the latter describing the (anisotropic) deviatoric effect of the matric suction. We noted that in this way Bishop’s χ function occurs as the isotropic part of the capillary-suction-induced intergranular stress tensor.

To illustrate the characteristics of the capillary-suction-induced macroscale measures, some synthetic relations were composed by combining the numerically calculated macroscale results of monotonic deviatoric deformation of the packing of spheres, the analytically derived evolution of the capillary-suction-induced interparticle forces, and some experimental data on unsaturated granular materials, mainly relating the degree of saturation to the capillary matric suction. The applied numerical results of the packing of spheres involve both the evolution of the deviatoric part of the fabric tensor and the fabric stress tensor at peak loading. Finally, the corresponding stress measures of unsaturated granular materials were derived, leading to both the combined intergranular stress tensor **σ*** and the combined pore fluid stress tensor **p**.

The constitutive relations of granular materials, composed of elastic or rigid particles with mainly frictional interaction, were considered to remain independent of the degree of saturation. Under such restrictive conditions, the effects of unsaturation on the mechanical behavior of such granular skeletons should be described by simply replacing the intergranular stress tensor and the isotropic pore fluid stress for the saturated case with the combined intergranular stress tensor and the combined pore fluid stress tensor, respectively.

Consequently, the proposed potential approach of inverse modeling for the determination of unsaturated capillary characteristics, involving mainly the capillary-suction-induced intergranular stress, remains limited to the specified class of granular materials, including quartz sands and silts, namely with elastic (or rigid) particles and mainly frictional particle interaction.

For clays and peats, with potential physical-chemical and bio-chemical actions and double porosity and/or fibrous microstructures, the constitutive models can be expected to be physically more complicated, thus involving more physically relevant parameters. Hence, clays and peats must be considered to fall outside the scope of the proposed model framework.

This research was supported by the Dutch Technology Foundation STW, which is part of the Netherlands Organization of Scientific Research (NWO) and which is partly funded by the Ministry of Economic Affairs, Agriculture and Innovation. We are grateful for the support of Royal Boskalis Westminster N.V. for this research.