Introduction

Recently, a large number of infrastructure projects such as railways, highways, airports and water conservancy hubs have been built in areas with seasonal and permanent frozen soils. With the rise of temperature, slopes in permafrost degradation areas and seasonal permafrost areas may be affected by freeze–thaw action, which may lead to their unique freeze–thaw collapse1. Slope protection by vegetation is to maintain the stability of rock and soil slope, reduce slope erosion, prevent soil and water loss, and protect the ecological environment by using the principle of vegetation drainage and soil fixation. Meanwhile, it also has the advantages of low cost and superior economy2. Therefore, it is necessary to study the influence of freeze–thaw cycle and root system on soil mechanical properties.

Freeze–thaw cycles can lead to significant changes in soil engineering properties. Many scholars’3,4,5,6,7,8,9,10,11 studies about fine grained soils (clay and silt) show that: freeze–thaw cycles will affect soil structure, density, porosity, water redistribution, and the strength, compressibility and pore water pressure. Freeze–thaw cycles can also cause changes in soil microstructure, resulting in changes in soil properties in many ways, such as permeability changes, often up to several levels of magnitude12. Although many scholars have made fruitful achievements in the research of frozen soil mechanics in recent years, the results are not uniform. In terms of stress–strain relationship, Alkire et al.13 and Ono et al.14 respectively conducted experiments on reshaped silty soil and reshaped clay, and the stress–strain curve of soil undergoing a freeze–thaw cycle was above that of unfreeze-thaw soil, indicating that freeze–thaw action enhanced soil's resistance to shear stress. On the contrary, tests conducted by Graham et al.15 show that the soil compression curve after freeze–thaw is located below unbroken soil without freeze–thaw cycles. The test results of Wang et al.16 on Qinghai-Tibet clay show that, at low confining pressure, the stress–strain relationship of samples after freeze–thaw cycles gradually changes from strain-softening type to strain-hardening type, that is, the breakage mode of samples changes from brittle breakage to plastic breakage due to freeze–thaw action, but if at high confining pressure, the strain hardening was observed before and after freeze–thaw. The results of freeze–thaw cycles effect on soil strength are different. Broms17, Leroueil18 and Qin et al.19 found that soil strength decreased after freeze–thaw cycles. Yong20 found that the strength of soil increased after freeze–thaw cycles. However, Bondarenko et al.21 and Swan et al.22 found that freeze–thaw cycles had little influence on soil strength, and the strength was basically unchanged before and after freeze–thaw cycles. At present, the different mechanical behaviors of soil under freeze–thaw action may be closely related to the test methods used as well as the type and initial state of soil mass. Due to the complex changes of soil physical properties after freeze–thaw cycles, the study on the influence of freeze–thaw cycles on soil mechanical properties and its mechanism is not enough, and no unified conclusion and theory has been formed.

The soil fixation mechanism of plant roots is mainly reflected in three aspects: the anchorage effect of deep roots, the reinforcement effect of shallow roots, and the reduction of slope pore water pressure23,24,25. It is generally believed that the tensile properties of roots have the greatest influence on the mechanical properties of soil containing roots. Some scholars have made in-depth studies on the factors affecting the tensile properties of plant roots, for example, Genet26 and Tosi et al.27 studied the tensile strength of roots of different trees and shrubs in the Mediterranean region, and showed that there were great differences in the tensile strength of roots of different tree species. The shear stress on the soil is converted into the pull on the roots through the surface friction of the roots28. In order to quantify the effect of roots on soil, Wu and Waldron29,30 believed that soil reinforcement by roots was mainly reflected in the improvement of soil shear strength, and proposed the Wu-Waldron model of soil shear strength, which considers the influence of plant roots mechanical properties. This model has few parameters and is easy to understand. Most subsequent studies are based on this model to revise and improve. Pollen31 argues that one of the assumptions of the Wu-Waldron model (simultaneous and transient fracture of all roots) can overestimate the reinforcement effect of roots on soil, and proposes the Fiber Bundle Model that considers progressive root breakage, and a more reasonable root cohesion is obtained. But this model is affected by the change of root density. Schwarz32 proposes the Root Bundle Model, assuming that the friction between the root and loam varies linearly from 0.1 to 10 kPa with the soil volume moisture content. In recent years, studies on atmosphere-vegetation-soil coupling have been developed33,34 and Ng et al.35 proposed a new model of coupled migration of groundwater seepage and surface runoff considering the influence of plant root shape.

Although many previous scholars have studied freeze–thaw cycles and rooted soil respectively, few experimental studies and theoretical models have taken these two factors into account. In order to explore this issue, the paper firstly studied the mechanical properties of rooted soil under freeze–thaw cycles, and analyzed the influences of confining pressure, freeze–thaw cycle and root system. Then based on the binary medium model36,37,38,39,40,41,42,43,44,45, the paper extends it to three-phase mixed materials, derives the extended model of binary media, and considers the influence of freeze–thaw cycle. Finally, the validity and applicability of the model are verified by experimental results.

Test material and method

Test material

The test soil samples were taken from Hailuogou, Sichuan Province, China at an altitude of about 3000 m, and the depth of the soil samples was within 0.5–1.0 m below the surface. The grain distribution curve of the test soil sample is shown in Fig. 1, and the physical properties are shown in Table 1. According to the classification standard of soil, the soil belongs to silty sand. According to the coefficient of inhomogeneity and the coefficient of curvature, the soil samples were judged to be well graded. The root system of Sorbus was selected. The length of the root was about 5 cm, the diameter of the root was about 5 mm, the weight of each root was about 1.4 g, the mass ratio was 0.292% (dry bulk density), the volume ratio was 0.262%, and the shear area ratio was 0.254%.

Figure 1
figure 1

The grain distribution curve of tested soils.

Full size image
Table 1 Physical parameters of test soil samples.
Full size table

Sample preparation

After retrieving the test soil, let it air dry naturally, and pass through a 2 mm sieve after grinding. The soil were mixed evenly at the dry density of 1.28 g/cm3, and then packed into the sample preparation cylinder and compacted in four layers, with a size of 61.8 mm in diameter and 125 mm in height. For the samples containing roots, live Sorbus roots were selected, as shown in Fig. 2, then cut out part of the root diameter for about 5 mm and the straight part for 5 cm, and placed them in the middle of the second and the third layer of soil during sample preparation. Since Sorbus roots is generally buried shallowly in the natural production process and grows horizontally, as shown in Fig. 3. Therefore, the root system of the sample was arranged horizontally and placed in the middle of the sample, as shown in Fig. 4.

Figure 2
figure 2

Sorbus root system.

Full size image
Figure 3
figure 3

Natural root state of Sorbus chinensis.

Full size image
Figure 4
figure 4

Root layout diagram.

Full size image

After the sample is prepared, it is vacuumized for 2 h and then water is released to immerse the sample in hydraulic saturation for more than 12 h. Then take out the saturated sample and cover it with rubber film and permeable stone as soon as possible. Rubber rings are attached to the positions of the permeable stone at both the upper and lower ends to prevent the sample from deforming too much in the process of freeze–thaw cycles. Then, triaxial test is carried out immediately for samples that do not need freeze–thaw cycle. For samples subjected to freeze–thaw cycles, the sample is quickly placed in water to remove air bubbles, and transferred in water to a transparent plastic tank without a lid. The transparent tank with samples was numbered and put into the automatic low-temperature freeze–thaw machine to start the freeze–thaw cycle, as shown Fig. 5. In the process of freeze–thaw cycle, the sample can be replenished by permeable stones at both ends. A freeze–thaw cycle consists of lowering room temperature to – 15 °C for 12 h and then rising to 20 °C for 12 h.

Figure 5
figure 5

Samples subjected to freeze–thaw cycles.

Full size image

Test method

After the set number of freeze–thaw cycles, the sample was taken out and installed on the triaxial tester to perform consolidated drained triaxial test at room temperature. The axial loading speed was 0.1 mm/min and the confining pressure was 25, 50, 100, 200 kPa. If the axial deformation is up to 15%, the sample will be considered damaged, the test would stop. The specific test scheme is shown in Table 2, and N in the figures and tables in this paper represents the number of freeze–thaw cycles, and r represents the number of roots contained. The soil after the test is pulverized and reconstituted again, which is called remolded soil, and the sample can be prepared under rootless and non-freeze–thaw conditions for triaxial test, which is used for the determination of some parameters.

Table 2 Test plan table.
Full size table

Ethical approval

Since this study did not recruit any human and/or animal subjects, this section does not apply.

Consent to participate

Since this study did not recruit any human subjects, this section does not apply.

Experimental analysis of mechanical properties of rooted soil under freeze–thaw cycles

Figures 6, 7, 8, 9, 10, 11 respectively describe the relationship curves of deviatoric stress, volumetric strain and axial strain under different confining pressures, different freeze–thaw cycles and different root number. \({\varepsilon }_{a}\) in the figure is for axial strain, \({\varepsilon }_{v}\) is for volumetric strain,\({\sigma }_{1}-{\sigma }_{3}\) is for deviatoric stress. It can be seen from the figure that the soil containing roots is greatly affected by confining pressure and the number of freeze–thaw cycles, and the volume ration of roots to sample has strengthening effect on the sample to some degree.

Figure 6
figure 6figure 6

Relations between axial strain and deviatoric stress under different confining pressures.

Full size image
Figure 7
figure 7figure 7

Relations between axial strain and volume strain under different confining pressures.

Full size image
Figure 8
figure 8

Relations between axial strain and deviatoric stress with different freeze–thaw cycles.

Full size image
Figure 9
figure 9

Relations between axial strain and volume strain with different freeze–thaw cycles.

Full size image
Figure 10
figure 10

Relationships between axial strain and deviatoric stress with different root number.

Full size image
Figure 11
figure 11

Relationships between axial strain and volume strain with different root number.

Full size image

Influence of confining pressure

Figures 6 and 7 show that the stress–strain curves of samples are significantly affected by confining pressure. The specimen exhibits strain softening and volumetric dilatancy under low confining pressure. When the confining pressure is 25 kPa, the deviatoric stress of the sample decreases rapidly after reaching the peak strength, and then gradually flattens out, and finally gradually approaches the residual strength. The volumetric strain of the sample expands after a short shrinkage, and the volumetric strain gradually becomes negative (that is, exceeds the initial volume of the sample), and continues to expand. At the confining pressure of 50 kPa, strain softening and volumetric dilatancy still occur, but the phenomenon is weakened. When the deviatoric stress reaches its peak, the strength decreases slowly and finally approaches the residual strength gradually. The volumetric strain of the sample is also the first bulk shrinkage, but the volume shrinkage is larger, and then the bulk expansion is slow, and the volume expansion is relatively small, and the negative bulk change will not be generated eventually. With the increase of confining pressure, the phenomenon of strain softening and swelling disappear gradually. Under high confining pressure, the mechanical properties of the samples are constant strain hardening and volume shrinkage. As the confining pressure increases, the breakage mode of the sample changes from shear breakage to swelling breakage, as shown in Fig. 12. The strength of the sample increases with the increase of confining pressure.

Figure 12
figure 12

Breakage forms of samples under different confining pressures.

Full size image

The mechanical properties of rooted soil with respect to confining pressure variation are consistent with those of typical sandy soils, as explained in the article. This is because during the consolidation process, the confining pressure damages the bonding and structural integrity between particles, and the higher the confining pressure, the greater the degree of disruption. At low confining pressures, the bonding and structural integrity between soil particles remain relatively intact during consolidation. When the specimen is sheared, the external load is mainly borne by the bonding elements, and the interlocking between particles in the weak zones of stress concentration is completely destroyed, resulting in strain softening and shear failure. At high confining pressures, the bonding and structural integrity between soil particles undergo varying degrees of damage during consolidation. Some of the bonding elements have transformed into frictional elements, and the external load is mainly borne by the frictional elements. During the shear process, the resistance to sliding and friction between particles predominantly come into play.

Influence of the number of freeze–thaw cycles

As can be seen from Figs. 8 and 9, with the increase of the number of freeze–thaw cycles, both stress and volumetric strain curves are located below, which means that the deviatoric stress decreases while the volumetric strain increases. At low confining pressure, with the increase of freeze–thaw cycles, the phenomenon of strain softening and volume expansion will be weaker, the deviatoric stress will decrease more gently after reaching the peak, the volume shrinkage will be larger, and the volume dilatancy will be more gentle. At the confining pressure of 25 kPa, the volume change will not reach negative value after the freeze–thaw cycles have reached 5 times. Under high confining pressure, with the increase of freeze–thaw cycles, the bearing capacity of the sample decreases, which is manifested as the deviatoric stress decreases gradually and the volume shrinkage value increases gradually. In general, with the increase of the number of freeze–thaw cycles, the mechanical properties and strength of the samples decreased gradually.

This is because the soil before freeze–thaw is structural and there is cementation among the particles, while in the freeze–thaw cycle, the soil particles bear the frost heaving force. In the freeze–thaw cycle, this bonding between the particles is destroyed, and the sample changes from a relatively dense state to a relatively loose state (‘tight sand’ becomes ‘loose sand’). In the original soil, there existed structural interlocking and bonding between soil particles. During the freezing process, the original bonding structure is disrupted due to the significant frost heave forces acting on the soil particles. However, the generated ice crystals bond the soil particles together, providing new and stronger bonding effects. This is why frozen soil often exhibits increased strength. During the thawing process, the ice crystals transition from a solid to a liquid state, causing the bonding effects of the ice crystals to disappear. After a complete freeze–thaw cycle, a portion of the soil particles undergo damage, transitioning from bonding elements to frictional elements, leading to a decrease in the mechanical properties and bearing capacity of the samples. As the number of freeze–thaw cycles increases, the proportion and degree of soil particle damage increase, resulting in a greater reduction in the mechanical properties and bearing capacity of the samples.

Effect of volume ration of roots to sample

As can be seen from Figs. 10 and 11, with the increase of volume ration of roots to sample, the deviatoric stress increases, the deformation decreases, and the strength gradually increases. On one hand, According to the homogenization theory of composite materials, plant roots have a certain strengthening effect on the mechanical properties of the samples. The modulus and strength of plant roots in the samples are higher than that of soil, and the bearing capacity of the samples with roots is slightly higher than that without roots. On the other hand, the combination of roots and soil allows them to leverage their respective advantages, with roots providing tensile strength and soil providing compressive strength. When subjected to external loads, both the soil and roots undergo deformation. Due to the significant difference in their deformation moduli, there will be relative displacement or a tendency for relative displacement between the soil and roots. Stress within the soil redistributes, resulting in surface frictional forces between the roots and soil particles. This leads to tensile forces in the roots, which, in turn, hinder soil deformation through frictional resistance, thereby enhancing the overall strength of the soil. But the rising level is not high, because the root system occupies a small proportion of the volume. In actual soil mass, the strength of soil mass is improved to some extent due to the reinforcement effect of shallow roots and the anchorage effect of deep roots, which can be used as a safety reserve for design.

Strength

When the stress–strain curve is strain softening, the peak value is taken as the shear strength of the sample. When the stress–strain curve is continuous hardening, the corresponding deviatoric stress of and the axial strain is 15% is the shear strength of the sample. The shear strength of all samples is summarized in Tables 3, 4, 5. It can be seen from the table that the strength of the sample increases with the increase of confining pressure (it is also obtained by Tian et al.46), decreases with the increase of freeze–thaw cycles, and increases with the increase of the number of roots contained. Comparing specimens with 1, 5, and 15 freeze–thaw cycles to samples without freeze–thaw cycles, the average reduction in strength is 4.66%, 9.08%, and 12.11%, respectively. It is because (1) with the increase of confining pressure, the soil becomes more compact and the shear modulus is larger; (2) In the process of freeze–thaw cycle, the soil structure is gradually destroyed; and (3) the root has reinforcing effect on the soil.

Table 3 Summary of shear strength of samples without root (unit: kPa).
Full size table
Table 4 Summary of shear strength of samples with one root (unit: kPa).
Full size table
Table 5 Summary of shear strength of samples with three roots (unit: kPa).
Full size table

Extended binary medium model of rooted soil

Formulation of the constitutive model

The traditional binary medium model abstracts the soil as a mixture composed of structural blocks and weak zones (called bonded elements and frictional elements, respectively). The load is shared by the bonded elements and frictional elements that will be damaged and transformed. During the loading process, the bonded element is gradually destroyed and transformed into frictional element. Rooted soil is regarded as a heterogeneous material composed of soil and roots. Soil can be regarded as the component of damaged transformed bonded element and friction element, while root can be regarded as another similar bonded elements material that will not be damaged. It is assumed that the strength and modulus of the root are large enough to prevent the root from being damaged during the loading process. The schematic diagram of the model is shown in Fig. 13.

Figure 13
figure 13

Schematic diagram of the model.

Full size image

Take a representative volume element (RVE), and the stress and strain on the bonded element, frictional element and root are respectively:

$${\sigma }_{ij}^{b}=\frac{1}{{V}^{b}}\int {\sigma }_{ij}^{loc}d{V}^{b}$$
(1)
$${\varepsilon }_{ij}^{b}=\frac{1}{{V}^{b}}\int {\varepsilon }_{ij}^{loc}d{V}^{b}$$
(2)
$${\sigma }_{ij}^{f}=\frac{1}{{V}^{f}}\int {\sigma }_{ij}^{loc}d{V}^{f}$$
(3)
$${\varepsilon }_{ij}^{f}=\frac{1}{{V}^{f}}\int {\varepsilon }_{ij}^{loc}d{V}^{f}$$
(4)
$${\sigma }_{ij}^{r}=\frac{1}{{V}^{r}}\int {\sigma }_{ij}^{loc}d{V}^{r}$$
(5)
$${\varepsilon }_{ij}^{r}=\frac{1}{{V}^{r}}\int {\varepsilon }_{ij}^{loc}d{V}^{r}$$
(6)

where \({\sigma }_{ij}^{loc}\) and \({\varepsilon }_{ij}^{loc}\) are the local microscopic stress and strain in the RVE. For the RVE, we have

$$V={V}^{b}+{V}^{f}+{V}^{r}$$
(7)

Based on the homogenization theory of heterogeneous media, the following formula can be derived:

$${\sigma }_{ij}=\frac{1}{V}\int {\sigma }_{ij}^{loc}d\mathrm{V}$$
(8)
$${\varepsilon }_{ij}=\frac{1}{V}\int {\varepsilon }_{ij}^{loc}d\mathrm{V}$$
(9)

Substituting Eqs. (1)–(7) into (8) and (9), we get:

$${\sigma }_{ij}=\frac{{V}^{b}}{V}{\sigma }_{ij}^{b}+\frac{{V}^{f}}{V}{\sigma }_{ij}^{f}+\frac{{V}^{r}}{V}{\sigma }_{ij}^{r}$$
(10)
$${\varepsilon }_{ij}=\frac{{V}^{b}}{V}{\varepsilon }_{ij}^{b}+\frac{{V}^{f}}{V}{\varepsilon }_{ij}^{f}+\frac{{V}^{r}}{V}{\varepsilon }_{ij}^{r}$$
(11)

Definition \({\lambda }_{v}^{f}=\frac{{V}^{f}}{V}\) represents volume breakage ratio; \({\lambda }_{v}^{r}=\frac{{V}^{r}}{V}\) represents the volume ration of roots to sample. λ represents the stress sharing ratio, and Eqs. (10) and (11) are transformed into:

$${\sigma }_{ij}=\left(1-{\lambda }^{r}-{\lambda }^{f}\right){\sigma }_{ij}^{b}+{\lambda }^{f}{\sigma }_{ij}^{f}+{\lambda }^{r}{\sigma }_{ij}^{r}$$
(12)
$${\varepsilon }_{ij}=\left(1-{\lambda }^{r}-{\lambda }^{f}\right){\varepsilon }_{ij}^{b}+{\lambda }^{f}{\varepsilon }_{ij}^{f}+{\lambda }^{r}{\varepsilon }_{ij}^{r}$$
(13)

Even if the stress–strain relationship of the bonded element and root system adopts linear elastic constitutive relationship, the stress–strain relationship of the frictional element is still nonlinear elastic–plastic. Therefore, the previous formula cannot simulate the actual stress–strain relationship, so it is necessary to use the linear method by replacing the elastic method to conduct full differentiation of Eqs. (12) and (13), and change them into incremental expressions, as follows:

$$d{\sigma }_{ij}=\left(1-{\lambda }^{r}-{\lambda }^{f}\right){d\sigma }_{ij}^{b}+{\lambda }^{f}d{\sigma }_{ij}^{f}+{\lambda }^{r}{d\sigma }_{ij}^{r}+d{\lambda }^{f}\left({\sigma }_{ij}^{f0}-{\sigma }_{ij}^{b0}\right)+d{\lambda }^{r}\left({\sigma }_{ij}^{r0}-{\sigma }_{ij}^{b0}\right)$$
(14)
$${d\varepsilon }_{ij}=\left(1-{\lambda }^{r}-{\lambda }^{f}\right)d{\varepsilon }_{ij}^{b}+{\lambda }^{f}d{\varepsilon }_{ij}^{f}+{\lambda }^{r}d{\varepsilon }_{ij}^{r}+d{\lambda }^{f}\left({\varepsilon }_{ij}^{f0}-{\varepsilon }_{ij}^{b0}\right)+d{\lambda }^{r}\left({\varepsilon }_{ij}^{r0}-{\varepsilon }_{ij}^{b0}\right)$$
(15)

It can be seen from Eqs. (14) and (15) that the average stress–strain increment of the representative element is composed of four parts, including the influence of bonded element, frictional element, root and breakage ratio.

Two-parameter binary medium model

The stress is decomposed into the spherical stress and deviatoric stress, i.e. mean stress and generalized shear stress. The strain is decomposed into the volume strain and shear strain.

$${\sigma }_{m}=\frac{1}{3}{\sigma }_{kk}$$
(16)
$${\sigma }_{s}=\sqrt{\frac{3}{2}\left({\sigma }_{ij}-\frac{1}{3}{\delta }_{ij}{\sigma }_{kk}\right)\left({\sigma }_{ij}-\frac{1}{3}{\delta }_{ij}{\sigma }_{kk}\right)}$$
(17)
$${\varepsilon }_{v}={\varepsilon }_{kk}$$
(18)
$${\varepsilon }_{s}=\sqrt{\frac{2}{3}\left({\varepsilon }_{ij}-\frac{1}{3}{\delta }_{ij}{\varepsilon }_{kk}\right)\left({\varepsilon }_{ij}-\frac{1}{3}{\delta }_{ij}{\varepsilon }_{kk}\right)}$$
(19)

Since it is assumed that the root will not be destroyed in the loading process and \({\lambda }^{r}\) is considered to be constant, and \({d\lambda }^{r}=0\), then:

$${d\sigma }_{m}=\left(1-{\lambda }_{v}^{r0}-{\lambda }_{v}^{f0}\right)d{\sigma }_{m}^{b}+{\lambda }_{v}^{f0}{d\sigma }_{m}^{f}+{\lambda }_{v}^{r0}d{\sigma }_{m}^{r}+d{\lambda }_{v}^{f}\left({\sigma }_{m}^{f0}-{\sigma }_{m}^{b0}\right)$$
(20)
$${d\varepsilon }_{v}=\left(1-{\lambda }_{v}^{r0}-{\lambda }_{v}^{f0}\right){d\varepsilon }_{v}^{b}+{\lambda }_{v}^{f0}{d\varepsilon }_{v}^{f}+{\lambda }_{v}^{r0}d{\varepsilon }_{v}^{r}+d{\lambda }_{v}^{f}\left({\varepsilon }_{v}^{f0}-{\varepsilon }_{v}^{b0}\right)$$
(21)
$${d\sigma }_{s}=\left(1-{\lambda }_{s}^{r0}-{\lambda }_{s}^{f0}\right)d{\sigma }_{s}^{b}+{\lambda }_{s}^{f0}{d\sigma }_{s}^{f}+{\lambda }_{s}^{r0}d{\sigma }_{s}^{r}+d{\lambda }_{s}^{f}\left({\sigma }_{s}^{f0}-{\sigma }_{s}^{b0}\right)$$
(22)
$$d{\varepsilon }_{s}=\left(1-{\lambda }_{s}^{r0}-{\lambda }_{s}^{f0}\right){d\varepsilon }_{s}^{b}+{\lambda }_{s}^{f0}{d\varepsilon }_{s}^{f}+{\lambda }_{s}^{r0}d{\varepsilon }_{s}^{r}+d{\lambda }_{s}^{f}\left({\varepsilon }_{s}^{f0}-{\varepsilon }_{s}^{b0}\right)$$
(23)

Similarly, \({\lambda }_{s}^{f}\) is defined as the area breakage ratio during cutting process, and \({\lambda }_{s}^{r}\) is the area ratio of roots to sample.

From Eqs. (20) to (23), we can get:

$${d\sigma }_{m}^{f}=\frac{1}{{\lambda }_{v}^{f0}}\left[{d\sigma }_{m}-\left(1-{\lambda }_{v}^{r0}-{\lambda }_{v}^{f0}\right)d{\sigma }_{m}^{b}-{\lambda }_{v}^{r0}d{\sigma }_{m}^{r}-d{\lambda }_{v}^{f}\left({\sigma }_{m}^{f0}-{\sigma }_{m}^{b0}\right)\right]$$
(24)
$${d\sigma }_{s}^{f}=\frac{1}{{\lambda }_{s}^{f0}}\left[{d\sigma }_{s}-\left(1-{\lambda }_{s}^{r0}-{\lambda }_{s}^{f0}\right)d{\sigma }_{s}^{b}-{\lambda }_{s}^{r0}d{\sigma }_{s}^{r}-d{\lambda }_{s}^{f}\left({\sigma }_{s}^{f0}-{\sigma }_{s}^{b0}\right)\right]$$
(25)
$${d\varepsilon }_{v}^{f}=\frac{1}{{\lambda }_{v}^{f0}}\left[{d\varepsilon }_{v}-\left(1-{\lambda }_{v}^{r0}-{\lambda }_{v}^{f0}\right){d\varepsilon }_{v}^{b}-{\lambda }_{v}^{r0}d{\varepsilon }_{v}^{r}-d{\lambda }_{v}^{f}\left({\varepsilon }_{v}^{f0}-{\varepsilon }_{v}^{b0}\right)\right]$$
(26)
$${d\varepsilon }_{s}^{f}=\frac{1}{{\lambda }_{s}^{f0}}\left[d{\varepsilon }_{s}-\left(1-{\lambda }_{s}^{r0}-{\lambda }_{s}^{f0}\right){d\varepsilon }_{s}^{b}-{\lambda }_{s}^{r0}d{\varepsilon }_{s}^{r}-d{\lambda }_{s}^{f}\left({\varepsilon }_{s}^{f0}-{\varepsilon }_{s}^{b0}\right)\right]$$
(27)

The constitutive relation of the bonded elements and root is assumed to be linear elasticity:

$${d\varepsilon }_{v}^{b}=\frac{1}{{K}^{b}}d{\sigma }_{m}^{b}$$
(28)
$${d\varepsilon }_{s}^{b}=\frac{1}{3{G}^{b}}d{\sigma }_{s}^{b}$$
(29)
$${d\varepsilon }_{v}^{r}=\frac{1}{{K}^{r}}d{\sigma }_{m}^{r}$$
(30)
$${d\varepsilon }_{s}^{r}=\frac{1}{3{G}^{r}}d{\sigma }_{s}^{r}$$
(31)

The frictional element adopts the nonlinear elastic–plastic constitutive, and the volume strain and shear strain are the coupling relations between the spherical stress and deviator stress:

$${d\varepsilon }_{v}^{f}={K}_{v}^{f}{d\sigma }_{m}^{f}+{K}_{s}^{f}{d\sigma }_{s}^{f}$$
(32)
$${d\varepsilon }_{s}^{f}={G}_{v}^{f}{d\sigma }_{m}^{f}+{G}_{s}^{f}{d\sigma }_{s}^{f}$$
(33)

Due to the non-uniformity of deformation of cement element, frictional elements and root in the RVE, local stress concentration coefficient is introduced to represent the relationship between local stress of bonded element, root and average stress of the REV.

$${\sigma }_{m}^{b}={C}_{m}^{b}{\sigma }_{m}$$
(34)
$${\sigma }_{s}^{b}={C}_{s}^{b}{\sigma }_{s}$$
(35)
$${\sigma }_{m}^{r}={C}_{m}^{r}{\sigma }_{m}$$
(36)
$${\sigma }_{s}^{r}={C}_{s}^{r}{\sigma }_{s}$$
(37)

The incremental forms of Eqs. (34) to (37) are:

$${d\sigma }_{m}^{b}={{\sigma }_{m}^{0}d{C}_{m}^{b}+C}_{m}^{b0}d{\sigma }_{m}={B}_{m}^{b}d{\sigma }_{m}$$
(38)
$${d\sigma }_{s}^{b}={{\sigma }_{s}^{0}d{C}_{S}^{b}+C}_{s}^{b0}d{\sigma }_{s}={B}_{S}^{b}d{\sigma }_{s}$$
(39)
$${d\sigma }_{m}^{r}={{\sigma }_{m}^{0}d{C}_{m}^{r}+C}_{m}^{r0}d{\sigma }_{m}={B}_{m}^{r}d{\sigma }_{m}$$
(40)
$${d\sigma }_{s}^{r}={{\sigma }_{s}^{0}d{C}_{S}^{r}+C}_{s}^{r0}d{\sigma }_{s}={B}_{s}^{r}d{\sigma }_{s}$$
(41)
$${B}_{m}^{b}={C}_{m}^{b0}+{\frac{\partial {C}_{m}^{b}}{\partial {\sigma }_{m}}}\sigma_{m}^{0}$$
(42)
$${B}_{s}^{b}={C}_{s}^{b0}+{\frac{\partial {C}_{s}^{b}}{\partial {\sigma }_{s}}}\sigma_{s}^{0}$$
(43)
$${B}_{m}^{r}={C}_{m}^{r0}+{\frac{\partial {C}_{m}^{r}}{\partial {\sigma }_{m}}}\sigma_{m}^{0}$$
(44)
$${B}_{s}^{b}={C}_{s}^{b0}+{\frac{\partial {C}_{s}^{r}}{\partial {\sigma }_{s}}}\sigma_{s}^{0}$$
(45)

Then, two parameters, volume breakage ratio and area breakage ratio, were introduced to describe the breakage characteristics of bonded elements and frictional elements. The expression is:

$${\lambda }_{v}^{f}={f}_{1}({\sigma }_{m},N)$$
(46)
$${\lambda }_{s}^{f}={f}_{2}({\sigma }_{s},N)$$
(47)

where \({f}_{1}\) is the function of the spherical stress and the number of freeze–thaw cycles, \({f}_{2}\) is the function of deviatoric stress and the number of freeze–thaw cycles.

The total differentiation of Eqs. (46) and (47) can be obtained as follows:

$${d\lambda }_{v}^{f}=\frac{\partial {f}_{1}}{\partial {\sigma }_{m}}{d\sigma }_{m}+\frac{\partial {f}_{1}}{\partial N}dN$$
(48)
$${d\lambda }_{s}^{f}=\frac{\partial {f}_{2}}{\partial {\sigma }_{s}}{d\sigma }_{s}+\frac{\partial {f}_{2}}{\partial N}dN$$
(49)

According to Eqs. (12), (13), (34–37), we can get:

$${\sigma }_{m}^{f0}-{\sigma }_{m}^{b0}=\frac{\left(1-{C}_{m}^{b0}\right)}{{\lambda }_{v}^{f0}}{\sigma }_{m}^{0}+\frac{{\lambda }_{v}^{r0}}{{\lambda }_{v}^{f0}}\left({C}_{m}^{b0}-{C}_{m}^{r0}\right){\sigma }_{m}^{0}$$
(50)
$${\sigma }_{s}^{f0}-{\sigma }_{s}^{b0}=\frac{\left(1-{C}_{s}^{b0}\right)}{{\lambda }_{s}^{f0}}{\sigma }_{s}^{0}+\frac{{\lambda }_{s}^{r0}}{{\lambda }_{s}^{f0}}\left({C}_{s}^{b0}-{C}_{s}^{r0}\right){\sigma }_{s}^{0}$$
(51)
$${\varepsilon }_{v}^{f0}-{\varepsilon }_{v}^{b0}=\frac{1}{{\lambda }_{v}^{f0}}\left[\left({\varepsilon }_{v}^{0}{-\varepsilon }_{v}^{b0}\right)+{\lambda }_{v}^{r0}\left({\varepsilon }_{v}^{b0}-{\varepsilon }_{v}^{r0}\right)\right]$$
(52)
$${\varepsilon }_{s}^{f0}-{\varepsilon }_{s}^{b0}=\frac{1}{{\lambda }_{s}^{f0}}\left[\left({\varepsilon }_{s}^{0}{-\varepsilon }_{s}^{b0}\right)+{\lambda }_{s}^{r0}\left({\varepsilon }_{s}^{b0}-{\varepsilon }_{s}^{r0}\right)\right]$$
(53)

By combining Eqs. (20)–(49), we can finally get:

$${d\varepsilon }_{v}={A}_{1}d{\sigma }_{m}+{B}_{1}d{\sigma }_{s}+{C}_{1}d\mathrm{N}$$
(54)
$${d\varepsilon }_{s}={A}_{2}d{\sigma }_{m}+{B}_{2}d{\sigma }_{s}+{C}_{2}d\mathrm{N}$$
(55)

Among them:

$${A}_{1} =\left[{K}_{v}^{f}+\left(1-{\lambda }_{v}^{r0}-{\lambda }_{v}^{f0}\right)\left(\frac{1}{{K}^{b}}-{K}_{v}^{f}\right){B}_{m}^{b} +{\lambda }_{v}^{r0}\left(\frac{1}{{K}^{r}}-{K}_{v}^{f}\right){B}_{m}^{r}\right]+\frac{1}{{\lambda }_{v}^{f0}}\left\{\left[\left({\varepsilon }_{v}^{0}{-\varepsilon }_{v}^{b0}\right)+{\lambda }_{v}^{r0}\left({\varepsilon }_{v}^{b0}-{\varepsilon }_{v}^{r0}\right)\right]-{K}_{v}^{f}\left[\left(1-{C}_{m}^{b0}\right){\sigma }_{m}^{0}+{\lambda }_{v}^{r0}\left({C}_{m}^{b0}-{C}_{m}^{r0}\right){\sigma }_{m}^{0}\right]\right\}\frac{\partial {f}_{1}}{\partial {\sigma }_{m}}$$
(56)
$${B}_{1}=\frac{{\lambda }_{v}^{f0}}{{\lambda }_{s}^{f0}}{K}_{s}^{f}\left\{1-\left(1-{\lambda }_{s}^{r0}-{\lambda }_{s}^{f0}\right){B}_{S}^{b}-{\lambda }_{s}^{r0}{B}_{s}^{r}-\left[\frac{\left(1-{C}_{s}^{b0}\right)}{{\lambda }_{s}^{f0}}{\sigma }_{s}^{0}+\frac{{\lambda }_{s}^{r0}}{{\lambda }_{s}^{f0}}\left({C}_{s}^{b0}-{C}_{s}^{r0}\right){\sigma }_{s}^{0}\right]\frac{\partial {f}_{2}}{\partial {\sigma }_{s}}\right\}$$
(57)
$${C}_{1}=\frac{1}{{\lambda }_{v}^{f0}}\left\{\left[\left({\varepsilon }_{v}^{0}{-\varepsilon }_{v}^{b0}\right)+{\lambda }_{v}^{r0}\left({\varepsilon }_{v}^{b0}-{\varepsilon }_{v}^{r0}\right)\right]-{K}_{v}^{f}\left[\left(1-{C}_{m}^{b0}\right){\sigma }_{m}^{0}+{\lambda }_{v}^{r0}\left({C}_{m}^{b0}-{C}_{m}^{r0}\right){\sigma }_{m}^{0}\right]\right\}\frac{\partial {f}_{1}}{\partial N}-\frac{{\lambda }_{v}^{f0}}{{\left({\lambda }_{s}^{f0}\right)}^{2}}{K}_{s}^{f}\left[{\left(1-{C}_{s}^{b0}\right)}\sigma_{s}^{0}+{\lambda }_{s}^{r0}\left({C}_{s}^{b0}-{C}_{s}^{r0}\right){\sigma }_{s}^{0}\right]\frac{\partial {f}_{2}}{\partial N}$$
(58)
$${A}_{2}=\frac{{\lambda }_{s}^{f0}}{{\lambda }_{v}^{f0}}{G}_{v}^{f}\left[1-\left(1-{\lambda }_{v}^{r0}-{\lambda }_{v}^{f0}\right){B}_{m}^{b}-{\lambda }_{v}^{r0}{B}_{m}^{r}\right]-\frac{{\lambda }_{s}^{f0}}{{\left({\lambda }_{v}^{f0}\right)}^{2}}{G}_{v}^{f}\left[\left(1-{C}_{m}^{b0}\right){\sigma }_{m}^{0}+{\lambda }_{v}^{r0}\left({C}_{m}^{b0}-{C}_{m}^{r0}\right){\sigma }_{m}^{0}\right]\frac{\partial {f}_{1}}{\partial {\sigma }_{m}}$$
(59)
$${B}_{2}=\left[{G}_{s}^{f}+\left(1-{\lambda }_{s}^{r0}-{\lambda }_{s}^{f0}\right)\left(\frac{1}{3{G}^{b}}-{G}_{s}^{f}\right) {B}_{S}^{b}+{\lambda }_{s}^{r0}\left(\frac{1}{3{G}^{r}}-{G}_{s}^{f}\right){B}_{s}^{r}\right]+\frac{1}{{\lambda }_{s}^{f0}}\left\{\left({\varepsilon }_{s}^{0}{-\varepsilon }_{s}^{b0}\right)+{\lambda }_{s}^{r0}\left({\varepsilon }_{s}^{b0}-{\varepsilon }_{s}^{r0}\right)-\left[\left(1-{C}_{s}^{b0}\right){\sigma }_{s}^{0}+{\lambda }_{s}^{r0}\left({C}_{s}^{b0}-{C}_{s}^{r0}\right){\sigma }_{s}^{0}\right]{G}_{s}^{f}\right\}\frac{\partial {f}_{2}}{\partial {\sigma }_{s}}$$
(60)
$${C}_{2}=\frac{1}{{\lambda }_{s}^{f0}}\left\{\left({\varepsilon }_{s}^{0}{-\varepsilon }_{s}^{b0}\right)+{\lambda }_{s}^{r0}\left({\varepsilon }_{s}^{b0}-{\varepsilon }_{s}^{r0}\right)-\left[\left(1-{C}_{s}^{b0}\right){\sigma }_{s}^{0}+{\lambda }_{s}^{r0}\left({C}_{s}^{b0}-{C}_{s}^{r0}\right){\sigma }_{s}^{0}\right]{G}_{s}^{f}\right\}\frac{\partial {f}_{2}}{\partial N}-\frac{{\lambda }_{s}^{f0}}{{\left({\lambda }_{v}^{f0}\right)}^{2}}{G}_{v}^{f}\left[\left(1-{C}_{m}^{b0}\right){\sigma }_{m}^{0}+{\lambda }_{v}^{r0}\left({C}_{m}^{b0}-{C}_{m}^{r0}\right){\sigma }_{m}^{0}\right]\frac{\partial {f}_{1}}{\partial N}$$
(61)

When N = 0, there is no dN and it becomes a binary medium model containing root soil.

When \({\lambda }^{r}={\lambda }_{v}^{r}={\lambda }_{s}^{r}=0\), it becomes a rootless freeze–thaw cycle binary medium model, and the coefficient is taken as \({A}_{1}^{{{\prime}}}\), \({B}_{1}^{{{\prime}}}\), \({C}_{1}^{{{\prime}}}\), \({A}_{2}^{{{\prime}}}\), \({B}_{2}^{{{\prime}}}\), \({C}_{2}^{{{\prime}}}\).

$${A}_{1}^{{{\prime}}}={K}_{v}^{f}+\left(1-{\lambda }_{v}^{f0}\right)\left(\frac{1}{{K}^{b}}-{K}_{v}^{f}\right){B}_{m}^{b}+\frac{1}{{\lambda }_{v}^{f0}}\left[\left({\varepsilon }_{v}^{0}{-\varepsilon }_{v}^{b0}\right)-{K}_{v}^{f}\left(1-{C}_{m}^{b0}\right){\sigma }_{m}^{0}\right]\frac{\partial {f}_{1}}{\partial {\sigma }_{m}}$$
(62)
$${B}_{1}^{{{\prime}}}=\frac{{\lambda }_{v}^{f0}}{{\lambda }_{s}^{f0}}{K}_{s}^{f}\left[1-\left(1-{\lambda }_{s}^{f0}\right){B}_{S}^{b}\right]-\frac{{\lambda }_{v}^{f0}}{{\left({\lambda }_{s}^{f0}\right)}^{2}}{K}_{s}^{f}\left(1-{C}_{s}^{b0}\right){\sigma }_{s}^{0}\frac{\partial {f}_{2}}{\partial {\sigma }_{s}}$$
(63)
$${C}_{1}^{{{\prime}}}=\frac{1}{{\lambda }_{v}^{f0}}\left[\left({\varepsilon }_{v}^{0}{-\varepsilon }_{v}^{b0}\right)-{K}_{v}^{f}\left(1-{C}_{m}^{b0}\right){\sigma }_{m}^{0}\right]\frac{\partial {f}_{1}}{\partial N}-\frac{{\lambda }_{v}^{f0}}{{\left({\lambda }_{s}^{f0}\right)}^{2}}{K}_{s}^{f}\left(1-{C}_{s}^{b0}\right){\sigma }_{s}^{0}\frac{\partial {f}_{2}}{\partial N}$$
(64)
$${A}_{2}^{{{\prime}}}=\frac{{\lambda }_{s}^{f0}}{{\lambda }_{v}^{f0}}{G}_{v}^{f}\left[1-\left(1-{\lambda }_{v}^{f0}\right){B}_{m}^{b}\right]-\frac{{\lambda }_{s}^{f0}}{{\left({\lambda }_{v}^{f0}\right)}^{2}}{G}_{v}^{f}\left(1-{C}_{m}^{b0}\right){\sigma }_{m}^{0}\frac{\partial {f}_{1}}{\partial {\sigma }_{m}}$$
(65)
$${B}_{2}^{{{\prime}}}={G}_{s}^{f}+\left(1-{\lambda }_{s}^{f0}\right)\left(\frac{1}{3{G}^{b}}-{G}_{s}^{f}\right) {B}_{S}^{b}+\frac{1}{{\lambda }_{s}^{f0}}\left[\left({\varepsilon }_{s}^{0}{-\varepsilon }_{s}^{b0}\right)-\left(1-{C}_{s}^{b0}\right){{\sigma }_{s}^{0}G}_{s}^{f}\right]\frac{\partial {f}_{2}}{\partial {\sigma }_{s}}$$
(66)
$${C}_{2}^{{{\prime}}}=\frac{1}{{\lambda }_{s}^{f0}}\left[\left({\varepsilon }_{s}^{0}{-\varepsilon }_{s}^{b0}\right)-\left(1-{C}_{s}^{b0}\right){{\sigma }_{s}^{0}G}_{s}^{f}\right]\frac{\partial {f}_{2}}{\partial N}-\frac{{\lambda }_{s}^{f0}}{{\left({\lambda }_{v}^{f0}\right)}^{2}}{G}_{v}^{f}\left(1-{C}_{m}^{b0}\right){\sigma }_{m}^{0}\frac{\partial {f}_{1}}{\partial N}$$
(67)

When \({\lambda }^{r}={\lambda }_{v}^{r}={\lambda }_{s}^{r}=0\), then N = 0, it is simplified to a general two-parameter binary medium model without dN term, and the coefficient is taken as \({A}_{1}^{{{\prime}}}\), \({B}_{1}^{{{\prime}}}\), \({A}_{2}^{{{\prime}}}\), \({B}_{2}^{{{\prime}}}\).

Model parameter determination

Constitutive relation of the bonded elements

It is assumed that the bonded elements is an ideal elastic-brittle material. Before the stress state reaches the breakage level, it is characteristic of the linear elastic constitutive model. As followed by Hooke’s law, it is presented in Eqs. (28) and (29). \({K}^{b} \; \mathrm{and } \; {G}^{b}\) can be tested by in-situ test or laboratory test of geotechnical materials. For example, \({K}^{b }\mathrm{and }{G}^{b}\) can be identified by the initial slope of curves \({\varepsilon }_{v}-{\sigma }_{m}\) and \({{\varepsilon }_{s}-\sigma }_{s}\) of conventional consolidated drained triaxial test. Or the elastic modulus E is approximated by the deformation modulus of 0.2% axial strain, and then from \(K=\frac{E}{3(1-2\nu )}\) and \(G=\frac{E}{2(1+\nu )}\), it can derive that \(\nu\) is the Poisso’s ratio.

Constitutive relation of the frictional elements

When the external stress of the bonded elements reaches the breakage level, the cementing bond breaks and the bonded elements transforms into the frictional elements to continue to bear the load. Frictional elements can be considered as isotropic elastoplastic materials. In order to better simulate the strain-softening and volumetric dilatancy phenomena of samples, the double-hardening constitutive model with wide applicability, which is proposed by Shen, is adopted in this paper47,48.

According to classical elastoplastic theory, the total strain is composed of elastic strain and plastic strain. Hooke’s law is adopted for elastic strain and plastic theory is adopted for plastic strain. The strain increment of frictional elements can be expressed as:

$${d\varepsilon }_{v}^{f}={d\varepsilon }_{v}^{fe}+{d\varepsilon }_{v}^{fp}$$
(68)
$${d\varepsilon }_{s}^{f}={d\varepsilon }_{s}^{fe}+{d\varepsilon }_{s}^{fp}$$
(69)

The elastic stress–strain relationship of frictional elements is:

$${d\varepsilon }_{v}^{fe}=\frac{{d\sigma }_{m}^{f}}{{K}^{f}}$$
(70)
$${d\varepsilon }_{s}^{fe}=\frac{{d\sigma }_{s}^{f}}{{3G}^{f}}$$
(71)

The double-hardening constitutive model includes two hardening parameters, plastic volumetric strain and plastic shear strain, and its yield function is:

$${f}^{f}=\frac{{\sigma }_{m}^{f}}{1-{({\eta }^{f}/\alpha )}^{n}}-H$$
(72)
$$\alpha ={\alpha }_{m}\left[1.0-{c}_{1}exp\left({\varepsilon }_{s}^{fp}/{c}_{2}\right)\right]$$
(73)
$$H={p}_{0}exp\left(\beta {\varepsilon }_{v}^{fp}\right)$$
(74)

where:

$${\eta }^{f}=\frac{{\sigma }_{s}^{f}}{{\sigma }_{m}^{f}}$$
(75)
$$\beta =\frac{1+{e}_{0}}{\lambda {\prime}-{\rm K}}$$
(76)

In order to better reflect the deformation characteristics of frictional elements, the non-associated flow law is adopted, that is, the yield function and the plastic potential function are the same. The plastic potential function is:

$${g}^{f}=\frac{{\sigma }_{m}^{f}}{1-{({\eta }^{f}/\alpha )}^{{n}_{1}}}-H$$
(77)

According to the flow rule:

$${d\varepsilon }_{v}^{fp}=d\Lambda \frac{\partial {g}^{f}}{\partial {\sigma }_{m}^{f}}$$
(78)
$${d\varepsilon }_{s}^{fp}=d\Lambda \frac{\partial {g}^{f}}{\partial {\sigma }_{s}^{f}}$$
(79)

in which \(d\Lambda >0\).

The consistency condition is applied Eq. (72), i.e.

$$\frac{\partial {f}^{f}}{\partial {\sigma }_{m}^{f}}{d\sigma }_{m}^{f}+\frac{\partial {f}^{f}}{\partial {\sigma }_{s}^{f}}{d\sigma }_{s}^{f}+\frac{\partial {f}^{f}}{\partial \alpha }\frac{\partial \alpha }{\partial {\varepsilon }_{s}^{fp}}d{\varepsilon }_{s}^{fp}+\frac{\partial {f}^{f}}{\partial H}\frac{\partial H}{\partial {\varepsilon }_{v}^{fp}}{d\varepsilon }_{v}^{fp}=0$$
(80)

We get:

$$d\Lambda =\frac{1}{h}\left(\frac{\partial {f}^{f}}{\partial {\sigma }_{m}^{f}}{d\sigma }_{m}^{f}+\frac{\partial {f}^{f}}{\partial {\sigma }_{s}^{f}}{d\sigma }_{s}^{f}\right)$$
(81)
$$h=-\frac{\partial {f}^{f}}{\partial H}\frac{\partial H}{\partial {\varepsilon }_{v}^{fp}}\frac{\partial {g}^{f}}{\partial {\sigma }_{m}^{f}}-\frac{\partial {f}^{f}}{\partial \alpha }\frac{\partial \alpha }{\partial {\varepsilon }_{s}^{fp}}\frac{\partial {g}^{f}}{\partial {\sigma }_{s}^{f}}$$
(82)

In which:

$$\frac{\partial {f}^{f}}{\partial {\sigma }_{m}^{f}}=\frac{1-(1+n){({\eta }^{f}/\alpha )}^{n}}{{\left[1-{({\eta }^{f}/\alpha )}^{n}\right]}^{2}}$$
(83)
$$\frac{\partial {g}^{f}}{\partial {\sigma }_{m}^{f}}=\frac{1-(1+{n}_{1}){({\eta }^{f}/\alpha )}^{{n}_{1}}}{{\left[1-{({\eta }^{f}/\alpha )}^{{n}_{1}}\right]}^{2}}$$
(84)
$$\frac{\partial {f}^{f}}{\partial {\sigma }_{s}^{f}}=\frac{n{({\eta }^{f}/\alpha )}^{n-1}}{{\alpha \left[1-{({\eta }^{f}/\alpha )}^{n}\right]}^{2}}$$
(85)
$$\frac{\partial {g}^{f}}{\partial {\sigma }_{s}^{f}}=\frac{{n}_{1}{({\eta }^{f}/\alpha )}^{{n}_{1}-1}}{{\alpha \left[1-{({\eta }^{f}/\alpha )}^{{n}_{1}}\right]}^{2}}$$
(86)
$$\frac{\partial {f}^{f}}{\partial H}=-1$$
(87)
$$\frac{\partial H}{\partial {\varepsilon }_{v}^{fp}}={p}_{0}\beta exp\left(\beta {\varepsilon }_{v}^{fp}\right)$$
(88)
$$\frac{\partial {f}^{f}}{\partial \alpha }=-\frac{n{\sigma }_{m}^{f}{({\eta }^{f}/\alpha )}^{n}}{\alpha {\left[1-{({\eta }^{f}/\alpha )}^{n}\right]}^{2}}$$
(89)
$$\frac{\partial \alpha }{\partial {\varepsilon }_{s}^{fp}}=-\frac{{\alpha }_{m}{c}_{1}}{{c}_{2}}exp\left({\varepsilon }_{s}^{fp}/{c}_{2}\right)$$
(90)

The plastic volumetric strain and shear strain of frictional elements are:

$${d\varepsilon }_{v}^{fp}=\frac{1}{h}\left(\frac{\partial {f}^{f}}{\partial {\sigma }_{m}^{f}}\frac{\partial {g}^{f}}{\partial {\sigma }_{m}^{f}}{d\sigma }_{m}^{f}+\frac{\partial {f}^{f}}{\partial {\sigma }_{s}^{f}}\frac{\partial {g}^{f}}{\partial {\sigma }_{m}^{f}}{d\sigma }_{s}^{f}\right)$$
(91)
$${d\varepsilon }_{s}^{fp}=\frac{1}{h}\left(\frac{\partial {f}^{f}}{\partial {\sigma }_{m}^{f}}{\frac{\partial {g}^{f}}{\partial {\sigma }_{s}^{f}}d\sigma }_{m}^{f}+\frac{\partial {f}^{f}}{\partial {\sigma }_{s}^{f}}\frac{\partial {g}^{f}}{\partial {\sigma }_{s}^{f}}{d\sigma }_{s}^{f}\right)$$
(92)

The total stress–strain formula of frictional elements is as follows:

$${d\varepsilon }_{v}^{f}=\frac{{d\sigma }_{m}^{f}}{{K}^{f}}+\frac{1}{h}\left(\frac{\partial {f}^{f}}{\partial {\sigma }_{m}^{f}}\frac{\partial {g}^{f}}{\partial {\sigma }_{m}^{f}}{d\sigma }_{m}^{f}+\frac{\partial {f}^{f}}{\partial {\sigma }_{s}^{f}}\frac{\partial {g}^{f}}{\partial {\sigma }_{m}^{f}}{d\sigma }_{s}^{f}\right)$$
(93)
$${d\varepsilon }_{s}^{f}=\frac{{d\sigma }_{s}^{f}}{{3G}^{f}}+\frac{1}{h}\left(\frac{\partial {f}^{f}}{\partial {\sigma }_{m}^{f}}{\frac{\partial {g}^{f}}{\partial {\sigma }_{s}^{f}}d}\sigma_{m}^{f}+\frac{\partial {f}^{f}}{\partial {\sigma }_{s}^{f}}\frac{\partial {g}^{f}}{\partial {\sigma }_{s}^{f}}{d\sigma }_{s}^{f}\right)$$
(94)

By comparing formula (32) and (33), we can obtain:

$${K}_{v}^{f}=\frac{1}{{K}^{f}}+\frac{1}{h}\frac{\partial {f}^{f}}{\partial {\sigma }_{m}^{f}}\frac{\partial {g}^{f}}{\partial {\sigma }_{m}^{f}}$$
(95)
$${K}_{s}^{f}=\frac{1}{h}\frac{\partial {f}^{f}}{\partial {\sigma }_{s}^{f}}\frac{\partial {g}^{f}}{\partial {\sigma }_{m}^{f}}$$
(96)
$${G}_{v}^{f}=\frac{1}{h}\frac{\partial {f}^{f}}{\partial {\sigma }_{m}^{f}}\frac{\partial {g}^{f}}{\partial {\sigma }_{s}^{f}}$$
(97)
$${G}_{s}^{f}=\frac{1}{{3G}^{f}}+\frac{1}{h}\frac{\partial {f}^{f}}{\partial {\sigma }_{s}^{f}}\frac{\partial {g}^{f}}{\partial {\sigma }_{s}^{f}}$$
(98)

The samples are all composed of frictional elements after complete damage. It can be approximately considered that the soil after the completion of the first test has been completely destroyed. Therefore, the parameters of frictional elements can be obtained through the consolidated drained triaxial test of the reshaped samples.\({K}^{f}\) and \({G}^{f}\) were approximately determined by the initial slope of \({\varepsilon }_{v}^{f}-{\sigma }_{m}^{f}\) and \({\varepsilon }_{s}^{f}-{\sigma }_{s}^{f}\) curves of the reshaped samples. \(\lambda {\prime}\) andκ were determined by the e-lnp curves of the reconstructed samples. \({\alpha }_{m}\), c1, c2, n, n1 were determined by curve fitting of the triaxial test of the reshaped sample.

Constitutive relation of roots

Assuming that the root is also an ideal elastic-brittle material, the linear elastic constitutive model is adopted and Hooke's law is followed, as shown in Eqs. (30) and (31). And the strength and modulus of the root are large enough that the root will not be destroyed before the sample is destroyed. \({K}^{r}\) and \({G}^{r}\) can be measured directly or indirectly through laboratory tests.

Structural parameters

Breakage ratio

Λ is an internal variable, as a structural mesoscopic parameter of soil, and is related to the strength of the bonded elements, confining pressure, external load, stress path, so it can not be measured directly. The breakage ratio caused by spherical stress is volume breakage ratio \({\lambda }_{v}^{f}\), and the breakage ratio caused by deviatoric stress is area breakage ratio \({\lambda }_{s}^{f}\). According to the definition and physical meaning of the breakage ratio, it can be seen that the breakage ratio is an undecreasing function varying in the range of 0 to 1. When there are roots, the breakage ratio is limited to \(1-{\lambda }_{v}^{r}\). When the specimen is strain softening, the breakage ratio does not increase after reaching the peak pressure and remains unchanged. Therefore, the definition similar to damage variable can be adopted, Weibull distribution function is used to describe the change rule of breakage ratio, and it is assumed that \({\lambda }_{v}^{f}\) and \({\lambda }_{s}^{f}\) conform to the change rule of the following expression:

$${\lambda }_{v}^{f}={f}_{1}\left({\sigma }_{m},N\right)=1-{\lambda }_{v}^{r}-\left(1-{\lambda }_{v}^{r}\right)exp\left\{-{\alpha }_{v}{\left(\frac{{\sigma }_{m}}{{P}_{a}}\right)}^{{\theta }_{v}}-{\beta }_{v}{N}^{{\omega }_{v}}\right\} \;\; \text{when} \;\;{d\sigma }_{m}>0$$
(99)
$${\lambda }_{v}^{f}={\lambda }_{vmax }^{f} \;\; \text{when} \;\;d{\sigma }_{m}\le 0$$
(100)
$${\lambda }_{s}^{f}={f}_{2}\left({\sigma }_{s},N\right)=1-{\lambda }_{s}^{r}-\left(1-{\lambda }_{s}^{r}\right)exp\left\{-{\alpha }_{s}{\left(\frac{{\sigma }_{s}}{{P}_{a}}\right)}^{{\theta }_{s}}-{\beta }_{s}{N}^{{\omega }_{s}}\right\} \;\; \text{when} \;\; {d\sigma }_{s}>0$$
(101)
$${\lambda }_{s}^{f}={\lambda }_{smax }^{f}\;\; \text{when} \;\; d{\sigma }_{m}\le 0$$
(102)

in which \({P}_{a}\) is the atmospheric pressure, \({\alpha }_{v}, {\beta }_{v}\), \({\theta }_{v}, {\omega }_{v}\), \({\alpha }_{s}\), \({\beta }_{s}\), \({\theta }_{s}\), \({\omega }_{s}\) is the test parameter, fitted by the triaxial test curve of the sample.

From Eqs. (99) to (102), we can get:

$$\frac{\partial {f}_{1}}{\partial {\sigma }_{m}}={\frac{1-{\lambda }_{v}^{r}}{{P}_{a}}}\alpha_{v}{{\theta }_{v}\left(\frac{{\sigma }_{m}}{{P}_{a}}\right)}^{{\theta }_{v}-1}exp\left\{-{\alpha }_{v}{\left(\frac{{\sigma }_{m}}{{P}_{a}}\right)}^{{\theta }_{v}}-{\beta }_{v}{N}^{{\omega }_{v}}\right\} \;\; \mathrm{when } \;\; {d\sigma }_{m}>0$$
(103)
$$\frac{\partial {f}_{1}}{\partial {\sigma }_{m}}=0, \;\; \mathrm{ when } \;\; d{\sigma }_{s}\le 0$$
(104)
$$\frac{\partial {f}_{1}}{\partial N}=\left(1-{\lambda }_{v}^{r}\right){\beta }_{v}{{\omega }_{v}N}^{{\omega }_{v}-1}exp\left\{-{\alpha }_{v}{\left(\frac{{\sigma }_{m}}{{P}_{a}}\right)}^{{\theta }_{v}}-{\beta }_{v}{N}^{{\omega }_{v}}\right\}$$
(105)
$$\frac{\partial {f}_{2}}{\partial {\sigma }_{s}}={\frac{1-{\lambda }_{s}^{r}}{{P}_{a}}}\alpha_{s}{\theta }_{s}{\left(\frac{{\sigma }_{s}}{{P}_{a}}\right)}^{{\theta }_{s}-1}exp\left\{-{\alpha }_{s}{\left(\frac{{\sigma }_{s}}{{P}_{a}}\right)}^{{\theta }_{s}}-{\beta }_{s}{N}^{{\omega }_{s}}\right\} \;\; \mathrm{when } \;\; {d\sigma }_{s}>0$$
(106)
$$\frac{\partial {f}_{2}}{\partial {\sigma }_{s}}=0 \;\; \mathrm{when} \;\; d{\sigma }_{s}\le 0$$
(107)
$$\frac{\partial {f}_{2}}{\partial N}=\left(1-{\lambda }_{s}^{r}\right){\beta }_{s}{\omega }_{s}{N}^{{\omega }_{s}-1}exp\left\{-{\alpha }_{s}{\left(\frac{W}{{P}_{a}}\right)}^{{\theta }_{s}}-{\beta }_{s}{N}^{{\omega }_{s}}\right\}$$
(108)
Coefficient of local stress concentration

The local stress concentration coefficient establishes the relationship between the stress of cemented elements and root and the stress of representative elements. It is also an unmeasurable internal variable, which is related to the external load on the specimen and its internal stress–strain coordination ability. Although the bonded elements and the root have similar constitutive relationship, their local stress concentration coefficients are different because of their different elastic modulus and different strain coordination ability. Similar to the definition of breakage ratio, the local stress concentration coefficient can be divided into spherical stress concentration coefficient and deviatoric stress concentration coefficient. The stress concentration coefficient of bonded elements meets the requirement that when the breakage ratio is 0, the stress concentration coefficient is 1. Ignoring the influence of root, the simplified assumption is:

$${C}_{m}^{b}=1+{m}_{1}{\lambda }_{v}^{f}$$
(109)
$${C}_{s}^{b}=1+{m}_{2}{\lambda }_{s}^{f}$$
(110)

The volume of the root is constant during the stress loading process and the stress concentration coefficient of the root is simply assumed to be a variable related to the confining pressure.

Model verification

Under triaxial test conditions, the mean stress, generalized shear stress, volumetric strain and shear strain have the following relations:

$${\varepsilon }_{v}={\varepsilon }_{1}+2{\varepsilon }_{3}$$
(111)
$${\varepsilon }_{s}=\frac{2}{3}\left({\varepsilon }_{1}-{\varepsilon }_{3}\right)$$
(112)
$${\sigma }_{m}=\frac{1}{3}\left({\sigma }_{1}+2{\sigma }_{3}\right)$$
(113)
$${\sigma }_{s}={\sigma }_{1}-{\sigma }_{3}$$
(114)

The actual number of freeze–thaw cycles is a discrete quantity, so it is verified under the same number of freeze–thaw cycles, that is, dN = 0. The above relation can be obtained by substituting it into Eqs. (54) and (55):

$${d\varepsilon }_{v}=d\left({\varepsilon }_{1}+2{\varepsilon }_{3}\right)={A}_{1}d\frac{1}{3}\left({\sigma }_{1}+2{\sigma }_{3}\right)+{B}_{1}d\left({\sigma }_{1}-{\sigma }_{3}\right)$$
(115)
$${d\varepsilon }_{s}=\frac{2}{3}d\left({\varepsilon }_{1}-{\varepsilon }_{3}\right)={A}_{2}d\frac{1}{3}\left({\sigma }_{1}+2{\sigma }_{3}\right)+{B}_{2}d\left({\sigma }_{1}-{\sigma }_{3}\right)$$
(116)

Due to \({d\sigma }_{3}=0\), we can get:

$$d{\varepsilon }_{3}=\frac{\left(\frac{{A}_{1}}{3}+{B}_{1}\right)-\frac{3}{2}\left(\frac{{A}_{2}}{3}+{B}_{2}\right)}{\left(\frac{{A}_{1}}{3}+{B}_{1}\right)+3\left(\frac{{A}_{2}}{3}+{B}_{2}\right)}{d\varepsilon }_{1}$$
(117)
$${d\sigma }_{1}=\frac{3}{\left(\frac{{A}_{1}}{3}+{B}_{1}\right)+3\left(\frac{{A}_{2}}{3}+{B}_{2}\right)}d{\varepsilon }_{1}$$
(118)

According to the triaxial test results, the basic parameters of the bonded elements, frictional elements and root were respectively determined by the above method, and then the internal state parameters were determined by the inversion of the triaxial test results. The parameters of the revised model of binary medium for soil samples containing roots under different confining pressures and different freeze–thaw cycles were determined as shown in Table 6.

Table 6 Parameter determination of extended binary medium model.
Full size table

The parameters determined in Table 6 were put into the stress–strain formula of the model for calculation, the test results were predicted, and curves \({\varepsilon }_{a}-\left({\sigma }_{1}-{\sigma }_{3}\right)\) and \({\varepsilon }_{a}-{\varepsilon }_{v}\) under different confining pressures and different times of freeze–thaw cycles with different volume ratio of roots were simulated. The simulation curves of part of the prediction are compared with the test curves, and the results are shown in Figs. 14, 15, 16, 17, 18, 19, where T represents the test and C represents the simulation calculation. It can be seen from the figure that the predicted results of the model are in good agreement with the triaxial test results of the samples, especially the strain softening and volumetric dilatancy phenomena of the samples can be simulated relatively well.

Figure 14
figure 14figure 14

Comparison of \({\varepsilon }_{a}-\left({\sigma }_{1}-{\sigma }_{3}\right)\) curves under different confining pressures between test and simulation.

Full size image
Figure 15
figure 15figure 15

Comparison of \({\varepsilon }_{a}-{\varepsilon }_{v}\) curves under different confining pressures between test and simulation.

Full size image
Figure 16
figure 16

Comparison of \({\varepsilon }_{a}-\left({\sigma }_{1}-{\sigma }_{3}\right)\) curves under different freeze–thaw cycles between test and simulation.

Full size image
Figure 17
figure 17

Comparison of \({\varepsilon }_{a}-{\varepsilon }_{v}\) curves under different freeze–thaw cycles between test and simulation.

Full size image
Figure 18
figure 18

Comparison of \({\varepsilon }_{a}-\left({\sigma }_{1}-{\sigma }_{3}\right)\) curves with different roots between test and simulation.

Full size image
Figure 19
figure 19

Comparison of \({\varepsilon }_{a}-{\varepsilon }_{v}\) curves with different roots between test and simulation.

Full size image

Conclusions

In this paper, the stress–strain characteristics of rooted soil under the action of freeze–thaw cycle were experimentally studied, and the variation law of rooted soil with confining pressure, the number of freeze–thaw cycles and the percentage of rooted soil was obtained. And on the basis of the test, the rooted soil was abstracted into the bonded elements, frictional elements and root of the three-phase composite material, the extended binary medium model was derived and verified. The main conclusions of this paper are as follows:

  1. 1.

    Laboratory triaxial tests show that rooted soil exhibits strain softening and volumetric dilatancy at low confining pressures, and strain hardening and volumetric shrinkage at high confining pressures. During the freeze–thaw cycles, the cement structure in the soil is damaged by frost heave force. With the increase of freeze–thaw times, the bearing capacity of the sample gradually decreases, the deviatoric stress decreases, and the volumetric strain increases. The greater the root volume, the more the bearing capacity of the sample, and the smaller the root volume, so the strengthening effect is limited. Under the same other conditions, the sample strength increases with the increase of confining pressure, decreases with the increase of freeze–thaw times, and increases with the increase of root volume.

  2. 2.

    Based on the theory of breakage mechanics for geological materials, this paper abstracts soil into bonded elements that can be destroyed and frictional elements, and roots into bonded elements that will not be destroyed to jointly bear external load and deformation, and deduces a extended binary medium model of rooted soil under the action of freeze–thaw cycles. In this way, the binary medium model can be extended to the three-phase composite materials, which can adopt different constitutive relations. In this paper, linear elastic constitutive model (Hooke’s law) and double-hardening constitutive model are adopted respectively according to their respective characteristics. The simulation results and experimental comparison show that the model can simulate the mechanical properties of rooted soil under freeze–thaw cycle.

  3. 3.

    The extension from binary medium to three-phase material inevitably brings about an increase in the number of model parameters, so the parameters are not easy to determine and the error may increase, which requires more accurate and careful tests and as much as possible to establish more relationships to reduce parameters. In addition, the characteristics of concrete materials are more in line with the characteristics of this model. Aggregate is regarded as bonded elements that will not be destroyed, and cement mortar gradually breaks under external load and transforms from bonded elements to frictional elements. The applicability of this model to concrete can be further studied in the later stage.