1 Introduction

For constitutive modeling of unsaturated soils, a somewhat modified definition of the Bishop stress is proposed. In the literature, several alternatives to Bishop stress \(\sigma ^B\) have been proposed, see [3] p. 17 for a list. The novel stress \(\sigma ^E\) is a macroscopic description of the grain-to-grain forces. Similar to the original Bishop stress \(\sigma ^B\), the new \(\sigma ^E\) alone is insufficient [5] for a constitutive description of unsaturated soil, unless an explicit dependence on water suction is introduced. Larger suction, for example, stiffens the unsaturated soil during virgin compression at a given \(\sigma ^B \), and decrease in suction may cause irreversible deformations (collapse), see the Barcelona Basic ModelFootnote 1 (BBM) [1] for a constitutive description of these effects.

There is no convincing reason to treat \(\sigma ^B\) as the ”effective stress” because \(\sigma ^B\) cannot be employed without suction. The advantage of \(\sigma ^E\) is that it has a sound physical and microscopic basis. This is in contrast with \(\sigma ^B\) that was originally used in modeling of strength and deformation alone (without suction). Furthermore, \(\sigma ^E\), unlike \(\sigma ^B\), satisfies the geostatic equilibrium with buoyant weight \(\gamma '\).

Grain-to-grain forces \(F_i^{K}\) from all contact points K on the border S surrounding a representative material element of volume V can be used to determine the effective stress (denoted as \(\sigma \) and not \(\sigma '\) here). The effective stress tensor is represented in the index notation by the average value of dyadic products

$$\begin{aligned} \sigma _{ij} = - \frac{1}{V} \sum _{K} F_i^{(K)} x_j^{(K)}\,, \end{aligned}$$
(1)

where \(x_j^{(K)} \) refers to the position of the contact point K. The positive compression sign convention is employed. The relationship (1) is often employed in the DEM postprocessing. Similarly, the average divergency of stress is

$$\begin{aligned} e \frac{1}{V} \int _V \sigma _{ij,j} \mathrm{d}V = \frac{1}{V} \int _S \sigma _{ij} n_j \mathrm{d}S = \frac{1}{V} \int _S -t_i \mathrm{d}S = -\frac{1}{V} \sum _K F_i^{(K)} \end{aligned}$$
(2)

with the effective traction \(t_i\) and corresponds to the resulting contact force per unit volume. The formulae above are commonly used for saturated soils. Here, they will be adopted to define the effective stress for unsaturated soils. This is preferable to ambiguousFootnote 2 and non-objective definitions based on estimations of how well the stress would perform, if unsaturated soils were calculated as saturated.

2 Scope of application of \(\sigma ^E\)

The degree of saturation S is used to divide the unsaturated soils into three groups:

  • \(0 \le S < S_C \) isolated water regions or lenses (pendular regime)

  • \(S_C< S < S_B\) interconnected water and air (funicular regime)

  • \(S_B < S \le 1 \) isolated air bubbles

with the bounds \(S_C \approx 0.2\) and \(S_B \approx 0.8\) being tentative estimates, see [8] and Sect. 6 for microstructural considerations. The novel stress \(\sigma ^E\) is applicable to soils with all grains entirely surrounded by either free or capillary (and hence interconnected) pore water. It must be able to transfer the hydrostatic pressure so it cannot be just a ”film” water. The proposed stress \(\sigma ^E\) is not suitable for the pendular regime. The proposed stress has the same form

$$\begin{aligned} \sigma ^E = \sigma ^\mathrm{tot}- \chi ^E u \text {with pore pressure} u = u_w - u_\mathrm{atm}\end{aligned}$$
(3)

as the Bishop stress \(\sigma ^B = \sigma ^\mathrm{tot}- \chi (u_w - u_\mathrm{atm}) \), where \(u_w\) and \(u_\mathrm{atm}\) denote the absolute pressures in water and air. The difference lies in factor \(\chi ^E \ne \chi \). The new definition should be used under static conditions and without seepage forces. The suction \(s = -u = u_\mathrm{atm}- u_w\) will be indispensable in modern constitutive models employing \(\sigma ^E\).

3 Divergency of stress in unsaturated soil

Factor \(\chi ^E(S,n)\) for \(\sigma ^E\) will be derived from the geostatic equilibrium in a 1D column of soil, Fig. 2. The divergency of the total stress follows from the static equilibrium

$$\begin{aligned} \partial _x \sigma ^\mathrm{tot}= \gamma \text { with }\gamma = (1-n) \gamma _s + n S \gamma _w\,, \end{aligned}$$
(4)

where x = coordinate axis that points downwards and starts at the ground surface; \(\sigma ^\mathrm{tot}\) = vertical total stress, compression positive; \(\gamma _s \approx 26.5 \) kN/m3 intrinsic unit weight of grains; \(\gamma _w = 10 \) kN/m3 intrinsic unit weight of water; \(n = e/(1+e)\) porosity; S = degree of saturation; \(\gamma ,\gamma _d = (1-n) \gamma _s \) moist and dry unit weights

Because the capillary pore water is interconnected, its pressure u has a vertical gradient

$$\begin{aligned} \partial _x u = \gamma _w \end{aligned}$$
(5)

The interconnected water channels resemble a system of water-filled pipes. Pascal’s law holds for such a system: in a fluid at rest, a pressure change in one section of a fluid at rest is transferred without loss to every portion of the fluid. As a result, (5) can be applied also to unsaturated soils with interconnected capillary water channels. The pore pressure \(u(x) = (x - h_c) \gamma _w \) may be simply calculated using the gradient (5) and starting from \(u(h_c)=0\) at the water level \(x=h_c\), where \(h_c\) is the height of the capillary rise. This relationship holds true both below \((x > h_c)\) and above \( (x< h_c)\) the water tableFootnote 3 irrespectively of saturation, see also [7] p. 21. The derivation of \(\sigma ^E\) described in Sect. 5 is based solely on (4) and (5). According to Terzaghi (1936):

If the voids of the soil are filled with water under a stress u, the total principal stresses consist of two parts. One part, u, acts in the water and in the solid in every direction with equal intensity. It is called the neutral stress (or pore water pressure). The balance \(\sigma = \sigma ^\mathrm{tot}- u\) represents an excess over the neutral stress u, and it has its seat exclusively in the solid phase of the soil. This fraction of the total principal stresses will be called the effective principal stresses.

Such stress should represent the contact forces between grains. Only three types of forces act on a solid grain surrounded by water, Fig. 1. Counting per unit volume of soil, they are:

  • The resultant of all contact forces (or \( - \partial _x \sigma ^E\) macroscopically)

  • The unit self-weight \(\gamma _d = (1-n)\gamma _s\) and

  • The unit buoyant force \( b= -(1-n) \gamma _w\) (from \(-\partial _x u (1-n)\))

These forces must be in static equilibrium and hence

$$\begin{aligned} -\partial _x \sigma ^E + \gamma ' = 0 \text { with }\gamma ' = (1-n)(\gamma _s - \gamma _w) \end{aligned}$$
(6)
Fig. 1
figure 1

Unit self-weight \(\gamma _d = (1-n) \gamma _s\) and buoyancy \( b = -(1-n) \gamma _w\) must be balanced by the macroscopic equivalent \(- \partial _x \sigma ^E\) of resultant contact force, \(- \partial _x \sigma ^E + (1-n) (\gamma _s -\gamma _w) = 0\)

For isolated air bubbles entrapped between grains, the buoyant unit weight in (6)\(_2\) may need a minor modification. If \(U_b\) is the fractionFootnote 4 of pore volume occupied by bubbles, then (6)\(_2\) has the form

$$\begin{aligned} \gamma ' = (1-n)(\gamma _s - \gamma _w) - n U_b \gamma _w \,, \end{aligned}$$
(7)

The extra term \(- n U_b \gamma _w\) in (7) results from the buoyancy of air bubbles pushing the soil skeleton upwards.

In this research, the geostatic equilibrium condition (6) is crucial. To the author’s knowledge, \(\sigma ^B\) has previously been neglected in that it violates (6), namely

$$\begin{aligned} \partial _x \sigma ^B = \partial _x ( \sigma ^\mathrm{tot}- S u ) = \gamma - S \gamma _w \ne \gamma ' \end{aligned}$$
(8)

for \(S < 1\) and assumingFootnote 5\(S(x) = \mathrm{const}\).

4 Weight carried by capillary surface tension

Let us start with a layer of saturated capillary water of thickness \(h_c\) located directly below the ground surface \(x=0\) and above the water table \(x = h_c\). As shown in Fig. 2c, the coordinate x is heading downwards. Two mathematically equivalent formulae can be used to calculate the effective vertical stress

$$\begin{aligned} \sigma (x)= & {} \sigma ^\mathrm{tot}- u = \gamma _r x - \gamma _w ( x- h_c ) \text { or }\nonumber \\ \sigma (x)= & {} q + \gamma ' x \text { with }q = \gamma _w h_c \,, \end{aligned}$$
(9)

where \(\gamma _r\) is the unit weight of the saturated soil and \(\gamma '\) is the buoyant unit weight (both are assumed constant for simplicity). The surface load q consists of microscopic forces exerted on skeleton by capillary water ”hanging” on skeleton. They are equally distributed at \(x=0\), as shown in Fig. 2a.

In the literature, if the capillary zone is not saturated, \(S < 1\), (9)\(_1\) is frequently replaced by \(\sigma ^B(x) = \sigma ^\mathrm{tot}- \chi u \). However, one might also try modifying (9)\(_2\). Doing so, the buoyant unit weight \(\gamma '\) remains unchanged, but the expression for q is different. Let the degree of saturation, \(S < 1\), be constant throughout the capillary zone for simplicity.

The surface tension must carry not only the

Fig. 2
figure 2

Weight \(n S \gamma _w \) of the capillary water but also the reaction \(-b\) to the buoyancy that acts on water downwards

At \(x=0\), the following additional load is transferred from the water to the skeleton

$$\begin{aligned} q = n S \gamma _w h_c + (1-n) \gamma _w h_c \end{aligned}$$
(10)

Two contributions to q in (10) (both acting downwards) are:

  • The weight of capillary water \(n S \gamma _w h_c \)

  • The reaction to buoyant force \(-b = (1-n) \gamma _w h_c\) , see Fig. 2b.

It turns out that the derivation of \(\sigma ^E\) from Sect. 5 is in agreement with (10) in the sense that \(\sigma ^E(0) = q\). For saturated soils, (10) and (9) produce the same \(q=\gamma _w h_c\).

5 Derivation of \(\sigma ^E\) from its divergency

Figure 3 shows a soil layer with a capillary zone directly beneath the ground level.

Stress

Fig. 3
figure 3

\(\sigma ^E\) in the capillary zone \(x \in (0,h_c)\)

The difference between divergencies from (4) and (6) is

$$\begin{aligned} \partial _x ( \sigma ^\mathrm{tot}- \sigma ^E ) = \gamma - \gamma ' \end{aligned}$$
(11)

Starting from the water level \(x=h_c\) (where \(u=0\) and \( \sigma ^\mathrm{tot}= \sigma ^E \)), this difference may be conveniently integrated with respect to x. If the upper limit \(\bar{x}\) of integration lies in the capillary zone \( 0 \le \bar{x} \le h_c\) then

$$\begin{aligned} \left. (\sigma ^\mathrm{tot}- \sigma ^E ) \right| _{\bar{x}}= & {} \int _{h_c}^{\bar{x}} (\gamma - \gamma ') \mathrm{d}x \end{aligned}$$
(12)

Application of (12) for the case of isolated air bubbles first leads to the formula

$$\begin{aligned} \left. ( \sigma ^\mathrm{tot}- \sigma ^E ) \right| _{\bar{x}}= & {} \int _{h_c}^{\bar{x}} \gamma _w \mathrm{d}x = u(\bar{x}) \text { or }\sigma ^E = \sigma ^\mathrm{tot}-u \end{aligned}$$
(13)

which is simply obtained inserting \(\gamma \) from (4) and \(\gamma '\) from (7) with \(U_b = 1-S\) into (12). The integral \(\int _{h_c}^{\bar{x}} \gamma _w \mathrm{d}x = \gamma _w(\bar{x} - h_c )\) has been expressed by the pore water pressure \(u(\bar{x}) < 0\) in (13), to demonstrate that the novel stress \(\sigma ^E\) is identical with the Terzaghi’s effective stress formula. It is so, if pore water contains isolated air bubbles only.

Next, the funicular regime with interconnected air and water channels (no bubbles or water lenses) is examined. If \(\gamma \) from (4) and \(\gamma '\) from (6) are substituted into (12), the result is

$$\begin{aligned} \left. ( \sigma ^\mathrm{tot}- \sigma ^E)\right| _{\bar{x}}= & {} \int _{h_c}^{\bar{x}} (1- n + nS ) \gamma _w \mathrm{d}x = (1- n + nS )^\mathrm{av}u(\bar{x})\,, \end{aligned}$$
(14)

where \(()^\mathrm{av}\) denotes the average value over the range \(x\in (\bar{x} , h_c)\). For constant porosity and constant degree of saturation, (14) can be simplified to

$$\begin{aligned} \sigma ^E = \sigma ^\mathrm{tot}- (1-n + nS ) u \text { for }~~ n(x)= \mathrm{const}, ~~S(x) = \mathrm{const} \end{aligned}$$
(15)

Comparing \(\sigma ^E\) from (15) with the Bishop stress [2] \( \sigma ^B = \sigma ^\mathrm{tot}- \chi u \), one obtains

$$\begin{aligned} \chi = \chi ^E = (1-n + nS ) \end{aligned}$$
(16)

The empirical estimates of the Bishop parameter \(\chi \) are usually different from \(\chi ^E\), and hence, they violate (6), i.e., \(\partial _x \sigma ^B \ne \gamma ' \).

Definition (15) is not valid, unless all grains are entirely surrounded by capillary water. The averaging of the factor \( \chi ^E = (1-n + nS )\) is required in the case of spatial variability \(n(x) \ne \mathrm{const}\) or \(S(x) \ne \mathrm{const}\) , see (14). A retention curve S(u) could be useful in this case.

6 Micro-mechanical considerations

A single water lens of diameter \(d_L\) between two spherical grains is considered, Fig. 4. The water lens glues a pair of grains by the force

$$\begin{aligned} \frac{1}{4} ( u_\mathrm{atm}- u_w ) \pi d_L^2 + T \pi d_L\,, \end{aligned}$$
(17)

where \(T \approx 73 \cdot 10^{-6}\) kN/m is the surface tension or interfacial (air–water) energy per unit area. The second summand of the aforesaid attractive force is an estimation based on an infinitesimally small separation of grains \(\delta x\) which is supposed to increase the interfacial surface area by \( \pi d_L \delta x \) and the interfacial energy by \( \delta E = T \pi d_L \delta x \). As a result, \( \delta E/\delta x = T \pi d_L\) is the force per water lens.

Under the funicular regime, Fig. 5, water is enclosed by a single but geometrically complicated air–water membrane (interface surface). The area of the membrane may somewhat grow due to a small separation \(\delta x\) of soil portions on opposite sides of a planar cross section. This increase in area is expected to be much smaller than that under pendular regime.

Fig. 4
figure 4

A grain under the pendular regime. Neighbouring grains are glued by the suction \(s = u_\mathrm{atm}- u_w\) and by the surface tension T

For example, let a cross section of a monodisperse silt with uniform grain diameter \(d = 0.03\) mm have a simple cubic packing with \(d^{-2}\) grain-to-grain contacts per unit area. Let each contact be glued by an isolated water lens of diameter \(d_L = 0.02\) mm. For \( d^{-2}\) lenses with the perimeter \(\pi d_L \) each, the contribution of the surface tension T to the compressive stress is

$$\begin{aligned}&d^{-2}\cdot T \cdot \, \pi \cdot d_L\nonumber \\&\quad = ( 0.03 \cdot 10^{-3} )^{-2} \, \cdot 73 \cdot 10^{-6}\cdot 3.14 \cdot \, 0.02 \cdot 10^{-3} \approx 5 \, \text {kPa} \end{aligned}$$
(18)

Under the funicular regime, the increase in the surface of the air–water membrane due to the separation \(\delta x\) is much smaller than that of isolated water lenses. The effect of surface tension on the effective stress is practically insignificant.

Fig. 5
figure 5

A cross section of unsaturated soil in funicular regime. Neighbouring grains aggregates are glued by the surface tension T on the air–water membrane

Nonetheless, under the funicular regime, unsaturated soils respond substantially differently from saturated soils under comparable conditions. The presence of the air–water membrane appears to be significant for the force chain stability and must be considered in constitutive models (but not as a part of the effective stress). Because all neighboring grains (or grain aggregates) are glued together, the pressure (grain-to-grain contact forces) due to surface tension is qualitatively different than the the same (in average) pressure applied directly to the soil skeleton at the boundary of the sample. The latter pressure is transferred through fewer and longer force chains that are more prone to buckling. Hence, the presence of the air–water membrane with surface tension has an important supporting effect on long force chains but has little influence on the average grain-to-grain forces and consequently on the effective stress.

7 Generalization to 3D

In Sect. 5, Eq. (15) for \(\chi ^E\) was derived for the vertical direction only. In this section, this equation is demonstrated to be valid also in the general case, such as horizontal direction. Because the horizontal stress component is statically indeterminate, the equilibrium principle cannot be applied. Instead, a micro-mechanical argument, Fig. 6, will be used. A unit plane cross section through a capillary zone is considered with water under tension, \(u = u_w - u_\mathrm{atm}<0\). The forces acting on such cross section can be calculated using the surface fractions of solid and water. The 2D porosity, i.e., the fraction of the cross section passing through pores, is identical as the usual porosity \(n = V_\mathrm{pores}/V\), according to the Delesse’s (1843) stereology principle. This principle is supposed to apply to the degree of saturation S as well, i.e., the fraction of the cross section passing through water is nS. The total traction vector on the cross section results from the following sources:

  • The contact forces acting through grains only

  • Pore pressure acting through the grains (with fraction \(1-n\))

  • Pore pressure acting directly through the pore water (with fraction nS)

  • Surface tension T (negligible for the contact forces but critical for the stability of force chains)

Adding all portions of pore pressure to the contact forces (and neglecting T), one obtains

$$\begin{aligned}t^\mathrm{tot}_i = t^E_i - n_i (1-n) u - n_i n S u \text { or }\sigma ^\mathrm{tot}_{ij}\nonumber\quad = \sigma ^E_{ij} + (1-n + S n) u \delta _{ij}\,, \end{aligned}$$
(19)

where \(n_i\) is the unit outer normal vector to the cross section, \(t_i^E = -\sigma _{ij}^E n_j^E\) is the effective stress vector (traction), and \(\delta _{ij}\) is the Kronecker symbol. The expression (19) is in agreement with the 1D derivation (15).

The absolute atmospheric pressure \(u_\mathrm{atm}\) acts on the cross section:

  • Though water on grains with surface fraction \(1-n\)

  • On water with surface fraction nS

  • Directly with surface fraction \(n (1-S)\)

i.e., on the whole area \( 1-n + n S + n (1-S) = 1\) of the cross section. However, the atmospheric pressure is excluded from the total net stress \(\sigma _{ij}^\mathrm{tot}\) by definition, and hence, it is absent in (19).

Fig. 6
figure 6

Water forces on a vertical unit cross section. Pore water acts directly through the surface nS and indirectly via grain through the surface \((1-n)\)

8 Example of an initial equilibrium

Modern constitutive models for unsaturated soils use suction \(s = -u\) explicitly in combination either with the total net stress \(\sigma _{ij}^\mathrm{tot}\) or with \(\sigma _{ij}^B\). Given \(\sigma _{ij}^E\) and u, one can use an existing constitutive code as follows: (1) convert \(\sigma _{ij}^E \rightarrow \sigma _{ij}^B\) (2) update \(\sigma _{ij}^B\) as usual (3) convert \(\sigma _{ij}^B \rightarrow \sigma _{ij}^E\). It would be more elegant to reformulate the constitutive equations using \(\sigma _{ij}^E\), of course.

As an example, let us consider the initial geostatic equilibrium in a column of soil with the total depth of \(H = 15\) m and with the water level at the depth of \(h_c = 6\) m, see Fig. 3 left. The simplified fields of void ratio \(e(x) =1\) (which corresponds to \(n=0.5\)) and the degree of saturation

$$\begin{aligned} S(x)= {\left\{ \begin{array}{ll} 0.4 \text { for }&{} x < h_c = 6 \\ 1.0 \text { for }&{} x > 6 \end{array}\right. } \end{aligned}$$
(20)

are assumed. The unit weights are \(\gamma _d = 14\), \(\gamma _s = 28\), \( \gamma _w = 10\), \(\gamma ' = 9\), \(\gamma = (1-n) \gamma _s + n S \gamma _w = 16 \) and \(\gamma _r = (1-n) \gamma _s + n \gamma _w = 19 \) kN/m\(^3\). The pore pressure is hydrostatic varying linearly from − 60 kPa on the ground level to 90 kPa at the bottom. The initial stress field is in equilibrium with the self-weight. It will be calculated analytically and with two popular commercial software packages Plaxis and Abaqus.

8.1 Analytical solution \(\sigma ^E \)

According to (15) with \(\sigma ^\mathrm{tot}= 0\) and \(u=-60\) kPa on the ground surface one obtains

$$\begin{aligned} \sigma ^E(0)= & {} -(1 - n + n S ) u = (1 - 0.5 + 0.5 \cdot 0.4 ) \cdot 60 = 42 \text { and }\\ \sigma ^E(x)= & {} \sigma ^E(0) + \gamma ' x = 42 + 9 \cdot x \end{aligned}$$

The stress \(\sigma ^E\) increases linearly with depth x as given in (6). This stress field is valid in the whole column, above and below the water table.

8.2 Defining initial stress for \(\sigma ^B = \sigma ^{\mathrm tot}- S u\) Abaqus

The initial equilibrium is calculated with the option *SOILS using 300 CPE4P elements. According to the manual, the geostatic initial effectiveFootnote 6 stress field has to be defined as

$$ \sigma ^{B} (x) = \left\{ {\begin{array}{*{20}l} {\sigma ^{{{\text{tot}}}} - Su{\text{ with }}\sigma ^{{{\text{tot}}}} = \gamma x = 16 \cdot x} \hfill & {{\text{for }}x < h_{c} } \hfill \\ {\sigma ^{{{\text{tot}}}} - u{\text{ with }}\sigma ^{{{\text{tot}}}} = \gamma h_{c} + \gamma _{r} (x - h_{c} ) = 19 \cdot x - 18} \hfill & {{\text{for }}x > h_{c} } \hfill \\ \end{array} } \right. $$

wherein the pore pressure is \(u(x) = -60 + 10 x\) kPa. A linear interpolation between the following values can be prescribed:

\(\sigma ^B(0) = 60 \cdot 0.4 = 24 \) and

\(\sigma ^B(6) = 16 \cdot 6 = 96 \) and

\(\sigma ^B(15) = (19 \cdot 15 - 18) - (-60 + 10 \cdot 15) = 177 \) kPa

The initial degree of saturation S(x) was defined via the discrete nodal values. Due to the jump from \(S=0.4\) to \(S=1.0\) at \(x=6\) m, the intermediate initial value \(S(6) = 0.7\) was chosen. In order to enforce \(S(x)=0.4\) everywhere in the capillary zone, the following (artificial) retention curve

$$\begin{aligned} S(u) = \left\{ \begin{array}{ll} 0.4, &{} \text { for }u < 0\\ 1, &{} \text { for }u \ge 0 \end{array} \right. \end{aligned}$$
(21)

was defined. Using a linear elastic material with \(E=10\) MPa and \(\nu =0\), the desired zero displacement field was obtained from initial equilibrium using a tiny adjustment \(\sigma ^B(0) \approx 27.2\) kPa (discretization with 30 elements) or \(\sigma ^B(0) \approx 24.15\) (discretization with 300 elements). Above the ground water level, the gradient \(\partial _x \sigma ^B = 12\) is larger than \(\gamma '=9\) kN/m\(^3\).

8.3 Calculation of initial stress \(\sigma = \sigma ^{\mathrm{tot}}- u\) with Plaxis

In Plaxis, the total unit weights \(\gamma = 16\) and \(\gamma _r = 19\) kN/m\(^3\) need to be prescribed above and below the water surface, respectively. The initial stress is calculated internally. Displacements to establish the static equilibrium are all zero as required, but the resulting ”effective” stress varies linearly according to

$$\begin{aligned} \sigma (x) = \sigma ^\mathrm{tot}(x) - u(x) \end{aligned}$$
(22)

in the whole column. In particular, \(\sigma (0) = 60\) kPa seems to be strongly overestimated. The equilibrium condition (6) is violated by Plaxis. Using some sophisticated options [4], the solution given by Abaqus can also be reproduced by Plaxis.

9 Summary

A novel effective stress \(\sigma ^E\) is defined by using \(\chi ^E\) instead of Bishop parameter \(\chi \). Regardless of the degree of saturation, this stress represents the grain-to-grain contact forces. Such physically based definition is beneficial for the constitutive modeling. The new stress, unlike Bishop stress \(\sigma ^B\), is in the geostatic equilibrium with the buoyant weight \(\gamma '\). The ”effective” stresses from Abaqus \(\sigma ^B \approx \sigma ^\mathrm{tot}- S u \) or from Plaxis \( \sigma = \sigma ^\mathrm{tot}- u\) are not in equilibrium with \(\gamma '\).

However, \(\sigma ^E\) is still just a definition and not a comprehensive constitutive model. It may become a part of a constitutive model and contribute to its clarity. The proposed \(\sigma ^E\) as such cannot simulate or explain experiment results. Existing constitutive models for unsaturated models that employ Bishop stress \(\sigma ^B\), and suction s can be adapted or modified to work with \(\sigma ^E\) and s instead. The validity of \(\chi ^E\) requires that all grains must be surrounded by capillary water.