Factor’s Influence Analysis of Frost Heave during the Twin-Tunnel Construction Using Artificial Horizontal Ground Freezing Method

In order to analyze the inﬂuence of diﬀerent twin-tunnel parameters on the frost heave of the ground, six tunnel clear distances (0.25 D, 1.00 D, 1.50 D, 2.00 D, 2.50 D, and 3.00 d), three tunnel buried depths (6 m, 12 m, and 18 m), and two freezing orders (simultaneous freezing and sequential freezing) are selected to establish the two-dimensional numerical calculation models, using ABAQUS ﬁnite-element program and the user subroutine of frost heaving deformation, and considering the orthotropic deformation characteristics of frozen soil. Numerical simulation results show that the in-teraction of twin-tunnel freezing is weakened with the increase in tunnel clear distance. Correspondingly, the heaving displacement of the ground surface also decreases. Besides, the heaving displacement curve of the ground surface gradually changes from the unimodal N-type to bimodal M-type as the tunnel clear distance increases. However, the trough of the bimodal M-type curve gradually disappears as the tunnel buried depth increases. SPSS mathematical analysis results show that tunnel clear distance has the highest signiﬁcance, tunnel buried depth ranks the second, and freezing order ranks the last. It is important to select the appropriate tunnel clear distance during the design of the twin-tunnel construction using the artiﬁcial horizontal ground freezing method.


Introduction
e artificial horizontal ground freezing method has been widely applied to subway tunnel construction due to its temporary support capacity and sealing capacity. At present, many successful cases of artificial freezing tunnel construction are available for research. In China, the geological survey report of Tianhe passenger station of Guangzhou metro line 3 shows that the tunnel construction section passes through the sandy clay layer, which has poor compactness and rich water that often lead to quicksand and landslides [1]. To this end, the solution that combines the full section horizontal freezing method with CRD excavation method [2] was proposed and successfully applied to solve this engineering problem. e same situation also occurred in the tunnel project from Minzu University station to convention center station of Nanning rail transit line 1. e stratum of tunnel crossing is round gravel layer, which has large permeability that often leads to water flush accident, and the artificial horizontal freezing method works perfectly [3]. Hong Kong-Zhuhai-Macau Bridge, which is the key construction project in China, had attracted the attention of domestic and foreign engineers. Gongbei tunnel, an important part of its connecting line, also faced severe geological conditions. Many researchers [4][5][6][7] put forward different design schemes from the aspects of controlling the ground deformation, limiting the flow direction of groundwater and protecting the surrounding buildings. At last, the solution that combines the pipe curtain method with artificial freezing method [8] had been adopted with various efforts and made this tunnel engineering a good quality to resist natural disasters [9].
In view of this, the numerical simulation prediction accuracy of the ground frost heave has been further improved in the subway tunnel construction using artificial freezing technology based on the above successful cases. Konrad [10] first proposed the frost heave of segregation potential model and developed a two-dimensional numerical calculation program (SP-2D). is model revealed the frost heave law of the ground during the tunnel freezing construction of Kobe in Japan. Cheng and Zhang [11] used ADINA finite-element program to study the coupling law of the ground temperature and displacement based on a tunnel freezing project. Qiao and Tao [12] used FLAC finite-element program to study the ground displacement in the process of tunnel freezing construction. Yang et al. [13] used MARC finite-element program to establish a two-dimensional moisture-heat-stress model, which simulated the ground frost heaving displacement law in the connecting passage construction of Nanjing metro line 1. Yan and Fang [14] used GEO-SLOPE program to simulate the freezing and excavation process of connecting passage, based on the coupled thermo-solid theory considering the phase transition. Gao et al. [15] used COMSOL finite-element program to analyze the distribution of the temperature field, seepage field, and displacement field in the horizontal freezing engineering of connecting passage. Cai et al. [16] used ABAQUS finite-element program combining the subroutine of the frost heaving deformation to establish the thermal mechanical coupling numerical analysis method of the ground displacement during horizontal freezing of tunnel. e results of the above numerical calculation models are all consistent with the actual project, which also shows that the numerical simulation is an effective means to predict the ground change law in the tunnel freezing construction.
Based on the above point of view, it is meaningful to research the frost heaving phenomenon caused by the different twin-tunnel parameters during the horizontal freezing process of the twin tunnels through ABAQUS finite-element program, and the frost heaving law obtained by the numerical simulation will be of great guiding significance to the twintunnel design scheme using the artificial freezing technology.

Mathematical Model of Freezing Temperature Field.
Freezing temperature field of tunnel construction is a transient heat conduction problem with soil-water phase change. According to heat transfer theory [17,18], the governing partial differential equation for 2D freezing temperature can be written as where C * is called the equivalent volumetric heat capacity (kcal/(m 3 · C)), T is the soil temperature (°C), t is time (day), and k * is called the equivalent thermal conductivity (kcal/(m·d·°C)).
C * and k * are also given by where C f and C u denote the volumetric heat capacity (kcal/ (m 3 ·°C)) of frozen soil and unfrozen soil, and C f � ρ f ·c f and C u � ρ u ·c u , where ρ f and ρ u are the density of the frozen soil and unfrozen soil (kg/m 3 ), c f and c u are the mass specific heat of the frozen soil and unfrozen soil (kcal/(kg·°C)). k f and k u are, respectively, the thermal conductivity (kcal/(m·d·°C)) of the frozen soil and unfrozen soil. T d is the freezing temperature (°C) of soil; T r is the thawing temperature (°C) of soil. L is the latent heat of phase change for one cubic meter of soil (kcal/m 3 ).

Initial and Boundary Conditions of Freezing Temperature
Field. e initial condition of freezing temperature field is where T 0 is the initial temperature (°C) of soil, assuming that the formation temperature distribution is uniform. In the calculation model, each freeze pipe could be simplified into one point, so the boundary condition for each freeze pipe is expressed as where (x p , y p ) is the Cartesian coordinate of each point (freeze pipe) and T c (t) is the surface temperature (°C) of the freezing pipes. e boundary condition of soil in far field is given by On the surface of the ground, the boundary condition of convective heat transfer between atmosphere and soil is given by

Advances in Civil Engineering
where α 1 is the convective heat transfer coefficient between the atmosphere and soil (kcal/(m 2 ·d·°C)), T a is the atmospheric temperature (°C), and n 1 is the normal direction vector of the surface.

Orthotropic Frost Heave Deformation of Frozen Soil.
Existing studies show that [19][20][21][22] there is a big difference in the mechanical properties between the freezing front and the heat flow direction during the freezing of soil, and the frost heave deformation of frozen soil exhibits orthotropic characteristic. is means that the heat flow direction is the main principle direction, while the other two main directions are perpendicular to the heat flow direction, and the main frost heave deformation occurs along the heat flow direction.
When the temperature of soil reaches the freezing temperature (freezing point), the frost heave ratio of soil, which can be obtained from the indoor frozen soil test, is defined as the ratio of a volume increase to the initial volume [23]. erefore, the frost heave ratio can be regarded as the temperature-dependent instantaneous volumetric strain, which is symbolized by ε T f here, and the theoretical calculation method is shown as follows: where ε T f is the frost heave ratio of soil, Δh f is the longitudinal frost heave of sample, and h 0 is the original height of sample.
As shown in Figure 1, in the Cartesian coordinate system (x, y, z), a local coordinate system (1, 2, 3) is established along the heat flow direction, and 1, 2, and 3 represent the three principal axes of the local coordinate system, respectively.
When the soil temperature is colder than the freezing point, the total strain of the frozen soil consists of two parts, namely, the strain caused by load and the instantaneous volumetric strain caused by temperature. In the local coordinate system, the instantaneous thermal strain components can be obtained by the following equation: where η is an orthotropic parameter. When η � 1/3, the frost heave deformation of the frozen soil exhibits isotropic characteristic; when η � 1, the frost heave deformation of frozen soil exhibits unidirectional characteristic, and its direction is along the temperature gradient.
In the local coordinate system, the total strain component of frozen soil is calculated according to the elastic Hooke theorem, as follows: where E is elastic modulus of frozen soil and μ is Poisson's ratio of frozen soil. For the plane strain problem, there is ε 33 � ε 13 � ε 23 � 0, so that By substituting equation (10) into (9), the total strain component of the frozen soil in the local coordinate system can be expressed as e instantaneous thermal strain components can be written as

Advances in Civil Engineering
As shown in Figure 1, the heat flow angle θ is the angle between the local coordinate system and the Cartesian coordinate system. e instantaneous thermal strain components in the Cartesian coordinate system can be calculated by Substituting equation (12) into (13), we can get erefore, considering orthotropic frost heave deformation of frozen soil, the frost heave strain components (i.e., the instantaneous thermal strain components) can be obtained by equation (14) when the soil temperature is colder than the freezing point.

User Subroutine of Freeze-aw Deformation of Frozen
Soil.
e developed frost heave model considering orthotropic frost heave deformation of frozen soil is implemented by ABAQUS software. ABAQUS employs user-defined subroutines to extend the capabilities of the applications.
e developed frost heave model is programmed by three user subroutine interfaces (USDFLD, UEXPAN, and GETVRM). e flowchart of the user-defined subroutine is shown in Figure 2.
At each incremental step, GETVRM is called by USDFLD and the heat flow components of the integration point are used as state variables and transmitted to UEX-PAN. In the user subroutine interface UEXPAN, the heat flow angle in Figure 1 is calculated according to the heat flow components, and then the thermal strain components of each integration point are calculated according to equation (14). According to temperature determination of integral points, the frost heave deformation can be determined whether it occurs or not at the current incremental step.
In this paper, ABAQUS finite-element program and user subroutine are used to simulate the frost heave of the ground during the tunnel freezing construction. e analysis type is coupled temperature/displacement analysis. e main program of numerical simulation can simulate the whole construction process of tunnel freezing, excavation, and thawing.

Numerical Calculation Model
In order to verify the reliability of the numerical simulation method proposed in this paper, Cai [1] applied the method to the tunnel freezing project of Guangzhou Metro Line 3. e soil parameters input in ABAQUS finite-element software were obtained from the test report of Guangzhou artificial frozen soil. Eventually, the frost heaving displacement of the ground surface acquired from the numerical simulation is in good agreement with that of the field measurement. In view of this, using the same model and soil parameters to analyze the influence of tunnel parameters on the frost heaving displacement of the ground surface will produce certain engineering application value.

Implementation of the Frost Heave Model.
e length (xaxis) of the model is set as 80 m, and the width (y-axis) of the model is set as 50 m for eliminating the boundary effect according to modeling rules proposed by Cai et al. [24]. In addition, the cross section of each tunnel is circular with the outside diameter (D) of 6 m, and the centers of the two tunnels are on the same horizontal line. Correspondingly, the radius of horizontal freezing pipes layout circle (R d ) is set as 4m in order to form a frozen wall with an effective thickness of 2 m after freezing 95 days. In this paper, three tunnel buried depths (6 m, 12 m, and 18 m), six tunnel clear distances (0.25 D, 1.00 D, 1.50 D, 2.00 D, 2.50 D, and 3.00 D), and two freezing orders (simultaneous freezing and sequential freezing) are selected to establish two-dimensional models, and as a result, 18 models are established through a random combination between the tunnel clear distances and buried depths, which are shown in Table 1. In model 6-1, model 12-1, and model 18-1, the tunnel clear distance is 0.25 D, namely, 1.5 m, which is less than the effective thickness of frozen wall (2 m); therefore, three freezing pipes are shared in the middle of the two tunnels, and the total freezing pipes are 37, as shown in Figure 3(a). Meanwhile, the total freezing pipes are 40 in the remaining models, as shown in Figure 3 Finally, the ground frost heaving phenomenon is analyzed under the 18 conditions by ABAQUS finite-element program and the user subroutine of the frost heaving deformation based on the above two-dimensional geometric models. Due to the large number of models, only two representative meshing diagrams are shown here. Model 6-1 consists of 5384 elements shown in Figure 4(a), while model 12-3 consists of 6005 elements shown in Figure 4(b). e quadrilateral element (CPE4RT) with good stability is selected as the element type.

Calculation Method of ermodynamic Parameters.
Unfrozen water content is not only an important index to evaluate the characteristics of water migration in frozen soil, but also a commonly used parameter in the thermal calculation of frozen soil. Xu et al. [25] put forward a rapid method for determining unfrozen water content, who thinks that there is a dynamic equilibrium relation between the initial water content of soil and the freezing temperature, as shown in the following equation: where w 0 is the initial water content of soil (%), T d is the freezing temperature of soil (°C), and a and b are the empirical constants related to the soil properties.
Research [26] shows that the empirical constant b can be determined when the soil properties and the salt content of the soil are determined. Ulteriorly, Xu and Deng [27] put forward that the unfrozen water content can be calculated as follows: where w u is the water content of unfrozen soil (%), and T is the temperature of frozen soil (°C). e specific heat calculation equation of frozen soil is as follows: where C s is the specific heat of soil particles, it is taken as 0.18 kcal/(kg·°C) in this paper, C i is the specific heat of ice, it is taken as 0.5 kcal/(kg·°C) in this paper, and C w is the specific heat of water, it is taken as 1.0 kcal/(kg·°C) in this paper. e calculation equation of thermal conductivity of frozen soil is as follows: where Q is the heat quantity, ΔT/Δh is the temperature gradient (°C/m), ΔA is the area (m 2 ), and t is the time. e calculation equation of latent heat of frozen soil is as follows: where L w is the latent heat of phase change per unit mass of water, it is taken as 79.6 kcal/kg in this paper, and ρ d is the dry density of soil (kg/m 3 ). e calculation results of thermophysical parameters of soil are shown in Table 2 based on the above calculation equations.
e freezing temperature of soil is often colder than 0°C as the soil contains a certain amount of salt. Hou et al. [28] had effectively identified the variation characteristics of soil temperature, moisture, and salinity during freezing through a laboratory unidirectional freezing experiment. e salt in soil inhibited water migration, which obviously reduced the freezing point of soil. In this paper, the phase transition temperature range is considered to reflect this process. e freezing point of the soil selected is −0.5°C referring to the frozen soil test report [29]. e values of "Solidus Temp" and "Liquidus Temp" [30] are needed to be input in ABAQUS calculation program. erefore, the Solidus Temp is set as −0.5°C, and the Liquidus Temp is set as 0°C in the simulations. e mechanics parameters [27] of soil are shown in Table 3. T is the negative temperature of soil, and the frost heave ratio of soil is 1%. Besides, the freeze-thaw deformation characteristics of orthotropic soil are considered during the process of numerical calculation, where the deformation characteristic coefficient (η) is 0.9 [31].

Initial and Boundary Conditions.
e initial temperature of soil is set as 20°C according to the field measured ground temperature report [29]. e "ammonia-brine" cycle refrigeration system is widely adopted in practical projects,  which can stabilize the temperature of brine at −25°C [32]. erefore, the surface temperature of the freezing pipes is assumed to be −25°C. Besides, the temperature of frozen soil is much colder than that of air, which will inevitably lead to heat exchange between soil and air. erefore, the convective heat transfer coefficient between soil and air is assumed to be 175 kcal/(m 2 ·d·°C) [13].
For the bottom boundary of the calculation model, vertical and horizontal displacements are constrained. Horizontal displacement is not allowed on the left and right boundaries of calculation model. Only the upper boundary of the calculation model is free [24].

Effect of Tunnel Clear Distance
(1) Heaving Displacement Field of the Ground. After simultaneous freezing for 95 days, the distributions of the vertical displacement field are shown in Figure 5 with the tunnel buried depth of 6 m. Figure 5 shows that the vertical displacement field is symmetrically distributed along the center line of twin tunnels after simultaneous freezing for 95 days. In addition, the soil at the bottom of tunnel produces the downward frost heave displacement, while the soil at the top of tunnel produces the upward frost heave displacement. e vertical displacement isopleth of soil inside tunnel is distributed horizontally. It is caused by the fact that the soil inside tunnel is squeezed by the uniform frost heaving force in the movement of the frozen internal front. With the tunnel clear distances increasing from 0.25 D to 3.00 D, the maximum vertical displacements of the tunnel top soil are 60.5, 61.5, 56.3, 54.2, 52.2, and 50.9 mm, while the maximum vertical displacements of the tunnel bottom soil are −11.8, −11.2, −11.9, −12.9, −14.0, and −14.4 mm, respectively.
(2) Heaving Displacement of the Ground Surface. e heaving displacements of the ground surface with three tunnel buried depths are shown in Figure 6, which directly reflect the frost heaving phenomenon caused by the change of tunnel clear distance during freezing. Figure 6 shows that there are two types of the frost heave curves of the ground surface, which are unimodal N-type and bimodal M-type. e number of bimodal M-type curves decreases gradually with the increase of the tunnel buried depth under six tunnel clear distances; correspondingly, the number of the unimodal N-type curves increases gradually. Besides, the frost heave curves of the ground surface are all unimodal N-type under three tunnel buried depths as the tunnel clear distance is 0.25 D. Finally, the maximum vertical displacement of the ground surface decreases with the increase of the tunnel clear distance from 1.00 D to 3.00 D, which all appear at the center line of the tunnel.

Effect of Tunnel Buried Depth
(1) Heaving Displacement Field of the Ground. After simultaneous freezing for 95 days, the distributions of the vertical displacement field are shown in Figure 7 with the tunnel clear distance of 1.00 D. Figure 7 shows that they have the same laws that the vertical displacement field is symmetrically distributed along the center line of the twin tunnels after simultaneous freezing for 95 days. In addition, the soil at the bottom of the tunnel produced the downward frost heave displacement, while the soil at the top of the tunnel produced the upward frost heaving displacement. e vertical displacement isopleth of the soil inside the tunnel is distributed horizontally. It is caused by the fact that the soil inside the tunnel is squeezed by the uniform frost heaving force  line of twin tunnels. Besides, it can be seen from the other figures that the slope of the rising section of the vertical displacement curve of the ground surface decreases as the tunnel buried depth increases from 6 m to 18 m. It is easy to be found that the vertical displacement curve of the twin-tunnel freezing is composed of two normal curves [33].

Effect of Freezing Order
(1) Heaving Displacement Field of the Ground. e distributions of the vertical displacement field are shown in Figure 9 after sequential freezing for 190 days.         (2) Heaving Displacement of the Ground Surface. e heaving displacement curves of the ground surface are shown in Figure 10 after simultaneous freezing for 95 days and sequential freezing for 190 days, respectively. Figure 10 shows that the vertical displacement curve of the ground surface after simultaneous freezing for 95 days is included in the vertical displacement curve of the ground surface after sequential freezing for 190 days when the tunnel clear distance is 0.25 D and 1.00 D. However, two curves are basically coincident as the tunnel clear distance increases from 1.50 D to 3.00 D. Besides, the vertical displacement of the ground surface at the center line of the left line tunnel is slightly less than that of the right line tunnel. erefore, regardless the perspective of the construction period or the vertical displacement of the ground surface, simultaneous freezing is more economical and safer than sequential freezing [34].

Factor Significance Analysis
(1) Factor Influence Level. e vertical displacement of the ground surface is accumulated by the frost heave of the ground [35], which is one of the items that must be tested in the field. To this end, the vertical displacements of the ground surface corresponding to six tunnel clear distances, three tunnel buried depths, and two freezing orders have been imported to SPSS mathematical analysis program for the factor influence levels analysis, and the analysis results can be seen in Table 4.
It is assumed that the tunnel clear distance, tunnel buried depth, and freezing order are all significantly impacting the vertical displacement of the ground surface; namely, the value of significance is less than 0.05 [36]. Table 4 shows that the significance values of the tunnel clear distance, tunnel buried depth, and freezing order are 0, 0.069, and 0.410. Only the significance value of the tunnel clear distance is less than 0.05, which means that only tunnel clear distance has accepted the original hypothesis. Finally, the influence level of the three factors ranks as follows: tunnel clear distance > tunnel buried depth > freezing order.  Tables 5 and 6. It needs to be explained here that the number of the subfactors of the freezing order is less than 3, which is not satisfied to the condition of the post hoc test. Cai et al. [24] analyzed the freezing sequence of the uplink tunnel and downlink tunnel by model test and numerical simulation based on the actual twin-tunnel freezing project of the Shanghai Metro and found that the difference of the maximum ground surface frost heaving displacement in simultaneous freezing mode and sequential freezing mode is only 15%. Compared with the simultaneous freezing mode, the sequential freezing method extends the freezing time and increases the cost, but it can prevent frost heave damage and ensure project safety. So, the freezing order is not discussed below again. Table 5 shows that the influence level of the tunnel clear distance ranks as follows: 3.00 D < 2.50 D < 2.00 D < 1.50 D < 0.25 D < 1.00 D, which shows that the tunnel clear distance of less than or equal to 1.00 D should be avoided as far as possible during the design process of horizontal freezing engineering of twin tunnels. Table 6 shows that the influence level of the tunnel buried depth ranks as follows: 18 m < 12 m < 6 m, which shows that the frost heave of the ground surface caused by the artificial horizontal freezing method in the construction process of shallow buried twin tunnels cannot be ignored.

Conclusions
In this paper, six tunnel clear distances, three tunnel buried depths, and two freezing orders have been selected to study the ground deformation laws during the construction process of twin tunnels using the artificial horizontal freezing method, and the conclusions are as follows: (1) e ground displacement field is symmetrically distributed along the center line of the twin tunnels in the process of simultaneous freezing, while it is not symmetrically distributed in the process of sequential freezing, using ABAQUS finite-element program and the user subroutine of frost heave deformation.
(2) e vertical displacement curve of the ground surface gradually changes from unimodal N-type to bimodal M-type as the tunnel clear distance increases; however, the trough of the bimodal M-type curve gradually disappears as the tunnel buried depth increases.    (3) e influence level of the three factors ranks as follows: tunnel clear distance > tunnel buried depth > freezing order. It is also recommended that the tunnel clear distance of less than or equal to 1.00 D should be avoided as far as possible during the design process of the horizontal freezing engineering of twin tunnels, and it is necessary to strengthen the monitoring of the frost heave of the ground surface in the shallow buried twin tunnels. (4) Simultaneous freezing is more economical and safer than sequential freezing in the construction of twin tunnels.
Data Availability e datasets generated and analyzed during the current study are available from the corresponding author upon reasonable request.

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this paper.