Continuous modeling of soil morphology —thermomechanical behavior of embankment dams

Abstract

Macroscopic modeling of soils is based on a number of properties that refer to the mesoscopic morphology. The most fundamental parameters of this art are: 1) coupling parameters between partial stresses of components and deformations of components, 2) porosities, 3) saturation, and 4) permeability and diffusivity, tortuosity.

The main aim of this paper is to present in juxtaposition continuous one-, two-, and three-component models of geomaterials appearing in construction of embankment dams. In particular, the above mentioned features, especially saturation with water and seepage problems, modeling of fluidization yielding piping, and generalizations of the Darcy law and changes of porosity, are presented.

Keywords

thermomechanical modeling / soil morphology / saturation / porosity

Introduction

The inspection of textbooks and manuals for civil and geotechnical engineers reveals that the design of embankment dams and levees is still based on two issues. It is either a stability analysis based on the one-dimensional Mohr-Coulomb relation

(1)\tau =c+\sigma \tan \phi ,

where τ is the shear strength, σ denotes the normal effective stress on the failure plane, and c, ϕ denote the cohesion intercept and the friction angle, respectively, or these are flow nets and streamlines obtained by a graphical, for instance, Schmidt’s, method (e.g., see Fig. 1).

Sometimes it is supplemented by Darcy’s law for the estimation of seepage. The rest of those books contain hundreds of examples of existing constructions, which is a description of their behavior under various loading conditions and failures. Based on this empirical knowledge, some heuristic hints for designers are formulated.

This is very different from engineering books in other branches of civil engineering where the design is based on theoretical models that have replaced a sheer collection of observations. This yields as well a very fruitful development of software for computer aided design, such as CAD, Novapoint, etc. Vaníček et al. [1] wrote in their book: “When we look back on the whole process through which the geotechnical engineer has to go, we arrive at the conclusion that the degree of accuracy is significantly lower than for steel or concrete structure, where the differences in the design can be in the order of a few percent, while for earth structures, these differences can be in the order of tens percent. That is why on one side, the excellent knowledge of soil behavior and treatment of soil as construction material can bring significant savings against conventional design, but on the other one, disregarding this can lead to failures of earth structures.”

However, the situation is slowly changing to the better because the research in the field of soil mechanics has made a tremendous progress, and many theoretical issues, such as failure criteria, fluid flow in porous and granular media, heat transfer in soils, micro-macro transitions in theoretical modeling that incorporate porosity changes, saturation, phase changes, dynamics, and, particularly, thixotropy or sound wave propagation in soils and rocks, were successfully developed.

A choice of theoretical descriptions of aquifers, embankments, and many other geotechnical structures depends on the class of phenomena that we want to embrace and on conditions in which the construction or its part should work. For instance, a mechanical loading of a granular dry material yields fragmentation and abrasion. The same mechanical loading of a water saturated granular material yields diffusion fragmentation but much less abrasion. Hence, in the first case, we may expect considerable changes of porosity, while in the second case, changes of permeability, piping, particle segregation etc., play an important role. This means that the water content in a geomaterial may essentially influence the choice of the theoretical description, which is needed.

In this work, we present a comparison of fundamental approaches to the theoretical description of the thermomechanical behavior of geotechnical materials. First, we sketch a one-component model with additional internal variables. This may be appropriate not only for the description of plastic behavior of geomaterials abrasion but, in many cases of practical interest, also for the description of diffusion. Second, we present a two-component model of a saturated granular material. This model contains a number of additional variables that are able to describe such phenomena as diffusion with variable permeability, localization of deformation (e.g., on filters), and internal erosion processes. This yields a theoretical description of the backward erosion, concentrated leak, and suffusion, which are, in turn, the main reasons for piping. Finally, also, a three-component description for unsaturated granular materials is mentioned.

Thermodynamic modeling of soils developed as a brunch of the theory of immiscible mixtures. This modern continuous approach to systems with multi-component structures whose one component is solid has been initiated by the works of Bowen [2]. Numerous nonlinear effects in porous and granular materials yield in such a modeling a necessity of application of rather complex mathematical tools. This complexity is then hidden in modern computer programs whose applications in engineering do not require high mathematical skills from software users.

However, such multi-component thermodynamic models of soils require as well an identification of various quantities that do not usually appear in classical continuous models. This concerns, in particular, geometrical properties of such systems, true properties of their components, and some process variables characteristic for systems with microstructure. We call this identification the modeling of soil morphology. To name a few examples, one has to identify the porosity, permeability, and, for three-component materials, the saturation, moisture, capillary pressure, and many others. This leads very often to confusion and misinterpretation of results. For this reason, in this work, we show fundamental quantities of thermodynamic models of immiscible mixtures appearing in description of soils and their relation to quantities commonly used by soil engineers.

One-component modeling of geomaterials

The origin of the one-component models of geomaterials stems from the classical model of elastoplastic materials. They belong to two groups: one describing dry granular materials driven by elastic properties and frictional interactions of grains, and the other one describing fully saturated granular materials in which viscosity rather than friction contributes to the mechanical response of the system. Both classes of models contain the macroscopic deformation \mathit{B}=B_{ij}{\mathit{e}}_i\otimes {\mathit{e}}_j (the left Cauchy-Green deformation tensor), the velocity BoldItalic = v_iBoldItalic_i, and the temperature θ as unknown functions of the position BoldItalic and time t.

In linear models, one uses often the displacement vector, u= u_kBoldItalic_k, related to the Almansi-Hamel deformation tensor, BoldItalic, and the velocity, BoldItalic, by the following relations:

(2)\mathit{e}=e_{kl}{\mathit{e}}_k\otimes {\mathit{e}}_l,e_{kl}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\displaystyle \frac{\partial u_k}{\partial x_l}}\right)+\left({\displaystyle \frac{\partial u_l}{\partial x_k}}\right),\mathit{v}=v_k{\mathit{e}}_k,v_k={\displaystyle \frac{\partial u_k}{\partial t}}.

This deformation measure is related to the left Cauchy-Green tensor by the equation:

(3)\mathit{e}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}(\mathrm{1}-{\mathit{B}}^{-\mathrm{1}}).

However, these models differ in the set of unknown microstructural variables. The first class contains only the roughness a whose time derivative \dot{a} is called the abrasion, while the second class may contain the abrasion, but it must contain also the porosity n and the pore pressure p. For these quantities-fields, we have to construct additional equations.

The classical approach is based on the set of conservation laws of mass, momentum, and energy. In Cartesian reference frame, they have the following form:

(4){\displaystyle \frac{\partial \rho }{\partial t}}+{\displaystyle \frac{\partial (\rho v_i)}{\partial x_i}}=\mathrm{0},

(5){\displaystyle \frac{\partial (\rho v_i)}{\partial t}}+{\displaystyle \frac{\partial }{\partial x_k}}(\rho v_k{v}_i-{\sigma }_{ik})=\rho b_i,

(6){\displaystyle \frac{\partial }{\partial t}}\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\rho v_k{v}_k+\rho ϵ\right)+{\displaystyle \frac{\partial }{\partial x_k}}\left[\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\rho v_i{v}_i+\rho ϵ\right)v_k+q_k-{\sigma }_{ki}{v}_i\right]=\rho b_k{v}_k,

where ρ is the bulk mass density, σ_ik are components of the Cauchy stress tensor, \mathit{T}={\sigma }_{ik}{\mathit{e}}_i\otimes {\mathit{e}}_k, b_i are the body forces (e.g., gravitational or centrifugal), BoldItalic = b_kBoldItalic_k, ϵ is the specific internal energy, and q_k are components of the heat flux vector, BoldItalic = q_kBoldItalic_k. It is often assumed that the real grains of the material are incompressible. This means that the true mass density ρ^SR, ρ = (1-n)ρ^SR (S for “solid” and R for “real” or “true”; in soil mechanics, one denotes sometimes ρ^SR = γ) is constant. The porosity n is related to the void ratio e, 0≤e<∞, often used in soil mechanics, by the simple relation n = e/(1+ e), 0≤n≤1.

All these arguments are related to the two-component model that we discuss in Section 3. Therefore, the changes of porosity following from the incompressibility assumption cannot be consistently incorporated in a one-component model. The following relation is used in this model:

(7)n=n_{\mathrm{0}}\left(\mathrm{1}+{\displaystyle \frac{\mathrm{1}-n_{\mathrm{0}}}{n_{\mathrm{0}}}}e\right),

where n is the current porosity, n₀ the initial porosity, and e = tr(1-BoldItalic^-1) is the volume changes of the solid frame (caution: this e is not the same as the void ratio mentioned above) that must be considered as an additional assumption of the one-component model.

We return to the microscopic interpretation of these quantities in Section 4.1 of the present work.

To obtain the equations of the model of dry granular materials, we have to specify constitutive relations for stress tensor σ_ik, internal energy ϵ, heat flux q_k, and the abrasion \dot{a}.

Experience shows that granular materials behave plastically. Consequently, the classical Mohr-Coulomb relation has been extended to relate the stress tensor and the deformation tensor. As the so-called hardening effects play an important role in such models, one had to introduce additional internal variables (the so-called back-stress). The result is the cam-clay model commonly used in the literature on soil mechanics (e.g., Refs. [3-5]). As an alternative, a so-called hypoplasticity was introduced by Bauer [6], Wolffersdorff [7], and Kolymbas [8,9]). In contrast to the cam-clay model, the hypoplasticity is rate dependent, which means that the rate of deformation D_{ij}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left[{\displaystyle \frac{\partial v_i}{\partial x_j}}+{\displaystyle \frac{\partial v_j}{\partial x_i}}\right] (it is related to the time derivative of the deformation tensor B_ij) has an influence on the current values of the stress. The general constitutive relation for stresses in this model has the form

(8){\stackrel{\circ }{\sigma }}_{ij}=f_{ij}[{\sigma }_{ij},D_{ij},n],

where {\stackrel{\circ }{\sigma }}_{ij} denotes an objective time derivative of the stress tensor.

We do not need to go into details of these models in the present work. Many of them can be found in Ref. [10].

It remains to specify the internal energy ϵ, the heat flux q_k, and the abrasion \dot{a}. For the first two quantities, one usually assumes that the classical Fourier model of heat conduction is valid. This may be questionable in some fast processes, but for the thermomechanical description of embankments, it seems to be sufficient. For the purpose of our analysis, it is sufficient to point out the ways in which energy is transported in the medium. They are specified in Eq. (6) by contributions under the div operator (i.e.,{\displaystyle \frac{\partial }{\partial x_k}}). The first one is the convection. The second one, described by BoldItalic, consists of two parts: the conduction, BoldItalic_c and the diffusion, BoldItalic_d with BoldItalic = BoldItalic_c + BoldItalic_d. The latter means that the energy is transported by the relative motion of components. Such a mechanism cannot be described by a one-component model. The conduction is related to the transfer of energy due to the temperature gradient. Finally, the last contribution, BoldItalic = σ_iv_iBoldItalic_k, is the bulk working of stresses that is also of no interest in this work.

Conduction in isotropic materials is usually described by the linear Fourier law

(9){\mathit{q}}_{\mathrm{c}}=-\lambda \mathrm{g}\mathrm{r}\mathrm{a}\mathrm{d}\theta ,

where θ is the absolute temperature. The coefficient λ, the heat conductivity, is for soils heavily dependent on the morphology. This dependence cannot be reflected by a one-component model either. We return to this problem in Section 4.3.

The form of the equation for abrasion has a long history, and it goes back to the work of Goodman et al. [11] However, in this pioneering work, the equation was proposed rather for changes of volume fraction than for the abrasion. It was first the series of works of Wang et al. (e.g., Refs. [12,13]) and the PhD Thesis of Kirchner [14], where this equation was thermodynamically justified. Its form follows from the assumption that the microstructural changes of the configuration caused by the abrasion must be accompanied by the so-called configurational forces. Then, the abrasion \dot{a}={\displaystyle \frac{\partial a}{\partial t}}, corresponding to the classical notion of momentum, must satisfy a balance law, which is assumed to have the form

(10)\rho k{\displaystyle \frac{{\partial }^{\mathrm{2}}a}{\partial t^{\mathrm{2}}}}={\displaystyle \frac{\partial h_i}{\partial x_i}}+\rho (l+f),

where k is the material parameter describing the resistance of the material to changes of its internal surface. According to Kirchner [14], “change of surface properties includes the smoothening of initially rough grain surface (that is, ) as well as the roughening of initially smooth grain surfaces (that is, \dot{a}>\mathrm{0}).” h_i is the surface stress of abrasion. l and f are supply and production, respectively, and the latter must be given by a constitutive law of its own.

The above sketched one-component model of dry geomaterials is often extended by an equation describing the flow of water through the saturated granular material. All equations described above are assumed to remain unchanged. The seepage process through the saturated material is supposed to satisfy some additional balance law, which is justified experimentally. Such a justification goes back on works of Darcy [15], and it has been incorporated in soil mechanics by von Terzaghi [16]. In the local form, this law can be written in the form:

(11)Q_i=-{\displaystyle \frac{k_{ij}}{\mu }}{\displaystyle \frac{\partial p}{\partial x_j}},

where Q_i is the so-called specific discharge, i.e., a relative velocity of water with respect to the skeleton, k_ij is the matrix of permeability that reduces to a scalar k for isotropic materials, and μ is the kinematic viscosity of water. It has been shown that such a relation holds for small relative velocities that are connected with a laminar flow of the water (low Reynolds number, Re<1/10 [17]). In the case of fast flows yielding turbulence (high Reynolds numbers), the Darcy law does not hold. Most likely, Forchheimer [18] was the first who proposed nonlinear corrections to Eq. (11) in order to describe such flows. They play a particularly important role in processes of piping, commonly appearing in embankment dams. A more rational procedure of description of seepage is proposed by theories of multicomponent systems.

Two-component modeling of geomaterials

The thermomechanical model of a one-component geomaterial can be considerably improved when one applies a theory of immiscible mixtures. We shall do so first for fully saturated materials. In this case, the extension yields a better physical insight, but for purposes of geotechnics, it is not necessary. However, for partially saturated materials, such extensions are unavoidable, and, simultaneously, they are similar to two-component models in many technical details. We will see this in the following sections on multi-component modeling.

In the case of two components, one has to describe the partial macroscopic fields for each component. These are partial mass densities ρ^S, ρ^F with the bulk mass density ρ = ρ^S + ρ^F, partial velocities BoldItalic^S, BoldItalic^F with the bulk (barycentric) velocity BoldItalic = (ρ^S/ρ)BoldItalic^S + (ρ^F/ρ)BoldItalic^F, partial Cauchy stresses {\mathit{T}}^{\mathrm{S}}={\sigma }_{ij}^{\mathrm{S}}{\mathit{e}}_i\otimes {\mathit{e}}_j, {\mathit{T}}^{\mathrm{F}}={\sigma }_{ij}^{\mathrm{F}}{\mathit{e}}_i\otimes {\mathit{e}}_j with the bulk stress BoldItalic≈BoldItalic^S + BoldItalic^F. These quantities must satisfy balance laws

(12){\displaystyle \frac{\partial {\rho }^{\alpha }}{\partial t}}+{\displaystyle \frac{\partial }{\partial x}}({\rho }^{\alpha }{v}_i^{\alpha })=\mathrm{0},\alpha =\mathrm{S},\mathrm{F},

(13){\displaystyle \frac{\partial ({\rho }^{\alpha }{v}_i^{\alpha })}{\partial t}}+{\displaystyle \frac{\partial }{\partial x_k}}({\rho }^{\alpha }{v}_i^{\alpha }\otimes v_k^{\alpha })={\displaystyle \frac{\partial {\sigma }_{ik}^{\alpha }}{\partial x_k}}+{\widehat{p}}_i^{\alpha }+{\rho }^{\alpha }{b}_i^{\alpha },{\widehat{p}}_i^{\mathrm{S}}+{\widehat{p}}_i^{\mathrm{F}}=\mathrm{0.}

It is easy to check that the conservation laws of a one-component model are then identically satisfied provided; we neglect quadratic terms in relative velocities BoldItalic^a - BoldItalic that seems to be well justified for processes in soils far from the structural loss of stability, such as fluidization. The momentum source {\widehat{\mathit{p}}}^{\mathrm{S}}={\widehat{p}}_i^{\mathrm{S}}{\mathit{e}}_i=-{\widehat{\mathit{p}}}^{\mathrm{F}}=-{\widehat{p}}_i^{\mathrm{F}}{\mathit{e}}_i is related to the diffusion force. In the case of an isotropic model linear in relative velocities, it can be written in the form:

(14){\widehat{\mathit{p}}}^{\mathrm{S}}=\pi ({\mathit{v}}^{\mathrm{F}}-{\mathit{v}}^{\mathrm{S}}),

(where π is the permeability coefficient. In the case of water, one can assume that the partial stress tensor BoldItalic^F is spherical, i.e., it reduces to the partial pressure p^F. Then, the momentum balance for the fluid written in Cartesian coordinates has the form:

(15){\rho }^{\mathrm{F}}\left({\displaystyle \frac{\partial v_i^{\mathrm{F}}}{\partial t}}+v_k^{\mathrm{F}}{\displaystyle \frac{\partial v_i^{\mathrm{F}}}{\partial x_k}}\right)=-{\displaystyle \frac{\partial p^{\mathrm{F}}}{\partial x_i}}+\pi (v_i^{\mathrm{F}}-v_i^{\mathrm{S}})+{\rho }^{\mathrm{F}}{b}_i^{\mathrm{F}}.

This equation yields Darcy’s law for processes with small changes of porosity and negligible inertial forces. In such case, it follows that

(16)p^{\mathrm{F}}=n_{\mathrm{0}}p,(v_i^{\mathrm{F}}-v_i^{\mathrm{S}})=-{\displaystyle \frac{n_{\mathrm{0}}}{\pi }}{\displaystyle \frac{\partial p}{\partial x_i}},

where n₀ is the initial porosity and n≈n₀. This is identical with Eq. (11) for isotropic materials with an appropriate definition of the permeability k/µ = n₀/π. Hence, the two-component model yields the one-component model as a particular case.

On the other hand, the general form of partial momentum balance (12) admits also nonlinear contributions that yield the loss of stability of the fluid motion. If the changes of porosity and relative velocity are not small, one can introduce the following relation for the source of momentum [19]:

(17){\widehat{p}}_i^{\mathrm{S}}=\pi (v_i^{\mathrm{F}}-v_i^{\mathrm{S}})-\left(p+{\rho }^{\mathrm{S}}{\displaystyle \frac{\partial {\psi }^{\mathrm{S}}}{\partial n}}\right){\displaystyle \frac{\partial n}{\partial x_i}},

where p is the pore pressure, and ψ^S is the Helmholtz free energy of the solid component dependent on deformations, porosity, and relative velocity. This form of the source is justified by thermodynamic considerations that we shall not discuss in this work. The simplest choice of the dependence on the relative velocity that yields piping is given as follows:

(18){\rho }^{\mathrm{S}}{\displaystyle \frac{\partial {\psi }^{\mathrm{S}}}{\partial n}}={\displaystyle \frac{\mathit{\Gamma }}{\sqrt{\mathrm{2}}}}\left(\mathrm{1}+{\displaystyle \frac{W-Y}{|W-Y|}}\right)\sqrt{W},\mathit{\Gamma },Y>\mathrm{0},W={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}(v_i^{\mathrm{F}}-v_i^{\mathrm{S}})(v_i^{\mathrm{F}}-v_i^{\mathrm{S}}),

where Γ is a material parameter, and Y is the threshold velocity. As shown in Ref. [19], this model yields a quantitative agreement with experiments on sands.

Changes of porosity

Changes of porosity may be described by relation (7), following from the assumption on incompressibility of grains. However, there is an evidence stemming from poroacoustics that such an assumption eliminates an important P2-wave from the model [20]. Many other problems must also account for changes of porosity. For instance, within soil mechanics, these are large plastic deformations [6], damage (for instance, in freezing and thawing), or, for granular materials, combustion problems of solid fuels.

Generally, there are a few different approaches to changes of porosity. We mention here the five most commonly appearing in the literature:

1) constitutive assumption, e.g., equilibrium changes coupled on volume changes of skeleton (Gassmann [21]);

2) incompressibility assumption (Bowen [2]);

3) evolution equation (Bowen [2]);

4) second-order equation based on a principle of equilibrated pressure (Goodman and Cowin [11]), extended by Passman et al. [23] and Faria et al. [24];

5) balance equation(Wilmanski [25,26]).

Gassmann’s model follows from a simplified micro-macro description and results in the relation

(19)n=n_{\mathrm{0}}(\mathrm{1}+\delta e),

where e denotes small volume changes of the skeleton, and δ is a material parameter related to compressibilities of components. n₀ is the initial porosity. We show some properties of this model further in this section.

The assumption of incompressibility that is essential for Bowen’s approach [2] states that

(20){\rho }^{\mathrm{S}\mathrm{R}}=\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{s}\mathrm{t}.,

which yields the following form of the macroscopic mass balance equation for the skeleton

(21){\rho }^{\mathrm{S}\mathrm{R}}{\displaystyle \frac{\partial }{\partial t}}(\mathrm{1}-n)+{\rho }^{\mathrm{S}\mathrm{R}}\mathrm{d}\mathrm{i}\mathrm{v}[(\mathrm{1}-n){\mathit{v}}^{\mathrm{S}}]=\mathrm{0},

provided; there is no mass exchange with other components.

Easy integration of this equation yields for small deformations

(22)n=n_{\mathrm{0}}\left(\mathrm{1}+{\displaystyle \frac{\mathrm{1}-n_{\mathrm{0}}}{n_{\mathrm{0}}}}e\right),

which is, obviously, Eq. (7) of the one-component model, and it reminds Gassmann’s relation (19), but there is no relation to material parameters. In this case, δ = (1-n₀)/ n₀. We will not discuss the next two models and mention only that the evolution equation proposed by Bowen is a particular case of the balance equation of porosity, at least for small deformations, when one neglects the influence of diffusion. The Goodman and Cowin proposition is related to some microstructural considerations that have a bearing in the case of combustion problems for powders. Hutter et al. have shown that some extensions of this model describe well the behavior of snow avalanches and mud flows.

The porosity balance of the last model, mentioned above, is thermodynamically admissible and yields a consistent model for large deformations [27]. As we will see below, in the case of small deformations of soils, this equation can be immediately solved. Its linear version for two-component systems has the form

(23){\displaystyle \frac{\partial {\mathit{\Delta }}_n}{\partial t}}+\mathit{\Phi }\mathrm{d}\mathrm{i}\mathrm{v}({\mathit{v}}^{\mathrm{F}}-{\mathit{v}}^{\mathrm{S}})=-{\displaystyle \frac{{\mathit{\Delta }}_n}{\tau }},{\mathit{\Delta }}_n=n-n_{\mathrm{E}},n_{\mathrm{E}}=n_{\mathrm{0}}(\mathrm{1}+\delta e),

where BoldItalic^F is the macroscopic velocity of the fluid, and τ, δ, and Φ are material constants. In the case of soils, one can usually neglect relaxation effects described by the right-hand side of the equation. This means that one can take the limit τ→∞. Then, the equation can be integrated, and the following relation follows:

(24)n=n_{\mathrm{0}}(\mathrm{1}+\delta e+\gamma (e-ϵ)),\gamma ={\displaystyle \frac{\mathit{\Phi }}{n_{\mathrm{0}}}},

where ϵ is the volume change of the fluid component. The quantity ζ = (e - ϵ) /n₀ is called the increment of fluid content. The above relation can be easily extended to three components [28]. We shall not present this equation in this work. Similar to the Gassmann relation, it can be shown that material parameters δ, γ can be identified by means of compressibilities of components [29]. Simultaneously, they coincide with coefficients of porosity changes predicted by Biot’s model [30]. In Fig. 2, we present a few examples of the behavior of these parameters. The left panel shows the dependence of the coefficient δ (equilibrium changes of porosity) as a function of the initial porosity, n₀. The curves correspond to the compressibility of the solid skeleton 48 GPa and 35 GPa, respectively. The two lower curves are plotted for water, and the two coinciding upper curves are plotted for air. The middle panel shows the illustration of δ for an incompressible skeleton (relation (22); in this case, γ ≡ 0). There is almost no difference between this curve and the curves for air in the left panel. Finally, the right panel shows the behavior of γ, which reflects the influence of diffusion on changes of porosity. Obviously, this influence is rather small, and the values for the saturation with air are so small that the curves are not visible in the scale in Fig. 2.

Tortuosity

The notion of tortuosity contributes to numerous confusions in modeling porous media. It is still disputable in what way one should include this measure of complexity of microstructure. The tortuosity, τ, is taken as the ratio of the average pore length to the macroscopic characteristic length along, for instance, the major flow (i.e., the longest local streamline).

An example of different values of tortuosity for a medium with the same local porosity is exemplified in Fig. 3. The left situation yields the smallest value of τ and the right situation–the largest value of τ among these three cases. The value of τ is the ratio of the distance of points A and B to the length of channels indicated by the arrows.

As shown by Epstein [31], the tortuosity influences the intensity of diffusion by entering a material parameter relating the pressure gradient and the diffusion velocity. We return to this coefficient in the next section. However, it should be mentioned that it is the square of τ that appears in this relation. The error of the linear dependence made by Kozeny [32] has been corrected in many works, and the quadratic dependence seems to be well established.

In many papers on acoustics of porous media, it is claimed that the tortuosity enters the model by the so-called added mass effect. This has been introduced by Biot [33] in the form of an off-diagonal contribution ρ₁₂ to the partial mass matrix. It can be understood as a coupling of components through inertial forces. One can show in a simple thermodynamic analysis [34] that such a coupling is nondissipative. This means that the tortuosity cannot have an influence on the damping of acoustic waves. This is, of course, a nonsensical conclusion. This contribution can be indeed introduced to poroelastic models after some nonlinear corrections, but it cannot be interpreted as an influence of tortuosity.

Permeability

The notion of permeability of porous media is usually related with the Darcy law, which expresses a total discharge of fluid in a one-dimensional flow in terms of the pressure difference. However, there is still a bit of confusion in the terminology. We use in the model of porous materials the term of a coefficient of permeability or coefficients of permeability if we deal with more than one fluid component. We proceed to introduce these notions in a systematic way.

The fundamental notion of intrinsic permeability, κ, is related solely to the morphology of a porous skeleton, and it is independent of the kind of fluid in the channels. It may be related to the effective diameter of the pores, d, by the relation

(25)\kappa =C{d}^{\mathrm{2}},

where C is a dimensionless constant. The units of κ are m² but in practical applications often the unit 1 [darcy] ≡ 1 D≈10^-12 m² is used.

In the second column in Table 1, we quote some typical values of this parameter.

To characterize the flow of a particular fluid, one has to account also for properties of this fluid. If the kinematic viscosity is µ and its mass density is ρ^FR, then the parameter

(26)K={\displaystyle \frac{\kappa {\gamma }^{\mathrm{F}\mathrm{R}}}{\mu }},{\gamma }^{\mathrm{F}\mathrm{R}}={\rho }^{\mathrm{F}\mathrm{R}}g

is called the hydraulic conductivity; g is the earth acceleration. A few typical values for water (pressure 10⁵ Pa, temperature 20°C, µ = 1.002 × 10^-3 Pa·s) are quoted also in Table 1. Then, the typical form of the Darcy law is given as follows:

(27)Q({\rho }^{\mathrm{F}\mathrm{R}}g)=KA{\displaystyle \frac{\mathrm{\Delta }p}{L}},

where Q (m³/s) is the total discharge through the surface A, ∆p= p_a - p_b is the pressure difference, and L is the distance between two faces a and b.

In addition, for horizontal aquifers, one uses the notion of transmissivity, T, which is the product of the hydraulic conductivity, K, and the thickness, d, of the aquifer.

If we refer the diffusive flow to the difference of the concentration rather than to the difference of the pressure—one speaks then about the first Fick law—then the diffusion flux and the gradient of concentration are connected by the diffusion coefficient, D. This notion is, of course related to the hydraulic conductivity, K, by a simple change of variables.

Finally, let us remark that the constant C of the relation (25) is dependent on the porosity and on the tortuosity (a quadratic dependence as we argued before).

Theories of porous media are based on the model of immiscible mixtures. Then, Darcy’s law does not enter the model at all. It is replaced by partial momentum equations. In the case of a linear model of two components, the simplest form of this equation is shown as follows:

(28){\rho }^{\mathrm{F}}{\displaystyle \frac{\partial {\mathit{v}}^{\mathrm{F}}}{\partial t}}=-\mathrm{g}\mathrm{r}\mathrm{a}\mathrm{d}{p}^{\mathrm{F}}+\pi ({\mathit{v}}^{\mathrm{F}}-{\mathit{v}}^{\mathrm{S}}).

If we neglect the acceleration, it becomes a precursor of Darcy’s law. The coefficient π, which is inversely proportional to the hydraulic conductivity, is called the coefficient of permeability.

The microstructural justification of the above relations for parameters of permeability is difficult. Such laws are known from the kinetic theory of mixture of gases, and in the case of granular materials, they are mimicked by assuming that grains have a very big mass in comparison to the gas in pores and, consequently, can be assumed to be immobile. One obtains the so-called gas-dust diffusivity [35], which is the counterpart of the hydraulic conductivity. Some results, also experimental, are obtained for regular geometries of channels. Otherwise, one has to rely on purely macroscopic observations.

We leave out considerations concerning heat conduction. This problem shall be presented in some details in Section 4 on three-component models.

One can conclude that the above remarks that the two-component model may play an important role in geotechnics for processes in which nonlinearities are essential. This concerns large deformations and, consequently, large changes of porosity and permeability. Of particular importance are, however, large relative (seepage) velocities that yield the loss of stability and piping. Otherwise, one-component models seem to be acceptable for both dry and wet granular materials of geotechnical bearing.

Three-component modeling of geomaterials

As indicated in the previous section, for three-component continuous models besides the porosity, at least one more microstructural variable appears, namely, the saturation, i.e., the fraction of the volumes of one of the pore fluids to this of the void space. Then, in the thermodynamical analysis, the following fields have to be considered:

1) partial mass densities, ρ^S, ρ^F, ρ^G;

2) velocities of components, BoldItalic^S, BoldItalic^F, BoldItalic^G;

3) common temperature of components, T;

4) porosity, n;

5) saturation, S.

The physical significance of the first and second group will be explained in the following section.

Micro-macro transitions for porous materials

The construction of thermodynamic models of materials with small randomly distributed voids requires always a smearing-out procedure, which transforms functions on complicated and different domains to a common domain, say, B₀, which, in turn, is identified with a continuous porous medium. Such constructions for multi-component systems with a relative motion (diffusion) yield the so-called immiscible mixtures. Quantities appearing in such models are not directly measurable and require, as we indicated in the introduction to this work, certain identification rules. In some cases, these rules can be simply constructed by volume averaging, and in some other cases, they require more sophisticated methods of identification. We begin our presentation with a few typical quantities which follow from the former procedure.

The volume averaging is performed by means of the so-called representative elementary volume (REV, e.g., Ref. [17]) consisting of points BoldItalic. REV is small enough to be replaced by a material point BoldItalic of B₀ at which a set of macroscopic quantities replaces real (true) quantities of the microdomain REV. For three-component systems that we consider in this section, the following quantities are defined by averaging over REV(BoldItalic, t) prescribed to the point BoldItalic∈B₀ at the instant of time t:

1) Porosity

(29)\begin{array}{l}n(\mathit{X},t)=\mathrm{1}-{\displaystyle \frac{\mathrm{1}}{\mathrm{v}\mathrm{o}\mathrm{l}(\mathrm{R}\mathrm{E}\mathrm{V})}}{\displaystyle \underset{\mathrm{R}\mathrm{E}\mathrm{V}(\mathit{X},t)}{\int }{\chi }^S(\mathit{Z},t)\mathrm{d}{V}_{\mathit{Z}}},\hfill \\ V_{\mathrm{S}}={\displaystyle \underset{\mathrm{R}\mathrm{E}\mathrm{V}(\mathit{X},t)}{\int }{\chi }^{\mathrm{S}}(\mathit{Z},t)\mathrm{d}{V}_{\mathit{Z}},}\hfill \end{array}

where χ^S(·, t) is the characteristic function of the skeleton, i.e., it is equal to one if the point BoldItalic is occupied by a particle of the skeleton (solid phase) at the instant of time t, and zero, otherwise, and vol(REV) = V is the volume of REV. In continuous models, the porosity is usually identical with the porosity available for the transport in pores. This means that, for instance, contributions of dead-end channels are not included in n. The latter contribute to changes of the true (real) effective properties of the skeleton.

2) Mass density of the skeleton

(30){\rho }^{\mathrm{S}}(\mathit{X},t)={\displaystyle \frac{\mathrm{1}}{\mathrm{v}\mathrm{o}\mathrm{l}(\mathrm{R}\mathrm{E}\mathrm{V})}}{\displaystyle \underset{\mathrm{R}\mathrm{E}\mathrm{V}(\mathit{X},t)}{\int }{\rho }^{\mathrm{S}\mathrm{R}}(\mathit{Z},t){\chi }^{\mathrm{S}}(\mathit{Z},t)\mathrm{d}{V}_{\mathit{Z}}},

where ρ^SR(BoldItalic, t) is the real (true) mass density of the skeleton at the point BoldItalic and at the instant of time t.

3) Mass density of the fluid

(31)\begin{array}{l}{\rho }^{\mathrm{F}}(\mathit{X},t)={\displaystyle \frac{\mathrm{1}}{\mathrm{v}\mathrm{o}\mathrm{l}(\mathrm{R}\mathrm{E}\mathrm{V})}}{\displaystyle \underset{\mathrm{R}\mathrm{E}\mathrm{V}(\mathit{X},t)}{\int }{\rho }^{\mathrm{F}\mathrm{R}}(\mathit{Z},t){\chi }^{\mathrm{F}}(\mathit{Z},t)\mathrm{d}{V}_{\mathit{Z}}},\hfill \\ V_{\mathrm{F}}={\displaystyle \underset{\mathrm{R}\mathrm{E}\mathrm{V}(\mathit{X},t)}{\int }{\chi }^{\mathrm{F}}(\mathit{Z},t)\mathrm{d}{V}_{\mathit{Z}}},\hfill \end{array}

where ρ^FR(BoldItalic, t) is the real (true) mass density of the fluid at the point BoldItalic and at the instant of time t, and χ^F(·, t) is the characteristic function of the fluid, i.e., it is equal to one if the point BoldItalic is occupied by a particle of the fluid at the instant of time t and zero, otherwise.

4) Mass density of the gas

(32)\begin{array}{l}{\rho }^{\mathrm{G}}(\mathit{X},t)={\displaystyle \frac{\mathrm{1}}{\mathrm{v}\mathrm{o}\mathrm{l}(\mathrm{R}\mathrm{E}\mathrm{V})}}{\displaystyle \underset{\mathrm{R}\mathrm{E}\mathrm{V}(\mathit{X},t)}{\int }{\rho }^{\mathrm{G}\mathrm{R}}(\mathit{Z},t){\chi }^{\mathrm{G}}(\mathit{Z},t)\mathrm{d}{V}_{\mathit{Z}}},\hfill \\ V_{\mathrm{G}}={\displaystyle \underset{\mathrm{R}\mathrm{E}\mathrm{V}(\mathit{X},t)}{\int }{\chi }^{\mathrm{G}}(\mathit{Z},t)\mathrm{d}{V}_{\mathit{Z}}},\hfill \end{array}

where ρ^GR(BoldItalic, t) is the real (true) mass density of the gas at the point BoldItalic and at the instant of time t, and χ^G(·, t) is the characteristic function of the gas, i.e., it is equal to one if the point BoldItalic is occupied by a particle of the gas at the instant of time t and zero, otherwise.

In many cases of practical importance, one assumes a microhomogeneity of the microstructure, which means that true mass densities are approximately constant on REV. Then

(33){\rho }^{\mathrm{S}}=(\mathrm{1}-n){\rho }^{\mathrm{S}\mathrm{R}},{\rho }^{\mathrm{F}}=nS{\rho }^{\mathrm{F}\mathrm{R}},{\rho }^{\mathrm{G}}=n(\mathrm{1}-S){\rho }^{\mathrm{G}\mathrm{R}},n={\displaystyle \frac{V_{\mathrm{F}}+V_{\mathrm{G}}}{V}},

where S denotes the saturation

(34)S={\displaystyle \frac{V_{\mathrm{F}}}{V_{\mathrm{F}}+V_{\mathrm{G}}}},

and ρ^SR, ρ^FR, and ρ^GR are evaluated in an arbitrary point BoldItalic₀∈REV (BoldItalic, t) (microhomogeneity!).

Sometimes the product

(35)\theta =nS\equiv {\displaystyle \frac{V_{\mathrm{F}}}{V}},V=V_{\mathrm{S}}+V_{\mathrm{F}}+V_{\mathrm{G}},

is called the volumetric water content. This should be distinguished from the gravimetric water content (a moisture fraction):

(36)\begin{array}{l}w\equiv m_{\mathrm{W}}={\displaystyle \frac{M_{\mathrm{F}}}{M_{\mathrm{S}}}}={\displaystyle \frac{{\rho }^{\mathrm{F}\mathrm{R}}{V}_{\mathrm{F}}}{{\rho }^{\mathrm{S}\mathrm{R}}{V}_{\mathrm{S}}}},\hfill \\ V_{\mathrm{L}}=nSV,V_{\mathrm{S}}=(\mathrm{1}-n)V,\hfill \end{array}

i.e., it is the fraction of the mass of fluid to the mass of the (real, true, dry) skeleton. Certainly, bearing relations (33) in mind, it can be written in the form

(37)m_{\mathrm{W}}={\displaystyle \frac{{\rho }^{\mathrm{F}}}{{\rho }^{\mathrm{S}}}}.

Consequently, in contrast to porosity, n, and saturation, S, neither the volumetric water content, θ, nor the moisture fraction, m_w, are independent microstructural quantities. However, the latter is sometimes used as a measure of compactness of soils reflected by the partial mass density of the skeleton, ρ^S.

The construction of volume averages illustrated above fails in the case of transport coefficients. This concerns both classical coefficients, such as heat conductivity or hydraulic conductivity, as well as partial stresses. In the latter case, one makes sometimes simplifying assumptions

(38)p^{\mathrm{S}}=(\mathrm{1}-n)p^{\mathrm{S}\mathrm{R}},p^{\mathrm{F}}=nS{p}^{\mathrm{F}\mathrm{R}},p=n(\mathrm{1}-S)p^{\mathrm{G}\mathrm{R}},p^{\mathrm{L}}=p^{\mathrm{F}}+p^{\mathrm{G}}=n{p}^{\mathrm{L}\mathrm{R}},p^{\mathrm{L}\mathrm{R}}=S{p}^{\mathrm{F}\mathrm{R}}+(\mathrm{1}-S)p^{\mathrm{G}\mathrm{R}},

where p^S = -1/3trBoldItalic^S is the partial pressure in the skeleton (BoldItalic^S denotes the partial Cauchy stress tensor in the skeleton), p^SR denotes the true (real) pressure in the skeleton, p^F and p^FR are partial and, respectively, real pressure in the fluid, p^G and p^GR in the gas, and p^LR is the pore pressure. One can show that relations (38) are a combination of volume averaging under very restrictive conditions and Dalton’s law for fluid components. However, one should bear in mind that they may not hold, for instance, near boundaries or for dynamical processes, such as shock waves. The following quantity

(39)p_{\mathrm{c}}=p^{\mathrm{G}\mathrm{R}}-p^{\mathrm{F}\mathrm{R}}

is identified with the capillary pressure and used in constitutive relations for the saturation, S.

Fields of velocities of components are more difficult to interpret in terms of microstructural quantities. It is obvious that deviations of real velocities of the skeleton and of the fluid components from their macroscopic counterparts, both with respect to directions and the magnitude, are usually very large. Such fluctuations cannot be easily estimated by averaging. Hence, the terminology, diffusion velocities BoldItalic^F-BoldItalic^S, BoldItalic^G-BoldItalic^F, with respect to the skeleton, or frequently appearing in soil mechanics, seepage velocity or velocity of filtration should be understood as macroscopic notions, and their interpretation may change from one process to the other. For instance, in cases of catastrophic phenomena, such as liquefaction, the filter velocity cannot be identified with any of the macroscopic average velocities. Macroscopic models based on the above notions must be then correspondingly extended.

At the first glance, the problem of averaging of velocities can be replaced by the averaging of momenta. These are volume densities, and consequently, one could write

(40)\begin{array}{l}{\rho }^{\mathrm{S}}{\mathit{v}}^{\mathrm{S}}={\displaystyle \frac{\mathrm{1}}{\mathrm{v}\mathrm{o}\mathrm{l}(\mathrm{R}\mathrm{E}\mathrm{V})}}{\displaystyle \underset{\mathrm{R}\mathrm{E}\mathrm{V}(\mathit{X},t)}{\int }{\rho }^{\mathrm{S}\mathrm{R}}(\mathit{Z},t){\mathit{v}}^{\mathrm{S}\mathrm{R}}(\mathit{Z},t){\chi }^{\mathrm{S}}(\mathit{Z},t)\mathrm{d}{V}_{\mathit{Z}}},\hfill \\ {\rho }^{\mathrm{F}}{\mathit{v}}^{\mathrm{F}}={\displaystyle \frac{\mathrm{1}}{\mathrm{v}\mathrm{o}\mathrm{l}(\mathrm{R}\mathrm{E}\mathrm{V})}}{\displaystyle \underset{\mathrm{R}\mathrm{E}\mathrm{V}(\mathit{X},t)}{\int }{\rho }^{\mathrm{F}\mathrm{R}}(\mathit{Z},t){\mathit{v}}^{\mathrm{F}\mathrm{R}}(\mathit{Z},t){\chi }^{\mathrm{F}}(\mathit{Z},t)\mathrm{d}{V}_{\mathit{Z}}},\hfill \\ {\rho }^{\mathrm{G}}{\mathit{v}}^{\mathrm{G}}={\displaystyle \frac{\mathrm{1}}{\mathrm{v}\mathrm{o}\mathrm{l}(\mathrm{R}\mathrm{E}\mathrm{V})}}{\displaystyle \underset{\mathrm{R}\mathrm{E}\mathrm{V}(\mathit{X},t)}{\int }{\rho }^{\mathrm{G}\mathrm{R}}(\mathit{Z},t){\mathit{v}}^{\mathrm{G}\mathrm{R}}(\mathit{Z},t){\chi }^{\mathrm{G}}(\mathit{Z},t)\mathrm{d}{V}_{\mathit{Z}}},\hfill \end{array}

where BoldItalic^SR, BoldItalic^FR, and BoldItalic^GR denote real (true) velocities in channels. However, this is not much of the help. First of all, the macroscopic momentum balance equations that yield field equations for BoldItalic^S, BoldItalic^F, and BoldItalic^G must contain source terms, such as diffusive forces reflecting interactions through walls of channels or frictional forces between fluid components. Such terms are not present in microstructural momentum balance equations and must be introduced by some additional surface or line integrals on REV, which, of course, has nothing to do with volume averaging. Second, average momentum densities do not describe strong deviations (fluctuations) of true momenta from average values anyway. These would have to be introduced additionally to averaged momentum balance equations. One of the microstructural variables designed for this purpose is the tortuosity described in Section 3.2. Consequently, in spite of a few attempts to develop this procedure, it seems to be easier to deal directly with purely macroscopic models.

The assumption on a common temperature of components, T, is related with difficulties of thermodynamic procedures for mixtures with multiple temperatures. Such procedures are scarce in literature and usually not directly applicable to porous media [36,37].

The field of porosity, n, has been discussed in detail in the last section about two-component models. In the followingsection, the notion of saturation, S, is discussed in some detail.

Capillary pressure, saturation, and retention curves

As mentioned earlier, the construction of a relation for the saturation as one of the fields for a three-component immiscible mixture is based on considerations concerning mechanical properties of the microstructure. They are still limited to isothermal conditions.

Between two immiscible fluids, of which one may be gaseous, e.g., water and air, a discontinuity in pressure exists across the interface separating them. The difference is called capillary pressure (compare Eq. (39)):

(41)p_{\mathrm{c}}=p^{\mathrm{G}\mathrm{R}}-p^{\mathrm{F}\mathrm{R}},

where p^GR is the true pressure in the nonwetting phase (gas in our case), and p^FR is the corresponding value in the wetting fluid (water in our case). The fluid phase whose molecules or atoms preferentially are adsorbed on a solid surface is called the wetting fluid, while the superseded material is denoted as a nonwetting fluid.

In a fluid-gas-mixture, the Young-Laplace equation describes the capillary pressure difference due to the phenomenon of surface tension σ and relates it to the radius of the bubbles r:

(42)p_{\mathrm{c}}=\mathrm{2}{\displaystyle \frac{\sigma }{r}}.

For a porous medium, the capillary pressure is a measure of the tendency to suck in the wetting fluid or to repel the nonwetting phase. In soil science, the negative of the capillary pressure (expressed as the pressure head) is called suction. The radius r then is of the order of magnitude of the pore or grain size. The capillary pressure, thus, depends on the geometry of the void space, on the nature of the solids and fluids and on the degree of saturation. As we have already pointed out, in soils, the geometry of the void space is extremely irregular and complex. Hence, only an idealized model may be adopted (e.g., capillary tubes, spheres of constant radius, or a bundle of parallel circular rods). Laboratory experiments are probably the only method to derive the relationship p_c = p_c(S).

In the experimentally determined capillary pressure curves a hysteresis occurs. This means that different capillary pressures may be obtained for a certain degree of saturation, depending on whether a sample is initially saturated with a wetting or with a nonwetting fluid. In both cases, the fluid initially saturating the sample is slowly displaced by the other fluid. When the sample is initially saturated with a wetting fluid, the process is called drainage, otherwise imbibition. Figure 4, which is taken from Ref. [17], shows a typical capillary pressure wetting fluid saturation relationship (kerosene and water in a sandstone) including the effect of hysteresis. In theoretical approaches, the hysteresis is neglected for simplification. The capillary pressure curve is known also under several other denotations, e.g., it is also called retention curve, pF-curve, or soil-water characteristic curve.

Inspection in Fig. 4 shows that a certain quantity of wetting fluid remains in the sample even at high capillary pressures. The value of the water saturation at this point (Bear denotes it by S_w0, otherwise, it is often described by S^rF-r stands for residual) is called irreducible saturation of the wetting fluid. Something similar appears if one looks at the imbibition curve. It is observed that at zero capillary pressure, there remains a certain amount of the nonwetting fluid—the residual saturation of the nonwetting fluid (Bear denotes it by S_nw0, otherwise, it is often described by S^rG). It indicates the amount of entrapped air in the pores. This leads to the introduction of the term effective saturation (for more information, see Ref. [38]):

(43)S_{\mathrm{e}}:={\displaystyle \frac{S-S^{\mathrm{r}\mathrm{F}}}{\mathrm{1}-S^{\mathrm{r}\mathrm{F}}-S^{\mathrm{r}\mathrm{G}}}}\mathrm{o}\mathrm{r}{S}_{\mathrm{e}}:={\displaystyle \frac{S-S^{\mathrm{r}\mathrm{F}}}{\mathrm{1}-S^{\mathrm{r}\mathrm{F}}}}.

The measured capillary pressure curves for several soil types differ considerably. Particularly, the range of the saturation, which can appear in a special soil, is rather distinct. A rough overview for sands, silts, and clays is given in Fig. 5, which is a modification of a figure in Ref. [39]. The entire curve, according to Ref. [40], can only be obtained by combining several measuring techniques because the range of capillary pressures can cover up to seven orders of magnitude, while each technique is applicable only up to three orders of magnitude.

There are some attempts to put the once experimentally measured curves into formulas. The most common approaches are those of Brooks & Corey [41] and van Genuchten [42]. According to the latter, the relationship is described by

(44)p_{\mathrm{c}}={\displaystyle \frac{\mathrm{1}}{{\alpha }_{\mathrm{v}\mathrm{G}}}}{[S_{\mathrm{e}}^{(-\mathrm{1}/m_{\mathrm{v}\mathrm{G}})}-\mathrm{1}]}^{\mathrm{1}/n_{\mathrm{v}\mathrm{G}}},

where α_vG, m_vG, and n_vG are parameters that depend on the type of the soil, and S_e is the effective water saturation (43). For simplicity, mostly, it is assumed that S_e ≡ S. Obviously, the van Genuchten equation is a nonlinear relation in the saturation S.

It is clear that the relation (44) can be incorporated in a macroscopic model only under the condition that we know how to transfer true pressures p^GR and p^FR to the macroscopic level. The simplest possibility is to use relations (38), but they may be too simplistic as we have already indicated. They may give reasonable results in linear models, but in general, the problem is still open.

Heat conduction

In Section 2, on one-component models, the Fourier law (9) has been already introduced. It was mentioned that the coefficient λ, the heat conductivity, is for soils heavily dependent on the morphology and that this dependence cannot be reflected by a one-component model. Therefore, we return here to this point in order to describe it by a three-component model. In Fig. 6 (compare Ref. [43]), we show an example of a nomogram in which the dependence of λ is demonstrated for various moisture contents (i.e., for various mass densities of the fluid if the mass density γ_d = ρ^S of the skeleton is fixed), saturations (i.e., the volume fraction of the gas to the fluid component), and mass densities of the skeleton. As we see, in this example, λ varies between 0.1 W/mK to 1.4 W/mK.

Unfortunately, the heat conductivity, λ, cannot be derived by means of any averaging procedure from microscopic conductivities of components. For this reason, we have to rely on empirical relations. These were proposed for soils around 30 years ago.

Recently developed experimental equipment, such as TP O2 probe, allow to make nonsteady-state measurements of heat conductivity [44] for various morphologies of soils. In a work, published in 2008, Chen [45] has proposed the following empirical relation for the conductivity:

(45)\lambda ={\lambda }_{\mathrm{0}}^{\mathrm{1}-n}{\lambda }_{\mathrm{w}}^n{[(\mathrm{1}-b)S+b]}^{cn},

where S is the saturation, λ₀ is the grain heat conductivity, λ_w = 0.61 W/mK is heat conductivity of water, and b and c are fitting parameters. For example, for sandy soils, λ₀ = 7.5 W/mK, b = 0.0022, and c = 0.78. In an implicit way, this relation accounts as well for a dependence on temperature through λ₀ and λ_w.

The results for the heat conductivity and corresponding theoretical one-component and multicomponent models play a particularly important role in description of freezing and frost heaving of soils.

Field equations for the above fields follow either from balance equations or from additional constitutive assumptions (saturation). Balance equations contain fluxes and sources and these, in turn, require additional information on material parameters.

Concluding remarks

Three important issues of theoretical modeling should be mentioned.

The first one is the formulation of boundary conditions. In one-component models, these are classical and extensively discussed in elasticity or plasticity. In multi-component models, the problem is more complicated because one has to formulate additional conditions for the extended set of partial differential equations. Even in the case of impermeable boundaries and such are phreatic surfaces of contact between saturated and dry domains of soils, one has to formulate an equation of motion of the surface itself. Such moving boundaries yield the boundary value problems with free boundaries, and these are usually ill-posed and create big mathematical problems. A physical presentation of this problem can be found in the book of Bear [17]. The situation is even worse on permeable surfaces. A part of the conditions on such surfaces has been formulated by von Terzaghi who had shown that the external loading must be taken over by the whole stress vector BoldItalic = σ_ikn_k, where BoldItalic = n_kBoldItalic_k is the unit normal vector of the boundary. The second condition was extensively discussed in the literature, and it concerns the flow through the boundary. This boundary condition for inviscid fluids relates the pressure difference and the velocity of flow through the surface. It contains an additional material parameter, the so-called surface permeability. It plays a very important role on contact surfaces between different layers saturated with water and on the external surface, which is the seepage face.

The second issue appears if the transition zones of not fully saturated soils appear. They are created, for example, by infiltration processes. In such processes one has to account for the capillary effects and an appropriate theoretical model must describe more than one fluid component. In the last few years, such models are being developed. The presentation of a linear three-component model with capillary effects in applications to poroacoustics can be found in the book of Albers [28].

The third important issue is the development of software for geotechnical engineers, which would account for all those theoretical problems that we have mentioned above. Such computational packages do not exist yet. It is only very recent that the research in this direction has been intensified. In particular, it concerns the formulation of some macroscopic constitutive laws in terms of microscopic material properties in which the micro-macro transition would be done in a numerical way. Mattsson et al. [46] formulate this problem in their extensive survey work from the year 2008 in the following manner: “The main objective with this literature survey is to elucidate the state of the art of internal erosion in embankment dams in order to be able to formulate a research program for numerical modeling of internal erosion in a physically sound manner.”

References

[1] Vaníček M, Vaníček I. Earth Structures in Transport, Water and Environmental Engineering (Geotechnical, Geological, and Earthquake Engineering). New York: Springer, 2008

[2] Bowen R M. Incompressible porous media models by use of the theory of mixtures. International Journal of Engineering Science, 1980, 18(9): 1129–1148

[3] Muir Wood D. Soil Behaviour and Critical State Soil Mechanics. Cambridge, UK: Cambridge University Press, 1991

[4] Muir Wood D. Geotechnical Modelling. Oxfordshire: Spon Press, 2004

[5] Lancellotta R. Geotechnical Engineering. Rotterdam: Balkema A. A., 1995

[6] Bauer E. Constitutive modeling of critical states in hypoplasticity. In: Pande, Pietruszczak, eds. Proceedings of the Fifth International Symposium on Numerical Models in Geomechanics. Rotterdam: Balkema A. A., 1995, 1520

[7] von Wolffersdorff P A. A hypoplastic relation for granular materials with a predefined limit state surface. Mechanics of Cohesive-frictional Materials, 1996, 1(3): 251–271

[8] Kolymbas D. A ratedependent constitutive equation for soils. Mechanics Research Communications, 1977, 4(6): 367–372

[9] Kolymbas D. Introduction to Hypoplasticity. Rotterdam: A. A. Balkema, 2000

[10] Bauer E, Tantono S F, Zhu Y, Liu S, Kast K. Modeling rheological properties of materials for rockfill dams. In: Zhu Y, eds, Long Time Effects and Seepage Behavior of Dams. Nanjing: Hohai University Press, 2008, 73–80

[11] Goodman M A, Cowin S C. A continuum theory for granular materials. Archive for Rational Mechanics and Analysis, 1972, 44(4): 249–266

[12] Wang Y, Hutter K. A constitutive model of multiphase mixtures and its applications in shearing flows of saturated solid-fluid mixtures. Granular Matter, 1999, 1(4): 163–181

[13] Wang Y, Hutter K. Shearing flows in a Goodman-Cowin type granular material-theory and numerical results. Particulate Science and Technology, 1999, 17(1-2): 97–124

[14] Kirchner N. Thermodynamics of Structured Granular Materials. Dissertation for the Doctoral Degree. Berichte aus der Thermodynamik. Aachen: Shaker Verlag, 2001

[15] Darcy H. Les Fontaines Publiques de la Ville de Dijon. Paris: Dalmont, 1856

[16] von Terzaghi K. Erdbaumechanik auf bodenphysikalischer Grundlage. Leipzig: Deuticke, 1925

[17] Bear J. Dynamics of Fluids in Porous Media. Dover: Dover Publication, 1972

[18] Forchheimer P. Wasserbewegung durch Boden. Z. Ver. Deutsch. Ing., 1901, 45: 1782–1788

[19] Wilhelm T, Wilmanski K. On the onset of flow instabilities in granular media due to porosity inhomogeneities. International Journal of Multiphase Flow, 2002, 28(12): 1929–1944

[20] Wilmanski K. A few remarks on Biot’s model and linear acoustics of poroelastic saturated materials. Soil Dynamics and Earthquake Engineering, 2006, 26(6-7): 509–536

[21] F. Gassmann; Über die Elastizität poröser Medien. In: Vierteljahresschrift der Naturforschenden Gesellschaft in Zürich 96. 1951, 1: 1–23

[22] Bowen R M. Compressible porous media models by use of the theory of mixtures. International Journal of Engineering Science, 1982, 20(6): 697–735

[23] Passman S L, Nunziato J W, Walsh E K. A theory of multiphase mixtures. In: Truesdell C, ed. Rational Thermodynamics, Berlin: Springer, 1984

[24] Faria S H, Hutter K, Kirchner N, Wang Y. Continuum Description of Granular Materials. Heidelberg: Springer, 2009

[25] Wilmanski K. A Thermodynamic model of compressible porous materials with the balance equation of porosity. Transport in Porous Media, 1998, 32(1): 21–47

[26] Wilmanski K. Thermomechanics of Continua. Heidelberg: Springer, 1998

[27] Wilmanski K. Continuum Thermodynamics. Part I: Foundations. Singapore: World Scientific, 2008

[28] Albers B. Modeling and Numerical Analysis of Wave Propagation in Saturated and Partially Saturated Porous Media. Aachen: Shaker Verlag, 2009

[29] Wilmanski K. On microstructural tests for poroelastic materials and corresponding Gassmanntype relations. Geotechnique, 2004, 54(9): 593–603

[30] Biot M A, Willis D G. The elastic coefficients of the theory of consolidation. Journal of Applied Mechanics, 1957, 24: 594–601

[31] Epstein N. On tortuosity and the tortuosity factor in ow and diffusion through porous media. Chemical Engineering Science, 1989, 44(3): 777–779

[32] Kozeny J. Über kapillare Leitung des Wassers im Boden (Aufstieg, Versickerung und Anwendung auf die Bewässerung). Sber. Akad. Wiss. Wien, 1927, 136 (Abt. IIa): 271–306

[33] Biot M A. Theory of propagation of elastic waves in a fluid saturated porous solid. I. Low frequency range. Journal of the Acoustical Society of America, 1956, 28(2): 168–178

[34] Wilmanski K. Tortuosity and objective relative accelerations in the theory of porous materials. Proceedings of the Royal Society, Mathematical, Physical and Engineering Sciences, 2005, 461(2057): 1533–1561

[35] Wesselingh J A, Krishna R. Mass Transfer in Multicomponent Mixtures. Delft: Delft University Press, 2000

[36] Ruggeri T, Simič S. Mixture of gases with multi-temperature: Maxwellian iteration. In: Ruggeri T, Sammartino M, eds. Asymptotic Methods in Nonlinear Wave Phenomena. New Jersey: World Scientific, 2007, 186–194

[37] Leşan D. On a theory of micromorphic elastic solids with microtemperatures. Journal of Thermal Stresses, 2001, 24(8): 737–752

[38] Helmig R. Multiphase Flow and Transport Processes in the Subsurface. Heidelberg: Springer, 1997

[39] Hartge K H, Horn R. Einführung in die Bodenphysik. Stuttgart: Schweizerbart, 1999

[40] Schick P. Ein quantitatives Zwei-Komponenten-Modell der PorenwasserBindekräfte in teilgesättigten Böden. Habilitationthesis. München: Universität der Bundeswehr, 2003

[41] Brooks R H, Corey A T. Hydraulic properties of porous media. In: Hydrology Papers, Vol. 3. Fort Collins: Colorado State University, 1964

[42] van Genuchten M Th. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Science Society of America, 1980, 44: 892–898

[43] Farouki O T. Ground thermal properties. In: Krzewinski T G, Tart R G, eds. Thermal Design Considerations in Frozen Ground Engineering. New York: ASCE, 1985: 186–203

[44] Gontaszewska A. Thermophysical Properties of Soils in Vicinity of Zielona Gora in Relation to Soil Frost Depth. Dissertation for the Doctoral Degree. Poznan: University of Poznan, 2006

[45] Chen S X. Thermal conductivity of sands. Heat and Mass Transfer, 2008, 44(10): 1241–1246

[46] Mattsson H, Hellström J G I, Lundström S. On internal erosion in embankment dams: a literature survey of the phenomenon and the prospect to model it numerically. Report No. 2008:I14, 2008

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag Berlin Heidelberg