1923–2023: One Century since Formulation of the Effective Stress Principle, the Consolidation Theory and Fluid–Porous-Solid Interaction Models

Abstract

In 1923, Karl Terzaghi developed the theory of soil consolidation in which he introduced the concept of effective stress (ES). Over the past century, various theoretical aspects have been unraveled regarding the Effective Stress Principle (ESP) and the fluid–porous-medium interaction in deformable permeable media; nevertheless, some aspects have been debated for a long time, and some perplexities are still perceived among scientists and professionals. By way of example, in the study of flow in deformable permeable media, particularly in fractured porous systems, some problems are still open. This review is aimed at providing an overview of the progress achieved over the past century in the theoretical and experimental treatment of ESP—with particular reference to saturated porous media—and of the geomechanical aspects of fluid flow and fluid–rock interaction, trying to answer to some common questions among professionals, such as what is the correct expression for the ES to be used in applications and why there are various formulations? Additionally, we try to answer questions related to the modeling of fluid flow in fractured porous media. Therefore, this review paper is divided into two main sections, “Effective Stress Principle” and “Fluid Flow, Consolidation, and Fluid–Rock Interaction”. In the first section, the basic concepts and the theory underlying the ESP are preliminarily illustrated, with a simple but rigorous theoretical proof, and, subsequently, historical remarks are provided. The second illustrates the different adopted theoretical approaches to fluid flow, starting from Terzaghi’s theory of one-dimensional consolidation up to the recent dual- and multiple-porosity models.

1. Introduction

In 1936, Karl Terzaghi enunciated the Effective Stress Principle (ESP) at the First International Conference on Soil Mechanics at Harvard University [1]. However, he provided the first application of this concept as early as 1923, as part of his studies on the phenomenon of the consolidation of deformable permeable porous media [2]. This work represents the starting point of practically all the studies concerning both the role of effective stress (ES) and the consolidation phenomena and, in general, the fluid–rock interaction underlying the transient hydraulic behavior of deformable permeable porous media.

Terzaghi proposed the following expression for the ES in one-dimensional form for the normal stress:

σ′ = σ − p

(1)

where σ′ denotes the effective stress, σ the total stress, and p the pore water pressure or neutral stress. In terms of stress tensor components, the equivalent form is the following:

σ′_ij = σ_ij − p δ_ij

(2)

where δ_ij denotes the Kronecker delta, and the usual convention used by geologists of a positive sign for compressive stresses is adopted (here and in the following text). According to Terzaghi, only the fraction of the total stress, σ, exceeding the pore-water pressure, p, which he called effective stress, σ′, is responsible for soil and rock strain and/or failure. In other words, pore-water pressure (or neutral stress) magnitude changes, at the same effective stress conditions, have practically no influence on stress–strain and/or strength behavior of saturated porous media.

He introduced this principle as an experimental law, and over the past century, several underlying theoretical facets have been explored. Nevertheless, some aspects remained controversial for a long time before being clarified, and some questions are still open. Some perplexities are still perceived by engineers, geologists, and geophysicists who currently apply ESP but who are not specialized in the complex theories of porous media and poroelasticity. A common question is, for example, why a porous medium follows different ES laws for strength and stress–strain behavior or why different media, e.g. rock, soil and metal powders, follow the same ES law for strength.

Several ES laws have been proposed by various authors in the following formulation (in one-dimensional form):

σ′ = σ − η p

(3)

where η is a parameter which generally depends on the properties of the porous medium, as well as on σ and p values. Note that, in Equation (3), the linear dependence of σ′ on σ and p is only apparent since, in general, the term η can be a strongly nonlinear function of these variables.

Different formulations have been proposed over time for the coefficient η. Some of them are as follows [3,4,5,6,7,8].

For the stress–strain behavior of soil and/or rock, we have the following:

η = 1; (Terzaghi)

σ′ = σ − n p; (Fillunger)

(4)

σ′ = σ − (1 − a_c) p; (Skempton)

(5)

σ′ = σ − (1 − K/K_s) p; (Skempton, Biot, Nur, and Byerlee)

(6)

σ′ = σ − (1 − n_s K/K_s) p; (Suklje)

(7)

For strength behavior, we have the following:

η = 1; (Terzaghi)

σ′ = σ − (1 − a_c tan ψ/tan ϕ′) p; (Skempton)

(8)

where n denotes the porosity, n_s = (1 − n) is the solid fraction; a_c denotes the contact area ratio (i.e., the area of contact between the particles, per unit gross area of the material); ϕ is the angle of shearing resistance; ψ is the angle of intrinsic friction [5]; and K_s and K are the bulk moduli of the constituent material (i.e., nonporous, non-fractured rock/clast) and of rock or soil aggregate, respectively.

Furthermore, various physical interpretations and definitions were attributed to the ES, and this may have caused significant confusion about its meaning. The main interpretations of the ES were (i) that of equivalent stress, dependent on the applied stress, σ, and pore pressure, p, capable of producing the same effect (strain, failure, etc.) on the same porous medium in dry conditions; (ii) the mean stress acting over the solid phase alone; and (iii) the intergranular stress (as discussed by Skempton, 1960, [5]). Definition (iii) had similar aspects to (ii), although it related only to granular media and led to a different formulation [5]. The currently most widely adopted definition is (i); nevertheless, this is substantially conventional (see Refs. [9,10,11]).

Regarding the theory of consolidation, the one-dimensional model developed by Terzaghi [2,12,13] is still currently used in civil-engineering design. It represents the first simple model of fluid–rock interaction in deformable permeable porous media. Later, more realistic models were developed, giving rise to autonomous disciplines such as poroelasticity and flow models in fractured media. In the study of flow processes in deformable permeable media, one may be mainly interested in the evolution of internal stress and deformations, e.g., in applications such as foundation design, subsidence, etc., or mainly interested in the quantification of fluid circulation and flow rates, e.g., in the field of oil recovery, diffusion of underground pollutants, and various hydrogeological or environmental studies involving ground fluids.

In the field of fluid-flow modeling in deformable permeable porous media, several questions are still open, particularly in studies involving fractured rock. This review is mainly aimed to provide an overview of different proposed models and theoretical approaches to the ESP and to flow modeling, over the past century. Given the vast literature concerning ESP in unsaturated soils, this topic has not been dealt with in this paper, and the reader is directed to the recent and comprehensive review by Khalili et al. [14] or to other relevant papers [15,16,17]. Regarding fluid-flow modeling and fluid–rock interaction in permeable deformable porous and/or fractured media, the main theoretical approaches are discussed, with particular reference to its study at several scales of observation and dual- and multiple-porosity models.

2. Effective Stress Principle

2.1. Definition of ES and Attributed Physical Meaning

Several physical interpretations and definitions have been attributed to the ES, and this may have caused significant confusion about its meaning. The currently most adopted definition of ES is well described, for example, by Lade and De Boer [18], who defined the ES as the stress, depending on the applied tension, σ_ij, and pore pressure, p, which controls the strain or strength behavior of soil and rock (or a generic porous body) for whatever pore-pressure value. In other words, this can be defined as the stress which, when applied over a dry porous body (i.e., at p = 0), provides the same strain or strength behavior which is observed at p ≠ 0 [9]. The usefulness of an appropriate ESP formulation consists in allowing us to predict the behavior of a porous body for whatever pore-pressure value on the basis of experiments carried out at a known pore-pressure value (e.g., at p = 0, in case of dry samples).

Other definitions of ES that have been adopted in the past were that of mean stress in the solid phase (which might be viewed as stress transmitted through the solid skeleton) and of intergranular stress (for granular media, see Reference Skempton [5]). Although these two definitions were conceptually similar, they led to very different formulations for the ES. In particular, the first was associated with the expression in Equation (4) and the second to Equation (5).

The above given ES definition may be extended to whatever property of a porous medium, e.g., seismic velocities, density, permeability, etc.; that is, the ES can be defined as the external stress that, if applied on a dry sample, would produce the same effect on that property (porosity, permeability, etc.) as a combination of tension, σ_ij, and pore pressure, p (see Refs. [10,11]). Therefore, it should always be taken into account that different properties are associated with different ES laws. A thorough and interesting discussion about the conventionality of ES definition was provided by Robin [10], who pointed out that such a definition may be, in some cases, misleading or also useless, whereas it may be useful and of practical use in other (frequent) cases, when it reduces to the simple Terzaghi’s formulation. In their recent paper, Hampton and Boitnott [11], based also on Robin’s studies, denote the term “effective stress” as a misnomer in order to emphasize how it may have the misleading connotation of equivalent stress, and they recommend its use only with specific reference to a particular property. The topic discussed in this section has been thoroughly illustrated by Zimmerman [19]; theoretical derivations about different ES laws, which can be associated with diverse kinds of elastic processes, are provided.

2.2. A Simple Example Aimed to Clarify the Conventional Nature of ES Concept

Let us imagine carryingout an unjacketed test on a porous sample made up of an isotropic linear elastic material, such as that described by Cheng [20] (Section 1.2.4) and also Refs. [5,6,10], in which pore pressure, P, and confining stress, Pc, are raised by the same amount. Although the constituting material is isotropic linear elastic, the aggregate, which can be a multi-connected solid including interconnected pores or constituted by grains in contact, is, in general, nonlinear elastic or inelastic, according to voids (pores and/or fractures) spatial distribution, orientation, and shape. In such a condition, this sample is transformed into a geometrically similar smaller one (cubic strain). According to the above given definition, the ES is the stress which, when applied on the dry sample, provides the same here-observed volume strain, i.e., for which C Pc’ = C_s p, where Pc’ denotes the confining ES, and C and C_s denote the compressibility of the aggregate and of the constituting material, respectively.

It should be noted that, although the sample in dry condition undergoing such stress exhibits the same volume change observed in the unjacketed test, these two strain processes are geometrically and mechanically substantially different. The unjacketed test produces a cubic volume strain with no pore-shape change, whereas the dry sample shows volume strain involving relevant pore-geometry change, pore collapse, etc.

By means of this simple example, we can understand that the ES definition is purely conventional, and the choice between more definitions is substantially dictated by the opportunity in a certain work context. By way of example, in the case of volume-variation studies of soils, rocks, or other porous materials, the choice of the most appropriate definition depends on the characteristics we are interested in. If we are interested in comparing overall volume changes associated with different pore pressures (e.g., in studying large-scale subsidence in different stress scenarios), then the most suitable law is the Skempton–Nur–Byerlee’s (Equation (6)). If, on the other hand, we are interested in volume variations associated with variations in pore geometry (which, for the same constituent material, controls all mechanical and hydraulic properties of the porous medium), then the Terzaghi formulation (Equation (1)) provides excellent approximations for the ES and is preferable to that of Skempton–Nur–Byerlee.

2.3. What Do Rock, Soil and Other Porous Media Have in Common from the Viewpoint of Mechanical Behavior? Strain Mechanics of Porous Media at Several Scales of Observation

Rock, soil, and other granular materials (e.g., metallic powders) show different mechanical behaviors which are associated with different micromechanics. Nevertheless, in terms of strength behavior, they have in common that this is mainly controlled by concentrated stresses in certain areas of the material at different scales of observation, such as fracture tips (e.g., rock), grain contact area (soil), or crystal lattice dislocations (e.g., metallic aggregates). Here, the term “concentrated stress” denotes a local stress tensor, within the solid material, in which at least one component is greater than all the components of the stress tensor acting on the surface of representative volume (viewed as a remote stress) by any order of magnitude.

2.3.1. Rock

The stress–strain behavior of shallow rocks is ordinarily characterized by nonlinear elasticity, although the constituent material may often display linear elastic behavior over a wide range of stress values. The rock macroscopic strain is mainly controlled by the local strain of elongated pores and/or fractures and pore collapse [6,19,21,22,23,24,25,26,27,28,29,30,31,32].

The strength behavior of rock is associated with the concurrent phenomena of brittle and ductile strain in different proportions, according to several variables (thermodynamics, strain rate, chemical environment, etc.). The wide literature on the field observations and theoretical considerations points out that fractures in rock start at cavities and microcracks (see Refs. [33,34,35,36,37,38,39,40,41] and many others). Tensional and shear fracture growth are mainly affected by tensions acting in areas of high stress concentration, such as fracture tips. In the framework of brittle strain, both mode I and II or III fractures, spread when the local stress is sufficient to produce detachment between adjacent crystals or between two portions of the same crystal associated to bond breaking at the crystal lattice scale. According to the leading models of crack growth of Griffith [42,43,44,45] and of Barenblatt (see Ref. [46]), as well as of Nonlocal Continuum Field Theory [47,48], such a phenomenon is characterized by stresses at the fracture tips (microscale stress) greater than the remote stress (that acting over a representative volume or a sample of rock) of several orders of magnitude. For the kind of glass used in his experiments, Griffith found a theoretical tensional strength of ca. 11 GPa (1.6 × 106 psi in his original paper). The theoretical tensile strength of common crystals, such as quartz or calcite, is of the same order of magnitude of that of Griffith’s glass (ca. 10 GPa), whereas tensile rock strength commonly falls in the range 1–5 Mpa. This is also valid for crystal plasticity because, also in correspondence of dislocations, there are considerable stress concentrations which are responsible for various phenomena of dislocation creep, sliding, etc., as evidenced by Burgers [49] and references therein, as well as later theoretically discussed by Eringen (Refs. [47,48] and Ref. [50] Section 6), through criteria of the Nonlocal Continuum Field Theory. Therefore, rock strength behavior is always controlled by concentrated stresses at fracture tips or, at crystal lattice scale, in correspondence of dislocations.

2.3.2. Soil

Soils, unlike rocks, always exhibit inelastic behavior. The differentiation between stress–strain and strength behavior for soils is substantially conventional, as the strain process is always irreversible. Furthermore, as observed by Lade and de Boer [18], a variety of porous media exists in nature with intermediate characteristics, e.g., pyroclastites, weakly cemented sediments, etc., for which rock and soil represent two extreme situations of a continuum of possible cases.

Soil strain is a consequence of reciprocal displacements between component particles, which depend on their elastic strain both in terms of small variation in shape and of variation of forces on the contact areas. When the local tangential force exceeds the local friction on contact between two particles, producing sliding, and then the system reorganizes itself, reaching a new spatial arrangement. The reciprocal displacement between grains (which are associated with the macroscopic strain of aggregate) can occur in the absence of significant grain deformations, as in the case of soils under load conditions in usual geotechnical applications, or accompanied by grain fracturing (e.g., experiment on gypsum sand under high pressure by Lade and de Boer; see Ref. [18]). Secondarily, strain of grains may be associated with plastic strain. In all of these cases, the soil strain is always strictly dependent on concentrated stresses, at the contact areas between particles or at fracture tips (within grains), or in correspondence of dislocations (at crystal lattice scale).

2.3.3. Metallic Powders and Lead Shots

The stress–strain and strength behaviors of metallic powders and granular media made up of lead shots (which were involved in a main experiment by Skempton [5]) are mainly associated to ductile strain of particles. Therefore, these are substantially controlled by concentrated stresses in correspondence of dislocations at the crystal lattice scale. It is commonly assumed that there is no significant interaction between crystal plasticity phenomena and local hydrostatic pressure and namely that (1) pressure variations do not significantly hinder/favor phenomena, such as dislocation sliding; and (2) that, following or simultaneously with dislocation gliding, the material always reacts to local (i.e., at a scale smaller than grains) hydrostatic pressure in the same way, i.e., according to elastic behavior. The hypothesis of point (1) is guaranteed, as metals exhibit negligible frictional behavior showing intrinsic friction angle values near to zero (e.g., see Ref. [5]). The hypothesis (2) implies an approximation, as local plasticization can lead to local anisotropies and heterogeneity within single clasts, due to dislocation propagation/extension/formation. However, it appears reasonable and in agreement with the experimental evidence in which nonporous crystalline aggregates, which do not show significant iso-orientations of crystals, show macroscopic linear elastic isotropic response at hydrostatic pressures, although local deviatoric stresses may occur at the interface between neighboring crystals. Anyway, the assumption that (non-porous) particles behave elastically in hydrostatic compression and that dislocation-related and viscoelastic effects occur only in shear is commonly adopted (e.g., see Ref. [51], Chapter 9).

2.4. Theory

2.4.1. Problem Setting

The here-studied porous-medium model consists of a solid with interconnected voids or a granular medium, i.e., consisting of elastic bodies in contact. In Ref. [9], it was illustrated how the elastic analysis involving a porous elastic body with interconnected voids can be extended to porous granular media. A basic here-adopted assumption consists in considering the porous medium as constituted by a linear elastic, homogeneous, and isotropic material. It should be noted that, despite such an assumption, the porous aggregate can be nonlinear elastic (e.g., rock) or inelastic (e.g., soil), anisotropic according to pore and fracture orientation, heterogeneous according to void spatial distribution (i.e., as a consequence of spatially variable porosity). The suitability of such assumption was discussed in detail by Zimmerman et al. [19,22,23,51]. Therefore, although this hypothesis reduces the general validity of such model, this model covers, with good approximation, a wide variety of real cases and shows the benefit of a relatively simple theoretical treatment.

To study the characteristics of the stress field within a porous medium, an analysis on a smaller scale than pores and/or clasts is here carried out. Specifically, here a representative volume of porous medium, including many pores and voids, is considered. The elastic problem within this representative volume will involve elementary volumes which are smaller than pores and/or clasts. Here, σ_ij denotes the stress tensor acting on a cubic representative volume of the two-phase medium; for example, along a certain axis, the normal component, σ, is given by the ratio between the force, F, and the surface area, A, orthogonal to this axis, where F is the force acting on such surface, which is resisted by the solid and the liquid phases, i.e., F = F_solid + F_fluid (Figure 1). Denoted by σ_s the stress acting on the solid fraction only, i.e., σ_s = F_solid/(n_s A), results in F = σ_s n_s A + p n A (the symbols associated with the different variables are illustrated in Figure 1). It should be noted that the term σ − p = F/A − p can be put in the form σ − p = (σ_s − p) n_s:

σ − p = σ_s n_s + p n − p = σ_s n_s − p ∙ (1 – n) = (σ_s − p) n_s

(9)

Therefore, the term (σ − p) is proportional to (σ_s − p); that is, they are in one-to-one correspondence. This consideration will be useful in the following sections because, in considering the boundary conditions of the solid phase only (and not of the two-phase medium), the corresponding term (σ − p) is employed in the equations instead of (σ_s − p). This choice is due to the fact that the experimentally measured quantity, e.g., in triaxial tests, is the force, F, and consequently, given the specimen section area, A, the stress σ and not the σ_s.

The stress acting on an elementary volume smaller than pores is denoted by ψ_ij (Figure 2). Therefore, σ_sij (and so σ_ij) can be viewed as a boundary condition for the elastic problem having as an unknown function the stress field, ψ_ij. The boundary of the considered continuum is a surface, Σ, which is composed by the external boundary surface, Σ₁, and the pore surface, Σ₂ (Figure 2).

2.4.2. Equations of Elasticity

Here, some recall about the well-known elastic problem and related equations, is provided. The stress field inside a linear elastic body in absence of body forces is governed by the following linear elliptic differential system (see Ref. [42]):

(1 + ν) ∇²ψ_ij + ∂²ψ_kk/∂x_i∂x_j = 0

(10)

where ν is the Poisson ratio, and ∇² is the Laplace operator. Furthermore, the usual Einstein convention about repeated indexes is adopted (e.g., ψ_kk = ψ₁ + ψ₂ + ψ₃).

The boundary conditions are defined by the external stress acting over the (deformed) boundary surface, Σ, i.e., over the outer surface (Σ₁) and pore surface (Σ₂):

Ψ_ij n_j = t_i

(11)

where n_j denotes the local boundary surface normal vector, and t_i denotes the tension acting on the boundary at a point. Specifically, the solid phase of a representative volume is subjected to a stress σs_ij acting over Σ₁ and a pressure, p, acting on Σ₂. The pore pressure at the moment is assumed here to be uniform within the representative volume.

The equation system (10) provides the basic relations of the Theory of Elasticity and is derived by combining local equilibrium equations with the compatibility conditions between infinitesimal strain components, as well as the elastic relations, as illustrated by Refs. [24,46].

The equilibrium conditions of elementary volumes in static case and in absence of body forces (whose role in our treatment can be neglected) are given by the following three-equation system:

∂ ψ_ij/∂ x_j = 0

(12)

The term on the left side of each equation represents one component, along the x_j axis, of the resulting force acting over the surface of a generic elementary volume.

The compatibility equations express the condition that strain components cannot be taken arbitrarily as functions of x_j, as these are a function of only three independent variables, i.e., the displacement u_j [42]. According to the hypothesis of linear elasticity, the strain components are linear functions of those of stress.

The achieved system of six elliptic differential equations, together with appropriate boundary conditions, allows us to determine the unknown stress field, ψ_jk(x). As an alternative, the system of equilibrium Equation (10) can be expressed in terms of strain and then of displacement components. This criterion allows us to obtain a system of three differential equations, having as unknown functions the three displacement components, u_j(x), and boundary conditions in terms of displacements. The choice between this latter system and that given by Equation (10) depends on the kind of boundary conditions, i.e., whether these are in terms of imposed external forces or imposed external displacements. In this study, the former kind of conditions is adopted, being primarily interested in studying the relation between the internal local stress field, ψ_jk(x); the (representative volume) boundary stress, σ_jk; and the pore pressure, p.

2.4.3. Theoretical Proof of the ESP

A simple theoretical proof of the ESP, based on the theory of elasticity, can be carried out according to the following steps: (1) decomposition of the boundary stress, σ_ij, in two contributions, one proportional to (σ − p) and another given by pore pressure, p; (2) assessment of local stress field features which are associated to these boundary stresses; (3) proof of applicabilityconditions of the superposition principle even when the porous medium shows nonlinear elastic or inelastic behavior (Appendix A); (4) determination of the local stress field and its effect on macroscopic strain; and (5) concluding reflections about the role of stresses (σ − p) and p in stress–strain and strength behaviors of different porous media.

• The stress, σ_ij, acting over the boundary surface of a representative volume of the two-phase medium, can be decomposed as follows (Figure 3):

  σ_ij = (σ_ij – p δ_ij) + p δ_ij

  (13)

This is equivalent to the following decomposition, involving the only solid phase, whose applied stress on the outer surface, Σ₁, is σs_ij:

σs_ij = (σs_ij – p δ_ij) + p δ_ij

which, for Equation (9), can be written as follows:

σs_ij = (σ_ij−pδ_ij)/n_s + p δ_ij

(14)

This is only a simple algebraic operation which is always allowed in mechanics (with no generality loss); analogous decompositions are widely adopted in Continuum Mechanics and Theory of Elasticity (e.g., in Ref. [42], in several theoretical proofs along the whole textbook) in order to decompose the boundary conditions in more simple systems whose associate solutions are already known or easily manageable.

2.

The stress field associated with the second term on the right side of Equation (14), (p δ_ij), which is now applied over Σ₁ and Σ₂, i.e., the whole boundary surface, is a well-known solution consisting of a uniform isotropic pressure equal to p (see Ref. [42]) for whatever boundary or pore surface shape, whereas the stress at a point associated with the first term (σ_ij − p δ_ij) is an unknown tensor, ψ0_ij = ψ0_ij(σ_kl − p δ_kl), which depends only on the differential stress (σ_kl − p δ_kl) (Figure 3). We are not interested in calculating the solution ψ_ij for each point but at achieving information about its features, namely at comparing the local stress and strain fields in a porous medium subjected to a pore pressure, with those into an identical one in dry condition.

3.

Although the superposition principle is generally not valid in the case of nonlinear elastic material (specifically, nonlinear elasticity of the model defined in Section 2.4.1. is related to a varying boundary pore surface), it can be proved that, in the special case in which one of the two stress systems consists of a uniform hydrostatic pressure, applied to the whole boundary, i.e., on the external and pore surfaces, such a principle is applicable. Preliminary, it is recalled that two geometrically similar elastic bodies, undergoing the same external stress, show the same local stress at corresponding points; in other words, the system of Equation (10) provides solutions that depend only on boundary surface shape and not on its size (see Ref. [9]).Let us consider a sample of porous medium, with an undeformed boundary surface ΣI. Imagine applying a stress σ_a over the outer surface and then subsequently unloading the sample and applying a pressure, p, on the whole boundary surface. Let ΣII denote the deformed surface after application of σ_a and by ΣIII that was obtained after the application of p alone (Figure 3). Let ψ_a(x) denote the stress field associated with σ_a and by ψ_p(x), which is associated with p. The solution associated with the latter boundary condition is the well-known constant isotropic tensor, ψ_pij(x) = p δ_ij. In this case, as the pressure, p, produces a cubic strain (which implies geometric similarity), it follows that the surface ΣIV is geometrically similar to ΣII, i.e., the one that would be obtained for p = 0. Therefore, the solution ψ_a(x) fulfils the boundary conditions (11) also over ΣIV. Taking into account that the solution ψ_pij = p δ_ij fulfils the boundary condition relating to a uniform pressure for whatever boundary shape, we obtain the following:

ψ_ij(σ_akl + p δ_kl) = ψ_aij(σ_akl) + p δ_ij

(15)

Moreover, both terms on the right side fulfil the boundary conditions on ΣIV.

4.

The solution of Equation (10), associated with the boundary stress of Equation (13), is given by the superposition principle, by the sum of ψ0_ij and p:

ψ_ij = ψ0_ij (σ_kl − p δ_kl) + p δ_ij

(16)

Equation (16) is not an ES law but rather a representation of the local stress field on a smaller scale than pores. This plays a key role in the here-developed model, as it allows us to evaluate the relevant properties of local stresses within the material and can clarify several significant aspects of ESP. The stress field inside the material that constitutes a saturated porous body includes the following:

• A term (p δ_kl) which produces only a small volume reduction (cubic strain), in accordance with the intrinsic bulk modulus of the solid, K_s, with no pore-shape change;
• A term ψ0ij, which depends on (σkl – p δkl) and which produces the often-more-evident strain at the macroscopic scale, in terms of aggregate volume and shape change, as well as pore surface shape deformation. It should be noted that this stress and strain depend exclusively on the differential stress (σkl − p δkl), i.e., Terzaghi’s ES. This provides a first explanation of the ESP for stress–strain behavior of rock. Two identical porous bodies, one subjected to pore pressure and one in dry condition but undergoing to the same differential stress (σkl − p δkl), show macroscopic strains which differ only by a small volume change associated with p.

5.

The role of the differential stress (σkl − p δkl) and pressure, p, in stress–strain and strength behavior of different porous media can be depicted as follows:

• Stress–strain behavior: The differential stress (σ − p) is responsible for the variation in the geometry of the pores and/or grain spatial arrangement, which, with the same material, controls all the mechanical and hydraulic characteristics of the porous medium itself. The pressure (p) produces a volume strain which, in many geotechnical engineering applications, is often modest or negligible but may be significant in low-porosity rocks subjected to high pressures in the Earth’s crust. Such a strain is associated with a geometrically similar transformation of the porous medium geometry with no pore/void shape change. In the study of stress–strain behavior of soils and rocks (and a great variety of other porous media), this volume variation is often taken into account by means of the Equation (6), as it is measured together with that due to the other stresses acting on the sample or on the generic representative volume.
• Strength (and inelastic) behavior: All of the above illustrated inelastic deformations (fracture in rocks and concrete, reciprocal displacements between particles in soils, ductile strain of grains in metal powders or lead shots) are practically exclusively controlled by forces acting on areas of high stress concentration, such as fracture tips in rock, contact areas in soils, and dislocations at the crystal lattice scale within metallic grains or lead spheres, in which pore pressure is negligible with respect to the maximum local stresstensor component, denoted by ψ0_Max:

  ψ0_Max >> p

  (17)
  This implies that, for the maximum local stress component, Equation (16) becomes as follows:

  ψ_Max ≈ ψ0_Max(σ_kl − p δ_kl)

  (18)
  Therefore, in areas of stress concentration, the maximum local stress component fully depends on Terzaghi’s ES. As a consequence, any variation in the geometry of the cracks and fractures, as well as in the arrangement of the particles in porous granular media, is controlled exclusively by Terzaghi ES (Equation (1)).
• Clarification on stress–strain behavior of granular media: The above drawn considerations about the strength behavior of porous media are also valid in those deformations in soils and granular media that are often improperly considered nonlinear elastic, i.e., those in the pre-failure phase, as they are in fact inelastic. The trend of a strain process and the reciprocal displacements between particles are controlled, with excellent approximation, exclusively by Terzaghi’s ES, (σ − p) (Equation (1)), even in the case of high pressures. However, in this latter case, it is often considered convenient to use the Skempton–Nur–Byerlee expression of the ES (Equation (6)), rather than that of Terzaghi (Equation (1)), in order to also take into account the volume variation due to the above considered term, p, as it is measured experimentally together with the volume variation due to particle rearrangement. Anyway, it should be taken into account that this choice is purely conventional (Section 2.2).

2.4.4. Choice between ESP Formulations

Numerous formulations of the ES have been provided in the past; some of them, e.g., Equations (4) and (7), refer to definitions of ES different from the one illustrated in Section 2.1., whereas others are in contrast with experimental data. By way of example, Equation (4) proposed by Fillunger [3,4] does not refer to the ES definition, as given in Section 2.1., but to the stress transmitted by the solid phase, and its use can be useful in balance or equilibrium equations [9]. Instead, for a long time, it was commonly believed that ES was associated with intergranular stress and its formulation was that in Equation (5); however, Skempton [5] proved experimentally (by means of oedometer test on lead shots) that it was wrong. In the same work [5], he proved that the formulation he proposed for ES for strength behavior, Equation (8)—once plausible values were entered for the parameters a_c and tan ψ, both for rocks and for metal powders—gave ES values experimentally indistinguishable from those obtained with the simple Terzaghi formulation (Equation (1)).

Therefore, among the various formulations illustrated in the Introduction, only those of Equations (1) and (6) are valid or useful. Experimental observations have always shown that, in studying the strength behavior of rocks, soils, and other saturated porous materials, Terzaghi’s formulation provides excellent agreement with the observed data. This is explained by the considerations made above in this section. Therefore, the correct formulation of the ES for strength behavior is that provided by Terzaghi, as shown inEquation (1).

The correct ES formulation to study the stress–strain behavior of rocks and granular media is that provided by Equation (6), as it also allows us to take into account the volume reduction associated with the pore pressure and is in agreement with the ES definition given in Section 2.1. However, it must be taken into account that Equations (6) and (1) provide different values only in the case of high pressures (e.g., in the deep Earth’s crust), whereas, for almost all geotechnical applications, these equations yield quite identical values.

2.5. Historical Remarks and Overview of Several Theoretical Approaches

2.5.1. The Early Works

The first formulation of ESP, as currently interpreted (i.e., according to definition in Section 2.1), was introduced in 1923 by Karl Terzaghi in a study on the consolidation of clay layers [2], and about a decade later, it was presented at the First International Conference on Soil Mechanics at Harvard [1]. The effect of pore pressure on the strength of materials had already been studied a few years earlier by Paul Fillunger(in 1915) [4] and Rudeloff and Panzerbieter(in 1912) [52]; however, these studies led to controversial results. According to Rudeloff and Panzerbieter, pore pressure showed a significant influence on the strength of porous media, whereas Fillunger sensed that only a fraction of the total stress controlled the mechanical behavior of such materials. Fillunger [3] proposed a formulation of the ES according to Equation (4), which was, however, associated with a different definition of ES from that currently shared (given in Section 2.1) that did not refer to an experimentally measurable quantity but rather to a variable (stress transmitted through the solid phase) which was useful in the mathematical models he developed in the field of soil consolidation theory [9]. Because the theory of soil and rock mechanics was something newat that time, it was not easy to interpret the subtle difference between stress effectively transmitted by the solid skeleton of a porous medium and ES as defined in Section 2.1. Terzaghi was more oriented to an experimental approach than Fillunger, and his definition of ESP was (and is) compatible with that given in Section 2.1, being aimed at predicting the behavior of saturated soils and rocks on the basis of experimental tests.

According to Skempton [53], who commented as follows, although Fillunger and other scientists had introduced the concept of ES in studying porous saturated materials, only Terzaghi has gasped the extent of this fundamental law: “Fillunger naturally concluded that the tensile strength does not vary with water pressure, at least within the pressure range and limits of accuracy of his experiments. This amounts to a corollary of the principle of effective stress, in the special case under consideration, yet neither he nor anyone else at the time realised the significance of these results”. He continued to say, “This result is similarly a direct consequence of the principle of effective stress. Nevertheless it is clear that the physical meaning of these tests was in no way understood, and it required the genius of Terzaghi to clarify and enunciate this basic law of the mechanical properties of porous materials”.

2.5.2. Successive Studies and Formulations

After the work of Terzaghi, the concept of ES, although widely used in professional practice, as well as in the scientific field, has been debated for a long time. In 1960, Skempton [5] opened its paper with the following sentence: “… it is of philosophical interest to examine the fundamental principles of effective stress, since it would seem improbable that an expression of the form (σ′ = σ − u), is strictly true”. Moreover, in relatively more recent works, some perplexities have been outlined. For example, in 1990, Oka said [54], “The question of why the effective stress is meaningful, however, still remains unanswered…”, and in 1996, De Buhan and Dormieux [55] referred to ESP as “Terzaghi’s postulate” and tried to clarify several theoretical aspects of ESP validity. In the same year, according to Bluhm and De Boer [56], “Although the phenomenon of effective stresses was known for a long time, the theoretical foundation has remained unsatisfactory until now”.

The main theoretical aspects underlying the ESP are now assessed, and in most cases, the appropriate ES formulation under given hypotheses has been clarified. In 1941, Biot, in developing the theory of poroelasticity, explained, according to a rigorous theory, various aspects about the role of fluids in porous media strain and of ESP [7]. Later, Biot and Willis [57] provided the formulation of the ES according to Equation (6), under the hypothesis of linear elasticity of the porous medium.

In 1960, Skempton carried out an extensive review of the available formulations and experimental data in the literature in order to make the theory clearer and to assess what would be the appropriate ES formulation to be used [5]. Besides providing a mechanical model for the different formulations known at that time, he proved, on the basis of experiments on lead shots carried out by means of a special oedometer, that the formulation by Equation (5) was wrong. At that time, a common opinion was that ES was actually the intergranular stress (Equation (5)) and, as the contact area ratio (a_c) was negligible compared to unity, Equation (5) was reduced to Equation (1). Nevertheless, the experiment involving lead shots under high pressure pointed out that, when the area ratio a_c approached unity, Equation (5) provided ES values that were in significant disagreement with the experimental data, whereas Terzaghi’s Equation (1) yielded an acceptable approximation, and Equation (6) (Biot–Skempton–Nur–Byerlee) provided the best stress–strain data fit. In the same work, Skempton derived the ES formulation by Equation (8) for the strength behavior of soil and rock. Nevertheless, the experimental data proved that Equations (1) and (8) provided indistinguishable ES values. Skempton attributed this result to the circumstance that reduced the values of a_c in soils, and tan ψ in lead shots made the term (a_c tan ψ/tan ϕ′) negligible compared to unity, whereas the a_c value remained uncertain for rocks. In all cases, the simple Terzaghi expression gave ES values that were in excellent agreement with the experimental data.

In 1971, Nur and Byerlee [6] theoretically derived Equation (6) for the stress–strain behavior of rock, highlighting that the main limitation of their proof was related to the assumption of rock linear elastic behavior. Later, Garg and Nur [25] provided a theoretical proof of Equation (6) based on the mixture theory, corroborating the conclusions achieved by Nur and Byerlee and extending their validity to the case of nonlinear elastic behavior of rock. Namely, they established that Equation (6) correctly describes the stress–strain behavior of rock, while Terzaghi’s ES (Equation (1)) adequately describes their strength behavior. Walsh [26,27,28] explained that the nonlinear elastic behavior of rocks is mainly controlled by the change in geometry of more elongated pores and voids. This concept was subsequently further clarified by Walsh and Grosenbaugh [29], as well as later by Zimmerman [19,21,22,23,30] and more recently by several other authors (e.g., Refs. [24,32]). In the late 1970s, Auriault and Sanchez-Palencia [58] provided a theoretical treatise based on a homogenization approach about saturated porous media that also proved the ESP; nevertheless, in their proof, the authors denoted the ES as intergranular stress (“contrainteintergranulaire” in the original paper language), thus causing some confusion about its meaning. Zimmerman [22], similar to Geertsma [59], tackled the poroelasticity problem in terms of compressibility-related variables, such as the volume of porous body and of pores, confining and pore pressures, providing relevant conclusions about stress–strain behavior of nonlinear elastic saturated porous media and on ESP, rigorously corroborating the validity of Equation (6) for stress–strain ES, in which the rock compressibility depends only on the Terzaghi ES (Equation (1)).

Since the late 1980s the study of porous media aroused growing interest, besides in the geological or geotechnical engineering disciplines, also in the field of porous materials of various nature such as metal powders, foams, composite materials, and biological tissues (see Refs. [60,61,62,63]). Consequently, the theory of Poroelasticity has undergone significant developments, also making use of the results of the homogenization criteria [55,58,64,65,66,67,68,69,70,71,72,73,74] and of the modern Theory of Porous Media [25,56,60,61,62,63,75,76,77,78,79,80,81,82,83,84,85,86], unravelling different aspects of porous media behavior such as constitutive equations and micromechanics [66,67,68,69,70,71,87], poroelasticity at finite strain [88,89,90], variational energy formulations [60,79,84,85,86], poroelastodynamics [64,65,72,73,82,91,92], poroviscoelasticity, porothermoelasticity, and porochemoelasticity [20,74,93,94,95].

2.5.3. Currently Open Issues about Fractured Rock and Fracturing Processes

In a recent work, Chen et al. (2020) [96] opened their paper with the following sentence: “The applicability of the classic Terzaghi or Biot effective stress laws to fractured rocks is not clear”, pointing out that some aspects of ESP, after about one century since its introduction, are still unclear. The mentioned work studies the quantification of appropriate effective modules for fractured rocks.

Although the study of the mechanical behavior of fractured rocks goes beyond the issues related to the definition of ES and might fall outside the scope of this work, in this section, we provide a brief overview of some of the main open issues regarding the role of ES in fracturing processes, redirecting the reader to the review paper by Guerriero and Mazzoli [9] for a specific and more exhaustive discussion. Hydraulic properties of fractured rocks are discussed in the sections below. The study of the mechanical behavior of fractured rocks and fracturing processes is problematic due to the notable local heterogeneity of the material, as well as of the pore pressure fields, at different observation scales [97]. Shear or tensional fractures in rock occur over several scales of observation, ranging from large faults to microfractures [98,99,100,101,102,103,104,105,106], and their occurrence significantly affects the hydraulic and mechanical properties of rock. On the other hand, the study of fractured porous media is of considerable interest both for what concerns the hydraulic behavior, e.g., in the field of oil/gas recovery or environmental issues relating to the circulation of ground fluids, and the fracturing processes, e.g., in seismology or structural geology.

The role of pore pressure in fracturing geological processes at several scales of observation, such as faulting and jointing, is still debated. Several works attribute a significant role to the pore pressure in triggering phenomena of earthquakes and co-seismic slip evolution [107,108,109,110,111], which, together with other weakening mechanisms, such as flash heating and frictional weakening (see Refs. [112,113,114]), is among the competing phenomena in fault slip evolution [111,115,116,117,118,119,120,121,122,123,124]. Nevertheless, it is not clear what its weight is and how it varies at different depths in the Earth’s crust [9]. During an earthquake, high temperatures develop on the fault plane, which cause an increase in the pore pressure associated with vaporization (see Refs. [107,108,109,110,111,117]). Such an increase in the co-seismic phase may considerably affect the evolution and speed of sliding; moreover, in the post-seismic phase, it can control the aftershock phenomenon [118], as the increase in pore pressure slowly propagates in the surrounding rock. At present, it is not clear whether this phenomenon is related to the statistical distribution of hypocentral depths observed by Marone and Scholz [119]. Even with reference to artificial fluid overpressure (e.g., by fluid injection in Earth crust; see Refs. [120,125,126] and references therein), the effect of fluid overpressure in enhancing faulting process instability (e.g., according to the Scholz model [127]) still appears controversial and still needs confirmation [128].

Within the framework of tensional joint network development, the role of pore pressure has been the focus of some debate and mechanics involving it remain still unclear to date. It is well-known that tensional joints commonly form under conditions of total compressive stress (with rare exceptions) and in which the pore pressure sufficiently exceeds the lower total stress [34,35,129]. However, fracture-growth models have been subject of debate; moreover, regarding the processes of joint network evolution, some kinds of their spatial distribution, e.g., exponential statistical distribution of spacings in certain categories of joint, remain unexplained to date.

Regarding fracture-growth processes, the first main models were introduced in the 1960s by Hubbert and Willis [129] and later by Secor [35]; nevertheless, Fyfe et al. questioned Secor’s model [33,130]. In the 1990s, Engelder and Lacazette illustrated, in detail, a model of hydraulic fracturing under conditions of constrained lateral expansion (uniaxial strain) [131]. In the same period, other important hydraulic fracturing models have been developed [132,133,134]. Note that the models proposed in Refs. [33,130,131] correctly illustrate hydraulic fracturing processes under uniaxial strain conditions; however, in nature, there are other recurrent situations in which the strain shows different features, such as heterogeneous pore-pressure increase (e.g., fluid injection) and tectonic-related total horizontal stress release (see Refs. [9] and [135]). In such cases, it is appropriate to apply the models described in Refs. [132,133,134].

Regarding the spatial distribution of joints, several statistical models have been proposed for different kinds of fracture [99,100,101,102,103,104,105,136,137,138,139,140]. In layered rocks, there are mainly two categories of tensional fracture: (i) stratabound joint, which entirely cuts the mechanical layer and shows regular spacing, which is proportional to bed thickness, regardless of mineralogy [136,137,141,142]; and (ii) non-stratabound joints, which have irregular spacing (often with exponential distribution) and opening values varying over several orders of magnitude [97,99,100,101,102,103,104,105]. Regarding the formation processes of fracture networks that can give rise to the observed spacings (joint filling processes), several models have been proposed [38,143,144,145,146,147,148]; however, they mainly explain the spacing observed for stratabound fractures but not for non-stratabound ones. The latter often exhibit a uniform random spatial distribution (associated with exponential distribution of the spacing values), which leads to the (possibly apparent) conclusion that, during the joint-filling process, new fracture nucleation appears to occur regardless of the occurrence of other existing joints in the surrounding area. Some models have been proposed in order to try to explain this kind of spatial distribution [37,149,150,151,152,153,154]; nevertheless, at present, it remains unexplained.

3. Fluid Flow, Consolidation, and Fluid–Rock Interaction

3.1. The Early Approaches by Terzaghi and Fillunger

According to Burland [155], Terzaghi’s work has provided the basis to soil and rock mechanics and, more generally, to geotechnical engineering, whereas De Boer and Ehlers [75,77,156] outlined as the studies by Fillunger posed the basis for a rigorous mathematical model of soil/rock consolidation, defining him as the precursor of the Theory of Porous Media.

In fact, Terzaghi had an engineering approach to the problem of subsidence and soil consolidation and provided simplified models that are still widely used in engineering practice today, whereas, on the other hand, Fillunger had a rigorous approach to the above problems and provided rigorous mathematical models that paid particular attention to the methods of averaging of the involved variables. Fillunger’s models were impractical and unsuitable for use by professionals; however, they provided the basis for advanced theoretical studies of particularly complex problems (Section 3.4). Due to the different approach to the problem of consolidation by the two scientists, a bitter scientific dispute arose between them, and this unfortunately led to a tragic ending [75,157].

Terzaghi’s consolidation model considered a one-dimensional condition in which both the deformation of the sample (or soil layer) and the fluidflow occurred only vertically.

Assuming that Darcy’s law is valid, we can obtain the following:

v = − K ∇ h

(19)

where v denotes the fluid velocity vector, h the hydraulic head, ∇ is the well-known differential vector operator Nabla, and K denotes the permeability coefficient.

Furthermore, Terzaghi adopted the following assumptions: constant permeability and linear elastic soil (i.e., constant oedometric modulus); incompressible solid particles and fluid (but compressible porous aggregate).

Denoting the incoming fluid flow per unit volume, the divergence theorem took the following form:

Φ = −∇∙v= K ∇² h

(20)

The ingoing/outgoing fluid rate, a cause of the above hypotheses, was proportional to the porosity time derivative, and then to effective stress derivative, which was, in turn, proportional to the piezometric head derivative:

Φ = constant∙∂h/∂t

Therefore, a well-known parabolic differential equation was achieved:

c_v∇² h = ∂h/∂t

(21)

where c_v denoted a constant depending on permeability, compressibility (oedometric modulus), and density of the porous medium. The solution of this equation was well-known, and this made the model of practical use both in the interpretation of oedometric tests and in the calculation of foundation settlements.

Fillunger, on the basis of the same Terzaghi’s assumptions about fluid and porous medium conditions, followed a different approach. He considered the local Newton’s law for both liquid and solid phases:

ρ_w dv_w/dt = −(Z + ∇P_w)
ρ_s dv_s/dt = (Z − ∇P_s)

(22)

where v_w and v_s denoted the fluid and solid velocities, respectively; Z denotes the interaction force between the porewater and the soil body; and P_w and P_s denote the partial pressure acting on the fluid and solid, respectively. According to the Fillunger model, the average pressure, P, acting over the two-phase medium was individually distributed on the fluid and solid phases as follows:

P_w = n P, P_s = n_s P

(23)

where, clearly, P_w + P_s = P.

It is recalled here that dv/dt = ∂v/∂t + (v ∙ ∇) v.

The continuity equation, both for solid and liquid, was in the following form:

∂ρ_w/∂t = −∇∙(ρ_w v_w)
∂ρ_s/∂t = −∇∙(ρ_s v_s)

(24)

In order to deal with these differential equations, Fillunger, in accordance with the theory of mixtures, substituted stress acting on the solid and pore pressure by the respective partial pressures. In the same way, density was also weighted by a volume fraction; that is, the reduced fluid density, here denoted by ρ_f, consists of an averaged value over the control volume, equal to ρ_w = n ρ_fluid, where ρ_fluid denotes the real fluid density, and n denotes the porosity.The solid reduced density was ρ_s = n_sρ_solid, where n_s was the solid fraction, i.e., (1 −n), and ρ_solid the real solid density. Therefore, solid and fluid phases were substituted by “smeared” continua that could be treated by methods of continuum mechanics. By substitution of the averaged variables in the motion and continuity equations (in one-dimension form), the following differential system was achieved:

∂v_w/∂t + v_w ∂v_w/∂z = −1/(n ρ_w) (Z + ∂(n P)/∂z)
∂v_s/∂t + v_s∂v_s/∂z = 1/(n_s ρ_s) (Z − ∂(n_s P)/∂z)
∂n/∂t = −∂(n v_w)/∂z
∂n/∂t = ∂(n_s v_s)/∂z

(25)

The model developed by Fillunger was very abstract and involved variables that were difficult to detect experimentally, and, therefore, it was not applicable to the study of real cases by engineers and/or designers. However, these models provided the basis for subsequent theoretical developments in the field of theory of porous media and poroelasticity.

3.2. BiotApproach

After Fillunger’s suicide, his theoretical results were forgotten for decades, whereas the methods proposed by Terzaghi found widespread diffusion among scientists and professionals. Only a few studies in the 1950s [158,159] resumed the methods proposed by Fillunger, and from the 1970s, the mixture theory was reused [75,160]. Biot, at the beginning of his career, followed Terzaghi’s approach, giving the basis to the modern Poroelasticity. Biot fully developed the three-dimensional soil consolidation theory [7,161], extending the one-dimensional model previously proposed by Terzaghi [2,13] to more general hypotheses and introducing the set of basic equations of poroelasticity. In this work, Biot explained various aspects about the role of fluids in porous media strain and of ESP. Afterward, Biot and Willis [57] provided the formulation of the ES according to Equation (6), under the hypothesis of linear elasticity of the porous medium. In the same year, Geertsma [59] rederived Biot’s poroelasticity equations, using a different set of variables, and he also proved Equation (6) for the ES, according to the hypothesis of linear elastic behavior. Biot then extended the consolidation theory to the case of anisotropic behavior [162] and later to the case of nonlinear elasticity [163], using second-order approximations of the elastic stress–strain relations. Furthermore, he has discussed in detail the general solutions of the poroelastic and consolidation model [164]. This model has also been adapted in order to study the behavior of saturated porous media in dynamic and seismic conditions [165,166,167,168]. The porous media strain was then analyzed from the thermodynamic point of view under very general hypotheses [169]. After Biot’s extensive work, since the 1980s, the theory of poroelasticity has undergone significant developments thanks to the extensive theoretical work carried out by several authors, such as Zimmerman and coauthors [19,21,22,23,24,30,31,170], as well as Coussy, Ehlers, and De Boer and coauthors, who have developed most of the modern theory of porous media [56,60,61,62,63,76,77,78,80,83,92,160,171].

3.3. Homogenization Approach

The homogenization theory has allowed for a deeper analysis of several problems concerning the porous media behavior. A theoretical analysis of the ESP by means of this approach was carried out by Auriault and Sanchez-Palencia [58] and, under hypothesis of inelastic behavior, by De Buhan and Dormieux [55]. Such an approach was used in order to build up micro-mechanical models for porous media behavior [66,67,68,69], as well as for theoretical approaches to Darcy law [69,70,71], poroelastodynamic models [64,65,72,73], and porochemoelasticity(see Ref. [74]). The homogenization approach assumes periodic pore geometry and then extends the achieved results about the behavior of a template elementary porous volume (representative volume) to an aggregate constituted by identical elementary volumes. Within this latter volume, solid and fluid phases occupy different spaces, i.e., variables such as density, stress, and pressure fields, etc., are defined within the pore or solid space for fluid and solid, respectively. In such spaces, the governing partial differential equations are defined. These are elasticity or other rheological models for solid and Navier–Stokes for fluid. Such equations, together with boundary conditions, are therefore tackled sometimes by means of simplification or linearization or considering second order terms. Such an approach makes extensive use of mathematics and different mathematical treatments may lead to different (and sometimes even inconsistent) results. Therefore, these results always need to be experimentally validated. Nevertheless, compared to classical approaches of poroelasticity, such as those of Terzaghi or Biot, the homogenization methods have the advantage of allowing us to obtain averaged quantities at the scale of the representative volume on the basis of rigorous hypotheses at the microscale (i.e., small compared to pores and/or grains).

3.4. Theory of Porous Media

After its involvement by Fillunger, the mixture theory had been poorly adopted in porous media studies until the 1970s [75]. At that time, Morland [79] used the mixture theory combined with volume fraction in order to build up the Theory of Interacting Continua, within the framework of porous media analysis. Based on such an approach, Garg and Nur [25] have critically re-examined the theoretical and experimental bases of several ES laws. The modern Theory of Porous Media is based on the mixture theory combined with the volume fraction concept [77,78]. Unlike the homogenization approach, in such models, solid and liquid phases are considered to be overlying continua within a porous control volume, whose associated characters, e.g., density, pressure, etc., are weighted by a volume fraction; as an example, the averaged fluid density, here denoted by ρ_w, consists of an averaged value over the control volume, equal to ρ_w = n ρ_fluid, where ρ_fluid denotes the real fluid density, and n denotes the porosity. Analogously, the solid reduced density is ρ_s = n_sρ_solid, where n_s is the solid fraction, i.e., (1 − n) and ρ_solid the real solid density. Therefore, solid and fluid phases are substituted by smeared continua that can be treated by methods of continuum mechanics [60,76,78]. The general form of an averaged variable (e.g., density, pressure, etc.) is according to the following example of pore pressure, where p_w denotes the volume averaged value of the pressure, p:

$${\mathrm{p}}_{\mathrm{w}}=\frac{{{\displaystyle \int }}_{\mathrm{V}}^{}\mathrm{p}\left(\mathrm{x}\right)\cdot \mathrm{g}\left(\mathrm{x}\right) \mathrm{dV}}{{{\displaystyle \int }}_{\mathrm{V}}^{}\mathrm{g}\left(\mathrm{x}\right) \mathrm{dV}}$$

(26)

where V denotes the representative volume, and g(x) determined as follows:

g(x) = 1, inside pore space; g(x) = 0, outside pore space

(27)

Moreover, this approach is largely mathematical and due to its abstractness, its results need to be carefully dealt with in order to avoid some misinterpretations. In the last decades this theory has allowed the analysis of porous media behavior under general hypotheses in several areas of science and engineering, going from geotechnics to biomechanics (see Refs. [60,61,62,63]). The theory of porous media has provided fundamental tools aimed at unravelling several aspects of macro- and micro-scale stress and strain field features, within saturated and unsaturated porous media [25,56,60,80,81,82,83], furthermore has allowed for the development of variational energy formulations [60,79,84,85,86], elastodynamic models for porous media [82,91,92], porochemoelasticity [93,94,95], mechanics of microporous media (see Ref. [87]), as well as dual- and multiple-porosity model aimed to study the behavior of fractured porous deformable media.

3.5. Hydraulic Behavior of Porous Fractured Media

3.5.1. Single-Porosity and Dual-Porosity Models and Their Main Limits

The equations of the theory of fluid flow in permeable media are set by imposing continuity, i.e., that the flow rate entering or leaving an elementary volume is equal to the variation in mass of the liquid inside it per time unit and on some assumptions about hydraulics (e.g., Darcy law) and the elastic behavior of fluid and permeable medium.

It is assumed that flow follows Darcy’s law in its generalized form (Equation (19)), where K is generally a second-order tensor, when permeability is anisotropic (which is a rather common condition in soil and rock), and variable within the studied space. The assumption of Darcian flow is not mandatory and can be replaced, if necessary, by a different formulation, e.g., a nonlinear law. Therefore, the ingoing flow per unit volume, as given by Equation (20), now becomes the following:

Φ = −∇∙v = ∇∙(K ∇h)

(28)

Note that, in several contexts (e.g., reservoir simulation), K is a second-order tensor whose components are generally variable in space (anisotropic heterogeneous medium), and, therefore, in this expression, they too may be subject to derivation. Furthermore, the operator Nabla is multiplied by the vector (K ∇h); for this reason, Equation (28) uses parentheses, although these are unnecessary.

Regarding the elastic behavior of fluid and the permeable medium (e.g., fractured rock), the most common assumption is that of linear elastic behavior. However, even this assumption may, if necessary, be replaced by a more appropriate one. In this regard, it is highlighted here that the compressibility of the rock aggregate, including fractures and pores, must be considered, which may differ significantly from the behavior of the constituting material.

The continuity equation assumes the following form:

Φ = ∇∙(K ∇h) = C_v ∂h/∂t

(29)

where C_v denotes a constant parameter which depends on the elastic moduli of fluid and fractured rock and on fluid density. In many practical problems, e.g., in the field of oil recovery or other studies involving ground fluids, it must be taken into account that fluids with different viscosities and specific densities (oil, gas, water, etc.) move within the fractured rock. The model described here is called singleporosity, as opposed to the more recent dualporosity described below, as it assimilates the fractured porous medium to a locally homogeneous permeable medium (i.e., in which each volume element is homogeneous). Despite the introduction of the more realistic dual-porosity model, which has taken place for several decades, the single-porosity model is currently still widely used.

The dual-porosity model was initially formulated by the mathematicians G. I. Barenblatt et al. [172] and subsequently developed for oil applications by Warren and Root [173]. It takes into account that, within rock, the fluid flow exhibits different features within the porous matrix (low permeability) and fractures (high permeability). This model considers the fractured porous rock as consisting of two overlying sub-systems: (i) the porous matrix, characterized by high porosity and low permeability (primary porosity); and (ii) the fracture network, characterized by low porosity and high permeability (secondary porosity). In such permeable subsystems, two overlying pore pressure fields, p_m and p_f, respectively, are defined. For each representative volume, the two pressure values relating to pore and fracture space can be averaged according to the criteria of the mixture theory, i.e., through Equation (26). Following the same line of reasoning as in the previous section, for each of the considered systems, the continuity equations can be written by taking into account their reciprocal interaction. In this case, in such a model, due to the low permeability of the non-fractured matrix, it is assumed that the exchange of fluids between adjacent representative volumes of rock occurs only through the fracture system [173]. Therefore, in each representative volume of rock, the outgoing flow from the porous matrix, which we indicate here as Φ_mf, is fully conveyed into the fracture system. Therefore, the first balance equation is written as follows:

Φ_mf = C_m n_m ∂p_m/∂t

(30)

where C_m denotes a constant which depends on the elastic moduli of fluid and of the non-fractured porous rock, on fluid density. The outgoing flow from a generic representative volume (exclusively through the fracture system), therefore, is given by the fluid released (or trapped) by the fractures due to a reduction (increase) of pore pressure, to which is added the flow coming from the porous matrix:

Φ_f = C_f n_f ∂p_f/∂t + Φ_mf

(31)

where C_f has the same meaning as the corresponding parameter in Equation (30), but referring to the fractured rock. Equations (30) and (31) are coupled by assuming that the fluid exchange between the matrix and fracture network (i.e., Φ_mf) is proportional to the difference in mean pressure between the two subsystems (steady-state hypothesis) [173]:

Φ_mf = λ∙(p_f − p_m)

(32)

where λ denotes a coefficient, whether to be experimentally determined or to be estimated theoretically, depending on the fluid characters, matrix permeability, and geometry of pores and fracture network. The condition expressed by Equation (32) can be replaced in cases where it is deemed unrealistic (unsteady state model; see Refs. [174,175,176]).

As in the case of the single-porosity model, also for the dual-porosity model the hypothesis that different fluids can circulate in the analyzed medium can be considered. Therefore, the model can be modified taking into account the characteristics of the considered fluids (dual- or multi-permeability model). The Equations (30)–(32) and the equivalent to Equation (29), related to the two pressure fields, provide a system of parabolic differential equations which, with opportune boundary conditions, can be solved by numerical methods. In numerous simple cases analytical solutions in closed form are available.

3.5.2. Multiple Porosity Models and Still Open Issues

One of the main limitations of the single-porosity model is related to the fact that it replaces each representative volume of fractured rock with an equivalent homogeneous (i.e., non-fractured) representative volume, showing the same hydraulic properties. Nevertheless, in Ref. [97], it was explained that representative volumes of fractured rock showing identical behavior under steady-state flow conditions (i.e., when ∂p/∂t = 0), but with different structural characteristics, can exhibit significantly different behavior and fluid-discharge time in an unsteady condition. The dual-porosity model should make it possible to overcome this obstacle. However, the numerical models developed in practice include only major fractures in the secondary porosity, whereas the matrix includes pores and all minor fractures in the primary porosity. Therefore, in such cases, the dual-porosity model may reproduce the same errors related to the single-porosity on a smaller scale.

In nature, rocks show fractures distributed over several scales of observation. These form a complex fracture network composed of sub-systems, with each one showing different geometry and features. In carbonate bedded rocks, usually four hierarchical permeable systems characterized by different size and permeability structure may be distinguished in building up a rock permeability model [97]: (i) the large-scale first-order structures consist of faults, with which high-permeability structures (damage zones) are associated; (ii) at the meter scale, stratabound joint systems, together with bedding joints bounding the mechanical layers, form a well-connected permeable network transporting fluids to the fault system; (iii) at a lower scale, down to crystal size, non-stratabound joints constitute a pervasive and capillary fracture network which conveys fluids within the rock mass; and (iv) the non-fractured host rock constitutes the lower permeability system (at the pore scale). Due to significant differences in permeability, such a system can be considered hierarchical; that is, each permeable system is assumed to convey fluids exclusively to that of immediately higher order. The numerical modeling of such a hierarchical system may be carried out by means of multiple-porosity models.

The concept of a system formed by superimposed permeable sub-elements, so as provided by the dual-porosity model, can be extended to the case in which more than two permeable structures within the fractured rock are considered, such as non-fractured matrix and different fracture networks characterized by different geometry and permeability features. In this case, it is possible to define multiple pressure fields, p_j, within the generic representative volume, associated with the different considered subsystems. Al-Ahmadi and Wattenbarger [177] illustrated different formulations of triple-porosity models, based on different working hypotheses. The equations for a hierarchical model are illustrated below, where p₁ is the pressure defined within the matrix; and p₂, p₃, … p_j indicate the pressure values relating to N − 1 different subsystems considered (e.g., various fracture networks with different characteristics), using the same notations of the previous section:

Φ₁₂ = C₁ n₁ ∂p₁/∂t
Φ₂₃ = C₂ n₂ ∂p₂/∂t + Φ₁₂
Φ_j,j+1 = C_jn_j ∂p_j/∂t + Φ_j−1,j
Φ_N = C_Nn_N ∂p_N/∂t + Φ_N−1,N

(33)

Following the same reasoning illustrated above for the dual-porosity model and using a similar notation, the following relations are obtained:

Φ_j,j+1 = λ_j,j+1∙ (p_j − p_j+1)
Φ_N = ∇∙K_N ∇h_N

(34)

where h_N denotes the hydraulic head associated with p_N.

The above-described equations govern the behavior of a hierarchical multiple-porosity system. These equations can be suitably modified in order to describe non-hierarchical systems (in which some sub-systems may act in parallel), as illustrated by Al-Ahmadi and Wattenbarger [177], or when more different fluids are considered within the fractured medium.

The first multiple-porosity model (a triple-porosity one) was developed in the early 1980s [178,179]. A few years later, such a model was involved in the petroleum literature [180,181]. More recent works employed multiple-porosity models in various fields, such as oil-reservoir simulation and nuclear-waste management [177,182,183,184,185,186,187,188,189,190,191]. Currently, multiple porosity models are still little used and are not available in most simulation software on the market.

4. Conclusions

Over the past century, the underlying mechanicsof the ESP and fluid–rock interaction in deformable porous fractured permeable media has been extensively studied and several aspects have been clarified; however, some problems are still open. Regarding the ESP in saturated media, an overview of the literature points out that, in the case of strength behavior of a great variety of porous media—rock, cement and concrete, soil, metal powders, etc.—Terzaghi’s expression (σ′ = σ − p) provides an excellent approximation for ES, and this is the correct formulation to use in these cases.

When studying the stress–strain behavior of soil and rock under high pressure, the Skempton–Nur–Byerlee formulation (σ′ = σ − (1 − K/K_s) p) allows us to take into account the often modest cubic volume strain associated with pore pressure, which it is superimposed to the (non-cubic) volume strain associated with the other acting stresses (i.e., those exceeding the pore pressure). In all cases, including those of porous media under high pressures, the geometry of the voids (pores and fractures in concrete and rocks; pores and grain spatial arrangement in granular media) is controlled, with excellent approximation, exclusively by the Terzaghi’s ES. It is important to keep this in mind because, with the same constituent material, the geometry of voids controls all the mechanical and hydraulic properties of porous (and/or fractured) media. To obtain more precise formulations, e.g., in the study of special situations or particular materials in which the classical formulation of the ES could significantly deviate from the experimental observations, more researched theories may be needed (e.g., Theory of Porous Media, poroelasticity, numerical models, etc.).

Regarding the phenomena of coupled fluid flow and stress–strain involving porous permeable deformable media, the methods originally proposed by Terzaghi (theory of one-dimensional consolidation) are still used (eventually modified) in current practice in civil engineering, whereas, in other fields of study (e.g., subsidence, seismology and fracturing, petroleum recovery, hydrogeology, environmental studies, and many others), more elaborate models are often used that were developed later on. In the scientific community, the current trend is moving toward the abandonment of single-porosity models in favor of dual- and multiple-porosity models and the study of flow characteristics at various observation scales, ranging from the pore scale to that of the entire area under consideration. In these areas, various problems are still open, particularly in the study of coupled fluid flow and stress–strain in fractured porous media.