Stress–Strain–Time Description and Analysis of Frozen–Thawed Silty Clay under Low Stress Level

Abstract

The construction of high-speed railways in cold regions needs to consider the effects of freeze–thaw cycles (FTHs) on the long-term deformation of subgrades. However, at present, research on the creep characteristics of frozen–thawed rocks and soils is not extensive. In the limited studies on frozen–thawed soil creep properties, current research focuses more on high stress–strain–time responses, but for the subgrades, the inner stress is usually low. This paper presents the results of triaxial compression creep tests on remolded, frozen–thawed silty clay sampled in the Yichun-Tieli area and describes its stress–strain–time relationship in an arctan function-based mathematical model. Each creep test condition is studied using three specimens. Frozen–thawed silty clay exhibits attenuation creep under low-level stress. In general, from 4 FTHs to 11 FTHs, the mean elasticity modulus decreases first, and then, increases. The exerted stress is higher than the yield stress; the more FTHs the specimens experience, the more time they need to be deformed stably under the same axial deviatoric stress (ADS). Under the same mean ADS, the mean stable strain of 7 FTHs exceeds the other two FTH conditions and, in general, the mean stable strain of 4 FTHs exceeds 11 FTHs. By dissecting the phenomena, it can be concluded that with FTHs increasing, moisture and voids reconstitute in the process; the elastic strain accounts for most of the total strain and significantly decides the extent of creep deformation; the arctan function-based model is basically able to describe, but not perfectly predict, the stress–strain–time relationship of frozen–thawed silty clay.

1. Introduction

In engineering practice in cold regions, freeze–thaw cycles (FTHs) cause inconvenient damage and destruction to engineering facilities, though sometimes the effect of FTH may have potential. At present, research on the mechanical properties of rock and soil under the action of FTHs is usually based on cohesion, Poisson’s ratio, the internal friction angle, the linear and nonlinear modulus, etc., focusing on the short-term performance of factors such as strength and strain. Research on the laws of rock and soil creep under FTHs is less abundant but significant, and includes preventing railway subgrade upon FTH active stratum or seasonal frozen soil due to uneven settlement, preventing slopes from sliding due to FTHs, keeping the pit wall safe during excavation using the Freezing Method, and keeping the tunnel wall stable during the operation period with years elapsing.

Combining the research results of repeated FTH experiments by Chamberlain, E.J. [1], Eigenbrod, K.D. [2], Viklander, P. [3], He, P., Yang, C.S., et al. [4], it can be concluded that FTHs have a densifying effect on loose soil and a loosening effect on ultra-compact soil, and the dry density or void ratio eventually tend to a certain value [5]. Generally, after about 3 to 18 cycles, the dry bulk density of soil tends to be at a fixed value, and the fixed value is related to the soil type and not related to the initial dry bulk density. Under the condition of an open water supply, the water content of the soil after an FTH is larger than the initial water content, and the increase in water content mainly occurs during the freezing process. When the initial water content is small, the water content of the entire sample after FTH action is greater than the initial water content, while when the initial water content is large, the water content of the sample after freeze–thaw action near the cold end is greater than the initial water content. The moisture content near the bottom of the sample is less than the initial moisture content when freezing occurs for the first time. Under the condition of an open water supply, after repeated FTHs, the soil body is replenished with water from the external environment. Generally, after 3~10 FTHs, the soil water content is greater than the initial water content. The variation in the soil water content of different soil types is not exactly the same. From the perspective of the changing speed of water content after FTHs, when other conditions are the same, generally, the speed of sandy soils is faster and more obvious than that of clay soils. Additionally, specimens with a large amount of water change faster and more significantly. During the freeze–thaw cycle of soil samples, moisture migrates to the cold end. Therefore, with sufficient time and water source, the soil will eventually reach a maximum frozen water content during the long-term freezing process. Generally, FTHs have the effect of strengthening and weakening the mechanical parameters of soil. In most engineering problems, the weakening effect is more prominent. The soil elastic modulus is a key variable in characterizing the ability of the soil structure to resist elastic deformation [6,7]. To sum up, if the FTH temperature is constant and the replacement law of the FTH is unchanged, the weakening effect generally occurs after the first FTH, and the strengthening effect occurs after 2 to 6 FTHs. After several FTHs, the mechanical properties of the soil tend to be stable again, and the elastic modulus tends to be a constant value [5].

Some research investigations about the creep behavior of frozen–thawed rocks and soils are summarized below in

Table 1. FTH creep test objectives and conditions of former researchers.

Source of Rock and Soil	Description of Rock and Soil	Remolded (Re) or Undisturbed (Un)	Water Content (%)	$\mathbf{Density},{\mathit{G}}_{\mathit{s}}$ and e Information	Temperature (°C) and Time (h) of F (Freeze)/TH (Thaw)	FTHs	Creep Test Type Uni/Bi/Triaxial	Size of Specimen (mm)	Creep Type of Describing Model
Dayan Zani Lake Strip Mine Clay, China [19].	Brown and sticky color, slightly powdery, basically composed of clay aggregates; fine, moist, soft plastic, with organic spots, iron streaks, and a small number of conglomerate particles in the lower part.	Re	15, 25, 35	ρ = 2.083 g/cm3 $G_s$= 2.674 e = 0.577	−47~+20	0, 1, 3, 6, 10, 15	Direct shearing creep test	H20/D61.8	Attenuation creep
Rock-like materials made of shale as raw rock, and made of cement, river sand and fly ash [11].	The fly ash content is m = 0%, 10%, 20%, 30%	Re	Saturated		−20~+20 6 h F/6 h TH	0, 5, 10, 20, 40, 60	Uniaxial creep test	H100/D50	Non-attenuation creep
Silty clay and clay from G25 highway subgrade in Jilin, Northeast China [22].	According to the evaluation of the Cu and Cc values of the soil samples, it belongs to well-graded soil. $I_p$ values are 11.25, 16.88 and 23.22, respectively.	Re	16, 17, 18, 19	Compaction degree: 0.8, 0.85, 0.9, 0.95 and the maximum dry density: 1.985, 1.985, 1.944. Two groups of conditions make permutations and combinations	−10~room temperature 24 h F/24 h TH	0, 3, 5, 7	Triaxial creep test	H80/D39.1	Non-attenuation creep
Red sandstone from Sichuan–Tibet railway FTH Area, China [13].		Un	0, 1.5, 2.4, 3.78		−20~+40 1 h F/1 h TH	0, 30, 60, 90, 120	Direct shearing creep test	L50/W50 /H50	Non-attenuation creep
Yellow sandstone mine in Zigong Area, Sichuan, China [14].	The lithology is dense and yellow in color. The main components are anorthite (68%), quartz (18%), andesite cuttings (3%), amphibole (3%), pyroxene (2%) and cement (6%), and the cements are mainly calcite and iron oxide.	Un	Saturated		−20~+20 12 h F/12 h TH	0, 15, 30, 45, 60	Triaxial creep test	H100/D50	Non-attenuation creep
Gneiss in Jilin Huibai Tunnel (located in the middle section of Longgang Mountains at the western foot of Changbai Mountain), China [9].	It is mainly composed of feldspar, quartz and various dark minerals (mica, amphibole, pyroxene, etc.), of which the content of feldspar and quartz is more than 50%.	Un	Saturated		−18~+22 6 h F/ 6 h TH	0, 10, 20, 40, 80	Triaxial creep test	H100/D50	Non-attenuation creep
Granite in the tunnel of Jilin Huibai (Huinan-Baishan Expressway Tunnel), China [10].	Gray-white granite; the main components are quartz, potassium feldspar and acid plagioclase, and the secondary minerals are biotite and amphibole.	Un	Saturated		−20~+20 12 h F/12 h TH	0, 10, 30, 50, 70	Shearing creep test	L200/W100 /H200	Non-attenuation creep
Road subgrade soil under construction in a seasonally frozen area in Liaoning, China [23].		Re			−20~room temperature 48 h F/48 h TH	0, 4, 8, 12	Triaxial creep test	Standard size	Non-attenuation creep
Quartz sandstone from Gansu province, China [15].		Un	Saturated	ρ = 2.14 g/cm3	−10~+204 h F/4 h TH	0, 5, 10	Uniaxial compression creep test	H100/D50	Non-attenuation creep
Gravel and soil collected from the Lianghekou Hydropower Station, China [21].		Re	8	Dry density = 2200 kg/m3	−10~+20 24 h F/12 h TH	0, 1, 2, 5, 10	Confined compression creep test	H200/D300	Non-attenuation creep
Fresh and unbroken red sandstone [16].		Un	Saturated		−20~+20 3 h F/3 h TH	5, 10, 20, 40	Triaxial compression creep test	H100/D50
Red sandstone from the Dafosi Coal Mine, the Binchang Mining Area, Jurassic Coalfield, Huanglong, Shaanxi Province, China [17].	The rock formations consist of mostly thick, weak, medium- to fine-grained sandstones. The sampled rock layer is light brownish-red fine- to medium-grained feldspar quartz sandstone, mainly composed of quartz, plagioclase, potash feldspar and calcite.	Un	Saturated, 5.31%	Dry density = 2.24 g/cm3 Saturated density = 2.35 g/cm3 Porosity = 9.6%	−20~+20 12 h F/12 h TH	0, 1, 5, 9, 13	Triaxial compression tests and multilevel loading and unloading triaxial creep tests	H100/D50
Qinghai–Tibetan clay, from the Beiluhe test site along the Qinghai–Tibetan highway [28].	Plastic limit = 16.5%, Liquid limit = 28.2%, Plastic index = 11.7.	Un	21.4	Maximum dry density = 1.82 g/cm3	−20~+20	0, 1, 2, 4, 6, 8, 10	Dynamic triaxial compression test	H100/D50
Soft clay collected from the construction site of the railway station of metro line 1 in Tianjin at a depth of 10–15 m [27].	Cohesive force c = 1.64 kPa, Plastic limit wp = 20.7%, Liquid limit wL = 34.2%.	Re	28.2 (saturated)	ρ = 1.82 g/cm3 Void ratio e = 1.17	Respectively, −3, −7, −10, −15, −20, −25, −30~+20 10 h F/24 h TH	0, 1	Dynamic triaxial compression test	H80/D39.1
Loess from Jiuzhoutai town of Lanzhou, China [25].	Liquid limit = 17.4%, Plastic limit = 25.7%.	Re	16.5 (initial)	Dry unit weight = 17.8 kg/m3	−6~+15 −12~+15 12 h F/12 h TH	0, 3, 6, 9, 12	Triaxial compression creep test	H125/D62

.

In the limited achieved research works, rocks are the main research object and are tested for different engineering demands (Li, X.J. [8], Yang, X.R. & Jiang, A.N. [9], Zhang, F.R., Jiang, A.N. et al. [10], Yang, X.R., Jiang, A.N. et al. Wu, Y. [11], Zhu, J. [12]).

Some research focuses on the FTH effect coupled with the mechanical or chemical effect on the rocks (Chen, G.Q. et al. [13], Feng, X., Qin, N., Cui, L.Z. et al. [14], Li, J.L. et al. [15], Zhang, H.M. et al. [16], Song, Y.J. et al. [17]).

Compared with rocks, there are fewer studies on the creep behaviour of frozen–thawed soils and the engineering backgrounds of most of the experiments are based on deep geomechanics, where the exerted loading stress is high. Wang, J.W. et al. [18] studied the evolution law of the microstructures and creep deformation characteristics of deep clay after an FTH. It is shown that the angle between the long axis and the vertical direction of soil particles in undisturbed soil, remodeled soil, undisturbed frozen–thawed soil and remolded frozen–thawed soil is mainly in the range of 80°~100°, and the ratio between the long and short axes is mainly in the range of 1~2. Under different FTH states, the pore size is mainly distributed in the <1000 nm range. Additionally, the FTH effect increases the maximum strain of each level of creep and reduces the long-term strength. The direct creep test on frozen–thawed clay performed by Guan, S. [19] exhibits attenuated creep behavior, that is, strain tends to a limit value with time and there is no critical turning point on the isochronous curve cluster; however, this test result was due to the preset shearing load and did not exceed the long-term strength of clay. Within the attenuated creep state, creep strain becomes larger with more FTHs or a larger initial water content. The variation in creep strain after the first FTH is more obvious, while the more FTHs occur, the smaller the disparity among strain–time curves tends to be. Guan, S. utilized power function form, logarithmic function form and a five-component Kelvin model to express the strain–time curves of the creep test. Cao, H.L. [20] carried out constant-axial-pressure unloading confining pressure creep tests on the deep clay that was consolidated and then frozen–thawed. The FTH in the intermittent filling process of the core wall of high earth-rockfill dams will change the deformation behavior of gravelly soil. For 300 m high earth-rockfill dams, the deformation of gravelly soil is a crucial factor in predicting the long-term deformation and settlement after construction. Feng, M. et al. [21] used a modified Burgers model to describe deformation under the influence of FTHs.

However, for the railway subgrades, creep research of frozen–thawed soils is more desired because most subgrades are constructed on the soil stratum, and not on the rock stratum directly. Shi, Y.H. [5] carried out dynamic triaxial tests on silty clay in permafrost region along the Qinghai–Tibet Railway, obtaining the variation law of classical mechanical parameters, and established the frozen soil subgrade calculation model. With this model, the distribution of temperature and water in the subgrade with the FTH effect was studied and an assessment method for subgrade stability was determined. Xia, M.L. [22] performed creep tests on subgrade soil going through different FTHs. The results show that along the depth direction, the influence of water content and freezing temperature is stronger than the influence of the compaction degree; meanwhile, the variation in the compaction degree is predictable along the depth, while that of freezing temperature or water content along the depth is less easy to predict. With FTHs increasing, the creep deformation part of total deformation tends to decrease, the creep duration gradually decreases, and the duration interval between two FTHs becomes larger; that is, the inner structural damage of soil specimens becomes aggravated with more FTHs, so the upper limit of exterior load also decreases and specimens quickly creep to destruction. FTHs also weaken the radial constraint of specimens; while experiencing 5~7 FTHs, both the axial and radial deformation of soil specimens generally become stable. Gao, Z. [23] performed a triaxial creep test on the subgrade soil of a highway project in a seasonally frozen area after FTHs, and analyzed the evolution law of the axial and radial strain of soil samples experiencing different FTHs. The confining pressure was applied by controlling force and the axial load was applied by controlling displacement. The results show that the more FTHs occur, the less the final creep of the specimens, the more gradual the decrease in the final creep duration, and the more serious the damage to the specimen. Upon exerting the first-level load, the instantaneous strain values are relatively large and continue to increase as the left load level continues to be applied. With FTH increasing, the final creep duration becomes closer, and according to these tests, when the sample undergoes 8 FTHs, the creep duration changes little, so it is considered that after 8 FTHs, the degree of damage to the samples is basically stable. Wang, L.N. [24] researched the effects of freezing temperature, FTH, dynamic stress amplitude, confining pressure, vibration frequency, and humidity on the cumulative plastic deformation, rebound strain and corresponding elastic modulus of Qinghai–Tibet railway silty clay; the cumulative dynamic permanent deformation of the subgrade under the long-term train load was also researched. Referring to the Cam-Clay model, based on certain assumptions, a set of creep models was established. Based on this established cumulative permanent-strain model of frozen–thawed silty clay under train load, a creep subroutine interfaced with the visco-elastic–plastic constitutive model was developed, and a cumulative permanent deformation prediction model of the subgrade for train running was established; based on this, the influence of train axle load, speed, equivalent vibration order and FTH state on the development of the cumulative permanent deformation of the subgrade was discussed. Zhou, Z.W. et al. [25] studied the thermal cycling effect in mechanical behaviors through a series of triaxial compression, creep and stress relaxation tests conducted on frozen loess, after experiencing different FTHs, and they proposed that the strength, stiffness and viscosity properties of frozen loess weakened gradually with the increase in the freeze–thaw procedure until the cycle number reached the critical value for the steady state of the soil sample. The freezing low temperature ranging from −6 °C to −12 °C has no evident effect on the testing results of the same freeze–thaw numbers. The strength distribution loci of frozen loess in the equivalent stress–mean stress space exhibits unidirectional shrinkage over the range of the freeze–thaw process studied. The triaxial strength characteristics of frozen loess experienced different freeze–thaw numbers, and are more in agreement with the Ma’s failure criterion for frozen soils [26]. The stiffness properties (elastic modulus, elastoplastic coupled degree and damage degree) and the viscosity properties (relaxing stress, stable relaxation time, stable creep rate and viscoelastic modulus) of frozen loess both depend on the freeze–thaw process. Lei, H.Y. et al. [27] investigated the strain behaviors of soft clay undergoing FTHs and discussed the differences between the soft clay undergoing and not undergoing FTHs under dynamic cyclic loading; this was based on a series of dynamic triaxial tests with consideration of different influencing factors such as dynamic stress amplitudes, loading frequencies, freezing temperatures, and number of cyclic loadings. An empirical model for predicting the accumulative plastic strain was proposed and validated. Additionally, the microscopic structural characters of the artificially freeze–thawed soft clay were explored via mercury intrusion porosimetry (MIP) tests. The results show that the total accumulative plastic strain of specimens increases as the dynamic stress amplitude increases, while it decreases as the loading frequency increases. The specimens undergoing FTHs produce larger plastic strain, and the lower the freezing temperature, the larger the total accumulative plastic strain. Xu, X.T. et al. [28] conducted dynamic triaxial tests on frozen Qinghai–Tibet clay under different conditions to investigate the effect of freeze–thaw on the accumulative deformation behaviors. The development characteristics of axial permanent strain for frozen clay are sensitive to the freeze–thaw process. The permanent strain curve gradually decreases with the progress of the freeze–thaw process until it reaches the critical cycle number for the steady state. The microstructure variations induced by the thermodynamic cycles were discussed and analyzed to study the micro-mechanism for property degeneration. The experimental results demonstrate that the normalized accumulative strain curve rises with the increase in the number of FTHs, the stress amplitude and the number of loading cycles, but decreases as the confining pressure increases. Finally, a modified accumulation model was presented to simulate the development process of plastic strain under cyclic loading conditions, in which the freeze–thaw procedure, the loading cyclic number, the confining pressure and the dynamic stress amplitude were considered and supplemented.

Although there has already been some research about soils, Wang, M. et al. [29] performed an investigation and suggested that the confining pressure in the previous tests was usually 0.3~3.5 MPa, which is suitable for permafrost areas, but in seasonally frozen soil areas and areas where the working space of roadbed is less than 10 m, the maximum depth of frozen soil is about 10 m. Therefore, if a high confining pressure is also used in the consolidation process of the specimen, the stress behavior of the subgrade cannot be representative of the true depth, and the confining pressure in the frozen soil test should be about 0.1 MPa. Wang, M. et al. [29] also proposed that in the previous tests, the amplitude of dynamic stress was usually 0.4~5.6 MPa, which is much larger than realistic dynamic loads. To sum up, the inner stress that subgrade soils experience is not as high as deep geomechanics or rock mechanics, usually in the magnitude of several kilo-pascals. Therefore, this research carries out a series of low-stress triaxial compression creep tests on frozen–thawed silty clay. The creep properties are analyzed, and a mathematical model is constructed to describe the creep stress–strain–time relationship under the FTH effect.

2. Materials and Experimental Methods

2.1. Source of the Soil

The soil was initially undisturbed permafrost bored from a field between Yichun and Tieli in Heilongjiang Province, where a railway which will serve the Belt and Road Initiative will be built. To ensure the test results were available for analysis of the factors affecting creep behavior characteristics, the permafrost (shown in Figure 1a–c) was thawed and remolded into triaxial specimens (shown in Figure 1d,e).

2.2. Parameters of the Specimens

The grain gradation is illustrated in Figure 2. The optimistic moisture content and maximum dry density were 22.38% and 1.55 g/cm³, respectively, as exhibited in Figure 3, the soil was compacted into triaxial specimens with the compaction factor being 0.95, so the dry density of the specimens ultimately was 1.47 g/cm³. The size of the specimens was generally in 121 mm (Height) × 61.5 mm (Diameter); the accurate sizes were measured with the deviation in specimen manufacture, and are presented in

Table 2. Frozen–thawed cycle, preset confining pressure, geometrical size, mass and moisture content information of the remolded triaxial specimens.

No.	Freeze–Thaw Cycles (FTHs)	$\mathbf{Preset}\mathbf{Confining}\mathbf{Pressure}\left({\mathit{\sigma }}_3\right)$/kPa	Initial Average Diameter (D)/mm	Initial Height (H)/mm	Mass (M)/g	Moisture Content ($\mathit{W}$)/%
11	4	150	61.585	122.50	651.5	23.45
12	4	150	61.605	122.06	652.9	22.23
23	4	150	61.320	121.76	653.0	23.79
15	7	150	61.570	122.66	652.9	22.13
16	7	150	61.540	120.40	654.2	22.57
24	7	150	61.245	120.60	653.6	22.97
25	11	150	61.440	121.68	652.3	23.55
26	11	150	61.405	122.86	653.3	23.25
28	11	150	61.310	122.34	652.9	23.40

.

2.3. Test Conditions and Procedure

Before the triaxial creep tests (${\sigma }_2={\sigma }_3=constant$), a group of classical triaxial shear tests (${\sigma }_2={\sigma }_3=constant$) were performed as a reference to determine the upper limit and progression of the axial deviatoric stress load to exert in the creep tests. The ${\rho }_d$ and W were kept consistent between the creep and shear tests and the specimens for the shear tests were generally 80 mm (H) $\times$ 39.1 mm (D). The axial deviatoric principal stress (${\sigma }_1-{\sigma }_3$) at the failure point under each preset confining pressure $\left({\sigma }_3\right)$ is displayed in

Table 3. (${\sigma }_1-{\sigma }_3$) at the failure point under each ${\sigma }_3$ in the triaxial shear tests.

No.	${\mathit{\sigma }}_3$/kPa	$({\mathit{\sigma }}_1-{\mathit{\sigma }}_3)$/kPa
1	100	227.63
2	150	250.40
3	200	264.15

.

The external axial load levels of stepwise loading in the triaxial creep tests were set as displayed in Figure 4 in the form of deviatoric stress.

After being compacted, the specimens went through different freeze–thaw cycles (FTHs), as shown in Table 2, before finally participating in the creep tests. Each FTH consisted of 12 h thawing at +20 °C subsequent to 12 h freezing at −20 °C.

Considering the natural hydraulic state and environment that permafrost experiences and exists in, and to ensure the effect of FTHs expressed in the long-term creep tests, the specimens were exerted at the 1st level of axial load while still frozen, without consolidation or drainage. As time elapsed, the specimens were thawed and consolidated. For the left levels of loading, the specimens were consolidated and drained. Considering that the difference was unclear between “creep” and “consolidation” in these tests, this manuscript will discuss the variation in the total strain over time instead of taking out the “creep strain” part to discuss alone. The criterion for soil to be deformed in a stable state was axial deformation $\le$ 0.06 mm within 24 h. This criterion was not only set based on a criterion proposed by a book [30] about the mechanics of frozen soil, but also combined with the processes of practical tests. In the book, $\frac{\mathrm{d}\epsilon }{\mathrm{d}t}\le 0.0005{\mathrm{h}}^{-1}$, that is, for a 125 mm-tall specimen, 1.5 mm/24 h is set as the criterion of stable deformation, while through primary tests, 0.06 mm/24 h is already a relatively strict criterion.

3. Results

Figure 5, Figure 6, Figure 7, Figure 8, Figure 9, Figure 10, Figure 11, Figure 12 and Figure 13 display the axial creep strain–time curves (${\epsilon }_1\left(t\right)$ curves) of all the tested specimens, already converted from the tested curves of stepwise loading. The appearances of all the specimens when finishing the tests are displayed in Figure 14, Figure 15 and Figure 16.

4. Analysis

4.1. Analysis of the Isochronous Stress–Strain Curves

The isochronous stress $\left({\sigma }_1-{\sigma }_3\right)$-strain $\left({\epsilon }_1\right)$ curves are exhibited in Figure 17, Figure 18, Figure 19, Figure 20, Figure 21, Figure 22, Figure 23, Figure 24 and Figure 25. In these figures, the legend “stress at certain min” means “waiting time at each stress step: certain min”.

Generally, it is believed that the trend of the isochronous stress $\left(\sigma \right)$–strain $\left(\epsilon \right)$ curve with time elapsed is a criterion determining whether the creep behavior occurring is attenuation creep or non-attenuation creep [31]. If the stress–strain curve cluster tends to be scattered, the creep is considered to be non-attenuation creep, and vice versa. Additionally, the progression of being gathered to scattered is explained as soil physical–mechanical state changes. That is, upon experiencing the progression, the soil density and structure vary so much that a specimen will present quite a different mechanical response towards external stress. Additionally, the scattered tendency also makes the curve of stress–strain no longer linear, which indicates elasticity modulus changes. In philosophy, before and after the progression, the specimen is no longer the original specimen. Take steel constitutive research as a reference; this progression stage of soil is also defined as the yield stage and the beginning of this stage is defined as the yield stress of soil $\left({\sigma }_s\right)$. Apparently, when $\sigma \le {\sigma }_s$, the creep is attenuation creep, and vice versa. However, different scholars may determine different values of ${\sigma }_s$ (specifically, ${\left({\sigma }_1-{\sigma }_3\right)}_s$) because the recognization of the final gathered point on the curve is affected by the subjective will of a scholar in some way. So, here, instead, we list the yield stress interval of all nine of these creep specimens in Table 3.

Actually, the yield and plastic development of deformation can also be observed in the strain–time curves. In Figure 5, Figure 6, Figure 7, Figure 8, Figure 9, Figure 10, Figure 11, Figure 12 and Figure 13, the specimens display dramatically increasing deformation under the stress levels listed in

Table 4. The yield stress interval and dramatically increasing deformation stress level of the creep specimens.

No.	$\mathbf{Interval}\mathbf{of}{\mathit{\sigma }}_{\mathit{s}}$/kPa	Dramatically Increasing Deformation Stress Level/kPa
11	109.32~168.00	-
12	109.32~192.97	192.97, 217.94
23	75.96~168.00	168.00, 192.97
15	109.32~168.00	-
16	75.96~168.00	217.94
24	75.96~168.00	168.00
25	109.32~168.00	168.00
26	109.32~192.97	217.94
28	109.32~192.97	-

, too.

The physical meaning of the isochronous stress–strain curve cluster is that at different moments, soil maintains an extremely similar ability to resist deformation from external forces, exhibiting a similar linear or nonlinear constitutive relationship. However, in Figure 19, Figure 24 and Figure 25, the isochronous stress–strain curve of No. 23, 26 and 28 at 2048 min deviated largely from the other curves. Under the same stress, these specimens were deformed. The deviation means that no matter how large the axial deviatoric stress, with enough time elapsed, the soil will present a different ability to resist deformation to before. The deviation may result from two reasons. One is in the test process, whereby moisture migrates and drains away partly; the other is that even though there is no moisture running off, time is also a factor affecting the real-time constitutive relationship of soil, and there should be a critical moment for the qualitative change to occur in the soil. The time-induced change seems to not be inseparable from stress. At present, we have not carried out more research about time-induced qualitative change (for example, what is its induced factor?), but it deserves more research. Some inspiration can be found in the famous Asphalt Drip Test. It may be inferred that there is certain structural or compositional variation; for example, with enough time elapsing, soil particles disintegrate. After all, all the specimens for this manuscript experienced FTHs and the FTH effect was impressed.

From the isochronous stress–strain curve clusters, the elasticity modulus (E) of each specimen can be obtained, as shown in

Table 5. The elasticity modulus of each specimen and the mean elasticity modulus of a same FTH.

No.	FTH	E/MPa	$\mathbf{Mean}\mathit{E}\left(\overline{\mathit{E}}\right)$/MPa
11	4	7.452	8.140
12	4	11.008
23	4	5.959
15	7	9.518	7.038
16	7	4.613
24	7	6.983
25	11	8.174	8.791
26	11	9.369
28	11	8.831

.

For each FTH condition, three specimens were tested. The mean elasticity modulus of such soil specimen is affected by FTH in Formula (1) with good fitting accuracy.

$$\overline{E}\left(\mathrm{MPa}\right)=0.114{\mathrm{FTH}}^2-1.617\mathrm{FTH}+12.703$$

(1)

In general, from 4 FTHs to 11 FTHs, the mean elasticity modulus decreases first and then, increases. This is because with FTHs increasing, moisture and voids reconstitute in the process. At 7 FTHs, the structure becomes unstable, while at 11 FTHs, the structure develops into a firmer state. A similar conclusion can be found in the Introduction of this manuscript and it should be noted that different soils may vary on the FTH needed for them to be stable.

4.2. Analysis of Accumulative Time Needed for Stable Deformation When Going through Different FTHs

Figure 26, Figure 27 and Figure 28 show the accumulative time for the specimens to be stably deformed when exerting each preset axial deviatoric stress via stepwise loading.

Here, the rate of time that specimens took to deform, which agrees with $\le 0.06\mathrm{mm}/1440\min$, was extracted from each specimen. The time t(min) correlates with the axial deviatoric stress $\left({\sigma }_1-{\sigma }_3\right)$ (kPa) in approximate functions after going through 4, 7 and 11 FTHs, respectively, as in Equations (2)–(4), with good fitting accuracy, where $t_0$ is unit time (min) and $P_0$ is the standard atmosphere pressure = 1 atm.

$$4\mathrm{FTHs} t/t_0=4.46\times {\left(\frac{101\left({\sigma }_1-{\sigma }_3\right)}{P_0}\right)}^{1.55} \mathrm{Adj.} \mathrm{R}\text{-}\mathrm{Square} = 0.97677$$

(2)

$$7\mathrm{FTHs} t/t_0=1.57\times {\left(\frac{101\left({\sigma }_1-{\sigma }_3\right)}{P_0}\right)}^{1.77} \mathrm{Adj.} \mathrm{R}\text{-}\mathrm{Square} = 0.95688$$

(3)

$$11\mathrm{FTHs} t/t_0=0.07\times \left(\frac{101\left({\sigma }_1-{\sigma }_3\right)}{P_0}\right) { }^{2.39} \mathrm{Adj.} \mathrm{R}\text{-}\mathrm{Square} = 0.98496$$

(4)

The comparison of mean accumulative time–mean axial deviatoric stress under different FTHs is shown in Figure 29.

Clearly, until $\left({\sigma }_1-{\sigma }_3\right)=142.68 \mathrm{kPa}$, the mean accumulative time differs slightly under different FTHs, while after 142.68 kPa, the accumulative time needed until stable deformation differs largely. The increasing disparity reflects that if $\left({\sigma }_1-{\sigma }_3\right)\ge 142.68\mathrm{kPa}$, the more FTHs the specimens experienced, the more time they needed to be deformed stably under the same deviatoric stress. The possible cases are: 1. the more FTHs occurred, the larger the deformation. With extremely similar ability against deformation among the specimens, the larger the deformation, the longer the time. 2. The more FTHs occurred, the stronger the ability of a specimen to resist deformation. To reach almost the same deformation, the more difficult the deformation, the longer the time needed. 3. The more FTHs occurred, the larger the deformation and the stronger the ability. Further analysis of the cases is discussed in Section 4.4.

4.3. Analysis of the Effect of Axial Deviatoric Stress (ADS) on Stable Strain under 150 kPa Confining Pressure and 4, 7 and 11 FTHs

The correlation between mean ADS and mean stable strain is displayed in Figure 30 under varied FTHs.

From Figure 30, at the same mean ADS, the mean stable strain of FTH = 7 exceeds the other two FTH conditions and, in general, the mean stable strain of FTH = 4 exceeds FTH = 11. This phenomenon is in good agreement with the analysis of the elasticity modulus with varied FTH. As mentioned before, the modulus is 11 FTHs $>$ 4 FTHs $>$ 7 FTHs, while the mean stable strain is 11 FTHs $<$ 4 FTHs $<$ 7 FTHs under a same mean ADS. So, the elastic strain accounts for most of the total strain and significantly decides the extent of creep deformation. With the highest mean elasticity modulus, the specimens going through 11 FTHs are the last to reach the equivalent failure strain of 15%.

4.4. Analysis of the FTH Effect on Stable Strain under 150 kPa Confining Pressure and Preset Axial Deviatoric Stress

The correlation of mean stable strain–FTH exerted varied mean ADS is illustrated in Figure 31.

The mean stable strain of 7 FTHs is always higher than it of 4 FTHs and 11 FTHs, while that of 4 FTHs is not always higher than that of 11 FTHs, such as when ADS = 168.00 kPa. Though there were only three FTHs set in the experiment, as it may not be perfect to assess the effect of comprehensive FTHs on the stable strain, the correlations of 4, 7 and 11 FTHs are still fitted and provided.

$$\mathrm{When}\mathrm{ADS}=37.98\mathrm{kPa},{\epsilon }_{1stable} \left(\%\right)=0.3267+0.35534\mathrm{FTH}-0.0251 {\mathrm{FTH}}^2, \mathrm{Adj.} \mathrm{R}\text{-}\mathrm{Square} = 0.99241;$$

$$\mathrm{When}\mathrm{ADS}=75.96\mathrm{kPa},{\epsilon }_{1stable} \left(\%\right)=0.8976+0.20613\mathrm{FTH}-0.01544{\mathrm{FTH}}^2, \mathrm{Adj.} \mathrm{R}\text{-}\mathrm{Square} = 0.98485;$$

$$\mathrm{When}\mathrm{ADS}=109.32\mathrm{kPa},{\epsilon }_{1stable} \left(\%\right)=0.3234+0.46212\mathrm{FTH}-0.03241{\mathrm{FTH}}^2, \mathrm{Adj.} \mathrm{R}\text{-}\mathrm{Square} = 0.99273;$$

$$\mathrm{When}\mathrm{ADS}=142.68\mathrm{kPa},{\epsilon }_{1stable} \left(\%\right)=0.7557+0.48475\mathrm{FTH}-0.03477{\mathrm{FTH}}^2, \mathrm{Adj.} \mathrm{R}\text{-}\mathrm{Square} = 0.99117;$$

$$\mathrm{When}\mathrm{ADS}=168.00\mathrm{kPa},{\epsilon }_{1stable} \left(\%\right)=-1.1352+1.55465\mathrm{FTH}-0.10053 {\mathrm{FTH}}^2, \mathrm{Adj.} \mathrm{R}\text{-}\mathrm{Square} = 0.99495;$$

$$\mathrm{When}\mathrm{ADS}=192.97\mathrm{kPa},{\epsilon }_{1stable} \left(\%\right)=4.2339+2.41242\mathrm{FTH}-0.21061{\mathrm{FTH}}^2, \mathrm{Adj.} \mathrm{R}\text{-}\mathrm{Square} = 0.99868;$$

$$\mathrm{When}\mathrm{ADS}=217.94\mathrm{kPa},{\epsilon }_{1stable }\left(\%\right)=2.1945+3.78993\mathrm{FTH}-0.27284 {\mathrm{FTH}}^2, \mathrm{Adj.} \mathrm{R}\text{-}\mathrm{Square} = 0.99481.$$

Comparing Figure 30 and Figure 31, the stable strain increment between ADS = 142.68~168.00 kPa of 7 FTHs > 11 FTHs > 4 FTHs, and the real-time deformation modulus between ADS = 142.68~168 kPa of 4 FTHs > 11 FTHs > 7 FTHs. Though 4 FTHs corresponds to the smallest deformation increment, the specimen also has the largest modulus, so the time needed to be stable might be short and, in the same way, the circumstance of 7 FTHs is similar. For the 11-FTH specimen, though the deformation increment is medium, its modulus is also medium, resulting a longer time needed for stability. However, through the comparison, it can be inferred that the modulus is the principal factor because the deformation increment does not vary too much.

4.5. Analysis of Stress–Strain–Time Relationship under Varied FTHs

Through the literature investigation, we can see that there exist two defining modes of “creep strain”. Some scholars directly treat the total strain obtained in the creep tests as “creep strain”, while others divide the total strain into elastic, plastic and creep strain. Limited to the effective stress measurement in these tests, it is hard to rightly extract the creep strain part from the total strain. So, in the manuscript, the model to be established focuses on the total strain description. However, still, the total strain is expressed as two parts: elastic strain and visco-plastic creep strain. Usually, the elastic strain is considered equal to the instantaneous strain upon the exertion of stress.

Following the classical research on elastic strain, the Hooke element is adopted to denote the elastically deformed behavior of the soil. For the visco-plastic creep strain, to fit the progression tendency of the strain–time curves in Figure 5, Figure 6, Figure 7, Figure 8, Figure 9, Figure 10, Figure 11, Figure 12 and Figure 13, firstly, it is assumed that the strain will eventually develop to a stable limit as long as the time is long enough; that is, it is attenuation creep. Even though the soil is rapidly deformed at an early period, with abundant time, the strain will always gradually reach the convergence level, such as in the strain–time curve in Figure 6, Figure 7, Figure 9, Figure 10, Figure 11 and Figure 12 under the unstable stress level listed in Table 3. Secondly, this assumption refers to the fact that the model is at least required to describe the tendency. From the perspective of mathematics, for the classical models such as the Nishihara Model, it is the $-e^{-at}$ adopted as the base function that enables the model to describe the nonlinear behavior of soil. Mathematically, the $-e^{-at}$ is a convergence function tending to zero; in this way, the developed model based on the $-e^{-at}$ is able to describe the total strain–time relationships obtained in the attenuation creep tests.

Through exploration and many trials on the fitting of the strain–time curves, the arctan(t) function is finally selected as the base function. Here, the initially designed empirical strain–time function is

$${\epsilon }_1\left(\sigma ,t,\mathrm{FTH}\right)={\epsilon }_{elastic}+{\epsilon }_{elastic\text{-}plastic\text{-}creep}=\frac{{\sigma }_1-{\sigma }_3}{\overline{E}}+Aarctan\left(t/\left(C{t}_0\right)-D\right)$$

(5)

In this function, t is time (min). A is designed as the Amplification Coefficient because with a large enough positive t, the arctan(t) function will eventually tend to $\mathsf{\pi }/2$. The main idea is to evaluate the amplification of strain caused by FTHs and other potential factors. If there are some regular correlations between the reaction and inducing factors, the empirical function (5) would be of certain significance. D, a time value, is introduced to achieve the fitting. Because the strain–time curves in Figure 5, Figure 6, Figure 7, Figure 8, Figure 9, Figure 10, Figure 11, Figure 12 and Figure 13 are converted from stepwise loading curves, actually, only the curve stair under 37.98 kPa begins from 0 min and 0% strain; the other stairs of the curve begin from nonzero real coordinates. If the designed function lacks D, at 0 min, the nonzero ${\epsilon }_1\left(\sigma ,t\right)$ will not be suitable to describe the strain–time curve under 37.98 kPa. All the fitting results of the experiment data show that D is a minor value. The arctan(t) function changes rapidly at the beginning stage and even minor D will cause a considerable change in strain. However, the creep research focuses more on the stable strain, which is decided by A. Roughly, at the design stage, C is a Discount Coefficient for time, which shortens the developing time until convergence, enhancing the fitting accuracy.

4.5.1. The Mathematical Significance of the Empirical Stress–Strain–Time Model

In math, within the rectangular coordinate system, if the curve function is $y=f\left(x\right)$, and $f\left(x\right)$ processes the secondary derivative, from the geometrical significance of the derivative, $y^{\prime }=tan\alpha$, where $\alpha$ is the obliquity of the curve to the abscissa axis. There exists

$$\alpha ={arctany}^{\prime }$$

(6)

and its differential,

$$\mathrm{d}\alpha =\frac{1}{1+{y^{\prime }}^2}\mathrm{d}{y}^{\prime }=\frac{y^{″}}{1+{y^{\prime }}^2}\mathrm{d}x$$

(7)

For a curve, there exists $\mathrm{d}s=\sqrt{1+{y^{\prime }}^2}\mathrm{d}x$. If $\mathrm{d}\alpha$ is divided by $\mathrm{d}s$, and we take the absolute value, the formula of curvature k will be obtained.

$$k=\left|\frac{y^{″}}{{\left(1+{y^{\prime }}^2\right)}^{\frac{3}{2}}}\right|$$

(8)

If on the curve, $k\ne 0$ at a point M, $R=\frac{1}{k}$ is defined as the radius of curvature at the point M.

Theoretically, for the designed empirical model,

$${\epsilon }_{visco\text{-}plastic\text{-}creep}={\epsilon }_1-{\epsilon }_{elastic}=Aarctan\left(t/\left(C{t}_0\right)-D\right)$$

(9)

$$\frac{{\epsilon }_{visco\text{-}plastic\text{-}creep}}{A}=arctan\left(t/\left(C{t}_0\right)-D\right)$$

(10)

$$\mathrm{d}\alpha =\mathrm{d}\left(\frac{{\epsilon }_{visco\text{-}plastic\text{-}creep}}{A}\right)=\frac{1}{1+{\left(\frac{t}{C{t}_0}-D\right)}^2}\mathrm{d}\left(\frac{t}{C{t}_0}-D\right)=\frac{\frac{1}{C{t}_0}}{1+{\left(\frac{t}{C{t}_0}-D\right)}^2}\mathrm{d}t$$

(11)

$$\mathrm{d}s=\sqrt{1+{\left(\frac{t}{C{t}_0}-D\right)}^2}\mathrm{d}t$$

(12)

$$\left|k\right|=\frac{\left|\mathrm{d}\alpha \right|}{\left|\mathrm{d}s\right|}=\frac{\frac{1}{C{t}_0}}{{\left[1+{\left(\frac{t}{C{t}_0}-D\right)}^2\right]}^{\frac{3}{2}}}$$

(13)

When t = 0 min,

$$\left|k\right|=\frac{1}{C{t}_0}$$

(14)

$$R=C{t}_0$$

(15)

That is, C is actually the radius of curvature at the ordinate of the $\frac{{\epsilon }_{visco\text{-}plastic\text{-}creep}}{A}-t$ curve. $\frac{{\epsilon }_{visco\text{-}plastic\text{-}creep}}{A}$ is actually the obliquity of the function curve to the abscissa axis, and the function’s derivative is $\left(\frac{t}{C{t}_0}-D\right)$.

When t = $\infty$, $\left|k\right|$ = 0, agreeing well with the convergence hypothesis. Additionally, out of this, the values of C and D are actually of little sense for a long-enough duration, and the final amplitude of strain is more significant, which is decided by A.

4.5.2. The Fitting Results of All the Strain–Time Curves

Limited to the length of the manuscript, here, only the fitting image of strain–time of No. 23 is displayed, as shown in Figure 32.

4.5.3. The Specific Expression of A

The fitting values of A of all the specimens under all the ADS level with Equation (5) are listed in

Table 6. The fitting A values of all the specimens under all the ADS levels.

$\left({\mathit{\sigma }}_1-{\mathit{\sigma }}_3\right)$/kPa	11 A	12 A	23 A	15 A	16 A	24 A	25 A	26 A	28 A
37.98	0.4454	0.4325	0.5914	0.4404	0.6048	0.7808	0.4221	0.4264	0.5309
75.96	0.2941	0.2945	0.3429	0.2734	0.2851	0.2180	0.2425	0.2943	0.3003
109.32	0.0966	0.1351	0.1863	0.1502	0.2147	0.1676	0.1289	0.1824	0.1580
142.68	0.1777	0.0610	0.2188	0.1333	0.2108	0.1135	0.1732	0.1347	0.1531
168.00	0.3367	0.0971	1.9545	0.4795	0.5967	3.2099	3.1692	0.2157	0.2890
192.97	3.9808	2.5552	8.9364	5.1686	1.9565	8.0530	4.2001	0.6187	1.0293
217.94	7.0337	6.4043		8.6997	6.6189		6.3806	6.4529	3.3383

. The correlations between A and ADS are illustrated in Figure 33.

From Figure 33, as evidenced, the correlation between A and $\left({\sigma }_1-{\sigma }_3\right)$ of all the specimens varies little before 0.14268 MPa confining pressure. From 0.14268 MPa, specimens produce evident bifurcation of A. The A of most specimens sharply increases from 0.168 MPa but already differs from 0.14268 Mpa. So, from the tendency distinction, 0.14268 MPa actually should be the ${\sigma }_s$ of most specimens. Compared with other specimens, the A values of No. 23, 24 and 25 show different variation, which may result from the specimens’ preparation because the small disparity in the initial silt-to-clay ratio will affect the compressibility [32]. However, in general, excluding these three specimens, the specimens experiencing 11 FTHs (No. 26, 28) present a wider ADS range in the yield stage than those experiencing 4 (No. 11, 12) and 7 FTHs (No. 15, 16).

Evidently, the variation in A follows two modes with the increase in ADS. So, the corrections between A and ADS were studied, respectively, towards the two parts interrupted at 0.14268 MPa. Eventually, the function

$$A=e^{-a\left({\sigma }_1-{\sigma }_3\right)}$$

(16)

is found to be suitable to fit the correlations at low-level $\mathrm{ADS}\le 0.14268\mathrm{MPa}$, shown in Figure 34.

When FTH = 4, a is fitted to be 16.53169 with Adj. R-Square being 0.80397. When FTH = 7, a is fitted to be 15.18503 with Adj. R-Square being 0.81007. When FTH = 11, a is fitted to be 17.31641 with Adj. R-Square being 0.81504.

As mentioned in the analysis of isochronous stress–strain curves, the mean elasticity modulus of FTH = 4, 7 and 11 is 8.140, 7.038 and 8.791 MPa. Roughly calculated, it is found that the fitting value of a is twice the value of $\overline{E}$ (MPa). Considering some potential deviation in the fitting of A, C and D, it is finally decided that the ratio of 2 will be used in the empirical model. That is,

$${\epsilon }_1\left(\sigma ,t,\mathrm{FTH}\right)=\frac{{\sigma }_1-{\sigma }_3}{\overline{E}}+e^{-\frac{2\times {101}^2·\overline{E}}{{10}^6·P_0^2}\left({\sigma }_1-{\sigma }_3\right)}arctan\left(t/\left(C{t}_0\right)-D\right),\left({\sigma }_1-{\sigma }_3\right)\le {\left({\sigma }_1-{\sigma }_3\right)}_s$$

(17)

where $\overline{E}$ refers to Equation (1). The units used for the parameters in this model are

$$\%=\frac{\mathrm{MPa}}{\mathrm{MPa}}\times 100+e^{-\frac{\mathrm{MPa}}{{\mathrm{Pa}}^2}\left(\mathrm{MPa}\right)}arctan\left(\min /\min -\min \right)$$

(18)

While the coefficient of A increases sharply at high $\mathrm{ADS}>0.14268\mathrm{MPa}$, to reduce the error of the fitting results, the function mode of

$$ln\left(A\right)=k\left({\sigma }_1-{\sigma }_3\right)+b$$

(19)

is adopted, as shown in Figure 35.

$$\mathrm{When}\mathrm{FTH}=4,lnA=59.29\left(\frac{101\times 1000\left({\sigma }_1-{\sigma }_3\right)}{P_0}\right)-10.36, \mathrm{with} \mathrm{Adj.} \mathrm{R}\text{-}\mathrm{Square} \mathrm{being} 0.73998.$$

$$\mathrm{When}\mathrm{FTH}=7,lnA=55.01\left(\frac{101\times 1000\left({\sigma }_1-{\sigma }_3\right)}{P_0}\right)-9.50, \mathrm{with} \mathrm{Adj.} \mathrm{R}\text{-}\mathrm{Square} \mathrm{being} 0.81988.$$

$$\mathrm{When}\mathrm{FTH}=11,lnA=45.58\left(\frac{101\times 1000\left({\sigma }_1-{\sigma }_3\right)}{P_0}\right)-8.33, \mathrm{with} \mathrm{Adj.} \mathrm{R}\text{-}\mathrm{Square} \mathrm{being} 0.81988.$$

To present the effect of FTHs, the correlation among k, b and FTH is fitted and shows a linear characteristic, as shown in Figure 36.

So, the empirical model under high-level $\mathrm{ADS}>0.14268\mathrm{MPa}$ can be expressed as

$${\epsilon }_1\left(\sigma ,t,\mathrm{FTH}\right)=\frac{{\sigma }_1-{\sigma }_3}{\overline{E}}+Aarctan\left(t/\left(C{t}_0\right)-D\right),\left({\sigma }_1-{\sigma }_3\right)>{\left({\sigma }_1-{\sigma }_3\right)}_s$$

(20)

where $\overline{E}$ refers to Equation (1). The units used in this model are

$$\%=\frac{\mathrm{MPa}}{\mathrm{MPa}}\times 100+Aarctan\left(\min /\min -\min \right)$$

(21)

while A is determined by

$$ln\left(A\right)=k\left(\frac{101\times 1000\left({\sigma }_1-{\sigma }_3\right)}{P_0}\right)+b$$

(22)

$$k=-1.85 \mathrm{FTH}+66.49$$

(23)

$$b=0.29 \mathrm{FTH}-11.52$$

(24)

The unit of $\left({\sigma }_1-{\sigma }_3\right)$ used in the determination of A is MPa.

4.5.4. The Physical Significance of the Coefficients in the Stress–Strain–Time Model

If we differentiate the base empirical model in Equation (5) towards time t,

$$\frac{\mathrm{d}\epsilon }{\mathrm{d}t}=\dot{\epsilon }=\frac{\mathrm{d}\left(\frac{{\sigma }_1-{\sigma }_3}{\overline{E}}+Aarctan\left(t/\left(C{t}_0\right)-D\right)\right)}{\mathrm{d}t}=\frac{AC{t}_0}{C^2{t}_0^2+C^2{D}^2{t}_0^2+t^2-2CD{t}_{0}t}$$

(25)

With this expression, the strain rate of Specimen 23 is taken as an example to illustrate the comparison of fitting theoretical values and practical values, as in Figure 37, Figure 38 and Figure 39.

From the three figures, the strain rate calculated from the theoretical formula Equation (25) is always lower than the practical strain rate. So, the empirical model is not accurate enough at predicting the strain rate. However, in general, the empirical model can predict the tendency of variation in the strain rate with the increase in time.

5. Discussion

(1) Section 4.2, “Analysis of accumulative time needed for stable deformation going through different FTHs”, involves the analysis of accumulative time; however, the accumulative time is extracted from the test process via stepwise loading, with the time itself containing the stress history effect. The more densely compressed the soil in the former loading period, the less time it needs to be deformed stably when exerting the next load level. It seems that time is better-obtained from separate loading; however, given that the soil foundation does not tolerate one upper load, but instead, tolerates gradual loading in the engineering construction procedure, the analysis of Section 4.2 is still of some significance.

(2) Confining pressure is also a vital factor in the creep behavior of soil, and there is already some research on its effects on creep, while the coordinating effect with FTH still remains to be studied further. In the present study, we performed a preliminary study on the effect of confining pressure on the stable creep stage of soil that has experienced one FTH. We tested Specimens 6, 14 and 17 with the same moisture content, dry density and axial deviatoric stress as in Section 2, but under preset confining pressure of 100 kPa, 150 kPa and 200 kPa. The time needed until stable creep (TNTSC) under the series of axial deviatoric stress (ADS) is illustrated in Figure 40.

With a same confining pressure, the TNTSC increases with the increase in ADS. With the same ADS, the TNTSC decreases with the increase in confining pressure. These tendencies are coordinated with the research mentioned in the Introduction.

Additionally, the effect of confining pressure on the mean stable creep axial strain is also illustrated in Figure 41.

From Figure 41, we can see that with the same confining pressure, the stable strain increases with the increase in ADS. With a same ADS, the stable strain decreases with the increase in confining pressure, and the drop from approximately 100 kPa confining pressure to approximately 150 kPa is larger than the drop from approximately 150 kPa to 200 kPa when ADS = 75.96 and 109.32 kPa; this reflects the variation in the confining pressure effect on reducing creep deformation.

Because only three confining pressure levels were tested, the number of specimens participating in the tests is not abundant. However, some fitting results are still presented below, and more supplements in later research are expected.

For FTH = 1, when ADS = 37.98 kPa, ${\epsilon }_{1stable}=4.849-0.018\times 101{\sigma }_3/P_0$; when ADS = 75.96 kPa, ${\epsilon }_{1stable}=0.676-\left(-81.229\right)\times {0.975}^{101{\sigma }_3/P_0}$; when ADS = 109.32 kPa, ${\epsilon }_{1stable}=0.393-\left(-232.792\right)\times {0.972}^{101{\sigma }_3/P_0}$. ${\epsilon }_{1stable}$—Axial stable creep strain, %; ${\sigma }_3$—Confining pressure, kPa. ADS—Axial deviatoric stress (${\sigma }_1-{\sigma }_3$), kPa.

(3) In the research history of creep, the definition and determination of the long-term strength of rock and soil were once discussed in the 1980s, and there were four genres: the strength will reduce, will improve, will not change or will change uncertainly [33]. In general, the definition of long-term strength is that when the shear stress applied to the soil is less than a certain critical shear strength of the soil, the speed of the shear strain becomes smaller and smaller, and finally, tends to be stable, and the soil sample will not be damaged; when the shear stress applied to the soil is greater than a certain critical strength value of the soil, the shear strain exhibits non-attenuating creep and, finally, causes failure. The critical shear strength value is usually referred to as the long-term strength value. As time goes on, many scholars accept the theory of Tan Tjong-kie [34,35,36,37,38,39]: Stepwise loading of one specimen, which produces isochronous stress–strain curves and the stress value of the inflection points (from linear to nonlinear) on the curves, is just the long-term strength.

For most rock and soil, attenuation creep occurs when the stress load is lower than the long-term strength; when the stress exceeds the long-term strength, the creep will accelerate. The description of accelerated creep is not easy compared with attenuation creep. In addition, in practical engineering, the field within the range of attenuation creep tends to be controlled to ensure safety, so researchers usually preset the upper limit of loading lower than the long-term strength. However, the long-term strength is obtained from creep tests, and is unknown before the creep test is carried out. Therefore, in order to take the long-term strength as a reference for the test design, the predecessors studied the relationship between the long-term creep strength and the triaxial shear strength, expecting to obtain an approximate value of long-term strength indirectly from the triaxial shear strength; moreover, on determining the upper limit of loading based on the long-term strength, the limit can be divided into several steps with appropriate intervals to apply to the specimens.

According to the Casagrande test data, the ultimate long-term strength of clay with small creep is equivalent to 80% of the standard strength, while the ultimate long-term strength of clay with large creep is only 40% of the standard strength. A similar experimental study was conducted by Sun, F.W. et al. [40]. The clay particle content of the test soil was 32.2% and the plasticity index was 10. Keeping the water content unchanged, an unconsolidated and undrained direct shear creep test was carried out according to the standard test and a certain horizontal load lasted for 48 h and 96 h. The results show that the standard strength c is 90 kPa, and the 48 h long-term strength is 53 kPa (equivalent to 59% of the standard strength) and 48 kPa, corresponding to 96 h long-term strength, respectively. Yang, R.H. [41] mentioned in his dissertation that the long-term strength of frozen soil is only about 20% of the instantaneous ultimate strength. From the results of this manuscript, the yield stress, which is also the long-term strength, is 142.68 kPa/250.40 kPa = 56.98% of the shear strength.

(4) In practical engineering, it is not safe or reasonable to take 15% strain of the subgrade as the criterion to stop the operation of a railway system. Actually, considering the upper safety limit of subgrade settlement, it is better for the end condition of the creep test to be set at 5% strain [29]. To ensure safety, usability and durability, the subgrade will be designed to be able to tolerate the upper loads and self-weight, protected from accelerated creep, and usually after a period of operation, the subgrade will enter a stabilization state. So, actually, in this manuscript, the low ADS level results before the yield stress are more significant to practical engineering.

6. Conclusions

In general, from 4 FTHs to 11 FTHs, the mean elasticity modulus decreases first, and, then increases. This is because with FTHs increasing, moisture and voids reconstitute in the process. The modulus is 11 FTHs $>$ 4 FTHs $>$ 7 FTHs, while the mean stable strain is 11 FTHs $<$ 4 FTHs $<$ 7 FTHs under a same mean ADS. So, the elastic strain accounts for most of the total strain and significantly decides the extent of creep deformation. With the highest mean elasticity modulus, the specimens going through 11 FTHs are the last to reach the equivalent failure strain of 15%. However, considering the upper safety limit of subgrade settlement, it is better for the end condition of the creep test to be set at 5% strain. To ensure safety, usability and durability, the subgrade should be designed to be able to tolerate the upper loads and self-weight and protected from accelerated creep, that is, the inner stress should be controlled less than the yield stress of the soil;in this way, usually after a period of operation, the subgrade will enter a stabilization state. From the creep test results, the yield stress of the frozen–thawed silty clay is about 56.98% of the shear strength.