Micromechanics-Based Strength Criterion for Root-Reinforced Soil

1. Conceptual Framework for Constructing the Strength Criterion

Test results show that the shear strength of root-reinforced soils is influenced not only by confining pressure but also by root content and the number of freeze–thaw cycles (Luo et al., 2023). If the rooted soil is treated as a single whole, the interaction between roots and soil cannot be captured. As shown in Figure 1, at the macroscopic scale the bonded element is taken as the matrix and the frictional element as the inclusion; a binary-medium representative volume element (RVE) is thus formed by bonded and frictional elements. At the mesoscopic scale, the bonded element itself consists of a cemented-soil matrix with root inclusions.

Through two homogenization steps, this model yields a strength criterion for rooted soils from the meso- to the macro-scale; a schematic of the homogenization procedure is given in Figure 2.

(i) Homogenization I: Homogenizing the cemented-soil matrix with the root inclusions gives the strength criterion F_b for the bonded element. Both the cemented-soil matrix and the roots are modeled as linear elastic; their strength criteria share the same form, and no plastic dissipation occurs during loading. Therefore, an elastic homogenization approach can be used directly to establish the strength criterion of the bonded element, which retains a linear character after homogenization.

(ii) Homogenization II: Homogenizing the bonded and frictional elements yields the macroscopic strength criterion for the rooted soil. The frictional element follows a nonlinear elastoplastic constitutive model; plastic dissipation occurs during loading, and its strength criterion F_f is nonlinear. In addition, during the homogenization of bonded and frictional elements, energy transfer caused by damage to the bonded elements must be considered; this includes both strain energy and dissipated energy. When some bonded elements fail and transform into frictional elements, additional energy dissipation occurs—analogous to the breakage of particle agglomerates—which is an energy-transfer process that must be accounted for in homogenization.

2. Strength Criterion

2.1. Strength Criterion for the Bonded Element (Homogenization Ⅰ)

Defining the mean stress $p$, generalized shear stress$q$, volumetric strain ${\epsilon }_v$ and generalized shear strain ${\epsilon }_s$ as follows:

$$\mathit{\sigma }=p\mathbb{I}+\mathit{s}$$

(1)

$$p=tr\mathit{\sigma }/3$$

(2)

$$q=\sqrt{3\left(\mathit{s}:\mathit{s}\right)/2}$$

(3)

$$\mathit{\epsilon }={\epsilon }_v\mathbb{I}+\mathit{e}$$

(4)

$${\epsilon }_v=tr\mathit{\epsilon }/3$$

(5)

$${\epsilon }_s=\sqrt{2\left(\mathit{e}:\mathit{e}\right)/3}$$

(6)

in which $\mathbb{I}$ is the unit tensor; $\mathit{\sigma }$, $\mathit{\epsilon }$ are respectively stress and strain tensor; $\mathit{s}$, $\mathit{e}$ are respectively deviatoric stress and deviatoric strain tensor.

Defining the strength for cemented soil matrix F^bm and roots inclusion F^r respectively as follows:

$$F^{bm}=q^{bm}-k^{bm}{p}^{bm}=0$$

(7)

$$F^r=q^r-k^r{p}^r=0$$

(8)

in which $q^{bm}$ and $p^{bm}$ denote the principal stress difference and mean stress, respectively, of the cemented-soil matrix at the failure state; $q^r$ and $p^r$ denote the principal stress difference and mean stress, respectively, of the root inclusions at the failure state; $k^{bm}$ and $k^r$ are the corresponding strength parameters. The strength parameters can be expressed as:

$$k^{bm}={\displaystyle \frac{q^{bm}}{p^{bm}}}={\displaystyle \frac{{3G}^{bm}}{K^{bm}}}{\displaystyle \frac{{\epsilon }_s^{bm}}{{\epsilon }_v^{bm}}}$$

(9)

$$k^r={\displaystyle \frac{q^r}{p^r}}={\displaystyle \frac{{3G}^r}{K^r}}{\displaystyle \frac{{\epsilon }_s^r}{{\epsilon }_v^r}}$$

(10)

where $K^{bm}$ and $G^{bm}$ are the bulk and shear moduli, respectively, of the cemented-soil matrix; ${\epsilon }_v^{bm}$ and ${\epsilon }_s^{bm}$ are the volumetric and shear strains, respectively, of the cemented-soil matrix at the failure state; $K^r$ and $G^r$ are the bulk and shear moduli, respectively, of the root inclusions; and ${\epsilon }_v^r$ and ${\epsilon }_s^r$ are the volumetric and shear strains, respectively, of the root inclusions at the failure state.

In a root–soil mixture, take a representative volume element (RVE). The volumes of its parts satisfy

$$V=V^b+V^f$$

(11)

$$V^b=V^{bm}+V^r$$

(12)

where $V$, $V^b$, $V^f$, $V^{bm}$, $V^r$ denote the volumes of the RVE, the bonded element, the frictional element, the cemented-soil matrix, and the root inclusion, respectively.

Define ${\mathsf{\lambda }}^f=V^f/V$ as the breakage ratio (volume fraction of the frictional element), and ${\mathsf{\lambda }}^r=V^r/V$ as the root volume fraction of the specimen (hereafter referred to as the root content). The volume fraction of roots within the bonded element, ${\lambda }^{br}$, is

$${\lambda }^{br}={\displaystyle \frac{V^r}{V^{bm}+V^r}}={\displaystyle \frac{{\lambda }^r}{1-{\lambda }^f}}$$

(13)

Within the bonded element, by the theory of volumetric homogenization, the stresses and strains of the cemented-soil matrix and the root inclusion satisfy:

$${\mathit{\sigma }}^b=\left(1-{\lambda }^{br}\right){\mathit{\sigma }}^{bm}+{\lambda }^{br}{\mathit{\sigma }}^r$$

(14)

$${\mathit{\epsilon }}^b=\left(1-{\lambda }^{br}\right){\mathit{\epsilon }}^{bm}+{\lambda }^{br}{\mathit{\epsilon }}^r$$

(15)

where ${\mathit{\sigma }}^b$ and ${\mathit{\epsilon }}^b$ are the stress and strain tensors of the bonded element; ${\mathit{\sigma }}^{bm}$ and ${\mathit{\epsilon }}^{bm}$ are those of the cemented-soil matrix; and ${\mathit{\sigma }}^r$ and ${\mathit{\epsilon }}^r$ are those of the root inclusion.

Treating the cemented soil as the matrix and the roots as inclusions within the bonded element, micromechanics and the Mori–Tanaka method give

$${\mathit{\epsilon }}^r={\mathbb{A}}^r:{\mathit{\epsilon }}^b$$

(16)

$${\mathbb{A}}^r={\left[\mathbb{I}+{\mathbb{S}}^{bm}:{\left({\mathbb{C}}^{bm}\right)}^{-1}:\left({\mathbb{C}}^r-{\mathbb{C}}^{\mathrm{bm}}\right)\right]}^{-1}$$

(17)

where ${\mathbb{C}}^{bm}$ and ${\mathbb{C}}^{bm}$ are the stiffness tensors of the cemented-soil matrix and the root inclusion, respectively. The Eshelby tensor of the bonded element with cemented soil as matrix is:

$${\mathbb{S}}^{\mathit{b}m}={\mathrm{S}}_1^{bm}\mathbb{J}+{\mathrm{S}}_2^{bm}\mathbb{K}$$

(18a)

$$S_1^{bm}={\displaystyle \frac{3K^{bm}}{3K^{bm}+4G^{bm}}}$$

(18b)

$$S_2^{bm}={\displaystyle \frac{6}{5}}{\displaystyle \frac{K^{bm}+2G^{bm}}{3K^{bm}+4G^{bm}}}$$

(18c)

Both the cemented-soil matrix and the root inclusion obey linear elastic constitutive laws (Hooke’s law). The constitutive relation and stiffness tensor ${\mathbb{C}}^{bm}$ of the cemented-soil matrix are:

$$d{\mathit{\sigma }}^{bm}={\mathbb{C}}^{bm}:d{\mathit{\epsilon }}^{bm}$$

(19)

$${\mathbb{C}}^{bm}=3K^{bm}\mathbb{J}+2G^{bm}\mathbb{K}$$

(20)

The constitutive relation and stiffness tensor ${\mathbb{C}}^r$ of the roots are

$$d{\mathit{\sigma }}^r={\mathbb{C}}^r:d{\mathit{\epsilon }}^r$$

(21)

$${\mathbb{C}}^r=3K^r\mathbb{J}+2G^r\mathbb{K}$$

(22)

Solving Eqs. (13)–(22) jointly yields：

$${\mathit{\sigma }}^{bm}={\mathbb{D}}^{bm}:{\mathit{\sigma }}^b$$

(23)

$${\mathbb{D}}^{bm}={\left[\left(1-{\lambda }^{br}\right)\mathbb{I}+{\lambda }^{br}{\mathbb{C}}^r:\mathbb{A}:{\left({\mathbb{C}}^{bm}\right)}^{-1}\right]}^{-1}$$

(24)

$${\mathit{\sigma }}^r={\mathbb{D}}^r:{\mathit{\sigma }}^b$$

(25)

$${\mathbb{D}}^r=\left[\mathbb{I}-\left(1-{\lambda }^{br}\right){\mathbb{D}}^{bm}\right]/{\lambda }^{br}$$

(26)

where ${\mathbb{D}}^{bm}$ and ${\mathbb{D}}^r$ are the local stress concentration tensors of the cemented-soil matrix and the roots, respectively.

From tensor algebra, any isotropic fourth-order tensor can be expressed as a linear combination of the spherical and deviatoric fourth-order identity tensors. By reducing the fourth-order tensor expressions to second-order matrix form, and since both the cemented-soil matrix and the root inclusions follow linear elastic constitutive models, the resulting second-order matrices are diagonal. After computing the above tensors and simplifying, they can be written as:

$${\mathbb{S}}^{\mathit{b}m}=\left(S_1^{bm}，S_2^{bm}\right)$$

(27)

$${\mathbb{C}}^{bm}=\left(K^{bm}，3G^{bm}\right)$$

(28)

$${\mathbb{C}}^r=\left(K^r，3G^r\right)$$

(29)

$${\mathbb{A}}^r=\left({\mathbb{A}}_1^r，{\mathbb{A}}_2^r\right)=\left({\displaystyle \frac{K^{bm}}{K^{bm}+S_1^{bm}\left({K^r-K}^{bm}\right)}}，{\displaystyle \frac{G^{bm}}{G^{bm}+S_2^{bm}\left({G^r-G}^{bm}\right)}}\right)$$

(30)

$${\mathbb{D}}^{bm}=\left({\mathbb{D}}_1^{bm}，{\mathbb{D}}_2^{bm}\right)$$

(31a)

$${\mathbb{D}}_1^{bm}={\displaystyle \frac{K^{bm}+S_1^{bm}\left({K^r-K}^{bm}\right)}{K^{bm}+\left(S_1^{bm}+{\lambda }^{br}-S_1^{bm}{\lambda }^{br}\right)\left({K^r-K}^{bm}\right)}}$$

(31b)

$${\mathbb{D}}_2^{bm}={\displaystyle \frac{G^{bm}+S_2^{bm}\left({G^r-G}^{bm}\right)}{G^{bm}+\left(S_2^{bm}+{\lambda }^{br}-S_2^{bm}{\lambda }^{br}\right)\left({G^r-G}^{bm}\right)}}$$

(31c)

$${\mathbb{D}}^r=\left({\mathbb{D}}_1^r，{\mathbb{D}}_2^r\right)=\left({\displaystyle \frac{1-\left(1-{\lambda }^{br}\right){\mathbb{D}}_1^{bm}}{{\lambda }^{br}}}，{\displaystyle \frac{1-\left(1-{\lambda }^{br}\right){\mathbb{D}}_2^{bm}}{{\lambda }^{br}}}\right)$$

(32)

At the critical state, Eq.(14) can be rewritten as the following expression through the mean stress and generalized shear stress:

$$p^b=\left(1-{\lambda }^{br}\right)p^{bm}+{\lambda }^{br}{p}^r$$

(33)

$$q^b=\left(1-{\lambda }^{br}\right)q^{bm}+{\lambda }^{br}{q}^r$$

(34)

Substituting Eq. (7) and Eq. (8) into Eq. (34), we can obtain

$$q^b=\left(1-{\lambda }^{br}\right)k^{bm}{p}^{bm}+{\lambda }^{br}{k}^r{p}^r$$

(35)

From Eq. (23) and Eq. (25), we have

$$p^{bm}={\mathbb{D}}_1^{bm}{p}^b$$

(36)

$$p^r={\mathbb{D}}_1^r{p}^b$$

(37)

Substituting the above two equations into Eq. (35), we have

$$q^b-\left[\left(1-{\lambda }^{br}\right)k^{bm}{\mathbb{D}}_1^{bm}+{\lambda }^{br}{k}^r{\mathbb{D}}_1^r\right]p^b=0$$

(38)

The above equation is the strength criterion of the bonded element, which can be rewritten as the following expression:

$$F^b=q^b-k^b{p}^b=0$$

(39)

$$k^b=\left(1-{\lambda }^{br}\right)k^{bm}{\mathbb{D}}_1^{bm}+{\lambda }^{br}{k}^r{\mathbb{D}}_1^r$$

(40)

It can be seen that the bonded element shares the same mathematical form as the cemented-soil matrix and the root inclusions. Since both the cemented-soil matrix and the root inclusions follow linear elastic constitutive models, the bonded element is a linear combination of the cemented-soil matrix and the root inclusions, which is consistent with linear superposition in elasticity.

The bonded element is a linear elastic material. Prior to damage, it undergoes only elastic deformation without plastic deformation. During loading, only the strain energy of the bonded element changes, with no plastic dissipation. Therefore, at the point of failure, the strain energy ${\phi }^b$ and dissipation energy ${\pi }^b$ of the bonded element are

$${\phi }^b={\displaystyle \frac{1}{2}}{\mathit{\epsilon }}^b:{\mathbb{C}}^b:{\mathit{\epsilon }}^b={\displaystyle \frac{1}{2}}{K}^b{\left({\epsilon }_v^b\right)}^2+{\displaystyle \frac{3}{2}}{G}^b{\left({\epsilon }_s^b\right)}^2$$

(41)

$${\pi }^b=0$$

(42)

Here, ${\epsilon }_v^b$ and ${\epsilon }_s^b$ are the volumetric and shear strains of the bonded element at failure; $K^b$ and $G^b$ are the bulk and shear moduli of the bonded element, which can be obtained via the Mori–Tanaka method or other mesomechanics-based approach:

$$K^b=K^{bm}\left[1+{\displaystyle \frac{{\lambda }^{br}(\frac{K^r}{K^{bm}}-1)}{1+(1-{\lambda }^{br})S_1^{bm}(\frac{K^r}{K^{bm}}-1)}}\right]$$

(43)

$$G^b=G^{bm}\left[1+{\displaystyle \frac{{\lambda }^{br}(\frac{G^r}{G^{bm}}-1)}{1+\left(1-{\lambda }^{br}\right)S_2^{bm}(\frac{G^r}{G^{bm}}-1)}}\right]$$

(44)

The linear elastic constitutive model of the bonded element is

$${\mathbb{C}}^b=\left(K^b，3G^b\right)$$

(45)

The nonlinear function $Y^b$ of the bonded element can be written as

$$Y^b=stat\left({\pi }^b-{\phi }^b\right)$$

(46)

Here, stat(·) denotes the stationarity operator, i.e., taking the extremum of the function in parentheses. By setting the partial derivatives of the variables to zero at the stationary point and solving the resulting coupled equations, the solution can be obtained.

2.2. Strength Criterion for the Frictional Element

From a micromechanical perspective, failure of granular soils is primarily caused by frictional sliding and shear slip (Zhao et al., 2018), exhibiting pronounced elastoplastic deformation characteristics. The strength of the frictional element is nonlinear, and the failure of soil particles can be well described by either the Mohr–Coulomb strength criterion or an elastoplastic strength criterion. For greater generality, this paper adopts an elliptical strength criterion to represent the strength of the frictional element, namely:

$$F^f=1-{\displaystyle \frac{{\left(p^f+S\right)}^2}{a}}+{\displaystyle \frac{{\left(q^f\right)}^2}{a{M}^2}}=0$$

(47)

where $q^f$ and $p^f$ are the principal stress difference and mean stress of the frictional element at the critical state; S and M are strength parameters; a is a regularization parameter. As a → 0, $F^f=q^f-M\left(p^f+S\right)$, and the elliptical strength criterion degenerates to the D–P criterion.

Because the frictional element undergoes nonlinear elastoplastic deformation during loading, in addition to the conversion of strain energy, plastic dissipation also occurs. The strain energy of the frictional element is nonlinear and difficult to solve directly. Here we employ the Linear Comparison Composite (LCC) method (Castaneda, 2002; Lopez-Pamies and Ponte Castaneda, 2004; Ortega et al., 2011; Wang et al., 2021; Zhou and Meschke,2018) to linearize the constitutive model and obtain an approximate solution for the strain energy. The LCC method refers to linearizing nonlinear features when the components of a composite share similar underlying structures; the nonlinear function Y is a “stationary estimate” of the difference between the microscale maximum dissipation energy $\pi$ and the strain energy ${\phi }^b$ of each constituent.

The constitutive model of the frictional element, after LCC-based approximation, is expressed as:

$${\mathit{\sigma }}^f={\mathit{\tau }}^f+C^f:{\mathit{\epsilon }}^f={\mathbb{C}}^f:{\mathit{\epsilon }}^f$$

(48)

where ${\mathit{\sigma }}^f$ and ${\mathit{\epsilon }}^f$ are the stress and strain tensors of the frictional element; ${\mathit{\tau }}^f$ is the prestress tensor of the frictional element; $C^f$ is the elastic stiffness tensor of the friction element, and ${\mathbb{C}}^f$ is the elastic stiffness of the friction element after linearization; they satisfy, respectively, the following:

$$C^f=\left(K^f，3G^f\right)$$

(49)

$${\mathbb{C}}^f=\left({\mathbb{K}}^f，3{\mathbb{G}}^f\right)$$

(50)

At the critical failure point, the prestress satisfies:

$${\mathit{\tau }}^f=C^{\tau }:{\mathit{\epsilon }}^f$$

(51)

where $C^{\tau }$ is the prestress stiffness tensor or the linearly corrected stiffness tensor of the friction element, hence:

$$C^{\tau }={\mathbb{C}}^f-C^f=\left(K^{\mathit{\tau }}，3G^{\mathit{\tau }}\right)$$

(52)

Regarding the prestress, some researchers directly assume a specific form, while others determine it via limit analysis by setting partial derivatives to zero. In this paper, it is obtained using the macro–meso energy equivalence principle and the Mori–Tanaka method.

Then the strain energy ${\phi }^f$ of the friction element can be expressed as:

$${\phi }^f={\displaystyle \frac{1}{2}}{\mathit{\epsilon }}^f:C^f:{\mathit{\epsilon }}^f+{\mathit{\tau }}^f:{\mathit{\epsilon }}^f={\displaystyle \frac{1}{2}}{K}^f{\left({\epsilon }_v^f\right)}^2+{\displaystyle \frac{3}{2}}{G}^f{\left({\epsilon }_s^f\right)}^2+{\tau }_v^f{\epsilon }_v^f+{\tau }_s^f{\epsilon }_s^f$$

(53)

where $K^f$ and $G^f$ denote the elastic bulk modulus and elastic shear modulus of the friction element; ${\epsilon }_v^f$ and ${\epsilon }_v^f$ are the volumetric strain and shear strain of the friction element at failure; ${\tau }_v^f$ and ${\tau }_s^f$ are the volumetric prestress and shear prestress of the friction element.

For the dissipated energy of the friction element, we first define a support function to represent, from a mesoscopic perspective, the maximum energy dissipation associated with the resistance of geomaterials to deformation (Dormieux et al., 2002, 2017; Maghous et al., 2009), namely:

$${\pi }^f=sup\left({\mathit{\sigma }}^f:{\mathit{\epsilon }}^{fp}\right)=p^f{\epsilon }_v^{fp}+q^f{\epsilon }_s^{fp}$$

(54)

where ${\epsilon }_v^{fp}$ and ${\epsilon }_s^{fp}$ are the plastic volumetric strain and plastic shear strain of the friction element at failure. The plastic strain of the friction element can be written as:

$${\mathit{\epsilon }}^{fp}={\mathit{\epsilon }}^f-{\mathit{\epsilon }}^{fe}={\left({\mathbb{C}}^f\right)}^{-1}:{\mathit{\sigma }}^f-{\left(C^f\right)}^{-1}:{\mathit{\sigma }}^f=\left[{\left({\mathbb{C}}^f\right)}^{-1}-{\left(C^f\right)}^{-1}\right]:{\mathit{\sigma }}^f$$

(55)

The nonlinear function $Y^f$ of the friction element is given by

$$Y^f=stat\left({\pi }^f-{\phi }^f\right)$$

(56)

2.3. Dissipated Energy During the Breakage Process

During the gradual failure of a bonded element and its transition into a friction element, energy is both converted and dissipated; this process can be characterized by a breakage stress (Wang et al., 2021). It can be considered that, in the transition from a bonded element to a friction element, energy dissipation predominates. Define the breakage stress ${\mathit{\sigma }}^{bf}$ as:

$${\mathit{\sigma }}^{bf}={\mathit{\sigma }}^f-{\mathit{\sigma }}^b$$

(57)

The bonded element deforms linearly elastically, and only the friction element undergoes plastic deformation during loading. Therefore, the plastic strain of the breakage process ${\mathit{\epsilon }}^{bfp}$ is taken to be the plastic strain of the friction element, i.e.

$${\mathit{\epsilon }}^{bfp}={\mathit{\epsilon }}^{fp}$$

(58)

The damage dissipated energy can be written as

$${\pi }^{bf}=sup\left({\mathit{\sigma }}^{bf}:{\mathit{\epsilon }}^{fp}\right)=\left(p^f-p^b\right){\epsilon }_v^{fp}+\left(q^f-q^b\right){\epsilon }_s^{fp}$$

(59)

Substituting the constitutive relation of the bonded element into the above gives:

$${\pi }^{bf}=\left(p^f-K^b{\epsilon }_v^b\right){\epsilon }_v^{fp}+\left(q^f-3G^b{\epsilon }_s^b\right){\epsilon }_s^{fp}$$

(60)

The nonlinear function for the damage process is

$$Y^{bf}=stat\left({\pi }^{bf}\right)$$

(61)

2.4. Macroscopic Strength Criterion (Homogenization II)

From a macroscopic perspective, the constitutive model of root-reinforced soil can also be treated using the LCC method, which approximately yields the macroscopic strain energy and dissipated energy. The linearized constitutive model of the RVE for root-reinforced soil is expressed as:

$$\mathit{\sigma }=\mathit{\tau }+C^{hom}:\mathit{\epsilon }={\mathbb{C}}^{hom}:\mathit{\epsilon }$$

(62)

in which $\mathit{\sigma }$, $\mathit{\epsilon }$, $\mathit{\tau }$ are the stress tensor, strain tensor, and prestress tensor of the RVE, respectively; $C^{hom}$ is the macroscopic elastic stiffness tensor of the RVE without considering plastic deformation; ${\mathbb{C}}^{hom}$ is the RVE macroscopic stiffness tensor corrected by the LCC method.

From homogenization theory and based on a binary medium model, we have

$$\mathit{\sigma }=\left(1-{\mathsf{\lambda }}^f\right){\mathit{\sigma }}^b+{\mathsf{\lambda }}^f{\mathit{\sigma }}^f$$

(63)

$$\mathit{\epsilon }=\left(1-{\mathsf{\lambda }}^f\right){\mathit{\epsilon }}^b+{\mathsf{\lambda }}^f{\mathit{\epsilon }}^f$$

(64)

Substituting the above into the constitutive relations of the RVE, the bonded element, and the friction element, we obtain:

$$\mathit{\tau }+C^{hom}:\mathit{\epsilon }=\left(1-{\mathsf{\lambda }}^f\right){\mathbb{C}}^b:{\mathit{\epsilon }}^b+{\mathsf{\lambda }}^f\left({\mathit{\tau }}^f+C^f:{\mathit{\epsilon }}^f\right)$$

(65)

By substituting Eq. (51) into Eq. (65) and comparing, we get:

$$C^{hom}:\mathit{\epsilon }=\left(1-{\mathsf{\lambda }}^f\right){\mathbb{C}}^b:{\mathit{\epsilon }}^b+{\mathsf{\lambda }}^f{C}^f:{\mathit{\epsilon }}^f$$

(66)

$$\mathit{\tau }={\mathsf{\lambda }}^f{C}^{\tau }:{\mathit{\epsilon }}^f$$

(67)

By computation via the Mori–Tanaka method, the macroscopic elastic stiffness tensor $C^{hom}$ can be expressed as follows:

$$C^{hom}=\left(K^{hom},3G^{hom}\right)$$

(68a)

$$K^{hom}=K^b\left[1+{\displaystyle \frac{{\lambda }^f(\frac{K^f}{K^b}-1)}{1+(1-{\lambda }^f)S_1^b(\frac{K^f}{K^b}-1)}}\right]$$

(68b)

$$G^{hom}=G^b\left[1+{\displaystyle \frac{{\lambda }^f(\frac{G^f}{G^b}-1)}{1+\left(1-{\lambda }^f\right)S_2^b(\frac{G^f}{G^b}-1)}}\right]$$

(68c)

in which $S_1^f$ and $S_2^f$ are the coefficients of the Eshelby tensor ${\mathbb{S}}^f$ for the case with the friction element as the matrix and the cemented element as the inclusion, satisfying the following:

$${\mathbb{S}}^b={\mathrm{S}}_1^b\mathbb{J}+{\mathrm{S}}_2^b\mathbb{K}$$

(69a)

$${\mathrm{S}}_1^b={\displaystyle \frac{3K^b}{3K^b+4G^b}}$$

(69b)

$${\mathrm{S}}_2^b={\displaystyle \frac{6}{5}}{\displaystyle \frac{K^b+2G^b}{3K^b+4G^b}}$$

(69c)

Similarly, the macroscopic stiffness tensor ${\mathbb{C}}^{hom}$ after linearization can be written as:

$${\mathbb{C}}^{hom}=\left({\mathbb{K}}^{hom},3G^{hom}\right)$$

(70a)

$${\mathbb{K}}^{hom}=K^b\left[1+{\displaystyle \frac{{\lambda }^f(\frac{{\mathbb{K}}^f}{K^b}-1)}{1+(1-{\lambda }^f)S_1^b(\frac{{\mathbb{K}}^f}{K^b}-1)}}\right]$$

(70b)

$${\mathbb{G}}^{hom}=G^b\left[1+{\displaystyle \frac{{\lambda }^f(\frac{{\mathbb{G}}^f}{G^b}-1)}{1+\left(1-{\lambda }^f\right)S_2^b(\frac{{\mathbb{G}}^f}{G^b}-1)}}\right]$$

(70c)

The strain tensors of the bonded and friction elements satisfy:

$${\mathit{\epsilon }}^b=\mathbb{A}:{\mathit{\epsilon }}^f$$

(71)

where A is the local strain concentration tensor, expressed as

$$\mathbb{A}=\mathbb{I}+{\mathbb{S}}^b:{\left({\mathbb{C}}^b\right)}^{-1}:\left({\mathbb{C}}^b-{\mathbb{C}}^{\mathrm{f}}\right)$$

(72)

The explicit expression of the local strain concentration tensor is

$$\mathbb{A}=\left({\mathbb{A}}_1，{\mathbb{A}}_2\right)=\left(1+{\mathrm{S}}_1^b\left(1-{\displaystyle \frac{{\mathbb{K}}^f}{K^b}}\right)，1+{\mathrm{S}}_2^b\left(1-{\displaystyle \frac{{\mathbb{G}}^f}{G^b}}\right)\right)$$

(73)

Substituting Eq.(71) into Eq.(64), we can obtain:

$$\mathit{\epsilon }=\left[\left(1-{\lambda }^f\right)\mathbb{A}+{\lambda }^f\mathbb{I}\right]:{\mathit{\epsilon }}^f=\mathbb{B}:{\mathit{\epsilon }}^f$$

(74)

The macroscopic strain energy is expressed as

$${\phi }^{hom}={\displaystyle \frac{1}{2}}\mathit{\epsilon }:C^{hom}:\mathit{\epsilon }+\mathit{\tau }:\mathit{\epsilon }$$

(75)

According to the macro–meso energy equivalence principle, the macroscopic and mesoscopic strain energy densities are equal, namely:

$${\phi }^{hom}=\left(1-{\lambda }^f\right){\phi }^b+{\lambda }^f{{\phi }^f+\lambda }^f{\phi }^{bf}$$

(76)

After linearization by the LCC method, the stiffness tensors and the stress/strain tensors of the bonded element, the friction element, and the RVE are all diagonal. Substituting the relevant expressions into Eqs. (75) and (76) reduces them to

$${\displaystyle \frac{1}{2}}\mathbb{B}:C^{hom}:\mathbb{B}+{\lambda }^f{C}^{\tau }:\mathbb{B}={\displaystyle \frac{1}{2}}\left(1-{\lambda }^f\right)\mathbb{A}:{\mathbb{C}}^{\mathrm{b}}:\mathbb{A}+{\lambda }^f\left({\displaystyle \frac{1}{2}}{C}^f+C^{\tau }\right)$$

(77)

The above is a quadratic equation in ${\mathbb{C}}^{\mathrm{f}}$ or $C^{\tau }$. It can be solved using Eq. (54), and we obtain

$$D_1{\left({\mathbb{C}}^b-{\mathbb{C}}^{\mathrm{f}}\right)}^2+D_2\left({\mathbb{C}}^b-{\mathbb{C}}^{\mathrm{f}}\right)+D_3=0$$

(78a)

$$D_1=\left[C^{hom}\left(1-{\lambda }^f\right)-{\mathbb{C}}^{\mathrm{b}}\right]\left(1-{\lambda }^f\right){\left({\displaystyle \frac{{\mathbb{S}}^b}{{\mathbb{C}}^b}}\right)}^2$$

(78b)

$$D_2=2\left[\left(C^{hom}-{\mathbb{C}}^{\mathrm{b}}\right)\left(1-{\lambda }^f\right){\displaystyle \frac{{\mathbb{S}}^b}{{\mathbb{C}}^b}}+{\lambda }^f\right]$$

(78c)

$$D_3=C^{hom}-\left(1+{\lambda }^f\right){\mathbb{C}}^b+{\lambda }^f{C}^f$$

(78d)

Solve for ${\mathbb{C}}^b-{\mathbb{C}}^{\mathrm{f}}$ using the quadratic formula. According to Eq. (37), the prestress tensor ${\mathit{\tau }}^f$ of the friction element and its stiffness tensor $C^{\tau }$ can then be determined. Based on the condition ${\mathbb{C}}^{\mathrm{f}}\ge 0$ and reasonable bounds, discard the root without physical meaning and use the physically meaningful root in subsequent calculations.

According to the LCC method and limit analysis theory, the maximum macroscopic dissipated energy of the homogenized composite is given by a stationary estimate of the macroscopic strain energy plus the volume-fraction-weighted sum of the nonlinear functions. The dissipated energy of the macroscopic RVE is

$${\pi }^{hom}=stat\left[{\phi }^{hom}+\left({1-\lambda }^f\right)Y^b+{\lambda }^f{Y}^f+{\lambda }^f{Y}^{bf}\right]$$

(79)

Regarding the nonlinear functions, some previous researchers computed the stationary function at each step separately, substituted the results into Eq. (79), and then performed another stationarity operation. This approach is questionable because each step’s stationary value corresponds to a particular stress–strain state, and these states may not be mutually consistent. A more reasonable approach is to omit the intermediate stepwise calculations and evaluate the extremum of the stationary functional in a single step using Eq. (79).

The above expression implicitly embodies the macro–meso energy equivalence principle. Substituting Eqs. (56), (61), and (76), which enforce equality of macro and meso strain energies, into Eq. (79), the right-hand side becomes the sum of the dissipated energies at the mesoscopic scale, namely:

$${\pi }^{hom}=stat\left({\lambda }^f{\pi }^f+{\lambda }^f{\pi }^{bf}\right)$$

(80)

Substituting Eqs.(54), (55) and (60) into Eq. (80), we can obtain the following by manipulation:

$${\pi }^{hom}={\lambda }^f\left[{\left(p^f\right)}^2\left({\displaystyle \frac{1}{{\mathbb{K}}^f}}-{\displaystyle \frac{1}{K^f}}\right)\left(2-{\displaystyle \frac{{{\mathbb{A}}_{1}K}^b}{{\mathbb{K}}^f}}\right)+{\left(q^f\right)}^2\left({\displaystyle \frac{1}{{3\mathbb{G}}^f}}-{\displaystyle \frac{1}{{3G}^f}}\right)\left(2-{\displaystyle \frac{{{\mathbb{A}}_{2}G}^b}{{\mathbb{G}}^f}}\right)\right]$$

(81)

The above expression must satisfy the strength criterion of the friction element; that is, find the stationary value of the expression subject to the condition in Eq. (32). Substituting Eq. (47) into Eq. (81), solve for the point where the partial derivative is zero.

$${\displaystyle \frac{{\partial \pi }^{hom}}{\partial p^f}}=0$$

(82)

Similarly, the dissipated energy at the macroscopic scale can be defined as:

$${\pi }^{hom}=\sup \left(\mathit{\sigma }:{\mathit{\epsilon }}^p\right)=P{\epsilon }_v^p+Q{\epsilon }_s^p$$

(83)

where P and Q are the mean stress and the principal stress difference of the RVE of root-reinforced soil at the critical failure state; ${\mathit{\epsilon }}^p$ is the plastic strain tensor at the critical failure state, and ${\epsilon }_v^p$ and ${\epsilon }_s^p$ are the plastic volumetric strain and plastic shear strain, which can be written as:

$${\mathit{\epsilon }}^p=\mathit{\epsilon }-{\mathit{\epsilon }}^p=\left[{\left({\mathbb{C}}^{hom}\right)}^{-1}-{\left(C^{hom}\right)}^{-1}\right]:\mathit{\sigma }$$

(84)

The macroscopic dissipated energy is expressed as

$${\pi }^{hom}=P^2\left({\displaystyle \frac{1}{{\mathbb{K}}^{hom}}}-{\displaystyle \frac{1}{K^{hom}}}\right)+Q^2\left({\displaystyle \frac{1}{{3\mathbb{G}}^{hom}}}-{\displaystyle \frac{1}{{3G}^{hom}}}\right)$$

(85)

According to the macro–meso energy equivalence principle, the dissipated energies at the macroscopic and mesoscopic scales are equal, from which the macroscopic strength criterion for root-reinforced soil is ultimately obtained:

$$F^{hom}={\displaystyle \frac{P^2}{X_1}}+{\displaystyle \frac{Q^2}{X_2}}-1=0$$

(86)

in which

$$X_1={\displaystyle \frac{{\mathbb{K}}^{hom}{K}^{hom}R}{K^{hom}-{\mathbb{K}}^{hom}}}$$

(87a)

$$X_2=3{\displaystyle \frac{{\mathbb{G}}^{hom}{G}^{hom}R}{G^{hom}-{\mathbb{G}}^{hom}}}$$

(87b)

$$R={{\lambda }^{f}M}^{2}N\left({\displaystyle \frac{{S^{2}M}^{2}N}{L+M^{2}N}}-a\right)$$

(87c)

$$L=\left({\displaystyle \frac{1}{{\mathbb{K}}^f}}-{\displaystyle \frac{1}{K^f}}\right)\left(2-{\displaystyle \frac{{{\mathbb{A}}_{1}K}^b}{{\mathbb{K}}^f}}\right)$$

(87d)

$$N=\left({\displaystyle \frac{1}{{3\mathbb{G}}^f}}-{\displaystyle \frac{1}{{3G}^f}}\right)\left(2-{\displaystyle \frac{{{\mathbb{A}}_{2}G}^b}{{\mathbb{G}}^f}}\right)$$

(87e)

3. Verification of Strength Criterion

3.1. Parameter Determination

To predict the strength of root-reinforced soil under freeze–thaw cycling, it is necessary to know the parameters of the roots, the cemented soil matrix, and the friction elements, as well as the breakage parameter. The bulk modulus $K^r$ and shear modulus $G^r$ of the roots can be determined directly or indirectly through laboratory tests. In the binary medium model, it can be considered that at the initial stage of loading the specimen is mainly composed of bonded elements, and the breakage ratio is low. Therefore, the bulk modulus $K^{bm}$ and shear modulus $G^{bm}$ of the cemented soil matrix can be approximately determined from the initial slopes of the ${\epsilon }_v-{\sigma }_m$ and ${{\epsilon }_s-\sigma }_s$ curves at 0.2% axial strain in conventional triaxial tests on specimens without roots and without freeze–thaw cycling.

After completion of the triaxial test, the bonded elements are essentially completely damaged; the remolded soil after freeze–thaw and triaxial testing consists predominantly of friction elements. The parameters of the friction elements can be obtained from triaxial tests on the remolded specimens after freeze–thaw and triaxial testing. The elastic bulk and shear moduli $K^f$ and $G^f$ of the friction elements can be approximately determined from the initial slopes of the ${\epsilon }_v^f-{\sigma }_m^f$ and ${\epsilon }_s^f-{\sigma }_s^f$ curves of these remolded specimens. The strength parameters a, S, and M of the friction elements can be obtained by fitting the strengths measured on the remolded specimens after freeze–thaw and triaxial testing.

For the breakage ratio, many studies have shown that its evolution follows a Weibull distribution function. We recommend the following equation:

$${\lambda }^f=1-exp\left[-\alpha {\left({\displaystyle \frac{P}{P_a}}\right)}^{\beta }\right]$$

(88)

where α and β are damage-fraction parameters that can be determined by inverse analysis, and Pa is the standard atmospheric pressure.

3.2. Experimental Verification

The test soil was extracted from Hailuogou, Sichuan Province, China, at an elevation of about 3000 m, from a depth of 0.5–1.0 m below the ground surface (Luo et al., 2023). The soil is a silty sand, and the roots used were from rowan (Sorbus). The roots were about 5 cm long with a diameter of about 5 mm; the mass of a single root was about 1.4 g, corresponding to a mass fraction of 0.292% (on a dry-mass basis), a volume fraction of 0.262%, and a shearing-area fraction of 0.254%. Specimens were prepared at a target dry density of 1.28 g/cm³, with a diameter of 61.8 mm and a height of 125 mm. After specimen preparation, they were vacuumed for 2 hours and then submerged to saturate in water for more than 12 hours. Specimens not subjected to freeze–thaw cycles were tested directly in triaxial compression; for specimens requiring freeze–thaw cycles, they were promptly placed in water to remove bubbles and transferred under water into lidless transparent plastic jars. The jars containing the specimens were labeled and placed in a fully automatic low-temperature freeze–thaw machine to begin cycling. During the freeze–thaw cycles, water could be replenished through the porous stones at the top and bottom ends of the specimen. According to the temperature of the field, one freeze–thaw cycle consisted of cooling from room temperature to −15°C and holding for 12 hours, then heating to 20°C and holding for 12 hours. After the prescribed number of cycles, the specimen was mounted on a triaxial apparatus and tested under room-temperature consolidated drained conditions, with an axial loading rate of 0.1 mm/min and confining pressures of 25, 50, 100, and 200 kPa. The test was terminated when the axial strain reached 15%, which was taken as failure. The test matrix included confining pressures of 25 kPa, 50 kPa, 100 kPa, and 200 kPa; freeze–thaw cycle counts of 0, 1, 5, 15; and root counts of none, 1 root, and 3 roots.

Based on the triaxial test results, the bulk and shear moduli of the roots, the cemented soil matrix, and the frictional element were determined. From the strengths of remolded specimens after freeze–thaw and triaxial testing, the strength parameters a, S, and M of the frictional element were fitted. A computational program was then developed to verify the strength criterion for root-reinforced soils. First, the parameters of the bonded element were obtained using the Mori–Tanaka method; next, by invoking strain-energy conservation, the prestress and the parameters of the frictional element were determined; then, from the dissipated energy, the parameters of the macroscopic strength criterion were derived; finally, the model predictions were obtained according to the relationships under triaxial testing conditions. Under triaxial stress, the mean stress P and the principal stress difference Q are:

$$P={\sigma }_3+{\displaystyle \frac{1}{3}}\left({\sigma }_1-{\sigma }_3\right)$$

(89)

$$Q={\sigma }_1-{\sigma }_3$$

(90)

Following the parameter-determination procedure described above, the parameters of the strength criterion for root-reinforced soil under freeze–thaw action were ultimately identified by fitting the test data and performing inverse analysis; these are listed in Table 1. Among them, the damage-rate parameter is affected by the number of freeze–thaw cycles: as the number of cycles increases, the parameter α increases and the damage rate rises accordingly, which corroborates the conclusion that freeze–thaw cycling exacerbates damage to the bonded element.

Based on the parameters in Table 1, the strength criterion model for root-reinforced soil was used to predict strengths under different confining pressures, numbers of freeze–thaw cycles, and root contents, as shown in Figure 3 and Figure 4. Figure 3 and Figure 4 indicate that the model-predicted trends of strength for root-reinforced soil under freeze–thaw action are consistent with the experimental results: strength increases with increasing confining pressure, decreases with increasing numbers of freeze–thaw cycles, and increases with higher root content. A comparison between the model predictions and the test data shows that the model can reproduce, with good accuracy, both the magnitude of strength and its variation under freeze–thaw cycling, with small prediction errors.

The proposed strength criterion model is derived on a micromechanical scale and can reasonably capture the micromechanical failure mechanisms of a three-constituent composite. The model requires relatively few parameters, and their physical meanings are clear; except for the damage-rate parameter, the parameters for the cemented soil matrix, the frictional element, and the roots can all be obtained from laboratory tests. In summary, the micromechanics-based strength criterion for root-reinforced soils, which accounts for the effects of freeze–thaw cycling, can provide a reliable basis for engineering practice.

4. Conclusions

Based on a micromechanical model, the RVE of root-reinforced soil is idealized as a three-constituent composite: a matrix of cemented elements—each composed of a cemented-soil matrix with root inclusions—and inclusions of frictional elements. A strength-criterion model for this three-phase composite was established. With the aid of the LCC method, the Mori–Tanaka scheme, and two-phase medium theory, a strength criterion for root-reinforced soil was developed while accounting for the number of freeze–thaw cycles. The model was then used to predict strengths and compared against experimental results. The main findings are summarized as follows:

(i) In the strength-criterion model for root-reinforced soil, the cemented element is treated as the matrix and the frictional element as the inclusion; the cemented element itself consists of a cemented-soil matrix with root inclusions. From a micromechanical perspective, this constructs a strength-criterion model for a three-constituent composite.

(ii) The constitutive model of the frictional element was linearized using the LCC method. Employing the two-phase medium model, micromechanics, macro–micro energy-equivalence principles, and limit analysis, a macroscopic strength criterion was derived for the root-reinforced soil RVE, thereby revealing the damage mechanisms at the microscale. In doing so, the work improves upon certain shortcomings in previous applications of limit analysis used to obtain the stationary function, and, using the Mori–Tanaka method together with macro–micro strain-energy equivalence, inversely determines the prestress and the moduli of the linearized frictional element.

(iii) The parameters of the strength criterion were identified through laboratory test data and inverse analysis. Using the proposed strength-criterion model, the strengths under different confining pressures, numbers of freeze–thaw cycles, and root contents were predicted with good accuracy, consistent with the strength trends observed in root-reinforced soil tests.