This study selected typical clay from China as the tested soil. According to the Chinese Specification of JTG 3430-2020, Test Methods of Soils for Highway Engineering, the main basic physical properties of the tested soils were listed in Table 1. On the basis of these parameters, the soils were determined as low liquid limit clay (CL) in accordance with the Unified Soil Classification System.

Limited by the difficulty of achieving undisturbed soils, and in order to ensure consistent humidity and density levels between specimens, the soils were remolded and tested in the laboratory. Firstly, the soils were crushed and collected after filtering by the 2-mm sieve, then dried at 105°C for 12 h. Subsequently, the dry soils were mixed with distilled water at the optimum moisture content and then stored in a sealed plastic bag for 24 h so that the moisture was uniform. The principle of this specimen preparation method is to simulate the actual humidity inside most embankments (Lin et al., 2017; Bozbey et al., 2018). At the same time, the dry density of the specimens was maintained at 1.73 g/cm³. Afterwards, the wet soil with the calculated mass was compacted by five layers into a cylinder that was 68 mm in diameter and 140 mm in height. The rationale behind this was that ASTM D2850 states the height-to-diameter ratio should be greater than 2. Finally, the specimens were wrapped to avoid water evaporation and waited for the next tasks.

The closed-system FT cycle experiments were applied in this study. This method is suitable for depicting the condition that there are no obvious differences in the in-situ moisture content all year round. Meanwhile, given the low permeability of clay, the rate of frost penetration is significantly faster than moisture transportation, so that there is no adequate time to achieve persistent replenishment of moisture to the frozen front during the freezing process (Güllü and Khudir, 2014). Thus, it is reasonable to adopt the closed FT type here.

As shown in Figure 1, the entire FT cycle tests were performed in the automatic refrigerator (the range was from −60 to 150°C) with internal temperature sensors (the precision was ±0.5°C). The prepared specimens were encased in preservative film to prevent moisture exchange. The specimens were first frozen at −20°C for 12 h, and following this period, they were thawed at 20°C for 12 h. The two processes were integrated into one FT cycle, and so on. In earlier pertinent reports, it can be found that these temperatures and durations were approved by some scholars (Orakoglu and Liu, 2017; Nguyen et al., 2019; Gowthaman et al., 2020; He et al., 2020; Lu et al., 2020). The numbers of FT cycles in this study were designed as 0, 1, 3, 6, 10, and 15.

DYNATRIAX-100/14 (Figure 2), a triaxial test apparatus, was employed to investigate the mechanical properties of clay in this study. Its axial load limit was 100 kN, the confining pressure could be adjusted from 0 to 1000 kPa, and the maximum axial distance of the actuator was 100 mm. The external transducer with the highest frequency of 1000 Hz collected real-time displacement and load data for the computer.

When the specimens had been subjected to the preset FT cycles, they were transferred to the triaxial chamber for UU tests. The axial loading rate was 0.8% per min and the values of confining pressures were set at 50, 100, 150, 200, 300, and 400 kPa, respectively. Once the axial strain of the specimens reached 20%, the test ended.

Figure 3 displays the stress-strain curves of the tested soils with various FT cycles and confining pressures (σ ₃). For strain softening soils, with the development of axial strain, the deviator stress first reaches its peak point and then declines gradually. While soils exhibit strain hardening characteristics, the deviator stress increases along the axial. As for strain stabilization curves, the deviator stress rise to a stable value with the accumulation of axial strain, and then continues to stabilize. Accordingly, three types of stress-strain responses can be observed in the experimental results. Moreover, at the identical confining pressures, it can be found that the more FT cycles, the smaller the peak value of softening curves, and after the 10th FT cycle, this change is not considerable.

Figure 3 also shows that, with the given FT cycles, the rise in confining pressures causes greater deviator stress when the soils achieve the same axial strain. This indicates that the constraint effect contributed by confining pressures strengthens the soils to better resist axial loads. In addition, the behaviors of strain hardening or stabilization occurred in the absence of FT cycles while the softening appeared with FT cycles. In other words, FT cycles can transform the types of stress-strain curves of the soils in this study. This phenomenon is consistent with the findings reported by Roustaei et al. (2015) and Orakoglu and Liu (2017), but different from Xian et al. (2019). One possible explanation is that the different clays with their own diverse properties were selected in this research, resulting in differences in the numbers of FT cycles and confining pressures corresponding to the critical conditions when the stress-strain curves transition.

Elastic modulus is a key indicator for describing the stiffness and anti-deformation ability of soils. For instance, in a recent study, Ling et al. (2020) discovered that railway coarse-grained materials that had undergone 10 FT cycles lost 68, 57, and 50% of their elastic modulus at confining pressures of 30, 60, and 90 kPa, respectively. For simplicity, the stress-strain curves were considered to be linear within the 2.0% axial strain, and the slope of this point was defined as the elastic modulus value (Wei et al., 2019; Liu et al., 2021), which can be determined in Eq. 1: E = ( σ 2.0 % − σ 0 ) / ( ε 2.0 % − ε 0 ) (1) where E denotes the elastic modulus; σ _2.0% is the deviator stress when the axial strain reaches 2.0% (ε _2.0%); and σ ₀ and ε ₀ are the initial stress and strain, respectively.

Figure 4 reveals the relationship between the elastic modulus and confining pressures with different FT cycles. Taking the state of 6 FT cycles as an instance, when the confining pressure ranged from 50 kPa to 100, 150, 200, 300, and 400 kPa, the elastic modulus increased by 20, 33, 51, 67, and 85%, respectively. This analysis result means that the confining pressure has a considerable positive effect on the elastic modulus of the studied clay. On the other hand, it can be observed that the accumulated FT cycles declined the elastic modulus. Furthermore, the value decayed sharply during the first FT cycle and stabilized after a certain number of cycles, in agreement with most research on clay in the past decade (Roustaei et al., 2015; Wang et al., 2007; Xian et al., 2019). Based on the foregoing analysis, to better understand the elastic modulus of the studied clay, this study suggested a novel empirical model through the multiple regression modeling of all 36 data points in Figure 4, especially adopting the modificatory logarithmic form to embody the effect of FT cycles, as shown in Eq. 2. The result is satisfactory in reference to the recommended criteria for evaluating the determination coefficient (R²) in geotechnical engineering (Zhang et al., 2020b; Zhang et al., 2021; Yao et al., 2022), namely, a fair fit is defined as 0.4 ≤ R² ≤ 0.69, a good fit is R² = 0.7–0.89 and an excellent fit is R² ≥ 0.9. E = 1.83 [ L n ( e + N ) ] − 0.15 ( σ 3 + 99.52 ) 0.49 ( R 2 = 0.95 ) (2) where N is the number of FT cycles.

This study ascertained the failure strength of clay in accordance with the instructions of ASTM D2166: for strain softening clay, the deviator stress corresponding to the peak of the stress-strain curve is the failure strength; and if strain hardening or stabilization occurs, the deviator stress at 15% of the axial strain is taken as the failure strength. The results of the failure strength due to various confining pressures and FT cycles are presented in Figure 5. It can be seen that the addition of confining pressures continuously improved the failure strength. Also, the failure strength weakened throughout the FT cycles and then basically tended to be stable after 10 FT cycles with an approximate attenuation rate of 6–22% compared to the state before FT cycles. This law may be produced as follows: the water in pores freezes into ice at a negative temperature, and the volume of ice is roughly 9% larger than liquid water (Zhou and Wei, 2020). The growth of ice crystals uplifted the soil particles, but when ice crystals thaw at positive temperatures, the lifted particles cannot be completely restored (Zhou et al., 2018). Thus, a complete process comprising freezing and thawing resulted in the expansion of pores between particles and the attenuation of failure strength. Furthermore, with increasing FT cycles, the contribution of ice crystals to the lifting of particles decreased gradually until the pores between the particles were expanded by repeated FT cycles were large enough to endure the volume growth due to the phase transitions in water, the failure strength was no longer attenuated. Given the above discussions, a new model can be conceived as Eq. 3. ( σ 1 − σ 3 ) lim = k 1 P a [ L n ( e + N ) ] k 2 ( σ 3 p a + 1 ) k 3 (3) where ( σ 1 − σ 3 ) lim is the failure strength, σ ₁ is the major principal stress; P _a = 101.3 kPa, which is the atmospheric pressure; k ₁ , k _2, and k ₃ are regression coefficients.

In the stability analysis of embankments, the failure strength of soils must be strictly considered (Zhao et al., 2020). Obviously, it is necessary and meaningful to verify the robustness and accuracy of Eq. 3. As for regression coefficients, k ₁ is the adjustment coefficient, that is, proportional to the failure strength, which can be ascertained to be a positive number because the failure strength cannot be negative. k ₂ as an exponent must be a negative number given that with the increase in the number of FT cycles, the failure strength decreases. In order to reflect the influence of confining pressures, k ₃ is a positive number as well as an exponent, since the higher the confining pressure, the larger the failure strength. This study adopted all 36 data points in Figure 5 and, in conjunction with previous data from other researchers, fitted Eq. 3 accordingly. The specific regression coefficients are summarized in Table 2, and the intuitionistic fitting effect is presented in Figure 6. The results demonstrated that the signs of regression coefficients in each group were consistent with the discussions, and the accuracies were mainly excellent. Simultaneously, the values of the fitted and experimented failure strength were close. The verification analysis concludes that the model established as Eq. 3 is suitable for different clays, although it was proposed only on the basis of test results from the low liquid limit clay in this study.

The shear strength of soils, specifically cohesion and internal friction angle, are the most popular parameters in embankment engineering, as they are the basic inputs necessary for calculating the slope safety factor. Currently, the common practice is to plot the Mohr semicircles and envelopes based on the failure strength determined by triaxial tests to obtain the shear strength (Zhou et al., 2019). Nonetheless, triaxial instruments are still far from universal due to their high cost. Furthermore, triaxial tests require experienced personnel to operate procedures and analyze test data, which is a complex and time-consuming process. Contrary to this, trying to propose a method for rapidly predicting the shear strength through performance-related soil basic parameters is advisable due to the measurements of these parameters being faster, more simple, reliable, and accurate.

In the above scenario, the weight percent of soils passing through 0.075-mm sieve (P _0.075), the plastic limit (w _p ) and the liquid limit (w _L ) were selected as physical parameters; relative humidity (RH, the ratio of actual moisture content to optimal moisture content) and relative compaction (RC, the ratio of actual dry density to maximum dry density) were taken as status parameters, because they had a considerable effect on the mechanical properties of soils (Nazzal and Mohammad, 2010). Corresponding to the various soils in Table 2, these detailed parameters were summarized in Table 3. Subsequently, an excellent multiple regression relationship between these basic parameters and the model regression coefficients (k ₁ , k ₂ and k ₃ ) of Eq. 3 was developed utilizing the JMP statistical analysis software, as presented in Eq. 4. { k 1 = − 12.44 + 0.08 P 0.075 − 0.94 w p + 0.41 w L + 0.01 R H + 0.14 R C ( R 2 =   0.97 ) k 2 = 1.46 − 0.01 P 0.075 + 0.07 w p − 0.01 w L − 0.002 R H − 0.02 R C ( R 2 =   0.92 ) k 3 = − 9.57 − 0.003 P 0.075 + 0.05 w p − 0.02 w L + 0.01 R H + 0.09 R C ( R 2 =   0.95 ) (4)

For the particular moment when the soils reach the failure state, the failure strength expressed by Eq. 3 can also be rewritten into Eq. 5: g ( σ 1 , σ 3 ) = σ 1 − σ 3 − k 1 P a [ L n ( e + N ) ] k 2 ( σ 3 p a + 1 ) k 3 = 0 (5) where g ( σ 1 , σ 3 ) is the function of σ ₁ and σ ₃ at the failure moment with FT cycles.

Most failures in geotechnical engineering can be attributed to shear action, and thus it is meaningful to describe the failure state by the normal stress and shear stress. This study adopted the widely used Mohr stress circle to reflect this relationship, as shown in Eq. 6: f ( σ 1 , σ 3 ) = ( σ − σ 1 + σ 3 2 ) 2 + τ 2 − ( σ 1 − σ 3 2 ) 2 = 0 (6) where f ( σ 1 , σ 3 ) is the function of σ ₁ and σ ₃ related to σ and τ, σ is the normal stress, and τ is the shear stress.

On the basis of the envelope theory, the g ( σ 1 , σ 3 ) and f ( σ 1 , σ 3 ) have the formulation in Eq. 7 (Yang et al., 2013): ∂ f ∂ σ 1 ∂ g ∂ σ 3 − ∂ f ∂ σ 3 ∂ g ∂ σ 1 = 0 (7)

Substituting Eq. 5, Eq. 6 into Eq. 7, the normal stress and shear stress on the shear plane can be derived as given in Eq. 8: { σ = σ 1 + σ 3 { 1 + k 1 k 3 [ L n ( e + N ) ] k 2 ( σ 3 P a + 1 ) ( k 3 − 1 ) } 2 + k 1 k 3 [ L n ( e + N ) ] k 2 ( σ 3 P a + 1 ) ( k 3 − 1 ) τ = ( σ 1 − σ 3 2 ) 2 − { σ 1 + σ 3 + k 1 k 3 σ 3 [ L n ( e + N ) ] k 2 ( σ 3 P a + 1 ) ( k 3 − 1 ) 2 + k 1 k 3 [ L n ( e + N ) ] k 2 ( σ 3 P a + 1 ) ( k 3 − 1 ) − σ 1 + σ 3 2 } 2 (8)

The rationality of Eq. 8 is logical, and the reasons are analyzed as follows. The regression model of failure strength established in this study, as shown in Eq. 3, has been verified to be reliable by multiple groups of clay data. Meanwhile, the Mohr stress circle and envelope theory have been proven to be successful in existing documents. Therefore, the normal stress and shear stress obtained by them together should be accurate. In order to have confidence in Eq. 8, the measured failure strength data of this study were used to draw the Mohr semicircles, and the set numbers of FT cycles (N), confining pressures (σ ₃), fitted k ₁ , k ₂ and k ₃ in Table 2, and Eq 5 and Eq. 8 of this study were adopted to calculate the envelope, as shown in Figure 7. It can be observed that the calculated envelopes were basically anastomotic with the Mohr semicircles, which indicates that the Eq. 8 obtained in this section can well describe the nonlinear relationship between the shear strength of clay and the numbers of FT cycles and confining pressures.

To sum up, the critical steps in the rapid prediction method for the shear strength of clay considering FT cycles were summarized and organized as follows:

1) Through a few simple laboratory basic physical tests to determine the values of P _0.075, w _p , w _L , RH and RC, and then substitute them into Eq. 4 to calculate k ₁ , k ₂ and k ₃ ;

Following the above steps, the results of the cohesion and internal friction angle were summarized in Figure 8. It can be observed that the cohesion of the tested clay decayed throughout a whole FT cycles, and the rate of attenuation gradually decreased until it became flat. This phenomenon has also been found in previous studies (Orakoglu and Liu, 2017; Wei et al., 2019). Figure 8 also showed that the internal friction angle of the tested clay was not sensitive to the numbers of FT cycle, that is, it fluctuated repeatedly around 23°, which was consistent with the studies from Liu et al. (2016), Orakoglu and Liu (2017) and Zhang et al. (2019), but unlike the findings from Zhou et al. (2019) and Xian et al. (2019). The difference was likely due to the effect of FT cycles do not cause particles rearrangement in the tested clay such as interlocking or slipping. Not long ago, a measurement on the micro-structure of pavement recycled materials seemed to indicate that the X-ray Micro-Computed Tomography technology might be an effective means of exploring this difference (Afshar et al., 2018), however it is beyond the scope of this study.

It is noteworthy that the method proposed in this work for rapidly predicting the shear strength of embankment clay under FT cycles is generally applicable. Concretely, when this method is applied to other clay, it is only necessary to carry out several simple basic physical tests and follow the above steps 1)–3) to achieve the purpose. For situations where the triaxial tests are limited or the shear strength of clay is urgently needed, this method has satisfactory engineering practicability.

The studied clay before FT cycles in Figure 3A exhibited obvious strain hardening or stabilization behavior. Here the Konder hyperbolic model as shown in Eq 9a, Eq. 9b was adopted to describe the stress-strain relationship, and the results were listed in Table 4. Although this model showed high accuracy, it could not easily reflect the influence of confining pressures. In other words, for different confining pressure levels, corresponding model coefficients need to be fitted separately, which increases the difficulty and complexity of operation. ( σ 1 − σ 3 ) = ε 1 a + b ε 1 (9a) ε 1 σ 1 − σ 3 = a + b ε 1 (9b) where ε ₁ represents the axial strain; 1/a is the initial tangent modulus E _tan; and 1/b is the ultimate deviator stress (σ ₁−σ ₃)_ult.

Casey and Germaine (2013) reported that the stress-strain relationship of soils under different conditions can be expressed by a normalized function. Based on this idea, multiply both sides of Eq. 9b by a normalized factor of N to obtain Eq. 10: ε 1 N σ 1 − σ 3 = a N + b N ε 1 = N E tan + N ( σ 1 − σ 3 ) ult ε 1 (10)

To implement the normalized function, the N should satisfy a linear relationship with the E _tan and (σ ₁−σ ₃)_ult, respectively. Through a series of attempts on the mechanical indicators obtained in Section 3, it was found that the elastic modulus before FT cycles (E ₀) met this normalized condition, as shown in Figure 9.

E ₀ was selected as the normalized factor N, based on Eq. 10, and Eq. 11 can be deduced: ε 1 E 0 σ 1 − σ 3 = E 0 1.6631 E 0 + E 0 0.0339 E 0 ε 1 = 0.601 + 29.496 ε 1 (11)

In addition, the E ₀ also had a relationship with confining pressures as shown in Figure 10 can be depicted by Eq. 12: E 0 = 39.941 σ 3 + 23662 ( R 2 = 0.95 ) (12)

Substituting Eq. 12 into Eq. 11, the stress-strain relationship of the studied clay before FT cycles can be expressed as Eq. 13: ( σ 1 − σ 3 ) 0 = ( 39.941 σ 3 + 23662 ) ε 1 0.601 + 29.496 ε 1 (13)

The studied clay under FT cycles in Figure 3B to Figure 4F demonstrated strain softening behavior, and thus the Konder model cannot be adopted to describe their stress-strain relationship directly. For this reason, this study introduced a new variable, namely, the deviator stress differences caused by FT cycles, and defined it in Eq. 14. The results of deviator stress differences with various FT cycles are plotted in Figure 11. Δ ( σ 1 − σ 3 ) = ( σ 1 − σ 3 ) 0 − ( σ 1 − σ 3 ) i (14) where Δ ( σ 1 − σ 3 ) is the deviator stress difference; ( σ 1 − σ 3 ) 0 and ( σ 1 − σ 3 ) i are respectively the deviator stress corresponding to the same axial strain before FT cycles and after the i-th FT cycles.

From Figure 11, the strain hardening behavior was found. Hence, the Konder model was applied to them, as expressed in Eq. 15, and the results are exhibited in Table 5. ε 1 Δ ( σ 1 − σ 3 ) = c + d ε 1 (15) where 1/c denotes the initial tangent modulus Δ E tan corresponding to the Δ ( σ 1 − σ 3 ) ; and 1/d is th15e ultimate deviator stress differences [ Δ ( σ 1 − σ 3 ) ] ult .

Following the processing method in Section 4.1, in order to obtain a normalized relationship between the deviator stress differences and axial strain, the normalized factor here needs to have a good linear relationship with both the Δ E tan and the [ Δ ( σ 1 − σ 3 ) ] ult . It can be discovered that the differences in elastic modulus compared with those before FT cycles ( Δ E 0 ) adapted to this condition, as illustrated in Figure 12.

∆E ₀ was chosen as the normalized factor, and Eq. 16 can be given: ε 1 Δ E 0 Δ ( σ 1 − σ 3 ) = Δ E 0 0.9739 Δ E 0 + Δ E 0 0.0717 E 0 ε 1 = 1.027 + 13.947 ε 1 (16)

To eliminate the intermediate variable Δ E 0 , this study was aided by the variation rate (V) of elastic modulus corresponding to various confining pressures and numbers of FT cycles, as stipulated in Eq. 17. The results were presented in Figure 13 can be formulated in Eq. 18. V = Δ E 0 E 0 (17) V = 1.38 σ 3 − 0.49 N 0.15 ( R 2 = 0.94 ) (18)

Combining Eqs 12–14, 16–18, the normalized empirical constitutive equation for the tested clay considering FT cycles can be obtained as: ( σ 1 − σ 3 ) i = ( 39.941 σ 3 + 23662 ) ε 1 0.601 + 29.496 ε 1 − 1.38 σ 3 − 0.49 N 0.15 ( 39.941 σ 3 + 23662 ) ε 1 1.027 + 13.947 ε 1 (19)

For the special case N = 0, the Eq. 19 is completely transformed to Eq. 13, and thus Eq. 19 is also applicable to the condition before FT cycles. Then Eq. 19 was selected to calculate the stress-strain curves with the corresponding confining pressures and numbers of FT cycles. The tested and calculated results are shown in Figure 14, demonstrating that the normalized Eq. 19 can preliminarily describe the stress-strain relationship of the tested clay affected by FT cycles, including strain hardening, stabilization, and softening forms.