#### Abstract

This paper addresses the difficult closure of a frozen wall in a coal mine shaft due to excessive seepage velocity in an aquifer when the aquifer is penetrated via the artificial freezing method. Based on hydrothermal coupling theory and considering the effect of decreased absolute porosity on seepage during the freezing process, a mathematical model of hydrothermal full-parameter coupling with a phase change is created. A shaft is used as a prototype, and COMSOL multiphysics finite element software is employed to perform a numerical simulation of the shaft freezing process at various stratum seepage velocities. The numerical simulation results are verified via a comparison with field measurement data. Based on the numerical simulation results, the impact of various underground water seepage velocities on the artificial frozen wall formation process with the seepage-temperature field coupling effect is analysed. Based on the analysis results, the recommended principles of the optimization design for a freezing plan are described as follows: first, the downstream area is closed to enable the water insulation effect, and second, the closure of the upstream area is expedited to reduce the total closure time of a frozen wall.

#### 1. Introduction

Due to core advantages, such as water insulation and temporary support, the artificial freezing method has been extensively applied in geotechnical engineering and become a common method for penetrating an aquifer in a coal mine shaft [1–3]. However, frozen wall closure failure due to excessive aquifer seepage velocity is common in engineering practice [4].

Many researchers have conducted extensive investigation of the freezing problem with high flow velocity. First, based on the Harlan [5] hydrothermal coupling model, a large amount of extension research was conducted. Li et al. [6], Liu et al. [7] and Kurylyk and Watanabe [8] elaborated differences among numerous studies. Vitel et al. [9] suggested that only two studies [9, 10] have performed experimental validation for high flow velocity. Second, based on hydrothermal coupling theory, Yang and Pi [11] created a mathematical model for single-hole freezing front development considering underground water flow and performed a comparative analysis of numerous influencing factors of single-pipe freezing via numeric calculation. Liu et al. [12] performed a qualitative analysis of the effect of horizontal underground water flow on horizontal freezing in a rectangular tunnel. Zhou and Lan [13] briefly analysed the variation of seepage and temperature fields during freezing via the numerical simulation of the artificial freezing process in a seepage stratum. Based on the proposed mathematical model, Marwan et al. [14, 15] employed an ant colony algorithm to optimize the tunnel freezing-hole layout in seepage conditions. However, the theoretical models of this numerical analysis were not validated by experiments, and the relationship between permeability variation and soil state during freezing was not explained. These numerical examples failed to summarize the common law for the effect of a seepage velocity increase on the formation of a shaft artificial frozen wall based on a detailed freezing plan.

Therefore, based on hydrothermal coupling theory and considering the effect of an absolute porosity decrease on seepage during freezing, a mathematical model for hydrothermal full-parameter coupling with a phase change is created. The combined programming of MATLAB and COMSOL is employed to provide a numerical solution for a coupling mathematical model. A mine shaft double-circle-piped freezing plan is employed as a prototype to perform a numerical simulation and investigate the effect of underground water seepage velocity on frozen wall formation. The focus is the analysis of the effect of the underground water seepage velocity on the seepage field, temperature field, frozen wall closure time, and maximum thickness during freezing. An optimization principle for a high-flow-velocity shaft freezing plan is proposed, which has a significant theoretical and engineering value for solving the problem of the difficult closure of a shaft frozen wall in a loose aquifer with high flow velocity.

In Section 1, the mathematical model for hydrothermal full-parameter coupling with a phase change is created. In Section 2, a mine shaft double-circle-piped freezing plan is employed as a prototype; a combined programming method of MATLAB and COMSOL is employed to obtain a numerical solution. The results are verified with field measurements. In Section 3, based on this model, development patterns of a seepage and temperature field during freezing with different initial seepage velocities (at the step of 1 m/d) are quantitatively investigated. The effect of initial seepage velocity on the frozen wall closure time and effective thickness is analysed.

#### 2. Mathematical Model

Generally, soil in freezing construction is a multiphase system that consists of a soil skeleton, water, air, and ice [16]. As this paper focuses on the effect of shaft underground water seepage on the freezing of a sandy soil layer, the soil is assumed to remain saturated during freezing. Considering the research focus and primary influencing factors, the following assumptions for the mathematical model derivation are presented:(1)The effect of temperature variation and solute transport-induced density gradient and temperature gradient on the transport of unfrozen water is disregarded. The heat transfer loss of a freezing pipe wall is disregarded.(2)Unfrozen water transport always satisfies the water flow continuity condition and Darcy’s law. The initial state is steady seepage.(3)A soil layer is a continuous, homogeneous and isotropic porous medium, which satisfies the basic assumption of mixture theory.(4)The soil skeleton gap is a constant that is unrelated to thermodynamics parameters, ignoring the influence of temperature and pressure on viscosity coefficient of fluid.

##### 2.1. Mathematical Expression for Content of Soil Composition during Freezing

The porosity of the unfrozen state is ; assume that the pore water volume is ; the volumes of soil compositions are listed in Table 1.

After soil freezes, it always contains some liquid water, and the content maintains a dynamic balance with temperature [16]. Based on considerable experimental analysis, the content is determined via formula (1) [16]. Therefore, the pore water volume is represented via formula (2):where is the initial water volume content in soil; is the volume content of unfrozen water in frozen soil; is the soil temperature, °C; is the initial soil freezing temperature, °C; and *b* is a constant determined by soil property.

##### 2.2. Creation of Hydrothermal Coupling Control Equation during Freezing

The temperature variation during freezing is a transient heat conduction problem. Based on the basic assumption, ice, water, and soil skeleton coexist in soil during freezing and exhibit different heat transfer characteristics. In addition to heat transfer by conduction, ice heat transfer also includes phase change heat transfer. Water heat transfer includes convective heat transfer, which is induced by liquid water flow. Based on the heat transfer and seepage theory, the temperature and seepage field control equation during freezing that considers a phase change is expressed as [17]where is the time; , , and are the specific heat capacities of ice, water, and the soil skeleton; , , are the thermal conductivities of ice, water, and the soil skeleton; , , are the densities of ice, water, and the soil skeleton; is the phase change latent heat of the unit mass of water; is the relative velocity vector of water; is the strength of the equivalent heat source; is the compression rate of water; is the equivalent compression rate; is the absolute porosity; is the seepage pressure; and is the flow rate.

Based on the relation between permeability and absolute porosity from the servo permeability test [18] and Darcy’s law, the hydraulic conductivity during freezing is represented aswhere is the hydraulic conductivity, is the hydraulic conductivity for the unfrozen state, and is the viscosity coefficient of water.

In saturated soil, the absolute porosity is equal to the volume of liquid water. Therefore, disregarding the variation of the water viscosity coefficient during the cooling period, formulae (4) and (6) are represented as

The temperature boundary at the water inlet is the temperature of the water inflow, i.e., equal to the initial soil temperature. The temperature boundary at the freezing pipe wall is the actual temperature measured at the external wall of the freezing pipe. Assume that the convective flux is at the water outlet. The boundary condition iswhere is the horizontal coordinate of the outlet boundary, and is a vector along the inner normal direction at the water outlet.

As the heat transfer of other boundaries is not considered, other boundaries are treated as adiabatic boundaries. The boundary condition iswhere are the boundary vertical coordinates and is a vector along the inner normal direction of the boundary.

In the calculation, assume that an initial water flow only has a flow velocity from left to right and perpendicular to the boundary, as shown in Figure 1. To form a directional and steady flow velocity in the initial state, the initial condition iswhere is the initial hydraulic pressure distribution, which is obtained via the seepage velocity formula.

Both the inlet and outlet are set to boundaries with a fixed pressure. The corresponding values are determined by the initial pressure distribution. Other boundaries are set as the boundaries without flow:where is a vector along the inner normal direction of the boundary.

Formulae (3), (5), (7), and (8) and the boundary conditions comprise the mathematical model for hydrothermal coupling during freezing.

#### 3. Numerical Example

The mathematical model considers the variation of each composition, absolute porosity, and permeability for the seepage-temperature coupling effect during freezing. This model also considers the time-dependence of the temperature boundary of the freezing pipe wall. Therefore, the control equation and boundary condition are highly nonlinear. In this paper, based on the interpolation function and the coupling equation developed by MATLAB script language and the combined programming of MATLAB and COMSOL, the transient state solution of the mathematical model for hydrothermal coupling during freezing is obtained by substituting the variables and boundary conditions in the field equation.

Extensive engineering experience reveals that the effect of seepage on freezing is substantially dependent on a freezing plan. Therefore, a main shaft freezing project of a mine in Shangdong is presented as an example in this paper. First, the freezing temperature field for this condition is analysed to validate the mathematical model and solution. Second, based on the effect of different seepage velocities on the freezing process, the effect of the seepage velocity and variation on the closure time and frozen wall thickness during freezing are quantitatively analysed to obtain the critical seepage velocity and common law for the freezing plan. Last, the optimal path of the freezing-hole layout is proposed based on the simulations.

##### 3.1. Engineering Background and Geometrical Model

In the main shaft of a mine in Shangdong, the artificial freezing method is employed to penetrate a 457.78 m thick alluvium; the shaft net diameter is 4.5 m, and the excavation diameter is 7.8 m [19]. The designed active freezing period is 150 days. The diameter of the primary pipe freezing hole is 159 mm; the circle diameter is 16.5 m; 40 holes exist; and the hole spacing is 1.295 m. The secondary pipe freezing-hole diameter is 133 mm; the circle diameter is 11 m; 10 holes exist; and the hole spacing is 3.45 m [19].

In shaft freezing construction, the geometrical model for frozen wall formation is three-dimensional. If the freezing pipe deflection, nonuniform axial distribution of the freezing pipe wall temperature, and soil heterogeneity during shaft freezing are disregarded, the model is simplified to a two-dimensional model. Therefore, a horizontal section at a depth of 342.5 m is selected for investigation. The shaft centre is defined as the origin; both the model length and the model width are set to 36 m; the seepage direction is perpendicular to the vertical boundary and extends from left to right. The constructed geometrical model and grid for calculation are shown in Figure 1. In the figure, 03, C1, C2, C3, and C4 represent the temperature measurement points with equivalent distance to the shaft centre. Of these points, 03 is the actual temperature measurement hole; C2 and C3 are the intersection points of the freezing-hole boundary with its axial plane for the primary pipe and secondary pipe, respectively; C1 and C4 are the two positions in the axis planes of two secondary pipe freezing holes along the boundary and close to the primary pipe; and C0 is the downstream velocity measurement point, which is the intersection of the axis and the boundary of two freezing holes in the downstream primary pipe.

##### 3.2. Calculation Parameters

Based on actual engineering hydrological data, the seepage velocity is set to 0.5 m/d in the calculations. The initial hydraulic conductivity is 3.656 m/d. The initial soil temperature is 22°C. Assume that the freezing pipe wall does not experience any heat transfer loss and the temperature at pipe outer wall is the brine temperature in the pipe. In the numeric calculation, therefore, the pipe wall temperature is based on the field measurement of the brine temperature. The average field measurement of the brine temperature at the outlet and inlet versus time is shown in Figure 2.

Based on corresponding physical test results and research conclusions in reference [16], the physical parameters are listed as follows: the phase change latent heat of unit mass of water is 334 ; the initial freezing temperature of saturated soil is −0.02°C; and the constant *b* determined by the property of saturated soil is 0.61 [16]. Values of other physical parameters in formulae (3) and (6) are listed in Table 2.

##### 3.3. Comparison and Verification of Numerical Calculation Result versus Field Measurement

To validate mathematical model and numerical solution method proposed in this paper, primary pipe closure time, effective frozen wall thickness after active freezing, and field temperature measurement are compared and analysed.

When primary shaft freezes for 150 d, actual measurement of frozen wall thickness is 6.6 m; actual measurement of primary pipe closure time is 22 days [19]. As shown in Figures 3 and 4, from numerical calculation result, the average thickness of the effective frozen wall is 6.54 m after the crude diameter is deducted, and the primary pipe closure time is 22 days, both match well with field measurement. Due to limitation of paper length, only 03 measurement point is selected to compare with actual temperature measurement. Comparison curve of actual measurement data versus numerical simulation is shown in Figure 5. Figure 5 shows that the simulation result and field measurement result have similar temperature variation pattern and amplitude.

Comparison results indicate that three numerical calculation results are consistent with the field measurements, which indicates that the numerical calculation results are accurate and valid for extension research.

#### 4. Analysis of Frozen Wall Formation Pattern under Seepage

To investigate the effect of seepage velocity on frozen wall formation, based on the model, the effect of seepage field, temperature field, seepage velocity on the frozen wall closure time, and the effective thickness for different initial seepage velocities (at a step of 1 m/d) during freezing is quantitatively investigated.

During freezing, with the formation of a frozen area, the roundabout flow induced by the water insulation effect of the frozen area along the freezing front will change the seepage velocity and direction. These changes in seepage characteristics affect the coupling with the temperature field. Coupling of the seepage and temperature fields is the main influencing factor of different characteristics of frozen wall formation. To investigate the effect of seepage velocity on frozen wall formation, the characteristics of two fields in seepage-temperature field coupling during freezing should be separately analysed.

##### 4.1. Analysis of Seepage Field during Freezing

Prior to freezing, with the exception of the freezing hole with the roundabout flow, other regions are located in the seepage field with uniform velocity. As freezing continues, the frozen area gradually forms around the freezing hole and the roundabout flow occurs along the freezing front. Then, the frozen regions between the midstream freezing holes are separately connected. In the hole deployment area, a “window” with a wide top and narrow bottom exists. Water flows from the wide “window” and aggregates at the narrow “window”, which causes an increase in the flow velocity at the point. When the downstream “window” closes, the seepage velocity at this point rapidly decreases to zero.

As shown in Figure 6, as freezing continues, the seepage velocity at the downstream velocity measurement point gradually increases. Then, the seepage velocity rapidly approaches the peak value. Subsequently, the velocity gradually increases to the peak value and significantly decreases to zero. With an increase in the initial seepage velocity, this pattern changes to a certain extent. As listed in Table 3, the flow velocity duration downstream increases; and both the peak seepage velocity and its duration at the downstream velocity measurement point increase. The rate of increase of the peak value gradually increases; however, the ratio between the peak seepage velocity to the initial seepage velocity gradually decreases. When initial seepage velocity is 10 m/d, a decrease in the seepage velocity in the calculation period is not exhibits.

As shown in Figure 7, when the downstream “window” closes, a semicircular water insulation frozen area is formed in the hole area. At this moment, the water flow by passes the freezing-hole area, and the seepage velocity in the freezing-hole area is almost zero. The seepage velocity in upstream and downstream areas significantly decreases, which indicates that the water insulation effect of freezing construction is formed.

However, in the midstream area along the tangential direction of flow velocity, the flow velocity exceeds the initial seepage velocity and is approximately two times the initial seepage velocity. The amplitude of the flow velocity increases with an increase in the initial seepage velocity, as listed in Table 4.

This phenomenon indicates that during freezing of a double-circle-piped shaft, the maximum seepage velocity downstream of the freezing-hole deployment area significantly exceeds the initial stratum seepage velocity. The downstream “window closure” and water flow that circumvents the freezing-hole deployment area almost simultaneously occur. After the downstream “window closure”, the effect of the seepage velocity on freezing is only reflected in the outward expansion rate of the midstream freezing front. The effect on frozen wall formation in the freezing-hole deployment area is insignificant.

##### 4.2. Effect of Seepage Velocity on Freezing Temperature Field

The frozen wall formation process is the frozen area expansion process induced by the variation in temperature field. The essence of the formation rule is the variation rule of the freezing temperature field. Therefore, an analysis of the effect of seepage velocity on frozen wall formation requires further analysis of the effect of seepage velocity on the freezing temperature field.

An analysis of the temperature field calculation results indicates that when seepage velocity is not considered, the frozen area expands along the freezing pipe in the outer normal direction in the form of a concentric circle prior to the frozen wall closure. After the frozen wall closure, the frozen area expands along the boundary until a homogeneous frozen wall is formed. Due to the seepage-temperature field coupling effect, as heat transfer between freezing holes changes to the superposition of fluid heat transfer and solid heat transfer, the freezing temperature field demonstrates inhomogeneity. Inhomogeneity increases with seepage flow velocity. As freezing continues, inhomogeneity gradually decreases after the frozen wall closure. At a particular moment when the seepage velocity increases, the frozen wall undergoes the following sequential states: formation of the entire frozen wall, midstream and downstream “window closure”, and frozen wall midstream closure.

To intuitively analyse the effect of the initial seepage velocity in different strata on the freezing temperature field, the calculation results for four key nodes with the same freezing duration are compared and analysed, as shown in Figure 8.

**(a)**

**(b)**

**(c)**

**(d)**

As shown in the nephogram in Figure 8, for the same freezing duration, when the initial stratum seepage velocity is 0 m/d, the primary pipe is closed and closure of the primary and secondary pipes starts, the temperature field is symmetric and homogeneous in the frozen area. When the initial stratum seepage velocity is 5 m/d, the frozen area demonstrates a pattern that displays a thin top and a thick bottom. The temperature of the lower part of midstream is significantly lower than the temperature of the other part. As the seepage velocity increases, the frozen area shrinks, the temperature in the hole area rises and inhomogeneity increases.

The temperature curve in Figure 8 demonstrates the following results. (1) When the flow velocity is zero, two measurement points in the midstream area have identical temperature variation patterns; the temperature variation patterns for the measurement points upstream and downstream are identical. As the measurement points are located in different positions, the temperature of two measurement points in the midstream area is lower than the temperature of the measurement points upstream and downstream. The difference gradually increases when the primary and secondary pipe closure starts and gradually disappears when closure is completed. (2) When the seepage velocity is a nonzero value, the temperature variation at each measurement point demonstrates four distinct phases as follows: (a) before formation of the frozen area around the freezing hole near the measurement point, the temperature rapidly decreases. (b) after formation of the frozen area around the freezing hole near measurement point and before closure of the frozen area, temperature at the measurement point gradually decreases. (c) after the downstream “window closure” and before total frozen wall formation, the temperature at the measurement point rapidly decreases. (d) after total frozen wall formation, the rate of temperature decrease at the measurement point decelerates; the difference between two temperatures at different measurement points is minimal; and the temperature at the downstream measurement point gradually becomes lower than the temperature at other measurement points. The temperatures at the measurement points maintain the following order: C4 < C3 < C1 < C2. (3) As the seepage velocity increases, the durations of phases b and c are prolonged. During a long period, the temperatures at different measurement points maintain the following order: C3 < C2 < C4 < C1. (4) The interval for the temperature at each measurement point to decrease to the phase change temperature increases. However, the increase for each interval at the upstream and downstream measurement points is insignificant. For different seepage velocities, following the order C3, C2, C4, C1, the patterns of temperature variation curve at each measurement point becomes more distinct.

The seepage and temperature field analysis results are described as follows: the seepage-temperature field coupling effect, the primary pipe at the lower part of the midstream area closes. The primary and secondary pipes also closed at the primary surface of the lower part of the midstream area, and the closure expands to the upstream and downstream area as freezing continues. After downstream closes, the impact of the seepage velocity on the freezing process in the freezing-hole deployment area is insignificant. The temperature in the frozen area rapidly decreases, and the difference among the temperatures at different points gradually decreases.

##### 4.3. Variation Pattern of Closure Time for Different Initial Seepage Velocities

Analysis of the individual characteristics of two fields with the seepage-temperature field coupling effect during freezing reveals that the seepage-temperature field coupling effect causes a loss of cooling energy and nonuniform temperature field distribution. As the initial stratum seepage velocity increases, the frozen wall inhomogeneity increases. Midstream freezing holes are gradually connected, and the closure expands towards the downstream and upstream areas until the frozen area between the upstream freezing holes closes to form the entire frozen area. Therefore, the upstream frozen wall closure time is used to represent the formation of the entire frozen area. The variation curve of the frozen wall closure time versus the initial stratum seepage velocity is shown in Figure 9.

As shown in Figure 9, when the initial stratum seepage velocity is within 7 m/d, the increased closure time is always less than 10% of the active freezing period for each 1 m/d increase in the initial stratum seepage velocity. However, when the initial stratum seepage velocity increases from 8 m/d to 9 m/d, the increased closure time is 27% of the active freezing period. When the initial stratum seepage velocity is 10 m/d, the frozen wall cannot close. The upstream and downstream frozen wall closure times increase with the initial stratum seepage velocity and both demonstrate a quasiexponential relation. However, the interval between the upstream closure time and the downstream closure time is not distinct.

The cause of this phenomenon is described as follows: after the midstream frozen area is formed, the downstream frozen area closes, and the resultant water insulation effect reduces the loss of cooling energy in the freezing-hole deployment area and facilitates closure of the entire frozen wall.

##### 4.4. Variation Pattern of Maximum Thickness of Frozen Wall during Active Freezing Period

Due to the underground water seepage effect, before the frozen wall closure, the upstream, midstream, and downstream frozen wall thicknesses along the seepage direction demonstrate a pattern with a thin top and a thick bottom. As the seepage velocity increases, the effective frozen walls in the upstream, midstream, and downstream areas become thinner, which causes an upstream closure failure, and the effective frozen wall cannot be formed. After the frozen wall closure, due to the water insulation effect by the frozen wall, the seepage velocities in the frozen wall interior area and the frozen wall exterior upstream and downstream areas are significantly lower than the initial stratum seepage velocity. However, the seepage velocity in the frozen wall exterior midstream area is approximately two times the initial stratum seepage velocity. Therefore, the freezing front expansion in the frozen wall interior area and the as frozen wall exterior upstream and downstream areas are not affected by the seepage velocity. The freezing front line at the frozen wall interior gradually forms a regular circle. The freezing front expansion in the frozen wall exterior upstream and downstream areas exhibits a similar pattern. However, the frozen wall exterior midstream area remains affected by high flow velocity and expansion decelerates.

Two completely different freezing front expansion patterns before and after the frozen wall closure cause a change in the frozen wall thickness growth trend with an increase in the seepage velocity and freezing time. To reflect the effect of the seepage velocity on freezing, the relation between the frozen wall maximum thickness without removing the crude diameter and the seepage velocity during the active freezing period are analysed in this paper, as shown in Figure 10.

As shown in Figure 10, when freezing progresses to 150 days, with the exception of a seepage velocity of zero, the average thickness of the entire frozen wall linearly decreases with an increase in the seepage velocity. The maximum thickness of the downstream frozen wall always exceeds the maximum thickness of the upstream and midstream frozen walls. The difference between the upstream thickness and the midstream thickness is maintained within 0.5 m. The difference between the downstream value and the upstream and midstream values range from 1.63 to 2.69 m. The maximum thickness of the downstream frozen wall initially increases with the seepage velocity and then decreases; the peak is attained when the seepage velocity is 2 m/d. The maximum thicknesses of the upstream and midstream frozen walls decrease with an increase in the seepage velocity. When the seepage velocity is less than 7 m/d, the thickness of the upstream wall is slightly greater than the thickness of the midstream wall. When the seepage velocity exceeds 7 m/d, the thickness of the midstream wall exceeds the thickness of the upstream wall.

The cause of this phenomenon is described as follows: during the same period, the cooling energy that is removed by the seepage effect is positively related to the seepage velocity. Prior to the frozen wall closure, due to the underground water seepage effect, the temperature of the water flow gradually reduces through upstream and midstream. The disturbance of low-temperature water in the downstream area is beneficial for increasing the frozen wall thickness. However, as the seepage velocity increases, the cooling energy removed from the frozen area increases; meanwhile, the decreasing trend of low-temperature water is also limited. Therefore, the beneficial effect gradually diminishes, and the seepage velocity eventually becomes an unfavourable factor for the downstream frozen wall thickness. An increase in the seepage velocity has the following effects: prior to the frozen wall closure, the frozen wall has a thin top and a thick bottom; after the frozen wall closure, the frozen wall closure time increased and the expansion period of the freezing front in the frozen wall interior and exterior in the upstream area decreases. However, when the seepage velocity exceeds 7 m/d, offsetting the difference between the thin top and the thick bottom prior to the frozen wall closure is difficult.

Due to the seepage-temperature field coupling effect, as the seepage velocity increases, the maximum thickness of the downstream frozen wall initially increases and then decreases; however, the average thickness of the frozen wall linearly decreases. In this freezing plan, the critical seepage velocity is 7 m/d.

#### 5. Conclusions and Suggestions

In this paper, the hydrothermal coupling process is analysed. Based on previous research, a hydrothermal full-parameter coupling equation is derived; its solution is provided via the combined programming method of MATLAB and COMSOL. First, a mine project is investigated to validate the mathematical model and solution. Second, the effect of the initial stratum seepage flow velocity on frozen wall formation is investigated and the following conclusions are obtained.(1)During freezing, the maximum seepage velocity that the downstream hole area should resist is significantly larger than the initial stratum seepage velocity. The maximum seepage velocity decreases with an increase in the initial seepage velocity. At the end of the active freezing period, the maximum seepage velocity at the midstream peak surface is approximately two times the initial seepage velocity.(2)Due to the seepage-temperature field coupling effect, at the lower part of the midstream area, the primary pipe closes, which is the same case for the primary and secondary pipes. After the downstream closure, the frozen area temperature rapidly decreases. The differences among the temperatures at different positions gradually decrease.(3)As the initial stratum seepage velocity increases, the freezing temperature field inhomogeneity increases. The maximum thickness of the downstream frozen wall initially increases and then decreases. However, the average thickness of the frozen wall linearly decreases. The upstream and downstream frozen wall closure time and initial stratum seepage velocity follow a quasiexponential relation and variation of closure time interval is not conspicuous. In this freezing plan, the critical seepage velocity is 7 m/d.

Therefore, to reduce the frozen wall closure time and improve the efficiency of the total freezing construction, the recommended design principles of the freezing plan are as follows: first, the downstream area should be closed to enable the water insulation effect to reduce the outflow of cooling energy; and second, the upstream closure should be expedited to reduce the closure time of the entire frozen wall.

#### Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

The work presented in this paper is financially supported by the National Natural Science Foundation of China (Nos. 51804006, 51778004, 51674006, 51474004 and 51374010), National Science and Technology Major Project of China (No. 2016ZX05068-001), and Youth Science Research Foundation of Anhui University of Science and Technology (Nos. QN2017211 and QN2017222). The authors gratefully acknowledge financial support of the above-mentioned agencies.