Analytical Solution of Three-Dimensional Heaving Displacement of Ground Surface during Tunnel Freezing Construction Based on Stochastic Medium Theory

Many studies indicate that the heaving displacement during tunnel freezing construction can be regarded as a plane strain problem. However, due to the fact that many of the frozen pipes are placed obliquely in practical engineering, it should be a three-dimensional problem. Therefore, considering stochastic medium theory, the analytical solution of three-dimensional heaving displacement of ground surface during tunnel freezing construction is established, and the analytical solution is applied to the project case in Beijing. The calculation results show that, considering the inclination angle between the freezing pipe and the center axis of the tunnel, the maximum heaving displacement of the ground surface occurs at the tunnel center. When the inclination angle is equal to 0 ° , the analytical solution of three-dimensional ground surface heaving displacement is in good agreement with that of the two-dimensional prediction model, which further veriﬁed the reliability of the established analytical solution of the three-dimensional frost heave prediction model.


Introduction
On the basis of the sandbox model test, the Polish scholar Litwiniszyn [1] proposed a stochastic medium theory to investigate the displacement of rock strata and ground surface during underground mining in the 1950s; in this theory, rock formation is regarded as a kind of "stochastic medium", and the deformation of rock formation and surface caused by underground mining is regarded as a stochastic process and the prediction of rock formation and surface deformation can be realized by using probability and statistics method. After the development and improvement of Chinese scholar Liu et al. [2,3], the application field of stochastic medium theory has developed from the analysis of surface deformation caused by underground mining to the analysis of surface deformation caused by open-pit mining, near surface excavation, and stratum drainage. Yang et al. [4] and Yang et al. [5] analyzed and discussed in detail the surface movement and deformation caused by urban tunnel construction by using stochastic medium theory, further expanding the application field of stochastic medium theory. Due to the deformation of the surrounding rock caused by shield construction of subway tunnel in soft stratum environment, surface subsidence and shield segment rupture problems occur frequently [6]. Shi et al. [7,8] considered the time-space development process of ground movement and deformation in tunnel excavation and, based on the stochastic medium theory, systematically researched the issues related to ground movement and deformation caused by mining method and shield method construction of urban tunnel. At present, the theory has become one of the practical methods to predict the surface deformation in urban metro tunnel engineering in China.
If ground frost heave induced by soil freezing is considered the inverse process of ground settlement [9] caused by soil excavation, then ground frost heave during tunnel freezing construction can also be predicted via the stochastic medium theory. Based on the stochastic medium theory and shear displacement method, Zeng and Huang [10] established a calculation model of soil deformation caused by double-O-tube shield tunnel construction. Yang et al. [4], Li [11], Tao et al. [12], and Zhou et al. [13] used stochastic medium theory to analyze the surface frost heaving deformation caused by the freezing method of subway tunnel construction and obtained some meaningful conclusions. However, the common deficiency is that freezing process is not considered in the analysis process. In the study of Li et al [14], considering the frost heaving force and frost heaving deformation of noncircular tunnel with lining support, the analytical solution of frost heaving force is obtained by the theory of complex function and continuity condition. In the study of Huang et al. [15], the analytical solution of circular tunnel frost heaving force is derived on the basis of the nonuniform frost heaving of tunnel surrounding rock. In the study of Cai [16][17][18] based on single-pipe freezing theory, flat-panel freezing theory, and calculation equation of frost heave ratio of soil, considering the formation process of the frozen wall and using stochastic medium theory, a twodimensional analytical solution of ground frost heave in horizontal freezing period of tunnel is established and the analytical solution is improved in this study by considering the freezing point of soil and assuming the constant surface temperature of the freezing pipe, the analytical solution methods for the radius of the freezing front and the outer radius of the expanded area are determined and applied to practical engineering cases, and then, the reliability of the analytical solution was verified by comparing with the field measurement results [19].
In the actual construction of the subway tunnel freezing method, the freezing pipe group is generally designed to be inclined and radial, which is not strictly horizontal freezing [20][21][22][23]. erefore, the research only on the three-dimensional frost heave deformation of the ground can reflect the actual situation of the project based on the analysis of the three-dimensional freezing temperature field of the ground.
In this paper, an analytical solution method of threedimensional frost heave deformation is established considering the inclined radial layout characteristics of the freezing pipe group in the actual tunnel construction. It is an improvement on the analytical solution of two-dimensional with the stochastic medium theory that the temperature field before the intersection of the frozen wall is calculated by single-pipe freezing theory and the temperature field after the intersection of the frozen wall is calculated by flat-panel freezing theory. Figure 1 shows that the length, width, and thickness of the microelement aredε, dξ, anddη, respectively. e distance between the center of the microelement and the surface is η. According to the stochastic medium theory, the vertical displacement of the soil at any point above the tunnel caused by the excavation of a single microelement can be obtained as

Ground Surface Settlement.
where r(z) is the convergence radius. e global coordinate system(X, Y, Z) is transformed into cylindrical coordinate system(R, θ, z). e equation for converting the two coordinate systems can be expressed as follows: erefore, in cylindrical coordinates, the vertical displacement of the upper soil caused by a single microelement is written as Based on the above assumptions, the volume of the microelement remains unchanged, that is, the volume of the soil in the whole excavation area remains unchanged, then where ε eR , ε eθ , and ε eZ is determined by the volumetric strain of the microelement in the cylindrical coordinate system, respectively. Substituting equations (3) and (4) into equation (5), the horizontal displacement U e can be obtained by the boundary conditions that R � 0and R ⟶ ± ∞, U e � 0.

Advances in Civil Engineering
In the Cartesian coordinate system, the horizontal displacement of the overlying stratum after the excavation of microelements is written as follows: Supposing r(z) is a proportional function of depth, the convergence radius can be obtained as where β is the main influence angle of the overlying soil layer and r(z) is the horizontal convergence radius in depth.
In actual engineering, excavation of the soil will make the soil around the tunnel contract inward. If the tunnel is excavated in the area Ω and the area after the excavation surface is contracted and becomes Ψ, then the size of the contraction deformation is Δ � Ω − Ψ.
In accordance with the superposition principle, the vertical displacement of the ground surface at any point will be the cumulative displacement of the upper soil sinking displacement, which is caused by excavation of each microunit body. It can be written as follows: e horizontal displacement of the upper soil settlement is expressed as follows: In accordance with equation (10), the three-dimensional model of the stochastic medium theory is utilized to predict the settlement of the upper soil after tunnel excavation. Two parameters need to be determined, namely, the main influence angle β of rock and soil layer and the deformation area Δ after tunnel excavation. Δ mainly affects the influence range of the soil layer, while β affects the shape of the settlement range.
From the above equations for calculating stratum settlement, there are many triple integral calculations involved, and the integral function of each triple integral is not integral. erefore, when solving these triple integrals, a numerical integration method should be adopted, which can be solved by using Mathematical and Maple software.

Analytical Solution for ree-Dimensional Ground Frost
Heave. Figure 2 shows that the angle α is located between the nth freezing pipe and the center axis of the tunnel. After the elapse of time t, as the cooling capacity is continuously transferred to the freezing pipe, the negative temperature in the freezing pipe will exchange heat with the surrounding soil layer, and the freezing will be diverged along the periphery of the freezing pipe, and the freezing radius of the freezing pipes gradually increases. e inner and outer freezing fronts of the frozen soil wall tend to rapidly smoothen. e frozen soil wall is assumed to expand uniformly from R(t) to R Δ (t) due to the frost heave characteristics of soil and the expanded area is Δ.

Advances in Civil Engineering
In accordance with the three-dimensional model of stochastic medium theory, when in the Cartesian coordinate system εξη, the vertical heaving displacement of the ground surface during the freezing time of tunnel can be expressed as follows: Since the integral region is a hollow cylinder, the space rectangular coordinate system is transformed four times in order to solve the integral equation. Figure 3(a) shows that the Cartesian coordinate system εξη(xyz) is converted into the Cartesian coordinate system ε 1 ξ 1 η 1 . e conversion equation of the two coordinate systems is Figure 3(b) shows that the coordinate system ε 1 ξ 1 η 1 is converted into the coordinate system ε 2 ξ 2 η 2 . And the center of the freezing pipe is connected with the center of the tunnel as an axis η 2 , the tangent of the freezing circle is taken as an axis ε 2 , and the axis ζ 2 coincides with the direction of ζ 1 . e angle between the axis η and the axis η 2 is set as c. e conversion equation of the two coordinate systems is expressed as follows: Figure 4(a) shows that the coordinate system ε 3 ξ 3 η 3 is obtained by rotating the coordinate system ε 2 ξ 2 η 2 around the axis ε 2 with angle α. And the conversion equation of the two coordinate systems can be written by Figure 4(b) shows that the integral region is a hollow cylinder, φ is the angle from oε 3 to om in the polar coordinate roφ, and the coordinate system ε 3 ξ 3 η 3 is converted into a cylindrical coordinate system rφζ for convenient calculation. e conversion equation of the two coordinate systems is as follows: From equations (12) to (16), the conversion equation from the Cartesian coordinate system εξη into the cylindrical coordinate system rϕζ can be expressed as follows: In accordance with the substitution equation of the triple integral, equation (12) may be written as follows:

Advances in Civil Engineering
e integral area is a hollow cylinder, the radius is from R(t) to R Δ (t), the length of the freezing pipe is l, and the vertical frost heave displacement of each point on the surface of the tunnel during the freezing period in the cylindrical coordinate system rφζ can be written as follows: where xyz are the coordinates of the points on the earth's surface in the Cartesian coordinate system (xyz). And X 1 Z 1 are the abscissa and vertical coordinates of the freezing pipe center in the Cartesian system (xyz), respectively. α is the angle between the freezing pipe and the central line of the tunnel, which is positive in elevation and negative in depression. β is the primary influence angle of ground heaving. rφζ are polar radius, polar angle, and vertical coordinates in cylindrical coordinates, respectively. And c is the angle between the connecting line and the axis η between the center of the freezing pipe and the center of the tunnel. l is the length of the freezing pipe.

Analytical Solution before the Closure of the Frozen Wall.
Suppose that a circle of freezing pipes is arranged around the tunnel, and the freezing pipes are arranged uniformly. Freezing pipes are arranged horizontally and the soil properties are the same. If the number of frozen pipes is n, in accordance with the superposition principle, the displacement of ground frost heave caused by tunnel freezing before the closure of the frozen wall is equal to the sum of ground frost heave displacement caused by n single-pipe freezing. en the freezing radius of all freezing pipes at t time is R(t), and the frost heave radius is R Δ (t).
e serial number of freezing pipes is recorded as 1, 2, ..., n. In the Cartesian coordinate system εξη(xyz), the coordinate of the center of the nth freezing pipe is defined as(x i , z i , 0), which can be obtained as where h is the depth between the tunnel center to the ground surface and R d is the radius of the freezing pipe layout circle. Before calculating the frost heave of the ground surface caused by the ith freezing pipe, the Cartesian coordinate system needs to be converted to the cylindrical coordinate system rφζ. e conversion equation is ε � r cos ϕ cos c +(ζ sin α + r sin ϕ cos α)sin c + X i , (21) ξ � r sin ϕ sin α − ζ cos α, Substituting equation (20) into equations (21)-(23), respectively, then η � r cos ϕ sin c − (ζ sin α + r sin ϕ cos α)sin c + Z i . (26) Hence,

Advances in Civil Engineering
Substituting equation (27) into equations (23)-(25), respectively, then ε � r cos ϕ cos i 2π n +(ζ sin α + r sin ϕ cos α)sin i 2π n From equation (19), under the condition of space problem, the vertical frost heave displacement of each point on the ground surface caused by the ith freezing pipe can be expressed as Substituting equations (28)-(30) into equation (31), respectively, then In accordance with the superposition principle, the displacement of ground frost heave before the closure of the frozen wall during tunnel freezing is equal to the sum of heaving displacement caused by n single pipes, which can be obtained by

Advances in Civil Engineering
Similarly, the horizontal frost heave displacement of the ground surface before the closure of frozen wall during tunnel freezing construction can be expressed by

Analytical Solution after the Closure of the Frozen Wall.
As the frozen soil columns gradually enlarge with prolonged freezing time, the adjacent frozen soil columns interconnect to form a continuous frozen soil wall around the tunnel. e inner and outer freezing fronts of the frozen soil wall tend to rapidly smoothen. Figure 5 shows that the frozen soil wall is assumed to expand uniformly from R(t) to R Δ (t) due to the frost heave characteristics of soil and the expanded area is Δ. Ground frost heave is induced by the expansion effect of the entire frozen wall after the closure of the frozen wall. e Cartesian coordinate system εξη is converted to the cylindrical coordinate system rθζ and θ is the angle of polar coordinate ro ' θ. e axis ξ is perpendicular to the outside of the paper, and the axis ζ is in the same direction as the axis ξ.
en the conversion equation is In accordance with the substitution equation of the triple integral, equation (12) may be written in the cylindrical coordinate system rφζ as When the angle between the freezing pipe and the horizontal direction is α, the integral region is approximately a hollow frustum, and the vertical heaving displacement of the ground surface is written by 8 Advances in Civil Engineering (41)

Solution for the Radius of the Freezing Front.
According to the paper written by Cai and Jiang et al. [17,24,25], it can be known that the single-pipe freezing theory can be applied to two-dimensional heat conduction problems, while the flat-panel freezing theory is applied to one-dimensional semi-infinite heat conduction problems. erefore, the radius of the freezing front is the radius of a single frozen soil column before the closure of the frozen wall, which can be solved using the single-pipe freezing theory.
After the closure of the frozen wall, the radii of the freezing front are the inner and outer radii of the frozen wall, which can be solved using the flat-panel freezing theory. Figure 6 shows that the temperature field can be divided into frozen and unfrozen areas according to the location of the freezing front. e temperature of the soil in the frozen area is T f and that in the unfrozen area is T u , both of which are functions of time t and radial coordinate r. e radius of the frozen soil column is r(t), and soil temperature at the freezing front is freezing point T d . e initial temperature of soil is T 0 and the outer radius of the freezing pipe is r 0 . e surface temperature of the freezing pipe is assumed constant and its value is T c .

Single-Pipe Freezing eory.
For this two-dimensional heat conduction problem, the differential equation of heat transmission in frozen and unfrozen areas can be written as follows: where α f and α u are the thermal diffusivity of frozen and unfrozen soils, respectively: where k f and k u are the thermal conductivity of frozen and unfrozen soils, respectively. c f and c u are the specific heat of frozen and unfrozen soils. ρ f and ρ u are the density of frozen and unfrozen soils. e initial conditions for the differential equations are as follows: e boundary conditions of the differential equations are as follows: Advances in Civil Engineering e surface temperature of the freezing pipe is easily measured using a temperature sensor. erefore, the surface temperature of the freezing pipe is assumed constant in this study. Equation (47) is selected as the boundary condition of the differential equation. e heat balance equation at the phase transition boundary is given by where L is the volumetric latent heat of soil, which is given by where L w is the latent heat of water, ρ d is the dry density of the soil, w 0 is the initial water content in the soil, and w u is the unfrozen water content in the frozen soil. Combining the above initial and boundary conditions to solve the partial differential equations, the temperature field distributions in frozen and unfrozen areas can be obtained as follows: where E i (x) is the exponential integral function, which is given by

Substituting equations (51) and (52) into equation (49) yields
where A is an undetermined constant that varies with time, which can be In accordance with single-tube freezing theory, the radius of the freezing front (i.e., the radius of the frozen soil column) may be solved by using equations (54) and (55) before the closure of the frozen wall when the surface temperature of the freezing pipe is constant. e adjacent frozen soil columns begin to interconnect, which indicates that the frozen wall has just closed, and the closure time of the frozen wall can be calculated by the following equation: where t i is the closure time of the frozen wall and l is the space between adjacent freezing pipes. Figure 7 shows that the flat panel is assumed to generate heat exchange only with the soil on the right side, the temperature of the flat panel is T C ′ , and the temperature of the soil in the freezing and unfrozen areas is T f and T u ′ , respectively, both of which are functions of time t and coordinate x. e distance between the freezing front and the flat panel is X(t), and the temperature of soil at the freezing front is the freezing point T d . e initial temperature of soil is T 0 .

Flat-Panel Freezing eory.
For the flat-panel freezing problem, when 0 ≤ x ≤ X(t), the differential equations of the frozen area can be written as (57) Similarly, when X(t) ≤ x < ∞., the differential equation is (58) e initial conditions for the differential equations are as follows: e boundary conditions of the differential equations are as follows: x � X(t), In addition to the above boundary conditions, the heat balance equation at the phase transition boundary X(t) can be expressed as Unfrozen area Frozen area T f r 0 Freezing pipe wall Freezing front T c r θ Figure 6: Schematic diagram of single tube freezing.
Combining the above initial and boundary conditions to solve the partial differential equations, the temperature field distributions in the frozen area can be obtained as follows (0 ≤ x ≤ X(t)): e temperature field distributions in the unfrozen area can be expressed as follows (X(t) ≤ x < ∞): where Φ(y) is Gaussian integral function, which is given by Substituting equations (64) and (65) into equation (63) yields where B is a constant to be determined, and the following condition exists: erefore, the flat panel is regarded as the freezing pipe circle after the closure of the frozen wall. e radius of the freezing fronts is calculated as follows: where R(d) is the radius of the freezing pipe circle, R(t) is the radius of the outer freezing front, and Rc(t) is the radius of the inner freezing front. e thickness of the frozen wall E(t) is obtained by In accordance with the flat-panel freezing theory, the radii of the freezing fronts after the closure of the frozen wall can be solved by using equations (67) to (70) when the thermophysical parameters of the frozen and unfrozen soils are determined. It is generally considered that the flat-pane cooling temperature T c ′ is equal to the freezing tube wall temperature T c .

Solution for the Outer Radius of the Expanded Area.
e frost heave rate in engineering is generally expressed as where ε f is the frost heave ratio of the soil, Δh f is the longitudinal frost heave of the sample, and h is the original height of the sample. Liu [26] thinks that the frost heave at a certain point during the freezing process of the soil is the integral of the frost heave ratio of the frozen soil layer along the thickness of the frozen layer, which is not directly related to the frost heave force. According to its principle, the frozen wall expands evenly outward, then Many experimental studies on frost heave of soils show that the frost heave of soils is restrained by loads, and the frost heave ratio of soils decreases with the increase of overlying loads in exponential function.
where ε f0 is the frost heave ratio of soil without load, P is the load, kPa, and b is empirical constant, kPa − 1 . Substituting equation (74) into equation (73) yields

Case Analysis
e return line tunnel of a passenger transport station in Beijing, China, with a length of 40.0 m, a width of 6.0 m, and a buried depth of 13.0 m, is a shallow-buried tunnel with a large cross section. Given that tunnel will pass beneath existing pipelines, artificial ground freezing and undermining methods, instead of the open cut method, are selected for tunnel construction. Figure 8 shows that the sand layer is frozen and reinforced, and the lower clay layer is closed, and the frozen wall is semicircular arched along the tunnel [27].
After the active freezing period of the tunnel for 35 days, a nearly semicircular arched horizontal frozen wall with a thickness of 1.2-1.6 m and an arch span of 6.0 m was finally formed. e frozen wall was a nonclosed structure, and the average temperature of the frozen wall was − 10°C. e layout of the freezing pipes in the tunnel is shown in Figure 8.

Advances in Civil Engineering
For the convenience of calculation, the simplified thermal and physical parameters of the tunnel are shown in Tables 1 and 2. e tangent value of the primary influence angle of the tunnel istan β � 0.8, and the convergence radius can be obtained by erefore, − 20 m ≤ x ≤ 20 m and − 10 m ≤ y ≤ 10 m are taken as the calculation area. In accordance with equation (54), the thermal diffusivity of frozen and unfrozen soils may be obtained by In accordance with equation (50), the volumetric latent heat of soil is e axis x is perpendicular to the freezing direction and the axis y is the direction along the freezing wall. And the frost heave radius and related parameters of different freezing time were introduced into the three-dimensional frost heave prediction model. Figure 9 shows that the inclination angle α exists between the freezing pipe and the center axis of the tunnel. Because of the horizontal freezing, the inclination angle is equal to 0°, and the length of the frozen area is 20 m. In accordance with the Legendre-Gauss integration method, the numerical programs for the displacement of the ground surface are compiled using Maple software, and the results are presented in MATLAB mathematical software. e heaving displacement of the ground surface due to freezing for 30 days, 45 days, 60 days, 75 days, and 90 days is shown in Figures 10-14, respectively.
However, if the inclination angle α is not 0°, three-dimensional frost heave prediction model will be more complex in the same case. Since the inclination angle given in reference [28] is 4 ∘ ∼ 14 ∘ , for the convenience of    Figures 15-19, respectively. As shown in Figures 10-14, when the inclination angle α between the freezing pipe and the center axis of the tunnel is 0°, the heaving displacement of the ground surface along the x-axis is in agreement with the results of the two-dimensional heaving displacement, which gradually decreases with an increase in distance from the tunnel center. Along the y-axis direction, that is, the tunnel axis direction, the heaving displacement of the ground surface in the middle part rarely changes, basically keeps horizontal. And the heaving displacement at both ends is obviously reduced. Since the length of the frozen pipe is 20 meters, the corresponding heaving displacement is the same, which can be regarded as the result of the translation of the same section. However, in two ends of the freezing pipe, the temperature of the freezing pipe is affected by the surrounding soil, the frost heave radius changes and the heaving displacement decreases.
While the inclination angle is equal to 10°, the maximum heaving displacement of the ground surface does not occur at the central axis of the tunnel. Due to the influence of the inclination angle of the freezing pipe, one end of the freezing pipe is relatively close to the ground surface, the frost heave phenomenon will be more obvious, and the heaving displacement of the ground surface is larger than that of the two-dimensional prediction model.       14 Advances in Civil Engineering