#### Abstract

With the high demand for construction of tunnels in China’s severe cold regions, the problem of frost heaving has become an important factor that endangers tunnel safety. This paper attempts to investigate the effect of frost heave of cavity water that widely exists in the tunneling engineering on the tunnel stability. According to the actual deformation of the surrounding rock of the tunnel, the viscoelastic behavior is considered to the surrounding rock. On the premise of the elastic solution of stagnant water frost heave, the viscoelastic solution of frost heaving pressure is deduced by Laplace transform using the generalized Kelvin model based on the elastic-viscoelastic correspondence principle. The frost heaving force is analyzed through a case study with variations in the size of the cavity defect as well as the constitutive model parameters. It is concluded that the frost heaving force increases with the cavity defect size; over time, the frost heaving force gradually increases, but it will eventually stabilize. It is found that when the frost heaving force reaches a certain level, the surrounding rock with low strength or the lining with insufficient strength will crack, and the frost heaving force will not continue to increase.

#### 1. Introduction

The fast development of China’s infrastructure projects has seen a large number of tunnels being built in severe cold and permafrost regions. However, due to the low temperature, all year round and seasonal freeze-thaw cycles, tunnels in these areas are threatened by frost heave disasters. When the tunnel is underexcavated and the lining and the surrounding rock are not properly fitted, there will be cavity defects between the lining and the surrounding rock [1]. After the water accumulates in the cavity, it will freeze and expand in a low temperature environment. The frost heaving force will act on the surrounding rocks and the lining, causing the surrounding rocks and the lining to crack and generate cracks [2]. When the temperature rises, the frost heave melts into water and continues to infiltrate. Cracks in the surrounding rocks and linings expand, making the leakage water increase seriously [3, 4]. In addition, when the frost heaving force in the cavity acts on the lining, the tunnel lining will be damaged to varying degrees, and ice hanging will often occur, which not only invades the tunnel building boundaries but also endangers the traffic safety in the tunnel [5]. Therefore, the frost damage caused by the frost heave problem is an urgent problem to be solved in cold regions.

At present, there are three types of theories on the study of tunnel frost heave [6]: the first is the theory of the whole frost heave of the freeze-thaw circle, the second theory is the theory of accumulated water frost heave, and the third is the theory of water-bearing weathering layer. Klein [7] proposed a finite element method for time-dependent problems of frozen soils and investigated the mechanical behavior of temporary frozen Earth support systems for tunnels using this kind of finite element. Lai et al. [8] derived the viscoelastic analytical solution of the frost heaving force of cold zone tunnels in Laplace space using the Poynting–Thomson’s model using the theory of the whole frost heave of the freeze-thaw circle. Wang and Hu [2] deduced the expression of frost heaving force by using the theory of standing water frost heave and put forward suggestions for lining construction and design. The Norwegian Public Roads Administration conducted a research study on two highway tunnels that produced frost heaving as a result of larger ares with poor drainage conditions. Liu et al. [9] studied the frost heaving pressure by considering an elastic model of the surrounding rock of a circular tunnel. Based on the fracture mechanics and mesodamage theory, the relationship between the elastic modulus of the surrounding rock and the number of freeze-thaw cycles was derived. Based on this, a frost heave failure model considering the combination of rock elastic modulus and porosity is proposed. Yang et al. [10] proposed a frost heave analysis model that is coupled with water freezing, temperature field, and stress field to predict and control the frost heave. The numerical prediction results were compared with the field measurement results to study the effects of overlying soil layer thickness, frozen soil wall thickness, excavation radius, and brine temperature on frost heave. Some investigators combined laboratory data to simulate the frost heaving in the soil column through finite element program and conducted sensitivity test for frost heaving to evaluate the transient frost heaving that may occur on specific roads in specific winter under different conditions [11]. In a tunnel jacking project, Ooi et al. [12] monitored the displacement of the support wall by freezing method, and compared the measured results with the finite element analysis. Selvadurai APS [13, 14] and Song et al. [15] studied the interaction between underground pipelines and the frost heaving ettects in cold regions. Fan et al. [16] used the equivalent elastic coefficient method to calculate the frost heaving force. Deng et al. [17, 18] proposed a frost heave pressure model similar to gas pressure and gave the corresponding theoretical formula. Huang et al. [19] deduced the analytical solution of the frost heaving force of circular tunnels in cold regions on the premise of considering uneven frost heave of surrounding rocks. Gao et al. [20] regarded frost heaving as a result of the combined effects of overall frost heaving and local frost heaving of surrounding rocks and deduced the formula for calculating frost heaving force by using the calculation method of elasticity.

In this paper, we use the viscoelastic model to derive the expression of frost heaving force based on the theory of accumulated water frost heave. In addition, we explore the influencing factors of frost heaving force, which may serve as useful indicators for making leakage prevention measures in the tunneling engineering practice.

#### 2. Constitutive Model and Basic Assumptions

The deformation properties of the surrounding rock of the tunnel are very important for the design of the tunnel structure. In the calculation of frost heaving force, different rheological properties have different effects on the frost heaving force. If the tunnel surrounding rock is elastic, the deformation and stress adjustment of the surrounding rock are all completed at the moment of tunnel excavation. In other words, if the stress at the moment of excavation does not exceed the allowable strength value of the surrounding rock, the tunnel will be stable without collapse, and the development of time will not affect the stability of the cavern, and it will be stable. However, in practice, many tunnels will crack or even collapse during operation after excavation, which indicates that the stress and deformation of the tunnel increase with time. Therefore, it is not reasonable to adopt the elastic constitutive model. When it is considered in terms of elastoplasticity, it may reflect the nonlinear properties of rock and soil. By doing so, the calculated deformation value and the size of the plastic circle of the surrounding rock are more practical than that by elasticity solution. In the current literature, the plastic zone development is commonly empirically determined by the elastoplastic deformation value converging to a certain stable state, which is inconsistent with the actual deformation of the surrounding rock [21, 22]. We can take the viscoelastic constitutive model, calculate the viscoelastic deformation value of the surrounding rock, and judge whether it converges to a certain value over time, so as to determine whether the surrounding rock can reach the self-stabilization state. The difference between this method and the elastoplastic method is that the stability of the surrounding rock can be judged by whether or not the viscoelastic deformation rate approaches zero, and it can timely alert when the surrounding rock needs to take support and reinforcement measures. Therefore, in this paper, the viscoelastic constitutive model will be used for calculation, and the following assumptions will be made on the calculation model:(1)Ignore initial ground stress(2)Surrounding rocks are isotropic and homogeneous(3)Plane strain problem

#### 3. Calculation of Frost Heaving Force

According to the elastic-viscoelastic correspondence principle, if the viscoelastic solution of the frost heaving force of the water in the cavity behind the lining is required, its elastic solution must be obtained first [23].

Ignoring the frost heaving of the rock, the frost heaving of the cavity water behind the tunnel lining is mainly considered. Therefore, the simplified calculation model is illustrated in Figure 1.

Surrounding rocks, ice bodies, and linings are constrained to meet the deformation coordination conditions. When the compressive stiffness of any of the constraints is 0, the frost heaving force is 0. Let the compressive stiffness of the lining be *K*_{1}, and the compressive stiffness of the surrounding rock restraints forming the cavity be *K*_{2} and *K*_{3}, respectively.

Let the frost heaving force be *P*, then

The following relationship exists between the compressive stiffness and the elastic modulus:

After the water body in the cavity frost heaves, the volume increase of the water body is equal to the sum of the product of the elastic deformation of the surrounding rock restraint and the lining restraint, as shown in Figure 2, and

Generally, *α* = 0.09.

From equations (1) and (3), the frost heaving force *P* can be expressed as follows:

For the simplicity of computation, let the shape of the frost heave be an equilateral triangle with side length *a*, then

According to the uniaxial creep test results of frozen soil specimens: Under the conditions of low stress level and medium stress level, the entire creep process is stable and belongs to stable creep; when the stress level is high, the entire creep process is unstable and belongs to accelerated creep. It is found through comparison that the form of the stable creep curve of frozen soil is similar to that of the generalized Kelvin model, so the generalized Kelvin model (i.e., Kelvin–Voigt model) can be used to characterize the stable creep characteristics of frozen soil. The viscoelastic deformation of surrounding rocks and linings under frost heave can be calculated using the generalized Kelvin model (as shown in Figure 3), and its constitutive equation is as follows:

Let , then equation (6) is simplified as follows:

For a linear viscoelastic body in a complex stress state, the stress state at one point is resolved into the spherical stress tensor *σ*_{ii} and the deviator stress tensor *S*_{ij}. The general form of the differential rheological constitutive equation is:

In equation (7), *P*_{1}, *Q*_{1}, *P*_{2}, and *Q*_{2} are polynomials of the operator , that is, ,

When the stress level is not high, the volumetric deformation caused by the spherical tensor can be elastic, while the deviator strain caused by the deviator stress tensor is described by the generalized Kelvin rheological equation. In this regard, equation (7) can be rewritten as follows:

From equation (8):

For a linear elastic body in a complex stress state, its constitutive equation is as follows:

According to the corresponding principle, the *K* and *G* modulus of the viscoelastic body can be obtained by comparing equations (6) and (9).

According to the relationship between *E*, *μ* and *K*, *G*, we get the following equation:

Performing Laplace transform on the above formula, the image function of the elastic modulus *E* in Laplace space is as follows:

Based on this, the elastic modulus solution of the corresponding material is substituted into the elastic equation (5) as follows:

From equation (5), substituting the Laplace inverse transformation solution of the elastic modulus *E* of each material property into it, the viscoelastic solution under the generalized Kelvin model can be obtained. The general solution of *E* obtained from the inverse Laplace transform is as follows:

Substituting equation (15) into equation (5) and performing inverse Laplace transform, the viscoelastic solution of the frost heaving force is as follows:

Among them,where *A*_{i}, *B*_{i}, and *C*_{i} are material constants as defined after equation (6).

#### 4. A Case Study

The proposed solution is illustrated by analyzing an example tunnel section selected from the literature. The cold zone tunnel project is at an altitude of 3800 m with a burial depth of *H* = 100 mm. The annual average temperature is −3°C, and the temperature in the cave is −15°C at the coldest time. Based on the thermophysical parameters, the freezing depth is 3.5 m. Among them, Poisson’s ratio *µ*_{1} of the lining is 0.2, and Poisson’s ratio values of the surrounding rocks are, respectively, *μ*_{2} = 0.3 and *μ*_{3} = 0.33. The viscoelastic mechanical parameters of the surrounding rock and lining samples measured by tests are shown in Table 1 [24]. By substituting these parameters into the viscoelastic solution (16) for calculation, the change of the frost heaving force *P* acting on the lining with time *t* was obtained, and discussed in detail as follows.

##### 4.1. The Effect of Viscosity Coefficient *η* on Frost Heaving Force

The cavity size was taken as 200 mm, and the curve of frost heaving force with time under different viscosity coefficients was plotted. It can be seen from Figure 4 that at the same time when the viscosity coefficient is increased from 200 to 600, the frost heaving force gradually decreases. As time progresses, the influence of the viscosity coefficient on the frost heaving force gradually decreases, and it tends to the same stable frost heaving force value, and the frost heaving force is 1.5 MPa at this time. Within a certain time range, the viscosity coefficient of the material will reduce the frost heaving force. After that, the different viscosity coefficients will no longer affect the frost heaving force.

##### 4.2. The Effect of Elastic Modulus *E*_{h} of Surrounding Rock Series Spring on Frost Heaving Force

Similarly, the cavity size is 200 mm, and the curve of the frost heaving force with time for different elastic modulus *E*_{h} of different surrounding rock series springs is plotted. According to Figure 5, the curves almost coincide at different *E*_{h}, indicating that *E*_{h} has little effect on the frost heaving force.

Lai et al. [8] deduced the elastic analytical solution of frost heaving force under the condition of uniform frost heaving of surrounding rock, and according to the corresponding principle of elasticity-viscoelasticity, the viscoelastic analytical solution of frost heaving force was obtained by Laplace transform.(1)If the effect of freezing on the mechanical parameters of surrounding rock is not considered; that is, the mechanical parameters of surrounding rock will not change before and after freezing, then no frost heaving force will be exerted on the lining.(2)Under normal circumstances, the elastic modulus of the surrounding rock increases after freezing, and Poisson’s ratio decreases. At this time, it is equivalent to the expansion of the frozen and expanding surrounding rock pulling the lining outward. In Lai's calculation model, pressure and tension can be transmitted between the lining and the surrounding rock. In actual engineering, there are waterproof plates between the lining and the surrounding rock. The rock and lining tend to disengage, and the lining will not be affected by frost heaving force.(3)In Lai's example, the elastic modulus of the surrounding rock decreases after freezing, which is inconsistent with the actual engineering situation. The smaller the elastic modulus value of frost heaving surrounding rock, the greater the calculated value of frost heaving force. The calculated frost heaving force is 4.2 MPa, which is obviously greater than the actual measured value on the spot.

##### 4.3. The Effect of Elastic Modulus *E*_{k} of Surrounding Rock Parallel Spring on Frost Heaving Force

When the cavity size is 200 mm, at the same time, as the elastic modulus *E*_{k} of the parallel spring of the surrounding rock increases, the frost heaving force increases uniformly, and the growth rate gradually slows down, as shown in Figure 6. Compared with the elastic modulus *E*_{h} of the series spring, the elastic modulus *E*_{k} of the parallel spring has a relatively large effect on the frost heaving force.

##### 4.4. The Effect of Cavity Size on Frost Heaving Force

Taking the side length *a* of the cavity behind the lining as 100 mm, 200 mm, 300 mm, 400 mm, and 500 mm, and substituting the generalized Kelvin model parameters in Table 1 into the viscoelastic solution equation (16), we can obtain the freezing pressure at different times, as shown in Figure 7. In general, as time *t* increases, the frost heaving force *P* gradually increases, and after a certain time, it tends to a certain stable value indefinitely. As the cavity size increases, the frost heaving force also increases. When *a* = 100 mm, the frost heaving force reaches its maximum at about 40 d, and the frost heaving force is about 0.7 MPa; when *a* = 500 mm, the frost heaving force reaches a stable value at about 60 d. The frost heaving force is about 3.8 MPa. Therefore, the design of tunnel lining in the cold area should consider the influence of frost heaving force. The larger the initial defects in the tunnel, the greater the effect on frost heave. That is, the larger the size of the cavity behind the lining, the greater the frost heaving force, and with time, the frost heaving force gradually increases. When it increases to a certain value, the frost heaving force will converge to a constant value. At this time, when the surrounding rock is cracked or the lining construction quality cannot be guaranteed, the lining cracks and the water leakage phenomenon may occur.

Based on the above analysis, the following measures are proposed to prevent tunnel leakage caused by frost heave:(i)Improve the construction quality, strictly supervise the construction of the lining, and ensure that the quality is acceptable.(ii)Ensure a close fit between the lining and the surrounding rocks, minimize the tunnel defect voids, and leave no space for water accumulation.(iii)Pay attention to measurements of water prevention and drainage so that the possibility of leakage can be greatly reduced. For instance, the waterproofing and drainage measures of the tunnels have a good effect on alleviating the leakage of water.

For rock tunnels, if the size of the cavity defect behind the lining is about 0.1 m, the frost heaving force generated is 0.6 MPa. Even if the size reaches 0.2 m, the frost heaving force is only 1.4 MPa. According to on-site measured data, the maximum frost heaving force of Kunlun Mountain Tunnel is 0.3 MPa [25], the maximum frost heaving force of Chaoyang No. 2 Tunnel of Nenlin Railway in Northeast Forest Region is 0.2 MPa [26], the maximum frost heaving force of Dabanshan Tunnel is 1.03 MPa [27], and the maximum frost heaving force of Qingshashan Tunnel is 1.3 MPa [28]. It can be seen from the comparison that, within a reasonable range of physical and mechanical parameters of surrounding rock lining, the results of the calculation model in this paper are close to the measured values of the frost heaving forces of various tunnel projects. It can be seen that the calculation results of the analytical solution in this paper are in good agreement with the actual situation.

#### 5. Conclusion

(1)By comparing the elastic, elastoplastic, and viscoelastic constitutive models, a generalized Kelvin model (viscoelasticity) that is closely related to the time factor is selected for calculation: The viscoelastic deformation value converges to a certain value with time, and the deformation rate tends to be null. The generalized Kelvin model is used, and according to the correspondence principle of elasticity and viscoelasticity, the viscoelastic solution of frost heaving force on Laplace image space is obtained by Laplace transform based on the elastic solution of frost heaving force. The inverse Laplace transform is used to obtain the formula for the change of frost heaving force with time.(2)Through a case study, it was found that when the size of the cavity defect is constant, the frost heaving force gradually increases with the development of time, and eventually reaches a stable value; when the time is constant, the frost heaving force increases with the increase of the size of the cavity defect. Also, the occurrence of frost heaving is greatly related to the defects of the early stage of the tunnel, such as the poor backfilling behind the lining, existence of voids, and poor construction quality.(3)Freezing-thawing is an important issue in tunnel construction in high altitude and cold regions. The occurrence of freeze-thaw cycles is closely related to temperature, of which the effect on frost heaving cannot be ignored. Note that the influence of temperature factors is not covered in this paper, and it should be further studied.

#### Abbreviations

K_{1}, K_{2}, and K_{3}: | Compressive stiffness of the lining and the surrounding rock restraints forming the cavity |

E: | Elastic modulus |

G: | Shear modulus |

E_{h} and E_{k}: | The elastic modulus of the series spring and the parallel spring, respectively |

μ: | Poisson’s ratio |

η: | Viscosity coefficient |

σ: | Frost heaving stress of the viscoelastic material |

ε: | Frost heaving strain of the viscoelastic material |

: | Rate of change of frost heaving stress with time |

: | Rate of change of frost heaving strain with time |

s_{ii}: | Spherical stress tensor |

S_{ij}: | Deviator stress tensor |

P: | Frost heaving force |

: | Compression deformation of the lining and the two restraint surfaces of the surrounding rock under the effect of frost heave, respectively |

α: | Volume expansion rate when the water body changes from the liquid phase to the solid phase |

a: | Side length of the rock cavity |

V: | Volume in the liquid phase of the frost heaving body |

S, S_{1}, S_{2}, and S_{3}: | The side length/area of the lining and the surrounding rock restraint body, respectively |

P_{1}, Q_{1}, P_{2}, and Q_{2}: | Polynomials of the operator |

L: | Laplace transformation parameter. |

#### Data Availability

Data are available from the corresponding author upon request.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

This work was supported by the Natural Science Foundation of China (NSFC Grant no. 51879091) and the Fundamental Research Funds for the Central Universities (2018B01214) and they are greatly appreciated.